lib/tool/modelfit_subs.pm

   1 # {{{ include
   2
   3 start include statements
   4       {
   5         # Perl libraries #
   6         use Config;
   7         use Cwd;
   8         use Data::Dumper;
   9         use File::Copy qw/cp mv/;
  10         use File::Path;
  11         use FindBin qw($Bin);
  12         # External #
  13         use Storable;
  14         use Math::Random;
  15         # PsN includes #
  16         use fcon;
  17         use grid::nordugrid::xrsl_file;
  18         if( $PsN::config -> {'_'} -> {'use_keyboard'} ) {
  19           require hotkey;
  20         }
  21         use nonmem;
  22         use POSIX ":sys_wait_h";
  23         use output;
  24         use OSspecific;
  25         use ui;
  26         use moshog_client;
  27         use Time::HiRes;
  28         use ext::IPC::Run3;
  29       }
  30 end include statements
  31
  32 # }}} include statements
  33
  34 # {{{ description, examples, synopsis, see also
  35
  36 # No method, just documentation
  37 start description
  38     #
  39     # In PsN versions < 2.0, the functionality for actually running
  40     # NONMEM on models and data PsN objects was provided by the model
  41     # class. As of PsN versions 2.0 and higher, this functinality has
  42     # been moved to the separate class I<modelfit> in order to make the
  43     # class responsibilities clearer.
  44     #
  45     # Fitting a model can be viewed as a special case of a more
  46     # general set of tools for population PK/PD. In PsN, the
  47     # modelfit class is therefore a specialization of a general PsN
  48     # tool class. The tool class itself is not capable of much at
  49     # all but to define a common structure for all PsN tools.
  50     #
  51     # All attributes and (class) methods specified for the general
  52     # tool class are inherited in the modelfit class. Some (class) methods
  53     # are defined in both classes (e.g. the L</run>) and in these
  54     # cases it is the modelfit version that will be used.
  55     #
  56     # =begin html
  57     #
  58     # <tr>Please look at the documentation for the <a
  59     # href="../tool.html">general tool class</a> as well as the <a
  60     # href="#examples">examples</a> section of this document for
  61     # descriptions of the setting that are available for all
  62     # tools.</tr>
  63     #
  64     # =end html
  65     #
  66     # =begin man
  67     #
  68     # Please look at the documentation for the general I<tool> class
  69     # as well as the L</examples> section of this document for
  70     # descriptions of the setting that are available for all tools.
  71     #
  72     # =end man
  73     #
  74 end description
  75
  76 start examples
  77     # The following code may be used to create a simple modelfit
  78     # object and to run it.
  79     #
  80     #   use lib qw(path to PsN installation/lib);
  81     #   use tool::modelfit;
  82     #   use model;
  83     #
  84     #   my $model = model -> new( filename => 'run1.mod' );
  85     #   my $modelfit = tool::modelfit -> new( models => [$model] );
  86     #   my %results = %{$modelfit -> run};
  87     #
  88     # To illustrate a more complex use of modelfit, we can run a
  89     # bootstrap of 200 resampled data sets of same size as the
  90     # original data set. In this example we assume that the modelfile
  91     # we use contains only one problem and one sub problem.
  92     #
  93     #   use lib qw(path to PsN installation/lib);
  94     #   use tool::modelfit;
  95     #   use data;
  96     #   use model;
  97     #   use output;
  98     #   Math::Random;       # For perturbation of initial estimates
  99     #                       # after unsuccessful minimizations
 100     #
 101     #   # set these to appropriate values for your own run
 102     #   my $samples = 200;
 103     #   my $modelfile = 'run1.mod';
 104     #   my $boot_sample_file = 'boot_ind.csv';
 105     #   my $boot_results_file = 'boot_results.csv';
 106     #
 107     #   # set the seed from a phrase (consistent over different hardware,
 108     #   # see Math::Random )
 109     #   random_set_seed_from_phrase('testing');
 110     #
 111     #   # ignore missing data and (especially) output files
 112     #   my $model = model -> new( filename => $modelfile,
 113     #                             ignore_missing_files => 1 );
 114     #
 115     #   my $data = $model -> datas -> [0];
 116     #
 117     #   # Create the bootstrap data sets. The default name for each
 118     #   # new resampled data file will be bs###.dta, where ### is a
 119     #   # number from 1 to $samples.
 120     #   my ( $dataref, $incl_ind_ref, $incl_keys_ref ) =
 121     #       $data -> bootstrap ( samples => $samples );
 122     #
 123     #   # Save the resampled ID numbers in a file, one row for each
 124     #   # bootstrap data set
 125     #   open ( LOG, ">$boot_sample_file" );
 126     #   foreach my $sample ( @{$incl_ind_ref} ) {
 127     #     print LOG join(';', @{$sample} ),"\n";
 128     #   }
 129     #
 130     #   # Create the boostrap models and set the correct resampled
 131     #   # data file in $DATA. Save them in @bs_models.
 132     #   my @bs_models  = ();
 133     #   for ( my $i = 1; $i <= $samples; $i++ ) {
 134     #     my $bs_mod = $model -> copy( filename    => 'bs'.$i.'.mod',
 135     #                                  copy_data   => 0,
 136     #                                  copy_output => 0);
 137     #     $bs_mod -> datafiles( new_names     => ['bs'.$i.'.dta'],
 138     #                           absolute_path => 1 );
 139     #     $bs_mod -> _write;
 140     #     push( @bs_models, $bs_mod );
 141     #   }
 142     #
 143     #   # Create a modelfit object with the bootstrap models as
 144     #   # input. Set the number of parallel runs to 2.
 145     #   my $mfitobj = tool::modelfit -> new ( models     => \@bs_models,
 146     #                                           threads    => 2 );
 147     #
 148     #   # Run the model fit. Since the bootstrap models are named
 149     #   # bs###.mod, the default names for the output files will be
 150     #   # bs###.lst.
 151     #   $mfitobj -> run;
 152     #
 153     #   # We'll save the OFV plus the theta, omega and sigma estimates
 154     #   # for each run in a file.
 155     #   open( RESULT, ">$boot_results_file" );
 156     #
 157     #   for ( my $i = 1; $i <= $samples; $i++ ) {
 158     #     my $bs_out = output -> new( filename => 'bs'.$i.'.lst' );
 159     #     my @ofv    = @{$bs_out -> ofv};
 160     #     my @thetas = @{$bs_out -> thetas};
 161     #     my @omegas = @{$bs_out -> omegas};
 162     #     my @sigmas = @{$bs_out -> sigmas};
 163     #     # We know that we only have one problem and one sub problem in this
 164     #     # example. ([0][0] comes from that fact)
 165     #     my @print_strings = ();
 166     #     push( @print_strings, $ofv[0][0] );
 167     #     push( @print_strings, @{$thetas[0][0]} );
 168     #     push( @print_strings, @{$omegas[0][0]} );
 169     #     push( @print_strings, @{$sigmas[0][0]} );
 170     #     print RESULT join( ';', @print_strings ),"\n";
 171     #   }
 172     #
 173     #   close( RESULT );
 174     #
 175     #   # We're done!
 176 end examples
 177
 178 start synopsis
 179     #   use tool::modelfit;
 180     #   use model;
 181     #
 182     #   my $model_obj = model -> new ( filename => 'run1.mod' );
 183     #
 184     #   my $modelfit_obj = tool::modelfit -> new ( models => [$model_obj] );
 185     #
 186     #   my $output_obj = $modelfit_obj -> run;
 187 end synopsis
 188
 189 start see_also
 190     # =begin html
 191     #
 192     # <a HREF="../data.html">data</a>, <a
 193     # HREF="../model.html">model</a> <a
 194     # HREF="../output.html">output</a>, <a
 195     # HREF="../tool.html">tool</a>
 196     #
 197     # =end html
 198     #
 199     # =begin man
 200     #
 201     # data, model, output, tool
 202     #
 203     # =end man
 204 end see_also
 205
 206 # }}}
 207
 208 # {{{ new
 209
 210 start new
 211     # Usage:
 212     #
 213     #   $modelfit_object = tool::modelfit -> new( models  => [$model_object],
 214     #                                               retries => 5 );
 215     #
 216     # This is the basic usage and it creates a modelfit object that
 217     # can later be run using the L</run> method. I<models> is an array
 218     # of PsN model objects.
 219     #
 220     #   $modelfitObject = $tool::modelfit -> new( 'retries' => 5 );
 221     #   $modelfitObject -> add_model( init_data => { filename => $modelFileName } );
 222     #
 223     # This way of using modelfit is suitable if you have a list with
 224     # filenames of modelfiles. "add_model> will create modelfitobject
 225     # for you.
 226     #
 227     # A more interresting attribute is I<threads> which sets how many
 228     # parallel executions of NONMEM that will run. Some tips are:
 229     # Setting the number of threads higher than the number of nodes in
 230     # your cluster/supercomputer can make your runs slower. The
 231     # biggest limiting factor is the amount of memory needed by
 232     # NONMEM. With smaller runs, just set the thread number to the
 233     # number of nodes available.
 234     #
 235     # The I<directory> is the folder where the tools stores
 236     # temporary data and runs subtools (or in the modelfit case,
 237     # runs NONMEM). Each NONMEM run will have its own sub directory
 238     # NM_run[X] where [X] is an index running from 1 to the number of
 239     # runs. If unsure of what this means, leave it undefined and a
 240     # default will be used, e.g. modelfit_dir3 or something.
 241     #
 242     # Next, the I<compress> and I<remove_temp_files> attributes are good
 243     # if you want to save some hard disk space. I<compress> set to 1
 244     # will put all NONMEM output in to an tar/gz archive named
 245     # I<nonmem_files.tgz> placed in the I<NM_run[X]> directory
 246     # described above. If I<remove_temp_files> is set to 1,  the NONMEM
 247     # files: 'FCON', 'FDATA', 'FSTREAM', 'PRDERR' will be removed.
 248     #
 249     # I<clean> is a stronger version of I<remove_temp_files>; it will also
 250     # remove I<NM_run[X]> and all that is in these.
 251     #
 252     # I<retries> is the number of times L</run> will alter initial
 253     # values and (re)execute NONMEM when executions fail. I<retries>
 254     # can either be an integer, specifying the number of retries for
 255     # all models, or it can be an array with the number of retries
 256     # specific for each modelfile as elements. The default value is
 257     # B<5>. The algorithm for altering the initial values works
 258     # roughly like this: For each each new try, a random new initial
 259     # value is drawn from a uniform distribution with limits +-n*10%
 260     # of the original intial estimate and where n i equal to the retry
 261     # number. I.e. the first retry, the borders of the distribution
 262     # are +-10%. The algorithm ensures that the new values are within
 263     # specified boundaries.
 264     #
 265     # =begin html
 266     #
 267     # For a full dexcription of the algorithm, see <a
 268     # href="../model/problem/record/init_option.html#set_random_init">set_random_init</a>
 269     # of the <a
 270     # href="../model/problem/record/init_option.html">init_option
 271     # class</a>.
 272     #
 273     # =end html
 274     #
 275     # =begin man
 276     #
 277     # For a full dexcription of the algorithm, see I<set_random_init>
 278     # of the I<init_option> class.
 279     #
 280     # =end man
 281     #
 282     # If I<picky> is set to 1, the output from NONMEM will be checked
 283     # more thoroughly. If any of the lines below are found in the
 284     # minimization message, a rerun is initiated.
 285     #
 286     #    COVARIANCE STEP ABORTED
 287     #    PROGRAM TERMINATED BY OBJ
 288     #    ESTIMATE OF THETA IS NEAR THE BOUNDARY AND
 289     #    PARAMETER ESTIMATE IS NEAR ITS BOUNDARY
 290     #    R MATRIX ALGORITHMICALLY SINGULAR
 291     #    S MATRIX ALGORITHMICALLY SINGULAR
 292     #
 293     # I<nm_version> is a string with the version number of NONMEM that
 294     # will be used. The installed versions of NONMEM must be specified
 295     # in OSspecific.pm, the class responsible for system specific
 296     # features settings.
 297     #
 298     # I<logfile> specifies the name of the logfile.
 299     #
 300     # if I<silent_logfile> is defined all NONMEM output will
 301     # be written to I<NM_run[X]/xxx>, where xxx is the defined value.
 302     #
 303     # I<extra_files> is an array of strings where each string is a
 304     # file needed for NONMEM execution. Those file will be moved
 305     # to the I<NM_run[X]> directory.
 306     #
 307     # I<seed> is just a way to set a random seed number.
 308     #
 309     # If I<run_on_nordugrid> is set to true modelfit will submit the nonmem
 310     # runs to a grid. A group of related parameters are also
 311     # specified.
 312
 313     # I<cpuTime> Is an estimated execution time for each individual
 314     # modelfile. It should preferably be a bit longer than reality. If
 315     # you specify a cpuTime that is to short, you risk that the grid
 316     # kills your jobs prematurely. The unit of I<cpuTime> is minutes.
 317
 318     # I<grid_cpuTime> is the time of the actual grid job. It should be
 319     # used to group modelfiles together. For example, if you set
 320     # I<cpuTime> to ten minutes, I<grid_cpuTime> to 60 minutes and the
 321     # number of modelfiles is 14 modelfit will create three grid jobs,
 322     # two with six model files each and one with two modelfiles.
 323
 324     # I<grid_adress> is the URL of the grid submission server,
 325     # e.g. hagrid.it.uu.se.
 326
 327
 328       {
 329         if ( defined $this -> {'logfile'} ) {
 330           my $dir;
 331           $this -> {'logfile'} = join('', OSspecific::absolute_path( $this -> {'directory'},
 332                                                                      $this -> {'logfile'}) );
 333         }
 334         if ( $this -> {'ask_if_fail'} ) {
 335           eval( 'use Tk' );
 336         }
 337
 338         $this -> {'run_local'} = 1;
 339
 340         if( $this -> {'run_on_lsf'} or
 341             $this -> {'run_on_ud'} or
 342             $this -> {'run_on_umbrella'} or
 343             $this -> {'run_on_sge'} ){
 344           $this -> {'run_local'} = 0;
 345         }
 346
 347         if( $this -> {'handle_msfo'} ){
 348           $this -> {'handle_crashes'} = 1;
 349         }
 350
 351         if( $this -> {'handle_crashes'} ){
 352           if( $this -> {'crash_restarts'} > 0 ){
 353             $this -> {'crash_restarts'}++;
 354           }
 355         }
 356
 357         $this -> calculate_raw_results_width();
 358
 359         $this -> {'raw_line_structure'} = ext::Config::Tiny -> new();
 360
 361       }
 362
 363 end new
 364
 365 # }}} new
 366
 367 # {{{ calculate_raw_results_width
 368
 369 start calculate_raw_results_width
 370   {
 371
 372     # 2008-01-24
 373     # This code comes largely from "prepare_raw_results" which should
 374     # be split over several functions to fit the serialized version of
 375     # PsN.
 376
 377     # Some column in "raw_results_header" are meta-columns, they
 378     # will be replaced by several columns. For example, the
 379     # 'theta' column will be replaced with TH1, TH2, TH3 in the
 380     # general case. (Actually it will be replaced with the
 381     # thetas labels from the model file. Here we search the meta
 382     # column of the raw_results_header to find the maximum
 383     # number of real columns.
 384
 385     my %max_hash;
 386
 387     foreach my $model ( @{$self -> {'models'}} ){
 388
 389       foreach my $category ( @{$self -> {'raw_results_header'}},'npomega' ) {
 390         if ( $category eq 'setheta' or $category eq 'seomega' or $category eq 'sesigma' ){
 391           next;
 392         } elsif ( $category eq 'theta' or $category eq 'omega' or $category eq 'sigma' or
 393                   $category eq 'npomega' or $category eq 'shrinkage_etas' or $category eq 'eigen') {
 394           my $numpar = 0;
 395           if( $category eq 'npomega' or $category eq 'shrinkage_etas' ) {
 396             my $nomegas = $model -> nomegas(with_correlations => 1);
 397             if( defined $nomegas ) {
 398               for( my $j = 0; $j < scalar @{$nomegas}; $j++ ) {
 399                 if( $category eq 'npomega' ) {
 400                   my $npar = $nomegas -> [$j];
 401                   $npar = $npar*($npar+1)/2;
 402                   $numpar = $numpar >= $npar ? $numpar : $npar;
 403                 } else {
 404                   $numpar = $numpar >= $nomegas -> [$j] ? $numpar : $nomegas -> [$j];
 405                 }
 406               }
 407             }
 408           } elsif( $category eq 'eigen') {
 409
 410             # This is an upper limit on the number of eigenvalues in
 411             # the output file. The accessors should not count "SAME"
 412             # but should count offdiagonals. It also counts "FIX"
 413             # which is not ideal.
 414
 415             my $max_sigmas = 0;
 416             foreach my $prob( @{$model -> nsigmas(with_correlations => 1 )} ) {
 417               if( $prob > $max_sigmas ){
 418                 $max_sigmas = $prob;
 419               }
 420             }
 421
 422             my $max_omegas = 0;
 423             foreach my $prob( @{$model -> nomegas(with_correlations => 1 )} ) {
 424               if( $prob > $max_omegas ){
 425                 $max_omegas = $prob;
 426               }
 427             }
 428
 429             $numpar = $model -> nthetas() + $max_omegas + $max_sigmas;
 430
 431           } else {
 432
 433             # Here we assume that $category is 'theta', 'omega' or
 434             # 'sigma'. We also asume that 'SAME' estimates have a zero
 435             # size.
 436
 437             my $accessor = 'n'.$category.'s';
 438             if( $category eq 'theta' ){
 439               $numpar = $model -> $accessor();
 440             } else {
 441               # Here accessor must be omega or sigma.
 442               $numpar = $model -> $accessor(with_correlations => 1);
 443             }
 444
 445             # omega and sigma is an array, we find the biggest member
 446
 447             if( ref($numpar) eq 'ARRAY' ){
 448               my $max = 0;
 449               foreach my $prob ( @{$numpar} ){
 450                 if( $prob > $max ){
 451                   $max = $prob;
 452                 }
 453               }
 454               $numpar = $max;
 455             }
 456
 457           }
 458           if ( $max_hash{ $category } < $numpar ) {
 459             $max_hash{ $category } = $numpar;
 460             $max_hash{ 'se'.$category } = $numpar;
 461           }
 462         } else {
 463           $max_hash{ $category } = 1;
 464         }
 465       }
 466     }
 467
 468     $self -> {'max_hash'} = \%max_hash;
 469
 470     if( 0 ) {
 471
 472       # {{{ Create a header
 473       my @params = ( 'theta', 'omega', 'sigma', 'npomega', 'shrinkage_etas', 'shrinkage_wres','eigen' );
 474       my @new_header;
 475       foreach my $category ( @{$self -> {'raw_results_header'}},'npomega' ){
 476         my @new_names = ();
 477         foreach my $param ( @params ) {
 478           if ( $category eq $param ) {
 479             if( defined $max_hash{$param} ){
 480               for ( my $i = 1; $i <= $max_hash{ $category }; $i++ ) {
 481                 push ( @new_names, uc(substr($category,0,2)).$i );
 482               }
 483             }
 484             last;
 485           }
 486           if ( $category eq 'se'.$param ) {
 487             if( defined $max_hash{$param} ){
 488               for ( my $i = 1; $i <= $max_hash{ $category }; $i++ ) {
 489                 push ( @new_names, uc(substr($category,2,2)).$i );
 490               }
 491               map ( $_ = 'se'.$_, @new_names );
 492             }
 493             last;
 494           }
 495         }
 496         if ( $#new_names >= 0 ) {
 497           push( @new_header, @new_names );
 498         } else {
 499           push( @new_header, $category );
 500         }
 501
 502       }
 503
 504       # }}}
 505
 506     }
 507
 508   }
 509 end calculate_raw_results_width
 510
 511 # }}} calculate_raw_results_width
 512
 513 # {{{ register_in_database
 514
 515 start register_in_database
 516   {
 517     if ( $PsN::config -> {'_'} -> {'use_database'} ) {
 518       my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
 519                                ";databse=".$PsN::config -> {'_'} -> {'project'},
 520                                $PsN::config -> {'_'} -> {'user'},
 521                                $PsN::config -> {'_'} -> {'password'},
 522                                {'RaiseError' => 1});
 523       my $sth;
 524
 525       # register modelfit (execute) tool
 526 #      print("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
 527 #           ".execute ( comment ) ".
 528 #           "VALUES ( '' )");
 529       $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
 530                              ".execute ( comment ) ".
 531                              "VALUES ( 'test' )");
 532       $sth -> execute;
 533       $self -> execute_id($sth->{'mysql_insertid'});
 534
 535       $sth -> finish;
 536       $dbh -> disconnect;
 537       tool::register_in_database( $self, execute_id => $self -> execute_id() );
 538
 539     }
 540   }
 541 end register_in_database
 542
 543 # }}} register_in_database
 544
 545 # {{{ prepare_raw_results
 546
 547 start prepare_raw_results
 548       {
 549         # As of PsN version 2.1.8, we don't handle problems and
 550         # subproblems in any of the tools but modelfit.
 551
 552         @{$self -> {'raw_results'}} = sort( {$a->[0] <=> $b->[0]} @{$self -> {'raw_results'}} );
 553
 554         my %max_hash = %{$self -> {'max_hash'}};
 555
 556         &{$self -> {'_raw_results_callback'}}( $self, \%max_hash )
 557             if ( defined $self -> {'_raw_results_callback'} );
 558
 559         # ---------------------------  Create a header  ----------------------------
 560
 561         # {{{ header
 562
 563         my %param_names;
 564         my @params = ( 'theta', 'omega', 'sigma' );
 565         foreach my $param ( @params ) {
 566           my $labels = $self -> models -> [0] -> labels( parameter_type => $param );
 567           $param_names{$param} = $labels -> [0] if ( defined $labels );
 568         }
 569
 570         my @new_header;
 571         foreach my $name ( @{$self -> {'raw_results_header'}} ) {
 572           my $success;
 573           my @params = ( 'theta', 'omega', 'sigma', 'npomega', 'shrinkage_etas', 'shrinkage_wres','eigen' );
 574           foreach my $param ( @params ) {
 575             if ( $name eq $param ){
 576               if ( $name eq 'shrinkage_etas' ){
 577                 for ( my $i = 1; $i <= $max_hash{ $name }; $i++ ) {
 578                   push ( @new_header, 'shrinkage_eta'.$i );
 579                 }
 580               } elsif ( $name eq 'shrinkage_wres' ){
 581                 push ( @new_header, 'shrinkage_wres' );
 582               } else {
 583                 for ( my $i = 1; $i <= $max_hash{ $name }; $i++ ) {
 584
 585                   my $label = undef;
 586                   if( defined $param_names{$name} -> [$i-1] ){
 587                     $label = $param_names{$name} -> [$i-1] ;
 588                   }
 589                   if( defined $label ){
 590                     push( @new_header, $label );
 591                   } else {
 592                     push( @new_header, uc(substr($name,0,2)).$i );
 593                   }
 594                 }
 595               }
 596               $success = 1;
 597               last;
 598             } elsif ( $name eq 'se'.$param ) {
 599               for ( my $i = 1; $i <= $max_hash{ $name }; $i++ ) {
 600                 my $label = undef;
 601
 602                 if( defined $param_names{$param} -> [$i-1] ){
 603                   $label = $param_names{$param} -> [$i-1] ;
 604                 }
 605                 if( defined $label ){
 606                   push( @new_header, 'se' . $label );
 607                 } else {
 608                   push( @new_header, 'se' . uc(substr($name,2,2)).$i );
 609                 }
 610               }
 611               $success = 1;
 612               last;
 613             }
 614           }
 615           unless( $success ){
 616             push( @new_header, $name );
 617           }
 618         }
 619
 620         $self -> {'raw_results_header'} = \@new_header;
 621
 622         # }}} header
 623       }
 624 end prepare_raw_results
 625
 626 # }}} prepare_raw_results
 627
 628 # {{{ print_raw_results
 629
 630 start print_raw_results
 631 {
 632
 633   ## print raw_results array
 634
 635   if ( defined $self -> {'raw_results'} ) {
 636
 637     my $raw_file;
 638
 639     if( ref $self -> {'raw_results_file'} eq 'ARRAY' ){
 640       $raw_file = $self -> {'raw_results_file'} -> [0];
 641     } else {
 642       $raw_file = $self -> {'raw_results_file'};
 643     }
 644
 645     my ($dir,$file) = OSspecific::absolute_path( $self -> {'directory'},
 646                                                  $raw_file );
 647
 648     my $append = $self -> {'raw_results_append'} ? '>>' : '>';
 649     open( RRES, $append.$dir.$file );
 650
 651     if( (not $self -> {'raw_results_append'}) and $PsN::output_header ){
 652       print RRES join(',',@{$self -> {'raw_results_header'}} ),"\n";
 653     }
 654
 655     for ( my $i = 0; $i < scalar @{$self -> {'raw_results'}}; $i++ ) {
 656       my @result_array = @{$self -> {'raw_results'} -> [$i]};
 657       map( $_ = defined $_ ? $_ : $PsN::out_miss_data, @result_array );
 658       print RRES join(',',@result_array ),"\n";
 659     }
 660     close( RRES );
 661   }
 662
 663   ## print raw_nonp_results
 664
 665   if( defined $self ->{'raw_nonp_results' } ) {
 666     my ($dir,$file) = OSspecific::absolute_path( $self -> {'directory'},
 667                                                  $self -> {'raw_nonp_file'} );
 668     my $append = $self -> {'raw_results_append'} ? '>>' : '>';
 669     open( RRES, $append.$dir.$file );
 670     for ( my $i = 0; $i < scalar @{$self -> {'raw_nonp_results'}}; $i++ ) {
 671       my @result_array;
 672       if ( defined $self -> {'raw_nonp_results'} -> [$i] ) {
 673         @result_array = @{$self -> {'raw_nonp_results'} -> [$i]};
 674         map( $_ = defined $_ ? $_ : $PsN::out_miss_data, @result_array );
 675       }
 676       print RRES join(',',@result_array ),"\n";
 677     }
 678     close( RRES );
 679   }
 680
 681   ## raw_line_structure should be printed to disk here for fast
 682   ## resumes. In the future
 683
 684   $self -> {'raw_line_structure'} -> write( 'raw_results_structure' );
 685
 686 }
 687 end print_raw_results
 688
 689 # }}} print_raw_results
 690
 691 # {{{ create_sub_dir
 692
 693 start create_sub_dir
 694     {
 695         my $file;
 696 #       ($tmp_dir, $file) = OSspecific::absolute_path( './' .
 697         ($tmp_dir, $file) = OSspecific::absolute_path( $self -> {'directory'}.'/' .
 698                                                        $subDir, '');
 699
 700         unless( -e $tmp_dir ){
 701             mkdir( $tmp_dir );
 702         }
 703     }
 704 end create_sub_dir
 705
 706 # }}}
 707
 708 # {{{ copy_model_and_input
 709
 710 start copy_model_and_input
 711     {
 712       # Fix new short names (i.e. No path)
 713       my @new_data_names;
 714
 715       if( defined $model -> datas ){
 716         foreach my $data ( @{$model -> datas} ) {
 717           my $filename = $data -> filename;
 718
 719           push( @new_data_names, $filename );
 720         }
 721       } else {
 722         debug -> die( message => 'No datafiles set in modelfile.' );
 723       }
 724
 725       # Fix new, short names (i.e. No Path)
 726       my @new_extra_data_names;
 727       my @problems = @{$model -> problems};
 728       for ( my $i = 1; $i <= $#problems + 1; $i++ ) {
 729         my $extra_data = $problems[$i-1] -> extra_data;
 730         if ( defined $extra_data ) {
 731           my $filename = $extra_data -> filename;
 732           push( @new_extra_data_names, $filename );
 733         }
 734       }
 735
 736       # Set the table names to a short version
 737       my @new_table_names = ();
 738       my @table_names = @{$model -> table_names( ignore_missing_files => 1 )};
 739       # Loop the problems
 740       for ( my $i = 0; $i <= $#table_names; $i++ ) {
 741         my @new_arr;
 742         # Loop the table files within each problem
 743         for ( my $j = 0; $j < scalar @{$table_names[$i]}; $j++ ) {
 744           my ( $dir, $filename ) = OSspecific::absolute_path( '.', $table_names[$i][$j] );
 745           push( @new_arr,  $filename );
 746         }
 747         push( @new_table_names, \@new_arr );
 748       }
 749
 750
 751       # Copy input files ( msfo, msfi, subroutines and extra at the
 752       # moment)
 753
 754       foreach my $file( @{$model -> input_files} ){
 755
 756         # $file is a ref to an array with two elements, the first is a
 757         # path, the second is a name.
 758
 759         cp( $file->[0] . $file -> [1], $file -> [1] );
 760
 761       }
 762
 763
 764       # Copy the model object. Set the new (shorter) data file names.
 765       # There's no need to physically copy these here since we took care of that above.
 766       $candidate_model = $model -> copy( filename              => 'psn.mod',
 767                                      data_file_names       => \@new_data_names,
 768                                      extra_data_file_names => \@new_extra_data_names,
 769                                      copy_data             => 1 );
 770
 771       $candidate_model -> register_in_database;
 772       $model -> flush_data;
 773
 774       if( $self -> {'handle_msfo'} ){
 775
 776         # Initialize sequence of msfi/msfo files.
 777
 778         my $msfo_names = $candidate_model -> msfo_names;
 779         my $msfi_names = $candidate_model -> msfi_names;
 780         my $msfi_in;
 781
 782         print "Model already has msfi: ". Dumper $msfi_names . "and msfo: " . Dumper $msfo_names;
 783
 784         if( defined $msfo_names ){
 785           $msfi_in = $msfo_names -> [0][0];
 786         } elsif ( defined $msfi_names ){
 787           $msfi_in = $msfi_names -> [0][0];
 788         }
 789
 790         if( -e $msfi_in ){
 791           mv( $msfi_in, 'psn_msfo-0' );
 792           $candidate_model->set_records(type=>'msfi',
 793                                         record_strings => ['psn_msfo-0']);
 794           $candidate_model->remove_records(type=>'theta');
 795           $candidate_model->remove_records(type=>'omega');
 796           $candidate_model->remove_records(type=>'sigma');
 797
 798         } else {
 799           # Probably an error, will be caught by NONMEM
 800         }
 801
 802         if( scalar @{$candidate_model -> record(record_name => 'estimation')} > 0 ){
 803           $candidate_model -> set_option( record_name => 'estimation',
 804                                           option_name => 'MSFO',
 805                                           option_value=> 'psn_msfo' );
 806         }
 807
 808       }
 809
 810       my $maxeval = $model -> maxeval -> [0][0];
 811       if ( $maxeval > 9999 ) {
 812         $candidate_model -> maxeval( new_values => [[9999]],
 813                                  problem_numbers => [1] );
 814       }
 815
 816       $candidate_model -> table_names( new_names            => \@new_table_names,
 817                                    ignore_missing_files => 1 );
 818 #      $candidate_model -> clean_extra_data_code;
 819       $candidate_model -> drop_dropped if ( $self -> {'drop_dropped'} );
 820       $candidate_model -> wrap_data if ( $self -> {'wrap_data'} );
 821       $candidate_model -> add_extra_data_code;
 822       $candidate_model -> write_get_subs;
 823       $candidate_model -> write_readers;
 824       $candidate_model -> _write( filename   => 'psn.mod' );# write_data => 1 );  #Kolla denna, den funkar inte utan wrap!!
 825       $candidate_model -> flush_data;
 826       $candidate_model -> store_inits;
 827
 828     }
 829 end copy_model_and_input
 830
 831 # }}}
 832
 833 # {{{ copy_model_and_output
 834
 835 start copy_model_and_output
 836 {
 837   my $outfilename = $model -> outputs -> [0] -> full_name;
 838
 839   my ($dir, $model_filename) = OSspecific::absolute_path($model -> directory,
 840                                                          $model -> filename );
 841
 842   # This is used with 'prepend_model_file_name'
 843   my $dotless_model_filename = $model_filename;
 844   $dotless_model_filename =~ s/\.[^.]+$//;
 845
 846   if( $self -> unwrap_table_files() ) {
 847     if( defined $final_model -> table_names ){
 848       foreach my $table_files( @{$final_model -> table_names} ){
 849         foreach my $table_file( @{$table_files} ){
 850
 851           open( TABLE, '<'.$table_file );
 852           my @table = <TABLE>;
 853           close( TABLE );
 854           open( TABLE, '>'.$table_file );
 855           my ( $j, $cont_column ) = ( 0, -1 );
 856           for( my $i = 0; $i <= $#table; $i++ ) {
 857             print TABLE $table[$i] and next if( $i == 0 );
 858             chomp($table[$i]);
 859             my @row = split(' ',$table[$i]);
 860             if( $i == 1 ) {
 861               for( my $k = 0; $k <= $#row; $k++ ) {
 862                 $cont_column = $k if( $row[$k] eq 'CONT' );
 863               }
 864             }
 865             for( my $k = 0; $k <= $#row; $k++ ) {
 866               next if( $k == $cont_column );
 867               print TABLE sprintf( "%12s",$row[$k] );
 868             }
 869             print TABLE "\n";
 870           }
 871           close( TABLE );
 872         }
 873       }
 874     }
 875     $final_model -> table_names( new_names => $model -> table_names );
 876   }
 877
 878   my @output_files = @{$final_model -> output_files};
 879
 880   foreach my $filename ( @output_files, 'compilation_output.txt' ){
 881
 882     my $use_name = $self -> get_retry_name( filename => $filename,
 883                                             retry => $use_run-1 );
 884
 885     # Copy $use_run files to files without numbers, to avoid confusion.
 886     cp( $use_name, $filename );
 887
 888     # Don't prepend the model file name to psn.lst, but use the name
 889     # from the $model object.
 890     if( $filename eq 'psn.lst' ){
 891       cp( $use_name, $outfilename );
 892       next;
 893     }
 894
 895     if( $self -> {'prepend_model_file_name'} ){
 896       cp( $use_name, $dir . $dotless_model_filename . '.' . $filename );
 897     } else {
 898       cp( $use_name, $dir .$filename );
 899     }
 900
 901   }
 902
 903   # TODO check if this is necessary
 904   my $final_output = output -> new( filename => $outfilename,
 905                                     model_id => $model -> model_id );
 906   $final_output -> register_in_database( model_id => $model -> model_id,
 907                                          force    => 1 ); # If we are here, the model has been run
 908   $final_model -> outputs( [$final_output] );
 909
 910
 911   # Keep files if debugging
 912
 913   if( 'debug' -> level == 0) {
 914     unlink 'nonmem', 'nonmem6', 'nonmem5',
 915            'nonmem.exe', 'nonmem5_adaptive','nonmem6_adaptive', 'nonmem_adaptive',
 916            'FDATA';
 917     # TODO
 918     # If we delete these, we can't resume /Lasse :
 919     # (pheno is not a good testing example for this)
 920     # unlink( @{$model -> datafiles}, @{$model -> extra_data_files} );
 921   }
 922
 923   if( $self -> {'clean'} >= 1 and 'debug' -> level == 0 ){
 924     unlink 'nonmem', 'nonmem'.$self -> {'nm_version'},
 925            'nonmem.exe','FDATA', 'FREPORT', 'FSUBS', 'FSUBS.f',
 926            'FSUBS.for', 'LINK.LNK', 'FSTREAM', 'FCON.orig', 'FLIB', 'FCON','PRDERR';
 927
 928
 929     if( defined $final_model -> extra_files ){
 930       foreach my $x_file( @{$final_model -> extra_files} ){
 931         my ( $dir, $filename ) = OSspecific::absolute_path( $final_model -> directory,
 932                                                             $x_file );
 933         unlink( $filename );
 934       }
 935     }
 936
 937
 938     if( $self -> {'clean'} >= 2 ){
 939       for ( my $i = 1; $i <= $self -> {'retries'}; $i++ ) {
 940         foreach my $filename ( @output_files ){
 941
 942           my $use_name = $self -> get_retry_name( filename => $filename,
 943                                                   retry => $i-1 );
 944           unlink( $use_name );
 945         }
 946
 947         unlink( "psn-$i.mod" );
 948         unlink( "compilation_output-$i.txt." );
 949
 950       }
 951       unlink( @{$model -> datafiles}, @{$model -> extra_data_files} );
 952     }
 953   }
 954
 955   if ( $self -> {'clean'} >= 3 ) {
 956     # Do nothing. "run_nonmem" will remove entire work directory
 957     # before returning.
 958   } else {
 959     system('tar cz --remove-files -f nonmem_files.tgz *')
 960         if ( $self -> {'compress'} and $Config{osname} ne 'MSWin32' );
 961     system('compact /c /s /q > NUL')
 962         if ( $self -> {'compress'} and $Config{osname} eq 'MSWin32' );
 963   }
 964 }
 965 end copy_model_and_output
 966
 967 # }}}
 968
 969 # {{{ get_retry_name
 970
 971 start get_retry_name
 972 {
 973   $retry++;
 974   unless( $filename =~ s/\.([^.]+)$/-$retry.$1/ ){
 975     $filename .= "-$retry";
 976   }
 977 }
 978 end get_retry_name
 979
 980 # }}}
 981
 982 # {{{ set_msfo_to_msfi
 983
 984 start set_msfo_to_msfi
 985     {
 986
 987       my $msfo = $self -> get_retry_name( 'filename' => 'psn_msfo',
 988                                           'retry' => $retry );
 989
 990       my $msfi;
 991
 992       if( $candidate_model -> outputs -> [0] -> msfo_has_terminated() ){
 993
 994         print "Msfo terminated\n";
 995
 996         $msfi = $msfo . '-' . ($queue_info -> {'crashes'}-1);
 997
 998         $candidate_model->remove_records( type => 'estimation' );
 999
1000       } else {
1001
1002         $msfi = $msfo . '-' . $queue_info -> {'crashes'};
1003
1004       }
1005
1006       print "$msfo -> $msfi\n";
1007       if( -e $msfo ){
1008         mv( $msfo, $msfi );
1009       }
1010
1011       $candidate_model->set_records(type=>'msfi',
1012                                     record_strings => [$msfi]);
1013
1014       $candidate_model->remove_records(type=>'theta');
1015       $candidate_model->remove_records(type=>'omega');
1016       $candidate_model->remove_records(type=>'sigma');
1017       $candidate_model->_write;
1018
1019     }
1020 end set_msfo_to_msfi
1021
1022 # }}}
1023
1024 # {{{ reset_msfo
1025
1026 start reset_msfo
1027     {
1028       my @data_ref = @{$candidate_model -> record( record_name => 'msfi' )};
1029       # Check if there is a msfi record and then delete it
1030       if (scalar(@data_ref)!=0) {
1031         $candidate_model->remove_records(type=>'msfi');
1032
1033         # Set the intial values + boundaries to the first  values (update theta, omega, sigma)
1034
1035         my @old_problems = @{$basic_model -> problems};
1036         my @new_problems = @{$candidate_model -> problems};
1037         for ( my $i=0; $i <= $#old_problems; $i++ ) {
1038           foreach my $param ( 'thetas', 'omegas', 'sigmas' ) {
1039             $new_problems[$i] -> $param( Storable::dclone( $old_problems[$i] -> $param ) );
1040           }
1041         }
1042
1043         $model_modified = 1;
1044       }
1045
1046     }
1047 end reset_msfo
1048
1049 # }}}
1050
1051 # {{{ cut_thetas
1052
1053 start cut_thetas
1054     {
1055       $candidate_model -> update_inits( from_output => $output_file,
1056                                     update_omegas => 1,
1057                                     update_sigmas => 1,
1058                                     update_thetas => 1);
1059
1060       foreach my $th_num ( @cutoff_thetas ) {
1061         my $init_val = $candidate_model ->
1062             initial_values( parameter_type    => 'theta',
1063                             parameter_numbers => [[$th_num]])->[0][0];
1064         if (abs($init_val)<=$self->{'cutoff'}) {
1065           $candidate_model->initial_values(parameter_type => 'theta',
1066                                        parameter_numbers => [[$th_num]],
1067                                        new_values =>[[0]]);
1068           $candidate_model->fixed(parameter_type => 'theta',
1069                               parameter_numbers => [[$th_num]],
1070                               new_values => [[1]] );
1071         }
1072       }
1073     }
1074 end cut_thetas
1075
1076 # }}}
1077
1078 # {{{ ask_user
1079
1080 start ask_user
1081     {
1082       # rand should not be used here. find some other way to create unique file names.
1083       my $num = rand;
1084       open( TMP, ">/tmp/$num" );
1085       print TMP "START MODEL FILE NAME\n";
1086       print TMP $basic_model -> filename,"\n";
1087       print TMP "END MODEL FILE NAME\n";
1088       foreach my $prob ( @reruns ) {
1089         my @theta_labels = @{$candidate_model -> labels( parameter_type => 'theta' )};
1090         my @omega_labels = @{$candidate_model -> labels( parameter_type => 'omega' )};
1091         my @sigma_labels = @{$candidate_model -> labels( parameter_type => 'sigma' )};
1092         my @theta_inits = @{$candidate_model -> initial_values( parameter_type => 'theta' )};
1093         my @omega_inits = @{$candidate_model -> initial_values( parameter_type => 'omega' )};
1094         my @sigma_inits = @{$candidate_model -> initial_values( parameter_type => 'sigma' )};
1095         print TMP "START PROBLEM NUMBER\n";
1096         print TMP $prob,"\n";
1097         print TMP "END PROBLEM NUMBER\n";
1098         print TMP "START MINIMIZATION_MESSAGE\n";
1099         print TMP $minimization_message -> [$prob-1][0],"\n";
1100         print TMP "END MINIMIZATION_MESSAGE\n";
1101         print TMP "START FINAL GRADIENT\n";
1102         print TMP join( " ",@{$output_file -> final_gradients -> [$prob-1][0]}),"\n";
1103         print TMP "END FINAL GRADIENT\n";
1104         print TMP "START OFV\n";
1105         print TMP $output_file -> ofv -> [$prob-1][0],"\n";
1106         print TMP "END OFV\n";
1107         print TMP "START INITIAL VALUES THETA\n";
1108         print TMP join(" ", @{$theta_inits[$prob-1]}),"\n";
1109         print TMP "END INITIAL VALUES THETA\n";
1110         print TMP "START INITIAL VALUES OMEGA\n";
1111         print TMP join(" ", @{$omega_inits[$prob-1]}),"\n";
1112         print TMP "END INITIAL VALUES OMEGA\n";
1113         print TMP "START INITIAL VALUES SIGMA\n";
1114         print TMP join(" ", @{$sigma_inits[$prob-1]}),"\n";
1115         print TMP "END INITIAL VALUES SIGMA\n";
1116         print TMP "START LABELS\n";
1117         print TMP join(" ", (@{$theta_labels[$prob-1]},@{$omega_labels[$prob-1]},
1118                              @{$sigma_labels[$prob-1]})),"\n";
1119         print TMP "END LABELS\n";
1120       }
1121       close( TMP );
1122       my $out = readpipe( "/users/lasse/PsN/Diagrams/test/scm_comm.pl $num ".$output_file->filename );
1123       my @out_per_prob = split("\n",$out);
1124       foreach my $prob ( @reruns ) {
1125         my ( $choice, $rest ) = split( ' ', shift( @out_per_prob ), 2 );
1126         if ( $choice == 0 ) {
1127           $retries = $tries + $rest;
1128         } elsif ( $choice == 1 ) {
1129           my ($theta_str,$omega_str,$sigma_str) = split(':',$rest);
1130           print "thstr $theta_str\n";
1131           print "omstr $omega_str\n";
1132           print "sistr $sigma_str\n";
1133           my @thetas = split( ' ', $theta_str );
1134           print "$prob: @thetas\n";
1135           $candidate_model -> initial_values( parameter_type  => 'theta',
1136                                           problem_numbers => [$prob],
1137                                           new_values      => [\@thetas],
1138                                           add_if_absent   => 0 );
1139           $retries = $tries+1;
1140         } else {
1141           last RETRY;
1142         }
1143         $candidate_model -> _write;
1144       }
1145       $return_value = $retries;
1146     }
1147 end ask_user
1148
1149 # }}}
1150
1151 # {{{ umbrella_submit
1152
1153 start umbrella_submit
1154     {
1155         if( $prepare_jobs ){
1156
1157             my $fsubs = join( ',' , @{$model -> subroutine_files} );
1158             my $nm_command = ($PsN::config -> {'_'} -> {'remote_perl'} ? $PsN::config -> {'_'} -> {'remote_perl'} : 'perl') . " -I" .
1159                 $PsN::lib_dir ."/../ " .
1160                 $PsN::lib_dir . "/nonmem.pm" .
1161                 " psn.mod psn.lst " .
1162                 $self -> {'nice'} . " ".
1163                 $nm_version . " " .
1164                 "1 " .
1165                 "1 " .
1166                 $fsubs . " " .
1167                 $self -> {'nm_directory'};
1168
1169             my $directory = $model -> directory();
1170
1171             push( @{$self -> {'__umbrella_insert'}}, [$nm_command,$directory] );
1172
1173             $queue_info->{$directory}{'candidate_model'} = $model;
1174
1175         } else {
1176
1177             my ($nonmem_insert, $tool_insert, $job_insert);
1178
1179             foreach my $row( @{$self -> {'__umbrella_insert'}} ){
1180                 $nonmem_insert .= ',' if( defined $nonmem_insert );
1181                 $nonmem_insert .= "('". $row -> [0] ."')";
1182             }
1183
1184             my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
1185                                      ";databse=".$PsN::config -> {'_'} -> {'project'},
1186                                      $PsN::config -> {'_'} -> {'user'},
1187                                      $PsN::config -> {'_'} -> {'password'},
1188                                      {'RaiseError' => 1});
1189
1190             my $project = $PsN::config -> {'_'} -> {'project'};
1191             $dbh -> do( "LOCK TABLES $project.nonmem WRITE,$project.tool WRITE,$project.job WRITE");
1192
1193             my $sth = $dbh -> prepare( "SELECT max(nonmem_id) FROM $project.nonmem" );
1194             $sth -> execute or 'debug' -> die( message => $sth -> errstr );
1195
1196             my $nonmem_old_max = $sth -> fetch;
1197
1198             # register nonmem tool
1199             $sth -> finish;
1200             my $sth = $dbh -> prepare( "INSERT INTO $project.nonmem ( command ) ".
1201                                        "VALUES $nonmem_insert");
1202
1203             $sth -> execute or 'debug' -> die( message => $sth->errstr ) ;
1204             $sth -> finish;
1205
1206             my $sth = $dbh -> prepare( "SELECT max(nonmem_id) FROM $project.nonmem" );
1207             $sth -> execute or 'debug' -> die( message => $sth -> errstr );
1208             my $nonmem_new_max = $sth -> fetch;
1209
1210             my $tool_insert;
1211
1212             for( my $i = ($nonmem_old_max -> [0] + 1); $i <= $nonmem_new_max -> [0]; $i++ ){
1213                 $tool_insert .= ',' if (defined $tool_insert);
1214                 $tool_insert .= "('".$self -> {'tool_id'}."','".$i."','".
1215                     $self -> {'__umbrella_insert'} -> [ $i-($nonmem_old_max->[0]+1) ]->[1]. "')";
1216             }
1217             $sth -> finish;
1218
1219             my $sth = $dbh -> prepare( "SELECT max(tool_id) FROM $project.tool" );
1220             $sth -> execute or 'debug' -> die( message => $sth -> errstr );
1221
1222             my $tool_old_max = $sth -> fetch;
1223
1224             # register generic tool
1225             $sth -> finish;
1226             my $sth = $dbh -> prepare( "INSERT INTO $project.tool ( parent_tool_id, nonmem_id, directory ) ".
1227                                        "VALUES $tool_insert" );
1228             $sth -> execute or 'debug' -> die( message => $sth->errstr ) ;
1229             $sth -> finish;
1230             my $sth = $dbh -> prepare( "SELECT max(tool_id) FROM $project.tool" );
1231             $sth -> execute or 'debug' -> die( message => $sth -> errstr );
1232
1233             my $tool_new_max = $sth -> fetch;
1234
1235             my $job_insert;
1236
1237             for( my $i = ($tool_old_max -> [0] + 1); $i <= $tool_new_max -> [0]; $i++ ){
1238                 $job_insert .= ',' if (defined $job_insert);
1239                 $job_insert .= "('" . $i . "','1')";
1240             }
1241
1242             # register job
1243
1244             $sth -> finish;
1245             my $sth = $dbh -> prepare( "INSERT INTO $project.job (tool_id, state) ".
1246                                        "VALUES $job_insert");
1247             $sth -> execute or 'debug' -> die( message => $sth->errstr ) ;
1248
1249
1250             $dbh -> do( "UNLOCK TABLES" );
1251
1252             $sth -> finish;
1253             $dbh -> disconnect;
1254         }
1255
1256     }
1257 end umbrella_submit
1258
1259 # }}}
1260
1261 # {{{ lsf_submit
1262
1263 start lsf_submit
1264     {
1265
1266       # This method will submit the nonmem.pm file as a script to the
1267       # LSF system.
1268
1269       my $fsubs = join( ',' , @{$model -> subroutine_files} );
1270
1271       # Check for vital lsf options.
1272       unless( $self -> {'lsf_queue'} ){
1273         if( $PsN::config -> {'_'} -> {'lsf_queue'} ){
1274           $self -> {'lsf_queue'} = $PsN::config -> {'_'} -> {'lsf_queue'};
1275         } else {
1276           'debug' -> die( message => 'No queue specified for lsf run' );
1277         }
1278       }
1279
1280       foreach my $lsf_option ( 'lsf_project_name', 'lsf_job_name', 'lsf_resources', 'lsf_ttl', 'lsf_options' ){
1281         unless( $self -> {$lsf_option} ){
1282           if( $PsN::config -> {'_'} -> {$lsf_option} ){
1283             $self -> {$lsf_option} = $PsN::config -> {'_'} -> {$lsf_option};
1284           }
1285         }
1286       }
1287
1288       my ($lsf_out, $lsf_err);
1289       for( my $i = 1; $i <= 5; $i++ ){
1290         my $str = "bsub -e stderr -o stdout " .
1291             "-q " . $self -> {'lsf_queue'} .
1292             ($self -> {'lsf_project_name'} ? " -P " . $self -> {'lsf_project_name'} : ' ') .
1293             ($self -> {'lsf_job_name'} ? " -J " . $self -> {'lsf_job_name'} : ' ') .
1294             ($self -> {'lsf_ttl'} ? " -c " . $self -> {'lsf_ttl'} : ' ') .
1295             ($self -> {'lsf_resources'} ? " -R " . $self -> {'lsf_resources'} : ' ') .
1296             $self -> {'lsf_options'} . " \"sleep 3 && " .
1297             ($PsN::config -> {'_'} -> {'remote_perl'} ? ' ' . $PsN::config -> {'_'} -> {'remote_perl'} : ' perl ') . " -I" .
1298             $PsN::lib_dir ."/../ " .
1299             $PsN::lib_dir . "/nonmem.pm" .
1300             " psn.mod psn.lst " .
1301             $self -> {'nice'} . " ".
1302             $nm_version . " " .
1303             1 . ' ' .
1304             1 . ' ' .
1305             $fsubs . " " .
1306             $self -> {'nm_directory'} . "\"";
1307
1308         run3( $str, undef, \$lsf_out, \$lsf_err );
1309         if ($lsf_out=~/Job \<(\d+)\> is submitted to queue/) {
1310           $jobId=$1;
1311         }
1312
1313         unless( $lsf_err =~ /System call failed/ or
1314                 $lsf_err =~ /Bad argument/ or
1315                 $lsf_err =~ /Request aborted by esub/ or
1316                 $lsf_err =~ /Bad user ID/ ) {
1317           last;
1318         }
1319
1320         print "$lsf_err\n";
1321         if( $lsf_err =~ /Bad argument/ or
1322             $lsf_err =~ /Request aborted by esub/ or
1323             $lsf_err =~ /Bad user ID/ ) {
1324           sleep(($i+1)**2);
1325         } else {
1326           chdir( $work_dir );
1327         }
1328         print "bsub command was not successful, trying ",(5-$i)," times more\n";
1329       }
1330       sleep(3);
1331     }
1332 end lsf_submit
1333
1334 # }}}
1335
1336 # {{{ lsf_monitor
1337
1338 start lsf_monitor
1339     {
1340       my ($stdout, $stderr);
1341       run3("bjobs $jobId",undef,\$stdout, \$stderr );
1342
1343       if ($stdout=~/DONE/m) {
1344         return $jobId; # Return the jobId found.
1345       }
1346
1347       return 0;
1348     }
1349 end lsf_monitor
1350
1351 # }}}
1352
1353 # {{{ ud_submit
1354
1355 start ud_submit
1356 {
1357   my $script;
1358   unless( defined $PsN::config -> {'_'} -> {'ud_nonmem'} ){
1359     if( $Config{osname} eq 'MSWin32' ) {
1360       $script = 'nonmem.bat';
1361     } else {
1362       $script = 'nonmem.sh';
1363     }
1364   } else {
1365     $script = $PsN::config -> {'_'} -> {'ud_nonmem'};
1366   }
1367
1368   if( system( "$script -s " . $model -> filename . "> nonmem_sh_stdout"  ) ){
1369     'debug' -> die( message => "UD submit script failed, check that $script is in your PATH.\nSystem error message: $!" );
1370   }
1371
1372   open(JOBFILE, "JobId") or 'debug' -> die( message => "Couldn't open UD grid JobId file for reading: $!" );
1373   $jobId = <JOBFILE>;
1374   close(JOBFILE);
1375
1376 }
1377 end ud_submit
1378
1379 # }}}
1380
1381 # {{{ ud_monitor
1382
1383 start ud_monitor
1384 {
1385 # unless( $self -> {'ud_native_retrieve'} ){
1386
1387   my $script;
1388   unless( defined $PsN::config -> {'_'} -> {'ud_nonmem'} ){
1389     if( $Config{osname} eq 'MSWin32' ) {
1390       $script = 'nonmem.bat';
1391     } else {
1392       $script = 'nonmem.sh';
1393     }
1394   } else {
1395     $script = $PsN::config -> {'_'} -> {'ud_nonmem'};
1396   }
1397
1398
1399   my $stdout; # Variable to store output from "nonmem.bat"
1400
1401   unless( run3( "$script -l " . $jobId, undef, \$stdout ) ){ # run3 will capture the output
1402     'debug' -> die( message => "UD submit script failed, check that $script is in your PATH.\nSystem error message: $!" );
1403   }
1404   my @output_lines = split( /\n/, $stdout ); # I'm splitting the output into an array for easier handling.
1405   debug -> warn( level => 2,
1406                  message => "$stdout" );
1407   foreach my $line( @output_lines ){ # loop over lines
1408     if( $line =~ /Job State:\s+Completed/ ){ # regexp to find finished jobs.
1409       debug -> warn( level => 1,
1410                      message => "Returning $jobId" );
1411       return $jobId; # Return the jobId found.
1412     }
1413   }
1414
1415   return 0;
1416
1417 # {{{ ud_native_retrieve
1418
1419 #    } else { # ud_native_retrieve
1420
1421 #     require SOAP::Lite;             # MGSI SOAP interface
1422 #     require LWP::UserAgent;         # MGSI Fileservice interface
1423
1424 #     my %confparams;
1425
1426 #     # MGSI constants
1427 #     uduserconf_read_file("./uduserconf", $ENV{HOME} ? ("$ENV{HOME}/uduserconf", "$ENV{HOME}/.uduserconf") : ());
1428 #     my $mgsisoapurl = $confparams{MGSI_SOAP_URL};
1429 #     my $mgsifilesvr = $confparams{MGSI_FILESVR_URL};
1430 #     my $mgsiuser    = $confparams{MGSI_USERNAME};
1431 #     my $mgsipwd     = $confparams{MGSI_PASSWORD};
1432
1433 #     # MGSI objects
1434 #     my $server = new SOAP::Lite
1435 #       -> uri('urn://ud.com/mgsi')
1436 #       -> proxy($mgsisoapurl);
1437 #     my $ua = new LWP::UserAgent; # mgsi filesvr HTTP object
1438
1439 #     ##############################
1440 #     ##  LOG IN TO MGSI SERVER   ##
1441 #     ##############################
1442 #     my $auth = soapvalidate($server->login($mgsiuser, $mgsipwd)); # mgsi authentication token
1443
1444 #     ################################
1445 #     ##  RETRIEVE JOB INFORMATION  ##
1446 #     ################################
1447 #     my $job = soapvalidate($server->getJobById($auth, $jobId));
1448
1449 #     my $jobstep;
1450 #     my $workunits;
1451 #     do {
1452 #       my $jobsteps = soapvalidate($server->getJobStepsByJob($auth, $$job{job_gid}));
1453
1454 #       $jobstep = @$jobsteps[1]; # this job only has one jobstep.
1455 #       my $workunits_array = soapvalidate($server->getWorkunits($auth,
1456 #                                                              {job_step_gid_match => [$$jobstep{job_step_gid}]},
1457 #                                                              "", 0, -1)); # no ordering, start at record 0, give me all wus
1458 #       $workunits = $$workunits_array{records};
1459 #       my $output_file = $$job{description};
1460
1461 #       #print "job $$job{job_gid}; jobstep $$jobstep{job_step_gid}\n" if $mgsidebug;
1462 #       #print "jobstep state is $$jobstep{state_id}\n" if $mgsidebug;
1463
1464 #       my $jobstepstatus = soapvalidate($server->getJobStepStatus($auth, $$jobstep{job_step_gid}));
1465
1466 #       sleep($self -> {'ud_sleep'});
1467 #     } while( $$jobstep{state_id} != 3 );
1468
1469 #     # Retrieve all results by going through every workunit in this job
1470
1471 #     foreach my $workunit (@$workunits) {
1472 #       my $results = soapvalidate($server->getResults($auth, {
1473 #       workunit_gid_match => [$$workunit{workunit_gid}], # filter by workunit guid
1474 #       success_active => 1, success_value => 1 # only retrieve successful results
1475 #                                                    },
1476 #                                                    "", 0, 1)); # no ordering, start at record 0, give me 1 result
1477 #       # if you want to actually compare redundant results to see if there is a quorum
1478 #       # here would be a good place to retrieve all results and 'diff' them
1479 #       # for now, I just retrieve 1 results through the queue and go with that
1480 #       if (not exists $$results{records}) {
1481 #       'debug' -> die( message => 'Found a Workunit without any successful Results, aborting retrieval.' );
1482
1483 #       }
1484 #       my $result = $$results{records}[0];
1485 #       my $tempfile = "package.tar"; #mktemp("tmpresXXXXXX");
1486 #       # open(PACKAGE, ">package.tar");
1487 #       my $resulturl = "$mgsifilesvr?auth=$auth&hash=$$result{result_hash}";
1488 #       my $request = HTTP::Request->new('GET', $resulturl);
1489 #       my $response = $ua->request($request, $tempfile);
1490 #       if ($response->is_error() ) {
1491 #       'debug' -> die( message => "Couldn't retrieve result file, server returned ".$response->status_line );
1492 #       } elsif ($response->header('Content-Length') != -s $tempfile) {
1493 #       'debug' -> die( message => "Incomplete file returned from server (expected ".$response->header('Content-Length')." but received ".(-s $tempfile).")." );
1494 #       }
1495
1496 #       require Archive::Tar;
1497
1498 #       my $tar = Archive::Tar->new;
1499
1500 #       $tar->read('package.tar',0, {extract => 1});
1501
1502 #       if( $Config{osname} ne 'MSWin32' ){
1503 #       cp('psn.LST','psn.lst');
1504 #       unlink('psn.LST');
1505 #       }
1506
1507 #       # add data to total frequency list
1508 #     }
1509
1510 #     # Optional deletion of job
1511 #     #if ($deletejob eq '-d') {
1512 #     #print "now deleting job $jobId...";
1513 #     #soapvalidate($server->deleteJob($auth, $$job{job_gid}));
1514 #     #}
1515
1516 #     # helper subroutines
1517
1518 #     # read in configuration
1519 #     sub uduserconf_read_file {
1520 #       my @files = @_;
1521
1522 #       my $file = undef;
1523 #       foreach (@files) {
1524 #       if (-f($_)) {
1525 #         $file = $_;
1526 #         last;
1527 #       }
1528 #       }
1529 #       defined($file) or 'debug' -> die(message => "Could not find any of: " . join(' ', @files) );
1530 #       open(FH,$file) or 'debug' -> die(message => "Could not open $file: $!");
1531
1532 #       # Main parsing loop for the file's contents.
1533 #       while (<FH>) {
1534 #       if (/^\s*(\w+)\s*=\s*\"([^\"]*)\"/ or /^\s*(\w+)\s*=\s*(\S+)\s*/) {
1535 #         $confparams{uc($1)} = $2;
1536 #       }
1537 #       }
1538
1539 #       close(FH);
1540
1541 #       foreach ("MGSI_FILESVR_URL",
1542 #              "MGSI_SOAP_URL",
1543 #              "MGSI_USERNAME",
1544 #              "MGSI_PASSWORD")
1545 #       {
1546 #       if (!defined($confparams{$_})) {
1547 #         'debug' -> die (message => "$_ must be defined in $file" );
1548 #       }
1549 #       }
1550 #     }
1551
1552 #     # soap response validation routine
1553 #     sub soapvalidate {
1554 #       my ($soapresponse) = @_;
1555 #       if ($soapresponse->fault) {
1556 #       'debug' -> die(message => "fault: ", $soapresponse->faultcode, " ", $soapresponse->faultstring );
1557 #       } else {
1558 #       return $soapresponse->result;
1559 #       }
1560 #     }
1561 #   }
1562
1563 # }}}
1564
1565 }
1566 end ud_monitor
1567
1568 # }}}
1569
1570 # {{{ ud_retrieve
1571
1572 start ud_retrieve
1573 {
1574   my $script;
1575   unless( defined $PsN::config -> {'_'} -> {'ud_nonmem'} ){
1576     if( $Config{osname} eq 'MSWin32' ) {
1577       $script = 'nonmem.bat';
1578     } else {
1579       $script = 'nonmem.sh';
1580     }
1581   } else {
1582     $script = $PsN::config -> {'_'} -> {'ud_nonmem'};
1583   }
1584
1585   my $subDir = "NM_run".($run_no+1);
1586   my ($tmp_dir, $file) = OSspecific::absolute_path( $self -> {'directory'}.'/' .
1587                                                     $subDir, '');
1588   if( system("$script -b -c -d ".$tmp_dir." -r $jobId > nonmem_bat_stdout") ){
1589     'debug' -> die( message => "UD submit script failed.\nSystem error message:$!" );
1590   }
1591
1592   if( $Config{osname} ne 'MSWin32' ){
1593     cp($tmp_dir.'/psn.LST',$tmp_dir.'/psn.lst');
1594     unlink($tmp_dir.'/psn.LST');
1595   }
1596 }
1597 end ud_retrieve
1598
1599 # }}}
1600
1601 # {{{ sge_submit
1602
1603 start sge_submit
1604 {
1605   my $fsubs = join( ',' , @{$model -> subroutine_files} );
1606   if( system( 'qsub -cwd -b y ' .
1607               ($self -> {'sge_resource'} ? '-l '.$self -> {'sge_resource'}.' ' : ' ') .
1608               ($self -> {'sge_queue'} ? '-q '.$self -> {'sge_queue'}.' ' : ' ') .
1609               ($PsN::config -> {'_'} -> {'remote_perl'} ? ' ' . $PsN::config -> {'_'} -> {'remote_perl'} : ' perl ') . " -I" .
1610               $PsN::lib_dir ."/../ " .
1611               $PsN::lib_dir . "/nonmem.pm" .
1612               " psn.mod psn.lst " .
1613               $self -> {'nice'} . " ".
1614               $nm_version . " " .
1615               1 . " " . # compilation
1616               1 . " " . # execution
1617               $fsubs . " " .
1618               $self -> {'nm_directory'} . ' > JobId' ) ){
1619     'debug' -> die( message => "Grid submit failed.\nSystem error message: $!" );
1620   }
1621
1622   open(JOBFILE, "JobId") or 'debug' -> die( message => "Couldn't open grid JobId file for reading: $!" );
1623   while( <JOBFILE> ){
1624     if( /Your job (\d+)/ ){
1625       $jobId = $1;
1626     }
1627   }
1628   close(JOBFILE);
1629 }
1630 end sge_submit
1631
1632 # }}}
1633
1634 # {{{ sge_monitor
1635
1636 start sge_monitor
1637 {
1638   my ($stdout, $stderr);
1639   run3("qstat -j $jobId 2> JobStat",undef,\$stdout, \$stderr );
1640   open(JOBFILE, "JobStat") or 'debug' -> die( message => "Couldn't open grid JobStat file for reading: $!" );
1641   while( <JOBFILE> ){
1642     if( /Following jobs do not exist:/ ){ # regexp to find finished jobs.
1643       close(JOBFILE);
1644       unlink( "JobStat" );
1645       return $jobId; # Return the jobId found.
1646     }
1647   }
1648   close(JOBFILE);
1649   unlink( "JobStat" );
1650
1651   return 0;
1652 }
1653 end sge_monitor
1654
1655 # }}}
1656
1657 # {{{ run_nonmem
1658
1659 start run_nonmem
1660     {
1661       my $candidate_model = $queue_info -> {'candidate_model'};
1662       my $tries = $queue_info -> {'tries'};
1663       # my $compile_only = $queue_info -> {'compile_only'};
1664       my $model = $queue_info -> {'model'};
1665
1666       # We do not expect any values of rerun lower than 1 here. (a bug otherwise...)
1667       if( not -e 'psn-' . ( $tries + 1 ) . '.lst' or $self -> {'rerun'} >= 2 ){
1668
1669         # {{{ Execution step
1670
1671         if( $self -> {'run_local'} ) {
1672
1673           # Normal local execution
1674
1675           my $fsubs = join( ',' , @{$candidate_model -> subroutine_files} );
1676
1677           my $command_line_options = " -I" .
1678               $PsN::lib_dir ."/../ " .
1679               $PsN::lib_dir . "/nonmem.pm" .
1680               " psn.mod psn.lst " .
1681               $self -> {'nice'} . " ".
1682               $nm_version . " " .
1683               1 . " " . # compilation
1684               1 . " " . # execution
1685               $fsubs . " " .
1686               $self -> {'nm_directory'} ;
1687
1688           if( $Config{osname} eq 'MSWin32' ){
1689
1690             # {{{ Windows execution
1691
1692             my $perl_bin = ($PsN::config -> {'_'} -> {'perl'} ? $PsN::config -> {'_'} -> {'perl'} : 'C:\Perl\bin\perl.exe');
1693
1694             require Win32::Process;
1695             require Win32;
1696             sub ErrorReport{ print Win32::FormatMessage( Win32::GetLastError() ); }
1697             my $proc;
1698             Win32::Process::Create($proc,$perl_bin,$perl_bin . $command_line_options,0,$Win32::Process::NORMAL_PRIORITY_CLASS,'.') || die ErrorReport();
1699
1700             $queue_info->{'winproc'}=$proc;
1701             $queue_map->{$proc->GetProcessID()} = $run_no;
1702
1703             # }}}
1704
1705           } else { #Asume *nix
1706
1707             # {{{ Unix execution
1708
1709             my $perl_bin = ($PsN::config -> {'_'} -> {'perl'} ? $PsN::config -> {'_'} -> {'perl'} : ' perl ');
1710
1711             my $pid = fork();
1712             if( $pid == 0 ){
1713               exec( $perl_bin . $command_line_options );
1714               exit; # Die Here if exec failed. Probably happens very rarely.
1715             }
1716             $queue_map->{$pid} = $run_no;
1717
1718             # }}}
1719
1720           }
1721
1722           $queue_info->{'start_time'} = time;
1723
1724         } elsif( $self -> {'run_on_lsf'} ) {
1725
1726           # lsf_submit will call the "nonmem" module that will figure
1727           # out that we want to run remotely. If we are also compiling
1728           # remotely, it will be done from here as well.
1729
1730           my $jobId = $self -> lsf_submit( model => $candidate_model,
1731                                            nm_version => $nm_version,
1732                                            work_dir   => $self -> {'directory'} . "/NM_run$run_no/");
1733
1734           $queue_map->{$jobId} = $run_no;
1735
1736         } elsif( $self -> {'run_on_umbrella'} ) {
1737
1738           # umbrella_submit will submit a request for execution
1739           # using the umbrella system.
1740           #sleep(1);
1741           $self -> umbrella_submit( queue_info => $queue_info,
1742                                     model      => $candidate_model,
1743                                     nm_version => $nm_version,
1744                                     prepare_jobs => 1);
1745         } elsif ( $self -> run_on_ud() ) {
1746           debug -> warn( level   => 1,
1747                          message => "Submitting to the UD system" );
1748           my $jobId = $self -> ud_submit( model => $candidate_model );
1749
1750           #$queue_info->{'compile_only'} = 0;
1751           $queue_map->{$jobId} = $run_no;
1752
1753         } elsif ( $self -> run_on_sge() ) {
1754           my $jobId = $self -> sge_submit( model => $candidate_model,
1755                                            nm_version => $nm_version );
1756
1757           #$queue_info->{'compile_only'} = 0;
1758           $queue_map->{$jobId} = $run_no;
1759         }
1760
1761         # }}}
1762
1763       } elsif( $self -> {'rerun'} >= 1 ){
1764
1765         # We are not forcing a rerun, but we want to recheck the
1766         # output files for errors. Therefore we put a fake entry in
1767         # queue_map to trigger "restart_nonmem()". We also need to
1768         # copy psn-x.lst to psn.lst to make sure that we don't stop
1769         # the next time we enter run_nonmem"
1770
1771         cp( 'psn-'.($tries+1).'.lst','psn.lst' );
1772         if( defined $model -> table_names ){
1773           foreach my $table_files( @{$model -> table_names} ){
1774             foreach my $table_file( @{$table_files} ){
1775               cp( $table_file.'-'.( $tries+1 ), $table_file );
1776             }
1777           }
1778         }
1779         if( defined $model -> extra_output() ){
1780           foreach my $file( @{$model -> extra_output()} ){
1781             cp( $file.'-'.( $tries+1 ), $file );
1782           }
1783         }
1784
1785         # $queue_info->{'compile_only'} = 0;
1786         $queue_map->{'rerun_'.$run_no} = $run_no; #Fake pid
1787       } # end of "not -e psn-$tries.lst or rerun"
1788
1789     }
1790 end run_nonmem
1791
1792 # }}}
1793
1794 # {{{ restart_needed
1795
1796 start restart_needed
1797 {
1798
1799   # -------------- Notes about automatic pertubation and retries -----------------
1800
1801   # Automatic pertubation of initial estimates are useful for two
1802   # purposes. One reason is when nonmem failes to produce a successful
1803   # minimization. In this case, we can try to direct the search by
1804   # selecting other estimates. It is also possible to get a successful
1805   # minimization in a local minima. In this case, we have no way of
1806   # knowing that it is a local minima without trying other initial
1807   # estimates. Two modelfit members govern the pertubation process;
1808   # "retries" which is a number giving the maximum number of retries
1809   # when nonmem failes and "min_retries" which is the number of runs
1810   # we want to do to get a global minima. If min_retries is 2 and
1811   # retries is 5 we will stop after 3 runs if we have reached a
1812   # successful minimization but continue until 5 if necessary.
1813
1814   # It is important to set $marked_for_rerun to 1 if $tries is
1815   # incremented!  Otherwise $tries can be incremented twice for
1816   # one run. The opposite is not true however, for instance a reset
1817   # of maxevals is not a retry but sets $marked_for_rerun to 1.
1818
1819   # We need the trail of files to select the most appropriate at the end
1820   # (see copy_model_and_output)
1821
1822   unless( defined $parm{'queue_info'} ){
1823     # The queue_info must be defined here!
1824     'debug' -> die( message => "Internal run queue corrupt\n" );
1825   }
1826   my $queue_info_ref = $parm{'queue_info'};
1827   my $run_results = $queue_info_ref -> {'run_results'};
1828   my $tries = \$queue_info_ref -> {'tries'};
1829   my $model = $queue_info_ref -> {'model'};
1830   my $candidate_model = $queue_info_ref -> {'candidate_model'};
1831   my $modelfile_tainted = \$queue_info_ref -> {'modelfile_tainted'};
1832
1833   my $lstsuccess = 0;
1834   for( my $lsttry = 1; $lsttry <= 5; $lsttry++ ){
1835     if( -e 'psn.lst' ){
1836
1837       my ( $output_file );
1838
1839       if( not -e 'psn-' . ( ${$tries} + 1 ) . '.lst' or $self -> {'rerun'} >= 2 ){
1840
1841         $output_file = $candidate_model -> outputs -> [0];
1842         $output_file -> abort_on_fail($self -> {'abort_on_fail'});
1843         $output_file -> _read_problems;
1844
1845         foreach my $filename ( @{$candidate_model -> output_files}, 'psn.mod', 'compilation_output.txt' ){
1846
1847           my $new_name = $self -> get_retry_name( filename => $filename,
1848                                                   retry => ${$tries} );
1849           mv( $filename, $new_name )
1850
1851         }
1852
1853       } else {
1854
1855         # This is rerun==1, i.e. re-evaluate the stuff that has been
1856         # run and (possibly) run extra runs to fix any problems left.
1857         # In this "else" statement it is true that psn-X.lst exists
1858         # and we copy it to psn.lst to make it the current version.
1859         cp( 'psn-'.(${$tries}+1).'.lst', 'psn.lst' );
1860
1861         $output_file = $candidate_model -> outputs -> [0];
1862         $output_file -> abort_on_fail($self -> {'abort_on_fail'});
1863         $output_file -> _read_problems;
1864
1865       }
1866
1867       # {{{ Create and write intermediate raw results
1868
1869       open( INTERMED, '>>intermediate_results.csv' );
1870       open( INTERMEDNONP, '>>intermediate_nonp_results.csv' );
1871
1872       my ($raw_results_row, $nonp_row) = $self -> create_raw_results_rows( max_hash => $self -> {'max_hash'},
1873                                                                            model => $candidate_model,
1874                                                                            model_number => $run_no + 1,
1875                                                                            raw_line_structure => $self -> {'raw_line_structure'} );
1876
1877       $queue_info_ref -> {'raw_results'} -> [${$tries}] = $raw_results_row;
1878       $queue_info_ref -> {'raw_nonp_results'} -> [${$tries}] = $nonp_row;
1879
1880       foreach my $row ( @{$raw_results_row} ){
1881         print INTERMED join( ',', @{$row} ), "\n";
1882       }
1883
1884       foreach my $row ( @{$nonp_row} ){
1885         print INTERMEDNONP join( ',', @{$nonp_row} ), "\n";
1886       }
1887
1888       close( INTERMED );
1889       close( INTERMEDNONP );
1890
1891       # }}}
1892
1893       # {{{ Check for minimization successfull an try to find out if lst file is truncated
1894
1895       my ( $minimization_successful, $minimization_message );
1896
1897       if( $output_file -> parsed_successfully() and
1898           not defined $output_file -> problems ){
1899         # This should not happen if we are able to parse the output file correctly
1900         $run_results -> [${$tries}] -> {'failed'} = 1;
1901         return(0);
1902       }
1903
1904       if( not $output_file -> parsed_successfully() ){
1905
1906         if( $self -> {'handle_crashes'} and $queue_info_ref -> {'crashes'} < $self -> crash_restarts() ) {
1907
1908           # If the output file could not be parsed successfully, this is
1909           # a sign of a crashed run. This is not a NONMEM error as such
1910           # but probably an operating system problem. To handle this, we
1911           # mark this for rerunning but do not increase the $tries
1912           # variable but instead increase $crashes and check whether
1913           # this value is below or equal to $crash_restarts.
1914           debug -> warn( level => 1,
1915                          message => "Restarting crashed run ".
1916                          $output_file -> full_name().
1917                          "\n".$output_file -> parsing_error_message() );
1918
1919
1920           $queue_info_ref -> {'crashes'}++;
1921
1922           my $message = "\nModel in NM_run".($run_no+1)." crashed, try nr ". ($queue_info_ref -> {'crashes'} );
1923           ui -> print( category => 'all',  message  => $message,
1924                        newline => 1);
1925
1926
1927           if( $self -> {'handle_msfo'} ){
1928             $self -> set_msfo_to_msfi( candidate_model => $candidate_model,
1929                                        retry => ${$tries},
1930                                        queue_info => $queue_info_ref);
1931           } else {
1932             cp( 'psn-'.(${$tries}+1).'.mod', 'psn.mod' );
1933           }
1934
1935           mv( 'psn-'.(${$tries}+1).'.lst', 'psn_crash-'.(${$tries}+1).'-'.$queue_info_ref -> {'crashes'}.'.lst' );
1936
1937           $output_file -> flush;
1938
1939           return(1); # Return a one (1) to make run() rerun the
1940                      # model. By returning here, we avoid the
1941                      # perturbation of the initial estimates later on in
1942                      # this method.
1943         } else {
1944           my $message = "\nModel in NM_run".($run_no+1)." crashed ".(($queue_info_ref -> {'crashes'}+1)." times. Not restarting." );
1945           ui -> print( category => 'all',  message  => $message,
1946                        newline => 1);
1947
1948           $output_file -> flush;
1949
1950           return(0);
1951         }
1952       }
1953
1954       # If the output file was parsed successfully, we (re)set the $crashes
1955       # variable and continue
1956
1957       $queue_info_ref -> {'crashes'} = 0;
1958
1959       $minimization_successful = $output_file -> minimization_successful();
1960       $minimization_message    = $output_file -> minimization_message();
1961
1962       unless( defined $minimization_successful ) {
1963         debug -> die( message => "No minimization status found in " . $output_file ->filename );
1964       }
1965
1966       # {{{ log the stats of this run
1967
1968       foreach my $category ( 'minimization_successful', 'covariance_step_successful',
1969                              'covariance_step_warnings', 'estimate_near_boundary',
1970                              'significant_digits', 'ofv' ){
1971         my $res = $output_file -> $category;
1972         $run_results -> [${$tries}] -> {$category} = defined $res ? $res -> [0][0] : undef;
1973       }
1974       $run_results -> [${$tries}] -> {'pass_picky'} = 0;
1975
1976       # }}}
1977
1978       # }}}
1979
1980       # {{{ Check if maxevals is reached and copy msfo to msfi
1981
1982       if ( not $marked_for_rerun and $handle_maxevals ) {
1983
1984         for ( @{$minimization_message -> [0][0]} ) {
1985           if ( /\s*MAX. NO. OF FUNCTION EVALUATIONS EXCEEDED\s*/) {
1986
1987             my $queue_info_ref -> {'evals'} += $output_file -> feval -> [0][0];
1988             my $maxeval = $model -> maxeval -> [0][0];
1989
1990             if( $maxeval > $queue_info_ref -> {'evals'} ){
1991               $self -> set_msfo_to_msfi( candidate_model => $candidate_model );
1992
1993               $candidate_model -> _write;
1994               ${$modelfile_tainted} = 1;
1995               $marked_for_rerun = 1;
1996               last;
1997             }
1998           }
1999         }
2000       }
2001
2002       # }}}
2003
2004       # {{{ Check for rounding errors and/or hessian_npd messages
2005
2006       if ( not $marked_for_rerun and $handle_rounding_errors || $handle_hessian_npd) {
2007         my $round_rerun;
2008         my $hessian_rerun;
2009         for ( @{$minimization_message -> [0][0]} ) {
2010           $round_rerun = 1 and last if ( /\s*DUE TO ROUNDING ERRORS\s*/);
2011           $hessian_rerun = 1 and last if ( /\s*NUMERICAL HESSIAN OF OBJ. FUNC. FOR COMPUTING CONDITIONAL ESTIMATE IS NON POSITIVE DEFINITE\s*/);
2012         }
2013
2014         if ( ($round_rerun && $handle_rounding_errors) or ($hessian_rerun && $handle_hessian_npd)) {
2015
2016           if( $self -> {'use_implicit_msfo'} ) {
2017             $self -> reset_msfo( basic_model => $model,
2018                                  candidate_model => $candidate_model );
2019           }
2020
2021           $self -> cut_thetas( candidate_model => $candidate_model,
2022                                cutoff_thetas => \@cutoff_thetas,
2023                                output_file => $output_file );
2024
2025           $candidate_model -> _write;
2026           ${$modelfile_tainted} = 1;
2027           #print "ROUND\n";
2028           $marked_for_rerun = 1;
2029
2030           ${$tries} ++; # This is a precaution, handle_rounding and handle_hessian should have
2031           # their own termination checks
2032         }
2033
2034       }
2035
2036       # }}}
2037
2038       # {{{ Check for failed problems and possibly check for picky errors.
2039
2040       if ( not $marked_for_rerun and $tweak_inits ) {
2041
2042         my @reruns;
2043         my @problems = @{$candidate_model -> problems};
2044         for ( my $problem = 1; $problem <= scalar @problems; $problem++ ) {
2045           unless( $candidate_model -> is_simulation( problem_number => $problem ) ){
2046             if ( $minimization_successful -> [$problem-1][0] ) {
2047               if ( $picky ) {
2048                 $run_results -> [${$tries}] -> {'pass_picky'} = 1;
2049                 for ( @{$minimization_message -> [$problem-1][0]} ) {
2050                   if ( /0COVARIANCE STEP ABORTED/ or
2051                        /0PROGRAM TERMINATED BY OBJ/ or
2052                        /0ESTIMATE OF THETA IS NEAR THE BOUNDARY AND/ or
2053                        /0PARAMETER ESTIMATE IS NEAR ITS BOUNDARY/ or
2054                        /0R MATRIX ALGORITHMICALLY SINGULAR/ or
2055                        /0S MATRIX ALGORITHMICALLY SINGULAR/ ) {
2056                     push( @reruns, $problem );
2057                     $run_results -> [${$tries}] -> {'pass_picky'} = 0;
2058                     last;
2059                   }
2060                 }
2061               }
2062             } else {
2063               my $significant_digits = $output_file -> significant_digits;
2064               if ( not ( $significant_digits_rerun and $significant_digits -> [$problem-1][0] > $significant_digits_rerun ) ) {
2065                 push( @reruns, $problem );
2066               }
2067             }
2068           }
2069         }
2070
2071         if( ${$tries} < ($retries -1) and scalar @reruns > 0 ) {
2072           $marked_for_rerun = 1;
2073           ${$tries} ++;
2074
2075           if( ${$tries} >= $self -> {'min_retries'} and $self -> {'verbose'} ){
2076             my $message = "R:".($run_no+1).":". (${$tries}+1) . " ";
2077             ui -> print( category => 'all',  message  => $message,
2078                          newline => 0);
2079           }
2080
2081           if( $self -> {'ask_if_fail'} ) {
2082             $retries = $self -> ask_user( basic_model => $model,
2083                                           candidate_model => $candidate_model,
2084                                           reruns => \@reruns,
2085                                           minimization_message => $minimization_message,
2086                                           output_file => $output_file,
2087                                           retries => $retries,
2088                                           tries => ${$tries} );
2089             $candidate_model->_write;
2090             ${$modelfile_tainted} = 1;
2091           } else {
2092
2093             # This code must be adjusted for multiple problems!!
2094             my $degree = 0.1*${$tries};
2095             if( $self -> {'handle_msfo'} ) {
2096               print "Ost på räka\n";
2097               $self -> reset_msfo( basic_model => $model,
2098                                    candidate_model => $candidate_model );
2099
2100               foreach my $prob ( @reruns ) {
2101                 $problems[$prob-1] -> set_random_inits ( degree => $degree );
2102               }
2103
2104               $candidate_model->_write;
2105               ${$modelfile_tainted} = 1;
2106
2107             } else {
2108               foreach my $prob ( @reruns ) {
2109                 $problems[$prob-1] -> set_random_inits ( degree => $degree );
2110               }
2111
2112               $candidate_model->_write;
2113               ${$modelfile_tainted} = 1;
2114               # The fcon code does not parse the dofetilide model correctly
2115 #               my $fcon = fcon -> new( filename => 'FCON.orig' );
2116 #               $fcon -> parse;
2117 #               $fcon -> pertubate_all( fixed_thetas => $final_model -> fixed( 'parameter_type' => 'theta' ),
2118 #                                       fixed_omegas => $final_model -> fixed( 'parameter_type' => 'omega' ),
2119 #                                       fixed_sigmas => $final_model -> fixed( 'parameter_type' => 'sigma' ),
2120 #                                       degree       => $degree );
2121 #               $fcon -> write( filename => 'FCON' );
2122             }
2123           }
2124         }
2125       }
2126
2127       # }}}
2128
2129       # {{{ Perturb estimates if min_retries not reached
2130
2131       # This "if" block should conceptually be last, since it is
2132       # something we should only do if the model succeeds. In
2133       # practise it could swap places with at least the tweak inits
2134       # block, but for simplicities sake, lets leave it at the
2135       # bottom.
2136
2137       if( not $marked_for_rerun and ${$tries} < $self -> {'min_retries'} ) {
2138         #Here we force pertubation when the model is successful.
2139
2140         ${$tries} ++;
2141         $marked_for_rerun = 1;
2142         my $degree = 0.1 * ${$tries};
2143         if( $self -> {'use_implicit_msfo'} and
2144             $self -> reset_msfo( basic_model => $model,
2145                                  candidate_model => $candidate_model ) ){
2146
2147           foreach my $prob ( @{$candidate_model -> problems} ) {
2148             $prob -> set_random_inits ( degree => $degree );
2149           }
2150
2151           $candidate_model->_write;
2152
2153         } else {
2154           foreach my $prob ( @{$candidate_model -> problems} ) {
2155             $prob -> set_random_inits ( degree => $degree );
2156           }
2157
2158           $candidate_model->_write;
2159         }
2160
2161         ${$modelfile_tainted} = 1;
2162       }
2163
2164       # }}}
2165
2166       $output_file -> flush;
2167
2168       $lstsuccess = 1; # We did find the lst file.
2169       last;
2170     } else {
2171       sleep(($lsttry+1)**2);
2172       # print " The lst-file is not present, trying ".(5-$lsttry)." times more\n";
2173     } # Did the lst file exist?
2174   } # The loop trying to read the lst file
2175   unless( $lstsuccess ) { # psn.lst doesn't exist.
2176     $run_results -> [${$tries}] -> {'failed'} = 1;
2177   }
2178
2179 }
2180 end restart_needed
2181
2182 # }}}
2183
2184 # {{{ select_best_model
2185
2186 start select_best_model
2187 {
2188   # -------------- Notes about Final model selection -----------------
2189
2190   # Since we have reruns with pertubation and now also forced (or
2191   # automatic) pertubation the final model is not equal to the
2192   # original model. We consider four implicit subsets. Those that pass
2193   # the picky test, those that don't pass the picky test but have
2194   # minimization successful, those that don't pass the minimization
2195   # step but produce an ofv and, finaly, those that doesn't produce an
2196   # ofv. The final model will be the model that passes the most tests
2197   # and have the lowest ofv value, and if no ofv value is produced, it
2198   # will be the basic model.
2199
2200   # Get all runs that passed the picky test (if picky is used)
2201
2202   # The order of categories is important. Highest priority goes last.
2203
2204   unless( defined $parm{'queue_info'} ){
2205     # The queue_info must be defined here!
2206     'debug' -> die( message => "Internal run queue corrupt\n" );
2207   }
2208   my $queue_info_ref = $parm{'queue_info'};
2209   my $run_results = $queue_info_ref -> {'run_results'};
2210   my $model = $queue_info_ref -> {'model'};
2211   my $candidate_model = $queue_info_ref -> {'candidate_model'};
2212
2213   my @selection_categories = ('really_bad','terminated','normal','picky');
2214   my ( %selection, $selected );
2215   foreach my $category ( @selection_categories ) {
2216     $selection{$category}{'best_significant_digits'} = 0;
2217   }
2218
2219   # First pass to get lowest OFV's
2220   $selection{'normal'}{'lowest_OFV'} = 999999999;
2221   $selection{'terminated'}{'lowest_OFV'} = 999999999;
2222   $selection{'really_bad'}{'lowest_OFV'} = 999999999;
2223   for(my $i = 0; $i < scalar @{$run_results}; $i++ ){
2224     if ( $run_results -> [$i] -> {'minimization_successful'} ) {
2225       if( defined( $run_results -> [$i] -> {'ofv'} ) ) {
2226         $selection{'normal'}{'lowest_OFV'} = $run_results -> [$i] -> {'ofv'} < $selection{'normal'}{'lowest_OFV'} ?
2227             $run_results -> [$i] -> {'ofv'} : $selection{'normal'}{'lowest_OFV'};
2228       }
2229     }
2230     if( defined( $run_results -> [$i] -> {'ofv'} ) ) {
2231       if( defined( $run_results -> [$i] -> {'significant_digits'} ) ) {
2232         $selection{'terminated'}{'lowest_OFV'} = $run_results -> [$i] -> {'ofv'} < $selection{'terminated'}{'lowest_OFV'} ?
2233             $run_results -> [$i] -> {'ofv'} : $selection{'terminated'}{'lowest_OFV'};
2234       }
2235       $selection{'really_bad'}{'lowest_OFV'} = $run_results -> [$i] -> {'ofv'} < $selection{'really_bad'}{'lowest_OFV'} ?
2236           $run_results -> [$i] -> {'ofv'} : $selection{'really_bad'}{'lowest_OFV'};
2237     }
2238   }
2239
2240   my $accepted_OFV_diff = 5;
2241
2242   # Loop through all categories, the order is not important here
2243   for(my $i = 0; $i < scalar @{$run_results}; $i++ ){
2244     foreach my $category ( @selection_categories ) {
2245
2246       if( $category eq 'picky' ) {
2247
2248         if ( $run_results -> [$i] -> {'pass_picky'} ) {
2249           if( $run_results -> [$i] -> {'significant_digits'} >
2250               $selection{$category}{'best_significant_digits'} and
2251               $run_results -> [$i] -> {'ofv'} < ($selection{'normal'}{'lowest_OFV'}+$accepted_OFV_diff) ){
2252             $selection{$category}{'selected'} = ($i+1);
2253             $selection{$category}{'best_significant_digits'} = $run_results -> [$i] -> {'significant_digits'};
2254           }
2255         }
2256       } elsif( $category eq 'normal' ) {
2257
2258         if ( $run_results -> [$i] -> {'minimization_successful'} ) {
2259           if( $run_results -> [$i] -> {'significant_digits'} >
2260               $selection{$category}{'best_significant_digits'} and
2261               $run_results -> [$i] -> {'ofv'} < ($selection{'normal'}{'lowest_OFV'}+$accepted_OFV_diff) ){
2262             $selection{$category}{'selected'} = ($i+1);
2263             $selection{$category}{'best_significant_digits'} = $run_results -> [$i] -> {'significant_digits'};
2264           }
2265         }
2266       } elsif( $category eq 'terminated' ) {
2267         if ( defined( $run_results -> [$i] -> {'ofv'} ) and
2268              $run_results -> [$i] -> {'ofv'} < ($selection{'terminated'}{'lowest_OFV'}+$accepted_OFV_diff) ) {
2269           if( defined( $run_results -> [$i] -> {'significant_digits'} ) ) {
2270             if ( $run_results -> [$i] -> {'significant_digits'} >
2271                  $selection{$category}{'best_significant_digits'} ){
2272               $selection{$category}{'selected'} = ($i+1);
2273               $selection{$category}{'best_significant_digits'} = $run_results -> [$i] -> {'significant_digits'};
2274             }
2275           }
2276         }
2277       } else {
2278         if ( defined( $run_results -> [$i] -> {'ofv'} ) and
2279              $run_results -> [$i] -> {'ofv'} < ($selection{'really_bad'}{'lowest_OFV'}+$accepted_OFV_diff) ) {
2280           $selection{$category}{'selected'} = ($i+1);
2281           $selection{$category}{'best_significant_digits'} = $run_results -> [$i] -> {'significant_digits'};
2282         }
2283       }
2284     }
2285   }
2286
2287   # Loop through all categories from less strict to strict and
2288   # replace the selected run as we find better runs. (I know that
2289   # this is a bit awkward but it is working.)
2290   foreach my $category ( @selection_categories ) {
2291     $selected = defined $selection{$category}{'selected'} ?
2292         $selection{$category}{'selected'} : $selected;
2293   }
2294   $selected = defined $selected ? $selected : 1;
2295
2296   open( STAT, '>stats-runs.csv' );
2297   print STAT Dumper \@{$run_results};
2298   print STAT "Selected $selected\n";
2299   close( STAT );
2300
2301   unless( $run_results -> [$selected-1] -> {'failed'} ){
2302
2303     my @raw_results_rows = @{$queue_info_ref -> {'raw_results'} -> [$selected-1]};
2304
2305     foreach my $row ( @raw_results_rows ){
2306       shift( @{$row} );
2307       unshift( @{$row}, $run_no+1 );
2308     }
2309
2310
2311     push( @{$self -> {'raw_results'}}, @raw_results_rows );
2312     push( @{$self -> {'raw_nonp_results'}}, @{$queue_info_ref -> {'raw_nonp_results'} -> [$selected-1]} );
2313
2314
2315     $self -> copy_model_and_output( final_model   => $candidate_model,
2316                                     model         => $model,
2317                                     use_run       => $selected ? $selected : '' );
2318
2319   }
2320
2321   if( $self -> {'nonparametric_etas'} and
2322       ( not -e 'np_etas.lst' or $self -> {'rerun'} >=2 )  ) {
2323
2324     # ---------------------  Create nonp eta model  -----------------------------
2325
2326     # {{{ nonp eta model
2327
2328     if( not -e 'np_etas.lst' or $self -> {'rerun'} >= 2 ){
2329
2330       ui -> print( category => 'execute',
2331                    message  => 'Creating NPETA model' );
2332
2333       my $np_eta_mod = $candidate_model ->
2334           copy( filename    => $self -> {'directory'}.
2335                 'NM_run'.($run_no+1).'/np_etas.mod',
2336                 target      => 'mem',
2337                 copy_data   => 0,
2338                 copy_output => 0);
2339
2340       my ( $msfo_ref, $junk ) = $candidate_model ->
2341           _get_option_val_pos( name            => 'MSFO',
2342                                record_name     => 'estimation' );
2343       # We should have an MSFO file here
2344       for( my $i = 0; $i < scalar @{$msfo_ref}; $i++ ) {
2345         my $msfi = $msfo_ref->[$i][0];
2346         $np_eta_mod -> set_records( problem_numbers => [($i+1)],
2347                                     type            =>'msfi',
2348                                     record_strings  => ["$msfi"]);
2349         $np_eta_mod -> set_records( problem_numbers => [($i+1)],
2350                                     type           => 'nonparametric',
2351                                     record_strings => [ 'ETAS UNCONDITIONAL '.
2352                                                         'MSFO=npmsfo'.($i+1) ] );
2353       }
2354       $np_eta_mod -> remove_option( record_name => 'estimation',
2355                                     option_name => 'MSFO' );
2356       $np_eta_mod -> remove_records( type => 'theta' );
2357       $np_eta_mod -> remove_records( type => 'omega' );
2358       $np_eta_mod -> remove_records( type => 'sigma' );
2359       $np_eta_mod -> remove_records( type => 'table' );
2360       my @nomegas = @{$candidate_model -> nomegas};
2361
2362       my @max_evals;
2363       for( my $i = 0; $i <= $#nomegas; $i++ ) {
2364         my $marg_str = 'ID';
2365         for( my $j = 1; $j <= $nomegas[$i]; $j++ ) {
2366           $marg_str = $marg_str.' ETA'.$j;
2367         }
2368         $marg_str = $marg_str.' FILE='.$model -> filename.'.nonp_etas'.
2369             ' NOAPPEND ONEHEADER NOPRINT FIRSTONLY';
2370         $np_eta_mod -> add_records( problem_numbers => [($i+1)],
2371                                     type            => 'table',
2372                                     record_strings  => [ $marg_str ] );
2373         $max_evals[$i] = [0];
2374       }
2375
2376 #       my @table_names = @{$np_eta_mod -> table_names};
2377 #       for ( my $i = 0; $i <= $#table_names; $i++ ) {
2378 #         foreach my $table ( @{$table_names[$i]} ) {
2379 #           $table = 'nonp_'.$table;
2380 #         }
2381 #       }
2382 #       $np_eta_mod -> table_names( new_names => \@table_names );
2383
2384       $np_eta_mod -> maxeval( new_values => \@max_evals );
2385       $np_eta_mod -> _write;
2386
2387       # }}} nonp eta model
2388
2389       # ----------------------  run nonp eta model  -------------------------------
2390
2391       # {{{ run eta model
2392
2393       ui -> print( category => 'execute',
2394                    message  => 'Running NPETA model' );
2395
2396       my $nonmem_object = nonmem -> new( adaptive => $self -> {'adaptive'},
2397                                       modelfile => 'np_etas.mod',
2398                                       version => $nm_version,
2399                                       nm_directory => $self -> {'nm_directory'},
2400                                       nice => $self -> {'nice'},
2401                                       show_version => not $self -> {'run_on_lsf'});
2402
2403       $nonmem_object -> fsubs( $np_eta_mod -> subroutine_files );
2404       unless( $nonmem_object -> compile() ){
2405         debug -> die( message => "NONMEM compilation failed:\n" .$nonmem_object -> error_message );
2406       }
2407       if( $nonmem_object -> nmtran_message =~ /WARNING/s and $self -> {'verbose'} ){
2408         ui -> print(category => 'all',
2409                     message => "NONMEM Warning: " . $nonmem_object -> nmtran_message );
2410       }
2411       $nonmem_object -> execute();
2412
2413       foreach my $table_files( @{$np_eta_mod -> table_names} ){
2414         foreach my $table_file( @{$table_files} ){
2415           my $dir = $model -> directory;
2416           cp( $table_file, $dir  );
2417         }
2418       }
2419
2420       unlink 'nonmem', 'nonmem6', 'nonmem5','nonmem.exe', 'nonmem6_adaptive', 'nonmem5_adaptive';
2421     }
2422
2423     # }}} run eta model
2424
2425   }
2426 }
2427 end select_best_model
2428
2429 # }}}
2430
2431 # {{{ print_finish_message
2432
2433 start print_finish_message
2434     {
2435       my $ui_text;
2436       # Log the run
2437       $ui_text .= sprintf("%3s",$run+1) . sprintf("%25s",$self -> {'models'} -> [$run] -> filename);
2438       my $log_text = $run+1 . ',' . $self -> {'models'} -> [$run] -> filename . ',';
2439       if( $self -> {'verbose'} or $self -> {'quick_summarize'} ){
2440         foreach my $param ( 'ofv', 'covariance_step_successful', 'minimization_message' ) {
2441           if( $param eq 'minimization_message' ){
2442             $ui_text .= "\n    ---------- Minimization Message ----------\n";
2443           }
2444           if( defined $candidate_model ){
2445             my $ests = $candidate_model -> outputs -> [0] -> $param;
2446             # Loop the problems
2447             for ( my $j = 0; $j < scalar @{$ests}; $j++ ) {
2448               if ( ref( $ests -> [$j][0] ) ne 'ARRAY' ) {
2449                 $ests -> [$j][0] =~ s/^\s*//;
2450                 $ests -> [$j][0] =~ s/\s*$//;
2451                 $log_text .= $ests -> [$j][0] .',';
2452                 #chomp($ests -> [$j][0]);
2453                 $ui_text .= sprintf("%10s",$ests -> [$j][0]);
2454               } else {
2455
2456                 # Loop the parameter numbers (skip sub problem level)
2457                 for ( my $num = 0; $num < scalar @{$ests -> [$j][0]}; $num++ ) {
2458                   #$ests -> [$j][0][$num] =~ s/^\s*//;
2459                   #$ests -> [$j][0][$num] =~ s/\s*$/\n/;
2460                   $log_text .= $ests -> [$j][0][$num] .',';
2461                   #chomp($ests -> [$j][0][$num]);
2462                   if( $param eq 'minimization_message' ){
2463                     $ui_text .= "    ";
2464                   }
2465                   $ui_text .= sprintf("%12s",$ests -> [$j][0][$num]);
2466                 }
2467               }
2468             }
2469           }
2470           if( $param eq 'minimization_message' ){
2471             $ui_text .= "    ------------------------------------------\n\n";
2472           }
2473         }
2474         ui -> print( category => 'all',
2475                      message  => $ui_text,
2476                      wrap     => 0,
2477                      newline => 0);
2478       }
2479
2480       open( LOG, ">>".$self -> {'logfile'} );
2481       print LOG $log_text;
2482       print LOG "\n";
2483       close LOG;
2484
2485     }
2486 end print_finish_message
2487
2488 # }}}
2489
2490 # {{{ run
2491
2492 start run
2493     {
2494
2495       my $cwd = getcwd();
2496       chdir( $self -> {'directory'} );
2497
2498       # {{{ sanity checks
2499
2500       my @models;
2501       if ( defined $self -> {'models'} ) {
2502         @models = @{ $self -> {'models'} };
2503       } else {
2504         debug -> die( message => "Have no models!" );
2505       }
2506
2507       my $threads = $self -> {'threads'};
2508       $threads = $#models + 1 if ( $threads > $#models + 1);
2509
2510       # Unless we store "run_no" in the database umbrella can't be parallel.
2511
2512       $threads = 1 if $self -> {'run_on_umbrella'};
2513
2514       # }}}
2515
2516       # {{{ print starting messages
2517
2518       ui -> print( category => 'all',
2519                    message  => 'Starting ' . scalar(@models) . ' NONMEM executions. '. $threads .' in parallel.' )
2520           unless $self -> {'parent_threads'} > 1;
2521       ui -> print( category => 'all',
2522                    message  => "Run number\tModel name\tOFV\tCovariance step successful." )  if $self -> {'verbose'};
2523
2524       # }}}
2525
2526       # {{{ Print model-NM_run translation file
2527
2528       open( MNT, ">model_NMrun_translation.txt");
2529       for ( my $run = 0; $run <= $#models; $run ++ ) {
2530         print MNT sprintf("%-40s",$models[$run]->filename),"NM_run",($run+1),"\n";
2531       }
2532       close( MNT );
2533
2534       # }}}
2535
2536       # {{{ Setup of moshog
2537
2538       my $moshog_client;
2539
2540       if( $self -> {'adaptive'} and $self -> {'threads'} > 1 ) {
2541
2542         # Initiate the moshog client
2543         $moshog_client = moshog_client -> new(start_threads => $self -> {'threads'});
2544
2545       }
2546       # }}}
2547
2548       # {{{ Local execution
2549
2550       # %queue_map is a mapping from nonmem.pm pID to run number.
2551
2552       my %queue_map;
2553
2554       # %queue_info is keyed on run number and contains information
2555       # about each nonmem run.
2556
2557       my %queue_info;
2558
2559       # @queue is an array of NM_run directory numbers. If "X" is in
2560       # @queue, then psn.mod in "NM_runX" is scheduled to be run. It
2561       # initialized to all models in the tool. Note that if X is in
2562       # the queue, it doesn't mean NM_runX exists.
2563
2564       my @queue = (0..$#models);
2565       my $all_jobs_started = 0;
2566
2567       # We loop while there is content in the queue (which shrinks and grows)
2568       # and while we have jobs running (represented in the queue_info)
2569
2570       my $poll_count = 0;
2571       while( @queue or (scalar keys %queue_map > 0) ){
2572         $poll_count++;
2573         # We start by looking for jobs that have been started and
2574         # finished. If we find one, we set $pid to that job.
2575
2576         my $pid = 0;
2577
2578         # {{{ Get potiential finished pid
2579
2580         foreach my $check_pid( keys %queue_map ){
2581
2582           if( $check_pid =~ /^rerun_/ ){
2583
2584             # A pid that starts with "rerun" is a rerun and is always
2585             # "finished".
2586
2587             $pid = $check_pid;
2588             last;
2589           }
2590
2591           # Diffrent environments have diffrent ways of reporting
2592           # job status. Here we check which environment we are in
2593           # and act accordingly.
2594
2595           if( $self -> {'run_on_ud'} ){
2596
2597             $pid = $self -> ud_monitor( jobId => $check_pid );
2598
2599             if( $pid ){
2600               $self -> ud_retrieve( jobId => $check_pid,
2601                                     run_no => $queue_map{$check_pid} );
2602             }
2603
2604           } elsif( $self -> {'run_on_sge'} ) {
2605
2606             $pid = $self -> sge_monitor( jobId => $check_pid );
2607
2608           } elsif( $self -> {'run_on_lsf'} ) {
2609
2610             $pid = $self -> lsf_monitor( jobId => $check_pid );
2611
2612           } else { # Local process poll
2613
2614             if( $Config{osname} eq 'MSWin32' ){
2615
2616               my $exit_code;
2617
2618               # GetExitCode is supposed to return a value indicating
2619               # if the process is still running, however it seems to
2620               # allways return 0. $exit_code however is update and
2621               # seems to be nonzero if the process is running.
2622
2623               $queue_info{$queue_map{$check_pid}}{'winproc'}->GetExitCode($exit_code);
2624
2625               if( $exit_code == 0 ){
2626                 $pid = $check_pid;
2627               }
2628
2629             } else {
2630
2631               $pid = waitpid($check_pid,WNOHANG);
2632
2633               # Waitpid will return $check_pid if that process has
2634               # finished and 0 if it is still running.
2635
2636               if( $pid == -1 ){
2637                 # If waitpid return -1 the child has probably finished
2638                 # and has been "Reaped" by someone else. We treat this
2639                 # case as the child has finished. If this happens many
2640                 # times in the same NM_runX directory, there is probably
2641                 # something wrong and we die(). (I [PP] suspect that we
2642                 # never/seldom reach 10 reruns in one NM_runX directory)
2643
2644                 my $run = $queue_map{$check_pid};
2645
2646                 $queue_info{$run}{'nr_wait'}++;
2647                 if( $queue_info{$run}{'nr_wait'} > 10 ){
2648                   debug -> die(message => "Nonmem run was lost\n");
2649                 }
2650                 $pid = $check_pid;
2651               }
2652
2653               # If the pid is not set, the job is not finished and we
2654               # check if it has been running longer than it should, and we
2655               # kill it.
2656               unless( $pid ){
2657                 if( $self -> {'max_runtime'}
2658                     and
2659                     time > $queue_info{$queue_map{$check_pid}}->{'start_time'} + $self -> {'max_runtime'}  ){
2660
2661                   # Only works on windows
2662
2663                   if( $Config{osname} ne 'MSWin32' ){
2664
2665                     print "Job ", $check_pid, " exceeded run time, trying to kill\n";
2666                     my $try=1;
2667                     while( kill( 0, $check_pid ) ){
2668                       if( $try <= 5 ){
2669                         kill( 'TERM', $check_pid );
2670                       } else {
2671                         kill( 'KILL', $check_pid );
2672                       }
2673                       waitpid( $check_pid, 0);
2674                       sleep(1);
2675                     }
2676                     print "Job ", $check_pid, " killed successfully on try nr $try\n";
2677                   }
2678                 }
2679               }
2680             }
2681           }
2682
2683           last if $pid;
2684
2685         }
2686
2687         # }}}
2688
2689         if( not $pid ){
2690
2691           # No process has finished.
2692
2693           # {{{ Return to polling if queue is empty or we have max number of jobs running.
2694
2695           if( (scalar @queue == 0) or scalar keys %queue_map >= $threads ){
2696
2697             # In that case we should not start another job. (Also
2698             # sleep to make polling less demanding).
2699
2700             if( defined $PsN::config -> {'_'} -> {'job_polling_interval'} and
2701                 $PsN::config -> {'_'} -> {'job_polling_interval'} > 0 ) {
2702               sleep($PsN::config -> {'_'} -> {'job_polling_interval'});
2703             } else {
2704               sleep(1);
2705             }
2706             my $make_request = 0;
2707             if ( defined $PsN::config -> {'_'} -> {'polls_per_moshog_request'} and
2708                  $self -> {'adaptive'} and
2709                  $self -> {'threads'} > 1 and
2710                  not ( $poll_count % $PsN::config -> {'_'} -> {'polls_per_moshog_request'} ) ) {
2711               $make_request = 1;
2712             }
2713
2714             if( $make_request and
2715                 scalar @queue > 0 ) {
2716
2717               # Make a moshog request.
2718
2719               $moshog_client -> request( request => scalar @queue );
2720               $threads += $moshog_client -> granted();
2721
2722             }
2723
2724             next; # Return to polling for finished jobs.
2725           }
2726
2727           # }}}
2728
2729           # This is where we initiate a new job:
2730
2731           my $run = shift( @queue );
2732
2733           # {{{ check for no run conditions. (e.g. job already run)
2734
2735           if ( -e $self -> {'models'} -> [$run] -> outputs -> [0] -> full_name
2736                and
2737                $self -> {'rerun'} < 1 ) {
2738
2739             if( not -e './NM_run' . ($run+1) . '/done' ){
2740               # here we have an .lst file, no done file and we are not
2741               # rerunning NONMEM. Which means we must create fake NM_run and
2742               # "done" files. (Usually this case occurs if we want to
2743               # use execute to create a summary or run table).
2744
2745               mkdir( "./NM_run" . ($run+1) );
2746               open( DONE, ">./NM_run". ($run+1) ."/done.1" );
2747               print DONE "This is faked\nseed: 1 1\n" ;
2748               close( DONE );
2749
2750               # TODO Must copy tablefiles if they exist.
2751
2752               # We use the existing .lst file as the final product.
2753               cp( $self -> {'models'} -> [$run] -> outputs -> [0] -> full_name, './NM_run' . ($run+1) .'/psn.lst' );
2754             }
2755
2756             # TODO Should check for tablefiles.
2757
2758             my $modulus = (($#models+1) <= 10) ? 1 : (($#models+1) / 10)+1;
2759
2760             if ( $run % $modulus == 0 or $run == 0 or $run == $#models ) {
2761               ui -> print( category => 'all', wrap => 0, newline => 0,
2762                            message  => 'D:'.( $run + 1 ).' .. ' )
2763                   unless( $self -> {'parent_threads'} > 1 or $self -> {'verbose'} );
2764             }
2765
2766             $queue_info{$run}{'candidate_model'} =
2767                 model -> new( filename             => "./NM_run".($run+1)."/psn.mod",
2768                               target               => 'disk',
2769                               ignore_missing_files => 1,
2770                               quick_reload         => 1,
2771                               cwres                => $models[$run] -> cwres() );
2772             $self -> print_finish_message( candidate_model => $queue_info{$run}{'candidate_model'},
2773                                            run => $run );
2774
2775             push( @{$self -> {'prepared_models'}[$run]{'own'}}, $queue_info{$run}{'candidate_model'} );
2776             push( @{$self -> {'ntries'}}, 0 );
2777             # push( @{$self -> {'evals'}}, \@evals );
2778
2779             next; # We are done with this model. It has already been run. Go back to polling.
2780
2781           }
2782
2783           # }}}
2784
2785           # {{{ delay code (to avoid overload of new processes)
2786
2787           if( $threads > 1 ){
2788
2789             if( $run > 1 ){
2790               my $start_sleep = Time::HiRes::time();
2791
2792               my ($min_sleep, $max_sleep); # min_sleep is in microseconds and max_sleep is in seconds.
2793
2794               if( defined $PsN::config -> {'_'} -> {'min_fork_delay'} ) {
2795                 $min_sleep = $PsN::config -> {'_'} -> {'min_fork_delay'};
2796               } else {
2797                 $min_sleep = 0;
2798               }
2799
2800               if( defined $PsN::config -> {'_'} -> {'max_fork_delay'} ) {
2801                 $max_sleep = $PsN::config -> {'_'} -> {'max_fork_delay'};
2802               } else {
2803                 $max_sleep = 0;
2804               }
2805
2806               if( $min_sleep > $max_sleep * 1000000 ){
2807                 $max_sleep = $min_sleep;
2808               }
2809
2810               # Dont wait for psn.lst if clean >= 3 it might have been
2811               # removed.
2812
2813               while( not( -e 'NM_run'.($run-1).'/psn.lst' ) and not $self -> {'clean'} >= 3
2814                      and
2815                      (Time::HiRes::time() - $start_sleep) < $max_sleep ){
2816                 Time::HiRes::usleep($min_sleep);
2817                 }
2818             }
2819
2820           }
2821
2822           # }}}
2823
2824           # {{{ Call to run_nonmem
2825
2826           # This will stop nasty prints from model, output and data
2827           # which are set to print for the scm.
2828           my $old_category = ui -> category();
2829
2830           ui -> category( 'modelfit' );
2831
2832           $self -> create_sub_dir( subDir => '/NM_run'.($run+1) );
2833           chdir( 'NM_run'.($run+1) );
2834
2835           # Initiate queue_info entry (unless its a restart)
2836
2837           unless( exists $queue_info{$run} ){
2838             $queue_info{$run}{'candidate_model'} = $self -> copy_model_and_input( model => $models[$run], source => '../' );
2839             $queue_info{$run}{'model'} = $models[$run];
2840             $queue_info{$run}{'modelfile_tainted'} = 1;
2841             $queue_info{$run}{'tries'} = 0;
2842             $queue_info{$run}{'crashes'} = 0;
2843             $queue_info{$run}{'evals'} = 0;
2844             $queue_info{$run}{'run_results'} = [];
2845             $queue_info{$run}{'raw_results'} = [];
2846             $queue_info{$run}{'raw_nonp_results'} = [];
2847             $queue_info{$run}{'start_time'} = 0;
2848
2849             # {{{ printing progress
2850
2851             # We don't want to print all starting models if they are
2852             # more than ten. But we allways want to print the first
2853             # and last
2854
2855             my $modulus = (($#models+1) <= 10) ? 1 : (($#models+1) / 10);
2856
2857             if ( $run % $modulus == 0 or $run == 0 or $run == $#models ) {
2858
2859               # The unless checks if tool owning the modelfit is
2860               # running more modelfits, in wich case we should be
2861               # silent to avoid confusion. The $done check is made do
2862               # diffrentiate between allready run processes.
2863
2864               ui -> print( category => 'all', wrap => 0, newline => 0,
2865                            message  => 'S:'.($run+1).' .. ' )
2866                   unless ( $self -> {'parent_threads'} > 1 or $self -> {'verbose'} );
2867             }
2868
2869             # }}}
2870           }
2871
2872           my %options_hash = %{$self -> _get_run_options(run_id => $run)};
2873
2874           $self -> run_nonmem ( run_no          => $run,
2875                                 nm_version      => $options_hash{'nm_version'},
2876                                 queue_info      => $queue_info{$run},
2877                                 queue_map       => \%queue_map );
2878
2879           chdir( '..' );
2880
2881           if( $self -> {'run_on_umbrella'} ){
2882             $queue_info{'NM_run'.($run+1)}{'run_no'} = $run;
2883           }
2884
2885           ui -> category( $old_category );
2886
2887           # }}}
2888
2889         } else {
2890
2891           # {{{ Here, a process has finished and we check for restarts.
2892
2893           my $run = $queue_map{$pid};
2894
2895           my $candidate_model = $queue_info{$run}{'candidate_model'};
2896
2897           my $work_dir = 'NM_run' . ($run+1) ;
2898           chdir( $work_dir );
2899
2900 #         if( $queue_info{$run}{'compile_only'} ){
2901 #           $queue_info{$run}{'compile_only'} = 0;
2902 #           $self -> run_nonmem ( run_no          => $run,
2903 #                                 nm_version      => $options_hash{'nm_version'},
2904 #                                 queue_info      => $queue_info{$run},
2905 #                                 queue_map       => \%queue_map );
2906 #
2907 #           delete( $queue_map{$pid} );
2908 #           chdir( '..' );
2909 #         } else {
2910             $self -> compute_cwres( queue_info => $queue_info{$run},
2911                                     run_no     => $run );
2912
2913             $self -> compute_iofv( queue_info => $queue_info{$run},
2914                                    run_no     => $run );
2915
2916             # Make sure that each process gets a unique random sequence:
2917             my $tmpseed = defined $self->seed() ? $self->seed() : random_uniform_integer(1,1,99999999);
2918             my $tmptry  = exists $queue_info{$run}{'tries'} ? $queue_info{$run}{'tries'} : 0;
2919             random_set_seed(($tmpseed+100000*($run+1)),($tmptry+1));
2920
2921             my %options_hash = %{$self -> _get_run_options(run_id => $run)};
2922
2923             my $do_restart = $self -> restart_needed( %options_hash,
2924                                                       queue_info        => $queue_info{$run},
2925                                                       run_no            => $run );
2926
2927             if( $do_restart ){
2928               unshift(@queue, $run);
2929               delete($queue_map{$pid});
2930               chdir( '..' );
2931             } else {
2932               $self -> select_best_model(run_no          => $run,
2933                                          nm_version      => $options_hash{'nm_version'},
2934                                          queue_info      => $queue_info{$run});
2935
2936               # {{{ Print finishing messages
2937
2938               if( scalar @queue == 0 ){
2939                 if( $all_jobs_started == 0 ){
2940
2941                   ui -> print( category => 'all',
2942                                message  => " done\n") if( $self -> {'parent_threads'} <= 1 and not $self -> {'verbose'} );
2943                   ui -> print( category => 'all',
2944                                message  => "Waiting for all NONMEM runs to finish:\n" ) if( $self -> {'parent_threads'} <= 1 and $threads > 1 and not $self -> {'verbose'} );
2945
2946                   $all_jobs_started = 1;
2947                 }
2948
2949                 my $modulus = (($#models+1) <= 10) ? 1 : (($#models+1) / 10)+1;
2950
2951                 if ( $run % $modulus == 0 or $run == 0 or $run == $#models ) {
2952                   ui -> print( category => 'all',
2953                                message  => 'F:'.($run+1).' .. ',
2954                                wrap => 0,
2955                                newline => 0)
2956                       unless ( $self -> {'parent_threads'} > 1 or $self -> {'verbose'} );
2957                 }
2958               }
2959
2960               # }}}
2961
2962               chdir( '..' );
2963
2964               # {{{ cleaning and done file
2965
2966               if( $self -> {'clean'} >= 3 ){
2967                 unlink( <$work_dir/*> );
2968                 unless( rmdir( $work_dir ) ){debug -> warn( message => "Unable to remove $work_dir directory: $! ." )};
2969                 sleep(2);
2970               } else {
2971
2972 #               unless( $queue_info{$run}{'tries'} > 0 ){
2973 #                 $self -> _print_done_file( run       => $run,
2974 #                                            tries     => $queue_info{$run}{'tries'},
2975 #                                            evals_ref => \@evals );
2976 #               }
2977
2978               }
2979
2980               # }}}
2981
2982               $self -> print_finish_message( candidate_model => $candidate_model,
2983                                              run => $run );
2984
2985               if( $threads <= 1 ){
2986                 push( @{$self -> {'prepared_models'}[$run]{'own'}}, $candidate_model );
2987                 push( @{$self -> {'ntries'}}, $queue_info{$run}{'tries'} );
2988               }
2989
2990               delete( $queue_info{$run} );
2991               delete( $queue_map{$pid} );
2992             }
2993
2994             next;
2995 #         }
2996
2997           # }}}
2998
2999         }
3000       } # queue loop ends here
3001
3002       if( $self -> {'run_on_umbrella'} ){
3003         # {{{ Umbrella code
3004
3005         $self -> umbrella_submit( prepare_jobs => 0 );
3006
3007         my $start_time = time();
3008
3009         my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
3010                                  ";databse=".$PsN::config -> {'_'} -> {'project'},
3011                                  $PsN::config -> {'_'} -> {'user'},
3012                                  $PsN::config -> {'_'} -> {'password'},
3013                                  {'RaiseError' => 1});
3014
3015         while( scalar keys %queue_info ){
3016
3017           # Get new "finished"
3018
3019           if( time() > ($start_time + $self -> {'umbrella_timeout'}) ) {
3020             debug -> die( message => "Waiting for NONMEM run timed out using the umbrella system" );
3021             last;
3022           }
3023           my $project = $PsN::config -> {'_'} -> {'project'};
3024           my $sth = $dbh -> prepare( "SELECT t.directory,j.tool_id, t.parent_tool_id, t.tool_id,j.state ".
3025                                      "FROM $project.job j, $project.tool t ".
3026                                      "WHERE j.tool_id=t.tool_id ".
3027                                      "AND t.parent_tool_id = ". $self -> tool_id() ." ".
3028                                      "AND j.state = 4" );
3029
3030           $sth -> execute or 'debug' -> die( message => $sth->errstr ) ;
3031
3032           my $select_arr = $sth -> fetchall_arrayref;
3033           $sth -> finish;
3034
3035           foreach my $row ( @{$select_arr} ){
3036             if( exists $queue_info{$row -> [0]} ){
3037
3038               my $work_dir = $row -> [0];
3039               my $id = $work_dir;
3040               chdir( $work_dir );
3041               my $modelfile_tainted = 0;
3042
3043               unless( defined $queue_info{$id}{'tries'} ){
3044                 $queue_info{$id}{'tries'} = 0;
3045               }
3046
3047               my %options_hash = %{$self -> _get_run_options(run_id => $queue_info{$id}{'run_no'})};
3048
3049               if( $self -> restart_needed( %options_hash,
3050                                            queue_info    => $queue_info{$id},
3051                                            run_no            => $id ) ){
3052
3053                 my $new_job_id;
3054
3055                 $self -> run_nonmem ( nm_version        => $options_hash{'nm_version'},
3056                                       model             => $models[$queue_info{$id}{'run_no'}],
3057                                       candidate_model   => $queue_info{$id}{'candidate_model'},
3058                                       work_dir          => $work_dir,
3059                                       tries             => \$queue_info{$id}{'tries'},
3060                                       job_id            => \$new_job_id,
3061                                       modelfile_tainted => \$modelfile_tainted);
3062
3063                 %{$queue_info{$new_job_id}} = %{$queue_info{$id}};
3064
3065                 delete $queue_info{$id};
3066               } else {
3067
3068                 $self -> select_best_model(run_no          => $queue_info{$id}{'run_no'},
3069                                            nm_version      => $options_hash{'nm_version'},
3070                                            queue_info      => $queue_info{$id});
3071                 delete $queue_info{$id};
3072
3073               }
3074
3075               chdir( '..' );
3076
3077             }
3078           }
3079         }
3080         $dbh -> disconnect;
3081
3082         # }}}
3083       }
3084
3085       # {{{ Print finish message
3086
3087       ui -> print( category => 'all',
3088                    message  => " done\n" ) if( $self -> {'parent_threads'} <= 1 and $threads > 1 and not $self -> {'verbose'});
3089
3090       # }}}
3091
3092       # }}} local execution
3093
3094       # {{{ Database code
3095
3096 #       my $dbh;
3097 #       if ( $PsN::config -> {'_'} -> {'use_database'} ) {
3098 #       $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.";databse=psn",
3099 #                             "psn", "psn_test",
3100 #                             {'RaiseError' => 1});
3101
3102 #       for( my $i = 1; $i <= scalar @{$self -> {'models'}}; $i++ ) {
3103 #         my $sth;
3104 #         foreach my $model ( @{$self -> {'prepared_models'}[$i-1]{'own'}} ) {
3105 #           my $model_id = $model -> model_id;
3106 #           if ( defined $model_id ) {
3107 #             $sth = $dbh -> prepare("INSERT INTO psn.tool_model (tool_id,".
3108 #                                    "model_id,prepared_model) ".
3109 #                                    "VALUES ('".$self -> {'tool_id'}."', ".
3110 #                                    $model_id.", 1 )");
3111 #             $sth -> execute;
3112 #           }
3113 #         }
3114 #         $sth -> finish if ( defined $sth );
3115 #       }
3116 #       $dbh -> disconnect;
3117 #       }
3118
3119       # }}}
3120
3121       $self -> prepare_raw_results();
3122
3123       $self -> print_raw_results();
3124
3125       chdir( $cwd );
3126
3127       # {{{ clean $self -> directory
3128       if( not $self -> {'top_tool'} and $self -> {'clean'} >= 3 ){
3129         my $dir = $self -> {'directory'};
3130         unlink( <$dir/NM_run*/*> );
3131         for( <$dir/NM_run*> ) {
3132           rmdir();
3133         }
3134         unlink( <$dir/*> );
3135         rmdir( $dir );
3136       }
3137       # }}}
3138
3139     }
3140 end run
3141
3142 # }}} run
3143
3144 # {{{ _get_run_options
3145
3146 start _get_run_options
3147 {
3148
3149   %options_hash = ( 'handle_rounding_errors' => undef, # <- Handle rounding errors a bit more intelligently
3150                     'handle_hessian_npd' => undef,        # <- Handle hessian not postiv definite a bit more intelligently
3151                     'handle_maxevals' => undef,        # <- Restart after maxeval
3152                     'cutoff_thetas' => 'ARRAY',
3153                     'tweak_inits' => undef,
3154                     'retries' => undef,
3155                     'picky' => undef,
3156                     'significant_digits_rerun' => undef,
3157                     'nm_version' => undef );
3158
3159   foreach my $option ( keys %options_hash ) {
3160
3161     # This loops allows run specific parameters to be
3162     # specified. We check that the parameter given is an
3163     # array we take out one element of that array and pass it
3164     # as a run specific parameter. If the option is specified
3165     # as being an "ARRAY" type in the has above, but there are
3166     # no subarrays, we send the entire array as an parameter.
3167
3168     if( ref( $self -> {$option} ) eq 'ARRAY' ) {
3169       if( ref( $self -> {$option} -> [$run_id] ) ne 'ARRAY' and $options_hash{$option} eq 'ARRAY' ) {
3170         $options_hash{$option} = $self -> {$option};
3171       } else {
3172         $options_hash{$option} = $self -> {$option} -> [$run_id];
3173       }
3174     } else {
3175       $options_hash{$option} = $self -> {$option};
3176     }
3177   }
3178 }
3179 end _get_run_options
3180
3181 # }}}
3182
3183 # {{{ summarize
3184
3185 start summarize
3186     {
3187
3188       my %raw_header_map;
3189
3190       for( my $i = 0 ; $i <= scalar @{$self -> {'raw_results_header'}}; $i++ ){
3191         $raw_header_map{$self -> {'raw_results_header'} -> [$i]} = $i;
3192       }
3193
3194       my %model_problem_structure;
3195
3196       foreach my $raw_row ( @{$self -> {'raw_results'}} ){
3197         my $model = $raw_row -> [$raw_header_map{'model'}];
3198         my $problem = $raw_row -> [$raw_header_map{'problem'}];
3199         my $subproblem = $raw_row -> [$raw_header_map{'subproblem'}];
3200
3201         $model_problem_structure{ $model }{ $problem }{ $subproblem } = 1;
3202       }
3203
3204       # ------------------------------------------------------------------#
3205
3206       # Model evaluation. Part one. Version 0.1
3207
3208       # This method checks a nonmem output file for flaws.  The parts
3209       # concerning parts of the covariance step obviously needs the
3210       # $COV # record. In addition, the conditioning number test needs
3211       # the PRINT=E # option of the $COV record. The limits for the
3212       # relative tests may be # changed below:
3213
3214       my $correlation_limit = $self -> {'correlation_limit'};      # All correlations above this number
3215                                                                    # will be listed.
3216       my $condition_number_limit = $self -> {'condition_number_limit'};    # An error will be rised if the
3217                                                                            # condition number is greater
3218                                                                            # than this number.
3219       my $near_bound_sign_digits = $self -> {'near_bound_sign_digits'};   # If a parameter estimate is equal
3220                                                                           # to a bound with this many
3221                                                                           # significant digits, a warning
3222                                                                           # will be printed.
3223       my $near_zero_boundary_limit = $self -> {'near_zero_boundary_limit'};  # When the bound is zero, the
3224                                                                              # check above is not valid. Use
3225                                                                              # this limit instead.
3226       my $sign_digits_off_diagonals =  $self -> {'sign_digits_off_diagonals'};   # The off-diagonal elements are
3227                                                                                  # checked against +-1 with this
3228                                                                                  # many significant digits.
3229       my $large_theta_cv_limit = $self -> {'large_theta_cv_limit'};   # Coefficients of variation larger
3230                                                                       # than these numbers will produce
3231                                                                       # warnings.
3232       my $large_omega_cv_limit = $self -> {'large_omega_cv_limit'};
3233       my $large_sigma_cv_limit = $self -> {'large_omega_cv_limit'};
3234
3235       my $confidence_level = $self -> {'confidence_level'};   # Percentage value;
3236       my $precision = $self -> {'precision'};   # Precision in output.
3237
3238       sub acknowledge {
3239         my $name    = shift;
3240         my $outcome = shift;
3241         my $file    = shift;
3242         my $l = (7 - length( $outcome ))/2;
3243         my $text = sprintf( "%-66s%2s%7s%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
3244         ui -> print( category => 'all', message => $text, wrap => 0 );
3245         print $file $text if defined $file;
3246       }
3247
3248       my $raw_line_structure = ext::Config::Tiny -> read( $self -> {'directory'} . "/raw_file_structure" );
3249       my $raw_res_index = 0;
3250
3251       for( my $model = 1; $model <= scalar keys %model_problem_structure; $model ++ ){
3252
3253 #       my $output = $self -> {'models'} -> [ $model - 1 ] -> outputs -> [0];
3254
3255 #       $output -> _read_problems;
3256
3257         open( my $log, ">test.log" );
3258
3259         my $eigen_run     = 1;
3260         my $corr_matr_run = 1;
3261
3262         my $raw_res = $self -> {'raw_results'};
3263
3264         for ( my $problem = 1; $problem <= scalar keys %{$model_problem_structure{ $model }} ; $problem++ ) {
3265           print "PROBLEM ",$problem,"\n" if scalar keys %{$model_problem_structure{ $model }} > 1;
3266           for ( my $subproblem = 1; $subproblem <= scalar keys %{$model_problem_structure{ $model }{ $problem }}; $subproblem++ ) {
3267             print "SUBPROBLEM ",$subproblem,"\n" if scalar keys %{$model_problem_structure{ $model }{ $problem }} > 1;
3268
3269             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'minimization_successful'}] == 1 ) { # Is in raw
3270               acknowledge( 'Successful minimization', 'OK', $log );
3271             } else {
3272               acknowledge( 'Termination problems', 'ERROR', $log );
3273             }
3274             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'rounding_errors'}] eq '0' ) { # Is in raw
3275               acknowledge( 'No rounding errors', 'OK', $log );
3276             } else {
3277               acknowledge( 'Rounding errors', 'ERROR', $log );
3278             }
3279
3280             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'zero_gradients'}] eq '0' ) { # Is in raw
3281               acknowledge( 'No zero gradients', 'OK', $log );
3282             } else {
3283               acknowledge( 'Zero gradients found '.
3284                            $raw_res -> [ $raw_res_index ][$raw_header_map{'zero_gradients'}].
3285                            ' times', 'ERROR', $log );
3286             }
3287
3288             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'final_zero_gradients'}] eq '0' ) { # Is in raw
3289               acknowledge( 'No final zero gradients', 'OK', $log );
3290             } else {
3291               acknowledge( 'Final zero gradients', 'ERROR', $log );
3292             }
3293
3294             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'hessian_reset'}] eq '0' ) { # Is in raw
3295               acknowledge( 'Hessian not reset', 'OK', $log );
3296             } else {
3297               acknowledge( 'Hessian reset '.
3298                            $raw_res -> [ $raw_res_index ][$raw_header_map{'hessian_reset'}].
3299                            ' times', 'WARNING', $log );
3300             }
3301
3302             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'estimate_near_boundary'}] eq '0' ){
3303               acknowledge( 'Estimate not near boundary', 'OK', $log );
3304             } else {
3305               acknowledge( 'Estimate near boundary','WARNING',$log);
3306             }
3307
3308             my ( $n_b, $f_b, $f_e );# =
3309 #               $output -> near_bounds( zero_limit         => $near_zero_boundary_limit,
3310 #                                       significant_digits => $near_bound_sign_digits,
3311 #                                       off_diagonal_sign_digits => $sign_digits_off_diagonals ); # Is in raw
3312
3313             if ( defined $n_b -> [$problem] and
3314                  defined $n_b -> [$problem][$subproblem] ) {
3315               my @near_bounds = @{$n_b -> [$problem][$subproblem]};
3316               my @found_bounds = @{$f_b -> [$problem][$subproblem]};
3317               my @found_estimates = @{$f_e -> [$problem][$subproblem]};
3318               if ( $#near_bounds < 0 ) {
3319                 acknowledge( 'No parameter near boundary', 'OK', $log );
3320               } else {
3321                 acknowledge( 'Parameter(s) near boundary', 'WARNING', $log );
3322                 for ( my $problem = 0; $problem <= $#near_bounds; $problem++ ) {
3323                   printf( "\t%-20s%10g   %10g\n", $near_bounds[$problem],
3324                           $found_estimates[$problem], $found_bounds[$problem] );
3325                   print $log sprintf( "\t%-20s%10g   %10g\n", $near_bounds[$problem],
3326                                       $found_estimates[$problem], $found_bounds[$problem] );
3327                   print "\n" if ( $problem == $#near_bounds );
3328                   print $log "\n" if ( $problem == $#near_bounds );
3329                 }
3330               }
3331             }
3332
3333             if ( $raw_res -> [ $raw_res_index ][ $raw_header_map{'covariance_step_run'} ] ) { # Is in raw
3334               if (  $raw_res -> [ $raw_res_index ][$raw_header_map{'covariance_step_successful'}] eq '0' ) {
3335                 acknowledge( 'Covariance step', 'ERROR', $log );
3336               } else {
3337                 if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'covariance_step_successful'}] eq '1' ) {
3338                   acknowledge( 'Covariance step ', 'OK', $log );
3339                 }else {
3340                   acknowledge( 'Covariance step', 'WARNING', $log );
3341                 }
3342
3343                 my ( $l_se, $f_cv );# =
3344 #                   $output -> large_standard_errors( theta_cv_limit => $large_theta_cv_limit,
3345 #                                                     omega_cv_limit => $large_omega_cv_limit,
3346 #                                                     sigma_cv_limit => $large_sigma_cv_limit ); # Is NOT in raw
3347                 if ( defined $l_se -> [$problem] and
3348                      defined $l_se -> [$problem][$subproblem] ) {
3349                   my @large_standard_errors = @{$l_se -> [$problem][$subproblem]};
3350                   my @found_cv = @{$f_cv -> [$problem][$subproblem]};
3351                   if ( $#large_standard_errors < 0 ) {
3352                     acknowledge( 'No large standard errors found', 'OK', $log );
3353                   } else {
3354                     acknowledge( 'Large standard errors found', 'WARNING', $log );
3355                     for ( my $problem = 0; $problem <= $#large_standard_errors; $problem++ ) {
3356                       printf( "\t%-20s%10g\n", $large_standard_errors[$problem], $found_cv[$problem] );
3357                       print $log sprintf( "\t%-20s%10g\n", $large_standard_errors[$problem], $found_cv[$problem] );
3358                       print "\n" if ( $problem == $#large_standard_errors );
3359                       print $log "\n" if ( $problem == $#large_standard_errors );
3360                     }
3361                   }
3362                 }
3363
3364                 my $eigens = $raw_res -> [ $raw_res_index ][$raw_header_map{'eigens'}]; # Is in raw
3365                 if ( $eigens ){# defined $eigens and
3366                     # defined $eigens -> [$problem][$subproblem] and
3367                     # scalar @{$eigens -> [$problem][$subproblem]} > 0 ) {
3368                   if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'condition_number'}] < $condition_number_limit ) {
3369                     acknowledge( 'Condition number ', 'OK', $log );
3370                   } else {
3371                     acknowledge( 'Large condition number ('.
3372                                  sprintf("%6.0f",
3373                                          $raw_res -> [ $raw_res_index ][$raw_header_map{'condition_number'}]).
3374                                  ')', 'WARNING', $log );
3375                   }
3376                 } else {
3377                   $eigen_run = 0;
3378                 }
3379
3380                 if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'s_matrix_singular'}] ne '1' ) { # Is in raw
3381                   acknowledge( 'S-matrix ', 'OK', $log );
3382                 } else {
3383                   acknowledge( 'S-matrix is singular', 'WARNING', $log );
3384                 }
3385
3386                 my ( $h_c, $f_c );# = $output high_correlations( limit => $correlation_limit ); # Is NOT in raw
3387                 if ( defined $h_c -> [$problem] and
3388                      defined $h_c -> [$problem][$subproblem] ) {
3389                   my @high_correlations = @{$h_c -> [$problem][$subproblem]};
3390                   my @found_correlations = @{$f_c -> [$problem][$subproblem]};
3391                   if ( $#high_correlations < 0 ) {
3392                     acknowledge( 'Correlations', 'OK', $log );
3393                   } else {
3394                     acknowledge( 'Large correlations between parameter estimates found', 'WARNING', $log );
3395                     for ( my $problem = 0; $problem <= $#high_correlations; $problem++ ) {
3396                       printf( "\t%-20s%10g\n", $high_correlations[$problem],
3397                               $found_correlations[$problem] );
3398                       print $log sprintf( "\t%-20s%10g\n", $high_correlations[$problem],
3399                                           $found_correlations[$problem] );
3400                       print "\n" if ( $problem == $#high_correlations );
3401                       print $log "\n" if ( $problem == $#high_correlations );
3402                     }
3403                   }
3404                 }
3405               }
3406             } else {
3407               print "No Covariance step run\n";
3408               print $log "No Covariance step run\n";
3409             }
3410
3411             # sumo script
3412
3413             my $confidence_interval = 0;
3414             my $csv = 0;
3415             my $form = '%.' . $precision . 'g';
3416
3417             my @output_matrix;
3418             my @output_matrix_sizes;
3419
3420             my %c_levels = ( '90'   => 1.6449,
3421                              '95'   => 1.96,
3422                              '99'   => 2.5758,
3423                              '99.9' => 3.2905 );
3424
3425             if( $confidence_interval ) {
3426               if( not defined $c_levels{$confidence_level} ) {
3427                 debug -> die( message =>  "Sorry, confidence intervals for level ".$confidence_level.
3428                               " can not be output. Valid levels are: ".join(',', keys %c_levels),
3429                               level => 1 );
3430
3431               }
3432             }
3433
3434             my ( %nam, %est, %cest, %ses );
3435
3436
3437             foreach my $estimate ( 'theta','omega','sigma' ){
3438               my @pos_length = split( /,/ , $raw_line_structure -> {($problem-1)}{$estimate});
3439               my @estimates = @{$raw_res -> [ $raw_res_index ]}[ $pos_length[0] .. $pos_length[0] + $pos_length[1] ];
3440               my @names;
3441               foreach my $num ( 0 .. $pos_length[1] - 1){
3442                 my $template = uc(substr($estimate,0,2));
3443                 push( @names, $template.($num+1));
3444               }
3445
3446               my @pos_length = split( /,/ , $raw_line_structure -> {($problem-1)}{'se'.$estimate});
3447               my @se_estimates = @{$raw_res -> [ $raw_res_index ]}[ $pos_length[0] .. $pos_length[0] + $pos_length[1] ];
3448
3449               my @pos_length = split( /,/ , $raw_line_structure -> {($problem-1)}{'cvse'.$estimate});
3450 #             my @cvse_estimates = @{$raw_res -> [ $raw_res_index ]}[ $pos_length[0] .. $pos_length[0] + $pos_length[1] ];
3451               my @cvse_estimates = ();
3452
3453               $nam{$estimate} = \@names;
3454               $est{$estimate} = \@estimates;
3455               $est{'cvse'.$estimate} = \@cvse_estimates;
3456               $ses{$estimate} = \@se_estimates;
3457             }
3458
3459             my $ofv = $raw_res -> [ $raw_res_index ][ $raw_header_map{'ofv'} ];
3460
3461             if ( defined $ofv ) {
3462               print "Objective function value: ",$ofv,"\n\n";
3463             } else {
3464               print "Objective function value: UNDEFINED\n\n";
3465             }
3466
3467             push( @output_matrix, ["","THETA","","","OMEGA","","","SIGMA", ""] );
3468             for( my $i = 0; $i <= $#{$output_matrix[0]}; $i++ ){
3469               if( $output_matrix_sizes[$i] < length( $output_matrix[0][$i] ) ){
3470                 $output_matrix_sizes[$i] = length( $output_matrix[0][$i] );
3471               }
3472             }
3473
3474             my $max_par = $#{$est{'theta'}};
3475             $max_par = $#{$est{'omega'}} if ( $#{$est{'omega'}} > $max_par );
3476             $max_par = $#{$est{'sigma'}} if ( $#{$est{'sigma'}} > $max_par );
3477
3478             for ( my $max = 0; $max <= $max_par; $max++ ) {
3479               my ( @row, %cis );
3480               if( $confidence_interval ) {
3481                 foreach my $param ( 'theta', 'omega', 'sigma' ) {
3482                   if ( defined $est{$param}[$max] ) {
3483                     my $diff = $c_levels{$confidence_level}*$ses{$param}[$max];
3484                     my ( $lo, $up );
3485                     if( defined $diff ) {
3486                       $lo = $est{$param}[$max]-$diff;
3487                       $up = $est{$param}[$max]+$diff;
3488                     }
3489                     my $cis = sprintf( "($form - $form)", $lo, $up );
3490                     push( @row, $nam{$param}[$max],
3491                           sprintf( $form, $est{$param}[$max] ),
3492                           $cis );
3493                   } else {
3494                     push( @row, '','','' );
3495                   }
3496                 }
3497               } else {
3498                 foreach my $estimate ( 'theta', 'omega','sigma' ){ # Loop over comegas instead
3499                   if ( defined $nam{$estimate}->[$max] ) {
3500                     push( @row, $nam{$estimate}->[$max], defined $est{$estimate}->[$max] ? sprintf( $form, $est{$estimate}->[$max] ) : '........',
3501                           $est{'cvse'.$estimate}->[$max] ? sprintf( "($form)", $est{'cvse'.$estimate}->[$max] ) : '(........)' );
3502                   } else {
3503                     push( @row, '','','' );
3504                   }
3505                 }
3506               }
3507
3508               push(@output_matrix, \@row);
3509               for( my $i = 0; $i <= $#row; $i++ ){
3510                 if( $output_matrix_sizes[$i] < length( $row[$i] ) ){
3511                   $output_matrix_sizes[$i] = length( $row[$i] );
3512                 }
3513               }
3514             }
3515
3516             foreach my $row ( @output_matrix ){
3517               for( my $i = 0; $i <= $#{$row}; $i++ ){
3518                 my $spaces = $output_matrix_sizes[$i] - length($row -> [$i]);
3519                 if( $csv ){
3520                   print $row -> [$i], ",";
3521                 } else {
3522                   print " " x $spaces, $row -> [$i], "  ";
3523                 }
3524               }
3525               print "\n";
3526             }
3527
3528             #$output -> flush;
3529             $raw_res_index ++;
3530           }
3531         }
3532         close( $log );
3533
3534       }
3535     }
3536 end summarize
3537
3538 # }}}
3539
3540 # {{{ compute_cwres
3541
3542 start compute_cwres
3543 my $queue_info_ref = defined $parm{'queue_info'} ? $parm{'queue_info'} : {};
3544 my $model = $queue_info_ref -> {'model'};
3545 if( defined $PsN::config -> {'_'} -> {'R'} ) {
3546   my $probs = $model -> problems();
3547   if( defined $probs ) {
3548     foreach my $prob ( @{$probs} ) {
3549       if( $prob -> {'cwres_modules'} ){
3550         my $sdno = $prob -> {'cwres_modules'} -> [0] -> sdno();
3551         my $sim = $prob -> {'cwres_modules'} -> [0] -> mirror_plots ? ',sim.suffix="sim"' : '';
3552 # Create the short R-script to compute the CWRES.
3553         open( RSCRIPT, ">compute_cwres.R" );
3554         print RSCRIPT "library(xpose4)\ncompute.cwres($sdno $sim)\n";
3555         close( RSCRIPT );
3556
3557 # Execute the script
3558         system( $PsN::config -> {'_'} -> {'R'}." CMD BATCH compute_cwres.R" );
3559       }
3560     }
3561   }
3562 }
3563
3564 end compute_cwres
3565
3566 # }}}
3567
3568 # {{{ compute_iofv
3569
3570 start compute_iofv
3571 {
3572   my $queue_info_ref = defined $parm{'queue_info'} ? $parm{'queue_info'} : {};
3573   my $model = $queue_info_ref -> {'model'};
3574
3575   if( defined $model -> iofv_modules ){
3576     $model -> iofv_modules -> [0] -> post_run_process;
3577   }
3578 }
3579 end compute_iofv
3580
3581 # }}}