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