Remaining switch to Math::Random::Free as default

UUPharmacometrics · Nov 29, 2023 · f7f23f1 · f7f23f1
1 parent 655342c
commit f7f23f1
Show file tree

Hide file tree

Showing 18 changed files with 438 additions and 439 deletions.
diff --git a/lib/data.pm b/lib/data.pm
@@ -4,7 +4,7 @@ use include_modules;
 use OSspecific;
 use File::Copy "cp";
 use Config;
-use Math::Random;
+use random;
 use Storable;
 use ui;
 use status_bar;

diff --git a/lib/model.pm b/lib/model.pm
@@ -9,7 +9,7 @@ use Config;
 use OSspecific;
 use Storable;
 use model::shrinkage_module;
-use Math::Random qw(random_multivariate_normal);
+use random qw(random_multivariate_normal);
 use model::iofv_module;
 use model::nonparametric_module;
 use model::annotation;

diff --git a/lib/random.pm b/lib/random.pm
@@ -4,20 +4,20 @@ use strict;
 use Digest::SHA qw(sha1_hex);
 use Exporter qw(import);
 use PsN;
+use linear_algebra;
 
-our @EXPORT = qw(random_set_seed_from_phrase random_uniform_integer random_uniform random_permuted_index random_get_seed random_set_seed);
+our @EXPORT = qw(random_set_seed_from_phrase random_uniform_integer random_normal random_uniform random_permuted_index random_get_seed random_set_seed);
+our @EXPORT_OK = qw(random_multivariate_normal);
 
 our $old;
 our @global_seed;
 
 
 if (PsN::use_old_math_random) {
     require Math::Random;
-    require Math::Random::Free;     # FIXME: Remove this
     $old = 1;
 } else {
     require Math::Random::Free;
-    require Math::Random;           # FIXME: Remove this
     $old = 0;
 }
 
@@ -27,7 +27,7 @@ sub random_set_seed
     if ($old) {
         Math::Random::random_set_seed(@_);
     } else {
-        my ($s1, $s2) = @_; 
+        my ($s1, $s2) = @_;
         my $seed = $s1 ^ $s2;
         $global_seed[0] = $seed;
         $global_seed[1] = 0;
@@ -50,17 +50,13 @@ sub random_set_seed_from_phrase
 {
     if ($old) {
         Math::Random::random_set_seed_from_phrase(@_); 
-        Math::Random::Free::random_set_seed_from_phrase(@_);    # FIXME: Remove this
     } else {
-        #Math::Random::Free::random_set_seed_from_phrase(@_);
-        Math::Random::random_set_seed_from_phrase(@_);      # FIXME: Remove this
+        my ($seed) = @_;
+        my $value = hex substr(sha1_hex($seed), 0, 6);
+        $global_seed[0] = $value;
+        $global_seed[1] = 0;
+        srand $value;
     }
-    my( $seed ) = @_;
-    # On 64-bit machine the max. value for srand() seems to be 2**50-1
-    my $value = hex substr( sha1_hex( $seed ), 0, 6 );
-    $global_seed[0] = $value;
-    $global_seed[1] = 0;
-    srand $value;
 }
 
 
@@ -84,6 +80,16 @@ sub random_uniform
 }
 
 
