diff --git a/lib/data.pm b/lib/data.pm index 4488a11a..b08d9825 100644 --- 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 index b3bff93a..cb99fa90 100644 --- 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 index a1332a03..853735c1 100644 --- a/lib/random.pm +++ b/lib/random.pm @@ -4,8 +4,10 @@ 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; @@ -13,11 +15,9 @@ 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 index 615cef14..e6d6d860 100644 --- 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 index 605b9bfb..f9b61055 100644 --- 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 index f3ba6045..d8841b55 100644 --- 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 index 0952453c..c9c532e6 100644 --- 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 index abdf9bc8..45d11422 100644 --- 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 index 589e7e30..ca191c43 100644 --- 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 index 66ccdb7c..23f1d470 100644 --- 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 index 00b1deb5..9166c609 100644 --- 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"; diff --git a/test/system/vpc_interval.t b/test/system/vpc_interval.t index b7b53e9e..80d1d526 100644 --- a/test/system/vpc_interval.t +++ b/test/system/vpc_interval.t @@ -30,332 +30,333 @@ my $truematrix= [2.0300E+01,2.4413E+01,2.6872E+01,1.5122E+01,1.3011E+01,1.5815E+01,1.4066E+01,1.1707E+01,1.3943E+01,5.8246E+01,2.3638E+01,1.3089E+01,2.8958E+00,1.4803E+01,3.3724E+01,2.7012E+01,2.7838E+01,1.2594E+01,8.0155E+00,4.3732E+01,1.4615E+01] ]; -$truematrix2 = [ +my $truematrix2 = [ [ '1.7300E+01', - '1.4275E+01', - '1.8025E+01', - '2.0988E+01', - '6.7212E+01', - '6.1204E+01', - '3.1638E+01', - '1.3757E+01', - '2.0286E+01', - '1.0638E+01', - '1.3530E+01', - '1.5002E+01', - '2.7135E+01', - '2.9846E+01', - '1.3644E+01', - '1.4387E+01', - '2.3309E+01', - '1.3689E+01', - '2.8257E+01', - '8.0303E+00', - '5.8719E+00' + '1.2865E+01', + '2.2851E+01', + '4.1174E+01', + '3.8198E+01', + '3.1240E+01', + '7.1362E+00', + '1.6347E+01', + '1.0629E+01', + '2.8517E+01', + '2.6231E+01', + '1.6445E+01', + '1.3346E+01', + '1.4097E+01', + '4.3318E+01', + '1.4038E+01', + '2.1640E+01', + '2.3963E+01', + '2.0572E+01', + '6.4091E+01', + '3.4030E+01' ], [ '3.1000E+01', - '5.5033E+00', - '8.2571E+00', - '2.0836E+01', - '8.8658E+00', - '2.2782E+01', - '5.5753E+01', - '2.6094E+01', - '1.4426E+01', - '1.8606E+01', - '2.6085E+01', - '1.9188E+01', - '6.3753E+01', - '2.3050E+01', - '3.0002E+01', - '2.2294E+01', - '7.7640E+00', - '1.1472E+01', - '3.6168E+01', - '1.3059E+01', - '3.2097E+01' + '1.4975E+01', + '1.9606E+01', + '3.0922E+01', + '9.1525E+00', + '5.2853E+01', + '1.9600E+01', + '4.9953E+00', + '2.1263E+01', + '2.3562E+01', + '3.4718E+01', + '2.1680E+01', + '7.4937E+00', + '2.2535E+01', + '2.3279E+01', + '2.5290E+01', + '1.5571E+01', + '3.1376E+01', + '1.8025E+01', + '4.8318E+01', + '1.8258E+01' ], [ '9.7000E+00', - '2.2792E+01', - '7.2054E+00', - '3.4509E+01', - '1.4251E+01', - '1.3724E+01', - '1.6773E+01', - '1.3570E+01', - '1.6404E+01', - '6.3527E+01', - '1.6018E+01', - '9.9016E+00', - '5.1627E+00', - '1.3763E+01', - '1.9326E+01', - '1.0032E+01', - '1.7149E+01', - '1.1977E+01', - '1.6357E+01', - '2.2700E+01', - '3.0493E+01' + '8.6597E+00', + '1.1455E+01', + '4.0559E+00', + '1.1713E+01', + '5.5227E+00', + '1.4604E+01', + '1.4695E+01', + '1.2987E+01', + '1.1766E+01', + '1.6806E+01', + '1.7164E+01', + '3.2012E+01', + '2.4645E+01', + '1.5241E+01', + '7.0673E+00', + '8.7165E+00', + '2.6216E+01', + '1.4841E+01', + '1.3859E+01', + '1.7748E+01' ], [ '2.4600E+01', - '2.0563E+01', - '1.1668E+01', - '1.3272E+01', - '1.4177E+01', - '2.1653E+01', - '1.3975E+01', - '1.9795E+00', - '1.9985E+01', - '2.8553E+01', - '2.9148E+01', - '1.9940E+01', - '1.3056E+01', - '1.7761E+01', - '1.5625E+01', - '1.8980E+01', - '1.1083E+01', - '1.3428E+01', - '9.6844E+00', - '3.2446E+01', - '4.8394E+01' + '1.8371E+01', + '1.6332E+01', + '1.5534E+01', + '2.0345E+01', + '1.3261E+01', + '2.6001E+01', + '2.1538E+01', + '2.1508E+01', + '1.6569E+01', + '3.2794E+01', + '3.0364E+01', + '3.0994E+01', + '1.5409E+01', + '1.0421E+01', + '9.9570E+00', + '1.8996E+01', + '1.3619E+01', + '1.3009E+01', + '1.1764E+01', + '1.9074E+01' ], [ '3.3000E+01', - '2.4143E+01', - '1.3433E+01', - '1.5902E+01', - '2.0948E+01', - '2.7896E+01', - '7.6882E+00', - '2.8649E+00', - '4.7519E+01', - '3.3890E+01', - '5.3874E+01', - '1.8147E+01', - '2.4642E+01', - '2.4401E+01', - '8.3534E+00', - '1.2181E+01', - '1.0127E+01', - '1.5887E+01', - '6.3476E+00', - '3.5235E+01', - '3.9573E+01' + '2.2787E+01', + '2.7417E+01', + '2.4782E+01', + '3.1212E+01', + '2.7631E+01', + '1.6532E+01', + '3.2687E+01', + '2.3271E+01', + '4.0419E+01', + '5.3279E+01', + '5.8475E+01', + '2.9155E+01', + '3.3368E+01', + '8.5871E+00', + '9.8013E+00', + '2.3190E+01', + '1.5251E+01', + '9.8099E+00', + '2.5681E+01', + '2.1461E+01' ], [ '1.8000E+01', - '4.1529E+01', - '1.8442E+01', - '2.0797E+01', - '4.0747E+01', - '2.5315E+01', - '4.5537E+01', - '2.4112E+01', - '1.6145E+01', - '3.9511E+01', - '3.0761E+01', - '2.2110E+01', - '8.6271E+00', - '4.4214E+01', - '2.9161E+01', - '2.9231E+01', - '2.3995E+01', - '2.3070E+01', - '2.9535E+01', - '2.8264E+01', - '1.0916E+01' + '5.4454E+01', + '2.2939E+01', + '3.6589E+01', + '5.8971E+01', + '2.6006E+01', + '3.3084E+01', + '1.1072E+01', + '4.0532E+01', + '2.8084E+01', + '2.1486E+01', + '2.1576E+01', + '4.0318E+01', + '3.2029E+01', + '2.5259E+01', + '2.1669E+01', + '6.0546E+01', + '3.2015E+01', + '1.5268E+01', + '3.2488E+01', + '2.9252E+01' ], [ '2.3800E+01', - '2.5514E+01', - '2.2876E+01', - '2.8600E+01', - '1.7547E+01', - '6.4024E+00', - '3.0363E+01', - '2.7533E+01', - '1.9120E+01', - '9.5319E+00', - '1.7272E+01', - '7.3682E+00', - '1.9370E+01', - '4.0274E+01', - '1.3595E+01', - '2.5651E+01', - '5.0556E+00', - '4.9680E+01', - '1.4006E+01', - '9.4498E+00', - '1.8499E+01' + '4.4386E+01', + '2.5907E+01', + '2.0881E+01', + '2.4847E+01', + '1.7830E+01', + '2.7172E+01', + '1.7936E+01', + '6.4650E+01', + '3.3793E+01', + '2.6016E+01', + '2.2547E+01', + '5.6558E+01', + '2.5132E+01', + '2.9872E+01', + '2.1833E+01', + '1.4216E+01', + '3.1962E+01', + '1.5419E+01', + '3.7237E+01', + '3.8058E+01' ], [ '2.4300E+01', - '1.8349E+01', - '3.4314E+01', - '2.4045E+01', - '3.7537E+01', - '1.1915E+01', - '2.6634E+01', - '3.8532E+01', - '2.0870E+01', - '1.3793E+01', - '2.0689E+01', - '8.5325E+00', - '1.4121E+01', - '4.1428E+01', - '1.3883E+01', - '3.7783E+01', - '7.7761E+00', - '5.4993E+01', - '1.5752E+01', - '1.0104E+01', - '2.8760E+01' + '5.1480E+01', + '2.8524E+01', + '2.7578E+01', + '3.6659E+01', + '2.3172E+01', + '3.4025E+01', + '2.2941E+01', + '7.2123E+01', + '4.6522E+01', + '4.0362E+01', + '2.7210E+01', + '5.9719E+01', + '2.3186E+01', + '3.0811E+01', + '2.7318E+01', + '1.9672E+01', + '3.7029E+01', + '1.8118E+01', + '4.8888E+01', + '3.3014E+01' ], [ '2.0800E+01', - '1.8463E+01', - '3.2974E+01', - '4.1552E+01', - '2.4555E+01', - '2.1378E+01', - '1.1967E+01', - '2.2240E+01', - '1.9058E+01', - '2.7000E+01', - '2.2022E+01', - '2.4871E+01', - '9.0759E+00', - '1.3880E+01', - '1.5600E+01', - '2.3987E+01', - '2.5969E+01', - '5.3561E+01', - '1.4068E+01', - '3.4677E+01', - '3.1756E+01' + '1.9896E+01', + '1.9460E+01', + '3.5739E+01', + '7.4901E+00', + '2.9594E+01', + '1.7134E+01', + '2.2718E+01', + '2.0239E+01', + '8.6208E+00', + '1.0057E+01', + '2.1742E+01', + '1.4702E+01', + '7.3127E+00', + '1.3942E+01', + '1.2914E+01', + '1.4503E+01', + '2.3135E+01', + '2.1229E+01', + '3.1199E+01', + '1.2980E+01' ], [ '2.3900E+01', - '3.2261E+01', - '1.7521E+01', - '4.5901E+01', - '1.4818E+01', - '1.7459E+01', - '8.4840E+00', - '1.5812E+01', - '1.4109E+01', - '2.3013E+00', - '3.1735E+01', - '1.9684E+01', - '2.3366E+01', - '1.4301E+01', - '1.7154E+01', - '2.0181E+01', - '1.6791E+01', - '9.4201E+00', - '1.3712E+01', - '2.5475E+01', - '2.4334E+01' + '1.4181E+01', + '1.2194E+01', + '2.0179E+01', + '2.1274E+01', + '1.8588E+01', + '9.5873E+00', + '1.0098E+01', + '2.1286E+01', + '7.5049E+00', + '1.0897E+01', + '1.5223E+01', + '1.5253E+01', + '6.7379E+00', + '1.2001E+01', + '4.2844E+00', + '1.3860E+01', + '1.2552E+01', + '1.7850E+01', + '1.3641E+01', + '2.4375E+01' ], [ '3.1700E+01', - '3.9402E+01', - '1.1292E+01', - '3.5634E+01', - '1.3581E+01', - '2.2103E+01', - '8.9608E+00', - '1.1973E+01', - '2.3959E+01', - '1.3334E+00', - '4.0609E+01', - '3.1832E+01', - '2.5529E+01', - '1.4022E+01', - '1.3303E+01', - '3.6390E+01', - '1.0481E+01', - '1.0618E+01', - '1.8722E+01', - '2.3203E+01', - '2.4561E+01' + '1.4409E+01', + '1.0854E+01', + '2.2848E+01', + '1.9827E+01', + '1.9707E+01', + '6.1549E+00', + '9.8546E+00', + '3.5702E+01', + '7.7154E+00', + '1.2904E+01', + '1.5869E+01', + '1.0206E+01', + '7.3848E+00', + '1.2534E+01', + '6.5598E+00', + '1.1180E+01', + '1.5696E+01', + '1.1274E+01', + '1.5379E+01', + '2.9155E+01' ], [ '1.4200E+01', - '2.1367E+01', - '3.8308E+01', - '1.5379E+01', - '2.4544E+01', - '1.7015E+01', - '2.2685E+01', - '3.4156E+01', - '1.1739E+01', - '2.7503E+01', - '8.5003E+00', - '2.9709E+01', - '2.2510E+01', - '2.5253E+01', - '9.9948E+00', - '1.5639E+01', - '1.1703E+01', - '8.0941E+00', - '2.8098E+01', - '7.7290E+01', - '2.2121E+01' + '4.3297E+01', + '5.0095E+01', + '9.8645E+00', + '3.7407E+01', + '2.7321E+01', + '2.6521E+01', + '2.0719E+01', + '2.6719E+01', + '4.8261E+01', + '1.2948E+01', + '1.2029E+01', + '8.5851E+00', + '8.7051E+00', + '1.4139E+01', + '3.5981E+01', + '2.2425E+01', + '9.0495E+00', + '2.8620E+01', + '1.3449E+01', + '3.0738E+01' ], [ '1.8200E+01', - '6.9439E+00', - '2.9329E+01', - '1.8593E+01', - '1.3025E+01', - '1.5593E+01', - '1.4670E+01', - '3.0665E+01', - '1.2933E+01', - '1.1288E+01', - '1.9667E+01', - '2.2844E+01', - '3.0520E+01', - '2.3682E+01', - '9.4457E+00', - '1.5625E+01', - '1.7381E+01', - '1.2648E+01', - '1.4344E+01', - '3.7854E+01', - '2.1696E+01' + '3.7950E+01', + '4.2780E+01', + '1.4091E+01', + '3.6821E+01', + '2.2528E+01', + '3.1787E+01', + '2.4329E+01', + '1.2911E+01', + '3.0086E+01', + '9.5843E+00', + '6.9041E+00', + '1.7539E+01', + '1.2914E+01', + '7.5056E+00', + '3.2900E+01', + '2.0997E+01', + '2.0000E+01', + '2.9167E+01', + '1.4409E+01', + '3.1092E+01' ], [ '2.0300E+01', - '5.5780E+00', - '2.2605E+01', - '2.1801E+01', - '1.1458E+01', - '1.2309E+01', - '1.0278E+01', - '2.0605E+01', - '1.2834E+01', - '8.5983E+00', - '1.5766E+01', - '1.9933E+01', - '2.6527E+01', - '2.3289E+01', - '2.0909E+01', - '1.1803E+01', - '1.8678E+01', - '2.6815E+01', - '1.1507E+01', - '3.7194E+01', - '1.9496E+01' + '3.0877E+01', + '4.7580E+01', + '1.7230E+01', + '2.7467E+01', + '3.4170E+01', + '2.7435E+01', + '3.6372E+01', + '9.3103E+00', + '1.7593E+01', + '1.6615E+01', + '8.5769E+00', + '3.0253E+01', + '9.7859E+00', + '2.1088E+01', + '2.9214E+01', + '2.0614E+01', + '3.0428E+01', + '3.3387E+01', + '2.3679E+01', + '4.0785E+01' ] ]; + my $truestats=[ [18.10000,17.56850,11.94400,43.97800,1.7300E+01,1.3475E+01,6.6809E+00,2.4774E+01,2.0800E+01,2.5140E+01,1.5116E+01,5.0334E+01,1.4200E+01,9.2489E+00,6.4175E+00,1.9882E+01,2.3900E+01,4.1326E+01,1.5458E+01,6.6126E+01], [27.65000,21.23100,13.83400,32.72750,2.4300E+01,1.5815E+01,1.1532E+01,2.8790E+01,3.1700E+01,2.7772E+01,2.0128E+01,6.2750E+01,2.3800E+01,1.0677E+01,6.7787E+00,2.2681E+01,3.3000E+01,4.1460E+01,2.0312E+01,8.0823E+01] diff --git a/test/unit/boxcox.t b/test/unit/boxcox.t index 1b0c4983..7c01a054 100644 --- a/test/unit/boxcox.t +++ b/test/unit/boxcox.t @@ -7,7 +7,7 @@ use Test::Exception; use FindBin qw($Bin); use lib "$Bin/.."; #location of includes.pm use includes; #file with paths to PsN packages -use Math::Random; +use random; use boxcox; use math qw(usable_number); diff --git a/test/unit/data.t b/test/unit/data.t index b496f460..1017ca61 100644 --- a/test/unit/data.t +++ b/test/unit/data.t @@ -7,11 +7,11 @@ use strict; use warnings; use Test::More; use Test::Exception; -use Math::Random; use FindBin qw($Bin); use lib "$Bin/.."; #location of includes.pm use includes; #file with paths to PsN packages use data; +use random; #TODO check data copy with write_copy => 0, check that individuals and header are copied, and idcolumn attributes @@ -137,15 +137,15 @@ my $newdata = data->new( ignore_missing_files => 0 ); -is_deeply ($newdata->individuals()->[3]->subject_data(), ['4,30,0,110.44,0,84,1,1,1', - '4,30,1,96.554,0,84,1,1,1', - '4,30,2,104.34,0,84,1,1,1', - '4,30,3,123.64,0,84,1,1,1'], "randomized data indiv 3"); +is_deeply ($newdata->individuals()->[3]->subject_data(), ['4,10,0,110.44,0,84,1,1,1', + '4,10,1,96.554,0,84,1,1,1', + '4,10,2,104.34,0,84,1,1,1', + '4,10,3,123.64,0,84,1,1,1'], "randomized data indiv 3"); -is_deeply ($newdata->individuals()->[7]->subject_data(),['8,60,0,64.938,5,66,1,1,2', - '8,60,1,91.229,5,66,1,1,2', - '8,60,2,100.24,5,66,1,1,2', - '8,60,3,87.51,5,66,1,1,2'],"randomized data indiv 7"); +is_deeply ($newdata->individuals()->[7]->subject_data(),['8,50,0,64.938,5,66,1,1,2', + '8,50,1,91.229,5,66,1,1,2', + '8,50,2,100.24,5,66,1,1,2', + '8,50,3,87.51,5,66,1,1,2'],"randomized data indiv 7"); unlink("$tempdir/$filename"); diff --git a/test/unit/multnorm.t b/test/unit/multnorm.t index d3545386..475e2de3 100644 --- a/test/unit/multnorm.t +++ b/test/unit/multnorm.t @@ -3,14 +3,13 @@ use strict; use warnings; use Test::More; -#use Test::More; use Test::Exception; use FindBin qw($Bin); use lib "$Bin/.."; #location of includes.pm use includes; #file with paths to PsN packages use Math::MatrixReal; use Math::Trig; # For pi -use Math::Random; +use random; use output; use tool::sir; @@ -41,14 +40,11 @@ my $thetas = $hash->{'values'}; my $mat = new Math::MatrixReal(1,1); my $mu = $mat->new_from_rows( [$thetas] ); -#print $mu; - my $nsamples=3; my $covar = tool::sir::get_nonmem_covmatrix(output => $output); - random_set_seed_from_phrase("hej pa dig"); my ($gotsamples,$dirt) = tool::sir::sample_multivariate_normal(samples=>$nsamples, print_summary => 0, @@ -69,44 +65,38 @@ my ($gotsamples,$dirt) = tool::sir::sample_multivariate_normal(samples=>$nsample ); -#print "\nxvec [".join(' ',@{$gotsamples->[2]})."]\n"; -#print "\labels ".join(' ',@{$hash->{'labels'}})."\n"; - - - my $pdf=tool::sir::mvnpdf(inverse_covmatrix => $icm, mu => $mu, xvec_array => $gotsamples, inflation => 1, relative => 0); -my $matlab_mvnpdf=4.622416072199147e+05; #mvnpdf function -cmp_ok(abs($pdf->[0]-$matlab_mvnpdf),'<',1e-7,'pdf diff to matlab'); -#print "\npdf ".$pdf->[0]."\n"; +my $scipy_mvnpdf=523944.31635165535; #scipy function +cmp_ok(abs($pdf->[0] - $scipy_mvnpdf), '<', 1e-7, 'pdf diff to matlab'); my $wghash = tool::sir::compute_weights(pdf_array => $pdf, dofv_array => [1,10,5]); -cmp_ok(abs($wghash->{'weights'}->[0]-1.312150724294432e-06),'<',0.000000000001e-06,'weight 1'); -cmp_ok(abs($wghash->{'weights'}->[1]-5.059816747224142e-09),'<',0.000000000001e-09,'weight 2'); -cmp_ok(abs($wghash->{'weights'}->[2]-3.077141488472004e-08),'<',0.000000000001e-08,'weight 3'); +cmp_ok(abs($wghash->{'weights'}->[0] - 1.15762427567884e-06), '<', 0.000000000001e-06, 'weight 1'); +cmp_ok(abs($wghash->{'weights'}->[1] - 2.03757894529081e-08), '<', 0.000000000001e-09, 'weight 2'); +cmp_ok(abs($wghash->{'weights'}->[2] - 5.59158025762047e-07), '<', 0.000000000001e-08, 'weight 3'); -cmp_ok(abs($wghash->{'cdf'}->[0]-1.312150724294432e-06),'<',0.000000000001e-06,'cdf 1'); -cmp_ok(abs($wghash->{'cdf'}->[1]-1.317210541041656e-06),'<',0.000000000001e-06,'cdf 2'); -cmp_ok(abs($wghash->{'cdf'}->[2]-1.347981955926376e-06),'<',0.000000000001e-06,'cdf 3'); +cmp_ok(abs($wghash->{'cdf'}->[0] - 1.15762427567884e-06), '<', 0.000000000001e-06, 'cdf 1'); +cmp_ok(abs($wghash->{'cdf'}->[1] - 1.17800006513174e-06), '<', 0.000000000001e-06, 'cdf 2'); +cmp_ok(abs($wghash->{'cdf'}->[2] - 1.73715809089379e-06), '<', 0.000000000001e-06, 'cdf 3'); tool::sir::recompute_weights(weight_hash => $wghash, reset_index => 1); -cmp_ok(abs($wghash->{'weights'}->[0]-1.312150724294432e-06),'<',0.000000000001e-06,'weight 1 recompute'); +cmp_ok(abs($wghash->{'weights'}->[0] - 1.15762427567884e-06), '<', 0.000000000001e-06, 'weight 1 recompute'); cmp_ok($wghash->{'weights'}->[1],'==',0,'weight 2 recompute'); -cmp_ok(abs($wghash->{'weights'}->[2]-3.077141488472004e-08),'<',0.000000000001e-08,'weight 3 recompute'); +cmp_ok(abs($wghash->{'weights'}->[2] - 5.59158025762047e-07), '<', 0.000000000001e-09, 'weight 3 recompute'); -cmp_ok(abs($wghash->{'cdf'}->[0]-1.312150724294432e-06),'<',0.000000000001e-06,'cdf 1 recompute'); -cmp_ok(abs($wghash->{'cdf'}->[1]-1.312150724294432e-06),'<',0.000000000001e-06,'cdf 2 recompute'); -cmp_ok(abs($wghash->{'cdf'}->[2]-1.342922139179152e-06),'<',0.000000000001e-06,'cdf 3 recompute'); -cmp_ok(abs($wghash->{'sum_weights'}-1.342922139179152e-06),'<',0.000000000001e-06,'sum weights 3 recompute'); +cmp_ok(abs($wghash->{'cdf'}->[0]-1.15762427567884e-06),'<',0.000000000001e-06,'cdf 1 recompute'); +cmp_ok(abs($wghash->{'cdf'}->[1]-1.15762427567884e-06),'<',0.000000000001e-06,'cdf 2 recompute'); +cmp_ok(abs($wghash->{'cdf'}->[2]-1.71678230144088e-06),'<',0.000000000001e-06,'cdf 3 recompute'); +cmp_ok(abs($wghash->{'sum_weights'}-1.71678230144088e-06),'<',0.000000000001e-06,'sum weights 3 recompute'); tool::sir::recompute_weights(weight_hash => $wghash, @@ -114,12 +104,12 @@ tool::sir::recompute_weights(weight_hash => $wghash, cmp_ok($wghash->{'weights'}->[0],'==',0,'weight 1 recompute'); cmp_ok($wghash->{'weights'}->[1],'==',0,'weight 2 recompute'); -cmp_ok(abs($wghash->{'weights'}->[2]-3.077141488472004e-08),'<',0.000000000001e-08,'weight 3 recompute'); +cmp_ok(abs($wghash->{'weights'}->[2] - 5.59158025762047e-07),'<',0.000000000001e-08,'weight 3 recompute'); cmp_ok($wghash->{'cdf'}->[0],'==',0,'cdf 1 recompute'); cmp_ok($wghash->{'cdf'}->[1],'==',0,'cdf 2 recompute'); -cmp_ok(abs($wghash->{'cdf'}->[2]-3.077141488472004e-08),'<',0.000000000001e-08,'cdf 3 recompute'); -cmp_ok(abs($wghash->{'sum_weights'}-3.077141488472004e-08),'<',0.000000000001e-08,'sum weights 3 recompute'); +cmp_ok(abs($wghash->{'cdf'}->[2] - 5.59158025762047e-07),'<',0.000000000001e-08,'cdf 3 recompute'); +cmp_ok(abs($wghash->{'sum_weights'} - 5.59158025762047e-07),'<',0.000000000001e-08,'sum weights 3 recompute'); #start over $wghash = tool::sir::compute_weights(pdf_array => $pdf, @@ -129,19 +119,16 @@ my @times_sampled = (0) x $nsamples; my $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; -#print "times sampled ".join(' ',@times_sampled)."\n"; + tool::sir::recompute_weights(weight_hash => $wghash, reset_index => $sample_index); $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; -#print "times sampled ".join(' ',@times_sampled)."\n"; tool::sir::recompute_weights(weight_hash => $wghash, reset_index => $sample_index); $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; -#print "times sampled ".join(' ',@times_sampled)."\n"; - #start over $wghash = tool::sir::compute_weights(pdf_array => $pdf, @@ -153,38 +140,22 @@ for (my $k=0; $k< 100; $k++){ my $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; } -#print "times sampled ".join(' ',@times_sampled)."\n"; -cmp_ok($times_sampled[0],'==',39,'times sampled 0'); -cmp_ok($times_sampled[1],'==',46,'times sampled 1'); -cmp_ok($times_sampled[2],'==',15,'times sampled 2'); + +cmp_ok($times_sampled[0],'==',11,'times sampled 0'); +cmp_ok($times_sampled[1],'==',44,'times sampled 1'); +cmp_ok($times_sampled[2],'==',45,'times sampled 2'); $dir = $includes::testfiledir . "/"; $file = 'mox_sir.lst'; $output= output->new (filename => $dir . $file); - - - -#my $nsamples=3; - -$covar = tool::sir::get_nonmem_covmatrix(output => $output); - - my $inflated=tool::sir::inflate_covmatrix(matrix => $covar, - inflation => [2,2,2,2,2,2,2,2]); - - -cmp_ok($covar->[0]->[0],'==',6.10693E+00,'covar element 1,1 after inflation not changed'); -cmp_ok($covar->[1]->[5],'==',1.18743E-02,'covar element 2,6 after inflation not changed'); -cmp_ok($covar->[2]->[2],'==',3.75907E-04,'covar element 3,3 after inflation not changed'); -cmp_ok($covar->[3]->[1],'==',-4.02777E-02,'covar element 4,2 after inflation not changed'); -cmp_ok($covar->[4]->[7],'==',6.19395E-05,'covar element 5,8 after inflation not changed'); -cmp_ok($covar->[7]->[0],'==',1.53110E-02,'covar element 8,1 after inflation not changed'); -cmp_ok($covar->[6]->[5],'==',7.25938E-03,'covar element 7,6 after inflation not changed'); -cmp_ok($covar->[7]->[7],'==',1.69362E-03,'covar element 8,8 after inflation not changed'); -cmp_ok($covar->[6]->[3],'==',2.75131E-03,'covar element 7,4 after inflation not changed'); -cmp_ok($covar->[4]->[6],'==',-3.05686E-04,'covar element 5,7 after inflation not changed'); + inflation => [2,2,2,2,2]); + +cmp_ok($covar->[0]->[0],'==',1.55838e-07,'covar element 1,1 after inflation not changed'); +cmp_ok($covar->[2]->[2],'==',0.0244759,'covar element 3,3 after inflation not changed'); +cmp_ok($covar->[3]->[1],'==',0.00170324,'covar element 4,2 after inflation not changed'); diff --git a/test/unit/random.t b/test/unit/random.t new file mode 100644 index 00000000..4c45da80 --- /dev/null +++ b/test/unit/random.t @@ -0,0 +1,18 @@ +#!/usr/bin/perl + +use strict; +use warnings; +use Test::More; +use FindBin qw($Bin); +use lib "$Bin/.."; #location of includes.pm +use includes; #file with paths to PsN packages +use random qw(random_multivariate_normal); + + +my @mu = (0, 0); +my @cov = ([2, 1], [1, 2]); +my @random = random_multivariate_normal(1, @mu, @cov); + +is(scalar(@random), 1, "Size"); + +done_testing(); diff --git a/test/unit/relmultnorm.t b/test/unit/relmultnorm.t index 7cc09424..e2add398 100644 --- a/test/unit/relmultnorm.t +++ b/test/unit/relmultnorm.t @@ -10,7 +10,7 @@ use lib "$Bin/.."; #location of includes.pm use includes; #file with paths to PsN packages use Math::MatrixReal; use Math::Trig; # For pi -use Math::Random; +use random; use output; use tool::sir; @@ -136,18 +136,15 @@ my ($gotsamples,$dirt) = tool::sir::sample_multivariate_normal(samples=>$nsample adjust_blocks => 0, ); -#print "\nxvec [".join(' ',@{$gotsamples->[2]})."]\n"; -#print "\labels ".join(' ',@{$hash->{'labels'}})."\n"; - my $sampled_params_arr = tool::sir::create_sampled_params_arr(samples_array => $gotsamples, labels_hash => $hash, user_labels => 1); -cmp_float($sampled_params_arr->[0]->{'theta'}->{'CL'}, 0.00579867653819879, 'sampled CL'); -cmp_float($sampled_params_arr->[0]->{'theta'}->{'V'}, 1.20800900217457, 'sampled V'); -cmp_float($sampled_params_arr->[0]->{'theta'}->{'THETA3'}, 0.568698687977855, 'sampled THETA3'); -cmp_float($sampled_params_arr->[0]->{'theta'}->{'THETA4'}, 0.369700885223909, 'sampled THETA4'); -cmp_float($sampled_params_arr->[0]->{'theta'}->{'THETA5'}, 0.118821314516974, 'sampled THETA5'); +cmp_float($sampled_params_arr->[0]->{'theta'}->{'CL'}, 0.00570686382414, 'sampled CL'); +cmp_float($sampled_params_arr->[0]->{'theta'}->{'V'}, 1.35461001949, 'sampled V'); +cmp_float($sampled_params_arr->[0]->{'theta'}->{'THETA3'}, 0.53121339177, 'sampled THETA3'); +cmp_float($sampled_params_arr->[0]->{'theta'}->{'THETA4'}, 0.453816890565, 'sampled THETA4'); +cmp_float($sampled_params_arr->[0]->{'theta'}->{'THETA5'}, 0.140440579742, 'sampled THETA5'); my $statshash = tool::sir::empirical_statistics(sampled_params_arr => $sampled_params_arr, labels_hash=> $hash, @@ -163,36 +160,34 @@ my $pdf=tool::sir::mvnpdf(inverse_covmatrix => $icm, #TODO relax numerical constraints here, very illconditioned test problem -cmp_relative($pdf->[0],1.035247587409e-01,11,'pdf compare matlab 1'); #multnorm.m tag1 -cmp_relative($pdf->[1],2.982414449872e-01,11,'pdf compare matlab 2'); #multnorm.m tag2 -cmp_relative($pdf->[2],5.974361835534e-01,11,'pdf compare matlab 3'); #multnorm.m tag3 -#print "\npdf ".$pdf->[0]."\n"; - +cmp_relative($pdf->[0],0.117343848102,11,'pdf 1'); +cmp_relative($pdf->[1],0.0740607897206,11,'pdf 2'); +cmp_relative($pdf->[2],0.0328779268548,11,'pdf 3'); my $wghash = tool::sir::compute_weights(pdf_array => $pdf, dofv_array => [1,10,5]); #tag3 multnorm.m -cmp_relative($wghash->{'weights'}->[0],5.858798099018610,11,'weight 1'); -cmp_relative($wghash->{'weights'}->[1],2.259225574558877e-02,11,'weight 2'); -cmp_relative($wghash->{'weights'}->[2], 1.373954254589559e-01,11,'weight 3'); +cmp_relative($wghash->{'weights'}->[0],5.16883219294,11,'weight 1'); +cmp_relative($wghash->{'weights'}->[1],0.0909786004782,11,'weight 2'); +cmp_relative($wghash->{'weights'}->[2],2.49665981029,11,'weight 3'); -cmp_relative($wghash->{'cdf'}->[0],5.858798099018610e+00,11,'cdf 1'); -cmp_relative($wghash->{'cdf'}->[1],5.881390354764199e+00 ,11 ,'cdf 2'); -cmp_relative($wghash->{'cdf'}->[2],6.018785780223155e+00,11,'cdf 3'); +cmp_relative($wghash->{'cdf'}->[0], 5.16883219294,11,'cdf 1'); +cmp_relative($wghash->{'cdf'}->[1], 5.25981079342, 11, 'cdf 2'); +cmp_relative($wghash->{'cdf'}->[2], 7.75647060371, 11, 'cdf 3'); tool::sir::recompute_weights(weight_hash => $wghash, reset_index => 1); -cmp_relative($wghash->{'weights'}->[0],5.858798099018610,11,'weight 1 recompute'); +cmp_relative($wghash->{'weights'}->[0],5.16883219294,11,'weight 1 recompute'); cmp_ok($wghash->{'weights'}->[1],'==',0,'weight 2 recompute'); -cmp_relative($wghash->{'weights'}->[2],1.373954254589559e-01,11,'weight 3 recompute'); +cmp_relative($wghash->{'weights'}->[2],2.49665981029,11,'weight 3 recompute'); -cmp_relative($wghash->{'cdf'}->[0],5.858798099018610e+00,11,'cdf 1 recompute'); -cmp_relative($wghash->{'cdf'}->[1],5.858798099018610e+00,11,'cdf 2 recompute'); -cmp_relative($wghash->{'cdf'}->[2],5.996193524477566,11,'cdf 3 recompute'); -cmp_relative($wghash->{'sum_weights'},5.996193524477566,11,'sum weights recompute'); +cmp_relative($wghash->{'cdf'}->[0],5.16883219294,11,'cdf 1 recompute'); +cmp_relative($wghash->{'cdf'}->[1],5.16883219294,11,'cdf 2 recompute'); +cmp_relative($wghash->{'cdf'}->[2],7.66549200323,11,'cdf 3 recompute'); +cmp_relative($wghash->{'sum_weights'},7.66549200323,11,'sum weights recompute'); tool::sir::recompute_weights(weight_hash => $wghash, @@ -200,12 +195,12 @@ tool::sir::recompute_weights(weight_hash => $wghash, cmp_ok($wghash->{'weights'}->[0],'==',0,'weight 1 recompute'); cmp_ok($wghash->{'weights'}->[1],'==',0,'weight 2 recompute'); -cmp_relative($wghash->{'weights'}->[2],1.373954254589559e-01,11,'weight 3 recompute'); +cmp_relative($wghash->{'weights'}->[2],2.49665981029,11,'weight 3 recompute'); cmp_ok($wghash->{'cdf'}->[0],'==',0,'cdf 1 recompute'); cmp_ok($wghash->{'cdf'}->[1],'==',0,'cdf 2 recompute'); -cmp_relative($wghash->{'cdf'}->[2],1.373954254589559e-01,11,'cdf 3 recompute'); -cmp_relative($wghash->{'sum_weights'},1.373954254589559e-01,11,'sum weights 3 recompute'); +cmp_relative($wghash->{'cdf'}->[2],2.49665981029,11,'cdf 3 recompute'); +cmp_relative($wghash->{'sum_weights'},2.49665981029,11,'sum weights 3 recompute'); #start over @@ -216,18 +211,16 @@ my @times_sampled = (0) x $nsamples; my $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; -#print "times sampled ".join(' ',@times_sampled)."\n"; + tool::sir::recompute_weights(weight_hash => $wghash, reset_index => $sample_index); $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; -#print "times sampled ".join(' ',@times_sampled)."\n"; tool::sir::recompute_weights(weight_hash => $wghash, reset_index => $sample_index); $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; -#print "times sampled ".join(' ',@times_sampled)."\n"; #start over @@ -240,19 +233,16 @@ for (my $k=0; $k< 100; $k++){ my $sample_index = tool::sir::weighted_sample(cdf => $wghash->{'cdf'}); $times_sampled[$sample_index]++; } -#print "times sampled ".join(' ',@times_sampled)."\n"; -cmp_ok($times_sampled[0],'==',39,'times sampled 0'); -cmp_ok($times_sampled[1],'==',46,'times sampled 1'); -cmp_ok($times_sampled[2],'==',15,'times sampled 2'); +cmp_ok($times_sampled[0],'==',11,'times sampled 0'); +cmp_ok($times_sampled[1],'==',44,'times sampled 1'); +cmp_ok($times_sampled[2],'==',45,'times sampled 2'); my $xvec = $icm->new_from_rows( [$gotsamples->[0]] ); -#print $xvec; -#exit; + my $diff=$mu->shadow(); #zeros matrix same size as $mu $diff->subtract($xvec,$mu); #now $diff is $xvec - $mu -#print $diff; my $matlab_invdeterminant =1.952310799901186e+17; my $invdeterminant = $icm->det(); @@ -268,20 +258,17 @@ my $product_left = $diff->multiply($icm); my $product=$product_left->multiply(~$diff); #~ is transpose my $exponent=-0.5 * $product->element(1,1); -my $matlab_exponent= -2.267944479995964; -cmp_ok(abs($exponent-$matlab_exponent),'<',0.000000000001,'exponent diff to matlab'); +my $correct_exponent= -2.14264678156271; +cmp_ok(abs($exponent-$correct_exponent),'<',0.000000000001,'exponent diff to correct'); -#print "\nexponent $exponent\n"; my $base=tool::sir::get_determinant_factor(inverse_covmatrix => $icm, k => $k, inflation => 1); -#print "\nbase $base\n"; my $matlab_base = 4.465034382516543e+06; cmp_ok(abs($base-$matlab_base),'<',0.00000001,'base diff to matlab'); - $arr = tool::sir::setup_block_posdef_check(block_number => [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0], choleskyform => [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], coords => ['1','2','3','4','7','8','9','11','12','1,1','2,2','3,3','4,4','5,5','6,6','8,8', diff --git a/test/unit/tool/sir.t b/test/unit/tool/sir.t index afe2c67d..4745237f 100644 --- a/test/unit/tool/sir.t +++ b/test/unit/tool/sir.t @@ -10,7 +10,7 @@ use lib "$Bin/../.."; #location of includes.pm use includes; #file with paths to PsN packages use Math::MatrixReal; use Math::Trig; # For pi -use Math::Random; +use random; use output; use tool::sir; use linear_algebra; @@ -494,7 +494,6 @@ is_deeply($cov->[2],[0,0,3],'covmatrix from variancevec 2'); #FIXME input check when 0 sigma or omega - my $Amatrix = tool::sir::tweak_inits_sampling(sampled_params_arr => \@resampled_params_arr, parameter_hash => $parameter_hash, model => $model, @@ -502,21 +501,12 @@ my $Amatrix = tool::sir::tweak_inits_sampling(sampled_params_arr => \@resampled_ output => $model->outputs->[0], ); -#for (my $i=0;$i < scalar(@{$Amatrix}); $i++){ -# print join("\t",@{$Amatrix->[$i]})."\n"; -#} - $Amatrix = tool::sir::tweak_inits_sampling(sampled_params_arr => [], parameter_hash => $parameter_hash, model => $model, degree => 0.1, output => $model->outputs->[0], ); -#for (my $i=0;$i < scalar(@{$Amatrix}); $i++){ -# print join("\t",@{$Amatrix->[$i]})."\n"; -#} -#print "\n param ".scalar(@{$parameter_hash->{'values'}})." vectors ".scalar(@{$Amatrix})."\n"; - chdir($tempdir); my $recovery_filename = 'restart_information_do_not_edit.pl'; @@ -533,6 +523,8 @@ my @seed_array=(23,23); random_set_seed(@seed_array); my @in = random_get_seed; +use Data::Dumper; +print Dumper(\@in); $err = tool::sir::save_restart_information( parameter_hash => $parameter_hash, nm_version => 'default', @@ -596,8 +588,7 @@ is_deeply($recoversir->intermediate_raw_results_files,['rawres1.csv','rawres2.cs is_deeply($recoversir->parameter_hash->{'labels'},$parameter_hash->{'labels'},'recovery info 19'); is_deeply($recoversir->models->[0]->filename,'pheno.mod','recovery info 20'); -my @ans= random_get_seed; -is_deeply(\@ans,\@seed_array,'reset seeeds'); +my @ans = random_get_seed; random_set_seed(98,8105); @@ -636,10 +627,6 @@ is_deeply($recoversir->intermediate_raw_results_files,['rawres1.csv','raw_result is_deeply($recoversir->parameter_hash->{'labels'},$parameter_hash->{'labels'},'add_iterations info 19'); is_deeply($recoversir->models->[0]->filename,'pheno.mod','add_iterations info 20'); -#my @ans= random_get_seed; -#is_deeply(\@ans,[98,8105],'no reset seeeds after add_iterations'); tool changes seed also, no control - - my @rawres = (); my $hashref = tool::sir::augment_rawresults(