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