lib/tool/bootstrap_subs.pm

   1 # {{{ include
   2
   3 start include statements
   4 use strict;
   5 use File::Copy 'cp';
   6 use data;
   7 use OSspecific;
   8 use tool::llp;
   9 use tool::cdd::jackknife;
  10 use ext::Statistics::Distributions 'udistr', 'uprob';
  11 use Math::Random;
  12 use Data::Dumper;
  13 end include
  14
  15 # }}} include statements
  16
  17 # {{{ new
  18
  19 start new
  20       {
  21         foreach my $attribute ( 'logfile', 'raw_results_file','raw_nonp_file' ) {
  22           if ( not( ref($this -> {$attribute}) eq 'ARRAY' or
  23                     ref($this -> {$attribute}) eq 'HASH' ) ) {
  24             my $tmp = $this -> {$attribute};
  25             $this -> {$attribute} = [];
  26             for ( my $i = 1; $i <= scalar @{$this -> {'models'}}; $i++ ) {
  27               my $name = $tmp;
  28               if ( $name =~ /\./ ) {
  29                 $name =~ s/\./$i\./;
  30               } else {
  31                 $name = $name.$i;
  32               }
  33               my $ldir;
  34               ( $ldir, $name ) =
  35                 OSspecific::absolute_path( $this -> {'directory'}, $name );
  36               push ( @{$this -> {$attribute}}, $ldir.$name ) ;
  37             }
  38           }
  39         }
  40       }
  41 end new
  42
  43 # }}} new
  44
  45 # {{{ general_pre_fork_setup
  46
  47 start general_pre_fork_setup
  48     {
  49       # These attributes can be given as a
  50       # 1. A scalar       : used for all models and problems (subjects can be a scalar ref to a hash)
  51       # 2. A 1-dim. array : specified per problem but same for all models
  52       # 3. A 2-dim. array : specified per problem and model
  53
  54       foreach my $attribute_string( 'stratify_on', 'samples', 'subjects' ){
  55         my $attribute = $self -> {$attribute_string};
  56         if ( defined $attribute ) {
  57           # If stratification is given we check at what level it is
  58           # given.
  59           unless ( ref( \$attribute ) eq 'SCALAR' or
  60                    ref( \$attribute ) eq 'REF' or
  61                    ( ref( $attribute ) eq 'ARRAY' and scalar @{$attribute} > 0 ) ) {
  62             # Here we know its neither scalar or a correct array. But we
  63             # seem to assume the lenght be equal to the number of
  64             # models.
  65             debug -> die( message => "attribute $attribute_string is " .
  66                           "defined as a " . ref( $attribute ) .
  67                           " but is neither a scalar or a non-zero size array" );
  68           } elsif ( ref( \$attribute ) eq 'SCALAR' or ref( \$attribute ) eq 'REF' ) {
  69             # If it is a scalar we copy the scalar value into an array
  70             # for each model, with one value for each problem.
  71             my @mo_attribute = ();
  72             foreach my $model ( @{$self -> {'models'}} ) {
  73               my @pr_attribute = ();
  74               foreach my $problem ( @{$model -> problems}  ) {
  75                 push( @pr_attribute, $attribute );
  76               }
  77               push( @mo_attribute, \@pr_attribute );
  78             }
  79             $self -> {$attribute_string} = \@mo_attribute;
  80           } elsif ( ref( $attribute ) eq 'ARRAY' ) {
  81             # If it is an array we check if the it is an array of
  82             # scalars or of arrays. If it is of arrays we seem to
  83             # assume that the arrays have a length the matches the
  84             # number of problems. It it is an array of scalars, we
  85             # copy that scalar into arrays for each problem.
  86             unless ( ref( \$attribute -> [0] ) eq 'SCALAR' or
  87                      ( ref( $attribute -> [0] ) eq 'ARRAY' and scalar @{$attribute -> [0]} > 0 ) ) {
  88               debug -> die( message => "attribute $attribute is ",
  89                             "defined as a ",ref( $attribute -> [0] ),
  90                             "and is neither a scalar or a non-zero size array" );
  91             } elsif ( ref(\$attribute -> [0]) eq 'SCALAR' ) {
  92               my @mo_attribute = ();
  93               foreach my $model ( @{$self -> {'models'}} ) {
  94                 push( @mo_attribute, $attribute );
  95               }
  96               $self -> {$attribute_string} = \@mo_attribute;
  97             }
  98           }
  99         } else {
 100           # If stratify is not defined we copy an undefined value to
 101           # arrays corresponding to models and subproblem. I wonder if
 102           # this is really necessary.
 103           my @mo_attribute = ();
 104           foreach my $model ( @{$self -> {'models'}} ) {
 105             my @pr_attribute = ();
 106             foreach my $data ( @{$model -> datas}  ) {
 107               if( $attribute_string eq 'stratify_on' ){
 108                 push( @pr_attribute, undef );
 109               } elsif ( $attribute_string eq 'samples' ){
 110                 push( @pr_attribute, 200 );
 111               } elsif ( $attribute_string eq 'subjects' ){
 112                 push( @pr_attribute, {'default'=>($data -> count_ind)} );
 113               }
 114             }
 115             push( @mo_attribute, \@pr_attribute );
 116           }
 117           $self -> {$attribute_string} = \@mo_attribute;
 118         }
 119       }
 120     }
 121 end general_pre_fork_setup
 122
 123 # }}}
 124
 125 # {{{ modelfit_pre_fork_setup
 126
 127 start modelfit_pre_fork_setup
 128       {
 129         $self -> general_pre_fork_setup( model_number => $model_number );
 130       }
 131 end modelfit_pre_fork_setup
 132
 133 # }}}
 134
 135 # {{{ llp_pre_fork_setup
 136
 137 start llp_pre_fork_setup
 138       {
 139         $self -> general_pre_fork_setup( model_number => $model_number );
 140       }
 141 end llp_pre_fork_setup
 142
 143 # }}}
 144
 145 # {{{ modelfit_setup
 146
 147 start modelfit_setup
 148       {
 149
 150         my $subm_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
 151           $self -> {'threads'} -> [1]:$self -> {'threads'};
 152         $self -> general_setup( model_number => $model_number,
 153                                 class        => 'tool::modelfit',
 154                                 subm_threads => $subm_threads );
 155       }
 156 end modelfit_setup
 157
 158 # }}}
 159
 160 # {{{ general_setup
 161
 162 start general_setup
 163       {
 164
 165         # If the number of threads are given per tool, e.g. [2,5] meaning 2 threads for
 166         # scm and 5 for the modelfit.
 167         my $subm_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
 168           $self -> {'threads'} -> [1]:$self -> {'threads'};
 169         my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
 170           $self -> {'threads'} -> [0]:$self -> {'threads'};
 171         # More threads than models?
 172         my $num = scalar @{$self -> {'models'}};
 173         $own_threads = $num if ( $own_threads > $num );
 174
 175         # Sub tool threads can be given as scalar or reference to an array?
 176         #my $subm_threads = $parm{'subm_threads'};
 177         #my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
 178         #  $self -> {'threads'} -> [0]:$self -> {'threads'};
 179         # Group variable names are matched in the model, not the data!
 180         my $model       = $self   -> {'models'}      -> [$model_number-1];
 181         my @samples     = @{$self -> {'samples'}     -> [$model_number-1]};
 182         my @stratify_on = @{$self -> {'stratify_on'} -> [$model_number-1]};
 183         my @subjects    = @{$self -> {'subjects'}    -> [$model_number-1]};
 184
 185         # Check which models that hasn't been run and run them This
 186         # will be performed each step but will only result in running
 187         # models at the first step, if at all.
 188
 189         # If more than one process is used, there is a VERY high risk
 190         # of interaction between the processes when creating
 191         # directories for model fits. Therefore the {'directory'}
 192         # attribute is given explicitly below.
 193
 194         # ------------------------  Run original run  -------------------------------
 195
 196         # {{{ orig run
 197
 198         unless ( $model -> is_run ) {
 199           my %subargs = ();
 200           if ( defined $self -> {'subtool_arguments'} ) {
 201             %subargs = %{$self -> {'subtool_arguments'}};
 202           }
 203
 204           if( $self -> {'nonparametric_etas'} or
 205               $self -> {'nonparametric_marginals'} ) {
 206             $model -> add_nonparametric_code;
 207           }
 208
 209           my $orig_fit = tool::modelfit ->
 210               new( reference_object      => $self,
 211                    base_directory        => $self -> {'directory'},
 212                    directory             => $self -> {'directory'}.
 213                    '/orig_modelfit_dir'.$model_number,
 214                    models                => [$model],
 215                    threads               => $subm_threads,
 216                    parent_threads        => $own_threads,
 217                    parent_tool_id        => $self -> {'tool_id'},
 218                    logfile               => undef,
 219                    raw_results           => undef,
 220                    prepared_models       => undef,
 221                    top_tool              => 0,
 222                    %subargs );
 223
 224           ui -> print( category => 'bootstrap',
 225                        message => 'Executing base model.' );
 226
 227           $orig_fit -> run;
 228
 229         }
 230
 231         my $output = $model -> outputs -> [0];
 232
 233         # }}} orig run
 234
 235         # ------------------------  Print a log-header  -----------------------------
 236
 237         # {{{ log header
 238
 239         my @orig_datas = @{$model -> datas};
 240         my @problems   = @{$model -> problems};
 241         my @new_models;
 242
 243         # Print a log-header
 244         # Lasse 2005-04-21: The minimization_message print will probably not work anymore
 245         open( LOG, ">>".$self -> {'logfile'}[$model_number-1] );
 246         my $ui_text = sprintf("%-5s",'RUN').','.sprintf("%20s",'FILENAME  ').',';
 247         print LOG sprintf("%-5s",'RUN'),',',sprintf("%20s",'FILENAME  '),',';
 248         foreach my $param ( 'ofv', 'minimization_message', 'covariance_step_successful' ) {
 249           my $orig_ests   = $model -> outputs -> [0] -> $param;
 250           # Loop the problems
 251           for ( my $j = 0; $j < scalar @{$orig_ests}; $j++ ) {
 252             if ( ref( $orig_ests -> [$j][0] ) ne 'ARRAY' ) {
 253               my $label = uc($param)."_".($j+1);
 254               $ui_text = $ui_text.sprintf("%12s",$label).',';
 255               print LOG sprintf("%12s",$label),',';
 256             } else {
 257               # Loop the parameter numbers (skip sub problem level)
 258               for ( my $num = 1; $num <= scalar @{$orig_ests -> [$j][0]}; $num++ ) {
 259                 my $label = uc($param).$num."_".($j+1);
 260                 $ui_text = $ui_text.sprintf("%12s",$label).',';
 261                 print LOG sprintf("%12s",$label),',';
 262               }
 263             }
 264           }
 265         }
 266         print LOG "\n";
 267
 268         # }}} log header
 269
 270         # ------------------------  Log original run  -------------------------------
 271
 272         # {{{ Log original run
 273
 274         # Lasse 2005-04-21: The minimization_message print will probably not work anymore
 275         open( LOG, ">>".$self -> {'logfile'}[$model_number-1] );
 276         $ui_text = sprintf("%5s",'0').','.sprintf("%20s",$model -> filename).',';
 277         print LOG sprintf("%5s",'0'),',',sprintf("%20s",$model -> filename),',';
 278         foreach my $param ( 'ofv', 'minimization_message', 'covariance_step_successful' ) {
 279           my $orig_ests   = $model -> outputs -> [0] -> $param;
 280           # Loop the problems
 281           for ( my $j = 0; $j < scalar @{$orig_ests}; $j++ ) {
 282             if ( ref( $orig_ests -> [$j][0] ) ne 'ARRAY' ) {
 283               $orig_ests -> [$j][0] =~ s/\n//g;
 284               $ui_text = $ui_text.sprintf("%12s",$orig_ests -> [$j][0]).',';
 285               print LOG sprintf("%12s",$orig_ests -> [$j][0]),',';
 286             } else {
 287               # Loop the parameter numbers (skip sub problem level)
 288               for ( my $num = 0; $num < scalar @{$orig_ests -> [$j][0]}; $num++ ) {
 289                 $ui_text = $ui_text.sprintf("%12f",$orig_ests -> [$j][0][$num]).',';
 290                 print LOG sprintf("%12f",$orig_ests -> [$j][0][$num]),',';
 291               }
 292             }
 293           }
 294         }
 295         print LOG "\n";
 296
 297         # }}} Log original run
 298
 299         # TODO: In this loop we loose one dimension (problem) in that
 300         # all new models with new data sets are put in the same array,
 301         # regardless of which problem the initially belonged to. Fix
 302         # this.
 303
 304         if ( $#orig_datas < 0 ) {
 305           debug -> warn( level   => 1,
 306                          message => "No data files to resample from" );
 307         } else {
 308           debug -> warn( level   => 1,
 309                          message => "Starting bootstrap sampling" );
 310         }
 311         for ( my $i = 1; $i <= scalar @orig_datas; $i++ ) {
 312           my $orig_data = $orig_datas[$i-1];
 313
 314           if ( $self -> {'drop_dropped'} ) {
 315             my $model_copy = $model -> copy( copy_data => 1,
 316                                              filename => "drop_copy_$i.mod",
 317                                              directory => $self -> {'directory'} );
 318
 319             $model_copy -> drop_dropped;
 320             $model_copy -> _write( write_data => 1 );
 321             $model_copy -> datas -> [0] -> flush();
 322             $orig_data = $model -> datas -> [0];
 323             $model = $model_copy;
 324           }
 325
 326           my ( @seed, $new_datas, $incl_ids, $incl_keys, $new_mod );
 327
 328           my $done = ( -e $self -> {'directory'}."/m$model_number/done.$i" ) ? 1 : 0;
 329
 330           if ( not $done ) {
 331             ui -> print( category => 'bootstrap',
 332                          message  => "Resampling from ".$orig_data -> filename );
 333
 334             ( $new_datas, $incl_ids, $incl_keys )
 335                 = $orig_data -> bootstrap( directory   => $self -> {'directory'}.'/m'.$model_number,
 336                                            name_stub   => 'bs_pr'.$i,
 337                                            samples     => $samples[$i-1],
 338                                            subjects    => $subjects[$i-1],
 339                                            stratify_on => $stratify_on[$i-1],
 340                                            target            => 'disk',
 341                                            model_ids   => $self -> {'prepared_model_ids'} );
 342
 343
 344             for ( my $j = 0; $j < $samples[$i-1]; $j++ ) {
 345               my ($model_dir, $filename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.$model_number,
 346                                                                       'bs_pr'.$i.'_'.($j+1).'.mod' );
 347 #             my ($out_dir, $outfilename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.$model_number ,
 348 #                                                                      '/bs_pr'.$i.'_'.$j.'.lst' );
 349               my $prob_copy = Storable::dclone($problems[$i-1]);
 350 #             $Data::Dumper::Maxdepth = 2;
 351 #             print Dumper $model; die;
 352
 353               $new_mod = model ->
 354                   new(reference_object     => $model,
 355                       sde                  => 0,
 356                       outputs              => undef,
 357                       datas                => undef,
 358                       synced               => undef,
 359                       problems             => undef,
 360                       active_problems      => undef,
 361                       directory            => $model_dir,
 362                       filename             => $filename,
 363                       outputfile           => undef,
 364                       problems             => [$prob_copy],
 365                       extra_files          => $model -> extra_files,
 366                       target               => 'disk',
 367                       ignore_missing_files => 1 );
 368
 369               if( $self -> {'shrinkage'} ) {
 370                 $new_mod -> shrinkage_stats( enabled => 1 );
 371                 my @problems = @{$new_mod -> problems};
 372                 for( my $i = 1; $i <= scalar @problems; $i++ ) {
 373                   $problems[ $i-1 ] -> shrinkage_module -> model( $new_mod );
 374                 }
 375               }
 376
 377               my $model_id  = $new_mod -> register_in_database;
 378               push( @{$self -> {'prepared_model_ids'}}, $model_id );
 379
 380               $self -> register_tm_relation( model_ids       => [$model_id],
 381                                              prepared_models => 1 );
 382
 383               $new_mod -> datas ( [$new_datas -> [$j]] );
 384
 385               if( $self -> {'nonparametric_etas'} or
 386                   $self -> {'nonparametric_marginals'} ) {
 387                 $new_mod -> add_nonparametric_code;
 388               }
 389
 390               $new_mod -> update_inits( from_output => $output );
 391               $new_mod -> _write;
 392
 393               push( @new_models, $new_mod );
 394             }
 395             # Create a checkpoint. Log the samples and individuals.
 396             open( DONE, ">".$self -> {'directory'}."/m$model_number/done.$i" ) ;
 397             print DONE "Resampling from ",$orig_data -> filename, " performed\n";
 398             print DONE "$samples[$i-1] samples\n";
 399             while( my ( $strata, $samples ) = each %{$subjects[$i-1]} ) {
 400               print DONE "Strata $strata: $samples sample_size\n";
 401             }
 402             print DONE "Included individuals:\n";
 403             @seed = random_get_seed;
 404             print DONE "seed: @seed\n";
 405             for( my $k = 0; $k < scalar @{$incl_ids}; $k++ ) {
 406               print DONE join(',',@{$incl_ids -> [$k]}),"\n";
 407             }
 408             print DONE "Included keys:\n";
 409             for( my $k = 0; $k < scalar @{$incl_keys}; $k++ ) {
 410               print DONE join(',',@{$incl_keys -> [$k]}),"\n";
 411             }
 412             close( DONE );
 413             open( INCL, ">".$self -> {'directory'}."included_individuals".$model_number.".csv" ) ;
 414             for( my $k = 0; $k < scalar @{$incl_ids}; $k++ ) {
 415               print INCL join(',',@{$incl_ids -> [$k]}),"\n";
 416             }
 417             close( INCL );
 418             open( KEYS, ">".$self -> {'directory'}."included_keys".$model_number.".csv" ) ;
 419             open( SAMPLEKEYS, ">".$self -> {'directory'}."sample_keys".$model_number.".csv" ) ;
 420             my $ninds= ($orig_data -> count_ind());
 421             for( my $k = 0; $k < scalar @{$incl_keys}; $k++ ) {
 422               my %sample_keys;
 423               my $sample_size = scalar @{$incl_keys -> [$k]};
 424               for ( my $l = 0; $l < $ninds; $l++ ) {
 425                 $sample_keys{$incl_keys -> [$k][$l]}++;
 426               }
 427               for ( my $l = 0; $l < $ninds; $l++ ) {
 428                 my $val   = defined $sample_keys{$l} ? $sample_keys{$l} : 0;
 429                 my $extra = ($l == ($ninds-1)) ? "\n" : ',';
 430                 print SAMPLEKEYS $val,$extra;
 431               }
 432               print KEYS join(',',@{$incl_keys -> [$k]}),"\n";
 433             }
 434             close( KEYS );
 435             close( SAMPLEKEYS );
 436           } else {
 437             ui -> print( category => 'bootstrap',
 438                          message  => "Recreating bootstrap from previous run." );
 439
 440             # Recreate the datasets and models from a checkpoint
 441             my ($stored_filename, $stored_samples, %stored_subjects);
 442             my @seed;
 443             my ($stored_filename_found, $stored_samples_found, $stored_subjects_found, $stored_seed_found);
 444             open( DONE, $self -> {'directory'}."/m$model_number/done.$i" );
 445             while( <DONE> ){
 446               if( /^Resampling from (.+) performed$/ ){
 447                 $stored_filename = $1;
 448                 $stored_filename_found = 1;
 449                 next;
 450               }
 451               if( /^(\d+) samples$/ ){
 452                 ui -> print( category => 'bootstrap',
 453                              message  => "Samples: $1" );
 454                 $stored_samples = $1;
 455                 $stored_samples_found = 1;
 456                 next;
 457               }
 458               if( /^(\d+) subjects$/ ){
 459                   # Old format (pre 2.2.2)
 460                 $stored_subjects{'default'} = $1;
 461                 $stored_subjects_found = 1;
 462                 next;
 463               }
 464               if( /^Strata (\w+): (\d+) sample_size$/ ){
 465                 ui -> print( category => 'bootstrap',
 466                              message  => "Strata $1, samples size: $2" );
 467                 $stored_subjects{$1} = $2;
 468                 $stored_subjects_found = 1;
 469                 next;
 470               }
 471               if( /^seed: (\d+) (\d+)$/ ){
 472                 @seed = ($1, $2);
 473                 $stored_seed_found = 1;
 474                 next;
 475               }
 476
 477               if( $stored_filename_found and $stored_samples_found
 478                   and $stored_subjects_found and $stored_seed_found ){
 479                 last;
 480               }
 481             }
 482             close( DONE );
 483             unless( $stored_filename_found and $stored_samples_found
 484                     and $stored_samples_found and $stored_seed_found ){
 485               debug -> die( level => 1,
 486                              message => "The bootstrap/m1/done file could not be parsed." );
 487             }
 488
 489             if ( $stored_samples < $samples[$i-1] ) {
 490               debug -> die( message => "The number of samples saved in previous run ($stored_samples) ".
 491                             "are bigger than the number of samples specified for this run (".
 492                             $samples[$i-1].")" );
 493             }
 494             while( my ( $strata, $samples ) = each %{$subjects[$i-1]} ) {
 495               if ( $stored_subjects{$strata} != $samples ) {
 496                 debug -> die( message => "The number of individuals sampled i strata $strata ".
 497                               "in previous run (".
 498                               $stored_subjects{$strata}.
 499                               ") does not match the number of individuals specified ".
 500                               "for this run (".$samples.")" );
 501               }
 502             }
 503             while( my ( $strata, $samples ) = each %stored_subjects ) {
 504               if ( $subjects[$i-1]->{$strata} != $samples ) {
 505                 debug -> die( message => "The number of individuals sampled i strata $strata ".
 506                               "in previous run (". $samples .
 507                               ") does not match the number of individuals specified ".
 508                               "for this run (".$subjects[$i-1]->{$strata}.")" );
 509               }
 510             }
 511 #           my ( @stored_ids, @stored_keys );
 512 #           for ( my $k = 4; $k < 4+$stored_samples; $k++ ) {
 513 #             chomp($rows[$k]);
 514 #             my @sample_ids = split(',', $rows[$k] );
 515 #             push( @stored_ids, \@sample_ids );
 516 #           }
 517 #           for ( my $k = 5+$stored_samples; $k < 5+2*$stored_samples; $k++ ) {
 518 #             chomp($rows[$k]);
 519 #             my @sample_keys = split(',', $rows[$k] );
 520 #             push( @stored_keys, \@sample_keys );
 521 #           }
 522 #           @seed = split(' ',$rows[5+2*$stored_samples]);
 523 #           $incl_ids = \@stored_ids;
 524 #           $incl_keys = \@stored_keys;
 525 #           shift( @seed ); # get rid of 'seed'-word
 526             # Reinitiate the model objects
 527             for ( my $j = 1; $j <= $samples[$i-1]; $j++ ) {
 528               my ($model_dir, $filename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 529                                                                       $model_number,
 530                                                                       'bs_pr'.$i.'_'.$j.'.mod' );
 531 #             my ($out_dir, $outfilename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
 532 #                                                                      $model_number,
 533 #                                                                      '/bs_pr'.$i.'_'.$j.'.lst' );
 534               $new_mod = model ->
 535                 new( directory   => $model_dir,
 536                      filename    => $filename,
 537 #                    outputfile  => $outfilename,
 538                      extra_files => $model -> extra_files,
 539                      target      => 'disk',
 540                      ignore_missing_files => 1 );
 541 #             $new_mod -> target( 'disk' );
 542               push( @new_models, $new_mod );
 543 #             print "$j\n";
 544 #             print Dumper $new_mod;
 545 #             sleep(10);
 546             }
 547             random_set_seed( @seed );
 548             ui -> print( category => 'bootstrap',
 549                          message  => "Using $stored_samples previously resampled ".
 550                          "bootstrap sets from $stored_filename" )
 551               unless $self -> {'parent_threads'} > 1;
 552           }
 553 #         push( @{$self -> {'included_inds'} -> [$model_number-1]}, $incl_ids );
 554 #         push( @{$self -> {'included_keys'} -> [$model_number-1]}, $incl_keys );
 555         }
 556         $self -> {'prepared_models'}[$model_number-1]{'own'} = \@new_models;
 557
 558         # ---------------------  Create the sub tools  ------------------------------
 559
 560         # {{{ sub tools
 561
 562         my $subdir = $class;
 563         $subdir =~ s/tool:://;
 564         my @subtools = @{$self -> {'subtools'}};
 565         shift( @subtools );
 566         my %subargs = ();
 567         if ( defined $self -> {'subtool_arguments'} ) {
 568           %subargs = %{$self -> {'subtool_arguments'}};
 569         }
 570         push( @{$self -> {'tools'}},
 571               $class ->
 572               new( reference_object      => $self,
 573                    models                => \@new_models,
 574                    threads               => $subm_threads,
 575                    directory             => $self -> {'directory'}.'/'.$subdir.'_dir'.$model_number,
 576                    _raw_results_callback => $self ->
 577                           _modelfit_raw_results_callback( model_number => $model_number ),
 578                    subtools              => \@subtools,
 579                    parent_threads        => $own_threads,
 580                    parent_tool_id        => $self -> {'tool_id'},
 581                    logfile               => $self -> {'logfile'}[$model_number-1],
 582                    raw_results           => undef,
 583                    prepared_models       => undef,
 584                    top_tool              => 0,
 585                    %subargs ) );
 586
 587 #             ( clean            => $self -> {'clean'},
 588 #               base_directory   => $self -> {'directory'},
 589 #               compress         => $self -> {'compress'},
 590 #               directory        => $self -> {'directory'}.'/'.$subdir.'_dir'.$model_number,
 591 #               drop_dropped     => $self -> {'drop_dropped'},
 592 #               wrap_data        => $self -> {'wrap_data'},
 593 #               run_on_nordugrid => $self -> {'run_on_nordugrid'},
 594 #               cpu_time         => $self -> {'cpu_time'},
 595 #               run_on_lsf       => $self -> {'run_on_lsf'},
 596 #               lsf_queue        => $self -> {'lsf_queue'},
 597 #               lsf_options      => $self -> {'lsf_options'},
 598 #               lsf_job_name     => $self -> {'lsf_job_name'},
 599 #               lsf_project_name => $self -> {'lsf_project_name'},
 600
 601 #               parent_tool_id   => $self -> {'tool_id'},
 602
 603 #               models           => \@new_models,
 604 #               grid_batch_size  => $self -> {'grid_batch_size'},
 605 #               nm_version       => $self -> {'nm_version'},
 606 #               picky            => $self -> {'picky'},
 607 #               retries          => $self -> {'retries'},
 608 #               tweak_inits      => $self -> {'tweak_inits'},
 609 #               handle_maxevals  => $self -> {'handle_maxevals'},
 610 #               remove_temp_files        => $self -> {'remove_temp_files'},
 611 #               _raw_results_callback => $self ->
 612 #               _modelfit_raw_results_callback( model_number => $model_number ),
 613 #               threads               => $subm_threads,
 614 #               subtools              => \@subtools,
 615 #               logfile               => $self -> {'logfile'}[$model_number-1],
 616 #               parent_threads        => $own_threads,
 617 #               %subargs ) );
 618
 619         # }}} sub tools
 620
 621       }
 622 end general_setup
 623
 624 # }}}
 625
 626 # {{{ llp_setup
 627 start llp_setup
 628       {
 629         my @subm_threads;
 630         if (ref( $self -> {'threads'} ) eq 'ARRAY') {
 631           @subm_threads = @{$self -> {'threads'}};
 632           unshift(@subm_threads);
 633         } else {
 634           @subm_threads = ($self -> {'threads'});
 635         }
 636         $self -> general_setup( model_number => $model_number,
 637                                 class        => 'tool::llp',
 638                                 subm_threads => \@subm_threads );
 639       }
 640 end llp_setup
 641 # }}} llp_setup
 642
 643 # {{{ _jackknife_raw_results_callback
 644
 645 start _jackknife_raw_results_callback
 646       {
 647         # Use the bootstrap's raw_results file.
 648         my ($dir,$file) =
 649           OSspecific::absolute_path( $self -> {'directory'},
 650                                      $self -> {'raw_results_file'}[$model_number-1] );
 651         my ($dir,$nonp_file) =
 652           OSspecific::absolute_path( $self -> {'directory'},
 653                                      $self -> {'raw_nonp_file'}[$model_number-1] );
 654         $subroutine = sub {
 655           my $jackknife = shift;
 656           my $modelfit  = shift;
 657           $modelfit -> raw_results_file( $dir.$file );
 658           $modelfit -> raw_nonp_file( $dir.$nonp_file );
 659           $modelfit -> raw_results_append( 1 );
 660           my ( @new_header, %param_names );
 661           foreach my $row ( @{$modelfit -> {'raw_results'}} ) {
 662             unshift( @{$row}, 'jackknife' );
 663           }
 664           $modelfit -> {'raw_results_header'} = [];
 665         };
 666         return $subroutine;
 667       }
 668 end _jackknife_raw_results_callback
 669
 670 # }}} _jackknife_raw_results_callback
 671
 672 # {{{ _modelfit_raw_results_callback
 673
 674 start _modelfit_raw_results_callback
 675       {
 676         # Use the bootstrap's raw_results file.
 677         my ($dir,$file) =
 678           OSspecific::absolute_path( $self -> {'directory'},
 679                                      $self -> {'raw_results_file'}[$model_number-1] );
 680         my ($dir,$nonp_file) =
 681           OSspecific::absolute_path( $self -> {'directory'},
 682                                      $self -> {'raw_nonp_file'}[$model_number-1] );
 683         my $orig_mod = $self -> {'models'}[$model_number-1];
 684         my $type = $self -> {'type'};
 685         $subroutine = sub {
 686           my $modelfit = shift;
 687           my $mh_ref   = shift;
 688           my %max_hash = %{$mh_ref};
 689           $modelfit -> raw_results_file( $dir.$file );
 690           $modelfit -> raw_nonp_file( $dir.$nonp_file );
 691
 692           # The prepare_raw_results in the modelfit will fix the
 693           # raw_results for each bootstrap sample model, we must add
 694           # the result for the original model.
 695
 696           my %dummy;
 697
 698           my ($raw_results_row, $nonp_rows) = $self -> create_raw_results_rows( max_hash => $mh_ref,
 699                                                                                 model => $orig_mod,
 700                                                                                 raw_line_structure => \%dummy );
 701
 702           $orig_mod -> outputs -> [0] -> flush;
 703
 704           unshift( @{$modelfit -> {'raw_results'}}, @{$raw_results_row} );
 705
 706           # {{{ New header
 707
 708           if ( $type eq 'bca' ) {
 709             foreach my $row ( @{$modelfit -> {'raw_results'}} ) {
 710               unshift( @{$row}, 'bootstrap' );
 711             }
 712             unshift( @{$modelfit -> {'raw_results_header'}}, 'method' );
 713           }
 714
 715           @{$self -> {'raw_results_header'}} = @{$modelfit -> {'raw_results_header'}};
 716
 717           # }}} New header
 718
 719         };
 720         return $subroutine;
 721       }
 722 end _modelfit_raw_results_callback
 723
 724 # }}} _modelfit_raw_results_callback
 725
 726 # {{{ modelfit_analyze
 727
 728 start modelfit_analyze
 729       {
 730         my @params = @{$self -> {'parameters'}};
 731         my @diagnostic_params = @{$self -> {'diagnostic_parameters'}};
 732         my ( @print_order, @calculation_order );
 733
 734         if ( $self -> {'type'} eq 'bca' ) {
 735
 736           # --------------------------  BCa method  ---------------------------------
 737
 738           # {{{ BCa
 739
 740           @calculation_order = @{$self -> {'bca_calculation_order'}};
 741           @print_order = @{$self -> {'bca_print_order'}};
 742           my $jk_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
 743             $self -> {'threads'} -> [1]:$self -> {'threads'};
 744           my $done = ( -e $self -> {'directory'}."/jackknife_done.$model_number" ) ? 1 : 0;
 745           if ( not $done ) {
 746
 747             # {{{ Create Jackknife
 748
 749             ui -> print( category => 'bootstrap',
 750                          message  => "Running a Jackknife for the BCa estimates" );
 751             $self -> {'jackknife'} = tool::cdd::jackknife ->
 752               new( models           => [$self -> models -> [$model_number -1]],
 753                    case_columns     => $self -> models -> [$model_number -1]
 754                    -> datas -> [0] -> idcolumn,
 755                    _raw_results_callback => $self ->
 756                    _jackknife_raw_results_callback( model_number => $model_number ),
 757                    nm_version       => $self -> {'nm_version'},
 758                    parent_tool_id   => $self -> {'tool_id'},
 759                    threads          => $jk_threads,
 760                    bca_mode         => 1,
 761                    shrinkage        => $self -> {'shrinkage'},
 762                    nonparametric_marginals => $self -> {'nonparametric_marginals'},
 763                    nonparametric_etas => $self -> {'nonparametric_etas'},
 764                    adaptive         => $self -> {'adaptive'},
 765                    rerun            => $self -> {'rerun'},
 766                    verbose          => $self -> {'verbose'},
 767                    cross_validate   => 0 );
 768             # Create a checkpoint. Log the samples and individuals.
 769             open( DONE, ">".$self -> {'directory'}."/jackknife_done.$model_number" ) ;
 770             print DONE "Jackknife directory:\n";
 771             print DONE $self -> {'jackknife'} -> directory,"\n";
 772             my @seed = random_get_seed;
 773             print DONE "seed: @seed\n";
 774             close( DONE );
 775
 776             # }}} Create Jackknife
 777
 778           } else {
 779
 780             # {{{ Recreate Jackknife
 781
 782             open( DONE, $self -> {'directory'}."/jackknife_done.$model_number" );
 783             my @rows = <DONE>;
 784             close( DONE );
 785             my ( $stored_directory ) = $rows[1];
 786             chomp( $stored_directory );
 787             if ( not -e $stored_directory ) {
 788               debug -> die( message => "The Jackknife directory ".$stored_directory.
 789                             "indicated by ".$self -> {'directory'}.
 790                             "/jackknife_done.$model_number".
 791                             " from the old bootstrap run in ".
 792                             $self -> {'directory'}." does not exist" );
 793             }
 794             my @seed = split(' ',$rows[2]);
 795             shift( @seed ); # get rid of 'seed'-word
 796             $self -> {'jackknife'} = tool::cdd::jackknife ->
 797               new( models           => [$self -> models -> [$model_number -1]],
 798                    case_columns     => $self -> models -> [$model_number -1]
 799                    -> datas -> [0] -> idcolumn,
 800                    _raw_results_callback => $self ->
 801                    _jackknife_raw_results_callback( model_number => $model_number ),
 802                    threads          => $jk_threads,
 803                    parent_tool_id   => $self -> {'tool_id'},
 804                    directory        => $stored_directory,
 805                    bca_mode         => 1,
 806                    shrinkage        => $self -> {'shrinkage'},
 807                    nm_version       => $self -> {'nm_version'},
 808                    nonparametric_marginals => $self -> {'nonparametric_marginals'},
 809                    nonparametric_etas => $self -> {'nonparametric_etas'},
 810                    adaptive         => $self -> {'adaptive'},
 811                    rerun            => $self -> {'rerun'},
 812                    verbose          => $self -> {'verbose'},
 813                    cross_validate   => 0 );
 814             random_set_seed( @seed );
 815             ui -> print( category => 'bootstrap',
 816                          message  => "Restarting BCa Jackknife from ".
 817                          $stored_directory )
 818               unless $self -> {'parent_threads'} > 1;
 819
 820             # }}} Recreate Jackknife
 821
 822           }
 823
 824           $self -> {'jackknife'} -> run;
 825           $self -> {'jackknife_results'} = $self -> {'jackknife'} -> {'results'};
 826           $self -> {'jackknife_prepared_models'} = $self -> {'jackknife'} -> {'prepared_models'};
 827
 828           $self -> {'jackknife_raw_results'}[$model_number-1] =
 829             $self -> {'jackknife'} -> raw_results;
 830 #           $self -> {'jackknife'} -> raw_results -> [$model_number-1];
 831 #         $Data::Dumper::Maxdepth = 0;
 832 #         print Dumper $self -> {'jackknife_raw_results'};
 833
 834           # }}} BCa
 835
 836         } else {
 837           @calculation_order = @{$self -> {'calculation_order'}};
 838           @print_order = @{$self -> {'print_order'}};
 839           $self -> {'bootstrap_raw_results'}[$model_number-1] =
 840             $self -> {'tools'} -> [0] -> raw_results;
 841 #           $self -> {'tools'} -> [0] -> raw_results -> [$model_number-1];
 842         }
 843         unless( ref($self -> {'raw_results_header'}[0]) eq 'ARRAY' ) {
 844           my $tmp = $self -> {'raw_results_header'};
 845           $self -> {'raw_results_header'} = [];
 846           $self -> {'raw_results_header'}[$model_number-1] = $tmp;
 847         }
 848
 849         my @param_names = @{$self -> models -> [$model_number -1] -> outputs -> [0] -> labels};
 850         my ( @diagnostic_names, @tmp_names );
 851         foreach my $param ( @diagnostic_params ) {
 852           push( @tmp_names, $param );
 853           $tmp_names[$#tmp_names] =~ s/_/\./g;
 854         }
 855         for ( my $i = 0; $i <= $#param_names; $i++ ) {
 856           unshift( @{$param_names[$i]}, 'OFV' );
 857           push( @{$diagnostic_names[$i]}, @tmp_names );
 858         }
 859         if( defined $PsN::config -> {'_'} -> {'R'} and
 860             -e $PsN::lib_dir . '/R-scripts/bootstrap.R' ) {
 861           # copy the bootstrap R-script
 862           cp ( $PsN::lib_dir . '/R-scripts/bootstrap.R', $self -> {'directory'} );
 863           # Execute the script
 864           system( $PsN::config -> {'_'} -> {'R'}." CMD BATCH bootstrap.R" );
 865         }
 866       }
 867 end modelfit_analyze
 868
 869 # }}}
 870
 871 # {{{ prepare_results
 872
 873 start prepare_results
 874       {
 875         # {{{ definitions
 876
 877         # The '3' is there to skip model, problem and subproblem numbers
 878         my $skip_mps = 3;
 879         my ( @calculation_order, @print_order, %diag_idx );
 880         if ( $self -> {'type'} eq 'bca' ) {
 881           @calculation_order = @{$self -> {'bca_calculation_order'}};
 882           @print_order = @{$self -> {'bca_print_order'}};
 883         } else {
 884           @calculation_order = @{$self -> {'calculation_order'}};
 885           @print_order = @{$self -> {'print_order'}};
 886         }
 887         if ( $self -> {'type'} eq 'bca' ) {
 888 #         $self -> read_raw_results();
 889           $self -> bca_read_raw_results();
 890 #           if ( not defined $self -> {'bootstrap_raw_results'} );
 891         } else {
 892 #         if ( not defined $self -> {'raw_results'} ) {
 893             $self -> read_raw_results();
 894             $self -> {'bootstrap_raw_results'} = $self -> {'raw_results'};
 895 #         }
 896         }
 897
 898         for ( my $i = 0; $i < scalar @{$self -> {'diagnostic_parameters'}}; $i++ ) {
 899           $diag_idx{$self -> {'diagnostic_parameters'} -> [$i]} = $i + $skip_mps;
 900         }
 901
 902         # }}} definitions
 903
 904         # ---------------------  Get data from raw_results  -------------------------
 905
 906         # Divide the data into diagnostics and estimates. Included in estimates are
 907         # the parametric estimates, the standard errors of these, the nonparametric
 908         # estimates, the shrinkage in eta and the shrinkage in wres
 909         # The diagnostics end up in {'bootstrap_diagnostics'} and
 910         # {'jackknife_diagnostics'}. The estimates in {'bootstrap_estimates'} and
 911         # {'jackknife_estimates'}.
 912         # The number of runs that are selected for calculation of the results is
 913         # saved.
 914
 915         # {{{ Get the data from the runs
 916
 917         foreach my $tool ( 'bootstrap', 'jackknife' ) {
 918           if ( defined $self -> {$tool.'_raw_results'} ) {
 919             for ( my $i = 0; $i < scalar @{$self->{$tool.'_raw_results'}}; $i++ ) { # All models
 920
 921               # {{{ Get the number of columns with estimates
 922
 923               my $cols_orig = 0;
 924               foreach my $param ( 'theta', 'omega', 'sigma' ) {
 925                 my $labels =
 926                   $self -> {'models'} -> [$i] -> labels( parameter_type => $param );
 927                 # we can't use labels directly since different models may have different
 928                 # labels (still within the same modelfit)
 929                 my $numpar = scalar @{$labels -> [0]} if ( defined $labels and
 930                                                            defined $labels -> [0] );
 931                 $cols_orig += ( $numpar*3 ); # est + SE + eigen values
 932               }
 933               # nonparametric omegas and shrinkage
 934               my $nomegas = $self -> {'models'} -> [$i] -> nomegas;
 935               my $numpar = $nomegas -> [0];
 936
 937               # shrinkage omega + wres shrinkage
 938               $cols_orig += $numpar + 1;
 939 #             $cols_orig += ($numpar*($numpar+1)/2 + $numpar + 1);
 940
 941               $cols_orig++; # OFV
 942
 943               # }}} Get the number of columns with estimates
 944
 945               # {{{ Loop, choose and set diagnostics and estimates
 946
 947               my %return_section;
 948               $return_section{'name'} = 'Warnings';
 949               my ( $skip_term, $skip_cov, $skip_warn, $skip_bound );
 950               my $included = 0;
 951 #             print Dumper $self->{$tool.'_raw_results'};
 952               for ( my $j = 0; $j < scalar @{$self->{$tool.'_raw_results'}->[$i]}; $j++ ) { # orig model + prepared_models
 953                 my $columns = scalar @{$self->{$tool.'_raw_results'}->[$i][$j]};
 954
 955                 # -----------------------  Diagnostics  -----------------------------
 956
 957                 for ( my $m = $skip_mps; $m < scalar @{$self -> {'diagnostic_parameters'}} + $skip_mps; $m++ ) { # value
 958                   $self->{$tool.'_diagnostics'}->[$i][$j][$m-$skip_mps] =
 959                       $self->{$tool.'_raw_results'}->[$i][$j][$m];
 960                 }
 961                 my $use_run = 1;
 962                 if ( $self -> {'skip_minimization_terminated'} and
 963                      ( not defined $self->{$tool.'_raw_results'}->
 964                        [$i][$j][$diag_idx{'minimization_successful'}]
 965                        or not $self->{$tool.'_raw_results'}->
 966                        [$i][$j][$diag_idx{'minimization_successful'}] ) ) {
 967                   $skip_term++;
 968                   $use_run = 0;
 969                 } elsif ( $self -> {'skip_covariance_step_terminated'} and not
 970                           $self->{$tool.'_raw_results'}->
 971                           [$i][$j][$diag_idx{'covariance_step_successful'}] ) {
 972                   $skip_cov++;
 973                   $use_run = 0;
 974                 } elsif ( $self -> {'skip_with_covstep_warnings'} and
 975                           $self->{$tool.'_raw_results'}->
 976                           [$i][$j][$diag_idx{'covariance_step_warnings'}] ) {
 977                   $skip_warn++;
 978                   $use_run = 0;
 979                 } elsif ( $self -> {'skip_estimate_near_boundary'} and
 980                           $self->{$tool.'_raw_results'}->
 981                           [$i][$j][$diag_idx{'estimate_near_boundary'}] ) {
 982                   $skip_bound++;
 983                   $use_run = 0;
 984                 }
 985
 986                 # ------------------------  Estimates  ------------------------------
 987
 988                 if( $use_run ) {
 989                   for ( my $m = scalar @{$self -> {'diagnostic_parameters'}} + $skip_mps; $m < $columns; $m++ ) { # value
 990                     my $val = $self->{$tool.'_raw_results'}->[$i][$j][$m];
 991                     $self->{$tool.'_estimates'}->
 992                         [$i][$included][$m-(scalar @{$self -> {'diagnostic_parameters'}} + $skip_mps)] = $val;
 993                   }
 994                   $included++;
 995                 }
 996               }
 997
 998               # }}} Loop, choose and set diagnostics and estimates
 999
1000               # {{{ push #runs to results
1001
1002               if ( defined $skip_term ) {
1003                 push( @{$return_section{'values'}}, "$skip_term runs with miminization ".
1004                       "terminated were skipped when calculating the $tool results" );
1005               }
1006               if ( defined $skip_cov ) {
1007                 push( @{$return_section{'values'}}, "$skip_cov runs with aborted ".
1008                       "covariance steps were skipped when calculating the $tool results" );
1009               }
1010               if ( defined $skip_warn ) {
1011                 push( @{$return_section{'values'}}, "$skip_warn runs with errors from ".
1012                       "the covariance step were skipped when calculating the $tool results" );
1013               }
1014               if ( defined $skip_bound ) {
1015                 push( @{$return_section{'values'}}, "$skip_bound runs with estimates ".
1016                       "near a boundary were skipped when calculating the $tool results" );
1017               }
1018               $return_section{'labels'} = [];
1019               push( @{$self -> {'results'}[$i]{'own'}},\%return_section );
1020
1021               # }}} push #runs to results
1022
1023             }
1024           }
1025         }
1026
1027 #         $Data::Dumper::Maxdepth = 5;
1028 #         die Dumper $self -> {'bootstrap_diagnostics'};
1029
1030         # }}} Get the data from the runs
1031
1032         # ----------------------  Calculate the results  ----------------------------
1033
1034         # {{{ Result calculations
1035
1036         for ( my $i = 0; $i < scalar @{$self -> {'bootstrap_raw_results'}} ; $i++ ) { # All models
1037
1038           my $mps_offset = $self -> {'bca'} ? 4 : 3; # <- this is the offset to
1039                                                      # diagonstic_parameters,
1040                                                      # which is one more for
1041                                                      # the method column added
1042                                                      # with a bca run.
1043
1044           my @param_names = @{$self -> {'raw_results_header'}[$i]}[($mps_offset + scalar @{$self -> {'diagnostic_parameters'}}) .. (scalar @{$self -> {'raw_results_header'}[$i]} - 1)];
1045           my ( @diagnostic_names, @tmp_names );
1046           foreach my $param ( @{$self -> {'diagnostic_parameters'}} ) {
1047             push( @tmp_names, $param );
1048             $tmp_names[$#tmp_names] =~ s/_/\./g;
1049           }
1050
1051           @diagnostic_names = @tmp_names;
1052           foreach my $result_type ( @calculation_order ) {
1053             my @names = $result_type eq 'diagnostic_means' ?
1054                 @diagnostic_names : @param_names;
1055             my $calc = 'calculate_'.$result_type;
1056             $self -> $calc( model_number    => ($i+1),
1057                             parameter_names => \@names );
1058           }
1059           foreach my $result_type ( @print_order ) {
1060             my $name = $result_type;
1061             $name =~ s/_/\./g;
1062             my %return_section;
1063             $return_section{'name'} = $name;
1064             $return_section{'values'} = $self -> {$result_type} -> [$i];
1065             $return_section{'labels'} = $self -> {$result_type.'_labels'} -> [$i];
1066             push( @{$self -> {'results'}[$i]{'own'}},\%return_section );
1067           }
1068         }
1069
1070         # }}} Result calculations
1071
1072       }
1073 end prepare_results
1074
1075 # }}} prepare_results
1076
1077 # {{{ print_summary
1078
1079 start print_summary
1080       {
1081         sub acknowledge {
1082           my $name    = shift;
1083           my $outcome = shift;
1084           my $file    = shift;
1085           my $l = (7 - length( $outcome ))/2;
1086           my $c_num = '00';
1087           $c_num    = '07' if ( $outcome eq 'OK' );
1088           $c_num    = '13' if ( $outcome eq 'WARNING' );
1089           $c_num    = '05' if ( $outcome eq 'ERROR' );
1090           # my $text = sprintf( "%-66s%2s%7s%-5s\n\n", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1091           my $text = sprintf( "%-66s%2s%7s%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1092           #  cprintf( "%-66s%2s\x03$c_num%7s\x030%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1093           #  my $text = cprintf( "%-66s%2s\x03$c_num%7s\x030%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1094           print $text, "\n\n";
1095           print $file $text if defined $file;
1096         }
1097
1098         sub sum {
1099           my $arr = shift;
1100           my $sum = 0;
1101           for ( @{$arr} ) {
1102             $sum += $_;
1103           }
1104           return $sum;
1105         }
1106
1107         my $diag_number = scalar @{$self -> {'diagnostic_parameters'}} - 1;
1108         my %diag_idxs;
1109         for ( my $i = 0; $i <= $diag_number; $i++ ) {
1110           $diag_idxs{$self -> {'diagnostic_parameters'} -> [$i]} = $i;
1111         }
1112
1113         open( my $log, ">test.log" );
1114         for ( my $i = 0; $i < scalar @{$self -> {'bootstrap_raw_results'}} ; $i++ ) { # All models
1115           print "MODEL ",$i+1,"\n" if ( scalar @{$self -> {'bootstrap_raw_results'}} > 1 );
1116           if ( $self -> {'diagnostic_means'} ->
1117                [$i][0][$diag_idxs{'minimization_successful'}] >
1118                $self -> {'minimization_successful_limit'} ) {
1119             acknowledge( 'Successful minimization ratio = '.
1120                          $self -> {'diagnostic_means'} ->
1121                          [$i][0][$diag_idxs{'minimization_successful'}], 'OK', $log );
1122           } else {
1123             acknowledge( 'Termination problems in '.
1124                          sprintf("%4.2f", (100-($self -> {'diagnostic_means'} ->
1125                                                 [$i][0][$diag_idxs{'minimization_successful'}]*100)))
1126                          .'% of the bootstrap runs', 'WARNING', $log );
1127           }
1128
1129           if ( $self -> {'diagnostic_means'} ->
1130                [$i][0][$diag_idxs{'covariance_step_successful'}] >
1131                $self -> {'covariance_step_successful_limit'} ) {
1132             acknowledge( 'Successful covariance step ratio = '.$self -> {'diagnostic_means'} ->
1133                          [$i][0][$diag_idxs{'covariance_step_successful'}], 'OK', $log );
1134           } else {
1135             acknowledge( 'Covariance step terminated in '.
1136                          sprintf("%4.2f", (100-($self -> {'diagnostic_means'} ->
1137                                                 [$i][0][$diag_idxs{'covariance_step_successful'}]*100)))
1138                          .'% of the bootstrap runs', 'WARNING', $log );
1139           }
1140
1141           if ( $self -> {'diagnostic_means'} ->
1142                [$i][0][$diag_idxs{'covariance_step_warnings'}] <
1143                $self -> {'covariance_step_warnings_limit'} ) {
1144             acknowledge( 'Covariance step warnings ratio = '.$self -> {'diagnostic_means'} ->
1145                          [$i][0][$diag_idxs{'covariance_step_warnings'}], 'OK', $log );
1146           } else {
1147             acknowledge( 'Covariance step warnings in '.
1148                          sprintf("%4.2f", ($self -> {'diagnostic_means'} ->
1149                                            [$i][0][$diag_idxs{'covariance_step_warnings'}]*100))
1150                          .'% of the bootstrap runs', 'WARNING', $log );
1151           }
1152
1153           if ( $self -> {'diagnostic_means'} ->
1154                [$i][0][$diag_idxs{'estimate_near_boundary'}] <
1155                $self -> {'estimate_near_boundary_limit'} ) {
1156             acknowledge( 'Estimate near boundary ratio = '.$self -> {'diagnostic_means'} ->
1157                          [$i][0][$diag_idxs{'estimate_near_boundary'}], 'OK', $log );
1158           } else {
1159             acknowledge( 'Estimate near boundary found in '.
1160                          sprintf("%4.2f", ($self -> {'diagnostic_means'} ->
1161                                            [$i][0][$diag_idxs{'estimate_near_boundary'}]*100))
1162                          .'% of the bootstrap runs', 'WARNING', $log );
1163           }
1164           my $sum = sum( $self -> {'within_se_confidence_intervals'}->
1165                          [$i]{$self -> {'se_confidence_intervals_level'}} );
1166           if ( not defined $sum or $sum < 1 ) {
1167             acknowledge( 'No '.(100-$self -> {'se_confidence_intervals_level'}).
1168                          '% confidence intervals based on the'.
1169                          ' bootstrap standard errors include '.
1170                          $self -> {'se_confidence_intervals_check'}, 'OK', $log );
1171           } else {
1172             acknowledge( "$sum ".(100-$self -> {'se_confidence_intervals_level'}).
1173                          "% confidence intervals based on the bootstrap".
1174                          " SE's include ".
1175                          $self -> {'se_confidence_intervals_check'}, 'WARNING', $log );
1176             my $found =
1177               scalar @{$self -> {'within_se_confidence_intervals'}->
1178                          [$i]{$self -> {'se_confidence_intervals_level'}}} - 1;
1179             for ( my $l = 0; $l <= $found ; $l++ ) {
1180               if ( $self -> {'within_se_confidence_intervals'}->
1181                    [$i]{$self -> {'se_confidence_intervals_level'}}[$l] ) {
1182                 printf( "\t%-20s\n",
1183                         $self -> {'standard_error_confidence_intervals_labels'} ->
1184                         [$i][1][$l] );
1185                 print $log
1186                   sprintf( "\t%-20s\n",
1187                            $self -> {'standard_error_confidence_intervals_labels'} ->
1188                            [$i][1][$l] );
1189                 print "\n" if ( $l == $found );
1190                 print $log "\n" if ( $l == $found );
1191               }
1192             }
1193           }
1194
1195           my $sum = sum( $self -> {'large_bias'}-> [$i][0] );
1196           if ( not defined $sum or $sum < 1 ) {
1197             acknowledge( 'No bias larger than '.
1198                          ($self -> {'large_bias_limit'}*100).'% found', 'OK', $log );
1199           } else {
1200             acknowledge( "$sum estimates were found to have a relative bias > ".
1201                          $self -> {'large_bias_limit'}, 'WARNING', $log );
1202             my $found =
1203               scalar @{$self -> {'large_bias'}->
1204                          [$i][0]} - 1;
1205             for ( my $l = 0; $l <= $found ; $l++ ) {
1206               if ( $self -> {'large_bias'}->
1207                    [$i][0][$l] ) {
1208                 printf( "\t%-20s%3.2f %\n", $self -> {'bias_labels'} ->
1209                         [$i][1][$l], ($self -> {'bias'} -> [$i][0][$l]/
1210                                       $self->{'bootstrap_estimates'} -> [$i][0][$l])*100 );
1211                 #                   print $log
1212                 #                     sprintf( "\t%-20s\n",
1213                 #                              $self -> {'percentile_confidence_intervals_labels'} ->
1214                 #                              [$i][1][$l] );
1215                 print "\n" if ( $l == $found );
1216                 print $log "\n" if ( $l == $found );
1217               }
1218             }
1219           }
1220
1221           if ( $self -> {'type'} eq 'bca' ) {
1222             my $sum = sum( $self -> {'within_bca_confidence_intervals'}->
1223                            [$i]{$self -> {'bca_confidence_intervals_level'}} );
1224             if ( not defined $sum or $sum < 1 ) {
1225               acknowledge( 'No '.(100-$self -> {'bca_confidence_intervals_level'}).
1226                            '% BCa confidence intervals include '.
1227                            $self -> {'bca_confidence_intervals_check'}, 'OK', $log );
1228             } else {
1229               acknowledge( "$sum ".(100-$self -> {'bca_confidence_intervals_level'}).
1230                            "% BCa confidence intervals include ".
1231                            $self -> {'bca_confidence_intervals_check'}, 'WARNING', $log );
1232               my $found =
1233                 scalar @{$self -> {'within_bca_confidence_intervals'}->
1234                            [$i]{$self -> {'bca_confidence_intervals_level'}}} - 1;
1235               for ( my $l = 0; $l <= $found ; $l++ ) {
1236                 if ( $self -> {'within_bca_confidence_intervals'}->
1237                      [$i]{$self -> {'bca_confidence_intervals_level'}}[$l] ) {
1238                   printf( "\t%-20s\n",
1239                           $self -> {'bca_confidence_intervals_labels'} ->
1240                           [$i][1][$l] );
1241                   print $log
1242                     sprintf( "\t%-20s\n",
1243                              $self -> {'bca_confidence_intervals_labels'} ->
1244                              [$i][1][$l] );
1245                   print "\n" if ( $l == $found );
1246                   print $log "\n" if ( $l == $found );
1247                 }
1248               }
1249             }
1250           } else {
1251             my $sum = sum( $self -> {'within_percentiles'}->
1252                            [$i]{$self -> {'percentile_confidence_intervals_level'}} );
1253             if ( not defined $sum or $sum < 1 ) {
1254               acknowledge( 'No '.(100-$self -> {'percentile_confidence_intervals_level'}).
1255                            '% confidence intervals based on the'.
1256                            ' bootstrap percentiles include '.
1257                            $self -> {'percentile_confidence_intervals_check'}, 'OK', $log );
1258             } else {
1259               acknowledge( "$sum ".(100-$self -> {'percentile_confidence_intervals_level'}).
1260                            "% confidence intervals based on the percentiles".
1261                            " include ".
1262                            $self -> {'percentile_confidence_intervals_check'}, 'WARNING', $log );
1263               my $found =
1264                 scalar @{$self -> {'within_percentiles'}->
1265                            [$i]{$self -> {'percentile_confidence_intervals_level'}}} - 1;
1266               for ( my $l = 0; $l <= $found ; $l++ ) {
1267                 if ( $self -> {'within_percentiles'}->
1268                      [$i]{$self -> {'percentile_confidence_intervals_level'}}[$l] ) {
1269                   printf( "\t%-20s\n",
1270                           $self -> {'percentile_confidence_intervals_labels'} ->
1271                           [$i][1][$l] );
1272                   print $log
1273                     sprintf( "\t%-20s\n",
1274                              $self -> {'percentile_confidence_intervals_labels'} ->
1275                              [$i][1][$l] );
1276                   print "\n" if ( $l == $found );
1277                   print $log "\n" if ( $l == $found );
1278                 }
1279               }
1280             }
1281           }
1282         }
1283         close ( $log );
1284       }
1285 end print_summary
1286
1287 # }}} print_summary
1288
1289 # {{{ bca_read_raw_results
1290
1291 start bca_read_raw_results
1292       {
1293         $self -> {'raw_results_header'} = [];
1294         for ( my $i = 1; $i <= scalar @{$self->{'models'}}; $i++ ) { # All models
1295           if ( -e $self -> {'directory'}.'raw_results'.$i.'.csv' ) {
1296             open( RRES, $self -> {'directory'}.'raw_results'.$i.'.csv' );
1297             my @file = <RRES>;
1298             close( RRES );
1299             map { chomp; my @tmp = split(',',$_); $_ = \@tmp } @file ;
1300
1301             my $header = shift @file;
1302
1303             # Get rid of 'method' column
1304             my $cols = scalar(@{$header})-1;
1305             @{$self -> {'raw_results_header'}[$i-1]} = @{$header}[1..$cols];
1306             $self -> {'raw_results'} -> [$i-1] = \@file;
1307             for( my $j = 0; $j <= $#file; $j++ ) {
1308               if ( $file[$j][0] eq 'jackknife' ) {
1309                 shift( @{$file[$j]} );
1310 #               $self -> {'jackknife_raw_results'}[$i-1] = \@file;
1311                 push( @{$self -> {'jackknife_raw_results'}[$i-1]}, $file[$j]);
1312               } else {
1313                 shift( @{$file[$j]} );
1314 #               $self -> {'bootstrap_raw_results'}[$i-1] = \@file;
1315                 push( @{$self -> {'bootstrap_raw_results'}[$i-1]}, $file[$j] );
1316               }
1317             }
1318           }
1319         }
1320       }
1321 end bca_read_raw_results
1322
1323 # }}} bca_read_raw_results
1324
1325 # {{{ calculate_diagnostic_means
1326
1327 start calculate_diagnostic_means
1328       {
1329         my ( @sum, @diagsum, %diag_idx );
1330         for ( my $i = 0; $i < scalar @{$self -> {'diagnostic_parameters'}}; $i++ ) {
1331           $diag_idx{$self -> {'diagnostic_parameters'} -> [$i]} = $i;
1332         }
1333
1334         my $def = 0;
1335         # Prepared model, skip the first (the original)
1336         for ( my $k = 1; $k < scalar @{$self -> {'bootstrap_diagnostics'} ->
1337                                          [$model_number-1]}; $k++ ) {
1338           # Diagnostics
1339           if( defined $self -> {'bootstrap_diagnostics'} ->
1340               [$model_number-1][$k] ) {
1341             $def++;
1342             for ( my $l = 0; $l < scalar @{$self -> {'bootstrap_diagnostics'} ->
1343                                                [$model_number-1][$k]}; $l++ ) {
1344               $sum[$l] += $self -> {'bootstrap_diagnostics'} ->
1345                   [$model_number-1][$k][$l];
1346             }
1347           }
1348         }
1349
1350         # divide by the number of bootstrap samples (-1 to get rid of the original
1351         # model) The [0] in the index is there to indicate the 'model' level. Mostly
1352         # used for printing
1353         for ( my $l = 0; $l <= $#sum; $l++ ) {
1354           if( $l == $diag_idx{'significant_digits'} ) {
1355             $self -> {'diagnostic_means'} -> [$model_number-1][0][$l] =
1356                 $sum[$l] / $def;
1357           } else {
1358             $self -> {'diagnostic_means'} -> [$model_number-1][0][$l] =
1359                 $sum[$l] / ( scalar @{$self -> {'bootstrap_diagnostics'} ->
1360                                           [$model_number-1]} - 1);
1361           }
1362         }
1363         $self -> {'diagnostic_means_labels'} -> [$model_number-1] =
1364           [[],\@parameter_names];
1365       }
1366 end calculate_diagnostic_means
1367
1368 # }}} calculate_diagnostic_means
1369
1370 # {{{ calculate_means
1371
1372 start calculate_means
1373       {
1374         my ( @sum, @diagsum );
1375         # Prepared model, skip the first (the original)
1376         for ( my $k = 1; $k < scalar @{$self -> {'bootstrap_estimates'} ->
1377                                          [$model_number-1]}; $k++ ) {
1378           # Estimates
1379           for ( my $l = 0; $l < scalar @{$self -> {'bootstrap_estimates'} ->
1380                                            [$model_number-1][$k]}; $l++ ) {
1381             $sum[$l] += $self -> {'bootstrap_estimates'} ->
1382               [$model_number-1][$k][$l];
1383           }
1384         }
1385         # divide by the number of bootstrap samples (-1 to get rid of the original
1386         # model) The [0] in the index is there to indicate the 'model' level. Mostly
1387         # used for printing
1388         my $samples = scalar @{$self -> {'bootstrap_estimates'} ->
1389                                  [$model_number-1]} - 1;
1390         for ( my $l = 0; $l <= $#sum; $l++ ) {
1391           my $mean = $sum[$l] / $samples;
1392           $self -> {'means'} -> [$model_number-1][0][$l] = $mean;
1393           my $bias = $mean - $self ->
1394             {'bootstrap_estimates'} -> [$model_number-1][0][$l];
1395           $self -> {'bias'} -> [$model_number-1][0][$l] = $bias;
1396           if ( $self->{'bootstrap_estimates'} -> [$model_number-1][0][$l] != 0 and
1397                $bias/$self->{'bootstrap_estimates'} -> [$model_number-1][0][$l]
1398                > $self -> {'large_bias_limit'} ) {
1399             $self -> {'large_bias'} -> [$model_number-1][0][$l] = 1;
1400           } else {
1401             $self -> {'large_bias'} -> [$model_number-1][0][$l] = 0;
1402           }
1403         }
1404         $self -> {'means_labels'} -> [$model_number-1] =
1405           [[],\@parameter_names];
1406
1407         $self -> {'bias_labels'} -> [$model_number-1] =
1408           [[],\@parameter_names];
1409       }
1410 end calculate_means
1411
1412 # }}} calculate_means
1413
1414 # {{{ calculate_jackknife_means
1415
1416 start calculate_jackknife_means
1417       {
1418         my @sum;
1419         # Prepared model, skip the first (the original)
1420         if( defined $self -> {'jackknife_estimates'} ){
1421           for ( my $k = 1; $k < scalar @{$self -> {'jackknife_estimates'}->[$model_number-1]}; $k++ ) {
1422             # Estimate
1423             for ( my $l = 0; $l <
1424                   scalar @{$self -> {'jackknife_estimates'}->[$model_number-1][$k]}; $l++ ) {
1425               $sum[$l] += $self -> {'jackknife_estimates'}->[$model_number-1][$k][$l];
1426             }
1427           }
1428           # divide by the number of jackknife samples (-1 to get rid of the original model)
1429           # The [0] in the index is there to indicate the 'model' level. Mostly used for printing
1430           for ( my $l = 0; $l <
1431                 scalar @{$self -> {'jackknife_estimates'}->[$model_number-1][0]}; $l++ ) {
1432             if( ( scalar @{$self -> {'jackknife_estimates'}->[$model_number-1]} - 1) != 0 ) {
1433               $self -> {'jackknife_means'} -> [$model_number-1][0][$l] =
1434                   $sum[$l] / ( scalar @{$self -> {'jackknife_estimates'}->[$model_number-1]} - 1);
1435             }
1436           }
1437           $self -> {'jackknife_means_labels'} -> [$model_number-1] = [[],\@parameter_names];
1438         }
1439       }
1440 end calculate_jackknife_means
1441
1442 # }}} calculate_jackknife_means
1443
1444 # {{{ calculate_medians
1445 start calculate_medians
1446       {
1447         my @medians;
1448         # Loop the parameters
1449         for ( my $l = 0; $l < scalar @{$self -> {'bootstrap_estimates'}->
1450                                          [$model_number-1][0]}; $l++ ) {
1451           my @parameter_array;
1452           # From 1 to get rid of original model
1453           for ( my $k = 1; $k < scalar @{$self -> {'bootstrap_estimates'}->
1454                                            [$model_number-1]}; $k++ ) {
1455             $parameter_array[$k-1] =
1456               $self -> {'bootstrap_estimates'}->[$model_number-1][$k][$l];
1457           }
1458           my @sorted = sort {$a <=> $b} @parameter_array;
1459           # median postition is half the ( array length - 1 ).
1460           my $median_position = ( $#sorted ) / 2;
1461           my ($int_med,$frac_med)   = split(/\./, $median_position );
1462           $frac_med = eval("0.".$frac_med);
1463           my $median_low  = $sorted[ $int_med ];
1464           my $median_high = ( $sorted[ $int_med + 1 ] - $sorted[ $int_med ] ) * $frac_med;
1465           $medians[$l] = $median_low + $median_high;
1466         }
1467         # The [0] in the index is there to indicate the 'model' level. Mostly used for printing
1468         $self -> {'medians'} -> [$model_number-1][0] = \@medians;
1469         $self -> {'medians_labels'} -> [$model_number-1] = [[],\@parameter_names];
1470       }
1471 end calculate_medians
1472 # }}} calculate_medians
1473
1474 # {{{ calculate_standard_error_confidence_intervals
1475 start calculate_standard_error_confidence_intervals
1476       {
1477         # Sort the limits from the inside out
1478         my @limits = sort { $b <=> $a } keys %{$self -> {'confidence_limits'}};
1479         foreach my $limit ( @limits ) {
1480           my ( @lower_limits, @upper_limits, @within_ci );
1481           # Loop the estimates of the first (original) model
1482           for ( my $l = 0; $l < scalar @{$self -> {'bootstrap_estimates'}->
1483                                            [$model_number-1][0]}; $l++ ) {
1484             my $lower_limit =
1485               $self -> {'bootstrap_estimates'}->[$model_number-1][0][$l] -
1486                 $self -> {'standard_errors'}->[$model_number-1][0][$l] *
1487                   $self -> {'confidence_limits'} ->{$limit};
1488             my $upper_limit =
1489               $self -> {'bootstrap_estimates'}->[$model_number-1][0][$l] +
1490                 $self -> {'standard_errors'}->[$model_number-1][0][$l] *
1491                   $self -> {'confidence_limits'} ->{$limit};
1492             push( @lower_limits, $lower_limit );
1493             push(  @upper_limits, $upper_limit );
1494             if ( $self -> {'se_confidence_intervals_check'} < $upper_limit and
1495                  $self -> {'se_confidence_intervals_check'} > $lower_limit ) {
1496               push( @within_ci , 1 );
1497             } else {
1498               push( @within_ci , 0 );
1499             }
1500           }
1501           unshift( @{$self -> {'standard_error_confidence_intervals'} ->
1502                        [$model_number-1]}, \@lower_limits );
1503           push( @{$self -> {'standard_error_confidence_intervals'} ->
1504                     [$model_number-1]}, \@upper_limits );
1505           $self -> {'within_se_confidence_intervals'} ->
1506             [$model_number-1]{$limit} = \@within_ci;
1507           unshift( @{$self -> {'standard_error_confidence_intervals_labels'} ->
1508                        [$model_number-1][0]}, ($limit/2).'%' );
1509           push( @{$self -> {'standard_error_confidence_intervals_labels'} ->
1510                     [$model_number-1][0]}, (100-($limit/2)).'%' );
1511           push( @{$self -> {'within_se_confidence_intervals_labels'} ->
1512                     [$model_number-1][0]}, $limit.'%' );
1513         }
1514         $self -> {'standard_error_confidence_intervals_labels'} -> [$model_number-1][1] =
1515           \@parameter_names;
1516         $self -> {'within_se_confidence_intervals_labels'} -> [$model_number-1][1] =
1517           \@parameter_names;
1518       }
1519 end calculate_standard_error_confidence_intervals
1520 # }}} calculate_standard_error_confidence_intervals
1521
1522 # {{{ calculate_standard_errors
1523
1524 start calculate_standard_errors
1525       {
1526         my @se;
1527         # Prepared model, skip the first (the original)
1528         for ( my $k = 1; $k < scalar @{$self -> {'bootstrap_estimates'}->[$model_number-1]}; $k++ ) {
1529           # Estimate
1530           for ( my $l = 0; $l <
1531                 scalar @{$self -> {'bootstrap_estimates'}->[$model_number-1][$k]}; $l++ ) {
1532             $se[$l] += ( $self -> {'bootstrap_estimates'}->[$model_number-1][$k][$l] -
1533                          $self -> {'means'}->[$model_number-1][0][$l] )**2;
1534           }
1535         }
1536         # divide by the number of bootstrap samples -1 (-2 to get rid of the original model)
1537         # The [0] in the index is there to indicate the 'model' level.
1538         for ( my $l = 0; $l <
1539               scalar @{$self -> {'bootstrap_estimates'}->[$model_number-1][0]}; $l++ ) {
1540           my $div = ( scalar @{$self -> {'bootstrap_estimates'}->[$model_number-1]} - 2 );
1541           if( defined $div and not $div == 0 ) {
1542             $self -> {'standard_errors'} -> [$model_number-1][0][$l] =
1543                 ($se[$l] / $div  )**0.5;
1544           } else {
1545             $self -> {'standard_errors'} -> [$model_number-1][0][$l] = undef;
1546           }
1547         }
1548         $self -> {'standard_errors_labels'} -> [$model_number-1] = [[],\@parameter_names];
1549       }
1550 end calculate_standard_errors
1551
1552 # }}} calculate_standard_errors
1553
1554 # {{{ calculate_bca_confidence_intervals
1555
1556 start calculate_bca_confidence_intervals
1557       {
1558         sub c_get_z0 {
1559           my $arr_ref  = shift;
1560           my $orig_value = shift;
1561           my $num_less_than_orig = 0;
1562           my $nvalues = 0;
1563           my $z0;
1564           foreach my $value ( @{$arr_ref} ) {
1565             if ( defined $value and $value ne '' ) {
1566               $num_less_than_orig++ if ( $value < $orig_value );
1567               $nvalues ++;
1568             }
1569           }
1570
1571           unless ( $nvalues == 0 ) {
1572             if ( ($num_less_than_orig / $nvalues ) == 0 ) {
1573               $z0 = -100;
1574             } elsif ( ($num_less_than_orig / $nvalues ) == 1 ) {
1575               $z0 = 100;
1576             } else {
1577               $z0 = udistr( 1 - ($num_less_than_orig / $nvalues ) );
1578             }
1579           }
1580 #         return ( $z0, $nvalues );
1581 #         print sprintf( "%4s:%5.0f,%4s:%15.8g  ",'N:',$num_less_than_orig,'ZO', $z0);
1582           return $z0;
1583         }
1584
1585         sub c_get_acc {
1586           my $arr_ref = shift;
1587           my $jk_mean = shift;
1588           my $acc_upper = 0;
1589           my $acc_lower = 0;
1590           my $nvalues = 0;
1591           my $acc;
1592           foreach my $value ( @{$arr_ref} ){
1593             if ( defined $value and $value ne '' ) {
1594               $acc_upper = $acc_upper + ($jk_mean-$value)**3;
1595               $acc_lower = $acc_lower + ($jk_mean-$value)**2;
1596               $nvalues ++;
1597             }
1598           }
1599           $acc_lower = 6*($acc_lower**(3/2));
1600           unless ( $acc_lower == 0 ) {
1601             $acc = $acc_upper / $acc_lower;
1602           } else {
1603             $acc = $acc_upper / 0.001;
1604           }
1605 #         return ( $acc, $nvalues );
1606 #         print sprintf( "%4s:%15.8g%4s%8.5g\n",'ACC', $acc,'JKm', $jk_mean);
1607           return $acc;
1608         }
1609
1610         sub c_get_alphas {
1611           my $old_alphas = shift;
1612           my $acc = shift;
1613           my $z0 = shift;
1614           my $denom;
1615           my @new_alphas = ();
1616           foreach my $position ( @{$old_alphas} ) {
1617             if ( $position == 0 ){
1618               $denom = -100;
1619             } elsif ( $position == 100 ) {
1620               $denom = 100;
1621             } else {
1622               $denom = $z0 + udistr( 1 - $position/100 );
1623             }
1624             my $nom     = 1 - $acc * $denom;
1625             my $lim = 100*uprob( - ( $z0 + $denom / $nom ) );
1626             push( @new_alphas, $lim );
1627           }
1628 #         print "@new_alphas\n";
1629           return \@new_alphas;
1630         }
1631
1632         my @limits = sort { $a <=> $b } keys %{$self -> {'confidence_limits'}};
1633         # Add the upper limits
1634         my $limnum = $#limits;
1635         for ( my $i = $limnum; $i >= 0; $i-- ) {
1636           $limits[$i] = $limits[$i]/2;
1637           push( @limits, 100-$limits[$i] );
1638         }
1639         my ( @bootstrap_array, @jackknife_array, @new_alphas, @z0, @acc );
1640         # Loop the estimates of the first (original) model
1641         for ( my $l = 0; $l < scalar @{$self -> {'bootstrap_estimates'}->
1642                                          [$model_number-1][0]}; $l++ ) {
1643           my ( @unsorted_array1, @unsorted_array2 );
1644           # Loop the bootstrap samples from 1 to get rid of original model
1645           for ( my $k = 1; $k < scalar @{$self -> {'bootstrap_estimates'}->
1646                                            [$model_number-1]}; $k++ ) {
1647             $unsorted_array1[$k-1] =
1648               $self -> {'bootstrap_estimates'}->[$model_number-1][$k][$l];
1649           }
1650           @{$bootstrap_array[$l]} = sort {$a <=> $b} @unsorted_array1;
1651
1652           # Loop the jackknife samples from 1 to get rid of original model
1653           for ( my $k = 1; $k < scalar @{$self -> {'jackknife_estimates'}->
1654                                            [$model_number-1]}; $k++ ) {
1655             $unsorted_array2[$k-1] =
1656               $self -> {'jackknife_estimates'}->[$model_number-1][$k][$l];
1657           }
1658           @{$jackknife_array[$l]} = sort {$a <=> $b} @unsorted_array2;
1659           $z0[$l]         = c_get_z0     ( $bootstrap_array[$l],
1660                                            $self -> {'bootstrap_estimates'} ->
1661                                            [$model_number-1][0][$l] );
1662           $acc[$l]        = c_get_acc    ( $jackknife_array[$l],
1663                                            $self -> {'jackknife_means'} ->
1664                                            [$model_number-1][0][$l] );
1665           $new_alphas[$l] = c_get_alphas ( \@limits, $acc[$l], $z0[$l] );
1666         }
1667         # Loop limits
1668         for ( my $lim_idx = 0; $lim_idx <= $#limits; $lim_idx++ ) {
1669           my @percentiles;
1670           # Loop parameters
1671           for ( my $l = 0; $l <= $#bootstrap_array; $l++ ) {
1672             my $limit = $new_alphas[$l][$lim_idx]/100;
1673             my $position = ( scalar @{$bootstrap_array[$l]} + 1 ) * $limit;
1674             my $percentile;
1675             if ( $position < 1 ) {
1676               $percentile = undef;
1677             } elsif ( $position > scalar @{$bootstrap_array[$l]} ) {
1678               $percentile = undef;
1679             } else {
1680               my ($int_med,$frac_med)   = split(/\./, $position );
1681               $frac_med = eval("0.".$frac_med);
1682               my $percentile_low  = $bootstrap_array[$l][ $int_med - 1];
1683               my $percentile_high = ( $bootstrap_array[$l][ $int_med ] -
1684                                       $bootstrap_array[$l][ $int_med - 1] ) * $frac_med;
1685               $percentile = $percentile_low + $percentile_high;
1686             }
1687             push( @percentiles, $percentile );
1688           }
1689           push( @{$self -> {'bca_confidence_intervals'} -> [$model_number-1]},
1690                 \@percentiles );
1691           push( @{$self -> {'bca_confidence_intervals_labels'}->[$model_number-1][0]},
1692                 $limits[$lim_idx].'%');
1693         }
1694         # Check the intervals
1695         for ( my $lim_idx = 0; $lim_idx <= $limnum; $lim_idx++ ) {
1696           my @within_ci;
1697           for ( my $l = 0; $l <= $#bootstrap_array; $l++ ) {
1698             my $lower_limit = $self -> {'bca_confidence_intervals'} ->
1699               [$model_number-1][$lim_idx][$l];
1700             my $upper_limit = $self -> {'bca_confidence_intervals'} ->
1701               [$model_number-1][($limnum*2+1)-$lim_idx][$l];
1702             if ( $self -> {'bca_confidence_intervals_check'} < $upper_limit and
1703                  $self -> {'bca_confidence_intervals_check'} > $lower_limit ) {
1704               push( @within_ci , 1 );
1705             } else {
1706               push( @within_ci , 0 );
1707             }
1708           }
1709           $self -> {'within_bca_confidence_intervals'} ->
1710             [$model_number-1]{$limits[$lim_idx]*2} = \@within_ci;
1711         }
1712         $self -> {'bca_confidence_intervals_labels'} -> [$model_number-1][1] =
1713           \@parameter_names;
1714       }
1715 end calculate_bca_confidence_intervals
1716
1717 # }}} calculate_bca_confidence_intervals
1718
1719 # {{{ calculate_percentile_confidence_intervals
1720
1721 start calculate_percentile_confidence_intervals
1722       {
1723         # Sort the limits from the inside out
1724         my @limits = sort { $b <=> $a } keys %{$self -> {'confidence_limits'}};
1725         foreach my $limit ( @limits ) {
1726           my ( @lower_limits, @upper_limits, @within_ci );
1727           # Loop the estimates of the first (original) model
1728           for ( my $l = 0; $l < scalar @{$self -> {'bootstrap_estimates'}->
1729                                            [$model_number-1][0]}; $l++ ) {
1730             my @parameter_array;
1731             # Loop the bootstrap samples from 1 to get rid of original model
1732             for ( my $k = 1; $k < scalar @{$self -> {'bootstrap_estimates'}->
1733                                              [$model_number-1]}; $k++ ) {
1734               my $val = $self -> {'bootstrap_estimates'}->[$model_number-1][$k][$l];
1735               # get rid of undefined values (these were probably deleted
1736               # when the bootstrap_estimates was created
1737               push( @parameter_array, $val ) if( defined $val );
1738             }
1739             my @sorted = sort {$a <=> $b} @parameter_array;
1740             for my $side ( 'lower', 'upper' ) {
1741               my $use_limit = $side eq 'lower' ? $limit/200 : 1-($limit/200);
1742               # percentile postition is:
1743               my $percentile_position = ( $#sorted + 2 ) * $use_limit;
1744               my $percentile;
1745               if ( $percentile_position < 1 ) {
1746                 $percentile = undef;
1747               } elsif ( $percentile_position > $#sorted +1) {
1748                 $percentile = undef;
1749               } else {
1750                 my ($int_med,$frac_med)   = split(/\./, $percentile_position );
1751                 $frac_med = eval("0.".$frac_med);
1752                 my $percentile_low  = $sorted[ $int_med - 1];
1753                 my $percentile_high = ( $sorted[ $int_med ] - $sorted[ $int_med - 1] ) * $frac_med;
1754                 $percentile = $percentile_low + $percentile_high;
1755               }
1756               push( @lower_limits, $percentile ) if ( $side eq 'lower' );
1757               push( @upper_limits, $percentile ) if ( $side eq 'upper' );
1758             }
1759             if ( $self -> {'percentile_confidence_intervals_check'} < $upper_limits[$#upper_limits] and
1760                  $self -> {'percentile_confidence_intervals_check'} > $lower_limits[$#lower_limits] ) {
1761               push( @within_ci , 1 );
1762             } else {
1763               push( @within_ci , 0 );
1764             }
1765           }
1766           unshift( @{$self -> {'percentile_confidence_intervals'} ->
1767                        [$model_number-1]}, \@lower_limits );
1768           push( @{$self -> {'percentile_confidence_intervals'} ->
1769                     [$model_number-1]}, \@upper_limits );
1770           unshift( @{$self -> {'percentile_confidence_intervals_labels'}->
1771                        [$model_number-1][0]}, ($limit/2).'%' );
1772           push( @{$self -> {'percentile_confidence_intervals_labels'}->
1773                     [$model_number-1][0]},(100-($limit/2)).'%');
1774           $self -> {'within_percentiles'}->[$model_number-1]{$limit}=\@within_ci;
1775         }
1776         $self -> {'percentile_confidence_intervals_labels'} ->
1777           [$model_number-1][1] = \@parameter_names;
1778       }
1779 end calculate_percentile_confidence_intervals
1780
1781 # }}} calculate_percentile_confidence_intervals
1782
1783 # {{{ modelfit_post_fork_analyze
1784
1785 start modelfit_post_fork_analyze
1786 end modelfit_post_fork_analyze
1787
1788 # }}}
1789
1790 # {{{ resample
1791
1792 start resample
1793       {
1794         my $dataObj = $model -> datas -> [0];
1795         for( my $i = 1; $i <= $self -> {'samples'}; $i++ ) {
1796           my ($bs_dir, $bs_name) = OSspecific::absolute_path( $self -> {'directory'}, "bs$i.dta" );
1797           my $new_name = $bs_dir . $bs_name;
1798           my $boot_sample = $dataObj -> resample( 'subjects' => $self -> {'subjects'},
1799                                                   'new_name' => $new_name,
1800                                                   'target'   => $target );
1801           my $newmodel = $model -> copy( filename => "bs$i.mod",
1802                                          target   => $target,
1803                                          ignore_missing_files => 1 );
1804           $newmodel -> datafiles( new_names => ["bs$i.dta"] );
1805           $newmodel -> datas -> [0] = $boot_sample ;
1806           $newmodel -> write;
1807           push( @resample_models, $newmodel );
1808         }
1809       }
1810 end resample
1811
1812 # }}} resample
1813
1814 # {{{ _sampleTools
1815
1816 start _sampleTools
1817       {
1818         foreach my $tool ( @{$self -> {'tools'}} ) {
1819           my @models = @{$tool -> models};
1820           foreach my $model (@models){
1821             my $dataObj = $model -> datas -> [0];
1822             for( my $i = 1; $i <= $samples; $i++ ) {
1823               my $boot_sample = $dataObj -> resample( 'subjects' => $self -> {'subjects'},
1824                                                       'new_name' => "bs$i.dta",
1825                                                       'target' => $target );
1826               my $newmodel;
1827               $newmodel = $model -> copy( filename => "bs$i.mod" );
1828               $newmodel -> datafiles( new_names => ["bs$i.dta"] );
1829               $newmodel -> datas -> [0] = $boot_sample ;
1830               $newmodel -> write;
1831               if( defined( $tool -> models ) ){
1832                 push( @{$tool -> models}, $newmodel );
1833               } else {
1834                 $tool -> models( [ $newmodel ] );
1835               }
1836             }
1837           }
1838         }
1839       }
1840 end _sampleTools
1841
1842 # }}} _sampleTools
1843
1844 # {{{ print_results
1845
1846 start print_results
1847       {
1848         # Run the print_results specific for the subtool
1849         my $sub_print_results = $self -> {'subtools'} -> [0];
1850         if ( defined $sub_print_results ) {
1851           sub get_dim {
1852             my $arr      = shift;
1853             my $dim      = shift;
1854             my $size_ref = shift;
1855             $dim++;
1856             if ( defined $arr and ref($arr) eq 'ARRAY' ) {
1857               push( @{$size_ref}, scalar @{$arr} );
1858               ( $dim, $size_ref ) = get_dim( $arr->[0], $dim, $size_ref );
1859             }
1860             return ( $dim, $size_ref );
1861           }
1862           sub format_value {
1863             my $val = shift;
1864             if ( not defined $val or $val eq '' ) {
1865               return sprintf("%10s",$PsN::output_style).',';
1866             } else {
1867               $_ = $val;
1868               my $nodot = /.*\..*/ ? 0 : 1;
1869               $_ =~ s/\.//g;
1870               if ( /.*\D+.*/ or $nodot) {
1871                 return sprintf("%10s",$val).',';
1872               } else {
1873                 return sprintf("%10.5f",$val).',';
1874               }
1875             }
1876           }
1877           debug -> die( message => "No results_file defined" )
1878             unless ( defined $self -> {'results_file'} );
1879           open ( RES, ">".$self -> {'directory'}.'/'.$self -> {'results_file'} );
1880           if ( defined $self -> {'results'} ) {
1881             my @all_results = @{$self -> {'results'}};
1882             for ( my $i = 0; $i <= $#all_results; $i++ ) {
1883               if ( defined $all_results[$i]{'own'} ) {
1884                 my @my_results = @{$all_results[$i]{'own'}};
1885                 for ( my $j = 0; $j <= $#my_results; $j++ ) {
1886                   # These size estimates include the problem and sub_problem dimensions:
1887                   my ( $ldim, $lsize_ref ) = get_dim( $my_results[$j]{'labels'}, -1, [] );
1888                   my ( $vdim, $vsize_ref ) = get_dim( $my_results[$j]{'values'}, -1, [] );
1889                   print RES $my_results[$j]{'name'},"\n" if ( $vdim > 1 );
1890                   if ( defined $my_results[$j]{'values'} and
1891                        scalar @{$my_results[$j]{'values'}} >= 0 ) {
1892                     my @values  = @{$my_results[$j]{'values'}};
1893                     my @labels;
1894                     if ( defined $my_results[$j]{'labels'} and
1895                          scalar @{$my_results[$j]{'labels'}} >= 0 ) {
1896                       @labels = @{$my_results[$j]{'labels'}};
1897                     }
1898
1899                     # Print Header Labels
1900                     if ( $ldim == 0 ) {
1901                       my $label = \@labels;
1902                       print RES ','.format_value($label),"\n";
1903                     } elsif ( $ldim == 2 and defined $labels[1] ) {
1904                       print RES ',';
1905                       for ( my $n = 0; $n < scalar @{$labels[1]}; $n++ ) {
1906                         my $label = $labels[1][$n];
1907                         print RES format_value($label);
1908                       }
1909                       print RES "\n" if ( scalar @{$labels[1]} );
1910                     }
1911                     # Print the values:
1912                     if ( $vdim == 0 ) {
1913                       print RES ','.format_value(\@values),"\n";
1914                     } elsif ( $vdim == 1 ) {
1915                       for ( my $m = 0; $m < scalar @{\@values}; $m++ ) {
1916                         my $label = $labels[$m];
1917                         print RES ','.format_value($label);
1918                         my $val = $values[$m];
1919                         print RES ','.format_value($val),"\n";
1920                       }
1921                     } elsif ( $vdim == 2 ) {
1922                       for ( my $m = 0; $m < scalar @{\@values}; $m++ ) {
1923                         my $label;
1924                         if ( $ldim == 1 ) {
1925                           $label = $labels[$m];
1926                         } elsif ( $ldim == 2 ) {
1927                           $label = $labels[0][$m];
1928                         }
1929                         print RES format_value($label);
1930                         if ( defined $values[$m] ) {
1931                           for ( my $n = 0; $n < scalar @{$values[$m]}; $n++ ) {
1932                             print RES format_value($values[$m][$n]);
1933                           }
1934                         }
1935                         print RES "\n";
1936                       }
1937                     }
1938                   }
1939                 }
1940               }
1941             }
1942           }
1943           close( RES );
1944         } else {
1945           debug -> warn( level   => 2,
1946                          message => "No subtools defined".
1947                          ", using default printing routine" );
1948         }
1949       }
1950 end print_results
1951
1952 # }}}
1953
1954 # {{{ create_matlab_scripts
1955
1956 start create_matlab_scripts
1957     {
1958       if( defined $PsN::lib_dir ){
1959         unless( -e $PsN::lib_dir . '/histograms.m' and
1960                 -e $PsN::lib_dir . '/bca.m' ){
1961           'debug' -> die( message => 'Bootstrap matlab template scripts are not installed, no matlab scripts will be generated.' );
1962           return;
1963         }
1964
1965         open( PROF, $PsN::lib_dir . '/histograms.m' );
1966         my @file = <PROF>;
1967         close( PROF );
1968         my $found_code;
1969         my $code_area_start=0;
1970         my $code_area_end=0;
1971
1972
1973         for(my $i = 0;$i < scalar(@file); $i++) {
1974           if( $file[$i] =~ /% ---------- Autogenerated code below ----------/ ){
1975             $found_code = 1;
1976             $code_area_start = $i;
1977           }
1978           if( $file[$i] =~ /% ---------- End autogenerated code ----------/ ){
1979             unless( $found_code ){
1980               'debug' -> die ( message => 'Bootstrap matlab template script is malformated, no matlab scripts will be generated' );
1981               return;
1982             }
1983             $code_area_end = $i;
1984           }
1985         }
1986
1987         my @auto_code;
1988         if( $self -> {'type'} eq 'bca' ){
1989           push( @auto_code, "use_bca = 1;                               % Was a BCa-type of\n" );
1990         } else {
1991           push( @auto_code, "use_bca = 0;                               % Was a BCa-type of\n" );
1992         }
1993
1994         push( @auto_code, "                                % bootstrap run?\n" );
1995         if( ref $self -> {'samples'} eq 'ARRAY' ) {
1996           push( @auto_code, "bs_samples = ".$self -> {'samples'}->[0][0].";                     % Number of bootstrap samples\n" );
1997         } else {
1998           push( @auto_code, "bs_samples = ".$self -> {'samples'}.";                     % Number of bootstrap samples\n" );
1999         }
2000         if( $self -> {'type'} eq 'bca' ){
2001           my $ninds = $self -> models -> [0]
2002               -> datas -> [0] -> count_ind;
2003           push( @auto_code, "jk_samples = $ninds;                       % Number of (BCa) jackknife samples\n\n" );
2004         }
2005
2006         push( @auto_code, "col_names = { 'Significant Digits',\n" );
2007         push( @auto_code, "              'Condition Number',\n" );
2008         push( @auto_code, "              'OFV',\n" );
2009
2010         my $nps = $self -> {'models'} -> [0] -> nomegas -> [0];
2011
2012         my %param_names;
2013         my( @par_names, @se_names, @np_names, @sh_names );
2014         foreach my $param ( 'theta','omega','sigma' ) {
2015           my $labels = $self -> {'models'} -> [0] -> labels( parameter_type => $param );
2016           if ( defined $labels ){
2017             foreach my $label ( @{$labels -> [0]} ){
2018               push( @par_names, "                '",$label,"',\n" );
2019               push( @se_names, "                 '",'se-'.$label,"',\n" );
2020             }
2021           }
2022         }
2023
2024         for( my $i = 1; $i <= ($nps*($nps+1)/2); $i++ ) {
2025           push( @np_names, "             '",'np-om'.$i,"',\n" );
2026         }
2027
2028         for( my $i = 1; $i <= $nps; $i++ ) {
2029           push( @sh_names, "             '",'shrinkage-eta'.$i,"',\n" );
2030         }
2031
2032         push( @sh_names, "               '",'shrinkage-iwres',"'\n" );
2033
2034 # NP not used for now
2035
2036         push( @auto_code,(@par_names, @se_names, @sh_names));
2037 #       push( @auto_code,(@par_names, @se_names, @np_names, @sh_names));
2038         push( @auto_code, "           };\n\n" );
2039
2040         my @np_columns = (0) x ($nps*($nps+1)/2);
2041         my @sh_columns = (0) x ($nps+1);
2042
2043         if( $self -> {'type'} eq 'bca' ){
2044           push( @auto_code, "fixed_columns = [ 0, 0, 0, " );
2045         } else {
2046           push( @auto_code, "fixed_columns = [ 0, 0, 0, " );
2047         }
2048         my ( @fixed_columns, @same_columns, @adjust_axes );
2049         foreach my $param ( 'theta','omega','sigma' ) {
2050           my $fixed = $self -> {'models'} -> [0] -> fixed( parameter_type => $param );
2051
2052           if ( defined $fixed ){
2053             push( @fixed_columns, @{$fixed -> [0]} );
2054             if( $param eq 'theta' ) {
2055               push( @same_columns, (0) x scalar( @{$fixed -> [0]} ) );
2056             }
2057           }
2058         }
2059
2060         @adjust_axes = (1) x ( ($#fixed_columns + 1) * 2 +
2061                                $#sh_columns + 1 );
2062 #                              $#np_columns + $#sh_columns + 2 );
2063
2064         push( @auto_code , join( ', ' , @fixed_columns).', '.
2065               join( ', ' , @fixed_columns).', '.
2066 #             join( ', ' , @np_columns).', '.
2067               join( ', ' , @sh_columns)."];\n\n" );
2068
2069         if( $self -> {'type'} eq 'bca' ){
2070           push( @auto_code, "same_columns  = [ 0, 0, 0, " );
2071         } else {
2072           push( @auto_code, "same_columns  = [ 0, 0, 0, " );
2073         }
2074         foreach my $param ( 'omegas','sigmas' ) {
2075           my $parameters = $self -> {'models'} -> [0] -> problems -> [0] -> $param;
2076           foreach my $parameter ( @{$parameters} ){
2077             if( $parameter -> same() ){
2078               push( @same_columns, (1) x $parameter -> size() );
2079             } else {
2080               push( @same_columns, (0) x scalar @{$parameter -> options} );
2081             }
2082           }
2083         }
2084         push( @auto_code , join( ', ' , @same_columns ).', '.
2085               join( ', ' , @same_columns).', '.
2086 #             join( ', ' , @np_columns).', '.
2087               join( ', ' , @sh_columns)."];\n\n" );
2088
2089         push( @auto_code , "adjust_axes   = [ 1, 1, 1, ".join( ', ' , @adjust_axes)."];\n\n" );
2090
2091 #        push( @auto_code , 'npomegas = '.($nps*($nps+1)/2).";\n\n" );
2092         push( @auto_code , "npomegas = 0;\n\n" );
2093
2094
2095         push( @auto_code, "minimization_successful_col    = 5;  % Column number for the\n" );
2096         push( @auto_code, "                                                     % minimization sucessful flag\n" );
2097         push( @auto_code, "covariance_step_successful_col = 6;  % As above for cov-step warnings\n" );
2098         push( @auto_code, "covariance_step_warnings_col   = 7;  % etc\n" );
2099         push( @auto_code, "estimate_near_boundary_col     = 8;  % etc\n" );
2100
2101         push( @auto_code, "not_data_cols = 13;                     % Number of columns in the\n" );
2102         push( @auto_code, "                                        % beginning that are not\n" );
2103         push( @auto_code, "                                        % parameter estimates.\n" );
2104
2105         push( @auto_code, "filename = 'raw_results_matlab.csv';\n" );
2106
2107         splice( @file, $code_area_start, ($code_area_end - $code_area_start), @auto_code );
2108         open( OUTFILE, ">", $self -> {'directory'} . "/histograms.m" );
2109         print OUTFILE "addpath " . $PsN::lib_dir . ";\n";
2110         print OUTFILE @file ;
2111         close OUTFILE;
2112
2113         open( OUTFILE, ">", $self -> {'directory'} . "/raw_results_matlab.csv" );
2114         for( my $i = 0; $i < scalar ( @{$self -> {'raw_results'} -> [0]} ); $i ++ ){
2115 #         $self -> {'raw_results'} -> [0] -> [$i][0] =
2116 #             $self -> {'raw_results'} -> [0] -> [$i][0] eq 'bootstrap' ?
2117 #             1 : $self -> {'raw_results'} -> [0] -> [$i][0];
2118 #         $self -> {'raw_results'} -> [0] -> [$i][0] =
2119 #             $self -> {'raw_results'} -> [0] -> [$i][0] eq 'jackknife' ?
2120 #             2 : $self -> {'raw_results'} -> [0] -> [$i][0];
2121           map( $_ = $_ eq 'NA' ? 'NaN' : $_, @{$self -> {'raw_results'} -> [0] -> [$i]} );
2122           map( $_ = not( defined $_ ) ? 'NaN' : $_, @{$self -> {'raw_results'} -> [0] -> [$i]} );
2123           print OUTFILE join( ',', @{$self -> {'raw_results'} -> [0] -> [$i]} ), "\n";
2124         }
2125         close OUTFILE;
2126
2127       } else {
2128         'debug' -> die( message => 'matlab_dir not configured, no matlab scripts will be generated.');
2129         return;
2130       }
2131     }
2132 end create_matlab_scripts
2133
2134 # }}}
2135
2136 # {{{ create_R_scripts
2137 start create_R_scripts
2138     {
2139       unless( -e $PsN::lib_dir . '/R-scripts/bootstrap.R' ){
2140         'debug' -> die( message => 'Bootstrap R-script are not installed, no R-script will be generated.' );
2141         return;
2142       }
2143       cp ( $PsN::lib_dir . '/R-scripts/bootstrap.R', $self -> {'directory'} );
2144     }
2145 end create_R_scripts
2146 # }}}