lib/output/problem/subproblem_subs.pm

   1 # {{{ include
   2
   3 start include statements
   4     use Config;
   5     use ext::Math::MatrixReal;
   6     use Data::Dumper;
   7     use debug;
   8 end include
   9
  10 # }}} include statements
  11
  12 # {{{ new
  13
  14 start new
  15       {
  16         if ( defined $this -> {'tablename'} ) {
  17           $this -> add_data( 'init_data' => {'filename' => $this -> {'tablename'},
  18                                              'idcolumn' => $this -> {'tableidcolumn'},
  19                                              'debug' => $this -> {'debug'}} );
  20         }
  21         $this -> _read_simulation;
  22         $this -> _read_iteration_path;
  23         $this -> _read_term;
  24         $this -> _read_minimization_message;
  25         $this -> _read_eval;
  26         $this -> _read_significant_digits;
  27         $this -> _read_ofv;
  28         $this -> _read_thomsi;
  29         $this -> _read_sethomsi;
  30         if ($this -> {'covariance_step_run'}) {
  31           $this -> _read_covmatrix;
  32           $this -> _read_eigen;
  33         }
  34         $this -> _read_npomegas;
  35         $this -> _compute_cvsetheta;
  36         $this -> _compute_cvseomega;
  37         $this -> _compute_cvsesigma;
  38
  39         $this -> _compute_comegas;
  40         $this -> _compute_csigmas;
  41         $this -> _set_thetanames;
  42         $this -> _set_omeganames;
  43         $this -> _set_sigmanames;
  44       }
  45 end new
  46
  47 # }}} new
  48
  49 # {{{ register_in_database
  50
  51 start register_in_database
  52     if ( $PsN::config -> {'_'} -> {'use_database'} ) {
  53       my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
  54                                ";databse=".$PsN::config -> {'_'} -> {'project'},
  55                                $PsN::config -> {'_'} -> {'user'},
  56                                 $PsN::config -> {'_'} -> {'password'},
  57                                {'RaiseError' => 1});
  58       my $sth;
  59       my $term_str;
  60       if ( defined $self->{'minimization_message'} ) {
  61         $term_str = join("\n",@{$self->{'minimization_message'}});
  62       }
  63       my @mod_str = ('','');
  64       if ( defined $model_id ) {
  65         @mod_str = ('model_id,',"$model_id,");
  66       }
  67       $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
  68                                ".subproblem ".
  69                              "(problem_id,output_id,".
  70                              $mod_str[0].
  71                              "condition_number,".
  72                              "covariance_step_successful,".
  73                              "covariance_step_warnings,".
  74                              "eval,".
  75                              "final_zero_gradients,".
  76                              "hessian_reset,".
  77                              "ofv,".
  78                              "rounding_errors,".
  79                              "s_matrix_singular,".
  80                              "significant_digits,".
  81                              "simulation,".
  82                              "minimization_successful,".
  83                              "minimization_message,".
  84                              "zero_gradients ) ".
  85                              "VALUES ( '$problem_id' ,".
  86                              "'$output_id' ,".
  87                              $mod_str[1].
  88                              "'$self->{'condition_number'}' ,".
  89                              "'$self->{'covariance_step_successful'}' ,".
  90                              "'$self->{'covariance_step_warnings'}' ,".
  91                              "'$self->{'eval'}' ,".
  92                              "'$self->{'final_zero_gradients'}' ,".
  93                              "'$self->{'hessian_reset'}' ,".
  94                              "'$self->{'ofv'}' ,".
  95                              "'$self->{'rounding_errors'}' ,".
  96                              "'$self->{'s_matrix_singular'}' ,".
  97                              "'$self->{'significant_digits'}' ,".
  98                              "'$self->{'simulation'}' ,".
  99                              "'$self->{'minimization_successful'}' ,".
 100                              "'$term_str' ,".
 101                              "'$self->{'zero_gradients'}' )");
 102       $sth -> execute;
 103       $self -> {'subproblem_id'} = $sth->{'mysql_insertid'};
 104       $self -> {'problem_id'} = $problem_id;
 105       $self -> {'output_id'} = $output_id;
 106       $self -> {'model_id'} = $model_id;
 107       foreach my $param ( 'theta', 'omega', 'sigma' ) {
 108         my $i = 1;
 109         foreach my $par_str ( @{$self -> {$param.'s'}} ) {
 110           $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
 111                                  ".estimate ".
 112                                  "(subproblem_id,problem_id,output_id,".
 113                                  $mod_str[0].
 114                                  "type,value, number, label) ".
 115                                  "VALUES ( '$self->{'subproblem_id'}' ,".
 116                                  "'$self->{'problem_id'}' ,".
 117                                  "'$self->{'output_id'}' ,".
 118                                  $mod_str[1].
 119                                  "'$param','$par_str', '$i', 'test_label')");
 120           $sth -> execute;
 121           push( @{$self -> {'estimate_ids'}}, $sth->{'mysql_insertid'} );
 122           $i++;
 123         }
 124       }
 125       $sth -> finish;
 126       $dbh -> disconnect;
 127     }
 128
 129 # if ( $PsN::config -> {'_'} -> {'use_database'} ) {
 130 #   my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
 131 #                              ";databse=".$PsN::config -> {'_'} -> {'project'},
 132 #                          $PsN::config -> {'_'} -> {'user'},
 133 #                               $PsN::config -> {'_'} -> {'password'},
 134 #                          {'RaiseError' => 1});
 135 #   my $sth = $dbh -> prepare( "UPDATE ".$PsN::config -> {'_'} -> {'project'}.
 136 #                              ".subproblem ".
 137 #                            "SET condition_number='".
 138 #                            $self -> {'condition_number'}."'".
 139 #                            "WHERE model_id='".$self -> {'model_id'}."'" );
 140 #   $sth -> execute or debug -> die( message => $sth->errstr ) ;
 141 #   $sth -> finish;
 142 #   $dbh -> disconnect;
 143 # }
 144 end register_in_database
 145
 146 # }}} register_in_database
 147
 148 # {{{ _compute_comegas_or_csigmas
 149
 150 start _compute_comegas_or_csigmas
 151 {
 152   # This method transforms omegas or sigmas.
 153
 154   # First we check if we have omegas or sigmas.
 155   if ( defined $self -> {'raw_' . $omega_or_sigma . 's'} ) {
 156
 157     my @raw_omegas_or_sigmas = @{$self -> {'raw_' . $omega_or_sigma . 's'}};
 158     # If the omega or sigma matrix has no offdiagonal elements, the
 159     # transformation is quite straightforward.
 160     if( $self -> { $omega_or_sigma . '_block_structure_type'} eq 'DIAGONAL' ){
 161
 162       # All diagonal omegas or sigmas will be in
 163       # @raw_omegas_or_sigmas, so we loop over it and take the square
 164       # root of non-negative, non-zero numbers (those become
 165       # undefined).
 166       my $row_size = 1;
 167       for( my $i = 0; $i <= $#raw_omegas_or_sigmas; $i+=$row_size ){
 168         my $omega_or_sigma_value = sqrt( $raw_omegas_or_sigmas[$i] );
 169         push( @{$self -> {'c' . $omega_or_sigma .'s'}}, $omega_or_sigma_value );
 170         $row_size ++;
 171       }
 172
 173     } elsif( $self -> {$omega_or_sigma . '_block_structure_type'} eq 'BLOCK' ) {
 174
 175       # If the omega or sigma matrix has block form, the transformation
 176       # becomes a bit more complex. The transformation for diagonal
 177       # elements is the same. But for offdiagonals it takes more into
 178       # account.
 179
 180       # $current_row_size is the number of omegas or sigmas in the row
 181       # in the matrix that we are currently transforming.
 182       my $current_row_size = 1;
 183
 184       # $diagonal_index is the index in @raw_omegas_or_sigmas where the current
 185       # on-diagonal omega or sigma we are transforming.
 186       my $diagonal_index = 0;
 187
 188       # $y_idx is the y-coordinate of the element in the matrix we
 189       # are transforming.
 190       my $y_idx = 1;
 191
 192       # First get all block sets sorted.
 193       my @keys = sort {$a <=> $b} keys( %{$self -> { $omega_or_sigma . '_block_sets'}} );
 194 #      print Dumper \@keys;
 195 #      print Dumper $self -> { $omega_or_sigma . '_block_sets'};
 196       foreach my $key ( @keys ) {
 197
 198         # If this is a block set $block_dimension will be > 1
 199         # otherwise it will be 1.
 200         my $block_dimension = $self -> { $omega_or_sigma . '_block_sets'}{$key}{'dimension'};
 201
 202         # The outer loop will be over the rows in the block
 203         # set.
 204         for( my $i = 0; $i < $block_dimension; $i ++ ){
 205
 206           # x_idx is the x-coordinate of the omega or sigma we are
 207           # transforming.
 208
 209           my $x_idx = $current_row_size - $i;
 210 #         print "X: $x_idx, Y: $y_idx\n";
 211
 212           # The inner loop is over the row elements in the block set.
 213           for( my $j = 0; $j <= $i; $j++ ){
 214
 215             # $omega_or_sigma will hold the transformed value.
 216             my $omega_or_sigma_value;
 217
 218             if( $x_idx == $y_idx ) {
 219
 220               # If we are on the diagonal, just take the square root
 221               # of the diagonal element.
 222               $omega_or_sigma_value = $raw_omegas_or_sigmas[ $diagonal_index ];
 223               $omega_or_sigma_value = sqrt( $omega_or_sigma_value );
 224
 225             } else {
 226               # If we are of the diagonal, we need to find two
 227               # on-diagonal omegas, one on the same column and one on
 228               # the same row.
 229
 230               # The equation y*(y+1)/2 is a geometric sum that gives
 231               # the index of the diagonal element on row y in the
 232               # compacted matrixarray @raw_omegas_or_sigmas. The equation
 233               # x*(x+1)/2 gives the diagonal element on column x in
 234               # the same array.
 235
 236               # We find the omega by finding the index of the diagonal
 237               # element on row y and subtracting the the difference
 238               # between the number of elements on the row ($i) and the
 239               # number of elements we have transformed ($j).
 240
 241               $omega_or_sigma_value = $raw_omegas_or_sigmas[(($y_idx*($y_idx+1))/2)-1 - ($i - $j)];
 242               my $denominator = $raw_omegas_or_sigmas[ (($x_idx*($x_idx+1))/2)-1 ]
 243                                  *
 244                                 $raw_omegas_or_sigmas[ (($y_idx*($y_idx+1))/2)-1 ];
 245               if( $denominator <= 0 ){ # To avoiding division by zero
 246                 $omega_or_sigma_value = undef;
 247               } elsif( $omega_or_sigma_value >= sqrt($denominator) ) {
 248                 # This rounding handles cases when the offdiagonals
 249                 # are greater or equal to one.
 250                 $omega_or_sigma_value = $omega_or_sigma_value/( int( 10000 * sqrt($denominator) )/10000 )
 251               } else {
 252                 $omega_or_sigma_value = $omega_or_sigma_value/sqrt($denominator);
 253               }
 254             }
 255             push( @{$self -> {'c' . $omega_or_sigma . 's'}}, $omega_or_sigma_value );
 256
 257             # Move the x-coordinate forwards.
 258             $x_idx ++;
 259           }
 260
 261           # Move the y-coordinate forwards.
 262           $y_idx ++;
 263           # Calculate the new row size.
 264           $current_row_size ++;
 265           # Skip one row and we get the index of the next diagonal.
 266
 267           $diagonal_index += $current_row_size;
 268         }
 269       }
 270     } else {
 271       'debug' -> warn( level => 1,
 272                        message => $omega_or_sigma . ' matrix has unknown form.' );
 273     }
 274   }
 275 }
 276 end _compute_comegas_or_csigmas
 277
 278 # }}} _compute_comegas_or_sigmas
 279
 280 # {{{ _compute_comegas
 281 start _compute_comegas
 282 {
 283   $self -> _compute_comegas_or_csigmas( omega_or_sigma => 'omega' );
 284 }
 285 end _compute_comegas
 286
 287 # This is the old algoritm for transforming omegas. It can't handle
 288 # of-diagonals omegas that are zero, if that should ever happen.
 289
 290 #     my $i;
 291 #     my @raw_omegas;
 292 #     my @computed_omegas;
 293 #     my @diags;
 294 #     my @indices;
 295
 296 #     @raw_omegas = @{$self -> {'raw_omegas'}};
 297 #     @computed_omegas = ();
 298
 299 #     if( $self -> {'omega_structure_type'} eq 'DIAGONAL' ){
 300 #       @diags = @raw_omegas;
 301 #     } else {
 302 #       ## Collect the diagonals
 303 #       foreach $i (0..$#raw_omegas) {
 304 #       if( $self -> _isdiagonal('index' => $i+1) ){
 305 #         push( @diags, $raw_omegas[$i] );
 306 #       }
 307 #       }
 308 #     }
 309
 310 #     foreach $i (0..$#raw_omegas){
 311 #       if($self -> _isdiagonal('index' => $i+1)) {
 312 #       if( $raw_omegas[$i] <= 0 ){
 313 #         push( @computed_omegas, $raw_omegas[$i] );
 314 #         next;
 315 #       }
 316 #       push @computed_omegas, sqrt($raw_omegas[$i]);
 317 #       next;
 318 #       } else {
 319 #       @indices = $self -> _rowcolind( index => $i+1 );
 320 #       if( ($raw_omegas[$i] == 0) or ($diags[$indices[0]-1]*$diags[$indices[1]-1]) <= 0 ) {
 321 #         next;
 322 #       } else {
 323 #         push @computed_omegas,
 324 #         $raw_omegas[$i]/sqrt($diags[$indices[0]-1]*$diags[$indices[1]-1]);
 325 #       }
 326 #       }
 327 #     }
 328 #   $self -> {'comegas'}= \@computed_omegas;}
 329
 330 # }}}
 331
 332 # {{{ _compute_csigmas
 333 start _compute_csigmas
 334 {
 335   $self -> _compute_comegas_or_csigmas( omega_or_sigma => 'sigma' );
 336 }
 337 end _compute_csigmas
 338 # This is the old algoritm for transforming sigmass. It can't handle
 339 # of-diagonals sigmass that are zero, if that should ever happen.
 340
 341 # {
 342 #     if (defined $self -> {'raw_sigmas'} ) {
 343 #       my $i;
 344 #       my ( @om, @com, @diags, @indices );
 345 #       @om = @{$self -> {'raw_sigmas'}};
 346 #       @com = ();
 347 #       ## Collect the diagonals
 348 #       foreach $i (0..$#om) {
 349 #           push @diags,$om[$i] if $self -> _isdiagonal('index' => $i+1);
 350 #       }
 351 #       foreach $i (0..$#om){
 352 #           if($self -> _isdiagonal('index' => $i+1)) {
 353 #             if( $om[$i] <= 0 ){
 354 #               next;
 355 #             }
 356 #             push @com, sqrt($om[$i]);
 357 #             next;
 358 #           } else {
 359 #             @indices = $self -> _rowcolind( index => $i+1 );
 360 #             if( ($om[$i] == 0) or ($diags[$indices[0]-1]*$diags[$indices[1]-1] <= 0) ) {
 361 #               next;
 362 #             } else {
 363 #               push @com,
 364 #               $om[$i]/sqrt($diags[$indices[0]-1]*$diags[$indices[1]-1]);
 365 #             }
 366 #           }
 367 #       }
 368 #       $self -> {'csigmas'}= \@com;
 369 #     }
 370 # }
 371
 372 # }}} _compute_csigmas
 373
 374 # {{{ _compute_cvseomega
 375
 376 start _compute_cvseomega
 377 {
 378     my @raw_omegas   = @{$self -> {'raw_omegas'}};
 379     my @raw_seomegas = @{$self -> {'raw_seomegas'}};
 380     my @cvseomega;
 381
 382     if ( scalar @raw_seomegas > 0) {
 383         my ($j,$i);
 384         my ( @init_om, @init_seom );
 385
 386         foreach $i (0.. $#raw_seomegas){
 387             if( $self -> _isdiagonal('index' => $i+1)) {
 388                 push @init_seom, $raw_seomegas[$i];
 389                 next;
 390             } else {
 391                 if(($raw_seomegas[$i] eq 'NA') or ($raw_seomegas[$i] == 0)) {
 392                     next;
 393                 } else {
 394                     push @init_seom, $raw_seomegas[$i];
 395                 }
 396             }
 397         }
 398
 399         foreach $i (0..$#raw_omegas){
 400             if( $self -> _isdiagonal('index' => $i+1)) {
 401                 push @init_om, $raw_omegas[$i];
 402                 next;
 403             } else {
 404                 if($raw_omegas[$i] == 0) {
 405                     next;
 406                 } else {
 407                     push @init_om, $raw_omegas[$i];
 408                 }
 409             }
 410         }
 411
 412         foreach my $i (0..$#init_om) {
 413           if( $init_seom[$i] ne 'NA' and $init_om[$i] ne 'NA' ){
 414             push @cvseomega,$init_seom[$i]/abs($init_om[$i]);
 415           } else {
 416             push @cvseomega,undef;
 417           }
 418 #           push @cvseomega,$init_seom[$i]/$init_om[$i]
 419 #               if ($init_seom[$i] && $init_om[$i]);
 420 #           push @cvseomega,"NA" unless ($init_seom[$i] && $init_om[$i]);
 421         }
 422
 423     }
 424
 425     $self -> {'cvseomegas'} = [@cvseomega];
 426 }
 427 end _compute_cvseomega
 428
 429 # }}} _compute_cvseomega
 430
 431 # {{{ _compute_cvsesigma
 432 start _compute_cvsesigma
 433 {
 434     my @raw_sigmas   = @{$self -> {'raw_sigmas'}};
 435     my @raw_sesigmas = @{$self -> {'raw_sesigmas'}};
 436     my @cvsesigma;
 437
 438     if ( scalar @raw_sesigmas > 0) {
 439         my ($j,$i);
 440         my ( @init_si, @init_sesi );
 441
 442         foreach $i (0.. $#raw_sesigmas){
 443             if($self -> _isdiagonal('index' => $i+1)) {
 444                 push @init_sesi, $raw_sesigmas[$i];
 445                 next;
 446             } else {
 447                 if(($raw_sesigmas[$i] eq 'NA') or ($raw_sesigmas[$i] == 0)) {
 448                     next;
 449                 } else {
 450                     push @init_sesi, $raw_sesigmas[$i];
 451                 }
 452             }
 453         }
 454
 455         foreach $i (0.. $#raw_sigmas){
 456             if($self -> _isdiagonal('index' => $i+1)) {
 457                 push @init_si, $raw_sigmas[$i];
 458                 next;
 459             } else {
 460                 if($raw_sigmas[$i] == 0) {
 461                     next;
 462                 } else {
 463                     push @init_si, $raw_sigmas[$i];
 464                 }
 465             }
 466         }
 467
 468         foreach my $i (0..$#init_si) {
 469           if( $init_sesi[$i] ne 'NA' and $init_si[$i] ne 'NA' ){
 470             push @cvsesigma,$init_sesi[$i]/abs($init_si[$i]);
 471           } else {
 472             push @cvsesigma,undef;
 473           }
 474         }
 475
 476     }
 477
 478     $self -> {'cvsesigmas'}= [@cvsesigma];
 479 }
 480 end _compute_cvsesigma
 481 # }}} _compute_cvsesigma
 482
 483 # {{{ _compute_cvsetheta
 484 start _compute_cvsetheta
 485 {
 486     my @thetas   = @{$self -> {'thetas'}};
 487     my @sethetas = @{$self -> {'sethetas'}};
 488     my @cvsethetas;
 489
 490     if ( scalar @sethetas > 0 ) {
 491         foreach my $i (0..$#thetas) {
 492             push @cvsethetas,$sethetas[$i]/abs($thetas[$i])
 493                 if ($sethetas[$i] && $thetas[$i]);
 494             push @cvsethetas,"NA" unless ($sethetas[$i] && $thetas[$i]);
 495         }
 496     }
 497
 498     $self -> {'cvsethetas'}= [@cvsethetas];
 499 }
 500 end _compute_cvsetheta
 501 # }}} _compute_cvsetheta
 502
 503 # {{{ _isdiagonal
 504 start _isdiagonal
 505 {
 506     my $previ = 1;
 507     my $j;
 508     return(1) if $index == 1;
 509     foreach my $j (2..100) {
 510         return(1) if $index == $previ+$j;
 511         $previ = $previ+$j;
 512         last if $index < $previ;
 513     }
 514 }
 515 end _isdiagonal
 516 # }}} _isdiagonal
 517
 518 # {{{ _read_matrixoestimates
 519 start _read_matrixoestimates
 520       {
 521         # Reads one matrix structure and returns the file handle at
 522         # the beginning of the next structure
 523         while ( $_ = @{$self -> {'lstfile'}}[ $pos++ ] ) {
 524           last if (/^\s*\*/);
 525           # Rewind one step if we find something that marks the end of
 526           # our structure
 527           $pos-- and last if ( /^ PROBLEM.*SUBPROBLEM/ or /^ PROBLEM NO\.:\s+\d/ );
 528           $pos-- and last if (/^[a-df-zA-DF-Z]/);
 529
 530           next if ( /^\s+TH/ or /^\s+OM/ or /^\s+SG/ );   # Header row
 531           next if ( /^1/ );                       # Those annoying 1's
 532
 533           chomp;                                # Get rid of line-feed
 534           my @row = split;
 535           shift( @row ) if ( $row[0] eq '+' );     # Get rid of +-sign
 536
 537           next if ( $#row < 0 );                           # Blank row
 538
 539           push( @subprob_matrix, @row );
 540         }
 541         $success = 1 if( scalar @subprob_matrix > 0 );
 542       }
 543 end _read_matrixoestimates
 544 # }}} _read_matrixoestimates
 545
 546 # {{{ _read_covmatrix
 547 start _read_covmatrix
 548       {
 549         my $matrix_nr = 0;
 550         my ( $t_success, $c_success, $corr_success, $i_success ) = (0,0,0,0);
 551         my $start_pos = $self -> {'lstfile_pos'}-1;
 552
 553         # {{{ sub clear dots
 554
 555         sub clear_dots {
 556           my $m_ref = shift;
 557           my @matrix = @{$m_ref};
 558           # get rid of '........'
 559           my @clear;
 560           foreach ( @matrix ) {
 561             push( @clear, $_ ) unless ( $_ eq '.........' );
 562           }
 563           #       print Dumper \@clear;
 564           return \@clear;
 565         }
 566
 567         # }}}
 568
 569         # {{{ sub make square
 570
 571         sub make_square {
 572           my $m_ref = shift;
 573           my @matrix = @{$m_ref};
 574           # Make the matrix square:
 575           my $elements = scalar @matrix; # = M*(M+1)/2
 576           my $M = -0.5 + sqrt( 0.25 + 2 * $elements );
 577           my @square;
 578           for ( my $m = 1; $m <= $M; $m++ ) {
 579             for ( my $n = 1; $n <= $m; $n++ ) {
 580               push( @{$square[$m-1]}, $matrix[($m-1)*$m/2 + $n - 1] );
 581               unless ( $m == $n ) {
 582                 push( @{$square[$n-1]}, $matrix[($m-1)*$m/2 + $n - 1] );
 583               }
 584             }
 585           }
 586           return \@square;
 587         }
 588
 589         # }}}
 590
 591         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 592           if (/T MATRIX/) {
 593             while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 594               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 595                 ( $start_pos, $self -> {'raw_tmatrix'}, $t_success )  = $self ->
 596                     _read_matrixoestimates( pos => $start_pos ) and last;
 597               }
 598             }
 599             last;                # No covariance matrix will be found!
 600           }
 601           if (/    COVARIANCE MATRIX OF ESTIMATE/) {
 602             while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 603               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 604                 ( $start_pos, $self -> {'raw_covmatrix'}, $c_success ) = $self ->
 605                     _read_matrixoestimates( pos => $start_pos ) and last;
 606               }
 607             }
 608           }
 609           if (/    CORRELATION MATRIX OF ESTIMATE/) {
 610             while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 611               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 612                 ( $start_pos, $self -> {'raw_cormatrix'}, $corr_success ) = $self ->
 613                     _read_matrixoestimates( pos => $start_pos ) and last;
 614               }
 615             }
 616           }
 617           if (/    INVERSE COVARIANCE MATRIX OF ESTIMATE/) {
 618             while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 619               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 620                 ( $start_pos, $self -> {'raw_invcovmatrix'}, $i_success ) = $self ->
 621                     _read_matrixoestimates( pos => $start_pos ) and last;
 622               }
 623             }
 624             last;                                       # Last matrix?
 625           }
 626         }
 627         foreach my $element ( @{$self -> {'raw_tmatrix'}} ) {
 628           push( @{$self -> {'t_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 629         }
 630         foreach my $element ( @{$self -> {'raw_covmatrix'}} ) {
 631           push( @{$self -> {'covariance_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 632         }
 633         foreach my $element ( @{$self -> {'raw_cormatrix'}} ) {
 634           push( @{$self -> {'correlation_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 635         }
 636
 637         if( defined $self -> {'raw_invcovmatrix'} ) {
 638             $self -> {'inverse_covariance_matrix'} = Math::MatrixReal ->
 639                 new_from_cols( make_square( clear_dots( $self -> {'raw_invcovmatrix'} ) ) );
 640         }
 641
 642 #       foreach my $element ( @{$self -> {'raw_invcovmatrix'}} ) {
 643 #         push( @{$self -> {'inverse_covariance_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 644 #       }
 645
 646         #If something has gone right!
 647         $self-> {'lstfile_pos'} = $start_pos if ( $t_success + $c_success + $corr_success + $i_success  );
 648       }
 649 end _read_covmatrix
 650 # }}} _read_covmatrix
 651
 652 # {{{ _read_eigen
 653
 654 start _read_eigen
 655       {
 656         my @eigens;
 657         my $start_pos = $self -> {'lstfile_pos'};
 658         my $eig_area  = 0;
 659         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 660           chomp;
 661           if ( /EIGENVALUES OF COR MATRIX OF ESTIMATE/ ) {
 662             $eig_area = 1;
 663             $start_pos = $start_pos + 4 ; # Jump forward to the index numbers
 664             debug -> warn( level   => 2,
 665                            message => "Found the eigenvalue area" );
 666           INNER: while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) { # Get rid of indexes
 667               last INNER if ( not /^\s+\d/ );
 668 #             $start_pos++ and last INNER if ( not /^\s+\d/ );
 669             }
 670           }
 671           if ( $eig_area ) {
 672             $start_pos-- and last if (/^[a-df-zA-DF-Z]/); #Rewind one step
 673             last if ( /^\s*\*/ or /^1/ );
 674             push( @eigens, split );
 675           }
 676           $start_pos-- and last if ( /^ PROBLEM.*SUBPROBLEM/ or /^ PROBLEM NO\.:\s+\d/ );
 677           $start_pos-- and last if (/^[a-df-zA-DF-Z]/); #Rewind one step
 678         }
 679         if ( scalar @eigens > 0 ) {
 680           my @list = sort { $a <=> $b } @eigens;
 681           $self -> {'condition_number'} = abs($list[$#list]/$list[0]) if ( $list[0] != 0 );
 682         }
 683         @{$self -> {'eigens'}} = @eigens;
 684       }
 685 end _read_eigen
 686
 687 # }}} _read_eigen
 688
 689 # {{{ _read_eval
 690 start _read_eval
 691 {
 692     my $start_pos = $self -> {'lstfile_pos'};
 693     my $success  = 0;
 694
 695     while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 696         if ( /MINIMUM VALUE OF OBJECTIVE FUNCTION/ or
 697              / NO. OF SIG. DIGITS IN FINAL EST./){
 698             debug -> warn( level   => 2,
 699                            message => "Hmmm, we have gone too far" );
 700             last;
 701         }
 702         if ( / NO. OF FUNCTION EVALUATIONS USED:\s*(\d*)/ ){
 703             $self -> {'feval'} = $1;
 704             $success = 1;
 705             last;
 706         }
 707     }
 708     unless ( $success ) {
 709         debug -> warn( level   => 2,
 710                        message => "rewinding to first position..." );
 711     } else {
 712         debug -> warn( level   => 2,
 713                        message => "Found number of Function evaluations" );
 714         $self -> {'lstfile_pos'} = $start_pos;
 715     }
 716 }
 717 end _read_eval
 718 # }}} _read_eval
 719
 720 # {{{ _read_iteration_path
 721 start _read_iteration_path
 722       {
 723         my $start_pos = $self -> {'lstfile_pos'};
 724         my $success = 0;
 725         my (@func_eval, @parameter_vectors,
 726             @gradient_vectors) = ((), (), (), (), ());
 727
 728         my $cumulative_evals = 0;
 729         my $zero_gradients = 0;
 730         my $hessian_reset = 0;
 731         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) { #Large reading loop
 732           if (/^0ITERATION NO/) {
 733             debug -> die( message => "Error in reading iteration path!\n$!" )
 734               unless (/0ITERATION NO\.:\s+(\d+)\s+OBJECTIVE VALUE:\s+(\S+)\s+NO\. OF FUNC\. EVALS\.:\s*(.+)/);
 735             $success = 1;
 736             push(@{$self -> {'iternum'}}, $1);
 737
 738             my $ofvpath = $2;
 739             unless( $ofvpath eq '**' ){ # If funcion evals are more than 10000, NONMEM will print out two stars.
 740                 push(@{$self -> {'ofvpath'}}, $ofvpath );
 741             } # If, in fact, we find stars, the number of evaluations are calculated below
 742
 743             my (@parameter_vector, @gradient_vector) = ((), ());
 744
 745             while ( $_ = @{$self -> {'lstfile'}}[ $start_pos ] ) {
 746               if (/^ CUMULATIVE NO\. OF FUNC\. EVALS\.:\s*(\d+)/) {
 747                 my $eval_path = $1;
 748
 749                 if( $ofvpath eq '**' ){
 750                     my $ofvpath = $eval_path - $cumulative_evals;
 751                     $cumulative_evals = $eval_path;
 752                     debug -> warn( level   => 2,
 753                                    message => "Calculated eval_path = $ofvpath" );
 754                     push(@{$self -> {'ofvpath'}}, $ofvpath );
 755                 }
 756
 757                 push(@{$self -> {'funcevalpath'}}, $eval_path);
 758
 759                 if (/RESET HESSIAN, TYPE (\w+)/) {
 760                   $hessian_reset++;
 761                 }
 762
 763                 $start_pos++;
 764               } elsif ( s/^ PARAMETER:\s*// ) {
 765                 do {
 766                   push(@parameter_vector, split);
 767                   $_ = @{$self -> {'lstfile'}}[ ++$start_pos ];
 768                 } while ( defined($_) and not /^ GRADIENT:\s*/ );
 769               } elsif (s/^ GRADIENT:\s*//) {
 770                 do {
 771                   push(@gradient_vector, split);
 772                   $_ = @{$self -> {'lstfile'}}[ ++$start_pos ]
 773                 } while ( defined($_)  and not /[A-D][F-X]/ );
 774               } else {
 775                 last;
 776               }
 777             }
 778             foreach my $grad ( @gradient_vector ) {
 779               $zero_gradients++ if $grad == 0;
 780             }
 781             $self -> {'initgrad'} = \@gradient_vector unless($self -> {'initgrad'});
 782             $self -> {'final_gradients'}  = \@gradient_vector;
 783             $self -> {'finalparam'} = \@parameter_vector;
 784             push(@parameter_vectors, \@parameter_vector);
 785             push(@gradient_vectors, \@gradient_vector);
 786
 787             last unless(/^0ITERATION NO/);
 788           }                     #End of if iteration no
 789         }                       #End of large reading loop
 790
 791         my ($kill_found, $file_end, $kill_message, $search_pos) = (0, 0, "", $start_pos);
 792         while ( $_ = @{$self -> {'lstfile'}}[ $search_pos++ ] ) { #Have a look, a few lines down...
 793           if( /kill/i ) {
 794             $kill_found = 1;
 795             $kill_message = $_;
 796             last;
 797           }
 798 #         $kill_found = 1 and $kill_message = $_ and last if(/kill/i);
 799           if( $search_pos + 1 == scalar @{$self -> {'lstfile'}} ) {
 800             $file_end = 1;
 801             $search_pos = $start_pos + 4;
 802           }
 803           last if( $search_pos > $start_pos + 3 )
 804 #         last if (--$kill_found < -2);
 805         }
 806         if (($kill_found == 1) or $file_end) { #Crash before last iteration
 807           $self -> {'minimization_message'} = ["PsN message:","The output file seems to have an abrupt ending," .
 808                                               " before the last","iteration has finished."];
 809           if ($kill_found == 1) {
 810             push(@{$self -> {'minimization_message'}}, "String found in output file: $kill_message" );
 811           } else {
 812             push(@{$self -> {'minimization_message'}}, " This is probably due to a crash");
 813           }
 814         }
 815
 816         unless ( $success ){
 817           debug -> warn( level   => 2,
 818                          message => "rewinding to first position..." );
 819         } else {
 820           $self -> {'lstfile_pos'} = $start_pos;
 821           $self -> {'parameter_path'}    = \@parameter_vectors;
 822           $self -> {'gradient_path'} = \@gradient_vectors;
 823           $self -> {'zero_gradients'} = $zero_gradients;
 824           my $final_zero_gradients = 0;
 825           foreach my $grad ( @{$self -> {'final_gradients'}} ) {
 826             $final_zero_gradients++ if $grad == 0;
 827           }
 828           $self -> {'final_zero_gradients'} = $final_zero_gradients;
 829           $self -> {'hessian_reset'} = $hessian_reset;
 830         }
 831       }
 832 end _read_iteration_path
 833 # }}} _read_iteration_path
 834
 835 # {{{ _read_npomegas
 836
 837 start _read_npomegas
 838 {
 839     my $start_pos = $self -> {'lstfile_pos'};
 840     my $success   = 0;
 841     my $npofvarea   = 0;
 842     my $nparea  = 0;
 843     my $npetabararea    = 0;
 844     my $npomegarea      = 0;
 845     my ( @npetabar, @npomega, @T, $i );
 846
 847     while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 848         $nparea      = 1   if /NONPARAMETRIC ESTIMATE/;
 849         last if ( /THERE ARE ERROR MESSAGES IN FILE PRDERR/ and $nparea );
 850         last if ( /^1/ and $nparea );
 851         last if ( /^1NONLINEAR/ and $nparea );
 852         last if ( /^[A-W]/ and $nparea );
 853
 854         if (/MINIMUM VALUE OF OBJECTIVE FUNCTION/ and $nparea ){ #Only nonmem6 version
 855             $npofvarea      = 1;
 856         }
 857         if ( /EXPECTED VALUE OF ETA/ and $nparea ) {
 858             $npetabararea = 1;
 859             $npofvarea = 0;
 860             $success   = 1;
 861         }
 862         if ( /COVARIANCE MATRIX OF ETA/ and $nparea ) {
 863             $npomegarea = 1;
 864             $npetabararea = 0;
 865         }
 866         if($npofvarea and /^\s+(-?\d*\.\d*)/) { #Assignment of attribute at the spot
 867             $self -> {'npofv'} = $1;
 868             $npofvarea = 0;
 869         }
 870         if($npetabararea and /^\s*-?\d*\.\d*/) {
 871             @T = split(' ',$_);
 872             for $i (0..(@T-1)) {$T[$i] = eval($T[$i]);}
 873             push(@npetabar,@T);
 874         }
 875         if($npomegarea and /^(\+|\s{2,})/) {
 876             next if /ET/;
 877             @T = split(' ',$_);
 878             shift @T if $T[0] eq '+';
 879             for  $i (0..(@T-1)) {$T[$i] = eval($T[$i]);}
 880             push(@npomega,@T);
 881         }
 882     }
 883     $self -> {'npetabar'} = [@npetabar];
 884     $self -> {'npomegas'} = [@npomega];
 885     unless ( $success ) {
 886         debug -> warn( level   => 2,
 887                        message => "rewinding to first position..." );
 888     } else {
 889         $self -> {'lstfile_pos'} = $start_pos;
 890     }
 891 #    print Dumper $self -> {'npomegas'};
 892 #    print Dumper $self -> {'npetabar'};
 893 }
 894 end _read_npomegas
 895
 896 # }}} _read_npomegas
 897
 898 # {{{ _read_ofv
 899
 900 start _read_ofv
 901 {
 902     my $start_pos = $self -> {'lstfile_pos'};
 903     my $success  = 0;
 904
 905     while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 906         if ( /^\s\*{50}\s+/ ) {
 907             (undef, my $ofvt, undef) = split(' ',$_,3);
 908             if ( $ofvt =~ /\*\*\*\*\*\*/ ) {
 909                 $self -> {'ofv'} = undef;
 910             } else {
 911                 $self -> {'ofv'} = $ofvt;
 912             }
 913             debug -> warn( level   => 2,
 914                            message => "OFV = $ofvt" );
 915             $success = 1;
 916             last;
 917         }
 918     }
 919     unless ( $success ) {
 920         debug -> warn( level   => 2,
 921                        message => "rewinding to first position..." );
 922     } else {
 923         debug -> warn( level   => 2,
 924                        message => "Found ofv" );
 925         $self -> {'lstfile_pos'} = $start_pos;
 926     }
 927 }
 928 end _read_ofv
 929
 930 # }}} _read_ofv
 931
 932 # {{{ _rowcolind
 933 start _rowcolind
 934       {
 935         my $i = 1;
 936         my $found = 0;
 937         while ( not $found ) {
 938           my $test = $index - ($i-1)*($i)/2;
 939           if ( $test <= $i ) {
 940             $row = $i;
 941             $col = $test;
 942             $found = 1;
 943           }
 944           $i++;
 945         }
 946 #       my ($i,$j);
 947 #       my @startind = (2,5,9,14,20,27,35,44);
 948 #       my $prevind;
 949 #       my $notfirst = 0;
 950
 951 #       OUTER: foreach $col (1..8) {
 952 #         $prevind = $startind[$col-1];
 953 #         foreach $row (($col+1)..8) {
 954
 955 #           ## If this is the first element of the column
 956 #           if($index == $prevind) {
 957 #             last OUTER;
 958 #           }
 959
 960 #           ## If it is a later element in the column
 961 #           if($index == $prevind+$row) {
 962 #       #      print "$index $col ",$row+1," ",$startind[$col-1],"\n" ;
 963 #             $row++;
 964 #             last OUTER;
 965 #           } else {
 966 #             $prevind = $prevind+$row;
 967 #             next;
 968 #           }
 969 #         }
 970 #       }
 971 #       print "$index ",$row," $col\n";
 972
 973       }
 974 end _rowcolind
 975 # }}} _rowcolind
 976
 977 # {{{ _read_significant_digits
 978 start _read_significant_digits
 979 {
 980     my $start_pos = $self -> {'lstfile_pos'};
 981     my $success  = 0;
 982
 983     while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 984         if ( /MINIMUM VALUE OF OBJECTIVE FUNCTION/ ){
 985             last;
 986         }
 987         if ( / NO. OF SIG. DIGITS IN FINAL EST.:\s*(-?\d*\.*\d*)/ ){
 988             $self -> {'significant_digits'} = $1;
 989             $success = 1;
 990             last;
 991         }
 992     }
 993     unless ( $success ) {
 994         debug -> warn( level   => 2,
 995                        message => "No Number of Significant digits found" );
 996     } else {
 997         debug -> warn( level   => 2,
 998                        message => "Found significant_digits" );
 999         $self -> {'lstfile_pos'} = $start_pos;
1000     }
1001 }
1002 end _read_significant_digits
1003 # }}} _read_significant_digits
1004
1005 # {{{ _read_sethomsi
1006 start _read_sethomsi
1007       {
1008         my $start_pos = $self -> {'lstfile_pos'};
1009         my $success  = 0;
1010
1011         my $searea   = 0;
1012         my $thetarea = 0;
1013         my $omegarea = 0;
1014         my $sigmarea = 0;
1015         my ( @setheta, @seomega, @sesigma, @T, $i, $tmp );
1016         my ( @raw_setheta, @raw_seomega, @raw_sesigma );
1017         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
1018           chomp;
1019           $searea   = 1 if /STANDARD ERROR OF ESTIMATE/;
1020           if ( $searea ) {
1021             if ( /THETA - VECTOR OF FIXED EFFECTS PARAMETERS/ ) {
1022               debug -> warn( level   => 2,
1023                              message => "Found standard erros of thetas" );
1024               $success  = 1;
1025               $thetarea = 1;
1026             }
1027             if ( /OMEGA - COV MATRIX FOR RANDOM EFFECTS - ETAS/ ) {
1028               debug -> warn( level   => 2,
1029                              message => "Found standard erros of omega" );
1030               $omegarea = 1;
1031               $thetarea = 0;
1032             }
1033             if ( /SIGMA - COV MATRIX FOR RANDOM EFFECTS - EPSILONS/ ) {
1034               debug -> warn( level   => 2,
1035                              message => "Found standard erros of sigmas" );
1036               $sigmarea = 1;
1037               $omegarea = 0;
1038               $thetarea = 0;
1039             }
1040             last if ( / COVARIANCE MATRIX OF ESTIMATE/ );
1041             last if ( /THERE ARE ERROR MESSAGES IN FILE PRDERR/ );
1042             # The row below checking for a one '1', crashes behaviour with MANY etas (>15)
1043 #           last if ( /^1/ );
1044             last if ( /^[A-W]/ and $searea );
1045             if ( $thetarea and /^\s*-?\d*\.\d*/ ) {
1046               @T = split(' ',$_);
1047               for $i (0..(@T-1)) {
1048                 $T[$i] = eval($T[$i]);
1049               }
1050               push(@setheta,@T);
1051             }
1052             if($omegarea and /^(\+|\s{2,})/) {
1053               next if /ET/;
1054               @T = split(' ',$_);
1055               shift @T if $T[0] eq '+';
1056               if ( not defined $self -> {'omega_block_structure'} ) {
1057                 for $i (0..(@T-1)) {
1058                   if(($T[$i] ne '.........') and (eval($T[$i]) != 0) ) {
1059                     push(@seomega,eval($T[$i]));
1060                   }
1061                 }
1062               }
1063               for  $i (0..(@T-1)) {
1064                 if($T[$i] ne '.........') {
1065                   $tmp = eval($T[$i]);
1066                 } else {
1067                   $tmp = 'NA';
1068                 }
1069                 $T[$i] = $tmp ;
1070               }
1071               push(@raw_seomega,@T);
1072             }
1073             if($sigmarea and /^(\+|\s{2,})/) {
1074               next if /EP/;
1075               @T = split(' ',$_);
1076               shift @T if $T[0] eq '+';
1077               if ( not defined $self -> {'sigma_block_structure'} ) {
1078                 for  $i (0..(@T-1)) {
1079                   if(($T[$i] ne '.........') and (eval($T[$i]) != 0) ) {
1080                     push(@sesigma,eval($T[$i]));
1081                   }
1082                 }
1083               }
1084               for $i (0..(@T-1)) {
1085                 if ($T[$i] ne '.........') {
1086                   $tmp = eval($T[$i]);
1087                 } else {
1088                   $tmp = 'NA';
1089                 }
1090                 $T[$i] = $tmp ;
1091               }
1092               push(@raw_sesigma,@T);
1093             }
1094           }
1095         }
1096
1097
1098         if ( defined $self -> {'omega_block_structure'} ) {
1099           my @omblock  = @{$self -> {'omega_block_structure'}};
1100           debug -> warn( level   => 2,
1101                          message => "OMEGA BLOCK DEFINED" ) if ( scalar @omblock > 0 );
1102           @seomega = ();
1103           my $i = 0;
1104           foreach my $row ( @omblock ) {
1105             foreach my $element ( @{$row} ) {
1106               push ( @seomega, $raw_seomega[$i] ) if ( $element );
1107               $i++;
1108             }
1109           }
1110         }
1111
1112         if ( defined $self -> {'sigma_block_structure'} ) {
1113           my @siblock  = @{$self -> {'sigma_block_structure'}};
1114           @sesigma = ();
1115           my $i = 0;
1116           debug -> warn( level   => 2,
1117                          message => "SIGMA BLOCK DEFINED" ) if ( scalar @siblock > 0 );
1118           foreach my $row ( @siblock ) {
1119             foreach my $element ( @{$row} ) {
1120               push ( @sesigma, $raw_sesigma[$i] ) if ( $element );
1121               $i++;
1122             }
1123           }
1124         }
1125
1126         $self -> {'sethetas'}            = [@setheta];
1127         $self -> {'raw_seomegas'}        = [@raw_seomega];
1128         $self -> {'raw_sesigmas'}        = [@raw_sesigma];
1129         $self -> {'seomegas'}            = [@seomega];
1130         $self -> {'sesigmas'}            = [@sesigma];
1131
1132         if ( scalar @setheta <= 0 ) {
1133           $self -> {'covariance_step_successful'} = 0;
1134         } else {
1135           $self -> {'covariance_step_successful'} = 1;
1136         }
1137
1138         unless ( $success ) {
1139           debug -> warn( level   => 2,
1140                          message => "No standard errors for thetas, sigmas or omegas." );
1141         } else {
1142           $self -> {'lstfile_pos'} = $start_pos;
1143         }
1144       }
1145 end _read_sethomsi
1146 # }}} _read_sethomsi
1147
1148 # {{{ _read_simulation
1149 start _read_simulation
1150     {
1151       my $start_pos = $self -> {'lstfile_pos'};
1152       while ( $_ = @{$self -> {'lstfile'}}[ $start_pos ++ ] ) {
1153         if ( /MINIMUM VALUE OF OBJECTIVE FUNCTION/ ) {
1154           last;
1155         }
1156         if ( /^\s*MONITORING OF SEARCH:/) {
1157           last;
1158         }
1159         if ( /\s*SIMULATION STEP PERFORMED/ ) {
1160           $self ->{'simulationstep'} = 1;
1161           last;
1162         }
1163       }
1164
1165       if ( $self -> {'simulationstep'} ) {
1166         $self -> {'lstfile_pos'} = $start_pos;
1167       }
1168     }
1169 end _read_simulation
1170 # }}} _read_simulation
1171
1172 # {{{ _read_term
1173
1174 start _read_term
1175       {
1176         my $start_pos = $self -> {'lstfile_pos'};
1177         my $success_pos;
1178         my $success  = 0;
1179
1180         $self -> {'minimization_successful'} = 0;
1181         if ( $self -> {'covariance_step_run'} ) {
1182           $self -> {'covariance_step_successful'} = 1;
1183           # If there are problems in the cov-step they will be caught later
1184         }
1185         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
1186           $self -> {'s_matrix_singular'} = 1 if ( /^0S MATRIX ALGORITHMICALLY SINGULAR/ );
1187           if ( /^0R MATRIX ALGORITHMICALLY SINGULAR/ or
1188                /^0S MATRIX ALGORITHMICALLY SINGULAR/ ) {
1189             $self -> {'covariance_step_warnings'} = 1;
1190             next;
1191           }
1192           if ( /^0ESTIMATE OF THETA IS NEAR THE BOUNDARY AND/ or
1193                /0PARAMETER ESTIMATE IS NEAR ITS BOUNDARY/ ) {
1194             $self -> {'estimate_near_boundary'} = 1;
1195           }
1196           if ( /ROUNDING ERRORS/ ) {
1197             $self -> {'rounding_errors'} = 1;
1198           }
1199           if ( /0COVARIANCE STEP ABORTED/ ) {
1200             $self -> {'covariance_step_successful'} = 0;
1201           }
1202
1203           if ( /^0MINIMIZATION SUCCESSFUL/ ) {
1204             $self -> {'minimization_successful'} = 1;
1205             $success = 1;
1206           }
1207           if ( /MINIMUM VALUE OF OBJECTIVE FUNCTION/ ) {
1208             debug -> warn( level   => 2,
1209                            message => "Hmmm, reached the OFV area" );
1210             last;
1211           }
1212           $success_pos = $start_pos unless ( $success );
1213         }
1214         if ($success) {
1215           debug -> warn( level   => 2,
1216                          message => "Found a minimization statement" );
1217           $self -> {'lstfile_pos'} = $success_pos; #Back to successline.
1218         } else {
1219           debug -> warn( level   => 2,
1220                          message => "No minimization statement found" ); #Back to starting line
1221         }
1222       }
1223 end _read_term
1224
1225 # }}} _read_term
1226
1227 # {{{ _read_minimization_message
1228
1229 start _read_minimization_message
1230       {
1231         my $success  = 0;
1232         my (@mess, @etabar,@pval);
1233         my $termarea = 0;
1234
1235         my $start_pos = $self -> {'lstfile_pos'};
1236         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
1237           if ( /^ PROBLEM NO\.:\s*\d+\s*SUBPROBLEM NO\.:\s*(\d+)/) {
1238             last unless ($1 == 1); #Continuition of the same subproblem
1239           }
1240           last if ( / PROBLEM NO\.:\s*\d+\n/);
1241           if ( not $termarea and ( /(^0MINIMIZATION.*$)/ or /(^0PROGRAM TERMINATED.*$)/ ) ) {
1242             $termarea = 1;
1243             my $tmp = $self -> {'lstfile_pos'};
1244             debug -> warn( level   => 2,
1245                            message => "Setting firstpos from $tmp to $start_pos" );
1246             $self -> {'lstfile_pos'} = $start_pos;
1247             $tmp = $self -> {'lstfile_pos'};
1248             debug -> warn( level   => 2,
1249                            message => "0MINIMIZATION found! and firstpos is: $tmp" );
1250           }
1251           if ( $termarea and /^1/ ) {
1252             debug -> warn( level   => 2,
1253                            message => "Found minimization area and reached the \'\^1\'" );
1254             $termarea = 0;
1255             $start_pos = $self -> {'lstfile_pos'};
1256             $success = 1;
1257             last;
1258           }
1259           if ( /MINIMUM VALUE OF OBJECTIVE FUNCTION/ ) {
1260             debug -> warn( level   => 2,
1261                            message => "Hmmm, reached the OFV area" );
1262             last;
1263           }
1264           if ( $termarea ) {
1265             push @mess,$_ ;
1266           }
1267         }
1268
1269         push( @{$self -> {'minimization_message'}}, @mess );
1270
1271 #       if ( $self -> {'minimization_message'} ) { # PsN message in sub _read_iteration_path
1272 #         $self -> {'minimization_message'} .= join "", @mess;
1273 #       } else {
1274 #         $self -> {'minimization_message'} = join "", @mess;
1275 #       }
1276
1277         my @temp;
1278         for (@mess) {
1279           if (s/\s+ETABAR:\s+//) {
1280             @temp = split;
1281             push @etabar, @temp;
1282 #           print "Etabar Push of: $_ \n";
1283             for (@mess) {       #If wrapped etabar-lines
1284               last if(/^1/);
1285               last if(/a-zA-DF-Z/);
1286               last unless (s/^\s+//);
1287               last unless (/\d/);
1288               @temp = split;
1289               push @etabar, @temp;
1290             }
1291             last;
1292           }
1293         }
1294         # Initialize the attribute only if we have found any data
1295         $self -> {'etabar'} = \@etabar if ( $#etabar > 0 );
1296
1297
1298         for (@mess) {
1299           if (s/\s+P VAL\.:\s+//) {
1300             @temp = split;
1301             push @pval, @temp;
1302 #           print "Pval Push of: $_ \n";
1303             for (@mess) {       #If wrapped etabar-lines
1304               last if(/^1/);
1305               last if(/a-zA-DF-Z/);
1306               last unless (s/^\s+//);
1307               last unless (/\d/);
1308               @temp = split;
1309               push @pval, @temp;
1310             }
1311             last;
1312           }
1313         }
1314         $self -> {'pval'} = \@pval;
1315
1316         unless ( $success ) {
1317           debug -> warn( level   => 2,
1318                          message => "No minimization message found" );
1319         } else {
1320           $self -> {'lstfile_pos'} = $start_pos;
1321         }
1322       }
1323 end _read_minimization_message
1324
1325 # }}} _read_minimization_message
1326
1327 # {{{ _read_thomsi
1328
1329 start _read_thomsi
1330       {
1331         my $start_pos = $self -> {'lstfile_pos'};
1332         my $success = 0;
1333
1334         my $estarea      = 0;
1335         my $thetarea     = 0;
1336         my $omegarea     = 0;
1337         my $sigmarea     = 0;
1338         my ( @theta, @omega, @raw_omega, @sigma, @raw_sigma, @T, $i, $tmp );
1339
1340         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
1341           chomp;
1342           $estarea  = 1  if /FINAL PARAMETER ESTIMATE/;
1343           last if ( /STANDARD ERROR OF ESTIMATE/ and $estarea );
1344           last if ( /^1NONLINEAR/ and $estarea );
1345           if ( /THETA - VECTOR OF FIXED EFFECTS PARAMETERS/ and $estarea ) {
1346             debug -> warn( level   => 2,
1347                            message => "Found the theta estimates area" );
1348             $thetarea = 1;
1349             $success  = 1;
1350           }
1351           if ( /OMEGA - COV MATRIX FOR RANDOM EFFECTS - ETAS/ and $estarea ) {
1352             debug -> warn( level   => 2,
1353                            message => "Found the omega estimates area" );
1354             $omegarea = 1;
1355             $thetarea = 0;
1356             $sigmarea = 0;
1357           }
1358           if ( /SIGMA - COV MATRIX FOR RANDOM EFFECTS - EPSILONS/ and $estarea ) {
1359             debug -> warn( level   => 2,
1360                            message => "Found the sigma estimates area" );
1361             $sigmarea = 1;
1362             $thetarea = 0;
1363             $omegarea = 0;
1364           }
1365           last if ( /THERE ARE ERROR MESSAGES IN FILE PRDERR/ and $estarea );
1366           # The row below checking for a one '1', crashes behaviour with MANY etas (>15)
1367 #       last if ( /^1/ and $estarea );
1368           last if ( /^\s+\*\*\*\*\*\*\*\*\*\*\*\*/ and ( $sigmarea or $thetarea or $omegarea ));
1369           last if ( /^[A-W]/ and $estarea );
1370
1371           if( $thetarea and /^\s*-?\d*\.\d*/ ) {
1372             @T = split(' ',$_);
1373             for $i (0..(@T-1)) {
1374               $T[$i] = eval($T[$i]);
1375             }
1376             push(@theta,@T);
1377           }
1378
1379           if($omegarea and /^(\+|\s{2,})/) {
1380             next if /ET/;
1381             @T = split(' ',$_);
1382             shift @T if $T[0] eq '+';
1383             if ( $self -> {'omega_block_structure_type'} eq 'DIAGONAL' ) {
1384               for  $i (0..(@T-1)) {
1385                 if(($T[$i] ne '.........') and (eval($T[$i]) != 0) ) {
1386                   push(@omega,eval($T[$i]));
1387                 }
1388               }
1389             }
1390             for  $i (0..(@T-1)) {
1391               if($T[$i] ne '.........') {
1392                 $tmp = eval($T[$i]);
1393               } else {
1394                 $tmp = 'NA';
1395               }
1396               $T[$i] = $tmp ;
1397             }
1398             push(@raw_omega,@T);
1399           }
1400
1401           if($sigmarea and /^(\+|\s{2,})/) {
1402             next if /EP/;
1403             next if /^\s*$/;
1404             @T = split(' ',$_);
1405             shift @T if $T[0] eq '+';
1406             if ( $self -> {'sigma_block_structure_type'} eq 'DIAGONAL' ) {
1407               for  $i (0..(@T-1)) {
1408                 if(($T[$i] ne '.........') and (eval($T[$i]) != 0) ) {
1409                   push(@sigma,eval($T[$i]));
1410                 }
1411               }
1412             }
1413             for  $i (0..(@T-1)) {
1414               if($T[$i] ne '.........') {
1415                 $tmp = eval($T[$i]);
1416               } else {
1417                 $tmp = 'NA';
1418               }
1419               $T[$i] = $tmp ;
1420             }
1421             push(@raw_sigma,@T);
1422
1423           }
1424         }
1425
1426         if ( $self -> {'omega_block_structure_type'} eq 'BLOCK' ) {
1427           my @omblock  = @{$self -> {'omega_block_structure'}};
1428           debug -> warn( level   => 2,
1429                          message => "OMEGA BLOCK DEFINED" ) if ( scalar @omblock > 0 );
1430           @omega = ();
1431           my $i = 0;
1432           foreach my $row ( @omblock ) {
1433             foreach my $element ( @{$row} ) {
1434               push ( @omega, $raw_omega[$i] ) if ( $element );
1435               $i++;
1436             }
1437           }
1438         }
1439
1440         if ( $self -> {'sigma_block_structure_type'} eq 'BLOCK' ) {
1441           my @omblock  = @{$self -> {'sigma_block_structure'}};
1442           debug -> warn( level   => 2,
1443                          message => "SIGMA BLOCK DEFINED" ) if ( scalar @omblock > 0 );
1444           @sigma = ();
1445           my $i = 0;
1446           foreach my $row ( @omblock ) {
1447             foreach my $element ( @{$row} ) {
1448               push ( @sigma, $raw_sigma[$i] ) if ( $element );
1449               $i++;
1450             }
1451           }
1452         }
1453
1454         $self -> {'thetas'}       = \@theta;
1455         $self -> {'nth'}          = $#theta + 1;
1456         $self -> {'raw_omegas'} = \@raw_omega;
1457         $self -> {'nrom'} = $#raw_omega + 1;
1458         $self -> {'omegas'}       = \@omega;
1459         $self -> {'nom'} = $#omega + 1;
1460         $self -> {'raw_sigmas'} = \@raw_sigma;
1461         $self -> {'sigmas'}       = \@sigma;
1462
1463
1464         # Gather the "true" estimates, i.e. for the parameters that are not fixed or SAME.
1465         foreach my $param ( 'theta', 'omega', 'sigma' ) {
1466           my @allests = eval( '@'.$param );
1467           my @estflags = @{$self -> {'estimated_'.$param.'s'}};
1468           my @ests;
1469 #         print "$param\n";
1470 #         print Dumper \@allests;
1471 #         print Dumper \@estflags;
1472
1473           die "Something is wrong: All $param"."s: ".($#allests+1)." and estimated $param"."s: ".
1474               ($#estflags+1)." do not match\n" unless
1475               ( $#allests == -1 or $#estflags == $#allests );
1476           my $defs = 0;
1477           for( my $i = 0; $i <= $#allests; $i++ ) {
1478             if( $estflags[$i] ) {
1479               if ( defined $allests[$i] ) {
1480                 push( @ests, $allests[$i]);
1481                 $defs++;
1482               } else {
1483                 push( @ests, 0 );
1484               }
1485             }
1486           }
1487           $self -> {'est_'.$param.'s'} = \@ests;
1488
1489 #         if( $#ests > -1 and $defs > 0 ) {
1490 #           $self -> {'est_'.$param.'s'} = Math::MatrixReal ->
1491 #               new_from_cols( [\@ests] );
1492 #         }
1493         }
1494
1495         unless ( $success ) {
1496           debug -> warn( level   => 2,
1497                          message => "No thetas, omegas or sigmas found" );
1498         } else {
1499 # Keep this code
1500 #         if ( $PsN::config -> {'_'} -> {'use_database'} and
1501 #              $self -> {'register_in_database'} ) {
1502 #           my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
1503 #                              ";databse=".$PsN::config -> {'_'} -> {'project'},
1504 #                                    $PsN::config -> {'_'} -> {'user'},
1505 #                               $PsN::config -> {'_'} -> {'password'},
1506 #                                    {'RaiseError' => 1});
1507 #           my $sth;
1508 #           my @mod_str = ('','');
1509 #           if ( defined $self -> {'model_id'} ) {
1510 #             @mod_str = ('model_id,',"$self->{'model_id'},");
1511 #           }
1512 #           foreach my $param ( 'theta', 'omega', 'sigma' ) {
1513 #             my $i = 1;
1514 #             foreach my $par_str ( eval('@'.$param) ) {
1515 #               $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
1516 #                              ".estimate ".
1517 #                                      "(subproblem_id,problem_id,output_id,".
1518 #                                      $mod_str[0].
1519 #                                      "type,value, number, label) ".
1520 #                                      "VALUES ( '$self->{'subproblem_id'}' ,".
1521 #                                      "'$self->{'problem_id'}' ,".
1522 #                                      "'$self->{'output_id'}' ,".
1523 #                                      $mod_str[1].
1524 #                                      "'$param','$par_str', '$i', 'test_label')");
1525 #               $sth -> execute;
1526 #               push( @{$self -> {'estimate_ids'}}, $sth->{'mysql_insertid'} );
1527 #               $i++;
1528 #             }
1529 #           }
1530 #           $sth -> finish;
1531 #           $dbh -> disconnect;
1532 #         }
1533
1534
1535           $self -> {'lstfile_pos'} = $start_pos-1;
1536         }
1537       }
1538 end _read_thomsi
1539
1540 # }}} _read_thomsi
1541
1542 # {{{ _set_omeganames
1543 start _set_omeganames
1544       {
1545         my @raw_omegas = @{$self -> {'raw_omegas'}};
1546
1547         unless ( scalar @raw_omegas > 0 ) {
1548           return undef;
1549         }
1550         ;
1551
1552         my ( @omeganames, %omeganameval, @omegas );
1553         my ($j,$ndiags);
1554         my @indices;
1555         foreach $j (1..scalar @raw_omegas) {
1556           push @omeganames, "OM".++$ndiags and next if $self -> _isdiagonal('index' => $j);
1557           if ($raw_omegas[$j-1] !=0) {
1558             @indices = $self -> _rowcolind( index => $j);
1559             push @omeganames,"OM".$indices[0]."\.".$indices[1];
1560           }
1561         }
1562         @omegas = @{$self -> {'omegas'}};
1563         for ( my $i = 0; $i <= $#omeganames; $i++ ) {
1564           $omeganameval{$omeganames[$i]} = $omegas[$i];
1565         }
1566         $self ->{'omeganameval'} = \%omeganameval;
1567         $self ->{'omeganames'} = \@omeganames;
1568       }
1569 end _set_omeganames
1570 # }}} _set_omeganames
1571
1572 # {{{ _set_sigmanames
1573 start _set_sigmanames
1574       {
1575         my @sigmas = @{$self -> {'raw_sigmas'}};
1576         unless ( scalar @sigmas > 0 ) {
1577           return undef;
1578         }
1579
1580         my ( @sigmanames, %sigmanameval );
1581         my ($j,$ndiags);
1582         my @indices;
1583         foreach $j (1..scalar @sigmas) {
1584           push @sigmanames, "SI".++$ndiags and next if $self -> _isdiagonal('index' => $j);
1585           if ($sigmas[$j-1] !=0) {
1586             @indices = $self -> _rowcolind( index => $j);
1587             push @sigmanames,"SI".$indices[0].$indices[1];
1588           }
1589         }
1590         @sigmas = @{$self -> {'sigmas'}};
1591         for ( my $i = 0; $i <= $#sigmanames; $i++ ) {
1592           $sigmanameval{$sigmanames[$i]} = $sigmas[$i];
1593         }
1594         $self ->{'sigmanameval'} = \%sigmanameval;
1595         $self ->{'sigmanames'} = \@sigmanames;
1596       }
1597 end _set_sigmanames
1598 # }}} _set_sigmanames
1599
1600 # {{{ _set_thetanames
1601 start _set_thetanames
1602       {
1603         my $nth = $self -> {'nth'};
1604         my ( @thetanames, %thetanameval, @thetas );
1605
1606         for (1..$nth) {
1607           push @thetanames, "TH$_";
1608         }
1609         @thetas = @{$self -> {'thetas'}};
1610         for ( my $i = 0; $i <= $#thetanames; $i++ ) {
1611           $thetanameval{$thetanames[$i]} = $thetas[$i];
1612         }
1613         $self ->{'thetanameval'} = \%thetanameval;
1614         $self ->{'thetanames'} = \@thetanames;
1615       }
1616 end _set_thetanames
1617 # }}} _set_thetanames
1618
1619