+sub random_normal
+{
+    if ($old) {
+        return Math::Random::random_normal(@_); 
+    } else {
+        return Math::Random::Free::random_normal(@_);
+    }
+}
+
+
 sub random_permuted_index
 {
     if ($old) {
@@ -94,4 +100,34 @@ sub random_permuted_index
 }
 
 
+sub random_multivariate_normal
+{
+    if ($old) {
+        return Math::Random::random_multivariate_normal(@_);
+    } else {
+        my $n = shift;
+        my @results;
+        my $size = scalar(@_) / 2;
+        my @mu = @_[0..$size-1];
+        my @cov = @_[$size..scalar(@_) - 1];
+        use Storable qw(dclone);
+        my @covcopy = @{dclone(\@cov)};
+        linear_algebra::cholesky(\@covcopy);
+        for (my $i = 0; $i < $n; $i++) {
+            my @standard_normal = random_normal($size, 0, 1);
+            my @vec;
+            for (my $row = 0; $row < $size; $row++) {
+                my $s = $mu[$row];
+                for (my $col = 0; $col <= $row; $col++) {
+                    $s += $covcopy[$col]->[$row] * $standard_normal[$col];
+                }
+                push @vec, $s;
+            }
+            push @results, \@vec;
+        }
+        return @results;
+    }
+}
+
+
 1;
diff --git a/lib/tool.pm b/lib/tool.pm
@@ -8,7 +8,7 @@ use File::Path qw(mkpath rmtree);
 use File::Spec;
 use Scalar::Util qw(blessed);
 use OSspecific;
-use Math::Random;
+use random;
 use Archive::Zip;
 use utils::file;
 use PsN;

diff --git a/lib/tool/bootstrap.pm b/lib/tool/bootstrap.pm
@@ -4,7 +4,7 @@ use include_modules;
 use strict;
 use File::Copy 'cp';
 use Statistics::Distributions 'udistr', 'uprob';
-use Math::Random;
+use random;
 use Mouse;
 use MouseX::Params::Validate;
 use array;

diff --git a/lib/tool/modelfit.pm b/lib/tool/modelfit.pm
@@ -9,7 +9,7 @@ use File::Path;
 use File::Glob;
 use File::Spec;
 use Storable;
-use Math::Random;
+use random;
 use nonmemrun;
 use nonmemrun::localunix;
 use nonmemrun::localwindows;

diff --git a/lib/tool/sir.pm b/lib/tool/sir.pm
@@ -6,7 +6,7 @@ use File::Copy 'cp';
 use data;
 use OSspecific;
 use tool::modelfit;
-use Math::Random qw(random_get_seed random_set_seed random_uniform random_multivariate_normal);
+use random qw(random_get_seed random_set_seed random_uniform random_multivariate_normal);
 use Mouse;
 use MouseX::Params::Validate;
 use Math::MatrixReal;
@@ -1426,7 +1426,6 @@ sub mvnpdf
             push(@pdf_array,$det_factor*exp(-0.5 * $product->element(1,1))); #matlab absolute mvnpdf possibly change back
         }
     }
-
     return \@pdf_array;
 }
 

diff --git a/lib/tool/sse.pm b/lib/tool/sse.pm
@@ -5,7 +5,7 @@ use model;
 use tool::cdd;
 use tool::modelfit;
 use model_transformations;
-use Math::Random;
+use random;
 use Config;
 use Mouse;
 use MouseX::Params::Validate;

