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       }
  70 end new
  71
  72 # }}} new
  73
  74 # {{{ parsing_error
  75 start parsing_error
  76     $self -> parsed_successfully( 0 );
  77 $self -> parsing_error_message( $message );
  78 end parsing_error
  79 # }}} parsing_error
  80
  81 # {{{ register_in_database
  82
  83 start register_in_database
  84     if ( $PsN::config -> {'_'} -> {'use_database'} ) {
  85       my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
  86                                ";databse=".$PsN::config -> {'_'} -> {'project'},
  87                                $PsN::config -> {'_'} -> {'user'},
  88                                 $PsN::config -> {'_'} -> {'password'},
  89                                {'RaiseError' => 1});
  90       my $sth;
  91       my $term_str;
  92       if ( defined $self->{'minimization_message'} ) {
  93         $term_str = join("\n",@{$self->{'minimization_message'}});
  94       }
  95       my @mod_str = ('','');
  96       if ( defined $model_id ) {
  97         @mod_str = ('model_id,',"$model_id,");
  98       }
  99       $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
 100                                ".subproblem ".
 101                              "(problem_id,output_id,".
 102                              $mod_str[0].
 103                              "condition_number,".
 104                              "covariance_step_successful,".
 105                              "covariance_step_warnings,".
 106                              "eval,".
 107                              "final_zero_gradients,".
 108                              "hessian_reset,".
 109                              "ofv,".
 110                              "rounding_errors,".
 111                              "s_matrix_singular,".
 112                              "significant_digits,".
 113                              "simulation,".
 114                              "minimization_successful,".
 115                              "minimization_message,".
 116                              "zero_gradients ) ".
 117                              "VALUES ( '$problem_id' ,".
 118                              "'$output_id' ,".
 119                              $mod_str[1].
 120                              "'$self->{'condition_number'}' ,".
 121                              "'$self->{'covariance_step_successful'}' ,".
 122                              "'$self->{'covariance_step_warnings'}' ,".
 123                              "'$self->{'eval'}' ,".
 124                              "'$self->{'final_zero_gradients'}' ,".
 125                              "'$self->{'hessian_reset'}' ,".
 126                              "'$self->{'ofv'}' ,".
 127                              "'$self->{'rounding_errors'}' ,".
 128                              "'$self->{'s_matrix_singular'}' ,".
 129                              "'$self->{'significant_digits'}' ,".
 130                              "'$self->{'simulation'}' ,".
 131                              "'$self->{'minimization_successful'}' ,".
 132                              "'$term_str' ,".
 133                              "'$self->{'zero_gradients'}' )");
 134       $sth -> execute;
 135       $self -> {'subproblem_id'} = $sth->{'mysql_insertid'};
 136       $self -> {'problem_id'} = $problem_id;
 137       $self -> {'output_id'} = $output_id;
 138       $self -> {'model_id'} = $model_id;
 139       foreach my $param ( 'theta', 'omega', 'sigma' ) {
 140         my $i = 1;
 141         foreach my $par_str ( @{$self -> {$param.'s'}} ) {
 142           $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
 143                                  ".estimate ".
 144                                  "(subproblem_id,problem_id,output_id,".
 145                                  $mod_str[0].
 146                                  "type,value, number, label) ".
 147                                  "VALUES ( '$self->{'subproblem_id'}' ,".
 148                                  "'$self->{'problem_id'}' ,".
 149                                  "'$self->{'output_id'}' ,".
 150                                  $mod_str[1].
 151                                  "'$param','$par_str', '$i', 'test_label')");
 152           $sth -> execute;
 153           push( @{$self -> {'estimate_ids'}}, $sth->{'mysql_insertid'} );
 154           $i++;
 155         }
 156       }
 157       $sth -> finish;
 158       $dbh -> disconnect;
 159     }
 160
 161 # if ( $PsN::config -> {'_'} -> {'use_database'} ) {
 162 #   my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
 163 #                              ";databse=".$PsN::config -> {'_'} -> {'project'},
 164 #                          $PsN::config -> {'_'} -> {'user'},
 165 #                               $PsN::config -> {'_'} -> {'password'},
 166 #                          {'RaiseError' => 1});
 167 #   my $sth = $dbh -> prepare( "UPDATE ".$PsN::config -> {'_'} -> {'project'}.
 168 #                              ".subproblem ".
 169 #                            "SET condition_number='".
 170 #                            $self -> {'condition_number'}."'".
 171 #                            "WHERE model_id='".$self -> {'model_id'}."'" );
 172 #   $sth -> execute or debug -> die( message => $sth->errstr ) ;
 173 #   $sth -> finish;
 174 #   $dbh -> disconnect;
 175 # }
 176 end register_in_database
 177
 178 # }}} register_in_database
 179
 180 # {{{ _compute_comegas_or_csigmas
 181
 182 start _compute_comegas_or_csigmas
 183 {
 184   # This method transforms omegas or sigmas.
 185
 186   # First we check if we have omegas or sigmas.
 187   if ( defined $self -> {'raw_' . $omega_or_sigma . 's'} ) {
 188     my @raw_omegas_or_sigmas = @{$self -> {'raw_' . $omega_or_sigma . 's'}};
 189     # If the omega or sigma matrix has no off-diagonal elements, the
 190     # transformation is quite straightforward.
 191     if( $self -> { $omega_or_sigma . '_block_structure_type'} eq 'DIAGONAL' ){
 192
 193       # All diagonal omegas or sigmas will be in
 194       # @raw_omegas_or_sigmas, so we loop over it and take the square
 195       # root of non-negative, non-zero numbers (those become
 196       # undefined).
 197       my $row_size = 1;
 198       for( my $i = 0; $i <= $#raw_omegas_or_sigmas; $i+=$row_size ){
 199         my $omega_or_sigma_value = undef;
 200         if( $raw_omegas_or_sigmas[$i] >= 0 ) {
 201           $omega_or_sigma_value = sqrt( $raw_omegas_or_sigmas[$i] );
 202         } else {
 203           ui -> print( category => 'all',
 204                        message  => "Warning: cannot take the square root of ".
 205                        "$omega_or_sigma with value ".$raw_omegas_or_sigmas[$i] );
 206         }
 207         push( @{$self -> {'c' . $omega_or_sigma .'s'}}, $omega_or_sigma_value );
 208         $row_size ++;
 209       }
 210
 211     } elsif( $self -> {$omega_or_sigma . '_block_structure_type'} eq 'BLOCK' ) {
 212
 213       # If the omega or sigma matrix has block form, the transformation
 214       # becomes a bit more complex. The transformation for diagonal
 215       # elements is the same. But for off-diagonals it takes more into
 216       # account.
 217
 218       # $current_row_size is the number of omegas or sigmas in the row
 219       # in the matrix that we are currently transforming.
 220       my $current_row_size = 1;
 221
 222       # $diagonal_index is the index in @raw_omegas_or_sigmas where the current
 223       # on-diagonal omega or sigma we are transforming.
 224       my $diagonal_index = 0;
 225
 226       # $y is the y-coordinate of the element in the matrix we
 227       # are transforming.
 228       my $y = 0;
 229
 230       # First get all block sets sorted.
 231       my @keys = sort {$a <=> $b} keys( %{$self -> { $omega_or_sigma . '_block_sets'}} );
 232
 233       foreach my $key ( @keys ) {
 234
 235         # If this is a block set $block_dimension will be > 1
 236         # otherwise it will be 1.
 237         my $block_dimension = $self -> { $omega_or_sigma . '_block_sets'}{$key}{'dimension'};
 238         my $block_size = $self -> { $omega_or_sigma . '_block_sets'}{$key}{'size'};
 239
 240 #       print "BS: $block_size\n";
 241 #       print "BD: $block_dimension\n";
 242
 243         # loop over the repeated blocks
 244         for( my $i = 1; $i <= $block_size; $i++ ) {
 245           # loop over the size of the block (i.e. the X in block(X))
 246           for( my $j = 1; $j <= $block_dimension; $j++ ) {
 247             $y++;
 248
 249             # The equation y*(y+1)/2 is a geometric sum that gives the
 250             # index of the diagonal element on row y in the compacted
 251             # matrixarray @raw_omegas_or_sigmas. -1 to get array
 252             # index.
 253             my $y_y_idx = ($y*($y+1)/2)-1;
 254
 255             my $start_x = ($y - $j + 1);
 256             for( my $x = $start_x; $x <= $y; $x++ ) {
 257 #             print "X: $x, Y: $y\n";
 258               my $omega_or_sigma_value;
 259               if( $i == 1 ) {
 260                 # This is the first of the repeated blocks, collect
 261                 # the values (otherwise, just leave
 262                 # omega_or_sigma_value undefined)
 263                 if( $x == $y ) {
 264
 265                   # If we are on the diagonal, just take the square root
 266                   # of the diagonal element.
 267                   $omega_or_sigma_value = $raw_omegas_or_sigmas[ $y_y_idx ];
 268                   $omega_or_sigma_value = sqrt( $omega_or_sigma_value );
 269 #                 print "y_y_idx: $y_y_idx\n";
 270
 271                 } else {
 272                   # If we are off the diagonal, we need to find two
 273                   # on-diagonal omegas, one on the same column and one on
 274                   # the same row.
 275
 276                   my $x_x_idx = ($x*($x+1)/2)-1; # See explanation for $y_y_idx above
 277
 278                   # Off diagonal index: $y_y_idx for the row before plus $x
 279                   my $x_y_idx = (($y-1)*$y/2+$x)-1;
 280
 281                   $omega_or_sigma_value = $raw_omegas_or_sigmas[$x_y_idx];
 282                   my $denominator = $raw_omegas_or_sigmas[ $x_x_idx ]
 283                       * $raw_omegas_or_sigmas[ $y_y_idx ];
 284                   if( $denominator <= 0 ){ # To avoiding division by zero
 285                     $omega_or_sigma_value = undef;
 286                   } elsif( $omega_or_sigma_value >= sqrt($denominator) ) {
 287                     # This rounding handles cases when the offdiagonals
 288                     # are greater or equal to one.
 289                     $omega_or_sigma_value = $omega_or_sigma_value/( int( 10000 * sqrt($denominator) )/10000 );
 290                   } else {
 291                     $omega_or_sigma_value = $omega_or_sigma_value/sqrt($denominator);
 292                   }
 293                 }
 294               }
 295               push( @{$self -> {'c' . $omega_or_sigma . 's'}}, $omega_or_sigma_value );
 296
 297             }
 298           }
 299         }
 300       }
 301 #       my $loop_limit = $block_dimension > $block_size ? $block_dimension : $block_size;
 302 #       print "loop limit: $loop_limit\n";
 303 #       for( my $i = 0; $i < $loop_limit; $i ++ ){
 304 #         # x_idx is the x-coordinate of the omega or sigma we are
 305 #         # transforming.
 306
 307 #         my $x_idx = $y_idx - $i;
 308 #         # current row size is equal to the y coordinate
 309 # #       my $x_idx = $current_row_size - $i;
 310
 311 #         # The innner_loop_limit must be 1 when we have SAME blocks,
 312 #         # which can be determined from the block_size variable.
 313 #         my $inner_loop_limit = $block_size > 1 ? 0 : $i;
 314 #         print "inner loop limit: $inner_loop_limit\n";
 315 #         # The inner loop is over the row elements in the block
 316 #         # set. (1 when we have SAME blocks)
 317 #         for( my $j = 0; $j <= $inner_loop_limit; $j++ ){
 318
 319 #           # $omega_or_sigma will hold the transformed value.
 320 #           my $omega_or_sigma_value;
 321
 322 #           print "X: $x_idx, Y: $y_idx\n";
 323 #           if( $x_idx == $y_idx ) {
 324
 325 #             # If we are on the diagonal, just take the square root
 326 #             # of the diagonal element.
 327 #             $omega_or_sigma_value = $raw_omegas_or_sigmas[ $diagonal_index ];
 328 #             $omega_or_sigma_value = sqrt( $omega_or_sigma_value );
 329
 330 #           } else {
 331 #             # If we are off the diagonal, we need to find two
 332 #             # on-diagonal omegas, one on the same column and one on
 333 #             # the same row.
 334
 335 #             # The equation y*(y+1)/2 is a geometric sum that gives
 336 #             # the index of the diagonal element on row y in the
 337 #             # compacted matrixarray @raw_omegas_or_sigmas. The equation
 338 #             # x*(x+1)/2 gives the diagonal element on column x in
 339 #             # the same array.
 340
 341 #             # We find the omega by finding the index of the diagonal
 342 #             # element on row y and subtracting the difference
 343 #             # between the number of elements on the row ($i) and the
 344 #             # number of elements we have transformed ($j).
 345
 346 #             $omega_or_sigma_value = $raw_omegas_or_sigmas[(($y_idx*($y_idx+1))/2)-1 - ($i - $j)];
 347 #             my $denominator = $raw_omegas_or_sigmas[ (($x_idx*($x_idx+1))/2)-1 ]
 348 #                                *
 349 #                               $raw_omegas_or_sigmas[ (($y_idx*($y_idx+1))/2)-1 ];
 350 #             if( $denominator <= 0 ){ # To avoiding division by zero
 351 #               $omega_or_sigma_value = undef;
 352 #             } elsif( $omega_or_sigma_value >= sqrt($denominator) ) {
 353 #               # This rounding handles cases when the offdiagonals
 354 #               # are greater or equal to one.
 355 #               $omega_or_sigma_value = $omega_or_sigma_value/( int( 10000 * sqrt($denominator) )/10000 )
 356 #             } else {
 357 #               $omega_or_sigma_value = $omega_or_sigma_value/sqrt($denominator);
 358 #             }
 359 #           }
 360 #           push( @{$self -> {'c' . $omega_or_sigma . 's'}}, $omega_or_sigma_value );
 361
 362 #           # Move the x-coordinate forwards.
 363 #           $x_idx ++;
 364 #         }
 365
 366 #         # Move the y-coordinate forwards.
 367 #         $y_idx ++;
 368 #         # Calculate the new row size.
 369 # #       $current_row_size ++;
 370 #         # Skip one row and we get the index of the next diagonal.
 371
 372 #         $diagonal_index += $y_idx;
 373 # #       $diagonal_index += $current_row_size;
 374 #       }
 375 #     }
 376     } else {
 377       'debug' -> warn( level => 1,
 378                        message => $omega_or_sigma . ' matrix has unknown form.' );
 379     }
 380   }
 381 }
 382 end _compute_comegas_or_csigmas
 383
 384 # }}} _compute_comegas_or_sigmas
 385
 386 # {{{ _compute_comegas
 387 start _compute_comegas
 388 {
 389   $self -> _compute_comegas_or_csigmas( omega_or_sigma => 'omega' );
 390 }
 391 end _compute_comegas
 392
 393 # This is the old algoritm for transforming omegas. It can't handle
 394 # of-diagonals omegas that are zero, if that should ever happen.
 395
 396 #     my $i;
 397 #     my @raw_omegas;
 398 #     my @computed_omegas;
 399 #     my @diags;
 400 #     my @indices;
 401
 402 #     @raw_omegas = @{$self -> {'raw_omegas'}};
 403 #     @computed_omegas = ();
 404
 405 #     if( $self -> {'omega_structure_type'} eq 'DIAGONAL' ){
 406 #       @diags = @raw_omegas;
 407 #     } else {
 408 #       ## Collect the diagonals
 409 #       foreach $i (0..$#raw_omegas) {
 410 #       if( $self -> _isdiagonal('index' => $i+1) ){
 411 #         push( @diags, $raw_omegas[$i] );
 412 #       }
 413 #       }
 414 #     }
 415
 416 #     foreach $i (0..$#raw_omegas){
 417 #       if($self -> _isdiagonal('index' => $i+1)) {
 418 #       if( $raw_omegas[$i] <= 0 ){
 419 #         push( @computed_omegas, $raw_omegas[$i] );
 420 #         next;
 421 #       }
 422 #       push @computed_omegas, sqrt($raw_omegas[$i]);
 423 #       next;
 424 #       } else {
 425 #       @indices = $self -> _rowcolind( index => $i+1 );
 426 #       if( ($raw_omegas[$i] == 0) or ($diags[$indices[0]-1]*$diags[$indices[1]-1]) <= 0 ) {
 427 #         next;
 428 #       } else {
 429 #         push @computed_omegas,
 430 #         $raw_omegas[$i]/sqrt($diags[$indices[0]-1]*$diags[$indices[1]-1]);
 431 #       }
 432 #       }
 433 #     }
 434 #   $self -> {'comegas'}= \@computed_omegas;}
 435
 436 # }}}
 437
 438 # {{{ _compute_csigmas
 439 start _compute_csigmas
 440 {
 441   $self -> _compute_comegas_or_csigmas( omega_or_sigma => 'sigma' );
 442 }
 443 end _compute_csigmas
 444 # This is the old algoritm for transforming sigmass. It can't handle
 445 # of-diagonals sigmass that are zero, if that should ever happen.
 446
 447 # {
 448 #     if (defined $self -> {'raw_sigmas'} ) {
 449 #       my $i;
 450 #       my ( @om, @com, @diags, @indices );
 451 #       @om = @{$self -> {'raw_sigmas'}};
 452 #       @com = ();
 453 #       ## Collect the diagonals
 454 #       foreach $i (0..$#om) {
 455 #           push @diags,$om[$i] if $self -> _isdiagonal('index' => $i+1);
 456 #       }
 457 #       foreach $i (0..$#om){
 458 #           if($self -> _isdiagonal('index' => $i+1)) {
 459 #             if( $om[$i] <= 0 ){
 460 #               next;
 461 #             }
 462 #             push @com, sqrt($om[$i]);
 463 #             next;
 464 #           } else {
 465 #             @indices = $self -> _rowcolind( index => $i+1 );
 466 #             if( ($om[$i] == 0) or ($diags[$indices[0]-1]*$diags[$indices[1]-1] <= 0) ) {
 467 #               next;
 468 #             } else {
 469 #               push @com,
 470 #               $om[$i]/sqrt($diags[$indices[0]-1]*$diags[$indices[1]-1]);
 471 #             }
 472 #           }
 473 #       }
 474 #       $self -> {'csigmas'}= \@com;
 475 #     }
 476 # }
 477
 478 # }}} _compute_csigmas
 479
 480 # {{{ _compute_cvseomega
 481
 482 start _compute_cvseomega
 483 {
 484   if( defined $self -> raw_omegas() and defined $self -> raw_seomegas() ) {
 485     my @raw_omegas   = @{$self -> raw_omegas()};
 486     my @raw_seomegas = @{$self -> raw_seomegas()};
 487     my @cvseomega;
 488
 489     if ( scalar @raw_seomegas > 0) {
 490       my ($j,$i);
 491       my ( @init_om, @init_seom );
 492
 493       foreach $i (0.. $#raw_seomegas){
 494         if( $self -> _isdiagonal('index' => $i+1)) {
 495           push @init_seom, $raw_seomegas[$i];
 496           next;
 497         } else {
 498           if(($raw_seomegas[$i] eq 'NA') or ($raw_seomegas[$i] == 0)) {
 499             next;
 500           } else {
 501             push @init_seom, $raw_seomegas[$i];
 502           }
 503         }
 504       }
 505
 506       foreach $i (0..$#raw_omegas){
 507         if( $self -> _isdiagonal('index' => $i+1)) {
 508           push @init_om, $raw_omegas[$i];
 509           next;
 510             } else {
 511               if($raw_omegas[$i] == 0) {
 512                 next;
 513               } else {
 514                 push @init_om, $raw_omegas[$i];
 515               }
 516             }
 517       }
 518
 519       foreach my $i (0..$#init_om) {
 520         if( $init_seom[$i] ne 'NA' and $init_om[$i] ne 'NA' ){
 521           if( ($init_seom[$i] == 'INF' or $init_seom[$i] == '-INF') and
 522               ($init_om[$i] == 'INF' or $init_om[$i] == '-INF') ) {
 523             push @cvseomega,undef;
 524           } elsif ( $init_seom[$i] == 'INF' ) {
 525             push @cvseomega,'INF';
 526           } elsif ( $init_seom[$i] == '-INF' ) {
 527             push @cvseomega,'-INF';
 528           } elsif ( $init_om[$i] == 'INF' or
 529                     $init_om[$i] == '-INF' ) {
 530             push @cvseomega,0;
 531           } else {
 532             push @cvseomega,$init_seom[$i]/abs($init_om[$i]);
 533           }
 534         } else {
 535           push @cvseomega,undef;
 536         }
 537 #           push @cvseomega,$init_seom[$i]/$init_om[$i]
 538 #               if ($init_seom[$i] && $init_om[$i]);
 539 #           push @cvseomega,"NA" unless ($init_seom[$i] && $init_om[$i]);
 540       }
 541
 542     }
 543
 544     $self -> cvseomegas([@cvseomega]);
 545   }
 546 }
 547 end _compute_cvseomega
 548
 549 # }}} _compute_cvseomega
 550
 551 # {{{ _compute_cvsesigma
 552 start _compute_cvsesigma
 553 {
 554   if( defined $self -> raw_sigmas() and defined $self -> raw_sesigmas() ) {
 555     my @raw_sigmas   = @{$self -> raw_sigmas()};
 556     my @raw_sesigmas = @{$self -> raw_sesigmas()};
 557     my @cvsesigma;
 558
 559     if ( scalar @raw_sesigmas > 0) {
 560       my ($j,$i);
 561       my ( @init_si, @init_sesi );
 562
 563       foreach $i (0.. $#raw_sesigmas){
 564         if($self -> _isdiagonal('index' => $i+1)) {
 565           push @init_sesi, $raw_sesigmas[$i];
 566           next;
 567         } else {
 568           if(($raw_sesigmas[$i] eq 'NA') or ($raw_sesigmas[$i] == 0)) {
 569             next;
 570           } else {
 571             push @init_sesi, $raw_sesigmas[$i];
 572           }
 573         }
 574       }
 575
 576       foreach $i (0.. $#raw_sigmas){
 577         if($self -> _isdiagonal('index' => $i+1)) {
 578           push @init_si, $raw_sigmas[$i];
 579           next;
 580         } else {
 581           if($raw_sigmas[$i] == 0) {
 582             next;
 583           } else {
 584             push @init_si, $raw_sigmas[$i];
 585           }
 586         }
 587       }
 588
 589       foreach my $i (0..$#init_si) {
 590         if( $init_sesi[$i] ne 'NA' and $init_si[$i] ne 'NA' ){
 591           if( ($init_sesi[$i] == 'INF' or $init_sesi[$i] == '-INF') and
 592               ($init_si[$i] == 'INF' or $init_si[$i] == '-INF') ) {
 593             push @cvsesigma,undef;
 594           } elsif ( $init_sesi[$i] == 'INF' ) {
 595             push @cvsesigma,'INF';
 596           } elsif ( $init_sesi[$i] == '-INF' ) {
 597             push @cvsesigma,'-INF';
 598           } elsif ( $init_si[$i] == 'INF' ) {
 599             push @cvsesigma,0;
 600           } else {
 601             push @cvsesigma,$init_sesi[$i]/abs($init_si[$i]);
 602           }
 603         } else {
 604           push @cvsesigma,undef;
 605         }
 606       }
 607
 608     }
 609
 610     $self -> cvsesigmas([@cvsesigma]);
 611   }
 612 }
 613 end _compute_cvsesigma
 614 # }}} _compute_cvsesigma
 615
 616 # {{{ _compute_cvsetheta
 617 start _compute_cvsetheta
 618 {
 619   if( defined $self -> thetas() and defined $self -> sethetas() ) {
 620     my @thetas   = @{$self -> thetas()};
 621     my @sethetas = @{$self -> sethetas()};
 622     my @cvsethetas;
 623
 624     if ( scalar @sethetas > 0 ) {
 625       foreach my $i (0..$#thetas) {
 626         if( defined $sethetas[$i] and
 627             defined $thetas[$i] ) {
 628           if( $thetas[$i] != 0 ) {
 629             push(@cvsethetas,$sethetas[$i]/abs($thetas[$i]));
 630           } elsif( $sethetas[$i] > 0 ) {
 631             push(@cvsethetas,'INF');
 632           } else {
 633             push(@cvsethetas,'-INF');
 634           }
 635         } else {
 636           push @cvsethetas,"NA";
 637         }
 638       }
 639     }
 640
 641     $self -> cvsethetas([@cvsethetas]);
 642   }
 643 }
 644 end _compute_cvsetheta
 645 # }}} _compute_cvsetheta
 646
 647 # {{{ _isdiagonal
 648 start _isdiagonal
 649 {
 650     my $previ = 1;
 651     my $j;
 652     return(1) if $index == 1;
 653     foreach my $j (2..100) {
 654         return(1) if $index == $previ+$j;
 655         $previ = $previ+$j;
 656         last if $index < $previ;
 657     }
 658 }
 659 end _isdiagonal
 660 # }}} _isdiagonal
 661
 662 # {{{ _read_matrixoestimates
 663 start _read_matrixoestimates
 664       {
 665         # Reads one matrix structure and returns the file handle at
 666         # the beginning of the next structure
 667         while ( $_ = @{$self -> {'lstfile'}}[ $pos++ ] ) {
 668           last if (/^\s*\*/);
 669           # Rewind one step if we find something that marks the end of
 670           # our structure
 671           $pos-- and last if ( /^ PROBLEM.*SUBPROBLEM/ or /^ PROBLEM NO\.:\s+\d/ );
 672           $pos-- and last if (/^[a-df-zA-DF-Z]/);
 673
 674           next if ( /^\s+TH/ or /^\s+OM/ or /^\s+SG/ );   # Header row
 675           next if ( /^1/ );                       # Those annoying 1's
 676
 677           chomp;                                # Get rid of line-feed
 678           my @row = split;
 679           shift( @row ) if ( $row[0] eq '+' );     # Get rid of +-sign
 680
 681           next if ( $#row < 0 );                           # Blank row
 682
 683           push( @subprob_matrix, @row );
 684         }
 685         $success = 1 if( scalar @subprob_matrix > 0 );
 686       }
 687 end _read_matrixoestimates
 688 # }}} _read_matrixoestimates
 689
 690 # {{{ _read_covmatrix
 691 start _read_covmatrix
 692       {
 693         my $matrix_nr = 0;
 694         my ( $t_success, $c_success, $corr_success, $i_success ) = (0,0,0,0);
 695         my $start_pos = $self -> {'lstfile_pos'}-1;
 696
 697         # {{{ sub clear dots
 698
 699         sub clear_dots {
 700           my $m_ref = shift;
 701           my @matrix = @{$m_ref};
 702           # get rid of '........'
 703           my @clear;
 704           foreach ( @matrix ) {
 705             unless ( $_ eq '.........' ) {
 706               if ( $_ eq 'NAN' ) {
 707                 push( @clear, 0 );
 708               } elsif ( $_ eq 'INF' ) {
 709                 push( @clear, 99999999999999 );
 710               } elsif ( $_ eq '-INF' ) {
 711                 push( @clear, -99999999999999 );
 712               } else {
 713                 push( @clear, $_ );
 714               }
 715             }
 716           }
 717           return \@clear;
 718         }
 719
 720         # }}}
 721
 722         # {{{ sub make square
 723
 724         sub make_square {
 725           my $m_ref = shift;
 726           my @matrix = @{$m_ref};
 727           # Make the matrix square:
 728           my $elements = scalar @matrix; # = M*(M+1)/2
 729           my $M = -0.5 + sqrt( 0.25 + 2 * $elements );
 730           my @square;
 731           for ( my $m = 1; $m <= $M; $m++ ) {
 732             for ( my $n = 1; $n <= $m; $n++ ) {
 733               push( @{$square[$m-1]}, $matrix[($m-1)*$m/2 + $n - 1] );
 734               unless ( $m == $n ) {
 735                 push( @{$square[$n-1]}, $matrix[($m-1)*$m/2 + $n - 1] );
 736               }
 737             }
 738           }
 739           return \@square;
 740         }
 741
 742         # }}}
 743
 744         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 745           if (/T MATRIX/) {
 746             while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 747               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 748                 ( $start_pos, $self -> {'raw_tmatrix'}, $t_success )  = $self ->
 749                     _read_matrixoestimates( pos => $start_pos ) and last;
 750               }
 751             }
 752             last;                # No covariance matrix will be found!
 753           }
 754           if (/    COVARIANCE MATRIX OF ESTIMATE/) {
 755             while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 756               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 757                 ( $start_pos, $self -> {'raw_covmatrix'}, $c_success ) = $self ->
 758                     _read_matrixoestimates( pos => $start_pos ) and last;
 759               }
 760             }
 761           }
 762           if (/    CORRELATION MATRIX OF ESTIMATE/) {
 763             while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 764               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 765                 ( $start_pos, $self -> {'raw_cormatrix'}, $corr_success ) = $self ->
 766                     _read_matrixoestimates( pos => $start_pos ) and last;
 767               }
 768             }
 769           }
 770           if (/    INVERSE COVARIANCE MATRIX OF ESTIMATE/) {
 771             while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 772               if (/^ TH (\d)/ or /^\s+TH (\d) \| /) { # Read matrix and get out of inner while loop
 773                 ( $start_pos, $self -> {'raw_invcovmatrix'}, $i_success ) = $self ->
 774                     _read_matrixoestimates( pos => $start_pos ) and last;
 775               }
 776             }
 777             last;                                       # Last matrix?
 778           }
 779         }
 780         foreach my $element ( @{$self -> {'raw_tmatrix'}} ) {
 781           push( @{$self -> {'t_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 782         }
 783         foreach my $element ( @{$self -> {'raw_covmatrix'}} ) {
 784           push( @{$self -> {'covariance_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 785         }
 786         foreach my $element ( @{$self -> {'raw_cormatrix'}} ) {
 787           push( @{$self -> {'correlation_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 788         }
 789
 790         if( defined $self -> {'raw_invcovmatrix'} ) {
 791             $self -> {'inverse_covariance_matrix'} = Math::MatrixReal ->
 792                 new_from_cols( make_square( clear_dots( $self -> {'raw_invcovmatrix'} ) ) );
 793         }
 794
 795 #       foreach my $element ( @{$self -> {'raw_invcovmatrix'}} ) {
 796 #         push( @{$self -> {'inverse_covariance_matrix'}}, eval($element) ) unless ( $element eq '.........' );
 797 #       }
 798
 799         #If something has gone right!
 800         $self-> {'lstfile_pos'} = $start_pos if ( $t_success + $c_success + $corr_success + $i_success  );
 801       }
 802 end _read_covmatrix
 803 # }}} _read_covmatrix
 804
 805 # {{{ _read_eigen
 806
 807 start _read_eigen
 808       {
 809         my @eigens;
 810         my $start_pos = $self -> {'lstfile_pos'};
 811         my $eig_area  = 0;
 812         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 813           chomp;
 814           if ( /EIGENVALUES OF COR MATRIX OF ESTIMATE/ ) {
 815             $eig_area = 1;
 816             $start_pos = $start_pos + 4 ; # Jump forward to the index numbers
 817             debug -> warn( level   => 2,
 818                            message => "Found the eigenvalue area" );
 819           INNER: while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) { # Get rid of indexes
 820               last INNER if ( not /^\s+\d/ );
 821 #             $start_pos++ and last INNER if ( not /^\s+\d/ );
 822             }
 823           }
 824           if ( $eig_area ) {
 825             $start_pos-- and last if (/^[a-df-zA-DF-Z]/); #Rewind one step
 826             last if ( /^\s*\*/ or /^1/ );
 827             push( @eigens, split );
 828           }
 829           $start_pos-- and last if ( /^ PROBLEM.*SUBPROBLEM/ or /^ PROBLEM NO\.:\s+\d/ );
 830           $start_pos-- and last if (/^[a-df-zA-DF-Z]/); #Rewind one step
 831         }
 832         if ( scalar @eigens > 0 ) {
 833           my @list = sort { $a <=> $b } @eigens;
 834           $self -> {'condition_number'} = abs($list[$#list]/$list[0]) if ( $list[0] != 0 );
 835         }
 836         @{$self -> {'eigens'}} = @eigens;
 837       }
 838 end _read_eigen
 839
 840 # }}} _read_eigen
 841
 842 # {{{ _read_eval
 843 start _read_eval
 844 {
 845     my $start_pos = $self -> {'lstfile_pos'};
 846
 847     while( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) {
 848       if ( / PROBLEM NO\.:\s*\d+\n/ ) {
 849         $self -> parsing_error( message => "Error in reading the number of ".
 850                                 "function evaluations!\nNext problem found" );
 851         return;
 852       }
 853
 854       if ( $start_pos >= scalar @{$self -> {'lstfile'}} ) {
 855         $self -> parsing_error( message => "Error in reading number of ".
 856                                 "function evaluations!\nEOF found\n" );
 857         return;
 858       }
 859
 860       if ( /^1/ or /MINIMUM VALUE OF OBJECTIVE FUNCTION/ ) {
 861         $self -> parsing_error( message => "Error in reading number of ".
 862                                 "function evaluations!\nOFV found" );
 863         return;
 864       }
 865
 866       if ( / NO. OF FUNCTION EVALUATIONS USED:\s*(\d*)/ ){
 867         $self -> feval($1);
 868         last;
 869       }
 870     }
 871 }
 872 end _read_eval
 873 # }}} _read_eval
 874
 875 # {{{ _read_iteration_path
 876 start _read_iteration_path
 877       {
 878         my $start_pos = $self -> {'lstfile_pos'};
 879         my $success = 0;
 880         my (@func_eval, @parameter_vectors,
 881             @gradient_vectors) = ((), (), (), (), ());
 882         my @numsigdig_vectors;
 883         my $cumulative_evals = 0;
 884         my $zero_gradients = 0;
 885         my $hessian_reset = 0;
 886         my $found_monitoring = 0;
 887         while ( $_ = @{$self -> {'lstfile'}}[ $start_pos++ ] ) { #Large reading loop
 888           if ( /MONITORING OF SEARCH/ ) {
 889             $found_monitoring = 1;
 890           }
 891
 892           if ( $found_monitoring and
 893                $start_pos >= scalar @{$self -> {'lstfile'}} ) {
 894             # This is probably (we cannot be really sure but it is
 895             # likely) a "NONMEM black hole"-situation. We stop parsing
 896             # but do not indicate that we had problems parsing the
 897             # file.
 898             $self -> parsing_error_message("Found \" MONITORING OF SEARCH:\" but no".
 899                                            " records of the iterations before the end".
 900                                            " of the output file. This is a sign of a".
 901                                            " \"NONMEM black hole\"-situation. We cannot ".
 902                                            "be 100% sure of this based on this type of".
 903                                            " output so please take a good look at the files\n");
 904             $self -> parsed_successfully(1);
 905             $self -> finished_parsing(1);
 906             return;
 907           }
 908
 909           if ( /0PROGRAM TERMINATED BY OBJ/ ) {
 910             # This is an error message which terminates NONMEM. We
 911             # return after reading the minimization message
 912             $self -> minimization_successful(0);
 913             $self -> _read_minimization_message();
 914             $self -> finished_parsing(1);
 915             return;
 916           }
 917
 918           if (/^0ITERATION NO/) {
 919             unless (/0ITERATION NO\.:\s+(\d+)\s+OBJECTIVE VALUE:\s+(\S+)\s+NO\. OF FUNC\. EVALS\.:\s*(.+)/) {
 920               debug -> warn( level => 1,
 921                              message => "Error in reading iteration path!\n$!" );
 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         for (@mess) {
1606           if (s/\s+ETABAR:\s+//) {
1607             @temp = split;
1608             push @etabar, @temp;
1609 #           print "Etabar Push of: $_ \n";
1610             for (@mess) {       #If wrapped etabar-lines
1611               last if(/^1/);
1612               last if(/a-zA-DF-Z/);
1613               last unless (s/^\s+//);
1614               last unless (/\d/);
1615               @temp = split;
1616               push @etabar, @temp;
1617             }
1618             last;
1619           }
1620         }
1621         # Initialize the attribute only if we have found any data
1622         $self -> {'etabar'} = \@etabar if ( $#etabar > 0 );
1623
1624
1625         for (@mess) {
1626           if (s/\s+P VAL\.:\s+//) {
1627             @temp = split;
1628             push @pval, @temp;
1629 #           print "Pval Push of: $_ \n";
1630             for (@mess) {       #If wrapped etabar-lines
1631               last if(/^1/);
1632               last if(/a-zA-DF-Z/);
1633               last unless (s/^\s+//);
1634               last unless (/\d/);
1635               @temp = split;
1636               push @pval, @temp;
1637             }
1638             last;
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