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