diff --git a/test/system/modelfit.t b/test/system/modelfit.t
@@ -55,9 +55,9 @@ my @command_line = (
 	get_command('execute') . " -seed=1 phenomaxeval10.mod -no-disp -clean=1 -no-tweak_inits -min_retries=0 -retries=0 -maxevals=9999 -handle_msfo -directory=dir0",
 	get_command('execute') . " -seed=1  phenomaxeval10.mod -no-disp -clean=1 -tweak_inits  -min_retries=0 -retries=0 -maxevals=0 -directory=dir1",
 	get_command('execute') . " -seed=1  phenomaxeval10.mod -no-disp -clean=1 -tweak_inits  -min_retries=0 -retries=1 -maxevals=0 -directory=dir2",
-	get_command('execute') . " -seed=1  pheno.mod -clean=1 -no-disp -tweak_inits  -min_retries=0 -retries=1 -reduced_model_ofv=70 -directory=dir3",
+	get_command('execute') . " -seed=52  pheno.mod -clean=1 -no-disp -tweak_inits  -min_retries=0 -retries=1 -reduced_model_ofv=70 -directory=dir3",
 	get_command('execute') . " -seed=1  phenomaxeval0.mod -no-disp -clean=1 -tweak_inits  -min_retries=0 -retries=1 -maxevals=0 -directory=dir4",
-	get_command('execute') . " -seed=1  phenomaxeval0.mod -no-disp -clean=1 -tweak_inits  -min_retries=1 -retries=0 -maxevals=0 -directory=dir5",
+	get_command('execute') . " -seed=154  phenomaxeval0.mod -no-disp -clean=1 -tweak_inits  -min_retries=1 -retries=0 -maxevals=0 -directory=dir5",
 	get_command('execute') . " -seed=1  create_tnpri_msf.mod  -no-disp",
 	get_command('execute') . " -seed=1  tnpri.mod -clean=1 -tweak_inits  -no-disp -min_retries=0 -retries=1 -extra_files=msf_tnpri -directory=dir6",
 	get_command('execute') . " -seed=1  tnpri.mod -clean=1 -tweak_inits  -no-disp -min_retries=1 -retries=0 -extra_files=msf_tnpri -maxevals=0 -directory=dir7",

diff --git a/test/system/psn_clean.t b/test/system/psn_clean.t
@@ -20,7 +20,7 @@ our $tempdir = create_test_dir('system_psn_clean');
 
 copy_test_files($tempdir, ["pheno.mod", "pheno.dta"]);
 chdir($tempdir);
-system(get_command('execute')." pheno.mod -dir=mydir -seed=1 -min_ret=1 -tweak -no-disp -clean=0");
+system(get_command('execute')." pheno.mod -dir=mydir -seed=52 -min_ret=1 -tweak -no-disp -clean=0");
 ok (-e "mydir/NM_run1", "NM_run1 ok");
 ok (-e "mydir/NM_run1/patab1", "patab1 ok");
 ok (((-e "mydir/NM_run1/psn-1.mod") or (-e "mydir/NM_run1/psn-1.ctl")) , "psn-1.mod ok");

diff --git a/test/system/rawres_input.t b/test/system/rawres_input.t
@@ -20,12 +20,12 @@ my $model_dir = $includes::testfiledir;
 
 my @commands = 
 	(
-	 get_command('sir') . " pheno.mod -rawres_input=raw_pheno_for_rawres_input.csv -samples=10 -offset_rawres=12 -auto_rawres=0.2 -seed=2000 -resamples=5 ",
-	 get_command('sir') . " pheno.mod -rawres_input=raw_pheno_for_rawres_input.csv -samples=10 -offset_rawres=1 -resamples=5 ",
-	 get_command('parallel_retries') . " pheno.mod -samples=2 -rawres_input=raw_pheno_for_rawres_input.csv -no-display",
-	 get_command('sse') . " pheno.mod -rawres_input=raw_pheno_for_rawres_input.csv -samples=5 -no-est ",
-	 get_command('sse') . " pheno.mod  -samples=20 -dir=$ssedir",
-	 get_command('vpc') . " pheno.mod -rawres_input=$ssedir/raw_results_pheno.csv -in_filter=problem.gt.0 -offset=0 -samples=20 -auto_bin=unique ");
+     get_command('sir') . " pheno.mod -rawres_input=raw_pheno_for_rawres_input.csv -samples=10 -offset_rawres=12 -auto_rawres=0.2 -seed=200 -resamples=5 ",
+     get_command('sir') . " pheno.mod -rawres_input=raw_pheno_for_rawres_input.csv -samples=10 -offset_rawres=1 -resamples=5 ",
+     get_command('parallel_retries') . " pheno.mod -samples=2 -rawres_input=raw_pheno_for_rawres_input.csv -no-display",
+     get_command('sse') . " pheno.mod -rawres_input=raw_pheno_for_rawres_input.csv -samples=5 -no-est ",
+     get_command('sse') . " pheno.mod  -samples=20 -dir=$ssedir",
+     get_command('vpc') . " pheno.mod -rawres_input=$ssedir/raw_results_pheno.csv -in_filter=problem.gt.0 -offset=0 -samples=20 -auto_bin=unique ");
 
 foreach my $command (@commands) {
 	print "Running $command\n";