bin/extended_grid

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