bin/extended_grid

   1 #!/usr/local/bin/perl
   2
   3 # Only for Development
   4 use FindBin qw($Bin);
   5 use lib "$Bin/../sfPsN/lib";
   6
   7 # Don't edit the line below, it must look exactly like this.
   8 # Everything above this line will be replaced #
   9
  10 # Perl includes #
  11 use strict;
  12 use Getopt::Long;
  13 # External modules #
  14 use Math::Random;
  15 # PsN includes #
  16 use PsN;
  17 use tool::modelfit;
  18 use model;
  19 use debug;
  20 use common_options;
  21 use ui;
  22 use File::Copy 'cp';
  23 use Data::Dumper;
  24
  25 ################################################################################
  26 #
  27 # This is the extended_gri. First the command line options are parsed
  28 # using Getopt::Long. PsN "common options" are used to fine tune the
  29 # running.
  30 #
  31 ################################################################################
  32
  33 # {{{ options handling
  34
  35 my $cmd_line = $0 . " " . join( " ", @ARGV );
  36
  37 my %options;
  38
  39 my %optional_options = ( 'ids_to_simulate:i' => undef,
  40                          'id_column:s' => undef,
  41                          'dv_column:s' => undef );
  42
  43 my $res = GetOptions( \%options,
  44                       @common_options::get_opt_strings,
  45                       keys(%optional_options) );
  46
  47 exit unless $res;
  48
  49 common_options::set_globals( \%options );
  50 common_options::get_defaults( \%options, 'extended_grid' );
  51 my %help_text;
  52 $help_text{Pre_help_message} = <<'EOF';
  53     <h3 class="heading1">extended_grid</h3>
  54
  55       Perl script for getting standard errors of etas.
  56
  57     <h3 class="heading1">Usage:</h3>
  58 EOF
  59
  60     $help_text{Description} = <<'EOF';
  61    <h3 class="heading1">Description:</h3>
  62
  63       extended_grid requires "MAXEVALS" to be abreviated to "MAX" if it
  64       exists in the model file.
  65
  66 EOF
  67     $help_text{Examples} = <<'EOF';
  68    <h3 class="heading1">Example:</h3>
  69
  70       <p class="style2">extended_grid pheno.mod</p>
  71
  72        This will create a new directory extended_grid_dirX where X is a
  73        number increased every time you run the program. Inside that
  74        directory it will create a new model file for each OMEGA in the
  75        file. You will be able to retrieve the results from the
  76        corresponding table file.
  77 EOF
  78
  79 common_options::online_help('extended_grid',\%options, \%help_text, {}, {});
  80
  81 if ( scalar( @ARGV ) < 1 ) {
  82   print "A model file must be specified. Use '$0 -h' for help.\n";
  83   exit;
  84 }
  85
  86 unless( defined $options{'seed'} ){
  87   $options{'seed'} = 1;
  88 }
  89
  90 unless( defined $options{'dv_column'} ){
  91   $options{'dv_column'} = 'DV';
  92 }
  93
  94 unless( defined $options{'id_column'} ){
  95   $options{'id_column'} = 'ID';
  96 }
  97
  98 my( $base_dir, $dummy) = OSspecific::absolute_path( '', '' );
  99
 100 my $directory;
 101
 102 if ( defined $options{'directory'} ) {
 103   my $dummy;
 104   ( $directory, $dummy ) = OSspecific::absolute_path( $options{'directory'} );
 105 } else {
 106   $directory = OSspecific::unique_path( 'extended_grid_dir' ,
 107                                         $base_dir );
 108 }
 109
 110 ui -> category( 'modelfit' );
 111
 112 random_set_seed_from_phrase( $options{'seed'} ) if ( defined $options{'seed'} );
 113
 114 debug -> level( $options{'debug'} );
 115 debug -> package( $options{'debug_package'} );
 116 debug -> subroutine( $options{'debug_subroutine'} );
 117 debug -> warn_with_trace( $options{'warn_with_trace'} );
 118
 119 my $eval_string = common_options::model_parameters(\%options);
 120
 121 unless( -d $directory ){
 122   mkdir( $directory );
 123
 124   open(CMD, ">", $directory . "/command.txt");
 125   print CMD $cmd_line, "\n";
 126   close(CMD);
 127 }
 128
 129 # }}}
 130
 131 ################################################################################
 132 #
 133 # The original model file is run and simulation started.
 134 #
 135 ################################################################################
 136
 137 # {{{ Orig + simulation
 138
 139 my $model;
 140
 141 $model = model -> new ( eval( $eval_string ),
 142                         filename                    => $ARGV[0],
 143                         ignore_missing_output_files => 1 );
 144
 145 my $datafile_name = $model -> datas -> [0] -> full_name;
 146
 147 chdir( $directory );
 148
 149 my $modelfit;
 150
 151 $modelfit = tool::modelfit ->
 152     new ( eval( $common_options::parameters ),
 153           directory => 'orig_modelfit_dir',
 154           models => [$model] );
 155
 156 $modelfit -> run;
 157
 158 my $model_copy = $model -> copy( filename => 'simulation.mod',
 159                                  copy_data => 1,
 160                                  data_file_names => ['simulation.dta']);
 161
 162 $model_copy -> update_inits( from_output => $model -> outputs -> [0] );
 163
 164 $model_copy -> remove_records( type => 'estimation' );
 165 $model_copy -> remove_records( type => 'covariance' );
 166 $model_copy -> remove_records( type => 'table' );
 167
 168 $model_copy -> set_records( type => 'simulation',
 169                             record_strings => ['('.$options{'seed'}.') ONLYSIM']);
 170
 171 my $nr_of_old_etas = $model -> nomegas -> [0];
 172 my $nr_of_old_sigmas = $model -> nsigmas -> [0];
 173 my @table_string;
 174 for( 1..$nr_of_old_etas ){
 175
 176   push( @table_string, "ETA$_" );
 177
 178 }
 179
 180 $model_copy -> set_records( type => 'table',
 181                             record_strings => [ "ID ",join( " ",@table_string ), " FIRSTONLY FILE=sim" ] );
 182
 183 $model_copy -> _write;
 184
 185 my $id_to_sim;
 186
 187 if( defined $options{'ids_to_simulate'} ){
 188   $id_to_sim = $options{'ids_to_simulate'};
 189
 190   my $data_set = $model_copy -> datas -> [0];
 191
 192   if( $id_to_sim <= scalar @{$data_set -> individuals} ){
 193     my $individuals = $data_set -> individuals;
 194     $data_set -> individuals( \@{$individuals}[0..$id_to_sim] );
 195   } else {
 196
 197     my $nr_of_new_ids = $id_to_sim - scalar @{$data_set -> individuals};
 198
 199     for( 1..$nr_of_new_ids ){
 200       my $new_individual = $data_set -> individuals -> copy;
 201
 202       push( @{$data_set -> individuals}, $new_individual );
 203
 204     }
 205
 206   }
 207
 208   $data_set -> _write;
 209
 210 } else {
 211   $id_to_sim = scalar @{$model_copy -> datas -> [0] -> individuals};
 212 }
 213
 214 $modelfit = tool::modelfit ->
 215     new ( eval( $common_options::parameters ),
 216           directory => 'simulation_dir',
 217           models => [$model_copy] );
 218
 219 $modelfit -> run;
 220
 221 # }}}
 222
 223 # {{{ create Dataset
 224
 225 my $extended_model = $model -> copy( filename => 'extended.mod',
 226                                      copy_data => 1,
 227                                      data_file_names => ['extended.dta']);
 228
 229 my $data_copy = $extended_model -> datas -> [0];
 230
 231 $extended_model -> update_inits( from_output => $model -> outputs -> [0] );
 232
 233 my $ids = $data_copy -> individuals;
 234
 235 #my $omegas = $models[0] -> problems -> [0] -> omegas();
 236 #my @L2_column;
 237 #my $L2_value = 0;
 238
 239 #foreach my $omega( @{$omegas} ){
 240 #  $L2_value++;
 241 #  for( 1..$omega -> size() ){
 242 #    push( @L2_column, $L2_value );
 243 #  }
 244 #}
 245
 246 foreach my $id( @{$ids} ){
 247   my @last_row_of_id = split(/,/,$id -> subject_data -> [$#{$id -> subject_data}]);
 248
 249   foreach my $row( @{$id -> subject_data} ){
 250     $row .= ",0";
 251   }
 252 }
 253
 254 my $id_template = $ids->[$#{$ids}];
 255
 256 my @last_row_of_id = split(/,/,$id_template -> subject_data -> [$#{$id_template -> subject_data}]);
 257 my $row_length = scalar @last_row_of_id;
 258
 259 my $id_column = $data_copy -> column_head_indices -> {$options{'id_column'}};
 260 my $dv_column = $data_copy -> column_head_indices -> {$options{'dv_column'}};
 261 my $evid_column = $data_copy -> column_head_indices -> {'EVID'};
 262 my $last_id = $id_template -> idnumber();
 263
 264 my $table_data = data -> new( filename => 'sim' );
 265
 266 my $table_ids = $table_data -> individuals();
 267
 268 foreach my $id_num ( 1..$id_to_sim ){
 269   my @new_id_rows;
 270   my @row_template = ('.') x $row_length;
 271   my $type = 1;
 272
 273   my $eta_from_sim = $table_ids -> [$id_num-1] -> subject_data -> [0];
 274
 275   foreach my $eta_num ( 1..$nr_of_old_etas ){
 276     my @eta_from_sim = split( /,/ , $eta_from_sim  );
 277     my @new_row = @row_template;
 278     $new_row[$id_column-1] = ($last_id + $id_num );
 279     $new_row[$dv_column-1] = $eta_from_sim[$eta_num]; # eta begins with 1, but there is an ID-column in $eta_from_sim
 280     $new_row[$evid_column-1] = 0;
 281     $new_row[$#new_row] = $type++;
 282     push( @new_id_rows, join(',',@new_row ) );
 283   }
 284   my $new_id = data::individual -> new( idcolumn => $last_id + $id_num,
 285                                         subject_data => \@new_id_rows,
 286                                         data_id => $data_copy -> {'data_id'} );
 287   push( @{$data_copy -> individuals}, $new_id );
 288 }
 289
 290 push( @{$data_copy -> {'header'}}, 'TYPE' );
 291 $data_copy -> synchronize;
 292
 293 $data_copy -> _write;
 294
 295 # }}}
 296
 297 # {{{ create the second model
 298
 299 # Get error block to replace EPS(n) with ETA(n + $nr_of_old_eta) and
 300 # replace ERR(n) with ERR(n+$nr_of_old_eta)
 301
 302 my $error = $extended_model -> record( record_name => 'error' );
 303
 304 my %seen_eps;
 305 my %seen_err;
 306
 307 my $error_has_if = 0;
 308
 309 my $loop_end = $#{$error -> [0]};
 310 for( my $row = 0 ; $row <= $loop_end; $row ++ ){
 311   $error -> [0] -> [$row] =~ s/
 312                               (.)EPS\((\d+)
 313                               (?{$seen_eps{$2}=1;$dummy=$2+$nr_of_old_etas;})\)
 314                               /$1ETA($dummy)/gx;
 315
 316   if( $error -> [0] -> [$row] =~ /IF\s*\(/ ){
 317     $error_has_if = 1;
 318   }
 319
 320 #  $error -> [0] -> [$row] =~ s/
 321 #                              (.)ERR\((\d+)
 322 #                             (?{$seen_err{$2}=1;$dummy=$2+$nr_of_old_etas;})\)
 323 #                             /$1ETA($dummy)/gx;
 324
 325
 326 #  if( $error -> [0] -> [$row] =~ /(.*)Y(\s*)=(.*)/ ){
 327 #
 328 #
 329 #    @{$error -> [0]} = ( @{$error -> [0]}[0..($row-1)],
 330 #                        '       WT=1/100000',
 331 #                        '       IF( TYPE.LE.0 ) ' . $error -> [0] -> [$row],
 332 #                        @{$error -> [0]}[$row+1..$#{$error->[0]}] );
 333 #    $row++;
 334 #
 335 #  }
 336
 337 }
 338
 339 unless( $error_has_if ){
 340   unshift( @{$error->[0]}, 'IF(TYPE.LE.0) THEN' );
 341   push( @{$error -> [0]}, 'ENDIF' );
 342 }
 343
 344 my $nr_of_added_etas = scalar keys %seen_eps;
 345
 346 my $input = $extended_model -> record( record_name => 'input' );
 347
 348 $input->[0]->[0] =~ s/\$INPUT//;
 349 chomp( $input->[0]->[ $#{$input->[0]} ] );
 350 $input->[0]->[ $#{$input->[0]} ] .= ' TYPE';
 351
 352 $extended_model -> set_records( type => 'input',
 353                                 record_strings =>  $input->[0]);
 354
 355
 356 foreach( 1..$nr_of_old_etas ){
 357   push(@{$error -> [0]}, '       IF( TYPE.EQ.'.($_).' ) Y = ETA('.($_).')+ERR(2)*W/100000');
 358 }
 359
 360
 361 $extended_model -> set_records( type => 'error',
 362                                 record_strings => $error -> [0] );
 363
 364 $extended_model -> initial_values( parameter_type => 'sigma',
 365                                    add_if_absent => 1,
 366                                    parameter_numbers => [[$nr_of_old_sigmas+1]],
 367                                    new_values => [[1]] );
 368
 369 $extended_model -> initial_values( parameter_type => 'omega',
 370                                    add_if_absent => 1,
 371                                    parameter_numbers => [[$nr_of_old_etas+1 .. $nr_of_old_etas+$nr_of_added_etas]],
 372                                    new_values => [[(1) x $nr_of_added_etas]] );
 373
 374 $extended_model -> fixed( parameter_type => 'sigma',
 375                           parameter_numbers => [[$nr_of_old_sigmas+1]],
 376                           new_values => [[1]] );
 377
 378 $extended_model -> set_records( type => 'estimation',
 379                                 record_strings => ['MAX=0 POSTHOC'] );
 380
 381 $extended_model -> set_records( type => 'nonparametric',
 382                                 record_strings => ['MARGINALS MSFO=MSF01 UNCONDITIONAL'] );
 383
 384 $extended_model -> remove_records( type => 'covariance' );
 385
 386 my $code_block;
 387 my $use_pk = 0;
 388 # Find PK or PRED block
 389 if( defined $extended_model -> pk ){
 390   $code_block = $extended_model -> pk;
 391   $use_pk = 1;
 392 } elsif( defined $extended_model -> pred ){
 393   $code_block = $extended_model -> pred;
 394 } else {
 395   die "Error: No \$PK or \$PRED found... \n";
 396 }
 397
 398 push( @{$code_block},'   JD = DEN_');
 399
 400 my @JD_string;
 401
 402 for( 1..$nr_of_old_etas ){
 403   push(@{$code_block},"   DN$_=CDEN_($_)" );
 404   push( @JD_string, "DN$_" );
 405 }
 406
 407 if( $use_pk ){
 408   $extended_model -> pk( new_pk => $code_block );
 409 } else {
 410   $extended_model -> pred( new_pred => $code_block );
 411 }
 412
 413 $extended_model -> set_records( type => 'table',
 414                                 record_strings => ['ID '.join( ' ', @table_string, 'JD', @JD_string,'NOPRINT NOAPPEND FIRSTONLY FILE=nptab' ) ] );
 415
 416 $extended_model -> extra_output( ['nptab'] );
 417
 418 $extended_model -> _write;
 419
 420 $modelfit = tool::modelfit ->
 421     new ( eval( $common_options::parameters ),
 422           prepend_model_file_name => 1,
 423           directory => 'extended_dir',
 424           models => [$extended_model] );
 425
 426 $modelfit -> run;
 427
 428 # }}}
 429
 430 # {{{ Results handling (Cut'n'Paste from nonp_bootstrap_v2)
 431
 432 ################################################################################
 433 #
 434 # Do one run per individual initial clearence and get fort.80 from
 435 # each run.
 436 #
 437 ################################################################################
 438
 439 # {{{ Fort.80 creation
 440
 441 my @parameter_P_sums;
 442 my @result_header;
 443 my $batch_nr = 0;
 444
 445 # {{{ create iofvcont.f
 446
 447 unless( -e 'iofvcont.f' ){
 448   open(IOFVCONT , '>', 'iofvcont.f');
 449
 450   print IOFVCONT <<'EOF';
 451       subroutine contr (icall,cnt,ier1,ier2)
 452       parameter (no=1000)
 453       common /rocm1/ y(no),data(no,3),nobs
 454       integer nobs, un
 455       double precision cnt,y
 456       DATA un /80/
 457       if (icall.le.1) return
 458       call ncontr (cnt,ier1,ier2,l2r)
 459 C     individual obj. funct. value for indiv. jj = cnt
 460       write(un,10) data(1,1),cnt
 461    10 FORMAT(1F6.2,1F19.4)
 462       return
 463       end
 464 EOF
 465    close( IOFVCONT );
 466 }
 467
 468 # }}}
 469
 470 # TODO update all initial values
 471
 472 my @new_ofv_models;
 473
 474 unless ( -e "ofv_dir" ){
 475   mkdir( "ofv_dir" );
 476 }
 477 chdir( "ofv_dir" );
 478 cp( '../iofvcont.f', './iofvcont.f' ) or print "Warning: unable to copy iofvcont.f\n";
 479
 480 my $table = table_file -> new( filename => "extended.nptab",
 481                                directory => '../' );
 482
 483 my $nomegas = $extended_model -> nomegas -> [0];
 484
 485 my @individual_etas;
 486 my @individual_etas_probabilities;
 487
 488 for( 0..($nomegas-1) ){
 489
 490   my $individual_eta = $table -> column_to_array( column => $_ + 1 ); # plus 1 is for ID column
 491   my $individual_eta_probabilities = $table -> column_to_array( column => ($nomegas-1) + $_ + 1 + 1 ); # first plus 1 is for ID column. Second is for JD
 492
 493   push( @individual_etas , $individual_eta );
 494   push( @individual_etas_probabilities , $individual_eta_probabilities );
 495
 496 }
 497
 498 my $number_of_individuals = @{$individual_etas[0]};
 499
 500 foreach my $id( 0..($number_of_individuals-1) ){
 501
 502   my $copy = $extended_model -> copy( filename => 'ofv_model_'.($id+1).'.mod' );
 503
 504   # Not needed ?
 505   # $copy -> update_inits( from_output => $extended_model -> outputs -> [0] );
 506
 507   my @new_values;
 508
 509   for my $parm ( 0..$nomegas-1 ){
 510
 511     push( @new_values, $individual_etas[$parm][$id] );
 512
 513   }
 514
 515 #  $copy -> lower_bounds( parameter_type => 'theta',
 516 #                        parameter_numbers => [[1..$nomegas]],
 517 #                        new_values => [[(-10000) x $nomegas]]) ;
 518
 519 #  $copy -> upper_bounds( parameter_type => 'theta',
 520 #                        parameter_numbers => [[1..$nomegas]],
 521 #                        new_values => [[(10000) x $nomegas]]) ;
 522
 523 #  $copy -> initial_values( parameter_type => 'theta',
 524 #                          parameter_numbers => [[1..$nomegas]],
 525 #                          new_values => [\@new_values]) ;
 526
 527   $copy -> add_records( type => 'contr',
 528                         record_strings => ['DATA=(ID)'] );
 529
 530   $copy -> add_option( record_name => 'subroutine',
 531                        option_name => 'CONTR',
 532                        option_value => 'iofvcont.f' );
 533
 534   #$copy -> add_option( record_name => 'DATA',
 535 #                      option_name => 'ACCEPT',
 536 #                      option_value => "ID.LE.$ids" );
 537
 538   $copy -> lower_bounds( parameter_type => 'omega',
 539                          parameter_numbers => [[1..$copy -> nomegas -> [0]]],
 540                          new_values => [[(-1) x $copy -> nomegas -> [0]]] );
 541
 542   $copy -> initial_values( parameter_type => 'omega',
 543                            parameter_numbers => [[1..$copy -> nomegas -> [0]]],
 544                            new_values => [[(0) x $copy -> nomegas -> [0]]] );
 545
 546 #  $copy -> initial_values( parameter_type => 'omega',
 547 #                          parameter_numbers => [[1..$copy -> nomegas -> [0]]],
 548 #                          new_values => [\@new_values] );
 549
 550   $copy -> fixed( parameter_type => 'theta',
 551                   new_values => [[(1) x $copy -> nthetas ]] );
 552
 553   $copy -> fixed( parameter_type => 'omega',
 554                   new_values => [[(1) x $copy -> nomegas -> [0] ]] );
 555
 556   $copy -> fixed( parameter_type => 'sigma',
 557                   new_values => [[(1) x $copy -> nsigmas -> [0] ]] );
 558
 559   $copy -> extra_output( ['fort.80'] );
 560
 561   $copy -> remove_records( type => 'table' );
 562   $copy -> remove_records( type => 'abbreviated' );
 563   $copy -> remove_records( type => 'nonparametric' );
 564
 565   $copy -> _write;
 566
 567   push( @new_ofv_models, $copy );
 568 }
 569
 570 my $new_modelfit = tool::modelfit ->
 571     new ( eval( $common_options::parameters ),
 572           prepend_model_file_name => 1,
 573           directory => 'modelfit',
 574           models => \@new_ofv_models );
 575
 576 $new_modelfit -> run;
 577
 578 # Post process fort.80
 579
 580 my $id = 1;
 581 my @all_ind_ofv;
 582
 583 foreach( 1..$number_of_individuals ){
 584   open (IND, '<', "ofv_model_$id.fort.80") || die("Couldn't open ofv_model_$id.mod.fort.80 for reading: $! \n");
 585   $id++;
 586   my @ind_ofv;
 587   while(<IND>) {
 588     chomp;
 589     my ($junk,$ind_ofv) = split(' ',$_);
 590     push @ind_ofv, $ind_ofv;
 591   }
 592   close IND;
 593   @ind_ofv = splice(@ind_ofv,-$number_of_individuals);
 594   push( @all_ind_ofv, \@ind_ofv );
 595 }
 596
 597 dump_csv( "ind_ofv.csv", \@all_ind_ofv );
 598
 599 dump_csv( "ind_etas.csv", \@individual_etas_probabilities );
 600
 601 my @etas_and_sums;
 602
 603 foreach my $parameter( 0..($nomegas - 1) ){
 604   my @result;
 605   my @row_sums = (0) x scalar @{$individual_etas[0]};
 606
 607   for( my $ind_ofv=0; $ind_ofv < $number_of_individuals ; $ind_ofv ++) {
 608     for( my $individual=0;$individual < $number_of_individuals; $individual++) {
 609       my $ofv = $all_ind_ofv[$ind_ofv] -> [$individual];
 610       my $p = $individual_etas_probabilities[$parameter] -> [$individual];
 611       my $LP = exp((-$ofv)/2)*$p;
 612       $result[$ind_ofv] -> [$individual] = $LP;
 613       $row_sums[$individual]+=$LP;
 614     }
 615   }
 616
 617   dump_csv( "L_times_P_ETA$parameter.csv", \@result );
 618
 619   for( my $ind_ofv=0; $ind_ofv < $number_of_individuals ; $ind_ofv ++) {
 620     for( my $individual=0;$individual <= $#{$result[$ind_ofv]}; $individual++ ) {
 621       $result[$ind_ofv] -> [$individual] = ($result[$ind_ofv] -> [$individual] / $row_sums[$individual]) / $number_of_individuals;
 622     }
 623   }
 624
 625   my $big_sum=0;
 626
 627   my @sums;
 628
 629   for( my $row = 0; $row < scalar @result; $row++ ){
 630     my $row_length = $#{$result[$row]};
 631     @{$result[$row]} = @{$result[$row]}[0..$id_to_sim];
 632
 633     my $scale = ($number_of_individuals/$id_to_sim);
 634
 635     my $sum = 0;
 636
 637     for(my $col = 0; $col < scalar @{$result[$row]}; $col++ ){
 638       $result[$row] -> [$col] *= $scale;
 639       $sum += $result[$row] -> [$col];
 640     }
 641
 642     $big_sum += $sum;
 643
 644     #push( @{$result[$row]}, $sum );
 645
 646     push( @sums, $sum );
 647
 648   }
 649
 650   my @order = 0..scalar(@{$individual_etas[$parameter]})-1 ;
 651
 652   @order = sort( { $sums[$a] <=> $sums[$b] } @order );
 653
 654   my @etas_reordered;
 655   my @sums_reordered;
 656   my $i = 0;
 657
 658   foreach my $index( @order ){
 659     $etas_reordered[$i] = $individual_etas[$parameter] -> [$index];
 660     $sums_reordered[$i] = $sums[$index];
 661     $i++;
 662   }
 663
 664   push( @etas_and_sums, \@etas_reordered );
 665   push( @etas_and_sums, \@sums_reordered );
 666   if( $big_sum > 1.01 or $big_sum < 0.99 ){
 667     print "Warning, sum of probabilties is outside of range [0.99..1.01]: $big_sum\n";
 668   }
 669
 670   dump_csv( "P_values_ETA$parameter.csv", \@result );
 671
 672 }
 673
 674 chdir( '..' );
 675
 676 dump_T_csv('results.csv', \@etas_and_sums);
 677
 678 # }}}
 679
 680 # {{{ helper routines
 681
 682 sub dump_csv {
 683   my $file_name = shift;
 684   my $arr = shift;
 685
 686   open( FILE, ">", $file_name );
 687
 688   foreach my $sub_arr( @{$arr} ){
 689     print( FILE join( ',', @{$sub_arr} ), "\n" );
 690   }
 691
 692   close( FILE );
 693
 694 }
 695
 696 sub dump_T_csv {
 697   my $file_name = shift;
 698   my $arr = shift;
 699
 700   open( FILE, ">", $file_name );
 701
 702   for( my $x = 0;$x < @{$arr->[0]}; $x++ ){
 703     for( my $y = 0;$y < @{$arr}; $y++ ){
 704       if( $y > 0 and $y < @{$arr} ){
 705         print( FILE ',' );
 706       }
 707       print( FILE $arr->[$y][$x]  );
 708     }
 709     print( FILE "\n" );
 710   }
 711
 712   close( FILE );
 713 }
 714
 715 # }}}
 716
 717 # }}}