lib/tool/sse_subs.pm

   1 # {{{ include statements
   2
   3 start include statements
   4 use tool::cdd;
   5 use tool::modelfit;
   6 use Math::Random;
   7 use Data::Dumper;
   8 end include statements
   9
  10 # }}} include statements
  11
  12 # {{{ new
  13
  14 start new
  15 foreach my $attribute ( 'logfile', 'raw_results_file' ) {
  16   if ( not( ref($this -> {$attribute}) eq 'ARRAY' or
  17             ref($this -> {$attribute}) eq 'HASH' ) ) {
  18     my $tmp = $this -> {$attribute};
  19     if ( not defined $tmp and $attribute eq 'logfile' ) {
  20         $tmp = 'mcs.log';
  21     }
  22     $this -> {$attribute} = [];
  23     for ( my $i = 1; $i <= scalar @{$this -> {'models'}}; $i++ ) {
  24       my $name = $tmp;
  25       if ( $name =~ /\./ ) {
  26         $name =~ s/\./$i\./;
  27       } else {
  28         $name = $name.$i;
  29       }
  30       my $ldir;
  31       ( $ldir, $name ) =
  32           OSspecific::absolute_path( $this -> {'directory'}, $name );
  33       push ( @{$this -> {$attribute}}, $ldir.$name ) ;
  34     }
  35   }
  36 }
  37
  38 if ( scalar (@{$this -> {'models'}->[0]-> problems}) != 1 ){
  39     'debug' -> die(message => 'The simulation model must contain exactly one problem.');
  40 }
  41
  42 foreach my $alt (@{$this -> {'alternative_models'}}){
  43     if ( scalar (@{$alt -> problems}) != 1 ){
  44         'debug' -> die(message => 'The alternative models must each contain exactly one problem.');
  45     }
  46 }
  47
  48 if ((scalar(@{$this -> {'alternative_models'}}) < 1) && (not $this->{'estimate_simulation'})){
  49   'debug' -> die(message => "Must estimate at least one model, ".
  50                  "either an alternative or the simulation model.");
  51 }
  52
  53 if ($this->{'estimate_simulation'} && (defined $this->{'ref_ofv'})){
  54   'debug'-> die(message => "Not allowed to use a reference ofv-value ".
  55                 "when estimating the simulation model.");
  56 }
  57 if ($this->{'samples'} < 1){
  58   'debug'-> die(message => "Must set -samples to at least 1.");
  59 }
  60
  61 end new
  62
  63 # }}}
  64
  65 # {{{ cdd_setup
  66
  67 start cdd_setup
  68
  69 my @subm_threads;
  70 if (ref( $self -> {'threads'} ) eq 'ARRAY') {
  71   @subm_threads = @{$self -> {'threads'}};
  72   shift(@subm_threads);
  73 } else {
  74   @subm_threads = ($self -> {'threads'});
  75 }
  76 $self -> general_setup( model_number => $model_number,
  77                         class        => 'tool::cdd',
  78                         subm_threads => \@subm_threads );
  79
  80 end cdd_setup
  81
  82 # }}} cdd_setup
  83
  84 # {{{ modelfit_setup
  85
  86 start modelfit_setup
  87
  88 my @subm_threads;
  89 if (ref( $self -> {'threads'} ) eq 'ARRAY') {
  90   @subm_threads = @{$self -> {'threads'}};
  91   shift(@subm_threads);
  92 } else {
  93   @subm_threads = ($self -> {'threads'});
  94 }
  95 $self -> general_setup( model_number => $model_number,
  96                         class        => 'tool::modelfit',
  97                         subm_threads => $subm_threads[$model_number-1] );
  98
  99 end modelfit_setup
 100
 101 # }}} modelfit_setup
 102
 103 # {{{ general_setup
 104
 105 start general_setup
 106 {
 107   my $model = $self -> {'models'} -> [$model_number-1];
 108   my @alternatives = @{$self -> {'alternative_models'}}; #may be empty
 109   my ( @seed, $new_datas, $skip_ids, $skip_keys, $skip_values );
 110   my ( @orig_est_models, @alt_est_models );
 111   my ( $sim_model, $est_alternative );
 112   my $alternative;
 113   my $time_in_input=0;
 114   my $datx_in_input=0;
 115
 116   my $done = ( -e $self -> {'directory'}."/m$model_number/done" ) ? 1 : 0;
 117
 118   if ( not $done ) {
 119
 120     # {{{ this is code for the first run , i.e. not a restart.
 121
 122     my @sim_models;
 123     my @table_header;
 124     my @all_simulated_files;
 125     my ( @orig_table_names, @alt_table_names );
 126
 127     for( my $sim_no = 1; $sim_no <= $self -> {'samples'}; $sim_no++ ) {
 128
 129       # {{{ Copy the model to new simulation models
 130
 131       my $est_original;
 132       my $sim_name = "mc-$sim_no.sim";
 133       my $orig_name = "mc-original-$sim_no.mod";
 134       my $out_name = "mc-$sim_no.lst";
 135       my $orig_out = "mc-original-$sim_no.lst";
 136       if( $sim_no == 1 ) {
 137         $sim_model = $model ->
 138             copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$sim_name,
 139                   target      => 'disk',
 140                   copy_data   => 1,
 141                   copy_output => 0);
 142         $sim_model -> drop_dropped;
 143
 144         if ($sim_model-> is_option_set(record=>'input',name=>'TIME')){
 145           #this assumes no synonym, and TIME is always option, not value.
 146           $time_in_input=1;
 147         }
 148         foreach my $col ('DATE','DAT1','DAT2','DAT3'){
 149           if ($sim_model-> is_option_set(record=>'input',name=>$col)){
 150             #this assumes no synonym, and name always options, not value.
 151             $datx_in_input=1;
 152             last;
 153           }
 154         }
 155
 156         if( $self -> {'estimate_simulation'} ) {
 157           $est_original = $model ->
 158               copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$orig_name,
 159                     target      => 'disk',
 160                     copy_data   => 0,
 161                     copy_output => 0);
 162           $est_original -> drop_dropped;
 163           $est_original -> remove_records( type => 'simulation' );
 164
 165           #ignore @ since simdata contains header rows. can skip old ignores since filtered
 166           $est_original -> set_option( record_name  => 'data',
 167                                        option_name  => 'IGNORE',
 168                                        option_value => '@');
 169
 170           #remove any DATX from est_original $INPUT
 171           if ($datx_in_input){
 172             foreach my $col ('DATE','DAT1','DAT2','DAT3'){
 173               $est_original -> remove_option(record_name => 'input',
 174                                              option_name => $col);
 175             }
 176             unless ($time_in_input){
 177               $est_original -> set_option(record_name => 'input',
 178                                           option_name => 'TIME');
 179             }
 180           }
 181         }
 182
 183         my $tbl_nm_ref = $model -> table_names();
 184         @orig_table_names = @{$tbl_nm_ref} if( defined $tbl_nm_ref );
 185
 186       } else {
 187         $sim_model = $sim_models[0] ->
 188             copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$sim_name,
 189                   target      => 'disk',
 190                   copy_data   => 0,
 191                   copy_output => 0);
 192
 193         if( $self -> {'estimate_simulation'} ){
 194           $est_original = $orig_est_models[0] ->
 195               copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$orig_name,
 196                     target      => 'disk',
 197                     copy_data   => 0,
 198                     copy_output => 0);
 199
 200         }
 201       }
 202
 203       $sim_model -> ignore_missing_files( 1 );
 204       $sim_model -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$out_name );
 205       $sim_model -> ignore_missing_files( 0 );
 206
 207       if( $self -> {'estimate_simulation'} ) {
 208         $est_original -> ignore_missing_files( 1 );
 209         $est_original -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$orig_out );
 210         $est_original -> ignore_missing_files( 0 );
 211
 212         my @new_orig_table_names;
 213
 214         for( my $pr = 0; $pr <= $#orig_table_names; $pr++ ) {
 215           for( my $tbl = 0; $tbl < scalar @{$orig_table_names[$pr]}; $tbl++ ) {
 216             $new_orig_table_names[$pr][$tbl] = $orig_table_names[$pr][$tbl].'-'.$sim_no;
 217           }
 218         }
 219
 220         if( $#new_orig_table_names >= 0 ) {
 221           $est_original -> table_names( new_names => \@new_orig_table_names );
 222         }
 223
 224       }
 225
 226       if( $self -> shrinkage() ) {
 227         $est_alternative -> shrinkage_stats( enabled => 1 );
 228
 229         if( $self -> {'estimate_simulation'} ){
 230           $est_original -> shrinkage_stats( enabled => 1 );
 231           my @problems = @{$est_original -> problems};
 232           for( my $i = 1; $i <= scalar @problems; $i++ ) {
 233             $problems[ $i-1 ] -> shrinkage_module -> model( $est_original );
 234           }
 235         }
 236       }
 237
 238       my @simulated_files;
 239       for( my $j = 0; $j < scalar @{$sim_model -> problems}; $j++ ) {
 240         my $prob = $sim_model -> problems -> [$j];
 241
 242         # set $SIMULATION record
 243
 244         my $sim_record = $sim_model -> record( problem_number => $j,
 245                                                record_name => 'simulation' );
 246
 247         if( scalar(@{$sim_record}) > 0 ){
 248
 249           my @new_record;
 250           foreach my $sim_line ( @{$sim_record -> [0]} ){
 251             my $new_line;
 252             while( $sim_line =~ /([^()]*)(\([^()]+\))(.*)/g ){
 253               my $head = $1;
 254               my $old_seed = $2;
 255               $sim_line = $3;
 256               $new_line .= $head;
 257
 258               while( $old_seed =~ /(\D*)(\d+)(.*)/ ){
 259                 $new_line .= $1;
 260                 $new_line .= random_uniform_integer( 1, 0, 1000000 ); # Upper limit is from nmhelp
 261                 $old_seed = $3;
 262               }
 263
 264               $new_line .= $old_seed;
 265
 266             }
 267
 268             push( @new_record, $new_line.$sim_line );
 269           }
 270
 271           $prob -> set_records( type => 'simulation',
 272                                 record_strings => \@new_record );
 273
 274         } else {
 275
 276           my $seed = random_uniform_integer( 1, 0, 1000000 ); # Upper limit is from nmhelp
 277           $prob -> set_records( type           => 'simulation',
 278                                 record_strings => [ '(' . $seed .
 279                                                     ') ONLYSIMULATION' ] );
 280         }
 281
 282         if( $sim_no == 1 ) {
 283           # remove $EST and $COV
 284
 285           if( $sim_model -> is_option_set( problem_number => $j,record => 'estimation',
 286                                            name => 'LIKELIHOOD',fuzzy_match => 1 )
 287               or
 288               $sim_model -> is_option_set( problem_number => $j,record => 'estimation',
 289                                            name => '-2LOGLIKELIHOOD',fuzzy_match => 1 )
 290               or
 291               $sim_model -> is_option_set( problem_number => $j, record => 'estimation',
 292                                            name => '-2LLIKELIHOOD',fuzzy_match => 1 )
 293               or
 294               $sim_model -> is_option_set( problem_number => $j, record => 'estimation',
 295                                            name => 'LAPLACIAN',fuzzy_match => 1 )
 296               ){
 297
 298             $sim_model -> set_option( problem_numbers => [$j], record_name => 'estimation',
 299                                       option_name => 'MAXEVALS',fuzzy_match => 1,
 300                                       option_value => 0 );
 301
 302           } else {
 303
 304             $prob -> remove_records(type => 'estimation');
 305
 306           }
 307
 308           $prob -> remove_records(type => 'covariance');
 309
 310           # set $TABLE record
 311           #when copying $INPUT to $TABLE: remove DATX
 312           #if not TIME present, remove TIME
 313           #if not TIME present but DATX then add back TIME at the end
 314
 315           if( defined $prob -> inputs and defined $prob -> inputs -> [0] -> options ) {
 316             foreach my $option ( @{$prob -> inputs -> [0] -> options} ) {
 317               push( @table_header, $option -> name ) unless
 318                   (($option -> name =~ /DAT(E|1|2|3)/) ||
 319                    ((not $time_in_input) && ($option -> name =~ /TIME/)));
 320             }
 321             if ((not $time_in_input) && ($datx_in_input )){
 322               push( @table_header, 'TIME');
 323             }
 324           } else {
 325             debug -> die( message => "Trying to construct table for monte-carlo simulation".
 326                           " but no headers were found in \$model_number-INPUT" );
 327           }
 328         }
 329         my $simulated_file = "mc-sim-".($j+1)."-$sim_no.dat";
 330         $prob -> set_records( type           => 'table',
 331                               record_strings => [ join( ' ', @table_header ).
 332                                                   ' NOPRINT NOAPPEND ONEHEADER FILE='.
 333                                                   $simulated_file ] );
 334         push( @simulated_files, $self -> {'directory'}.'m'.$model_number.'/'.
 335               $simulated_file );
 336       }
 337       $sim_model -> _write( write_data => $sim_no == 1 ? 1 : 0 );
 338       push( @sim_models, $sim_model );
 339       if( defined $est_original ){
 340         push( @orig_est_models, $est_original );
 341       }
 342       push( @all_simulated_files, \@simulated_files );
 343
 344       # }}}
 345
 346     }
 347
 348     my $threads = $self -> {'parallel_simulations'} ? $self -> {'parallel_simulations'} : 1;
 349
 350     my $mod_sim = tool::modelfit -> new( top_tool         => 0,
 351                                          reference_object => $self,
 352                                          models           => \@sim_models,
 353                                          base_directory   => $self -> {'directory'},
 354                                          directory        => $self -> {'directory'}.
 355                                          'simulation_dir'.$model_number,
 356                                          parent_tool_id   => $self -> {'tool_id'},
 357                                          retries          => 1,
 358                                          logfile                 => undef,
 359                                          raw_results           => undef,
 360                                          prepared_models       => undef,
 361                                          threads          => $threads );
 362
 363     $mod_sim -> run;
 364
 365     for( my $j = 0; $j <= $#orig_est_models; $j++ ) {
 366       $orig_est_models[$j] -> datafiles( new_names => $all_simulated_files[$j] );
 367       $orig_est_models[$j] -> _write;
 368       $orig_est_models[$j] -> flush_data();
 369     }
 370     @{$self -> {'mc_models'}} = ( @orig_est_models );
 371
 372     my $alternative_counter=0;
 373     foreach $alternative (@alternatives){
 374       $alternative_counter++;
 375       # {{{ create copies of the alternative models
 376       @alt_table_names = ();
 377       @alt_est_models =();
 378 #      my $filestem = 'mc-alt-' . $alternative -> filename();
 379 #      $filestem =~ s/\.[^.]*$//;
 380       my $filestem = 'mc-alternative_' . $alternative_counter;
 381       for( my $sim_no = 1; $sim_no <= $self -> {'samples'}; $sim_no++ ) {
 382         my $alt_name = "$filestem-$sim_no.mod";
 383         my $alt_out = "$filestem-$sim_no.lst";
 384
 385         if( $sim_no == 1 ) {
 386           $est_alternative = $alternative ->
 387               copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$alt_name,
 388                     target      => 'disk',
 389                     copy_data   => 0,
 390                     copy_output => 0);
 391           $est_alternative -> drop_dropped;
 392           #remove any DATX in $INPUT (drop_dropped does not)
 393           foreach my $col ('DATE','DAT1','DAT2','DAT3'){
 394             $est_alternative -> remove_option(record_name => 'input',
 395                                               option_name => $col);
 396           }
 397           #if added time then remove TIME (if present) and then add TIME (must be last in list)
 398           if ((not $time_in_input) && ($datx_in_input)){
 399             $est_alternative -> remove_option(record_name => 'input',
 400                                               option_name => 'TIME');
 401             $est_alternative -> set_option(record_name => 'input',
 402                                            option_name => 'TIME');
 403           }
 404
 405           $est_alternative -> remove_records( type => 'simulation' );
 406
 407           #ignore @ since simdata contains header rows. cannot skip other old ignores!!
 408
 409           my $ig_list = $est_alternative -> get_option_value( record_name  => 'data',
 410                                                               option_name  => 'IGNORE',
 411                                                               option_index => 'all');
 412           $est_alternative -> remove_option( record_name  => 'data',
 413                                              option_name  => 'IGNORE');
 414
 415           if (scalar (@{$ig_list})>0){
 416             foreach my $val (@{$ig_list}){
 417               unless ($val =~ /^.$/){
 418                 $est_alternative -> add_option( record_name  => 'data',
 419                                                 option_name  => 'IGNORE',
 420                                                 option_value => $val);
 421               }
 422             }
 423           }
 424           $est_alternative -> add_option( record_name  => 'data',
 425                                           option_name  => 'IGNORE',
 426                                           option_value => '@');
 427           ##done fixing ignore
 428
 429
 430
 431           my $tbl_nm_ref = $alternative -> table_names();
 432           @alt_table_names = @{$tbl_nm_ref} if( defined $tbl_nm_ref );
 433         } else {
 434           $est_alternative = $alt_est_models[0] ->
 435               copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$alt_name,
 436                     target      => 'disk',
 437                     copy_data   => 0,
 438                     copy_output => 0);
 439         }
 440         $est_alternative -> ignore_missing_files( 1 );
 441         $est_alternative -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$alt_out );
 442         $est_alternative -> ignore_missing_files( 0 );
 443         my @new_alt_table_names;
 444         for( my $pr = 0; $pr <= $#alt_table_names; $pr++ ) {
 445           for( my $tbl = 0; $tbl < scalar @{$alt_table_names[$pr]}; $tbl++ ) {
 446             $new_alt_table_names[$pr][$tbl] = $alt_table_names[$pr][$tbl].'-'.$sim_no;
 447           }
 448         }
 449         if( $#new_alt_table_names >= 0 ) {
 450           $est_alternative -> table_names( new_names => \@new_alt_table_names );
 451         }
 452         if( $self -> shrinkage() ) {
 453           my @problems = @{$est_alternative -> problems};
 454           for( my $i = 1; $i <= scalar @problems; $i++ ) {
 455             $problems[ $i-1 ] -> shrinkage_module -> model( $est_alternative );
 456           }
 457         }
 458         push( @alt_est_models, $est_alternative );
 459       } #end loop over sim_no
 460
 461       for( my $j = 0; $j <= $#alt_est_models; $j++ ) {
 462         $alt_est_models[$j] -> datafiles( new_names => $all_simulated_files[$j] );
 463         $alt_est_models[$j] -> _write;
 464         $alt_est_models[$j] -> flush_data();
 465       }
 466       push (@{$self -> {'mc_models'}}, @alt_est_models);
 467
 468       # }}}
 469
 470     } #end loop over alternatives
 471
 472     # Create a checkpoint.
 473     open( DONE, ">".$self -> {'directory'}."/m$model_number/done" ) ;
 474     print DONE "Simulation from ",$sim_models[0] -> filename," through ",
 475     $sim_models[$#sim_models] -> filename," performed\n";
 476     print DONE $self -> {'samples'}." samples\n";
 477     @seed = random_get_seed();
 478     print DONE "seed: @seed\n";
 479     close( DONE );
 480
 481     # }}}
 482
 483   } else {
 484
 485     # {{{ code for restarts
 486
 487     #$done=true
 488     # Recreate the datasets and models from a checkpoint
 489     ui -> print( category => 'mc',
 490                  message  => "Recreating models from a previous run" );
 491     open( DONE, $self -> {'directory'}."/m$model_number/done" );
 492     my @rows = <DONE>;
 493     close( DONE );
 494     my ( $junk, $junk, $stored_filename1, $junk,
 495          $stored_filename2, $junk ) = split(' ',$rows[0],4);
 496     my ( $stored_samples, $junk ) = split(' ',$rows[1],2);
 497     @seed = split(' ',$rows[2]);
 498     shift( @seed ); # get rid of 'seed'-word
 499     # Reinitiate the model objects
 500     for ( my $j = 1; $j <= $stored_samples; $j++ ) {
 501       my ($model_dir, $orig_name) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 502                                                                $model_number,
 503                                                                "mc-original-$j.mod" );
 504       my ($out_dir, $orig_out) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 505                                                             $model_number,
 506                                                             "mc-original-$j.lst" );
 507       my $est_original = model ->
 508           new( directory   => $model_dir,
 509                filename    => $orig_name,
 510                outputfile  => $orig_out,
 511                extra_files => $model -> extra_files,
 512                target      => 'disk',
 513                ignore_missing_files => 1,
 514                extra_output => $model -> extra_output(),
 515                cwres       => $model -> cwres(),
 516                quick_reload => 1 );
 517       push( @orig_est_models, $est_original );
 518       my $nl = $j == $stored_samples ? "" : "\r";
 519       ui -> print( category => 'mc',
 520                    message  => ui -> status_bar( sofar => $j+1,
 521                                                  goal  => $stored_samples+1 ).$nl,
 522                    wrap     => 0,
 523                    newline  => 0 );
 524     } #end loop over samples
 525     @{$self -> {'mc_models'}} = ( @orig_est_models);
 526
 527     #start alternatives section
 528     my $alternative_counter=0;
 529     foreach $alternative (@alternatives){
 530       $alternative_counter++;
 531 #      my $filestem = 'mc-alt-' . $alternative -> filename();
 532 #      $filestem =~ s/\.[^.]*$//;
 533       my $filestem = 'mc-alternative_' . $alternative_counter;
 534       @alt_est_models =();
 535       for ( my $j = 1; $j <= $stored_samples; $j++ ) {
 536         my ($model_dir, $alt_name) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 537                                                                 $model_number,
 538                                                                 "$filestem-$j.mod" );
 539         my ($out_dir, $alt_out) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 540                                                              $model_number,
 541                                                              "$filestem-$j.lst" );
 542         my $est_alternative = model ->
 543             new( directory   => $model_dir,
 544                  filename    => $alt_name,
 545                  outputfile  => $alt_out,
 546                  extra_files => $alternative -> extra_files,
 547                  target      => 'disk',
 548                  ignore_missing_files => 1,
 549                  extra_output => $model -> extra_output(),
 550                  cwres       => $model -> cwres(),
 551                  quick_reload => 1 );
 552         push( @alt_est_models, $est_alternative );
 553         my $nl = $j == $stored_samples ? "" : "\r";
 554         ui -> print( category => 'mc',
 555                      message  => ui -> status_bar( sofar => $j+1,
 556                                                    goal  => $stored_samples+1 ).$nl,
 557                      wrap     => 0,
 558                      newline  => 0 );
 559       } #end loop over samples
 560       push (@{$self -> {'mc_models'}}, @alt_est_models );
 561
 562     } #end loop over alternatives
 563
 564     ui -> print( category => 'mc', message  => " ... done." );
 565     random_set_seed( @seed );
 566     ui -> print( category => 'mc',
 567                  message  => "Using $stored_samples data sets, previously simulated ".
 568                  "from $stored_filename1 through $stored_filename2" )
 569         unless $self -> {'parent_threads'} > 1;
 570
 571     # }}}
 572
 573   } #end $done
 574
 575   my $subdir = $class;
 576   $subdir =~ s/tool:://;
 577   my @subtools = @{$self -> {'subtools'}};
 578   shift( @subtools );
 579   my %subargs = ();
 580   if ( defined $self -> {'subtool_arguments'} ) {
 581     %subargs = %{$self -> {'subtool_arguments'}};
 582   }
 583
 584   push( @{$self -> {'tools'}},
 585         $class -> new
 586         ( reference_object => $self,
 587           top_tool         => 0,
 588           logfile                => undef,
 589           raw_results      => undef,
 590           prepared_models  => undef,
 591           models         => $self -> {'mc_models'},
 592           nm_version     => $self -> {'nm_version'},
 593           parent_tool_id   => $self -> {'tool_id'},
 594           base_directory => $self -> {'directory'},
 595           directory      => $self -> {'directory'}.'/'.$subdir.'_dir'.$model_number,
 596           subtools       => $#subtools >= 0 ? \@subtools : undef,
 597           shrinkage      => $self -> {'shrinkage'},
 598           _raw_results_callback => $self ->
 599           _modelfit_raw_results_callback( model_number => $model_number ),
 600           %subargs ) );
 601 }
 602 end general_setup
 603
 604 # }}} general_setup
 605
 606 # {{{ cdd_analyze
 607 start cdd_analyze
 608     $self -> {'tools'}[0] -> print_results;
 609 end cdd_analyze
 610 # }}} cdd_analyze
 611
 612 # {{{ modelfit_analyze
 613
 614 start modelfit_analyze
 615 $self -> {'tools'}[0] -> print_results;
 616 end modelfit_analyze
 617
 618 # }}}
 619
 620 # {{{ _modelfit_raw_results_callback
 621
 622 start _modelfit_raw_results_callback
 623
 624 # Use the mc's raw_results file.
 625 my ($dir,$file) =
 626     OSspecific::absolute_path( $self -> {'directory'},
 627                                $self -> {'raw_results_file'}[$model_number-1] );
 628 my $orig_mod = $self -> {'models'}[$model_number-1];
 629 $subroutine = sub {
 630   my $modelfit = shift;
 631   my $mh_ref   = shift;
 632   my %max_hash = %{$mh_ref};
 633   $modelfit -> raw_results_file( $dir.$file );
 634
 635   # {{{ New header
 636
 637   if ( defined $modelfit -> raw_results() ) {
 638     my @rows = @{$modelfit -> raw_results()};
 639
 640     #the first 'samples' (a scalar) number of rows is null/simulation models
 641     if( $self -> {'estimate_simulation'} ){
 642       for( my $i = 0; $i < $self -> {'samples'}; $i++ ) {
 643         unshift( @{$rows[$i]}, 'simulation' );
 644       }
 645     }
 646     #remaining rows, if any, belong to alternatives model(s)
 647     my $offset;
 648     if( $self -> {'estimate_simulation'} ){
 649       $offset = $self -> {'samples'};
 650     } else {
 651       $offset = 0;
 652     }
 653     my $alternative_counter=0;
 654     foreach my $alternative (@{$self -> {'alternative_models'}}){
 655       $alternative_counter++;
 656       if (($offset + $self -> {'samples'}-1)> $#rows){
 657         debug -> die( message => "Too few rows in raw results file." );
 658       }
 659       my $altname = 'mc-alternative_'.$alternative_counter;
 660       for( my $i = 0; $i < $self -> {'samples'}; $i++ ) {
 661         unshift( @{$rows[($offset+$i)]}, $altname );
 662       }
 663       $offset += $self -> {'samples'};
 664     }
 665     if (($offset -1)< $#rows){
 666       debug -> die( message => "Too many rows in raw results file." );
 667     }
 668
 669     unshift( @{$modelfit -> {'raw_results_header'}}, 'hypothesis' );
 670
 671   }
 672   @{$self -> {'raw_results_header'}} = @{$modelfit -> {'raw_results_header'}};
 673   # }}} New header
 674
 675 };
 676 return $subroutine;
 677
 678 end _modelfit_raw_results_callback
 679
 680 # }}} _modelfit_raw_results_callback
 681
 682 # {{{ prepare_results
 683
 684 start prepare_results
 685 {
 686
 687   $self -> {'raw_results'} = $self -> {'tools'}[0] -> {'raw_results'};
 688
 689   my $model = $self -> {'models'} -> [0];
 690   my @alternatives = @{$self -> {'alternative_models'}};
 691
 692   my $samples = $self -> {'samples'};
 693
 694   my $raw_line_struct = $self -> {'tools'}[0] -> {'raw_line_structure'};
 695
 696   my %orig_results_section;
 697   my @alt_results_sections;
 698
 699   ### Retrieve initial parameter estimates from original model.
 700
 701   my %initial_values;
 702   my %n_initials;
 703   my (@thetalist,@omegalist,@sigmalist);
 704   my $tmp;
 705
 706   $n_initials{'theta'} = $model->nthetas(problem_number => 1);
 707   $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'theta');
 708   @thetalist=@{$tmp->[0]};
 709   'debug'->die(message=>'Error number of initial thetas.')
 710       unless (scalar(@thetalist) == $n_initials{'theta'});
 711   $initial_values{'theta'}=\@thetalist;
 712
 713   $tmp = $model->nomegas(problem_numbers => [1],
 714                          with_correlations => 1);
 715   $n_initials{'omega'} = $tmp->[0];
 716   $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'omega');
 717   @omegalist=@{$tmp->[0]};
 718   'debug'->die(message=>'Error number of initial omegas.')
 719       unless (scalar(@omegalist) == $n_initials{'omega'});
 720   $initial_values{'omega'}=\@omegalist;
 721
 722   $tmp = $model->nsigmas(problem_numbers => [1],
 723                          with_correlations => 1);
 724   $n_initials{'sigma'} = $tmp->[0];
 725   $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'sigma');
 726   @sigmalist=@{$tmp->[0]};
 727   'debug'->die(message=>'Error number of initial sigmas.')
 728       unless (scalar(@sigmalist) == $n_initials{'sigma'});
 729   $initial_values{'sigma'}=\@sigmalist;
 730
 731   $n_initials{'ofv'}=0;
 732
 733   ## Prepare general run info for output file
 734   my %return_section;
 735   $return_section{'name'} = 'SSE run info';
 736   my $modelname=$self-> {'models'}->[0] ->filename();
 737   $return_section{'labels'} = [[],['Date','samples','simulation model',
 738                                    'PsN version','NONMEM version']];
 739   my @datearr=localtime;
 740   my $the_date=($datearr[5]+1900).'-'.($datearr[4]+1).'-'.($datearr[3]);
 741   $return_section{'values'} = [[$the_date,$self->{'samples'},$modelname,
 742                                 $PsN::version,$self->{'nm_version'}]];
 743   push( @{$self -> {'results'}[0]{'own'}},\%return_section );
 744
 745
 746   ### Prepare result section with ofv reference value , if used, and true parameter values for sim
 747   my %reference_section;
 748   $reference_section{'name'} = 'Reference/True values';
 749   my @reference_headers=('ofv');
 750   my @reference_values;
 751   if (defined $self->{'ref_ofv'}){
 752     @reference_values=($self->{'ref_ofv'});
 753   }elsif ($self->{'estimate_simulation'}){
 754     @reference_values=('Sim. model result');
 755   } else{
 756     @reference_values=('Alt. model 1 result');
 757   }
 758   foreach my $measure ( 'theta','omega','sigma' ){
 759     for (my $i=1; $i<=$n_initials{$measure};$i++){
 760       push(@reference_headers, uc(substr($measure, 0,2))."_$i" );
 761       push(@reference_values,$initial_values{$measure}->[$i-1]);
 762     }
 763   }
 764
 765   $reference_section{'labels'}=[[],\@reference_headers];
 766   $reference_section{'values'}=[\@reference_values];
 767   push( @{$self -> {'results'}[0]{'own'}},\%reference_section );
 768
 769   ### Loop over the original model (the zero) and over all alternative
 770   ### models.
 771
 772   my $end_index = ( $self -> {'estimate_simulation'})? scalar(@alternatives): scalar(@alternatives) - 1;
 773
 774   for my $model_index ( 0..$end_index ){
 775
 776     ### Foreach such model, create a %results_section
 777
 778     my %results_section;
 779
 780     ### Each results_section contains statistics of ofv and
 781     ### parameters.
 782
 783     my @values;
 784     my @labels;
 785
 786     @{$labels[0]} = ('mean','median','sd','max','min','skewness',
 787                      'kurtosis','rmse','bias','rse');
 788     my @ci_list=(0.5,2.5,5,95,97.5,99.5);
 789     my @z_list=(-2.58,-1.96,-1.645,1.645,1.96,2.58);
 790     push(@{$labels[0]},('standard error CI'));
 791     foreach my $ci (@ci_list){
 792       push(@{$labels[0]},"$ci\%");
 793     }
 794
 795     if( $model_index == 0 and $self -> {'estimate_simulation'}){
 796       $results_section{'name'} = "\nSimulation model";
 797     } else {
 798       my $alt_index;
 799       if( $self -> {'estimate_simulation'} ){
 800         $alt_index = $model_index;
 801       } else {
 802         $alt_index = $model_index + 1;
 803       }
 804       $results_section{'name'} = "\nAlternative model ".$alt_index;
 805     }
 806
 807     $results_section{'labels'} = \@labels;
 808     $results_section{'values'} = \@values;
 809
 810     push( @{$self -> {'results'}[0]{'own'}},\%results_section );
 811
 812     MEASURE: foreach my $measure ( 'ofv','theta','omega','sigma' ){
 813
 814       # Create a header, we might use max_hash from modelfit for this
 815       # later on.
 816
 817       my $start=0;
 818       my $length=0;
 819
 820       for( 0..($samples-1) ){
 821         my ($s, $l) = split(/,/, $raw_line_struct -> {1+$model_index*$samples + $_} -> {$measure});
 822
 823         if( defined $s and $start != 0 and $s != $start ){
 824           'debug' -> die( message => "Internal structure corrupted at $model_index, ".
 825                           "$samples, $_, $measure, this is a bug" );
 826         } else {
 827           $start = $s;
 828         }
 829
 830         if( defined $l and $l > $length ){
 831           $length = $l;
 832         }
 833         if( $start == 0 ){
 834           for( 0..(10+scalar(@z_list)) ){
 835             # At least two reasons for this. Either the "measure" does
 836             # not exist or the model has not terminated.
 837             push(@{$values[$_]},'NA');
 838           }
 839           push( @{$labels[1]},$measure );
 840           next MEASURE;
 841         }
 842
 843       }
 844
 845
 846       if( $measure eq 'ofv' ){
 847         push( @{$labels[1]}, $measure );
 848       } else {
 849         for( 1..$length ){
 850           push( @{$labels[1]}, uc(substr($measure, 0,2))."_$_" );
 851         }
 852       }
 853
 854       my $initial_index=0;
 855       foreach my $col ( $start..($start + $length-1) ){
 856         my ($skew, $kurt, $mean, $stdev)
 857             = $self -> skewness_and_kurtosis( column_index => $col+1,
 858                                               start_row_index => $model_index*$samples,
 859                                               end_row_index => $model_index*$samples+$samples-1 );
 860         my ($maximum,$minimum)
 861             = $self -> max_and_min( column_index => $col+1,
 862                                     start_row_index => $model_index*$samples,
 863                                     end_row_index => $model_index*$samples+$samples-1 );
 864
 865
 866         push( @{$values[0]}, $mean );
 867         push( @{$values[2]}, $stdev );
 868         push( @{$values[3]}, $maximum );
 869         push( @{$values[4]}, $minimum );
 870         push( @{$values[5]}, $skew );
 871         push( @{$values[6]}, $kurt );
 872
 873         my $median = $self -> median( column_index => $col+1,
 874                                       start_row_index => $model_index*$samples,
 875                                       end_row_index => $model_index*$samples+$samples-1 );
 876
 877         push( @{$values[1]}, $median );
 878
 879         if( $initial_index < $n_initials{$measure}){
 880
 881           my $initial_value= $initial_values{$measure}->[$initial_index];
 882
 883           $initial_index++;
 884
 885           my $rmse = ' ';
 886           my $bias = ' ';
 887           my $rse = ' ';
 888
 889           if( defined $initial_value ){
 890             $rmse = $self -> rmse_percent( column_index => $col+1,
 891                                            start_row_index => $model_index*$samples,
 892                                            end_row_index => $model_index*$samples+$samples-1,
 893                                            initial_value => $initial_value );
 894
 895             $bias = $self -> bias_percent( column_index => $col+1,
 896                                            start_row_index => $model_index*$samples,
 897                                            end_row_index => $model_index*$samples+$samples-1,
 898                                            initial_value => $initial_value );
 899
 900             $rse = ($initial_value != 0)? 100*$stdev/$initial_value : 'NA';
 901           }
 902
 903           push( @{$values[7]}, $rmse );
 904           push( @{$values[8]}, $bias );
 905           push( @{$values[9]}, $rse );
 906
 907           my $label_index=11; #leave 10 empty
 908           foreach my $zval (@z_list){
 909             if( defined $initial_value ){
 910               push( @{$values[$label_index]},($bias+$zval*$rse/sqrt($samples)));
 911             } else {
 912               push( @{$values[$label_index]}, ' ' );
 913             }
 914             $label_index++;
 915           }
 916
 917         } else {
 918           #push placeholders so that values will stay in correct column
 919           push( @{$values[7]}, ' ' );
 920           push( @{$values[8]}, ' ' );
 921           push( @{$values[9]}, ' ' );
 922           my $label_index=11; #leave 10 empty
 923           foreach my $zval (@z_list){
 924             push( @{$values[$label_index]}, ' ' );
 925             $label_index++;
 926           }
 927         }
 928       }
 929     }
 930   }
 931
 932   ## Calculate OFV mean calculations
 933   #if at least two estimated models or ref_ofv given
 934   if( $end_index > 0 or defined $self -> {'ref_ofv'} ){
 935
 936     my @ofv_limits = (0,3.84,5.99,7.81,9.49);
 937
 938     my (%results_section_below, %results_section_above);
 939     my (@labels_below, @labels_above, @values_below, @values_above);
 940
 941     @{$labels_above[1]} = ('mean delta_OFV');
 942     @{$labels_below[1]} = ('mean delta_OFV');
 943
 944     foreach my $limit( @ofv_limits ){
 945       push( @{$labels_above[1]}, "$limit" );
 946       push( @{$labels_below[1]}, "$limit" );
 947     }
 948
 949     $results_section_above{'name'} = "\nOFV Statistics\nType I error rate\n,,Percent delta_OFV > N";
 950     $results_section_above{'labels'} = \@labels_above;
 951     $results_section_above{'values'} = \@values_above;
 952
 953     $results_section_below{'name'} = "\n1 - type II error rate (power)\n,,Percent delta_OFV < N";
 954     $results_section_below{'labels'} = \@labels_below;
 955     $results_section_below{'values'} = \@values_below;
 956
 957
 958     push( @{$self -> {'results'}[0]{'own'}},\%results_section_above );
 959     push( @{$self -> {'results'}[0]{'own'}},\%results_section_below );
 960
 961
 962     my $reference_string;
 963     my $alt_index;
 964     if ( defined $self -> {'ref_ofv'} ){
 965       $reference_string='reference ofv -';
 966       $alt_index = 1;
 967     } elsif( $self -> {'estimate_simulation'} ){
 968       $reference_string= 'simulation -';
 969       $alt_index=1;
 970     } else {
 971       $reference_string= 'alternative1 -';
 972       $alt_index = 2;
 973     }
 974
 975     # $end_index is already set.
 976     my $start_index = (defined $self -> {'ref_ofv'})? 0 : 1;
 977
 978     ALTMODEL: foreach my $model_index( $start_index..$end_index ){
 979
 980       push( @{$labels_above[0]}, "$reference_string alternative".($alt_index));
 981       push( @{$labels_below[0]}, "$reference_string alternative".($alt_index));
 982       $alt_index++;
 983
 984       my $start=0;
 985       my $delta_sum=0;
 986       my %nr_below;
 987       my %nr_above;
 988       my $nr_terminated=0;
 989       my $org_ofv;
 990
 991       foreach my $sample( 0..($samples-1) ){
 992
 993         my ($s, $l) = split(/,/, $raw_line_struct -> {1+$model_index*$samples + $sample} -> {'ofv'});
 994         my @dummy= split(/,/, $raw_line_struct -> {$model_index*$samples + $sample} -> {'ofv'});
 995         if( $s ne '' and $start != 0 and $s != $start ){
 996           'debug' -> die( message => "OFV: Internal structure corrupted at ".
 997                           "$model_index, $samples, $sample, this is a bug" );
 998         } else {
 999           $start = $s;
1000         }
1001         if($start eq '' ){
1002           # This model probably did not terminate.
1003           push(@{$values_below[$model_index-$start_index]},"Alternative model not terminated for all samples");
1004           push(@{$values_above[$model_index-$start_index]},"Alternative model not terminated for all samples");
1005           next ALTMODEL;
1006         }
1007
1008         if ( defined $self -> {'ref_ofv' } ){
1009           $org_ofv = $self -> {'ref_ofv'};
1010         } else {
1011           #check reference model terminated
1012           my ($s, $l) = split(/,/, $raw_line_struct -> {1+ $sample} -> {'ofv'});
1013           if( $s ne ''){
1014             $org_ofv = $self -> {'raw_results'} -> [$sample][$start+1];
1015           }else {
1016             next; #sample
1017           }
1018         }
1019
1020         my $alt_ofv = $self -> {'raw_results'} -> [$model_index*$samples + $sample][$start+1];
1021         $nr_terminated++;
1022         my $delta = $org_ofv - $alt_ofv;
1023         $delta_sum += $delta;
1024
1025         foreach my $limit ( @ofv_limits ){
1026
1027           if( $delta < -$limit ){
1028             $nr_below{$limit}++;
1029           }
1030
1031           if( $delta > $limit ) {
1032             $nr_above{$limit}++;
1033           }
1034
1035         }
1036
1037       } #end loop over samples
1038
1039       if ($nr_terminated > 0){
1040         push( @{$values_below[$model_index-$start_index]}, ($delta_sum / $nr_terminated) );
1041         foreach my $limit( @ofv_limits ){
1042           push( @{$values_below[$model_index-$start_index]}, ($nr_below{$limit} / $nr_terminated)*100 . "%" );
1043         }
1044
1045         push( @{$values_above[$model_index-$start_index]}, ($delta_sum / $nr_terminated) );
1046         foreach my $limit( @ofv_limits ){
1047           push( @{$values_above[$model_index-$start_index]}, ($nr_above{$limit} / $nr_terminated)*100 . "%" );
1048         }
1049       } else {
1050         push( @{$values_below[$model_index-$start_index]},"Reference model not terminated for any sample" );
1051         push( @{$values_above[$model_index-$start_index]},"Reference model not terminated for any sample");
1052       }
1053     }
1054   }
1055
1056 }
1057 end prepare_results
1058
1059 # }}}
1060
1061 # {{{ rmse_percent
1062
1063 start rmse_percent
1064 {
1065   #input is integers $column_index, $start_row_index, $end_row_index and scalar float $initial_value
1066   #output is scalar $rmse_percent
1067
1068   unless( $end_row_index ){
1069     $end_row_index = $#{$self -> {'raw_results'}};
1070   }
1071
1072   debug -> die( message => "Bad row index input") if ($start_row_index > $end_row_index);
1073   my $row_count = $end_row_index - $start_row_index +1;
1074   my $sum_squared_errors=0;
1075   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
1076     $sum_squared_errors += ($self->{'raw_results'}->[$i][$column_index] - $initial_value)**2;
1077   }
1078
1079   if ($initial_value == 0){
1080     debug -> warn( level => 2,
1081                    message => "Initial value 0, returning absolute rmse instead of relative.");
1082     $rmse_percent= sqrt($sum_squared_errors/$row_count);
1083   }else{
1084     $rmse_percent= (sqrt($sum_squared_errors/$row_count))*100/abs($initial_value);
1085   }
1086 }
1087 end rmse_percent
1088
1089 # }}} rmse_percent
1090
1091 # {{{ bias_percent
1092
1093 start bias_percent
1094 {
1095   #input is integers $column_index, $start_row_index, $end_row_index and scalar float $initial_value
1096   #output is scalar $bias_percent
1097
1098   unless( $end_row_index ){
1099     $end_row_index = $#{$self -> {'raw_results'}};
1100   }
1101
1102   debug -> die( message => "Bad row index input") if ($start_row_index > $end_row_index);
1103
1104   my $row_count = $end_row_index - $start_row_index +1;
1105   my $sum_errors=0;
1106   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
1107     $sum_errors += ($self->{'raw_results'}->[$i][$column_index] - $initial_value);
1108   }
1109
1110   if ($initial_value == 0){
1111     debug -> warn( level => 2,
1112                    message => "Initial value 0, returning absolute bias instead of relative.");
1113     $bias_percent= ($sum_errors/$row_count);
1114   }else{
1115     $bias_percent= ($sum_errors/$row_count)*100/abs($initial_value);
1116   }
1117 }
1118 end bias_percent
1119
1120 # }}} bias_percent
1121
1122 # {{{ skewness_and_kurtosis
1123
1124 start skewness_and_kurtosis
1125 {
1126   #input is integers $column_index, $start_row_index, $end_row_index
1127
1128   unless( $end_row_index ){
1129     $end_row_index = $#{$self -> {'raw_results'}};
1130   }
1131
1132   debug -> die( message => "Bad row index input") if ($start_row_index >= $end_row_index);
1133
1134   my $row_count = $end_row_index - $start_row_index +1;
1135   my $sum_values=0;
1136   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
1137     $sum_values += $self->{'raw_results'}->[$i][$column_index];
1138   }
1139
1140   $mean=$sum_values/$row_count;
1141
1142   my $error=0;
1143   my $sum_errors_pow2=0;
1144   my $sum_errors_pow3=0;
1145   my $sum_errors_pow4=0;
1146
1147   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
1148     $error = ($self->{'raw_results'}->[$i][$column_index] - $mean);
1149     $sum_errors_pow2 += $error**2;
1150     $sum_errors_pow3 += $error**3;
1151     $sum_errors_pow4 += $error**4;
1152   }
1153
1154   ## TODO fråga om missing values. och om SD
1155
1156   $stdev=0;
1157   unless( $sum_errors_pow2 == 0 ){
1158
1159     $stdev= sqrt ($sum_errors_pow2/($row_count-1));
1160     $skewness = $sum_errors_pow3/($row_count*($stdev**3));
1161     $kurtosis = -3 + $sum_errors_pow4/($row_count*($stdev**4));
1162
1163   }
1164 }
1165 end skewness_and_kurtosis
1166
1167 # }}} skewness_and_kurtosis
1168
1169 # {{{ max_and_min
1170
1171 start max_and_min
1172 {
1173   #input is integers $column_index, $start_row_index, $end_row_index
1174
1175   unless( $end_row_index ){
1176     $end_row_index = $#{$self -> {'raw_results'}};
1177   }
1178
1179   debug -> die( message => "Bad row index input") if ($start_row_index >= $end_row_index);
1180
1181   $maximum=$self->{'raw_results'}->[$start_row_index][$column_index];
1182   $minimum=$self->{'raw_results'}->[$start_row_index][$column_index];
1183   for (my $i=$start_row_index+1; $i<=$end_row_index; $i++){
1184     $maximum = $self->{'raw_results'}->[$i][$column_index] if
1185         ($self->{'raw_results'}->[$i][$column_index] > $maximum);
1186     $minimum = $self->{'raw_results'}->[$i][$column_index] if
1187         ($self->{'raw_results'}->[$i][$column_index] < $minimum);
1188   }
1189
1190 }
1191 end max_and_min
1192
1193 # }}} max_and_min
1194
1195
1196 # {{{ median
1197
1198 start median
1199 {
1200   #input is integers $column_index, $start_row_index, $end_row_index
1201
1202   unless( $end_row_index ){
1203     $end_row_index = $#{$self -> {'raw_results'}};
1204   }
1205
1206   debug -> die( message => "Bad row index input") if ($start_row_index >= $end_row_index);
1207
1208   my @temp;
1209
1210   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
1211     push( @temp, $self->{'raw_results'}->[$i][$column_index] );
1212   }
1213
1214   @temp = sort({$a <=> $b} @temp);
1215   if( scalar( @temp ) % 2 ){
1216     $median = $temp[$#temp/2];
1217   } else {
1218     $median = ($temp[@temp/2]+$temp[(@temp-2)/2]) / 2;
1219   }
1220
1221 }
1222 end median
1223
1224 # }}} median