lib/tool/mc_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 original 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 end new
  49
  50 # }}}
  51
  52 # {{{ cdd_setup
  53
  54 start cdd_setup
  55
  56 my @subm_threads;
  57 if (ref( $self -> {'threads'} ) eq 'ARRAY') {
  58   @subm_threads = @{$self -> {'threads'}};
  59   shift(@subm_threads);
  60 } else {
  61   @subm_threads = ($self -> {'threads'});
  62 }
  63 $self -> general_setup( model_number => $model_number,
  64                         class        => 'tool::cdd',
  65                         subm_threads => \@subm_threads );
  66
  67 end cdd_setup
  68
  69 # }}} cdd_setup
  70
  71 # {{{ modelfit_setup
  72
  73 start modelfit_setup
  74
  75 my @subm_threads;
  76 if (ref( $self -> {'threads'} ) eq 'ARRAY') {
  77   @subm_threads = @{$self -> {'threads'}};
  78   shift(@subm_threads);
  79 } else {
  80   @subm_threads = ($self -> {'threads'});
  81 }
  82 $self -> general_setup( model_number => $model_number,
  83                         class        => 'tool::modelfit',
  84                         subm_threads => $subm_threads[$model_number-1] );
  85
  86 end modelfit_setup
  87
  88 # }}} modelfit_setup
  89
  90 # {{{ general_setup
  91
  92 start general_setup
  93
  94 my $model = $self       -> {'models'} -> [$model_number-1];
  95 #my $alternative = $self -> {'alternative_models'} -> [$model_number-1];
  96 my @alternatives = @{$self -> {'alternative_models'}}; #may be empty
  97 my $subm_threads = $parm{'subm_threads'};
  98 my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
  99     $self -> {'threads'} -> [0]:$self -> {'threads'};
 100 my $done = ( -e $self -> {'directory'}."/m$model_number/done" ) ? 1 : 0;
 101 my ( @seed, $new_datas, $skip_ids, $skip_keys, $skip_values );
 102 my ( @orig_est_models, @alt_est_models );
 103 my ( $sim_model, $est_original, $est_alternative );
 104 my $alternative;
 105
 106 if ( not $done ) {
 107   my @sim_models;
 108   my @table_header;
 109   my @all_simulated_files;
 110   my ( @orig_table_names, @alt_table_names );
 111
 112 # Copy the model to a new simulation models
 113   for( my $sim_no = 1; $sim_no <= $self -> {'samples'}; $sim_no++ ) {
 114     my $sim_name = "mc-$sim_no.sim";
 115     my $orig_name = "mc-original-$sim_no.mod";
 116     my $out_name = "mc-$sim_no.lst";
 117     my $orig_out = "mc-original-$sim_no.lst";
 118     if( $sim_no == 1 ) {
 119       $sim_model = $model ->
 120           copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$sim_name,
 121                 target      => 'disk',
 122                 copy_data   => 1,
 123                 copy_output => 0);
 124       $sim_model -> drop_dropped;
 125       $est_original = $model ->
 126           copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$orig_name,
 127                 target      => 'disk',
 128                 copy_data   => 0,
 129                 copy_output => 0);
 130       $est_original -> drop_dropped;
 131       $est_original -> remove_records( type => 'simulation' );
 132
 133
 134       #ignore @ since simdata contains header rows. can skip old ignores since filtered
 135       $est_original -> set_option( record_name  => 'data',
 136                                    option_name  => 'IGNORE',
 137                                    option_value => '@');
 138
 139       my $tbl_nm_ref = $model -> table_names();
 140       @orig_table_names = @{$tbl_nm_ref} if( defined $tbl_nm_ref );
 141       } else {
 142       $sim_model = $sim_models[0] ->
 143           copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$sim_name,
 144                 target      => 'disk',
 145                 copy_data   => 0,
 146                 copy_output => 0);
 147       $est_original = $orig_est_models[0] ->
 148           copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$orig_name,
 149                 target      => 'disk',
 150                 copy_data   => 0,
 151                 copy_output => 0);
 152     }
 153     $sim_model -> ignore_missing_files( 1 );
 154     $sim_model -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$out_name );
 155     $sim_model -> ignore_missing_files( 0 );
 156     $est_original -> ignore_missing_files( 1 );
 157     $est_original -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$orig_out );
 158     $est_original -> ignore_missing_files( 0 );
 159     my @new_orig_table_names;
 160     for( my $pr = 0; $pr <= $#orig_table_names; $pr++ ) {
 161       for( my $tbl = 0; $tbl < scalar @{$orig_table_names[$pr]}; $tbl++ ) {
 162         $new_orig_table_names[$pr][$tbl] = $orig_table_names[$pr][$tbl].'-'.$sim_no;
 163       }
 164     }
 165     if( $#new_orig_table_names >= 0 ) {
 166       $est_original -> table_names( new_names => \@new_orig_table_names );
 167     }
 168     if( $self -> shrinkage() ) {
 169       $est_original -> shrinkage_stats( enabled => 1 );
 170       $est_alternative -> shrinkage_stats( enabled => 1 );
 171       my @problems = @{$est_original -> problems};
 172       for( my $i = 1; $i <= scalar @problems; $i++ ) {
 173         $problems[ $i-1 ] -> shrinkage_module -> model( $est_original );
 174       }
 175     }
 176
 177     my @simulated_files;
 178     for( my $j = 0; $j < scalar @{$sim_model -> problems}; $j++ ) {
 179       my $prob = $sim_model -> problems -> [$j];
 180
 181       # set $SIMULATION record
 182
 183       my $sim_record = $sim_model -> record( problem_number => $j,
 184                                              record_name => 'simulation' );
 185
 186       if( scalar(@{$sim_record}) > 0 ){
 187
 188         my @new_record;
 189         foreach my $sim_line ( @{$sim_record -> [0]} ){
 190           my $new_line;
 191           while( $sim_line =~ /([^()]*)(\([^()]+\))(.*)/g ){
 192             my $seed = random_uniform_integer( 1, 0, 1000000 ); # Upper limit is from nmhelp
 193
 194             my $head = $1;
 195             my $old_seed = $2;
 196             $sim_line = $3;
 197
 198             if( $old_seed =~ /\(\d+\s+(NORMAL|UNIFORM|NONPARAMETRIC)\)/ ){
 199               $new_line .= "$head($seed $1)";
 200             } else {
 201               $new_line .= "$head($seed)";
 202             }
 203
 204           }
 205
 206           push( @new_record, $new_line.$sim_line );
 207         }
 208
 209         $prob -> set_records( type => 'simulation',
 210                               record_strings => \@new_record );
 211
 212       } else {
 213
 214         my $seed = random_uniform_integer( 1, 0, 1000000 ); # Upper limit is from nmhelp
 215         $prob -> set_records( type           => 'simulation',
 216                               record_strings => [ '(' . $seed .
 217                                                   ') ONLYSIMULATION' ] );
 218       }
 219
 220       if( $sim_no == 1 ) {
 221         # remove $EST and $COV
 222
 223         if( $sim_model -> is_option_set( problem_number => $j,
 224                                          record => 'estimation',
 225                                          name => 'LIKELIHOOD',
 226                                          fuzzy_match => 1 )
 227             or
 228             $sim_model -> is_option_set( problem_number => $j,
 229                                          record => 'estimation',
 230                                          name => '-2LOGLIKELIHOOD',
 231                                          fuzzy_match => 1 )
 232             or
 233             $sim_model -> is_option_set( problem_number => $j,
 234                                          record => 'estimation',
 235                                          name => '-2LLIKELIHOOD',
 236                                          fuzzy_match => 1 )
 237             ){
 238
 239           $sim_model -> set_option( problem_numbers => [$j],
 240                                     record_name => 'estimation',
 241                                     option_name => 'MAXEVALS',
 242                                     fuzzy_match => 1,
 243                                     option_value => 0 );
 244
 245         } else {
 246
 247           $prob -> remove_records(type => 'estimation');
 248
 249         }
 250
 251         $prob -> remove_records(type => 'covariance');
 252
 253         # set $TABLE record
 254         if( defined $prob -> inputs and defined $prob -> inputs -> [0] -> options ) {
 255           foreach my $option ( @{$prob -> inputs -> [0] -> options} ) {
 256             push( @table_header, $option -> name );
 257           }
 258         } else {
 259           debug -> die( message => "Trying to construct table for monte-carlo simulation".
 260                         " but no headers were found in \$model_number-1NPUT" );
 261         }
 262       }
 263       my $simulated_file = "mc-sim-".($j+1)."-$sim_no.dat";
 264       $prob -> set_records( type           => 'table',
 265                             record_strings => [ join( ' ', @table_header ).
 266                                                 ' NOPRINT NOAPPEND ONEHEADER FILE='.
 267                                                 $simulated_file ] );
 268       push( @simulated_files, $self -> {'directory'}.'m'.$model_number.'/'.
 269             $simulated_file );
 270     }
 271     $sim_model -> _write( write_data => $sim_no == 1 ? 1 : 0 );
 272     push( @sim_models, $sim_model );
 273     push( @orig_est_models, $est_original );
 274     push( @all_simulated_files, \@simulated_files );
 275   }
 276   my $mod_sim = tool::modelfit -> new( reference_object => $self,
 277                                        models           => \@sim_models,
 278                                        base_directory   => $self -> {'directory'},
 279                                        directory        => $self -> {'directory'}.
 280                                        'simulation_dir'.$model_number,
 281                                        parent_tool_id   => $self -> {'tool_id'},
 282                                        retries          => 1,
 283                                        logfile           => undef,
 284                                        raw_results           => undef,
 285                                        prepared_models       => undef,
 286                                        threads          => 1 );
 287
 288   $mod_sim -> run;
 289   for( my $j = 0; $j <= $#orig_est_models; $j++ ) {
 290     $orig_est_models[$j] -> datafiles( new_names => $all_simulated_files[$j] );
 291     $orig_est_models[$j] -> _write;
 292     $orig_est_models[$j] -> flush_data();
 293   }
 294   @{$self -> {'mc_models'}} = ( @orig_est_models );
 295
 296   foreach $alternative (@alternatives){
 297       @alt_table_names = ();
 298       @alt_est_models =();
 299       my $filestem = 'mc-alt-' . $alternative -> filename();
 300       $filestem =~ s/\.[^.]*$//;
 301       for( my $sim_no = 1; $sim_no <= $self -> {'samples'}; $sim_no++ ) {
 302           my $alt_name = "$filestem-$sim_no.mod";
 303           my $alt_out = "$filestem-$sim_no.lst";
 304
 305           if( $sim_no == 1 ) {
 306               $est_alternative = $alternative ->
 307                   copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$alt_name,
 308                         target      => 'disk',
 309                         copy_data   => 0,
 310                         copy_output => 0);
 311               $est_alternative -> drop_dropped;
 312               $est_alternative -> remove_records( type => 'simulation' );
 313
 314               #ignore @ since simdata contains header rows. cannot skip other old ignores!!
 315
 316               my $ig_list = $est_alternative -> get_option_value( record_name  => 'data',
 317                                                                   option_name  => 'IGNORE',
 318                                                                   option_index => 'all');
 319               $est_alternative -> remove_option( record_name  => 'data',
 320                                                  option_name  => 'IGNORE');
 321
 322               if (scalar (@{$ig_list})>0){
 323                   foreach my $val (@{$ig_list}){
 324                       unless ($val =~ /^.$/){
 325                           $est_alternative -> add_option( record_name  => 'data',
 326                                                           option_name  => 'IGNORE',
 327                                                           option_value => $val);
 328                       }
 329                   }
 330               }
 331               $est_alternative -> add_option( record_name  => 'data',
 332                                               option_name  => 'IGNORE',
 333                                               option_value => '@');
 334               ##done fixing ignore
 335
 336               my $tbl_nm_ref = $alternative -> table_names();
 337               @alt_table_names = @{$tbl_nm_ref} if( defined $tbl_nm_ref );
 338           } else {
 339               $est_alternative = $alt_est_models[0] ->
 340                   copy( filename    => $self -> {'directory'}.'m'.$model_number.'/'.$alt_name,
 341                         target      => 'disk',
 342                         copy_data   => 0,
 343                         copy_output => 0);
 344           }
 345           $est_alternative -> ignore_missing_files( 1 );
 346           $est_alternative -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$alt_out );
 347           $est_alternative -> ignore_missing_files( 0 );
 348           my @new_alt_table_names;
 349           for( my $pr = 0; $pr <= $#alt_table_names; $pr++ ) {
 350               for( my $tbl = 0; $tbl < scalar @{$alt_table_names[$pr]}; $tbl++ ) {
 351                   $new_alt_table_names[$pr][$tbl] = $alt_table_names[$pr][$tbl].'-'.$sim_no;
 352               }
 353           }
 354           if( $#new_alt_table_names >= 0 ) {
 355               $est_alternative -> table_names( new_names => \@new_alt_table_names );
 356           }
 357           if( $self -> shrinkage() ) {
 358               my @problems = @{$est_alternative -> problems};
 359               for( my $i = 1; $i <= scalar @problems; $i++ ) {
 360                   $problems[ $i-1 ] -> shrinkage_module -> model( $est_alternative );
 361               }
 362           }
 363           push( @alt_est_models, $est_alternative );
 364       } #end loop over sim_no
 365       #is all_simulated_files still in scope? should be, declared inside if not $done
 366       for( my $j = 0; $j <= $#alt_est_models; $j++ ) {
 367           $alt_est_models[$j] -> datafiles( new_names => $all_simulated_files[$j] );
 368           $alt_est_models[$j] -> _write;
 369           $alt_est_models[$j] -> flush_data();
 370       }
 371       push (@{$self -> {'mc_models'}}, @alt_est_models);
 372
 373
 374   }#end loop over alternatives
 375
 376   # Create a checkpoint.
 377   open( DONE, ">".$self -> {'directory'}."/m$model_number/done" ) ;
 378   print DONE "Simulation from ",$sim_models[0] -> filename," through ",
 379   $sim_models[$#sim_models] -> filename," performed\n";
 380   print DONE $self -> {'samples'}." samples\n";
 381   @seed = random_get_seed;
 382   print DONE "seed: @seed\n";
 383   close( DONE );
 384 } else {
 385   #$done=true
 386   # Recreate the datasets and models from a checkpoint
 387   ui -> print( category => 'mc',
 388                message  => "Recreating models from a previous run" );
 389   open( DONE, $self -> {'directory'}."/m$model_number/done" );
 390   my @rows = <DONE>;
 391   close( DONE );
 392   my ( $junk, $junk, $stored_filename1, $junk,
 393        $stored_filename2, $junk ) = split(' ',$rows[0],4);
 394   my ( $stored_samples, $junk ) = split(' ',$rows[1],2);
 395   @seed = split(' ',$rows[2]);
 396   shift( @seed ); # get rid of 'seed'-word
 397   # Reinitiate the model objects
 398   for ( my $j = 1; $j <= $stored_samples; $j++ ) {
 399     my ($model_dir, $orig_name) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 400                                                             $model_number,
 401                                                             "mc-original-$j.mod" );
 402     my ($out_dir, $orig_out) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 403                                                           $model_number,
 404                                                           "mc-original-$j.lst" );
 405     my $est_original = model ->
 406         new( directory   => $model_dir,
 407              filename    => $orig_name,
 408              outputfile  => $orig_out,
 409              extra_files => $model -> extra_files,
 410              target      => 'disk',
 411              ignore_missing_files => 1,
 412              extra_output => $model -> extra_output(),
 413              cwres       => $model -> cwres(),
 414              quick_reload => 1 );
 415     push( @orig_est_models, $est_original );
 416     my $nl = $j == $stored_samples ? "" : "\r";
 417     ui -> print( category => 'mc',
 418                  message  => ui -> status_bar( sofar => $j+1,
 419                                                goal  => $stored_samples+1 ).$nl,
 420                  wrap     => 0,
 421                  newline  => 0 );
 422   } #end loop over samples
 423   @{$self -> {'mc_models'}} = ( @orig_est_models);
 424
 425   #start alternatives section
 426   foreach $alternative (@alternatives){
 427       my $filestem = 'mc-alt-' . $alternative -> filename();
 428       $filestem =~ s/\.[^.]*$//;
 429       @alt_est_models =();
 430       for ( my $j = 1; $j <= $stored_samples; $j++ ) {
 431           my ($model_dir, $alt_name) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 432                                                                   $model_number,
 433                                                                   "$filestem-$j.mod" );
 434           my ($out_dir, $alt_out) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 435                                                                $model_number,
 436                                                                "$filestem-$j.lst" );
 437           my $est_alternative = model ->
 438               new( directory   => $model_dir,
 439                    filename    => $alt_name,
 440                    outputfile  => $alt_out,
 441                    extra_files => $alternative -> extra_files,
 442                    target      => 'disk',
 443                    ignore_missing_files => 1,
 444                    extra_output => $model -> extra_output(),
 445                    cwres       => $model -> cwres(),
 446                    quick_reload => 1 );
 447           push( @alt_est_models, $est_alternative );
 448           my $nl = $j == $stored_samples ? "" : "\r";
 449           ui -> print( category => 'mc',
 450                        message  => ui -> status_bar( sofar => $j+1,
 451                                                      goal  => $stored_samples+1 ).$nl,
 452                        wrap     => 0,
 453                        newline  => 0 );
 454       } #end loop over samples
 455       push (@{$self -> {'mc_models'}}, @alt_est_models );
 456
 457   }#end loop over alternatives
 458
 459   ui -> print( category => 'mc', message  => " ... done." );
 460   random_set_seed( @seed );
 461   ui -> print( category => 'mc',
 462                message  => "Using $stored_samples data sets, previously simulated ".
 463                "from $stored_filename1 through $stored_filename2" )
 464       unless $self -> {'parent_threads'} > 1;
 465 } #end $done
 466 my $subdir = $class;
 467 $subdir =~ s/tool:://;
 468 my @subtools = @{$self -> {'subtools'}};
 469 shift( @subtools );
 470 my %subargs = ();
 471 if ( defined $self -> {'subtool_arguments'} ) {
 472   %subargs = %{$self -> {'subtool_arguments'}};
 473 }
 474
 475 push( @{$self -> {'tools'}},
 476       $class -> new
 477       ( reference_object => $self,
 478         logfile          => undef,
 479         raw_results      => undef,
 480         prepared_models  => undef,
 481         models         => $self -> {'mc_models'},
 482         nm_version     => $self -> {'nm_version'},
 483         parent_tool_id   => $self -> {'tool_id'},
 484         threads        => $subm_threads,
 485         base_directory => $self -> {'directory'},
 486         directory      => $self -> {'directory'}.'/'.$subdir.'_dir'.$model_number,
 487         subtools       => $#subtools >= 0 ? \@subtools : undef,
 488         parent_threads => $own_threads,
 489         shrinkage      => $self -> {'shrinkage'},
 490         _raw_results_callback => $self ->
 491         _modelfit_raw_results_callback( model_number => $model_number ),
 492         %subargs ) );
 493
 494 end general_setup
 495
 496 # }}} general_setup
 497
 498 # {{{ cdd_analyze
 499 start cdd_analyze
 500     $self -> {'tools'}[0] -> print_results;
 501 end cdd_analyze
 502 # }}} cdd_analyze
 503
 504 # {{{ modelfit_analyze
 505
 506 start modelfit_analyze
 507 $self -> {'tools'}[0] -> print_results;
 508 end modelfit_analyze
 509
 510 # }}}
 511
 512 # {{{ _modelfit_raw_results_callback
 513
 514 start _modelfit_raw_results_callback
 515
 516 # Use the mc's raw_results file.
 517 my ($dir,$file) =
 518     OSspecific::absolute_path( $self -> {'directory'},
 519                                $self -> {'raw_results_file'}[$model_number-1] );
 520 my $orig_mod = $self -> {'models'}[$model_number-1];
 521 $subroutine = sub {
 522   my $modelfit = shift;
 523   my $mh_ref   = shift;
 524   my %max_hash = %{$mh_ref};
 525   $modelfit -> raw_results_file( $dir.$file );
 526
 527   # {{{ New header
 528
 529   if ( defined $modelfit -> raw_results() ) {
 530     my @rows = @{$modelfit -> raw_results()};
 531
 532     #the first 'samples' (a scalar) number of rows is null models
 533     for( my $i = 0; $i < $self -> {'samples'}; $i++ ) {
 534       unshift( @{$rows[$i]}, 'null' );
 535     }
 536     #remaining rows, if any, belong to alternatives model(s)
 537     my $offset = $self -> {'samples'};
 538     foreach my $alternative (@{$self -> {'alternative_models'}}){
 539         if (($offset + $self -> {'samples'}-1)> $#rows){
 540           debug -> die( message => "Too few rows in raw results file." );
 541         }
 542         my $altname = 'alt-' . $alternative -> filename();
 543         $altname =~ s/\.[^.]*$//;
 544         for( my $i = 0; $i < $self -> {'samples'}; $i++ ) {
 545             unshift( @{$rows[($offset+$i)]}, $altname );
 546         }
 547         $offset += $self -> {'samples'};
 548     }
 549     if (($offset -1)< $#rows){
 550         debug -> die( message => "Too many rows in raw results file." );
 551     }
 552
 553     unshift( @{$modelfit -> {'raw_results_header'}}, 'hypothesis' );
 554
 555   }
 556   @{$self -> {'raw_results_header'}} = @{$modelfit -> {'raw_results_header'}};
 557   # }}} New header
 558
 559 };
 560 return $subroutine;
 561
 562 end _modelfit_raw_results_callback
 563
 564 # }}} _modelfit_raw_results_callback
 565
 566 # {{{ prepare_results
 567
 568 start prepare_results
 569 {
 570   # Mean
 571   # Median
 572   # SD
 573
 574   $self -> {'raw_results'} = $self -> {'tools'}[0] -> {'raw_results'};
 575
 576   my $model = $self -> {'models'} -> [0];
 577   my @alternatives = @{$self -> {'alternative_models'}};
 578
 579   my $samples = $self -> {'samples'};
 580
 581   my $raw_line_struct = $self -> {'tools'}[0] -> {'raw_line_structure'};
 582
 583   my %orig_results_section;
 584   my @alt_results_sections;
 585
 586   ### Retrieve initial parameter estimates from original model.
 587
 588   my %initial_values;
 589   my %n_initials;
 590   my (@thetalist,@omegalist,@sigmalist);
 591   my $tmp;
 592
 593   $n_initials{'theta'} = $model->nthetas(problem_number => 1);
 594   $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'theta');
 595   @thetalist=@{$tmp->[0]};
 596   'debug'->die(message=>'Error number of initial thetas.') unless (scalar(@thetalist) == $n_initials{'theta'});
 597   $initial_values{'theta'}=\@thetalist;
 598
 599   $tmp = $model->nomegas(problem_numbers => [1]);
 600   $n_initials{'omega'} = $tmp->[0];
 601   $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'omega');
 602   @omegalist=@{$tmp->[0]};
 603   'debug'->die(message=>'Error number of initial omegas.') unless (scalar(@omegalist) == $n_initials{'omega'});
 604   $initial_values{'omega'}=\@omegalist;
 605
 606   $tmp = $model->nsigmas(problem_numbers => [1]);
 607   $n_initials{'sigma'} = $tmp->[0];
 608   $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'sigma');
 609   @sigmalist=@{$tmp->[0]};
 610   'debug'->die(message=>'Error number of initial sigmas.') unless (scalar(@sigmalist) == $n_initials{'sigma'});
 611   $initial_values{'sigma'}=\@sigmalist;
 612
 613   $n_initials{'ofv'}=0;
 614
 615
 616   ### Loop over the original model (the zero) and over all alternative
 617   ### models.
 618
 619   for my $model_index ( 0..scalar(@alternatives) ){
 620
 621     ### Foreach such model, create a %results_section
 622
 623     my %results_section;
 624
 625     ### Each results_section contains statistics of ofv and
 626     ### parameters.
 627
 628     my @values;
 629     my @labels;
 630
 631     @{$labels[0]} = ('mean','median','sd','skewness','kurtosis','rmse','bias');
 632
 633     if( $model_index == 0 ){
 634       $results_section{'name'} = 'Original model';
 635     } else {
 636       $results_section{'name'} = 'Alternative model '.$model_index;
 637     }
 638
 639     $results_section{'labels'} = \@labels;
 640     $results_section{'values'} = \@values;
 641
 642     push( @{$self -> {'results'}[0]{'own'}},\%results_section );
 643
 644     foreach my $measure ( 'ofv','theta','omega','sigma' ){
 645
 646       # Create a header, we might use max_hash from modelfit for this
 647       # later on.
 648
 649       my $start=0;
 650       my $length=0;
 651
 652       for( 0..($samples-1) ){
 653         my ($s, $l) = split(/,/, $raw_line_struct -> {$model_index*$samples + $_} -> {$measure});
 654
 655         if( defined $s and $start != 0 and $s != $start ){
 656           'debug' -> die( message => "Internal structure corrupted at $model_index, $samples, $_, $measure, this is a bug" );
 657         } else {
 658           $start = $s;
 659         }
 660
 661         if( defined $l and $l > $length ){
 662           $length = $l;
 663         }
 664
 665       }
 666
 667       if( $start == 0 ){
 668         for( 0..6 ){
 669           push(@{$values[$_]},'Model not terminated');
 670         }
 671         push( @{$labels[1]},$measure );
 672         next;
 673       }
 674
 675       if( $measure eq 'ofv' ){
 676         push( @{$labels[1]}, $measure );
 677       } else {
 678         for( 1..$length ){
 679           push( @{$labels[1]}, uc(substr($measure, 0,2))."_$_" );
 680         }
 681       }
 682
 683       my $initial_index=0;
 684       foreach my $col ( $start..($start + $length-1) ){
 685         my ($skew, $kurt, $mean, $stdev) = $self -> skewness_and_kurtosis( column_index => $col+1,
 686                                                                            start_row_index => $model_index*$samples,
 687                                                                            end_row_index => $model_index*$samples+$samples-1 );
 688
 689         push( @{$values[0]}, $mean );
 690         push( @{$values[2]}, $stdev );
 691         push( @{$values[3]}, $skew );
 692         push( @{$values[4]}, $kurt );
 693
 694         my $median = $self -> median( column_index => $col+1,
 695                                       start_row_index => $model_index*$samples,
 696                                       end_row_index => $model_index*$samples+$samples-1 );
 697
 698         push( @{$values[1]}, $median );
 699
 700         if( $initial_index < $n_initials{$measure}){
 701
 702           my $initial_value= $initial_values{$measure}->[$initial_index];
 703
 704           $initial_index++;
 705
 706           my $rmse = $self -> rmse_percent( column_index => $col+1,
 707                                             start_row_index => $model_index*$samples,
 708                                             end_row_index => $model_index*$samples+$samples-1,
 709                                             initial_value => $initial_value );
 710
 711           push( @{$values[5]}, $rmse );
 712
 713           my $bias = $self -> bias_percent( column_index => $col+1,
 714                                             start_row_index => $model_index*$samples,
 715                                             end_row_index => $model_index*$samples+$samples-1,
 716                                             initial_value => $initial_value );
 717
 718
 719           push( @{$values[6]}, $bias );
 720         }else {
 721           #push placeholders so that values will stay in correct column
 722           push( @{$values[5]}, ' ' );
 723           push( @{$values[6]}, ' ' );
 724         }
 725
 726       }
 727     }
 728   }
 729
 730   ## Calculate OFV mean calculations
 731
 732   my ($s, $l) = split(/,/, $raw_line_struct -> {0} -> {'ofv'});
 733
 734   my @ofv_limits = (-3.84,-5.99,-7.81,-9.49);
 735
 736   my %results_section;
 737   my (@labels, @values);
 738
 739   @{$labels[1]} = ('mean delta_OFV');
 740
 741   foreach my $limit( @ofv_limits ){
 742     push( @{$labels[1]}, "delta_OFV < $limit" );
 743   }
 744
 745   $results_section{'name'} = 'OFV Statistics';
 746   $results_section{'labels'} = \@labels;
 747   $results_section{'values'} = \@values;
 748
 749   push( @{$self -> {'results'}[0]{'own'}},\%results_section );
 750
 751   foreach my $model_index( 1..scalar(@alternatives) ){
 752
 753     push( @{$labels[0]}, "original - alternative$model_index" );
 754
 755     my $start=0;
 756
 757     my $delta_sum=0;
 758
 759     my %nr_below;
 760
 761     foreach my $sample( 0..($samples-1) ){
 762       my ($s, $l) = split(/,/, $raw_line_struct -> {$model_index*$samples + $sample} -> {'ofv'});
 763
 764       if( $s ne '' and $start != 0 and $s != $start ){
 765         'debug' -> die( message => "OFV: Internal structure corrupted at $model_index, $samples, $sample, this is a bug" );
 766       } else {
 767         $start = $s;
 768       }
 769
 770       if( $start eq '' ){
 771         # This model probably did not terminate.
 772         next;
 773       }
 774
 775       my $alt_ofv = $self -> {'raw_results'} -> [$model_index*$samples + $sample][$s+1];
 776       my $org_ofv = $self -> {'raw_results'} -> [$sample][$s+1];
 777       my $delta = $org_ofv - $alt_ofv;
 778       $delta_sum += $delta;
 779
 780       foreach my $limit ( @ofv_limits ){
 781
 782         if( $delta < $limit ){
 783           $nr_below{$limit}++;
 784         }
 785
 786       }
 787
 788     }
 789
 790     push( @{$values[$model_index-1]}, ($delta_sum / ($samples * scalar(@alternatives))) . "%" );
 791     foreach my $limit( @ofv_limits ){
 792       push( @{$values[$model_index-1]}, ($nr_below{$limit} / ($samples * scalar(@alternatives)))*100 . "%" );
 793     }
 794   }
 795
 796
 797 }
 798 end prepare_results
 799
 800 # }}}
 801
 802 # {{{ rmse_percent
 803
 804 start rmse_percent
 805 {
 806   #input is integers $column_index, $start_row_index, $end_row_index and scalar float $initial_value
 807   #output is scalar $rmse_percent
 808
 809   unless( $end_row_index ){
 810     $end_row_index = $#{$self -> {'raw_results'}};
 811   }
 812
 813   debug -> die( message => "Bad row index input") if ($start_row_index > $end_row_index);
 814   my $row_count = $end_row_index - $start_row_index +1;
 815   my $sum_squared_errors=0;
 816   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
 817     $sum_squared_errors += ($self->{'raw_results'}->[$i][$column_index] - $initial_value)**2;
 818   }
 819
 820   if ($initial_value == 0){
 821     debug -> warn( level => 2,
 822                    message => "Initial value 0, returning absolute rmse instead of relative.");
 823     $rmse_percent= sqrt($sum_squared_errors/$row_count);
 824   }else{
 825     $rmse_percent= (sqrt($sum_squared_errors/$row_count))*100/abs($initial_value);
 826   }
 827 }
 828 end rmse_percent
 829
 830 # }}} rmse_percent
 831
 832 # {{{ bias_percent
 833
 834 start bias_percent
 835 {
 836   #input is integers $column_index, $start_row_index, $end_row_index and scalar float $initial_value
 837   #output is scalar $bias_percent
 838
 839   unless( $end_row_index ){
 840     $end_row_index = $#{$self -> {'raw_results'}};
 841   }
 842
 843   debug -> die( message => "Bad row index input") if ($start_row_index > $end_row_index);
 844
 845   my $row_count = $end_row_index - $start_row_index +1;
 846   my $sum_errors=0;
 847   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
 848     $sum_errors += ($self->{'raw_results'}->[$i][$column_index] - $initial_value);
 849   }
 850
 851   if ($initial_value == 0){
 852     debug -> warn( level => 2,
 853                    message => "Initial value 0, returning absolute bias instead of relative.");
 854     $bias_percent= ($sum_errors/$row_count);
 855   }else{
 856     $bias_percent= ($sum_errors/$row_count)*100/abs($initial_value);
 857   }
 858 }
 859 end bias_percent
 860
 861 # }}} bias_percent
 862
 863 # {{{ skewness_and_kurtosis
 864
 865 start skewness_and_kurtosis
 866 {
 867   #input is integers $column_index, $start_row_index, $end_row_index
 868
 869   unless( $end_row_index ){
 870     $end_row_index = $#{$self -> {'raw_results'}};
 871   }
 872
 873   debug -> die( message => "Bad row index input") if ($start_row_index >= $end_row_index);
 874
 875   my $row_count = $end_row_index - $start_row_index +1;
 876   my $sum_values=0;
 877   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
 878     $sum_values += $self->{'raw_results'}->[$i][$column_index];
 879   }
 880
 881   $mean=$sum_values/$row_count;
 882
 883   my $error=0;
 884   my $sum_errors_pow2=0;
 885   my $sum_errors_pow3=0;
 886   my $sum_errors_pow4=0;
 887
 888   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
 889     $error = ($self->{'raw_results'}->[$i][$column_index] - $mean);
 890     $sum_errors_pow2 += $error**2;
 891     $sum_errors_pow3 += $error**3;
 892     $sum_errors_pow4 += $error**4;
 893   }
 894
 895   ## TODO fråga om missing values. och om SD
 896
 897   $stdev=0;
 898   unless( $sum_errors_pow2 == 0 ){
 899
 900     $stdev= sqrt ($sum_errors_pow2/($row_count-1));
 901     $skewness = $sum_errors_pow3/($row_count*($stdev**3));
 902     $kurtosis = -3 + $sum_errors_pow4/($row_count*($stdev**4));
 903
 904   }
 905 }
 906 end skewness_and_kurtosis
 907
 908 # }}} skewness_and_kurtosis
 909
 910 # {{{ median
 911
 912 start median
 913 {
 914   #input is integers $column_index, $start_row_index, $end_row_index
 915
 916   unless( $end_row_index ){
 917     $end_row_index = $#{$self -> {'raw_results'}};
 918   }
 919
 920   debug -> die( message => "Bad row index input") if ($start_row_index >= $end_row_index);
 921
 922   my @temp;
 923
 924   for (my $i=$start_row_index; $i<=$end_row_index; $i++){
 925     push( @temp, $self->{'raw_results'}->[$i][$column_index] );
 926   }
 927
 928   @temp = sort({$a <=> $b} @temp);
 929   if( scalar( @temp ) % 2 ){
 930     $median = $temp[$#temp/2];
 931   } else {
 932     $median = ($temp[@temp/2]+$temp[(@temp-2)/2]) / 2;
 933   }
 934
 935 }
 936 end median
 937
 938 # }}} median