lib/fcon_subs.pm

   1 # 2004-12-08: FCON must be adjusted for multiple problems /Lasse
   2
   3 start include statements
   4 # A Perl module for parsing and manipulating NONMEM FCON files
   5 #use Math::Random;
   6 #use ui;
   7 end include statements
   8
   9 # {{{ _get_array
  10 start _get_array
  11 {
  12     my $items_collected = 0;
  13     while( $items_collected < $nr_items ){
  14         my $items_per_line = ($nr_items - $items_collected) < (72 / $tabSize) ? ($nr_items - $items_collected): (72 / $tabSize);
  15
  16         my $line = $fcon[$index++];
  17
  18         # strip head
  19         $line =~ s/^.{8}//;
  20
  21         for( ; $items_per_line > 0; $items_per_line -- ){
  22             $line =~ s/^(.{$tabSize})//;
  23             my $val = $1;
  24             unless( $val =~ /[\ ]{8}/ ){ # Check for blanks
  25                 $val = $val * 1; # Cheeky trick to make $val a number instead of a string
  26                 push(@values, $val);
  27             }
  28             $items_collected ++;
  29         }
  30     }
  31     push(@values, $index-1);
  32 }
  33 end _get_array
  34 # }}}
  35
  36 # {{{ _skip_lines
  37 start _skip_lines
  38 {
  39     while( not $fcon[$index+1] =~ /^[A-Z]{4}[\ ]{3}/ ){
  40         $index++;
  41         if( $index > $#fcon ){
  42             last;
  43         }
  44     }
  45 }
  46 end _skip_lines
  47 # }}}
  48
  49 # {{{ run
  50
  51 start run
  52 {
  53     for( my $i = $batchStart; $i <  $batchStart+$batchSize; $i++ ){
  54         for( my $try = 0; $try < $retries[$i - $batchStart]; $try++ ){
  55             system( "./NM_run".($i+1) . '_nonmem.sh' ) == 0 or die "Unable to execute NM_run".($i+1) . "_nonmem.sh $!\n";
  56             open( my $outfile, "<NM_run".($i+1) . '_psn.lst' );
  57             my $success = 1;
  58             while( <$outfile> ){
  59                 if((/^0MINIMIZATION TERMINATED/ or /^0PROGRAM TERMINATED/) # If hard termination
  60                    or # We are picky and check some softer rerun requirements
  61                    ($picky[$i - $batchStart] and
  62                     (/0ESTIMATE OF THETA IS NEAR THE BOUNDARY AND/ or
  63                      /0PARAMETER ESTIMATE IS NEAR ITS BOUNDARY/    or
  64                      /0R MATRIX ALGORITHMICALLY SINGULAR/          or
  65                      /0S MATRIX ALGORITHMICALLY SINGULAR/
  66                      )
  67                     )
  68                    ){
  69                     print "NONMEM failed. Preparing rerun number" .  $try+1 . " / " . $retries[$i - $batchStart] . "\n";
  70                     $success = 0;
  71                     $self -> whipe();
  72                     $self -> {'filename'} = "NM_run".($i+1) . '_FCON';
  73                     $self -> parse('filename' => "NM_run".($i+1) . '_FCON' );
  74                     $self -> pertubate_all('fixed_thetas' => $fixed_thetas[$i],
  75                                            'fixed_omegas' => $fixed_omegas[$i],
  76                                            'fixed_sigmas' => $fixed_omegas[$i]);
  77                     $self -> write('filename' => "NM_run".($i+1) . '_FCON' );
  78                 }
  79             }
  80             close( $outfile );
  81             if($success){ last };
  82         }
  83     }
  84 }
  85
  86
  87 end run
  88
  89 # }}}
  90
  91 # {{{ whipe
  92
  93 start whipe
  94 {
  95     $self -> {'fcon_lines'} = undef;
  96     $self -> {'theta_values'} = undef;
  97     $self -> {'have_omega_DIAG'} = undef;
  98     $self -> {'omega_values'} = undef;
  99     $self -> {'have_sigma_DIAG'} = undef;
 100     $self -> {'have_omega_BLST'} = undef;
 101     $self -> {'have_sigma_BLST'} = undef;
 102 }
 103 end whipe
 104
 105 # }}}
 106
 107 # {{{ parse
 108 start parse
 109 {
 110     my $fcon;
 111     open($fcon,"<",$self -> {'filename'}) or die "fcon_file: Unable to open ", $self -> {'filename'}, " : $!";
 112     my @fcon_lines =  <$fcon>;
 113     $self -> fcon_lines( \@fcon_lines );
 114
 115     close($fcon);
 116
 117     my $prob = -1;
 118
 119     my @have_FIND;
 120     my @have_initial_STRC;
 121
 122     my @have_sigma_STRC;
 123     my @have_omega_STRC;
 124
 125     my @expect_sigma_BLST;
 126     my @expect_omega_BLST;
 127
 128     my @have_sigma_BLST;
 129     my @have_omega_BLST;
 130     my @sigma_BLST_lengths;
 131     my @omega_BLST_lengths;
 132
 133     my @expect_sigma_DIAG;
 134     my @expect_omega_DIAG;
 135
 136     my @have_sigma_DIAG;
 137     my @have_omega_DIAG;
 138
 139     my @nr_thetas;
 140     my @expect_bounds;
 141
 142     for( my $i = 0; $i <= $#fcon_lines; $i++ ){
 143         my $line = $fcon_lines[$i];
 144         if( $line =~ /^FILE/ ){
 145             #print "$i " . $line;
 146         }elsif( $line =~ /^PROB/ ){
 147             $prob++;
 148             $self -> prob($prob);
 149             #print "$i " . $line;
 150         }elsif ( $line =~ /^DATA/ ){
 151             #print "$i " .$line;
 152         }elsif ( $line =~ /^ITEM/ ){
 153             #print "$i " .$line;
 154         }elsif ( $line =~ /^INDX/ ){
 155             #print "$i " .$line;
 156             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 157                                        'index' => $i );
 158         }elsif ( $line =~ /^LABL/ ){
 159             #print "$i " .$line;
 160             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 161                                        'index' => $i );
 162         }elsif ( $line =~ /^FORM/ ){
 163             #print "$i " .$line;
 164             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 165                                        'index' => $i );
 166         }elsif ( $line =~ /^FIND/ ){
 167             #print "$i " .$line;
 168             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 169                                        'index' => $i );
 170             $have_FIND[$prob] = 1;
 171         }elsif ( $line =~ /^STRC/ ){
 172             #print "$i " .$line;
 173             if( $have_FIND[$prob] ){
 174                 die "Malformated FCON (STRC unexpectedly found)\n";
 175             }
 176
 177             my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 178                                                 'index' => $i,
 179                                                 'tabSize' => 4,
 180                                                 'nr_items' => 9 )};
 181             $i = pop(@values);
 182
 183             unless( $have_initial_STRC[$prob] ){
 184                 $have_initial_STRC[$prob] = 1;
 185                 $nr_thetas[$prob] = $values[0];
 186                 #$expect_bounds = $values[3];
 187                 $expect_omega_DIAG[$prob] = $values[1];
 188                 $expect_sigma_DIAG[$prob] = $values[2];
 189                 $expect_omega_BLST[$prob] = $values[5] ? 0 : $values[6];
 190                 $expect_sigma_BLST[$prob] = $values[7] ? 0 : $values[8];
 191             } else {
 192                 unless( $have_omega_STRC[$prob] ){
 193                     if( $expect_omega_BLST[$prob] ){
 194                         @omega_BLST_lengths[$prob] = @values;
 195                     }
 196                     $have_omega_STRC[$prob] = 1;
 197                 } else {
 198                     if( $expect_sigma_BLST[$prob] ){
 199                         @sigma_BLST_lengths[$prob] = @values;
 200                     }
 201                     $have_sigma_STRC[$prob] = 1;
 202                 }
 203             }
 204
 205         }elsif ( $line =~ /^THTA/ ){
 206             #print "$i " .$line;
 207             if( $have_initial_STRC[$prob] ){
 208                 my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 209                                                     'index' => $i,
 210                                                     'tabSize' => 8,
 211                                                     'nr_items' => $nr_thetas[$prob] )};
 212                 $i = pop(@values);
 213                 $self -> {'theta_values'} = [] unless defined $self -> {'theta_values'};
 214                 $self -> theta_values -> [$prob] = \@values;
 215             } else {
 216                 die "Malformated FCON (THTA unexpectedly found)\n";
 217             }
 218
 219         }elsif ( $line =~ /^LOWR/ ){
 220             #print "$i " .$line;
 221             #if( $expect_bounds ){
 222                 my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 223                                                     'index' => $i,
 224                                                     'tabSize' => 8,
 225                                                     'nr_items' => $nr_thetas[$prob] )};
 226                 $i = pop(@values);
 227                 $self -> {'lower_bounds'} = [] unless defined $self -> {'lower_bounds'};
 228                 $self -> lower_bounds -> [$prob] = \@values;
 229             #} else { ## I'd like to uncomment this.. but most FCONs seams malformated in this way. And its easy to ignore
 230                 #die "Malformated FCON (LOWR unexpectedly found)\n";
 231             #}
 232         }elsif ( $line =~ /^UPPR/ ){
 233             #print "$i " .$line;
 234             #if( $expect_bounds ){
 235                 my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 236                                                     'index' => $i,
 237                                                     'tabSize' => 8,
 238                                                     'nr_items' => $nr_thetas[$prob] )};
 239                 $i = pop(@values);
 240                 $self -> {'upper_bounds'} = [] unless defined $self -> {'upper_bounds'};
 241                 $self -> upper_bounds -> [$prob] =  \@values;
 242             #} else { ## I'd like to uncomment this.. but most FCONs seams malformated in this way. And its easy to ignore
 243                 #die "Malformated FCON (UPPR unexpectedly found)\n";
 244             #}
 245         }elsif ( $line =~ /^DIAG/ ){
 246           #print "$i " .$line;
 247           if( $expect_omega_DIAG[$prob] and not $have_omega_DIAG[$prob] and not ($expect_omega_BLST[$prob] or $have_omega_BLST[$prob]) ){
 248             my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 249                                                 'index' => $i,
 250                                                 'tabSize' => 8,
 251                                                 'nr_items' => $expect_omega_DIAG[$prob] )};
 252             $i = pop(@values);
 253             $self -> {'omega_values'} = [] unless defined $self -> {'omega_values'};
 254             $self -> omega_values -> [$prob] -> [0] = \@values;
 255             $have_omega_DIAG[$prob] = 1;
 256             $self -> {'have_omega_DIAG'} = [] unless defined $self -> {'have_omega_DIAG'};
 257             $self -> have_omega_DIAG -> [$prob] = 1;
 258           } elsif( $expect_sigma_DIAG[$prob] and not $have_sigma_DIAG[$prob] and not ($expect_sigma_BLST[$prob] or $have_sigma_BLST[$prob]) ) {
 259             my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 260                                                 'index' => $i,
 261                                                 'tabSize' => 8,
 262                                                 'nr_items' => $expect_sigma_DIAG[$prob] )};
 263             $i = pop(@values);
 264             $self -> {'sigma_values'} = [] unless defined $self -> {'sigma_values'};
 265             $self -> sigma_values -> [$prob] -> [0] = \@values;
 266             $have_sigma_DIAG[$prob] = 1;
 267             $self -> {'have_sigma_DIAG'} = [] unless defined $self -> {'have_sigma_DIAG'};
 268             $self -> have_sigma_DIAG -> [$prob] = 1;
 269           }
 270         }elsif ( $line =~ /^BLST/ ){
 271             #print "$i " .$line;
 272             if( $expect_sigma_BLST[$prob] or $expect_omega_BLST[$prob] ){
 273                 if( $expect_omega_BLST[$prob] ){
 274                     $expect_omega_BLST[$prob] --;
 275
 276                     my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 277                                                         'index' => $i,
 278                                                         'tabSize' => 8,
 279                                                         'nr_items' => 9 )};
 280                     $i = pop(@values);
 281                     $self -> omega_values -> [$prob] -> [$have_omega_BLST[$prob]] = \@values;
 282
 283                     $have_omega_BLST[$prob]++;
 284                     $self -> {'have_omega_BLST'} = [] unless defined $self -> {'have_omega_BLST'};
 285                     $self -> have_omega_BLST -> [$prob] = 1;
 286                 } elsif( $expect_sigma_BLST[$prob] ){
 287                     $expect_sigma_BLST[$prob] --;
 288
 289                     my @values = @{$self -> _get_array( 'fcon' => \@fcon_lines,
 290                                                         'index' => $i,
 291                                                         'tabSize' => 8,
 292                                                         'nr_items' => 9 )};
 293                     $i = pop(@values);
 294
 295                     $self -> sigma_values -> [$prob] -> [$have_sigma_BLST[$prob]] = \@values;
 296
 297                     $have_sigma_BLST[$prob]++;
 298                     $self -> {'have_sigma_BLST'} = [] unless defined $self -> {'have_sigma_BLST'};
 299                     $self -> have_sigma_BLST -> [$prob] = 1;
 300                 }
 301             } else {
 302                 die "Malformated FCON (DIAG unexpectedly found)\n";
 303             }
 304         }elsif ( $line =~ /^ESTM/ ){
 305             #print "$i " .$line;
 306             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 307                                        'index' => $i );
 308         }elsif ( $line =~ /^COVR/ ){
 309             #print "$i " .$line;
 310             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 311                                        'index' => $i );
 312         }elsif ( $line =~ /^TABL/ ){
 313             #print "$i " .$line;
 314             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 315                                        'index' => $i );
 316         }elsif ( $line =~ /^SCAT/ ){
 317             #print "$i " .$line;
 318             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 319                                        'index' => $i );
 320         }else {
 321              print "Found something new: \n $line \n" if $self -> debug;
 322             $i = $self -> _skip_lines( 'fcon' => \@fcon_lines,
 323                                        'index' => $i );
 324         }
 325     }
 326 }
 327 end parse
 328 # }}}
 329
 330 # {{{ write
 331
 332 start write
 333 {
 334     my $skip_line = 0;
 335     my @have_omega_DIAG;
 336     my @omega_BLST_count;
 337     my @sigma_BLST_count;
 338
 339     my $prob = -1;
 340
 341     open( my $outfile, ">". $filename ) or die "fcon_file: Unable to open FCON file: $!";
 342
 343     foreach my $line (@{$self -> fcon_lines}){
 344         if( $line =~ /^PROB/ ){
 345             $prob ++;
 346             print $outfile $line;
 347         } elsif ( $line =~ /^(THTA.{4})/ ){
 348             print $outfile $1;
 349             my $i = 1;
 350             $skip_line = 1;
 351             foreach my $theta ( @{$self -> theta_values -> [$prob]} ){
 352                 print $outfile "\n        " unless( $i % 10 );
 353                 $i++;
 354                 printf $outfile ("%8s",$theta);
 355             }
 356             print $outfile "\n";
 357         } elsif ( $line =~ /^(DIAG.{4})/ ){
 358             print $outfile $1;
 359             my $i = 1;
 360             my @diag;
 361             $skip_line = 1;
 362             if( $self -> have_omega_DIAG -> [$prob] && !$have_omega_DIAG[$prob] ){
 363                 @diag = @{$self -> omega_values -> [$prob] -> [0]};
 364                 $have_omega_DIAG[$prob] = 1;
 365             } elsif ( $self -> have_sigma_DIAG -> [$prob]) {
 366                 @diag = @{$self -> sigma_values -> [$prob] -> [0]};
 367             }
 368             foreach my $diag ( @diag ){
 369                 print $outfile "\n        " unless( $i % 10 );
 370                 $i++;
 371                 printf $outfile ("%8s",$diag);
 372             }
 373             print $outfile "\n";
 374         } elsif ( $line =~ /^(BLST.{4})/ ) {
 375             print $outfile $1;
 376             my $i = 1;
 377             my @blst;
 378             $skip_line = 1;
 379             if( $self -> have_omega_BLST -> [$prob] && $omega_BLST_count[$prob] <= $#{$self -> omega_values -> [$prob]} ){
 380                 @blst = @{$self -> omega_values -> [$prob] -> [$omega_BLST_count[$prob]++]};
 381             } else {
 382                 @blst = @{$self -> sigma_values -> [$prob] -> [0]};
 383             }
 384
 385             foreach my $blst ( @blst ){
 386                 print $outfile "\n        " unless( $i % 10 );
 387                 $i++;
 388                 printf $outfile ("%8s",$blst);
 389             }
 390             print $outfile "\n";
 391         } elsif( not $line =~ /^[A-Z]{4}[\ ]{3}/ and $skip_line){
 392             next;
 393         } else {
 394             $skip_line = 0;
 395             print $outfile $line;
 396         }
 397     }
 398 }
 399 end write
 400
 401 # }}}
 402
 403 # {{{ pertubate_all
 404 start pertubate_all
 405 {
 406     for( my $i = 0; $i <= $self -> prob; $i++ ){
 407         $self -> _pertubate( 'estimates' => $self -> theta_values -> [$i],
 408                              'degree' => $degree,
 409                              'lobnds' => $self -> lower_bounds -> [$i],
 410                              'upbnds' => $self -> upper_bounds -> [$i],
 411                              'fixed' => $fixed_thetas[$i]);
 412
 413         foreach my $sigma_line ( @{$self -> sigma_values -> [$i]} ){
 414             $self -> _pertubate( 'estimates' => $sigma_line,
 415                                  'degree' => $degree,
 416                                  'fixed' => $fixed_sigmas[$i]);
 417         }
 418
 419         foreach my $omega_line ( @{$self -> omega_values -> [$i]} ){
 420             $self -> _pertubate( 'estimates' => $omega_line,
 421                                  'degree' => $degree,
 422                                  'fixed' => $fixed_omegas[$i]);
 423         }
 424     }
 425 }
 426 end pertubate_all
 427 # }}}
 428
 429 # {{{ pertubate
 430
 431 start _pertubate
 432 {
 433   for( my $i = 0; $i <= $#{$estimates}; $i++ ){
 434     unless( $fixed[$i] ){
 435       my ( $sign, $est, $form );
 436       my $lobnd;
 437       my $upbnd;
 438       my $init = $estimates -> [$i];
 439
 440       my $change = abs($degree*$init);
 441
 442       if ( defined $lobnds and defined $lobnds -> [$i] ){
 443         $lobnd = $lobnds -> [$i];
 444         $lobnd = $init-$change < $lobnd ? $lobnd : $init-$change;
 445       } else {
 446         $lobnd = $init-$change;
 447       }
 448
 449       if ( defined $upbnds and defined $upbnds -> [$i] ) {
 450         $upbnd = $upbnds -> [$i];
 451         $upbnd = $init+$change > $upbnd ? $upbnd : $init+$change;
 452       } else {
 453         $upbnd = $init+$change;
 454       }
 455
 456       $lobnd = 0.01 if ( ( $lobnd < 0.01 and $lobnd > -0.01)
 457                          and $upbnd >= 0.01001 );
 458       $upbnd = -0.01 if ( ( $upbnd < 0.01 and $upbnd > -0.01)
 459                           and $lobnd <= -0.01001 );
 460
 461
 462       if ( $lobnd <= -0.01 and $upbnd >= 0.01 ) {
 463         $est = $lobnd + 0.02 + rand($upbnd - ($lobnd+0.02));
 464         #$est = random_uniform(1, $lobnd + 0.02, $upbnd);
 465         $est = $est - 0.02 if ( $est <0.01 );
 466       } else {
 467         $est = $lobnd + rand($upbnd - $lobnd);
 468         #$est = random_uniform(1, $lobnd, $upbnd );
 469       }
 470       $form  = "%6.4f" if $est < 1000 and $est > -999;
 471       $form  = "%6.1f" if $est >= 1000 or $est <=-999;
 472
 473
 474       $estimates -> [$i] = sprintf $form,$est;
 475       if ( $estimates -> [$i] == 0 ) {
 476         $estimates -> [$i] = '0.0001';
 477       }
 478     }
 479   }
 480 }
 481 end _pertubate
 482
 483 # }}}