From 2772d322ed2bdfa27d94927d9f54427e75506cee Mon Sep 17 00:00:00 2001 From: pontus_pih Date: Wed, 11 Jun 2008 13:42:30 +0000 Subject: [PATCH] Fixed $PRIOR support and bugfix in shrinkage results handling --- Makefile | 1 + bin/extended_grid | 86 ++++++++++- bin/setup.pl | 2 +- diagrams/model.dia | 204 ++++++++++++++++++++++++++ lib/model/problem/record/theta_option_subs.pm | 2 +- lib/model/problem_subs.pm | 4 +- lib/model_subs.pm | 26 ++++ lib/output/problem_subs.pm | 24 +-- lib/tool/modelfit_subs.pm | 2 + 9 files changed, 330 insertions(+), 21 deletions(-) diff --git a/Makefile b/Makefile index a9c8c82..9217dfa 100644 --- a/Makefile +++ b/Makefile @@ -51,6 +51,7 @@ LIBFILES=debug.pm \ model/problem/nonparametric.pm \ model/problem/omega.pm \ model/problem/pk.pm \ + model/problem/prior.pm \ model/problem/problem.pm \ model/problem/pred.pm \ model/problem/scatter.pm \ diff --git a/bin/extended_grid b/bin/extended_grid index 62f0a2d..860e157 100644 --- a/bin/extended_grid +++ b/bin/extended_grid @@ -177,13 +177,14 @@ $model_copy -> set_records( type => 'table', $model_copy -> _write; my $id_to_sim; +my $nr_of_individuals = scalar @{$model_copy -> datas -> [0] -> individuals}; if( defined $options{'ids_to_simulate'} ){ $id_to_sim = $options{'ids_to_simulate'}; my $data_set = $model_copy -> datas -> [0]; - if( $id_to_sim <= scalar @{$data_set -> individuals} ){ + if( $id_to_sim <= $nr_of_individuals ){ my $individuals = $data_set -> individuals; my @new_individuals = @{$individuals}[0..($id_to_sim-1)]; @@ -439,7 +440,90 @@ $modelfit -> run; # }}} +my $nptab = table_file -> new( filename => "extended.nptab" ); + +my $JD_column = $nptab -> column_to_array( column => 1 + $nr_of_old_etas ); + +for( $nr_of_individuals .. ($nr_of_individuals + $id_to_sim-1) ){ + + $JD_column -> [$_] = $JD_column -> [$_] - (1/($nr_of_individuals + $id_to_sim)); + +} + +my $JD_sum = 0; +my $verification_sum = 0; + +for( 0 .. $#{$JD_column} ){ + + $JD_sum += $JD_column -> [$_]; + + $JD_column -> [$_] = $JD_column -> [$_] * ( ( $nr_of_individuals + $id_to_sim ) / $nr_of_individuals ); + + $verification_sum += $JD_column -> [$_]; + +} + +if( abs($verification_sum - 1) >= 0.01 ){ + print "Warning: sum of scaled joint densities is to inacurate: $verification_sum (should be between 0.99 and 1.01)n"; +} + +my @final_nptab; + +foreach my $eta_index ( 0..($nr_of_old_etas-1) ){ + + my $eta = $nptab -> column_to_array( column => $eta_index + 1 ); # plus 1 is for ID column + + my %index; + my $counter = 0; + foreach my $val ( @{$eta} ){ + $index{$counter} = $val; + $counter++; + } + + my @order = sort( {$index{$a} <=> $index{$b}} keys %index ); + + my @JD_etax; + my @sorted_eta; + + foreach my $index ( @order ){ + + push( @JD_etax, $JD_column -> [$index] ); + push( @sorted_eta, $eta -> [$index] ); + } + + my $previous_value = $JD_etax[0]; + + my @JD_etax_cumulative; + $JD_etax_cumulative[0] = $previous_value; + + for( 1 .. $#JD_etax ){ + + $previous_value = $previous_value + $JD_etax[$_]; + + push( @JD_etax_cumulative, $previous_value ); + + } + + unshift( @sorted_eta, "sorted_ETA" . ($eta_index+1) ); + unshift( @JD_etax_cumulative, "cumulative_JD". ($eta_index+1) ); + + push( @final_nptab, \@sorted_eta, \@JD_etax_cumulative ); + +} + +unshift( @{$JD_column}, 'NJD' ); + +unshift( @final_nptab, $JD_column ); + +dump_T_csv( 'final_nptab.csv', \@final_nptab ); + + + +################################################################################ +# # We are done for now, code below needs a specification and a review +# +################################################################################ exit; # {{{ Results handling (Cut'n'Paste from nonp_bootstrap_v2) diff --git a/bin/setup.pl b/bin/setup.pl index a20a856..62995d2 100644 --- a/bin/setup.pl +++ b/bin/setup.pl @@ -3,7 +3,7 @@ use Config; use CPAN; use File::Spec; -my $version = '2.2.6_rc1'; +my $version = '2.2.6_rc2'; my $name_safe_version = $version; $name_safe_version =~ s/\./_/g; diff --git a/diagrams/model.dia b/diagrams/model.dia index 2adbecb..8680123 100644 --- a/diagrams/model.dia +++ b/diagrams/model.dia @@ -14606,6 +14606,210 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + ## + + + ## + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ## + + + + + + + + ## + + + ## + + + + + + + + + + + ## + + + ## + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + #prior# + + + ## + + + ## + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/lib/model/problem/record/theta_option_subs.pm b/lib/model/problem/record/theta_option_subs.pm index e1373b0..d166ffd 100644 --- a/lib/model/problem/record/theta_option_subs.pm +++ b/lib/model/problem/record/theta_option_subs.pm @@ -36,7 +36,7 @@ start _read_option # use Data::Dumper; # print Dumper $PsN::config; $self -> {'lobnd'} = $PsN::config -> {'low_INF'}; - print $self -> {'lobnd'}; +# print $self -> {'lobnd'}; } if( defined $self -> {'upbnd'} and $self -> {'upbnd'} =~ /INF/ ){ diff --git a/lib/model/problem_subs.pm b/lib/model/problem_subs.pm index ae9879e..8a676c0 100644 --- a/lib/model/problem_subs.pm +++ b/lib/model/problem_subs.pm @@ -2,8 +2,8 @@ start include statements use Data::Dumper; -my @print_order = ('problem','abbreviated','input','data','msfi','contr','subroutine','model','tol','infn','omega','pk','aesinit','aes','des','error','pred','mix','theta','sigma','simulation','estimation','covariance','nonparametric','table','scatter'); -my @sde_print_order = ('problem','abbreviated','input','data','msfi','contr','subroutine','model','tol','infn','theta','omega','sigma','pk','aesinit','aes','des','error','pred','mix','simulation','estimation','covariance','nonparametric','table','scatter'); +my @print_order = ('problem','abbreviated','input','data','msfi','contr','subroutine','prior','model','tol','infn','omega','pk','aesinit','aes','des','error','pred','mix','theta','sigma','simulation','estimation','covariance','nonparametric','table','scatter'); +my @sde_print_order = ('problem','abbreviated','input','data','msfi','contr','subroutine','prior','model','tol','infn','theta','omega','sigma','pk','aesinit','aes','des','error','pred','mix','simulation','estimation','covariance','nonparametric','table','scatter'); my %abbreviations; # Here we intialize a hash used to find long names for abbreviated diff --git a/lib/model_subs.pm b/lib/model_subs.pm index 1941067..88271df 100644 --- a/lib/model_subs.pm +++ b/lib/model_subs.pm @@ -364,6 +364,32 @@ end shrinkage_stats # }}} shrinkage_stats +start shrinkage_modules +{ + if( defined $parm ){ + if( ref $parm ne 'ARRAY' + or + not ( scalar @{$parm} == scalar @{$self -> {'problems'}} ) ){ + 'debug' -> die( message => 'New number of shrinkage modules must be equal to number of problems' ); + } + + foreach my $prob( @{$self -> {'problems'}} ){ + my $new_module = shift( @{$parm} ); + $new_module -> model( $self ); + $prob -> shrinkage_module( shift( @{$parm} ) ); + + } + + } else { + my @return_array; + foreach my $prob( @{$self -> {'problems'}} ){ + push( @return_array, $prob -> shrinkage_module ); + } + return \@return_array; + } +} +end shrinkage_modules + # {{{ wres_shrinkage =head2 wres_shrinkage diff --git a/lib/output/problem_subs.pm b/lib/output/problem_subs.pm index b63857a..f263bf9 100644 --- a/lib/output/problem_subs.pm +++ b/lib/output/problem_subs.pm @@ -346,7 +346,6 @@ start _read_block_structures my $start_pos = $self -> {'lstfile_pos'}; my $success = 1; - my $tbarea = 0; my $obarea = 0; my $sbarea = 0; @@ -360,19 +359,10 @@ start _read_block_structures last; } - if ( /^0LENGTH OF THETA/ ) { - # If we find this we must find other stuff too. Set - # success = 0 and change this if we find the rest. Nope, - # this is not true. For example, you can supply only a - # $THETA if you use the LIKELIHOOD option for the - # $ESTIMATION block. - $tbarea = 1; - } - - if ( /^0INITIAL ESTIMATE/ ) { - # We want to find this if we are currently reading theta, - # omega or sigma block structures - $success = 1 if ( $tbarea or $sbarea or $obarea ); + if ( /^0INITIAL ESTIMATE/ or /^0DEFAULT OMEGA BOUNDARY TEST OMITTED:/ ) { + # We want to find this if we are currently reading omega + # or sigma block structures + $success = 1 if ( $sbarea or $obarea ); $start_pos --; last; } @@ -389,14 +379,12 @@ start _read_block_structures if(/0OMEGA HAS BLOCK FORM:/) { $self -> {'omega_block_structure_type'} = 'BLOCK'; $obarea = 1; - $tbarea = 0; $success = 1; next; } if(/0SIGMA HAS BLOCK FORM:/) { $self -> {'sigma_block_structure_type'} = 'BLOCK'; $sbarea = 1; - $tbarea = 0; $obarea = 0; $success = 1; next; @@ -718,6 +706,10 @@ while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) { return; } + if(/^ PRIOR SUBROUTINE USER-SUPPLIED/){ + $success = 1; + } + if(/^0ESTIMATION STEP OMITTED:\s*\b(.*)\b/) { $self -> estimation_step_initiated(1); $self -> estimation_step_run(0) if $1 eq 'YES'; diff --git a/lib/tool/modelfit_subs.pm b/lib/tool/modelfit_subs.pm index 7fa8c13..1532b08 100644 --- a/lib/tool/modelfit_subs.pm +++ b/lib/tool/modelfit_subs.pm @@ -769,6 +769,8 @@ start copy_model_and_input extra_data_file_names => \@new_extra_data_names, copy_data => 1 ); + $candidate_model -> shrinkage_modules( $model -> shrinkage_modules ); + $candidate_model -> register_in_database; $model -> flush_data; -- 2.11.4.GIT