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 'cp';
  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 ext::IPC::Run3;
  23         use output;
  24         use OSspecific;
  25         use ui;
  26         use moshog_client;
  27       }
  28 end include statements
  29
  30 # }}} include statements
  31
  32 # {{{ description, examples, synopsis, see also
  33
  34 # No method, just documentation
  35 start description
  36     #
  37     # In PsN versions < 2.0, the functionality for actually running
  38     # NONMEM on models and data PsN objects was provided by the model
  39     # class. As of PsN versions 2.0 and higher, this functinality has
  40     # been moved to the separate class I<modelfit> in order to make the
  41     # class responsibilities clearer.
  42     #
  43     # Fitting a model can be viewed as a special case of a more
  44     # general set of tools for population PK/PD. In PsN, the
  45     # modelfit class is therefore a specialization of a general PsN
  46     # tool class. The tool class itself is not capable of much at
  47     # all but to define a common structure for all PsN tools.
  48     #
  49     # All attributes and (class) methods specified for the general
  50     # tool class are inherited in the modelfit class. Some (class) methods
  51     # are defined in both classes (e.g. the L</run>) and in these
  52     # cases it is the modelfit version that will be used.
  53     #
  54     # =begin html
  55     #
  56     # <tr>Please look at the documentation for the <a
  57     # href="../tool.html">general tool class</a> as well as the <a
  58     # href="#examples">examples</a> section of this document for
  59     # descriptions of the setting that are available for all
  60     # tools.</tr>
  61     #
  62     # =end html
  63     #
  64     # =begin man
  65     #
  66     # Please look at the documentation for the general I<tool> class
  67     # as well as the L</examples> section of this document for
  68     # descriptions of the setting that are available for all tools.
  69     #
  70     # =end man
  71     #
  72 end description
  73
  74 start examples
  75     # The following code may be used to create a simple modelfit
  76     # object and to run it.
  77     #
  78     #   use lib qw(path to PsN installation/lib);
  79     #   use tool::modelfit;
  80     #   use model;
  81     #
  82     #   my $model = model -> new( filename => 'run1.mod' );
  83     #   my $modelfit = tool::modelfit -> new( models => [$model] );
  84     #   my %results = %{$modelfit -> run};
  85     #
  86     # To illustrate a more complex use of modelfit, we can run a
  87     # bootstrap of 200 resampled data sets of same size as the
  88     # original data set. In this example we assume that the modelfile
  89     # we use contains only one problem and one sub problem.
  90     #
  91     #   use lib qw(path to PsN installation/lib);
  92     #   use tool::modelfit;
  93     #   use data;
  94     #   use model;
  95     #   use output;
  96     #   Math::Random;       # For perturbation of initial estimates
  97     #                       # after unsuccessful minimizations
  98     #
  99     #   # set these to appropriate values for your own run
 100     #   my $samples = 200;
 101     #   my $modelfile = 'run1.mod';
 102     #   my $boot_sample_file = 'boot_ind.csv';
 103     #   my $boot_results_file = 'boot_results.csv';
 104     #
 105     #   # set the seed from a phrase (consistent over different hardware,
 106     #   # see Math::Random )
 107     #   random_set_seed_from_phrase('testing');
 108     #
 109     #   # ignore missing data and (especially) output files
 110     #   my $model = model -> new( filename => $modelfile,
 111     #                             ignore_missing_files => 1 );
 112     #
 113     #   my $data = $model -> datas -> [0];
 114     #
 115     #   # Create the bootstrap data sets. The default name for each
 116     #   # new resampled data file will be bs###.dta, where ### is a
 117     #   # number from 1 to $samples.
 118     #   my ( $dataref, $incl_ind_ref, $incl_keys_ref ) =
 119     #       $data -> bootstrap ( samples => $samples );
 120     #
 121     #   # Save the resampled ID numbers in a file, one row for each
 122     #   # bootstrap data set
 123     #   open ( LOG, ">$boot_sample_file" );
 124     #   foreach my $sample ( @{$incl_ind_ref} ) {
 125     #     print LOG join(';', @{$sample} ),"\n";
 126     #   }
 127     #
 128     #   # Create the boostrap models and set the correct resampled
 129     #   # data file in $DATA. Save them in @bs_models.
 130     #   my @bs_models  = ();
 131     #   for ( my $i = 1; $i <= $samples; $i++ ) {
 132     #     my $bs_mod = $model -> copy( filename    => 'bs'.$i.'.mod',
 133     #                                  copy_data   => 0,
 134     #                                  copy_output => 0);
 135     #     $bs_mod -> datafiles( new_names     => ['bs'.$i.'.dta'],
 136     #                           absolute_path => 1 );
 137     #     $bs_mod -> _write;
 138     #     push( @bs_models, $bs_mod );
 139     #   }
 140     #
 141     #   # Create a modelfit object with the bootstrap models as
 142     #   # input. Set the number of parallel runs to 2.
 143     #   my $mfitobj = tool::modelfit -> new ( models     => \@bs_models,
 144     #                                           threads    => 2 );
 145     #
 146     #   # Run the model fit. Since the bootstrap models are named
 147     #   # bs###.mod, the default names for the output files will be
 148     #   # bs###.lst.
 149     #   $mfitobj -> run;
 150     #
 151     #   # We'll save the OFV plus the theta, omega and sigma estimates
 152     #   # for each run in a file.
 153     #   open( RESULT, ">$boot_results_file" );
 154     #
 155     #   for ( my $i = 1; $i <= $samples; $i++ ) {
 156     #     my $bs_out = output -> new( filename => 'bs'.$i.'.lst' );
 157     #     my @ofv    = @{$bs_out -> ofv};
 158     #     my @thetas = @{$bs_out -> thetas};
 159     #     my @omegas = @{$bs_out -> omegas};
 160     #     my @sigmas = @{$bs_out -> sigmas};
 161     #     # We know that we only have one problem and one sub problem in this
 162     #     # example. ([0][0] comes from that fact)
 163     #     my @print_strings = ();
 164     #     push( @print_strings, $ofv[0][0] );
 165     #     push( @print_strings, @{$thetas[0][0]} );
 166     #     push( @print_strings, @{$omegas[0][0]} );
 167     #     push( @print_strings, @{$sigmas[0][0]} );
 168     #     print RESULT join( ';', @print_strings ),"\n";
 169     #   }
 170     #
 171     #   close( RESULT );
 172     #
 173     #   # We're done!
 174 end examples
 175
 176 start synopsis
 177     #   use tool::modelfit;
 178     #   use model;
 179     #
 180     #   my $model_obj = model -> new ( filename => 'run1.mod' );
 181     #
 182     #   my $modelfit_obj = tool::modelfit -> new ( models => [$model_obj] );
 183     #
 184     #   my $output_obj = $modelfit_obj -> run;
 185 end synopsis
 186
 187 start see_also
 188     # =begin html
 189     #
 190     # <a HREF="../data.html">data</a>, <a
 191     # HREF="../model.html">model</a> <a
 192     # HREF="../output.html">output</a>, <a
 193     # HREF="../tool.html">tool</a>
 194     #
 195     # =end html
 196     #
 197     # =begin man
 198     #
 199     # data, model, output, tool
 200     #
 201     # =end man
 202 end see_also
 203
 204 # }}}
 205
 206 # {{{ new
 207
 208 start new
 209     # Usage:
 210     #
 211     #   $modelfit_object = tool::modelfit -> new( models  => [$model_object],
 212     #                                               retries => 5 );
 213     #
 214     # This is the basic usage and it creates a modelfit object that
 215     # can later be run using the L</run> method. I<models> is an array
 216     # of PsN model objects.
 217     #
 218     #   $modelfitObject = $tool::modelfit -> new( 'retries' => 5 );
 219     #   $modelfitObject -> add_model( init_data => { filename => $modelFileName } );
 220     #
 221     # This way of using modelfit is suitable if you have a list with
 222     # filenames of modelfiles. "add_model> will create modelfitobject
 223     # for you.
 224     #
 225     # A more interresting attribute is I<threads> which sets how many
 226     # parallel executions of NONMEM that will run. Some tips are:
 227     # Setting the number of threads higher than the number of nodes in
 228     # your cluster/supercomputer can make your runs slower. The
 229     # biggest limiting factor is the amount of memory needed by
 230     # NONMEM. With smaller runs, just set the thread number to the
 231     # number of nodes available.
 232     #
 233     # The I<directory> is the folder where the tools stores
 234     # temporary data and runs subtools (or in the modelfit case,
 235     # runs NONMEM). Each NONMEM run will have its own sub directory
 236     # NM_run[X] where [X] is an index running from 1 to the number of
 237     # runs. If unsure of what this means, leave it undefined and a
 238     # default will be used, e.g. modelfit_dir3 or something.
 239     #
 240     # Next, the I<compress> and I<remove_temp_files> attributes are good
 241     # if you want to save some hard disk space. I<compress> set to 1
 242     # will put all NONMEM output in to an tar/gz archive named
 243     # I<nonmem_files.tgz> placed in the I<NM_run[X]> directory
 244     # described above. If I<remove_temp_files> is set to 1,  the NONMEM
 245     # files: 'FCON', 'FDATA', 'FSTREAM', 'PRDERR' will be removed.
 246     #
 247     # I<clean> is a stronger version of I<remove_temp_files>; it will also
 248     # remove I<NM_run[X]> and all that is in these.
 249     #
 250     # I<retries> is the number of times L</run> will alter initial
 251     # values and (re)execute NONMEM when executions fail. I<retries>
 252     # can either be an integer, specifying the number of retries for
 253     # all models, or it can be an array with the number of retries
 254     # specific for each modelfile as elements. The default value is
 255     # B<5>. The algorithm for altering the initial values works
 256     # roughly like this: For each each new try, a random new initial
 257     # value is drawn from a uniform distribution with limits +-n*10%
 258     # of the original intial estimate and where n i equal to the retry
 259     # number. I.e. the first retry, the borders of the distribution
 260     # are +-10%. The algorithm ensures that the new values are within
 261     # specified boundaries.
 262     #
 263     # =begin html
 264     #
 265     # For a full dexcription of the algorithm, see <a
 266     # href="../model/problem/record/init_option.html#set_random_init">set_random_init</a>
 267     # of the <a
 268     # href="../model/problem/record/init_option.html">init_option
 269     # class</a>.
 270     #
 271     # =end html
 272     #
 273     # =begin man
 274     #
 275     # For a full dexcription of the algorithm, see I<set_random_init>
 276     # of the I<init_option> class.
 277     #
 278     # =end man
 279     #
 280     # If I<picky> is set to 1, the output from NONMEM will be checked
 281     # more thoroughly. If any of the lines below are found in the
 282     # minimization message, a rerun is initiated.
 283     #
 284     #    COVARIANCE STEP ABORTED
 285     #    PROGRAM TERMINATED BY OBJ
 286     #    ESTIMATE OF THETA IS NEAR THE BOUNDARY AND
 287     #    PARAMETER ESTIMATE IS NEAR ITS BOUNDARY
 288     #    R MATRIX ALGORITHMICALLY SINGULAR
 289     #    S MATRIX ALGORITHMICALLY SINGULAR
 290     #
 291     # I<nm_version> is a string with the version number of NONMEM that
 292     # will be used. The installed versions of NONMEM must be specified
 293     # in OSspecific.pm, the class responsible for system specific
 294     # features settings.
 295     #
 296     # I<logfile> specifies the name of the logfile.
 297     #
 298     # if I<silent_logfile> is defined all NONMEM output will
 299     # be written to I<NM_run[X]/xxx>, where xxx is the defined value.
 300     #
 301     # I<extra_files> is an array of strings where each string is a
 302     # file needed for NONMEM execution. Those file will be moved
 303     # to the I<NM_run[X]> directory.
 304     #
 305     # I<seed> is just a way to set a random seed number.
 306     #
 307     # If I<run_on_nordugrid> is set to true modelfit will submit the nonmem
 308     # runs to a grid. A group of related parameters are also
 309     # specified.
 310
 311     # I<cpuTime> Is an estimated execution time for each individual
 312     # modelfile. It should preferably be a bit longer than reality. If
 313     # you specify a cpuTime that is to short, you risk that the grid
 314     # kills your jobs prematurely. The unit of I<cpuTime> is minutes.
 315
 316     # I<grid_cpuTime> is the time of the actual grid job. It should be
 317     # used to group modelfiles together. For example, if you set
 318     # I<cpuTime> to ten minutes, I<grid_cpuTime> to 60 minutes and the
 319     # number of modelfiles is 14 modelfit will create three grid jobs,
 320     # two with six model files each and one with two modelfiles.
 321
 322     # I<grid_adress> is the URL of the grid submission server,
 323     # e.g. hagrid.it.uu.se.
 324
 325
 326       {
 327         if ( defined $this -> {'logfile'} ) {
 328           my $dir;
 329           $this -> {'logfile'} = join('', OSspecific::absolute_path( $this -> {'directory'},
 330                                                                      $this -> {'logfile'}) );
 331         }
 332         if ( $this -> {'ask_if_fail'} ) {
 333           eval( 'use Tk' );
 334         }
 335       }
 336
 337 end new
 338
 339 # }}} new
 340
 341 # {{{ prepare_raw_results
 342
 343 start prepare_raw_results
 344       {
 345         # As of PsN version 2.1.8, we don't handle problems and
 346         # subproblems in any of the tools but modelfit.
 347         my @params = ( 'theta', 'omega', 'sigma', 'npomega', 'shrinkage_etas', 'shrinkage_wres','eigen' );
 348
 349         # ----------------------  Collect the parameter names  ----------------------
 350
 351         # {{{ param names
 352
 353         # Find the model with most parameter of each kind and put the
 354         # number of the parameter in "max_hash".
 355
 356         my ( %max_hash ); #, %category_names );
 357         for ( my $i = 0; $i < scalar @{$self -> {'models'}}; $i++ ) { # models level
 358
 359           # Some column in "raw_results_header" are meta-columns, they
 360           # will be replaced by several columns. For example, the
 361           # 'theta' column will be replaced with TH1, TH2, TH3 in the
 362           # general case. (Actually it will be replaced with the
 363           # thetas labels from the model file. Here we search the meta
 364           # column of the raw_results_header to find the maximum
 365           # number of real columns.
 366
 367           foreach my $category ( @{$self -> {'raw_results_header'}},'npomega' ) {
 368             if ( $category eq 'setheta' or $category eq 'seomega' or $category eq 'sesigma' ){
 369               next;
 370             } elsif ( $category eq 'theta' or $category eq 'omega' or $category eq 'sigma' or
 371                  $category eq 'npomega' or $category eq 'shrinkage_etas' or $category eq 'eigen') {
 372               my $numpar = 0;
 373               if( $category eq 'npomega' or $category eq 'shrinkage_etas' ) {
 374                 my $nomegas = $self -> {'models'} -> [$i] -> nomegas;
 375                 if( defined $nomegas ) {
 376                   for( my $j = 0; $j < scalar @{$nomegas}; $j++ ) {
 377                     if( $category eq 'npomega' ) {
 378                       my $npar = $nomegas -> [$j];
 379                       $npar = $npar*($npar+1)/2;
 380                       $numpar = $numpar >= $npar ? $numpar : $npar;
 381                     } else {
 382                       $numpar = $numpar >= $nomegas -> [$j] ? $numpar : $nomegas -> [$j];
 383                     }
 384                   }
 385                 }
 386               } elsif( $category eq 'eigen') {
 387                 foreach my $prob_eigens ( @{$self -> {'models'} -> [$i] -> outputs -> [0] -> eigens} ){
 388                   foreach my $subprob_eigens( @{$prob_eigens} ){
 389                     if( defined $subprob_eigens ){
 390                       $numpar = $numpar < scalar @{$subprob_eigens} ? scalar @{$subprob_eigens} : $numpar ;
 391                     }
 392                   }
 393                 }
 394               } else {
 395                 my $labels = $self -> {'models'} -> [$i] -> {'outputs'} -> [0] -> labels( parameter_type => $category );
 396
 397                 # we can't use labels directly since different models may have different
 398                 # labels (still within the same modelfit)
 399                 if( defined $labels ) {
 400                   for( my $j = 0; $j < scalar @{$labels}; $j++ ) {
 401                     if( defined $labels -> [$j] ) {
 402                       $numpar = $numpar >= scalar @{$labels -> [$j]} ? $numpar : scalar @{$labels -> [$j]};
 403                     }
 404                   }
 405                 }
 406               }
 407               if ( $max_hash{ $category } < $numpar ) {
 408                 $max_hash{ $category } = $numpar;
 409                 $max_hash{ 'se'.$category } = $numpar;
 410               }
 411             } else {
 412               $max_hash{ $category } = 1;
 413             }
 414           }
 415         }
 416
 417         # }}} param names
 418
 419         # --------------------------  Collect the values  ---------------------------
 420
 421         # {{{ values
 422
 423         my $pushed_rows = 0; # All rows pushed on result array in previous models
 424         my $raw_line_structure = ext::Config::Tiny -> new( );
 425         for ( my $i = 0; $i < scalar @{$self -> {'models'}}; $i++ ) { # models level
 426           my @probs = @{$self -> {'models'} -> [$i] -> {'outputs'} -> [0] -> problem_structure};
 427           my $np = scalar @probs; # #probs
 428           my $model_row = 0;
 429           # ------------  Push model, problem and sub-problem numbers  --------------
 430
 431           for( my $j = 0; $j < $np; $j++ ) {
 432             my $ns = $probs[$j]; # #subprobs
 433             for( my $k = 0; $k < $ns; $k++ ) {
 434               my $row = $pushed_rows+$model_row++;
 435               push( @{$self -> {'raw_results'} -> [$row]}, (($i+1),($j+1),($k+1)) );
 436             }
 437           }
 438
 439           # ---------------------  Loop all result categories  ----------------------
 440           my $position = 0;
 441           foreach my $category ( @{$self -> {'raw_results_header'}},'npomega' ){
 442             next if( $category eq 'model' or $category eq 'problem' or $category eq 'subproblem' );
 443             my $model_row = 0; # Need to mask previous definition of model_row
 444             my ( $accessor, $res );
 445
 446             # {{{ Get the values for the category
 447
 448             if ( $category eq 'theta' or $category eq 'omega' or $category eq 'sigma' or
 449                  $category eq 'setheta' or $category eq 'seomega' or $category eq 'sesigma' or
 450                  $category eq 'npomega' or $category eq 'eigen' ) {
 451               $accessor = $category.'s';
 452               $res = $self -> {'models'} -> [$i] -> {'outputs'} -> [0] -> $accessor;
 453
 454               if( defined $res and ref $res eq 'ARRAY' and $category ne 'npomega' ){
 455                 my $prob_num = 0;
 456                 foreach my $prob ( @{$res} ){
 457                   if( defined $prob and ref $prob eq 'ARRAY' and $prob -> [0] and ref $prob -> [0] eq 'ARRAY' ){
 458                     $raw_line_structure -> { $i+$prob_num }{ $category } = scalar @{$self -> {'raw_results'} -> [$i+$prob_num]}+$position
 459                                                                            . ",". scalar @{$prob -> [0]};
 460                     #$position += $#{$prob -> [0]};
 461                   }
 462                   $prob_num++;
 463                 }
 464               }
 465
 466             } elsif ( $category eq 'shrinkage_etas' ) {
 467               # Shrinkage does not work for subproblems right now.
 468               $res = $self -> {'models'} -> [$i] -> eta_shrinkage;
 469             } elsif ( $category eq 'shrinkage_wres' ) {
 470               # Shrinkage does not work for subproblems right now.
 471               # get ofv just to get the prob-subp structure
 472               $res = $self -> {'models'} -> [$i] -> wres_shrinkage;
 473             } else {
 474               $accessor = $category;
 475               $res = $self -> {'models'} -> [$i] -> {'outputs'} -> [0] -> $accessor;
 476             }
 477
 478             # }}} Get the values for the category
 479             my $results_file;
 480             if( $category eq 'npomega' ){
 481               $results_file = 'raw_nonp_results';
 482             } else {
 483               $results_file = 'raw_results';
 484             }
 485
 486             if( defined $res ) {
 487               for( my $j = 0; $j < $np; $j++ ) {
 488                 my $ns = $probs[$j]; # #subprobs
 489                 if( defined $res -> [$j] ) {
 490                   for( my $k = 0; $k < $ns; $k++ ) {
 491                     my $row = $pushed_rows+$model_row++;
 492                     if( ref $res -> [$j] eq 'ARRAY' ){
 493                       if( defined $res -> [$j][$k] ) {
 494                         if ( ref $res -> [$j][$k] eq 'ARRAY' ) {
 495                           push( @{$self -> {$results_file} -> [$row]}, @{$res -> [$j][$k]} );
 496                           push( @{$self -> {$results_file} -> [$row]},
 497                                 (undef) x ($max_hash{$category}- scalar @{$res -> [$j][$k]})  );
 498                         } else {
 499                           push( @{$self -> {$results_file} -> [$row]}, $res -> [$j][$k] );
 500                         }
 501                       } else {
 502                         push( @{$self -> {$results_file} -> [$row]},
 503                               (undef) x $max_hash{$category}  );
 504                       }
 505                     } else {
 506                       push( @{$self -> {$results_file} -> [$row]},
 507                             $res -> [$j] );
 508                     }
 509                   }
 510                 } else {
 511
 512                   # {{{ Push undefs for missing subprobs
 513
 514                   for( my $k = 0; $k < $ns; $k++ ) {
 515                     my $row = $pushed_rows+$model_row++;
 516                     push( @{$self -> {$results_file} -> [$row]},
 517                           (undef) x $max_hash{$category}  );
 518                   }
 519
 520                   # }}} Push undefs for missing subprobs
 521
 522                 }
 523
 524               }
 525             } else {
 526
 527               # {{{ Push undefs for missing probs/subprobs
 528
 529               for( my $j = 0; $j < $np; $j++ ) {
 530                 my $ns = $probs[$j]; # #subprobs
 531                 for( my $k = 0; $k < $ns; $k++ ) {
 532                   my $row = $pushed_rows+$model_row++;
 533                   push( @{$self -> {$results_file} -> [$row]},
 534                         (undef) x $max_hash{$category}  );
 535                 }
 536               }
 537
 538               # }}} Push undefs for missing probs/subprobs
 539
 540             }
 541           }
 542           $self -> {'models'} -> [$i] -> {'outputs'} -> [0] -> flush;
 543           map( $pushed_rows += $_, @probs );
 544         }
 545         $raw_line_structure -> write( 'raw_file_structure' );
 546         # }}} values
 547
 548         &{$self -> {'_raw_results_callback'}}( $self, \%max_hash )
 549             if ( defined $self -> {'_raw_results_callback'} );
 550
 551         # ---------------------------  Create a header  ----------------------------
 552
 553         # {{{ header
 554         my @new_header;
 555         foreach my $name ( @{$self -> {'raw_results_header'}} ) {
 556           my $success;
 557           foreach my $param ( @params ) {
 558             if ( $name eq $param ){
 559               for ( my $i = 1; $i <= $max_hash{ $name }; $i++ ) {
 560                 push ( @new_header, uc(substr($name,0,2)).$i );
 561               }
 562               $success = 1;
 563               last;
 564             } elsif ( $name eq 'se'.$param ) {
 565               my @new_names = ();
 566               for ( my $i = 1; $i <= $max_hash{ $name }; $i++ ) {
 567                 push ( @new_names, uc(substr($param,2,2)).$i );
 568               }
 569               map ( $_ = 'se'.$_, @new_names );
 570
 571               push( @new_header, @new_names );
 572               $success = 1;
 573               last;
 574             }
 575           }
 576           unless( $success ){
 577             push( @new_header, $name );
 578           }
 579         }
 580
 581         $self -> {'raw_results_header'} = \@new_header;
 582
 583         # }}} header
 584       }
 585 end prepare_raw_results
 586
 587 # }}} prepare_raw_results
 588
 589 # {{{ print_raw_results
 590
 591 start print_raw_results
 592       {
 593         if ( defined $self -> {'raw_results'} ) {
 594           my ($dir,$file) = OSspecific::absolute_path( $self -> {'directory'},
 595                                                        $self -> {'raw_results_file'} );
 596           my $append = $self -> {'raw_results_append'} ? '>>' : '>';
 597           open( RRES, $append.$dir.$file );
 598           print RRES join(',',@{$self -> {'raw_results_header'}} ),"\n"
 599             if( not $self -> {'raw_results_append'} );
 600           for ( my $i = 0; $i < scalar @{$self -> {'raw_results'}}; $i++ ) {
 601             my @result_array = @{$self -> {'raw_results'} -> [$i]};
 602             map( $_ = defined $_ ? $_ : $PsN::config -> {'_'} -> {'out_miss_data'}, @result_array );
 603             print RRES join(',',@result_array ),"\n";
 604           }
 605           close( RRES );
 606         }
 607         if( defined $self ->{'raw_nonp_results' } ) {
 608           my ($dir,$file) = OSspecific::absolute_path( $self -> {'directory'},
 609                                                        $self -> {'raw_nonp_file'} );
 610           my $append = $self -> {'raw_results_append'} ? '>>' : '>';
 611           open( RRES, $append.$dir.$file );
 612           # print RRES join(',',@{$self -> {'raw_results_header'}} ),"\n"
 613 #           if( not $self -> {'raw_results_append'} );
 614           for ( my $i = 0; $i < scalar @{$self -> {'raw_nonp_results'}}; $i++ ) {
 615             my @result_array = @{$self -> {'raw_nonp_results'} -> [$i]};
 616             map( $_ = defined $_ ? $_ : $PsN::config -> {'_'} -> {'out_miss_data'}, @result_array );
 617             print RRES join(',',@result_array ),"\n";
 618           }
 619           close( RRES );
 620         }
 621       }
 622 end print_raw_results
 623
 624 # }}} print_raw_results
 625
 626 # {{{ create_sub_dir
 627
 628 start create_sub_dir
 629     {
 630         my $file;
 631 #       ($tmp_dir, $file) = OSspecific::absolute_path( './' .
 632         ($tmp_dir, $file) = OSspecific::absolute_path( $self -> {'directory'}.'/' .
 633                                                        $subDir, '');
 634
 635         unless( -e $tmp_dir ){
 636             mkdir( $tmp_dir );
 637         }
 638     }
 639 end create_sub_dir
 640
 641 # }}}
 642
 643 # {{{ copy_model_and_input
 644
 645 start copy_model_and_input
 646     {
 647       # Copy data files to temporary dirctory and set short version names
 648       my $i = 1;
 649       my @new_data_names;
 650
 651       if( defined $model -> datas ){
 652         foreach my $data ( @{$model -> datas} ) {
 653           my $filename = $data -> filename;
 654
 655           # We always copy the datafile in case a crash corrupted the
 656           # copy.
 657
 658 #         cp( $data -> full_name, $filename );
 659
 660           push( @new_data_names, $filename );
 661           $i++;
 662         }
 663       } else {
 664         debug -> die( message => 'No datafiles set in modelfile.' );
 665       }
 666
 667       # Copy extra data files to temporary dirctory and set short version names
 668       my @new_extra_data_names;
 669       my @problems = @{$model -> problems};
 670       for ( my $i = 1; $i <= $#problems + 1; $i++ ) {
 671         my $extra_data = $problems[$i-1] -> extra_data;
 672         if ( defined $extra_data ) {
 673           my $filename = $extra_data -> filename;
 674           if ( not -e $filename ) {
 675 #           cp( $extra_data -> full_name, $filename );
 676           }
 677           push( @new_extra_data_names, $filename );
 678         }
 679       }
 680
 681       # Copy extra files needed to run NONMEM
 682       if( defined $model -> extra_files ){
 683         foreach my $x_file (@{$model -> extra_files}){
 684           debug -> warn( level   => 2,
 685                          message => "Copying extra file: $x_file" );
 686           my ( $dir, $filename ) = OSspecific::absolute_path( $model -> directory,
 687                                                               $x_file );
 688           cp( $dir.$filename, $filename );
 689         }
 690       }
 691
 692       # Set the table names to a short version
 693       my @new_table_names = ();
 694       my @table_names = @{$model -> table_names( ignore_missing_files => 1 )};
 695       # Loop the problems
 696       for ( my $i = 0; $i <= $#table_names; $i++ ) {
 697         my @new_arr;
 698         # Loop the table files within each problem
 699         for ( my $j = 0; $j < scalar @{$table_names[$i]}; $j++ ) {
 700           my ( $dir, $filename ) = OSspecific::absolute_path( '.', $table_names[$i][$j] );
 701           push( @new_arr,  $filename );
 702         }
 703         push( @new_table_names, \@new_arr );
 704       }
 705
 706       # Copy the model object. Set the new (shorter) data file names.
 707       # There's no need to physically copy these here since we took care of that above.
 708       $final_model = $model -> copy( filename              => 'psn.mod',
 709                                      data_file_names       => \@new_data_names,
 710                                      extra_data_file_names => \@new_extra_data_names,
 711                                      copy_data             => 1 );
 712       $final_model -> register_in_database;
 713       $model -> flush_data;
 714
 715       $final_model -> table_names( new_names            => \@new_table_names,
 716                                    ignore_missing_files => 1 );
 717 #      $final_model -> clean_extra_data_code;
 718       $final_model -> drop_dropped if ( $self -> {'drop_dropped'} );
 719       $final_model -> wrap_data if ( $self -> {'wrap_data'} );
 720       $final_model -> add_extra_data_code;
 721       $final_model -> write_get_subs;
 722       $final_model -> write_readers;
 723       $final_model -> _write( filename   => 'psn.mod',
 724                               write_data => 1 );  #Kolla denna, den funkar inte utan wrap!!
 725       $final_model -> flush_data;
 726       $final_model -> store_inits;
 727
 728       $self -> {'orig_filename'} = $model -> filename;
 729       $maxeval = $model -> maxeval -> [0][0];
 730       if ( $maxeval > 9999 ) {
 731         $final_model -> maxeval( new_values => 9999,
 732                                  problem_numbers => [1] );
 733       }
 734     }
 735 end copy_model_and_input
 736
 737 # }}}
 738
 739 # {{{ copy_model_and_output
 740
 741 start copy_model_and_output
 742 {
 743   my $outfilename = $model -> outputs -> [0] -> full_name;
 744
 745   cp( "psn-$use_run.lst", 'psn.lst' ) if ( $use_run ne '' );
 746   cp( "psn.lst", $outfilename );
 747   if( defined $model -> table_names ){
 748     foreach my $table_files( @{$model -> table_names} ){
 749       foreach my $table_file( @{$table_files} ){
 750         cp( $table_file.'-'.$use_run, $table_file );
 751       }
 752     }
 753   }
 754   my $final_output = output -> new( filename => $outfilename,
 755                                     model_id => $model -> model_id );
 756   $final_output -> register_in_database( model_id => $model -> model_id,
 757                                          force    => 1 ); # If we are here, the model has been run
 758   $final_model -> outputs( [$final_output] );
 759
 760   if( defined $model -> table_names ){
 761     foreach my $table_files( @{$model -> table_names} ){
 762       foreach my $table_file( @{$table_files} ){
 763         my ( $dir, $filename ) = OSspecific::absolute_path($model -> directory,
 764                                                            $table_file );
 765         cp( $filename, $dir.$filename  );
 766       }
 767     }
 768     $final_model -> table_names( new_names => $model -> table_names );
 769   }
 770
 771   my @problems = @{$model -> problems};
 772   my @final_problems = @{$final_model -> problems};
 773   for( my $i = 0; $i <= $#problems; $i++ ) {
 774     if( $problems[$i-1] -> shrinkage_module -> enabled ) {
 775       my ( $dir, $eta_filename ) =
 776           OSspecific::absolute_path($model -> directory,
 777                                     $problems[$i] -> shrinkage_module -> eta_tablename );
 778       cp( $final_problems[$i] -> shrinkage_module -> eta_tablename, $dir.$eta_filename  );
 779       my ( $dir, $wres_filename ) =
 780           OSspecific::absolute_path($model -> directory,
 781                                     $problems[$i] -> shrinkage_module -> wres_tablename );
 782       cp( $final_problems[$i] -> shrinkage_module -> wres_tablename, $dir.$wres_filename  );
 783     }
 784   }
 785
 786   unlink 'nonmem', 'nonmem6', 'nonmem5','nonmem.exe', 'nonmem6_adaptive', 'nonmem5_adaptive';
 787   if( $self -> {'clean'} >= 1){
 788     unlink 'nonmem', 'nonmem'.$self -> {'nm_version'}, 'nonmem.exe','FDATA', 'FREPORT', 'FSUBS', 'FSUBS.f', 'FSUBS.for', 'LINK.LNK', 'FSTREAM';
 789     foreach my $data ( @{$model -> datas} ){
 790       unlink $data -> filename;
 791     }
 792     for ( my $i = 1; $i <= $self -> {'retries'}; $i++ ) {
 793       unlink( 'FCON.'.$i );
 794       unlink( "psn-$i.lst" );
 795       unlink( "psn-$i.mod" );
 796       unlink( 'compilation_output.txt.'.$i );
 797     }
 798   }
 799
 800   if ( $self -> {'clean'} >= 2 ) {
 801     unlink 'FCON', 'PRDERR';
 802   }
 803
 804   if ( $self -> {'clean'} >= 3 ) {
 805     # Do nothing. "run_nonmem" will remove entire work directory
 806     # before returning.
 807   } else {
 808     system('tar cz --remove-files -f nonmem_files.tgz *')
 809         if ( $self -> {'compress'} and $Config{osname} ne 'MSWin32' );
 810     system('compact /c /s /q > NUL')
 811         if ( $self -> {'compress'} and $Config{osname} eq 'MSWin32' );
 812   }
 813 }
 814 end copy_model_and_output
 815
 816 # }}}
 817
 818 # {{{ _twoDimToString
 819
 820 start _twoDimToString
 821 {
 822     $string .= "[" ;
 823     for(my $j = 0; $j <= $#array; $j++ ){
 824         $string .= "[" ;
 825
 826         my $thetas = $array[$j];
 827
 828         for(my $i = 0;$i < $#{$thetas}; $i++ ){
 829             $string .= $thetas -> [$i] . ",";
 830         }
 831         $string .= $thetas -> [$#{$thetas}] . "]";
 832         if( $j < $#array ){
 833             $string .= ",";
 834         }
 835     }
 836     $string .= "]";
 837 }
 838 end _twoDimToString
 839
 840 # }}}
 841
 842 # {{{ _printGridScript
 843
 844 start _printGridScript
 845 {
 846     my $mainExecName = 'batch' . $batchStart . '.pl';
 847
 848     open( my $gridScript , ">>", $self -> directory . "/$mainExecName" );
 849     print $gridScript "#!/usr/bin/perl\n";
 850     print $gridScript 'use FindBin qw($Bin);', "\n";
 851     print $gridScript "use fcon;\n\n";
 852
 853     print $gridScript "my \$fcon = fcon -> new('filename' => 'FCON');\n";
 854     print $gridScript "\$fcon -> run('retries' => 5, 'batchStart' => $batchStart, 'batchSize' => $batchSize,\n";
 855     print $gridScript "'fixed_thetas' => [";
 856     $, = ",";
 857     print $gridScript @theta_fixed;
 858     $, = undef;
 859     print $gridScript "], 'fixed_omegas' => [";
 860     $, = ",";
 861     print $gridScript @omega_fixed;
 862     $, = undef;
 863     print $gridScript "], 'fixed_sigmas' => [";
 864     $, = ",";
 865     print $gridScript @sigma_fixed;
 866     $, = undef;
 867     print $gridScript "], 'retries' => [";
 868     $, = ",";
 869     print $gridScript @retries;
 870     $, = undef;
 871     print $gridScript "], 'picky' => [";
 872     $, = ",";
 873     print $gridScript @picky;
 874     $, = undef;
 875     print $gridScript "] )\n";
 876
 877     close( $gridScript );
 878 }
 879 end _printGridScript
 880
 881 # }}}
 882
 883 # {{{ set_msfo_to_msfi
 884
 885 start set_msfo_to_msfi
 886     {
 887       my ( $msfo_ref, $junk ) = $final_model ->
 888           _get_option_val_pos( problem_numbers => [0],
 889                                name            => 'MSFO',
 890                                record_name     => 'estimation' );
 891
 892       unless ( defined $msfo_ref ) {
 893         debug -> die( message => "handle_maxevals set, but no MSFO option " .
 894                       "defined in model file " . $final_model -> filename );
 895       }
 896
 897       my $msfo = @{$msfo_ref}->[0][0];
 898       $msfo =~ /(\d+)\Z/;
 899       my $msfi = $msfo;
 900       if ( defined ( $1 ) ) {
 901         cp( 'psn.lst', 'psn.lst.msfo'.($1) );
 902         cp( $final_model -> filename, $final_model -> filename.'.msfo'.($1) );
 903       }
 904       my $new_num = defined ( $1 ) ? $1+1 : 1;
 905       $msfo =~ s/\d*(\Z)/$new_num$1/;
 906
 907       $final_model->_option_val_pos( name             => 'MSFO',
 908                                      record_name      => 'estimation',
 909                                      problem_numbers  => [[0]],
 910                                      new_values       => [[$msfo]]);
 911
 912       $final_model->set_records(type=>'msfi',
 913                                 record_strings => ["$msfo"]); # $msfi or ????
 914     }
 915 end set_msfo_to_msfi
 916
 917 # }}}
 918
 919 # {{{ reset_msfo
 920
 921 start reset_msfo
 922     {
 923       my @data_ref = @{$final_model -> record( record_name => 'msfi' )};
 924       # Check if there is a msfi record and then delete it
 925       if (scalar(@data_ref)!=0) {
 926         $final_model->remove_records(type=>'msfi');
 927
 928         # Set the intial values + boundaries to the first  values (update theta, omega, sigma)
 929
 930         my @old_problems = @{$basic_model -> problems};
 931         my @new_problems = @{$final_model -> problems};
 932         for ( my $i=0; $i <= $#old_problems; $i++ ) {
 933           foreach my $param ( 'thetas', 'omegas', 'sigmas' ) {
 934             $new_problems[$i] -> $param( Storable::dclone( $old_problems[$i] -> $param ) );
 935           }
 936         }
 937
 938         $model_modified = 1;
 939       }
 940
 941       # Update the msfo file if defined
 942       my ( $msfo_ref, $junk ) = $final_model ->
 943           _get_option_val_pos( problem_numbers => [0],
 944                                name            => 'MSFO',
 945                                record_name     => 'estimation' );
 946
 947       if ( defined $msfo_ref ) {
 948         my $msfo = @{$msfo_ref}->[0][0];
 949         $msfo =~ /(\d+)\Z/;
 950
 951         my $new_num = defined ( $1 ) ? $1+1 : 1;
 952         $msfo =~ s/\d*(\Z)/$new_num$1/;
 953
 954         $final_model->_option_val_pos( name             => 'MSFO',
 955                                        record_name      => 'estimation',
 956                                        problem_numbers  => [[0]],
 957                                        new_values       => [[$msfo]]);
 958       }
 959     }
 960 end reset_msfo
 961
 962 # }}}
 963
 964 # {{{ cut_thetas
 965
 966 start cut_thetas
 967     {
 968       $final_model -> update_inits( from_output => $output_file,
 969                                     update_omegas => 1,
 970                                     update_sigmas => 1,
 971                                     update_thetas => 1);
 972
 973       foreach my $th_num ( @cutoff_thetas ) {
 974         my $init_val = $final_model ->
 975             initial_values( parameter_type    => 'theta',
 976                             parameter_numbers => [[$th_num]])->[0][0];
 977         if (abs($init_val)<=$self->{'cutoff'}) {
 978           $final_model->initial_values(parameter_type => 'theta',
 979                                        parameter_numbers => [[$th_num]],
 980                                        new_values =>[[0]]);
 981           $final_model->fixed(parameter_type => 'theta',
 982                               parameter_numbers => [[$th_num]],
 983                               new_values => [[1]] );
 984         }
 985       }
 986     }
 987 end cut_thetas
 988
 989 # }}}
 990
 991 # {{{ ask_user
 992
 993 start ask_user
 994     {
 995       # rand should not be used here. find some other way to create unique file names.
 996       my $num = rand;
 997       open( TMP, ">/tmp/$num" );
 998       print TMP "START MODEL FILE NAME\n";
 999       print TMP $basic_model -> filename,"\n";
1000       print TMP "END MODEL FILE NAME\n";
1001       foreach my $prob ( @reruns ) {
1002         my @theta_labels = @{$final_model -> labels( parameter_type => 'theta' )};
1003         my @omega_labels = @{$final_model -> labels( parameter_type => 'omega' )};
1004         my @sigma_labels = @{$final_model -> labels( parameter_type => 'sigma' )};
1005         my @theta_inits = @{$final_model -> initial_values( parameter_type => 'theta' )};
1006         my @omega_inits = @{$final_model -> initial_values( parameter_type => 'omega' )};
1007         my @sigma_inits = @{$final_model -> initial_values( parameter_type => 'sigma' )};
1008         print TMP "START PROBLEM NUMBER\n";
1009         print TMP $prob,"\n";
1010         print TMP "END PROBLEM NUMBER\n";
1011         print TMP "START MINIMIZATION_MESSAGE\n";
1012         print TMP $minimization_message -> [$prob-1][0],"\n";
1013         print TMP "END MINIMIZATION_MESSAGE\n";
1014         print TMP "START FINAL GRADIENT\n";
1015         print TMP join( " ",@{$output_file -> final_gradients -> [$prob-1][0]}),"\n";
1016         print TMP "END FINAL GRADIENT\n";
1017         print TMP "START OFV\n";
1018         print TMP $output_file -> ofv -> [$prob-1][0],"\n";
1019         print TMP "END OFV\n";
1020         print TMP "START INITIAL VALUES THETA\n";
1021         print TMP join(" ", @{$theta_inits[$prob-1]}),"\n";
1022         print TMP "END INITIAL VALUES THETA\n";
1023         print TMP "START INITIAL VALUES OMEGA\n";
1024         print TMP join(" ", @{$omega_inits[$prob-1]}),"\n";
1025         print TMP "END INITIAL VALUES OMEGA\n";
1026         print TMP "START INITIAL VALUES SIGMA\n";
1027         print TMP join(" ", @{$sigma_inits[$prob-1]}),"\n";
1028         print TMP "END INITIAL VALUES SIGMA\n";
1029         print TMP "START LABELS\n";
1030         print TMP join(" ", (@{$theta_labels[$prob-1]},@{$omega_labels[$prob-1]},
1031                              @{$sigma_labels[$prob-1]})),"\n";
1032         print TMP "END LABELS\n";
1033       }
1034       close( TMP );
1035       my $out = readpipe( "/users/lasse/PsN/Diagrams/test/scm_comm.pl $num ".$output_file->filename );
1036       my @out_per_prob = split("\n",$out);
1037       foreach my $prob ( @reruns ) {
1038         my ( $choice, $rest ) = split( ' ', shift( @out_per_prob ), 2 );
1039         if ( $choice == 0 ) {
1040           $retries = $tries + $rest;
1041         } elsif ( $choice == 1 ) {
1042           my ($theta_str,$omega_str,$sigma_str) = split(':',$rest);
1043           print "thstr $theta_str\n";
1044           print "omstr $omega_str\n";
1045           print "sistr $sigma_str\n";
1046           my @thetas = split( ' ', $theta_str );
1047           print "$prob: @thetas\n";
1048           $final_model -> initial_values( parameter_type  => 'theta',
1049                                           problem_numbers => [$prob],
1050                                           new_values      => [\@thetas],
1051                                           add_if_absent   => 0 );
1052           $retries = $tries+1;
1053         } else {
1054           last RETRY;
1055         }
1056         $final_model -> _write;
1057       }
1058       $return_value = $retries;
1059     }
1060 end ask_user
1061
1062 # }}}
1063
1064 # {{{ lsf_submit
1065
1066 start lsf_submit
1067     {
1068
1069       # This method will submit the nonmem.pm file as a script to the
1070       # LSF system. The nonmem.pm file will use the
1071       # "no_remote_compile" and "no_remote_execution" arguments what
1072       # to do on the remote host. It can either compile or execute or
1073       # both. If only execution is done remote, compilation will
1074       # already be done by a call to nonmem -> compile in
1075       # "run_nonmem".
1076
1077       my $fsubs = join( ',' , @{$model -> subroutine_files} );
1078
1079       # Check for vital lsf options.
1080       unless( $self -> {'lsf_queue'} ){
1081         if( $PsN::config -> {'_'} -> {'lsf_queue'} ){
1082           $self -> {'lsf_queue'} = $PsN::config -> {'_'} -> {'lsf_queue'};
1083         } else {
1084           'debug' -> die( message => 'No queue specified for lsf run' );
1085         }
1086       }
1087
1088       foreach my $lsf_option ( 'lsf_project_name', 'lsf_job_name', 'lsf_resources', 'lsf_ttl', 'lsf_options' ){
1089         unless( $self -> {$lsf_option} ){
1090           if( $PsN::config -> {'_'} -> {$lsf_option} ){
1091             $self -> {$lsf_option} = $PsN::config -> {'_'} -> {$lsf_option};
1092           }
1093         }
1094       }
1095
1096       run3( "bsub -e stderr -o stdout -K " .
1097             "-q " . $self -> {'lsf_queue'} .
1098             ($self -> {'lsf_project_name'} ? " -P " . $self -> {'lsf_project_name'} : ' ') .
1099             ($self -> {'lsf_job_name'} ? " -J " . $self -> {'lsf_job_name'} : ' ') .
1100             ($self -> {'lsf_ttl'} ? " -c " . $self -> {'lsf_ttl'} : ' ') .
1101             ($self -> {'lsf_resources'} ? " -R " . $self -> {'lsf_resources'} : ' ') .
1102             $self -> {'lsf_options'} .
1103             ($PsN::config -> {'_'} -> {'remote_perl'} ? ' ' . $PsN::config -> {'_'} -> {'remote_perl'} : ' perl ') . " -I" .
1104             $PsN::lib_dir ."/../ " .
1105             $PsN::lib_dir . "/nonmem.pm" .
1106             " psn.mod psn.lst " .
1107             $self -> {'nice'} . " ".
1108             $nm_version . " " .
1109             $self -> {'no_remote_compile'} . " " .
1110             $self -> {'no_remote_execution'} . " " .
1111             $fsubs . " " .
1112             $self -> {'nm_directory'} );
1113
1114     }
1115 end lsf_submit
1116
1117 # }}}
1118
1119 # {{{ run_nonmem
1120
1121 start run_nonmem
1122     {
1123       my $maxeval;
1124       my $evals = 0;
1125       my @run_results;
1126       ( $final_model, $maxeval ) = $self -> copy_model_and_input( model => $model, source => '../' );
1127
1128       # {{{ Experimental code for MSFO-resume
1129
1130 #       $final_model = $self -> copy_model_and_input( model => $model, source => $base_dir );
1131
1132 #       my ($msfo_ref, $junk) = $final_model -> _get_option_val_pos( problem_numbers => [0],
1133 #                                                                  name            => 'MSFO',
1134 #                                                                  record_name     => 'estimation' );
1135 #       if( defined $msfo_ref and -e @{$msfo_ref}->[0][0] ){
1136 #       $self -> set_msfo_to_msfi( final_model => $final_model );
1137
1138 #       $final_model->remove_records(type=>'theta');
1139 #       $final_model->remove_records(type=>'omega');
1140 #       $final_model->remove_records(type=>'sigma');
1141
1142 #       $final_model -> _write;
1143 #       }
1144
1145 #       $final_model = $self -> copy_model_and_input( model => $model, source => $base_dir );
1146
1147 # }}}
1148
1149       my $marked_for_rerun = 1;
1150       my $modelfile_tainted = 1;
1151       my $nonmem;
1152
1153       # Here we check wheter compilation and/or exection will be done
1154       # localy or remote. (Remote compilation is only supported on LSF
1155       # so far). So, unless we are doing everything remotely, we will
1156       # need a "nonmem" object. Which is created here.
1157
1158       if( not $self -> {'run_on_lsf'} or
1159           ( $self -> {'run_on_lsf'} and
1160             ( $self -> {'no_remote_compile'} or
1161               $self -> {'no_remote_execution'}))){
1162         $nonmem = nonmem -> new( adaptive => $self -> {'adaptive'},
1163                                  modelfile => 'psn.mod',
1164                                  version => $nm_version,
1165                                  nm_directory => $self -> {'nm_directory'},
1166                                  nice => $self -> {'nice'},
1167                                  show_version => not $self -> {'run_on_lsf'});
1168
1169         $nonmem -> fsubs( $final_model -> subroutine_files );
1170       }
1171
1172       my $fprob = $final_model -> problems;
1173       if( defined $fprob ) {
1174         my $nprob = scalar @{$fprob};
1175         for( my $i = 0; $i < $nprob; $i++ ) {
1176           my ( $model_msfos, $junk ) =
1177               $final_model -> _option_val_pos( name             => 'MSFO',
1178                                                record_name      => 'estimation',
1179                                                problem_numbers  => [$i] );
1180
1181           if( not scalar @{$model_msfos} > 0 ) {
1182             $final_model -> _option_val_pos( name             => 'MSFO',
1183                                              record_name      => 'estimation',
1184                                              problem_numbers  => [$i],
1185                                              new_values       => [["msfo$i"]]);
1186           }
1187         }
1188       }
1189       $final_model->_write();
1190
1191       # -------------- Notes about automatic pertubation and retries -----------------
1192
1193       # Automatic pertubation of initial estimates are usefull for two
1194       # purposes. On reason is when nonmem failed to produce a
1195       # successful minimization we can try to direct the search by
1196       # selection other estimates. But it is also possible to get a
1197       # successful minimization but in a local minima, in this case
1198       # we have now way of knowing that, but if we set new intial
1199       # estimates we can get lucky and find a global minima. Two
1200       # modelfit members govern the pertubation process "retries"
1201       # which is a number giving the maximum number of retries when
1202       # nonmem failes and "auto_perturbation" which is the number of
1203       # runs we want to do to get a global minima. So if
1204       # auto_pertubation is 3 and retries is 5 we will stop after 3
1205       # successful runs but continue until 5 if we nonmem fails to
1206       # minimize successful. The reason we don't do retries for each
1207       # auto_pertubation try is because the pertubation takes place in
1208       # the same space of possible estimates. In other words, if 5
1209       # retry pertubations of the original model fails then it will
1210       # likely fail for a retry of a auto pertubated model.
1211
1212       # It is important to set $marked_for_rerun to 1 if $tries is incremented!
1213       # Otherwise $tries can be incremented twice for one run. The
1214       # oposite is not true however, for instance a reset of maxevals
1215       # is not a retry but sets $marked_for_rerun to 1.
1216
1217       while( $marked_for_rerun and $tries < $retries ){
1218
1219         # {{{ Run nonmem and get minimization messages
1220
1221         # We do not expect any values of rerun lower than 1 here. (a bug otherwise...)
1222         if( not -e 'psn-' . ( $tries + 1 ) . '.lst' or $self -> {'rerun'} >= 2 ){
1223           if( $modelfile_tainted ){ # We need to recompile
1224             if( $self -> {'run_on_nordugrid'} or not $self -> {'run_on_lsf'} or
1225                 ( $self -> {'run_on_lsf'}       and $self -> {'no_remote_compile'} ) ){
1226
1227               # {{{ local compile requested
1228
1229               unless( $nonmem -> compile() ){
1230                 debug -> die( message => "NONMEM compilation failed:\n" . $nonmem -> error_message );
1231               }
1232
1233               unless ( -e 'FCON.orig' ) {
1234                 unless( $self -> {'clean'} ){
1235
1236                   # Instead of cleaning its better to never create
1237                   # garbage :) However, if FCON pertubation becomes
1238                   # reality, we will need FCON.orig.
1239
1240                   cp( 'FCON', 'FCON.orig' );
1241                 }
1242               }
1243
1244               if( $nonmem -> nmtran_message =~ /WARNING/s ){
1245                 debug -> warn( level => 1,
1246                                message => "NONMEM Warning: " . $nonmem -> nmtran_message );
1247               }
1248
1249               $modelfile_tainted = 0;
1250
1251               # }}}
1252
1253             } elsif ( $self -> {'run_on_lsf'} and $self -> {'no_remote_execution'} ) {
1254
1255
1256               # lsf_submit will call the "nonmem" module that will
1257               # figure out that we want to compile remotely. If we want
1258               # to compile AND run remotely that will be done in the
1259               # call to "lsf_submit" below.
1260
1261               $self -> lsf_submit( model => $final_model.
1262                                    nm_version => $nm_version );
1263
1264             }
1265           }
1266
1267           if( not $self -> {'run_on_nordugrid'} and not $self -> {'run_on_lsf'} or
1268               ( $self -> {'run_on_lsf'} and $self -> {'no_remote_execution'} ) ){
1269
1270             # Normal local execution
1271             $nonmem -> execute();
1272
1273           } elsif( $self -> {'run_on_lsf'} ) {
1274
1275             # lsf_submit will call the "nonmem" module that will figure
1276             # out that we want to run remotely. If we are also compiling
1277             # remotely, it will be done from here as well.
1278
1279             $self -> lsf_submit( model => $final_model,
1280                                  nm_version => $nm_version );
1281
1282           } elsif ( $self -> {'run_on_nordugrid'} ){
1283
1284             # {{{ submit and execute using nordugrid
1285
1286             my ( $ng_failed, $ng_tries ) = ( 1, 1 );
1287             while( $ng_failed and $ng_tries <= 5 ) {
1288               my $jobidfile = $self -> ng_submit( model => $final_model,
1289                                                   run   => $run_no);
1290               my ( $job_files, $job_number ) = $self -> ng_monitor( jobidfile => $jobidfile,
1291                                                                     run       => $run_no );
1292               $self -> ng_retrieve( job_number => $job_number,
1293                                     job_files => $job_files,
1294                                     jobidfile => $jobidfile,
1295                                     run       => $run_no );
1296               if ( -e $self -> {'directory'}.'NM_run'.($run_no+1).
1297                    "/$job_number/psn.log/failed" ) {
1298                 $ng_failed = 1;
1299                 $ng_tries++;
1300                 print "Job number ",($run_no+1)," with id $job_number failed at nordugrid\n",
1301                 "... rerunning\n";
1302               } else {
1303                 $ng_failed = 0;
1304               }
1305             }
1306
1307             # }}}
1308
1309           }
1310
1311           # We need the trail of files to select the most appropriate at the end
1312           # (see copy_model_and_output)
1313           cp( 'FCON', 'FCON.'.($tries+1) );
1314           cp( 'psn.lst', 'psn-'.($tries+1).'.lst' );
1315           cp( 'psn.mod', 'psn-'.($tries+1).'.mod' );
1316           if( defined $model -> table_names ){
1317             foreach my $table_files( @{$model -> table_names} ){
1318               foreach my $table_file( @{$table_files} ){
1319                 cp( $table_file, $table_file.'-'.($tries+1) );
1320               }
1321             }
1322           }
1323         } else {
1324
1325           # This is rerun==1, i.e. re-evaluate the stuff that has been
1326           # run and (possibly) run extra runs to fix any problems left.
1327           # In this "else" statement it is true that psn-X.lst exists
1328           # and we copy it to psn.lst to make it the current version.
1329           cp( 'psn-'.($tries+1).'lst', 'psn.lst' );
1330
1331         }
1332
1333         if ( $PsN::config -> {'_'} -> {'use_database'} ) {
1334           # Save the temporary run in the database
1335 #         for ( my $i = 1; $i <= $self -> {'retries'}; $i++ ) {
1336 #           if ( -e "psn-$i.mod" ) {
1337           my $tmp_model = model -> new( filename => 'psn-'.($tries+1).'.mod',
1338                                         datas    => [] ); # Avoid load of data
1339           my $model_id = $tmp_model -> register_in_database( force => 1 );
1340           if ( defined $model_id and defined $tmp_model -> outputs -> [0] ) {
1341             $tmp_model -> outputs -> [0] -> register_in_database( model_id => $model_id,
1342                                                                   force    => 1 );
1343             $self -> register_tm_relation( model_ids => [$model_id] );
1344           }
1345         }
1346 #         }
1347 #       }
1348         cp( 'compilation_output.txt', 'compilation_output.txt.'.($tries+1) );
1349
1350         $marked_for_rerun = 0;
1351
1352         my ( $minimization_successful, $significant_digits, $minimization_message );
1353         my ( $output_file );
1354
1355         if ( -e 'psn.lst' ) {
1356           $output_file = output -> new ( filename => 'psn.lst',
1357                                          model_id => $self -> {'model_id'} );
1358           $minimization_successful       = $output_file -> minimization_successful;
1359           $significant_digits            = $output_file -> significant_digits;
1360           $evals[$tries]              += $output_file -> feval -> [0][0];
1361
1362           # {{{ log the stats of this run
1363
1364           foreach my $category ( 'minimization_successful', 'covariance_step_successful',
1365                                  'covariance_step_warnings', 'estimate_near_boundary',
1366                                  'significant_digits', 'ofv' ){
1367             my $res = $output_file -> $category;
1368             $run_results[$tries] -> {$category} = defined $res ? $res -> [0][0] : undef;
1369           }
1370           $run_results[$tries] -> {'pass_picky'} = 0;
1371
1372           # }}}
1373
1374           unless( defined $minimization_successful ) {
1375             debug -> die( message => "No minimization status found in " . $output_file ->filename );
1376           }
1377
1378           $minimization_message = $output_file -> minimization_message;
1379           my $ofv = $output_file -> ofv;
1380           my $problems = $final_model -> problems;
1381
1382           # We should avoid testing this when we are doing a simulation
1383           unless( $#{$problems} == $#{$minimization_message} or
1384                   $#{$problems} == $#{$ofv} ){
1385             $marked_for_rerun = 1;
1386             debug -> warn( level => 1,
1387                            message => "The number of problems (".
1388                            scalar @{$problems}.
1389                            ") does not match the number of minimization messages (".
1390                            scalar @{$minimization_message}.
1391                            ").\nNONMEM was likely killed, restarting." );
1392             $tries++; # We need this here to avoid a possible infinite loop.
1393           }
1394
1395         } else {
1396           debug -> die( message => "NONMEM ran, but gave no lst-file" );
1397         }
1398
1399         # }}}
1400
1401         if ( not $marked_for_rerun and $handle_maxevals and $maxeval > $evals[$tries] ) {
1402
1403           # {{{ Check if maxevals is reached and copy msfo to msfi
1404
1405           for ( @{$minimization_message -> [0][0]} ) {
1406             if ( /\s*MAX. NO. OF FUNCTION EVALUATIONS EXCEEDED\s*/) {
1407
1408               if( $self -> {'use_implicit_msfo'} ) {
1409                 $self -> set_msfo_to_msfi( final_model => $final_model );
1410               }
1411
1412               $final_model->remove_records(type=>'theta');
1413               $final_model->remove_records(type=>'omega');
1414               $final_model->remove_records(type=>'sigma');
1415
1416               $final_model -> _write;
1417               $modelfile_tainted = 1;
1418               $marked_for_rerun = 1;
1419               last;
1420             }
1421           }
1422
1423           # }}}
1424
1425         }
1426
1427         if ( not $marked_for_rerun and $handle_rounding_errors || $handle_hessian_npd) {
1428
1429           # {{{ Check for rounding errors and/or hessian_npd messages
1430
1431           my $round_rerun;
1432           my $hessian_rerun;
1433           for ( @{$minimization_message -> [0][0]} ) {
1434             $round_rerun = 1 and last if ( /\s*DUE TO ROUNDING ERRORS\s*/);
1435             $hessian_rerun = 1 and last if ( /\s*NUMERICAL HESSIAN OF OBJ. FUNC. FOR COMPUTING CONDITIONAL ESTIMATE IS NON POSITIVE DEFINITE\s*/);
1436           }
1437
1438           if ( ($round_rerun && $handle_rounding_errors) or ($hessian_rerun && $handle_hessian_npd)) {
1439
1440             if( $self -> {'use_implicit_msfo'} ) {
1441               $self -> reset_msfo( basic_model => $model,
1442                                    final_model => $final_model );
1443             }
1444
1445             $self -> cut_thetas( final_model => $final_model,
1446                                  cutoff_thetas => \@cutoff_thetas,
1447                                  output_file => $output_file );
1448
1449             $final_model -> _write;
1450             $modelfile_tainted = 1;
1451             print "ROUND\n";
1452             $marked_for_rerun = 1;
1453
1454             $tries ++; # This is a precaution, handle_rounding and handle_hessian should have
1455                        # their own termination checks
1456           }
1457
1458           # }}}
1459
1460         }
1461
1462         if ( not $marked_for_rerun and $tweak_inits ) {
1463
1464           # {{{ Check for failed problems and possibly check for picky errors.
1465
1466           my @reruns;
1467           my @problems = @{$final_model -> problems};
1468           for ( my $problem = 1; $problem <= scalar @problems; $problem++ ) {
1469             unless( $final_model -> is_simulation( problem_number => $problem ) ){
1470               if ( $minimization_successful -> [$problem-1][0] ) {
1471                 if ( $picky ) {
1472                   $run_results[$tries] -> {'pass_picky'} = 1;
1473                   for ( @{$minimization_message -> [$problem-1][0]} ) {
1474                     if ( /0COVARIANCE STEP ABORTED/ or
1475                          /0PROGRAM TERMINATED BY OBJ/ or
1476                          /0ESTIMATE OF THETA IS NEAR THE BOUNDARY AND/ or
1477                          /0PARAMETER ESTIMATE IS NEAR ITS BOUNDARY/ or
1478                          /0R MATRIX ALGORITHMICALLY SINGULAR/ or
1479                          /0S MATRIX ALGORITHMICALLY SINGULAR/ ) {
1480                       push( @reruns, $problem );
1481                       $run_results[$tries] -> {'pass_picky'} = 0;
1482                       last;
1483                     }
1484                   }
1485                 }
1486               } else {
1487                 if ( not ( $significant_digits_rerun and $significant_digits -> [$problem-1][0] > $significant_digits_rerun ) ) {
1488                   push( @reruns, $problem );
1489                 }
1490               }
1491             }
1492           }
1493
1494           if( $tries < ($retries -1) and scalar @reruns > 0 ) {
1495             $marked_for_rerun = 1;
1496             $tries ++;
1497
1498             if( $tries >= $self -> {'min_retries'} and $self -> {'verbose'} ){
1499               my $message = "R:".($run_no+1).":". ($tries+1) . " ";
1500               ui -> print( category => 'all',  message  => $message,
1501                            newline => 0);
1502             }
1503
1504             if( $self -> {'ask_if_fail'} ) {
1505               $retries = $self -> ask_user( basic_model => $model,
1506                                             final_model => $final_model,
1507                                             reruns => \@reruns,
1508                                             minimization_message => $minimization_message,
1509                                             output_file => $output_file,
1510                                             retries => $retries,
1511                                             tries => $tries );
1512               $final_model->_write;
1513               $modelfile_tainted = 1;
1514             } else {
1515
1516               # This code must be adjusted for multiple problems!!
1517               my $degree = 0.1*$tries;
1518               if( $self -> {'use_implicit_msfo'} and
1519                   $self -> reset_msfo( basic_model => $model,
1520                                        final_model => $final_model ) ){
1521
1522                 foreach my $prob ( @reruns ) {
1523                   $problems[$prob-1] -> set_random_inits ( degree => $degree );
1524                 }
1525
1526                 $final_model->_write;
1527                 $modelfile_tainted = 1;
1528
1529               } else {
1530                 foreach my $prob ( @reruns ) {
1531                   $problems[$prob-1] -> set_random_inits ( degree => $degree );
1532                 }
1533
1534                 $final_model->_write;
1535                 $modelfile_tainted = 1;
1536                 # The fcon code does not parse the dofetilide model correctly
1537 #               my $fcon = fcon -> new( filename => 'FCON.orig' );
1538 #               $fcon -> parse;
1539 #               $fcon -> pertubate_all( fixed_thetas => $final_model -> fixed( 'parameter_type' => 'theta' ),
1540 #                                       fixed_omegas => $final_model -> fixed( 'parameter_type' => 'omega' ),
1541 #                                       fixed_sigmas => $final_model -> fixed( 'parameter_type' => 'sigma' ),
1542 #                                       degree       => $degree );
1543 #               $fcon -> write( filename => 'FCON' );
1544               }
1545             }
1546           }
1547
1548           # }}}
1549
1550         }
1551
1552         # This "if" block should conceptually be last, since it is
1553         # something we should only do if the model succeeds. In
1554         # practise it could swap places with at least the tweak inits
1555         # block, but for simplicities sake, lets leave it at the
1556         # bottom.
1557
1558         if( not $marked_for_rerun and $tries < $self -> {'min_retries'} ) {
1559
1560           # {{{ Perturb estimates
1561
1562           #Here we force pertubation when the model is successful.
1563
1564           $tries ++;
1565           $marked_for_rerun = 1;
1566           my $degree = 0.1 * $tries;
1567           if( $self -> {'use_implicit_msfo'} and
1568               $self -> reset_msfo( basic_model => $model,
1569                                    final_model => $final_model ) ){
1570
1571             foreach my $prob ( @{$final_model -> problems} ) {
1572               $prob -> set_random_inits ( degree => $degree );
1573             }
1574
1575             $final_model->_write;
1576
1577           } else {
1578             foreach my $prob ( @{$final_model -> problems} ) {
1579               $prob -> set_random_inits ( degree => $degree );
1580             }
1581
1582             $final_model->_write;
1583           }
1584
1585           $modelfile_tainted = 1;
1586
1587           # }}}
1588
1589         }
1590       }
1591
1592       # -------------- Notes about Final model selection -----------------
1593
1594       # Since we have reruns with pertubation and now also forced (or
1595       # automatic) pertubation the final model is not equal to the
1596       # original model. If we consider the models run in the algoritm
1597       # above, there are four implicit subsets. Those that pass the
1598       # picky test, those that don't pass the picky test but have
1599       # minimization successful, those that don't pass the
1600       # minimization step but produce ofv values and, finaly, those
1601       # that doesn't produce an ofv value. The final model will be the
1602       # model that passes the most tests and have the lowest ofv
1603       # value, and if no ofv value is produced, it will be the basic
1604       # model.
1605
1606
1607       # Get all runs that passed the picky test (if picky is used)
1608
1609       # The order of categories is important. Highes priority goes
1610       # first.
1611
1612       my @selection_categories = ('picky','normal');
1613       my ( %selection, $selected );
1614       foreach my $category ( @selection_categories ) {
1615         $selection{$category}{'best_significant_digits'} = 0;
1616       }
1617       for(my $i = 0; $i < scalar @run_results; $i++ ){
1618         foreach my $category ( @selection_categories ) {
1619
1620           if( $category eq 'picky' ) {
1621
1622             if ( $run_results[$i] -> {'pass_picky'} ) {
1623               if( $run_results[$i] -> {'significant_digits'} >
1624                   $selection{$category}{'best_significant_digits'} ){
1625                 $selection{$category}{'selected'} = ($i+1);
1626                 $selection{$category}{'best_significant_digits'} = $run_results[$i] -> {'significant_digits'};
1627               }
1628             }
1629           } else {
1630
1631             if ( $run_results[$i] -> {'minimization_successful'} ) {
1632               if( $run_results[$i] -> {'significant_digits'} >
1633                   $selection{$category}{'best_significant_digits'} ){
1634                 $selection{$category}{'selected'} = ($i+1);
1635                 $selection{$category}{'best_significant_digits'} = $run_results[$i] -> {'significant_digits'};
1636               }
1637             }
1638           }
1639         }
1640       }
1641
1642       foreach my $category ( @selection_categories ) {
1643         $selected = defined $selection{$category}{'selected'} ?
1644             $selection{$category}{'selected'} : $selected;
1645       }
1646       $selected = defined $selected ? $selected : 1;
1647
1648
1649       open( STAT, '>stats-runs.csv' );
1650       print STAT Dumper \@run_results;
1651       close( STAT );
1652
1653       $self -> copy_model_and_output( final_model   => $final_model,
1654                                       model         => $model,
1655                                       use_run       => $selected ? $selected : '',
1656                                       hit_max_tries => ( $retries > 1 and $tries > $retries ) );
1657
1658       if( $self -> {'nonparametric_etas'} and
1659           ( not -e 'np_etas.lst' or $self -> {'rerun'} >=2 )  ) {
1660
1661         # ---------------------  Create nonp eta model  -----------------------------
1662
1663         # {{{ nonp eta model
1664
1665         if( not -e 'np_etas.lst' or $self -> {'rerun'} >= 2 ){
1666
1667           ui -> print( category => 'execute',
1668                        message  => 'Creating NPETA model' );
1669
1670           my $np_eta_mod = $final_model ->
1671               copy( filename    => $self -> {'directory'}.
1672                     'NM_run'.($run_no+1).'/np_etas.mod',
1673                     target      => 'mem',
1674                     copy_data   => 0,
1675                     copy_output => 0);
1676
1677           my ( $msfo_ref, $junk ) = $final_model ->
1678               _get_option_val_pos( name            => 'MSFO',
1679                                    record_name     => 'estimation' );
1680           # We should have an MSFO file here
1681           for( my $i = 0; $i < scalar @{$msfo_ref}; $i++ ) {
1682             my $msfi = $msfo_ref->[$i][0];
1683             $np_eta_mod -> set_records( problem_numbers => [($i+1)],
1684                                         type            =>'msfi',
1685                                         record_strings  => ["$msfi"]);
1686             $np_eta_mod -> set_records( problem_numbers => [($i+1)],
1687                                         type           => 'nonparametric',
1688                                         record_strings => [ 'ETAS UNCONDITIONAL '.
1689                                                             'MSFO=npmsfo'.($i+1) ] );
1690           }
1691           $np_eta_mod -> remove_option( record_name => 'estimation',
1692                                         option_name => 'MSFO' );
1693           $np_eta_mod -> remove_records( type => 'theta' );
1694           $np_eta_mod -> remove_records( type => 'omega' );
1695           $np_eta_mod -> remove_records( type => 'sigma' );
1696           $np_eta_mod -> remove_records( type => 'table' );
1697           my @nomegas = @{$final_model -> nomegas};
1698
1699           my @max_evals;
1700           for( my $i = 0; $i <= $#nomegas; $i++ ) {
1701             my $marg_str = 'ID';
1702             for( my $j = 1; $j <= $nomegas[$i]; $j++ ) {
1703               $marg_str = $marg_str.' ETA'.$j;
1704             }
1705             $marg_str = $marg_str.' FILE='.$model -> filename.'.nonp_etas'.
1706                 ' NOAPPEND ONEHEADER NOPRINT FIRSTONLY';
1707             $np_eta_mod -> add_records( problem_numbers => [($i+1)],
1708                                         type            => 'table',
1709                                         record_strings  => [ $marg_str ] );
1710             $max_evals[$i] = [0];
1711           }
1712
1713 #       my @table_names = @{$np_eta_mod -> table_names};
1714 #       for ( my $i = 0; $i <= $#table_names; $i++ ) {
1715 #         foreach my $table ( @{$table_names[$i]} ) {
1716 #           $table = 'nonp_'.$table;
1717 #         }
1718 #       }
1719 #       $np_eta_mod -> table_names( new_names => \@table_names );
1720
1721           $np_eta_mod -> maxeval( new_values => \@max_evals );
1722           $np_eta_mod -> _write;
1723
1724           # }}} nonp eta model
1725
1726           # ----------------------  run nonp eta model  -------------------------------
1727
1728           # {{{ run eta model
1729
1730           ui -> print( category => 'execute',
1731                        message  => 'Running NPETA model' );
1732
1733           $nonmem = nonmem -> new( adaptive => $self -> {'adaptive'},
1734                                    modelfile => 'np_etas.mod',
1735                                    version => $nm_version,
1736                                    nm_directory => $self -> {'nm_directory'},
1737                                    nice => $self -> {'nice'},
1738                                    show_version => not $self -> {'run_on_lsf'});
1739
1740           $nonmem -> fsubs( $np_eta_mod -> subroutine_files );
1741           unless( $nonmem -> compile() ){
1742               debug -> die( message => "NONMEM compilation failed:\n" .$nonmem -> error_message );
1743           }
1744           if( $nonmem -> nmtran_message =~ /WARNING/s ){
1745               debug -> warn( level => 1,
1746                              message => "NONMEM Warning: " . $nonmem -> nmtran_message );
1747           }
1748           $nonmem -> execute();
1749
1750           foreach my $table_files( @{$np_eta_mod -> table_names} ){
1751             foreach my $table_file( @{$table_files} ){
1752               my $dir = $model -> directory;
1753               cp( $table_file, $dir  );
1754             }
1755           }
1756
1757           unlink 'nonmem', 'nonmem6', 'nonmem5','nonmem.exe', 'nonmem6_adaptive', 'nonmem5_adaptive';
1758         }
1759         # }}} run eta model
1760
1761       }
1762
1763       # {{{ Nonparametric stuff
1764
1765 #      if( $self -> {'nonparametric_marginals'} ) {
1766
1767         # ---------------------  Create nonp marginal model  -----------------------------
1768
1769 #       print "Creating NPMARG model\n";
1770 #       # {{{ nonp marginal model
1771
1772 #       my $np_marg_mod = $final_model ->
1773 #           copy( filename    => $self -> {'directory'}.
1774 #                 'NM_run'.($run_no+1).'/np_marginals.mod',
1775 #                 target      => 'mem',
1776 #                 copy_data   => 0,
1777 #                 copy_output => 0);
1778
1779 #       my @table_names = @{$np_marg_mod -> table_names};
1780 #       if ( $self -> {'nonparametric_etas'} ) {
1781 #         $np_marg_mod -> remove_records( type => 'table' );
1782 #       } else {
1783 #         for ( my $i = 0; $i <= $#table_names; $i++ ) {
1784 #           foreach my $table ( @{$table_names[$i]} ) {
1785 #             $table = 'nonp_'.$table;
1786 #           }
1787 #         }
1788 #         $np_marg_mod -> table_names( new_names => \@table_names );
1789 #       }
1790
1791 #       my ( $msfo_ref, $junk ) = $final_model ->
1792 #           _get_option_val_pos( problem_numbers => [0],
1793 #                                name            => 'MSFO',
1794 #                                record_name     => 'estimation' );
1795
1796 #       my $final_output = $final_model -> outputs -> [0];
1797 #       my $npomegas     = $final_output -> npomegas;
1798
1799         # We should have an MSFO file here
1800 #       for( my $i = 0; $i < scalar @{$msfo_ref}; $i++ ) {
1801 #         my $n_npomegas   = scalar @{$npomegas->[$i][0]};
1802 #         my $msfi = $msfo_ref->[$i][0];
1803 #         $np_marg_mod -> set_records( problem_numbers => [($i+1)],
1804 #                                      type            => 'msfi',
1805 #                                      record_strings  => [$msfi]);
1806 #         $np_marg_mod -> set_records( problem_numbers => [($i+1)],
1807 #                                      type            => 'nonparametric',
1808 #                                      record_strings  => [ 'MARGINALS UNCONDITIONAL '.
1809 #                                                          'MSFO=np_marg_msfo'.($i+1) ] );
1810 #         my $marg_str = 'ID';
1811 #         print "NPOMEGAS: ".$n_npomegas."\n";
1812 #         for( my $j = 0; $j <= $n_npomegas; $j++ ) {
1813 #           $marg_str = $marg_str.' COM('.($j+1).')=MG'.($j+1);
1814 #         }
1815 #         $marg_str = $marg_str.' FILE='.$final_model->filename.'.marginals'.
1816 #             ' NOAPPEND ONHEADER';
1817 #         $np_marg_mod -> add_records( problem_numbers => [($i+1)],
1818 #                                      type            => 'table',
1819 #                                      record_strings  => [ $marg_str ] );
1820 #         $np_marg_mod -> add_records( problem_numbers => [($i+1)],
1821 #                                      type            => 'abbreviated',
1822 #                                      record_strings  => [ "COMRES=".($n_npomegas+1) ] );
1823 #         $np_marg_mod ->
1824 #             add_marginals_code( problem_numbers => [($i+1)],
1825 #                                 nomegas         => [ $n_npomegas ] );
1826 #       }
1827
1828 #       $np_marg_mod -> remove_option( record_name => 'estimation',
1829 #                                      option_name => 'MSFO' );
1830 #       $np_marg_mod -> remove_records( type => 'theta' );
1831 #       $np_marg_mod -> remove_records( type => 'omega' );
1832 #       $np_marg_mod -> remove_records( type => 'sigma' );
1833
1834
1835 #       my @max_evals;
1836 #       @max_evals = @table_names;
1837 #       map( $_ = [0], @max_evals );
1838
1839 #       $np_marg_mod -> maxeval( new_values => \@max_evals );
1840 #       $np_marg_mod -> _write;
1841
1842 #       # }}} nonp marginal model
1843
1844 #       # ----------------------  run nonp marginal model  -------------------------------
1845
1846 #       print "Running NPMARG model\n";
1847 #       # {{{ run marginal model
1848
1849 #       $nonmem = nonmem -> new( modelfile => 'np_marginals.mod',
1850 #                                version => $nm_version,
1851 #                                nm_directory => $self -> {'nm_directory'},
1852 #                                nice => $self -> {'nice'},
1853 #                                show_version => not $self -> {'run_on_lsf'});
1854
1855 #       $nonmem -> fsubs( $np_marg_mod -> subroutine_files );
1856 #       unless( $nonmem -> compile() ){
1857 #         debug -> die( message => $nonmem -> error_message );
1858 #       }
1859 #       if( $nonmem -> nmtran_message =~ /WARNING/s ){
1860 #         debug -> warn( level => 1,
1861 #                        message => "NONMEM Warning: " . $nonmem -> nmtran_message );
1862 #       }
1863 #       $nonmem -> execute();
1864
1865 #       foreach my $table_files( @{$np_marg_mod -> table_names} ){
1866 #         foreach my $table_file( @{$table_files} ){
1867 #           my $dir = $model -> directory;
1868 #           cp( $table_file, $dir  );
1869 #         }
1870 #       }
1871
1872 #       # }}} run marginal model
1873
1874 #      }
1875
1876       # }}}
1877
1878     }
1879 end run_nonmem
1880
1881 # }}}
1882
1883 # {{{ _grid_proxy_init
1884
1885 start _grid_proxy_up
1886 {
1887     if( ($return_value = system( 'grid-proxy-info -exists &>/dev/null' )) ){
1888         if ($? == -1) {
1889             warn "Failed to execute \"grid-proxy-info\": $!\n";
1890         } else {
1891             if( ($return_value = system('grid-proxy-init -valid 168:0')) ){
1892                 if ($? == -1) {
1893                     warn "Failed to execute \"grid-proxy-init\": $!\n";
1894                 }
1895                 warn 'grid-proxy-init failed to setup a proxy\n';
1896             }
1897         }
1898     }
1899 }
1900 end _grid_proxy_up
1901
1902 # }}}
1903
1904 # {{{ ng_submit
1905
1906 start ng_submit
1907     {
1908       # The batch_size is the number of models that go into on
1909       # grid job.
1910
1911 #      $batchSize = $self -> {'grid_batch_size'};
1912
1913 #      my @allmodels = @{ $self -> {'models'} };
1914 #      my @xrsls;
1915 #      my $pm2;
1916       # Four parallel compilation jobs would be ok even on a one-processor machine
1917 #      $pm2 = ext::Parallel::ForkManager -> new (4);
1918 #      for ( my $run = 0; $run <= $#allmodels ; $run += $batchSize) {
1919 #       my $batchEnd = $run + $batchSize - 1 > $#allmodels ? $#allmodels : $run + $batchSize - 1;
1920 #       my @models = @allmodels[$run .. $batchEnd];
1921
1922         # {{{ Filename accumulation And Directory creation.
1923
1924         my (@inputFiles, @executables, @outputFiles );
1925
1926 #       push( @inputFiles, [$INC{'fcon.pm'}, 'fcon.pm'] );
1927
1928 #       push( @xrsls,'-f', $self -> directory . '/psn'. $run .'.xrsl' );
1929
1930 #       my $index = $run;
1931         my $mainExecName = 'batch' . $run . '.pl';
1932
1933 #       my @theta_fixed;
1934 #       my @omega_fixed;
1935 #       my @sigma_fixed;
1936 #       my @retries;
1937 #       my @picky;
1938
1939 #       if( ref( $self -> {'retries'} ) eq 'ARRAY' ){
1940 #         @retries = @{$self -> {'retries'}}[$run .. $batchEnd];
1941 #       }
1942 #       if(ref( $self -> {'picky'} ) eq 'ARRAY' ){
1943 #         @picky = @{$self -> {'picky'}}[$run .. $batchEnd];
1944 #       }
1945
1946 #       foreach my $model ( @models ){
1947 #         my $dirName = 'NM_run' . ($index+1) . '/';
1948
1949 #         push(@theta_fixed, $self -> _twoDimToString( 'array' => $model -> fixed( 'parameter_type' => 'theta' ) ) );
1950 #         push(@omega_fixed, $self -> _twoDimToString( 'array' => $model -> fixed( 'parameter_type' => 'omega' ) ) );
1951 #         push(@sigma_fixed, $self -> _twoDimToString( 'array' => $model -> fixed( 'parameter_type' => 'sigma' ) ) );
1952
1953           push(@inputFiles, ( 'nonmem.sh',
1954                               'nonmem'.$self -> {'nm_version'},
1955                               'FCON',
1956                               'FDATA',
1957                               'FSTREAM' ));
1958
1959           if ( defined $model -> extra_files ) {
1960             foreach my $extra_file( @{$model -> extra_files} ){
1961               if ( $extra_file =~ /^\// ) {
1962                 $extra_file =~ /\/([^\/]*)$/;
1963                 $extra_file = $1;
1964               }
1965               push( @inputFiles, $extra_file );
1966             }
1967           }
1968
1969 #         push(@inFilesToRename, 'FDATA', 'FSTREAM');
1970
1971           for( my $i = 0; $i < $model -> nproblems; $i++ ) {
1972             if( $model -> is_option_set( 'problem_number' => $i,
1973                                          'record' => 'subroutine',
1974                                          'name' => 'LIBRARY' ) ){
1975               push( @inputFiles, 'FLIB');
1976 #             push( @inFilesToRename, 'FLIB');
1977               last;
1978             }
1979           }
1980
1981           push(@executables, 'nonmem'.$self -> {'nm_version'} );
1982
1983           my @table_files = ();
1984           if( $model -> table_names ){
1985             foreach my $table_name_array ( @{$model -> table_names} ) {
1986               @table_files = (@table_files, @{$table_name_array});
1987 #             push(@outFilesToRename, @table_files);
1988             }
1989           }
1990
1991           for( my $i = 0; $i <= $#table_files; $i++ ){
1992             $table_files[$i] = $table_files[$i];
1993           }
1994
1995           push( @outputFiles, ( 'psn.lst',
1996                                 'PRDERR',
1997                                 @table_files ) );
1998
1999 #         unless( ref( $self -> {'retries'} ) eq 'ARRAY' ){
2000 #           push( @retries,$self -> {'retries'});
2001 #         }
2002 #         unless( ref( $self -> {'picky'} ) eq 'ARRAY' ){
2003 #           push( @picky, $self -> {'picky'} );
2004 #         }
2005
2006 #         $index ++;
2007 #       }
2008
2009
2010 #       $self -> _printGridScript( 'batchStart' => $run,
2011 #                                  'batchSize' => $batchSize,
2012 #                                  'theta_fixed' => \@theta_fixed,
2013 #                                  'omega_fixed' => \@omega_fixed,
2014 #                                  'sigma_fixed' => \@sigma_fixed,
2015 #                                  'retries' => \@retries,
2016 #                                  'picky' => \@picky );
2017
2018         # }}}
2019
2020         # {{{ xrsl creation
2021
2022         my $jobname = "nonmem_run_".($run+1);
2023
2024         if( $self -> {'cpu_time'} < 2 ){
2025           # Two minutes is the lowest allowed cpu_time.
2026           $self -> {'cpu_time'} = 2;
2027         }
2028
2029 #       my $total_cpuTime = $self -> {'cpu_time'} * scalar(@models) * 1.2; # We expect 20 % to fail
2030
2031 #       debug -> warn( level   => 2,
2032 #                      message => "grid_time: $total_cpuTime" );
2033
2034         my $xrslfile = grid::nordugrid::xrsl_file -> new( executable => 'nonmem.sh',
2035                                                           inputFiles => \@inputFiles,
2036                                                           executables => \@executables,
2037                                                           outputFiles => \@outputFiles,
2038                                                           ftpThreads => '1',
2039                                                           stdout => 'stdout',
2040                                                           stderr => 'stderr',
2041                                                           #in_base_url => 'gsiftp://joe.asdfbiof.uu.se/ftp/',
2042                                                           #out_base_url => 'gsiftp://joe.biof.uu.se/ftp/',
2043                                                           cpuTime => "\"" . $self -> {'cpu_time'} . " minutes\"",
2044                                                           gmlog => "\"psn.log\"",
2045                                                           jobName => "\"$jobname\"",
2046                                                           filename => 'psn'.($run+1).'.xrsl',
2047                                                           );
2048         $xrslfile -> write;
2049
2050         # }}}
2051
2052         # {{{ nonmem compilation
2053
2054 #       $index = $run;
2055 #       my $wdir = getcwd;
2056
2057 #       $pm2 -> start and next;
2058 #       foreach my $model ( @models ){
2059 #         my $dirName = 'NM_run' . ($index+1);
2060
2061 #         my $tmp_dir = $self -> create_sub_dir( 'subDir' => $dirName . '/' );
2062 #         chdir( $tmp_dir );
2063
2064 #         $self -> copy_model_and_input( 'model' => $model, 'source' => $self -> directory . '/../');
2065
2066 #         my $nm_version = ref( $self -> {'nm_version'} ) eq 'ARRAY' ?
2067 #             $self -> {'nm_version'} -> [$index] : $self -> {'nm_version'};
2068
2069 #         my $nonmem = 'nonmem'->new(modelfile => 'psn.mod',
2070 #                                    version => $nm_version,
2071 #                                    show_version => 0);
2072
2073 #         if( @{$model -> extra_data_files} > 0  ){
2074 #           $nonmem -> fsubs( ['get_sub*.f', 'reader*.f'] );
2075 #         }
2076
2077           my $filehandle;
2078           open( $filehandle, ">", 'nonmem.sh' );
2079           print $filehandle "#!/bin/bash \n";
2080
2081 #         foreach my $file(@inFilesToRename ){
2082 #           print $filehandle "cp $dirName"."_$file $file\n";
2083 #         }
2084 #         my $nmdir = $nonmem -> nm_directory;
2085           print $filehandle "./nonmem".$self -> {'nm_version'}." < FCON > psn.lst \n";
2086           print $filehandle "echo This file was created using the NONMEM version in directory [unknown] >> psn.lst\n";
2087           print $filehandle "echo `date` >> psn.lst\n";
2088           print $filehandle "rm FDATA FSTREAM\n";
2089 #         print $filehandle "mv PRDERR $dirName" . "_PRDERR\n";
2090
2091 #         foreach my $file(@outFilesToRename){
2092 #           print $filehandle "mv $file $dirName"."_$file\n";
2093 #         }
2094
2095         close $filehandle;
2096
2097 #         unless( $nonmem -> compile() ){
2098 #           debug -> die( message => $nonmem -> error_message );
2099 #         }
2100
2101 #         if( $nonmem -> nmtran_message =~ /WARNING/s ){
2102 #           debug -> warn( level => 1,
2103 #                          message => "NONMEM Warning: " . $nonmem -> nmtran_message );
2104 #         }
2105
2106 #         chdir( $wdir );
2107 #         $index ++;
2108 #       }
2109 #       $pm2 -> finish;
2110
2111         # }}}
2112
2113 #      $pm2 -> wait_all_children;
2114
2115       # {{{ grid submit
2116
2117         my ( $iddir, $idfile ) = OSspecific::absolute_path('.', 'job.id' );
2118         $jobidfile = $iddir.$idfile;
2119       unless( $self -> _grid_proxy_up ){
2120         my $clusters;
2121         if ( defined $PsN::config -> {'_'} -> {'ng_clusters'} ) {
2122           $clusters = '-c' . join(' -c', split(/,/, $PsN::config -> {'_'} -> {'ng_clusters'}));
2123         }
2124 #       my ( $out, $err );
2125
2126         unlink( $jobidfile );
2127         # run3 does not work here. We need the command to finish
2128         my $subtries = 0;
2129         my $oldhome = $ENV{'HOME'};
2130         $ENV{'HOME'} = './';
2131         while( $subtries < 6 ) {
2132           my $dbglvl = $subtries <= 1 ? $subtries : 1; # 3 is just too much
2133           my @command = ( 'ngsub', '-d', $dbglvl,'-f',  'psn'.($run+1).'.xrsl',
2134                           '-o', $idfile );
2135           if( $clusters ne '' ){
2136             push( @command, $clusters );
2137           }
2138
2139           open FILE, ">sub_wrapper.sh";
2140           print FILE "#!/bin/bash\n";
2141           print FILE join(" ", @command );
2142           close FILE;
2143           chmod 0755, "sub_wrapper.sh";
2144
2145           my ( $dir, $file ) = OSspecific::absolute_path( '.', "sub_wrapper.sh" );
2146           if( system( $dir.$file ) ){
2147             if ( $! eq 'Bad file descriptor' or
2148                  $! eq 'No such file or directory' ) {
2149               print "ngsub of run ",($run+1)," failed with $!, retrying in ",++$subtries," seconds\n";
2150               sleep($subtries);
2151             } else {
2152               debug -> die( message => "system call @command returned: $!" );
2153             }
2154           } else {
2155             $subtries = 6;
2156           }
2157         }
2158         $ENV{'HOME'} = $oldhome;
2159
2160       } else {
2161         'debug' -> die( message => "Could not submit gridjob, since grid-proxy-init failed" );
2162       }
2163
2164       # }}}
2165 }
2166 end ng_submit
2167
2168 # }}}
2169
2170 # {{{ ng_monitor
2171
2172 start ng_monitor
2173       {
2174         my $opentries = 0;
2175         while( $opentries < 1 ) {
2176           unless( open( JOB, "<" , $jobidfile ) ) {
2177             if ( $opentries == 0 ) {
2178               debug -> die( message => "unable to open jobid file for run ".($run+1) );
2179             } else {
2180               print "unable to open jobid file for run ",($run+1),
2181                 ", retrying in 5 seconds\n";
2182               sleep(5);
2183               $opentries++
2184             }
2185           } else {
2186             $opentries = 1;
2187           }
2188         }
2189
2190         my @jobs = <JOB>;
2191         close( JOB );
2192
2193         debug -> die( message => "More than one job id was found in $jobidfile" )
2194           if ( $#jobs > 0 );
2195
2196         my $job_id = $jobs[0];
2197         chomp( $job_id );
2198         my $status = 'SUBMITTED';
2199         my $job_name;
2200
2201         my $eager = 0;
2202         my $first_poll = 1; # Do some polling immediately
2203         my $time;
2204         while( not $status =~ /FINISHED/ and
2205                $status ne 'FAILED' ){
2206           if( $eager or $first_poll ){
2207             $first_poll = 0;
2208             $time = 60;
2209           } else {
2210             $time = defined $self -> {'grid_poll_interval'} ? $self -> {'grid_poll_interval'}*60 :
2211               int($self -> {'cpu_time'} * 1.2 * 60 * 0.5);
2212           }
2213           unless( $eager ){
2214             print "Polling T - " . $time . "\n"; #" (or press \'f\' to force an update)\n";
2215           }
2216
2217           my $slept = 0;
2218           while( $slept < $time ) {
2219             if( $PsN::config -> {'_'} -> {'use_keyboard'} ) {
2220               my $key = readkey();
2221               if( $key eq 'f' ){
2222                 print "Forced Poll\n";
2223                 last;
2224               }
2225             }
2226             $slept += sleep(10);
2227           }
2228
2229           my ( @info, $err );
2230
2231           open FILE, ">stat_wrapper.sh";
2232           print FILE "#!/bin/bash\n";
2233           print FILE "ngstat -i $jobidfile";
2234           close FILE;
2235           chmod 0755, "stat_wrapper.sh";
2236
2237           run3( "./stat_wrapper.sh", undef, \@info, \$err );
2238
2239           print "info: @info";
2240           print "err: $err\n";
2241           print "Status: $status\n";
2242
2243           my $became_eager = 0;
2244           foreach my $line (@info){
2245             if( $line =~ /Job (gsiftp:\/\/.*)/ ){
2246               next;
2247             } elsif( $line =~ /Jobname: (.*)/ ){
2248               $job_name = $1;
2249             } elsif( $line =~ /Status: (.*)/ ){
2250               my $new_status = $1;
2251               if( $new_status eq 'FINISHING' ){
2252                 print "Job [ $job_id ] is Finishing, polling every minute.\n"; #" (or press \'f\' to force an update)\n";
2253                 $eager = 1;
2254                 $became_eager = 1;
2255               }
2256               if( $status ne $new_status ) {
2257                 print "Job [ $job_id ] changed status from $status to $1\n";
2258               }
2259               $status = $new_status;
2260             } elsif( $line =~ /Error: (.*)/ ) {
2261               my $error = $1;
2262               if( $status =~ /FINISHED/ ){
2263                 debug -> warn( level   => 1,
2264                                message => "Job $job_id Failed: $error" );
2265                 $status = 'FAILED';
2266               }
2267             } else {
2268               # print "Found nothing: $line";
2269             }
2270           }
2271         }
2272
2273         $job_id =~ /\/([0-9]*)$/;
2274         $job_number = $1;
2275         my $dir = $self -> {'directory'}.'NM_run'.($run+1).'/';
2276         push( @job_files, [ $dir.$job_number.'/PRDERR',  $dir.'PRDERR' ]);
2277         push( @job_files, [ $dir.$job_number.'/psn.lst', $dir.'psn.lst' ]);
2278       }
2279 end ng_monitor
2280
2281 # }}}
2282
2283 # {{{ ng_retrieve
2284
2285 start ng_retrieve
2286     {
2287       my $dir = $self -> {'directory'}.'NM_run'.($run+1).'/';
2288       unless( $self -> _grid_proxy_up ){
2289         my $command = "ngget -i $jobidfile -dir $dir";
2290
2291         open FILE, ">get_wrapper.sh";
2292         print FILE "#!/bin/bash\n";
2293         print FILE $command;
2294         close FILE;
2295         chmod 0755, "get_wrapper.sh";
2296
2297         if( system( "./get_wrapper.sh" ) ){
2298           debug -> die( message => "system command $command returned: $!" );
2299         }
2300       } else {
2301         die "Could not get gridjob ",($run+1),", since grid-proxy-init failed\n";
2302       }
2303       my ( $out, $err );
2304       run3( "cp ".$dir.$job_number.'/* '.$dir, undef, \$out, \$err );
2305     }
2306 end ng_retrieve
2307
2308 # }}}
2309
2310 # {{{ print_finish_message
2311
2312 start print_finish_message
2313     {
2314       my $ui_text;
2315       # Log the run
2316       $ui_text .= sprintf("%3s",$run+1) . sprintf("%25s",$self -> {'models'} -> [$run] -> filename);
2317       my $log_text = $run+1 . ',' . $self -> {'models'} -> [$run] -> filename . ',';
2318       if( $self -> {'verbose'} or $self -> {'quick_summarize'} ){
2319         foreach my $param ( 'ofv', 'covariance_step_successful', 'minimization_message' ) {
2320           my $ests = $final_model -> outputs -> [0] -> $param;
2321           if( $param eq 'minimization_message' ){
2322             $ui_text .= "\n    ---------- Minimization Message ----------\n";
2323           }
2324
2325           # Loop the problems
2326           for ( my $j = 0; $j < scalar @{$ests}; $j++ ) {
2327             if ( ref( $ests -> [$j][0] ) ne 'ARRAY' ) {
2328               $ests -> [$j][0] =~ s/^\s*//;
2329               $ests -> [$j][0] =~ s/\s*$//;
2330               $log_text .= $ests -> [$j][0] .',';
2331               #chomp($ests -> [$j][0]);
2332               $ui_text .= sprintf("%10s",$ests -> [$j][0]);
2333             } else {
2334
2335               # Loop the parameter numbers (skip sub problem level)
2336               for ( my $num = 0; $num < scalar @{$ests -> [$j][0]}; $num++ ) {
2337                 #$ests -> [$j][0][$num] =~ s/^\s*//;
2338                 #$ests -> [$j][0][$num] =~ s/\s*$/\n/;
2339                 $log_text .= $ests -> [$j][0][$num] .',';
2340                 #chomp($ests -> [$j][0][$num]);
2341                 if( $param eq 'minimization_message' ){
2342                   $ui_text .= "    ";
2343                 }
2344                 $ui_text .= sprintf("%12s",$ests -> [$j][0][$num]);
2345               }
2346             }
2347           }
2348
2349           if( $param eq 'minimization_message' ){
2350             $ui_text .= "    ------------------------------------------\n\n";
2351           }
2352         }
2353         ui -> print( category => 'all',
2354                      message  => $ui_text,
2355                      wrap     => 0,
2356                      newline => 0);
2357       }
2358
2359       open( LOG, ">>".$self -> {'logfile'} );
2360       print LOG $log_text;
2361       print LOG "\n";
2362       close LOG;
2363
2364     }
2365 end print_finish_message
2366
2367 # }}}
2368
2369 # {{{ run
2370
2371 start run
2372     {
2373       # run:
2374       #
2375       # Arguments: None
2376       #
2377       # Returns: Like all tools modelfit returns a hash with two keys;
2378       # 'own' and 'subtool' but one where the 'subtool' doesn't have any
2379       # value.
2380       #
2381       # Usage:
2382       #
2383       # modelfit -> run implements the core functionality for
2384       # modelfit, it executes NONMEM on each specified
2385       # modelfile. And moves the outputfile to the directory
2386       # where the script was started from.
2387
2388       my $threads;
2389       my @models;
2390       if ( defined $self -> {'models'} ) {
2391         @models = @{ $self -> {'models'} };
2392       } else {
2393         debug -> die( message => "Have no models!" );
2394       }
2395       my $cwd = getcwd();
2396       chdir( $self -> {'directory'} );
2397
2398       $threads = $self -> {'threads'};
2399       $threads = $#models + 1 if ( $threads > $#models + 1);
2400
2401       # Currently parallel execution is not supported on windows
2402       # platforms. Because of regular expression issues with
2403       # Parallel::Forkmanager.
2404       $threads = 1 if( $Config{osname} eq 'MSWin32' );
2405
2406       ui -> print( category => 'all',
2407                    message  => 'Starting ' . scalar(@models) . ' NONMEM executions. '. $threads .' in parallel.' )
2408           unless $self -> {'parent_threads'} > 1;
2409       ui -> print( category => 'all',
2410                    message  => "Run number\tModel name\tOFV\tCovariance step successful." )  if $self -> {'verbose'};
2411
2412       # Print model-NM_run translation file
2413       open( MNT, ">model_NMrun_translation.txt");
2414       for ( my $run = 0; $run <= $#models; $run ++ ) {
2415         print MNT sprintf("%-40s",$models[$run]->filename),"NM_run",($run+1),"\n";
2416       }
2417       close( MNT );
2418
2419       # {{{ Setup of parallel forkmanager
2420
2421       my $moshog_client;
2422
2423       my $moshog_client;
2424       if( $self -> {'adaptive'} and $self -> {'threads'} > 1 ) {
2425
2426         # Make a moshog request. Answer will be read on run_on_finish
2427
2428         $moshog_client = moshog_client -> new(start_threads => $self -> {'threads'});
2429         # print "Making initial request for: ", scalar @models - $self -> {'threads'} , " threads\n";
2430         $moshog_client -> request( request => scalar @models - $self -> {'threads'} );
2431
2432         # print "Initial gran, got new threads: ", $self -> {'threads'} + $moshog_client -> granted(), "\n";
2433
2434         $threads = $self -> {'threads'} + $moshog_client -> granted();
2435
2436       }
2437
2438       my $pm = ext::Parallel::ForkManager -> new ($threads) if ( $threads > 1 );
2439
2440       # These are variables used in forkmanager callbacks.
2441       #
2442       # $run is the loop variable of the nonmem loop. Can be used to
2443       # keep track of progress. (Which we do in the "adaptive" code).
2444
2445       my $run = 0;
2446       my @children;
2447       my $keepforking = 1;
2448       my $aborting = 0;
2449       my $fork_finished = 0;
2450
2451       if ( $threads > 1 ) {
2452
2453         # This following call defines a function that will be
2454         # called when all children have exited. It detects
2455         # errors from children and kills remaining children.
2456         $pm -> run_on_finish(
2457                               sub { my ( $pid, $exit_code, $ident, $exit_signal ) = @_;
2458                                     my $tries = $exit_code;
2459
2460 #                                   if( $exit_signal ){
2461 #                                     debug -> die( message => "A NONMEM run failed" );
2462 #                                   }
2463
2464                                     if( $tries < 0 ) {
2465
2466 # !!!!!!!!!!!!!!!!!!!!!
2467 # tries starts on 0 so this will likely never occur
2468 # Talk to Pontus
2469
2470                                       if( $self -> {'abort_on_fail'} and not $aborting){
2471                                         $aborting = 1;
2472                                         debug -> die( message => "Some nonmem runs failed, aborting" );
2473                                         $keepforking = 0;
2474                                         foreach my $child ( @children ){
2475                                           if( $child ){
2476                                             if( kill( $child ) != 1 ){
2477                                               sleep(3);
2478                                               kill(9, $child);
2479                                             }
2480                                           }
2481                                         }
2482                                       } else {
2483                                         debug -> warn( level   => 0,
2484                                                        message => "One child died" );
2485                                       }
2486                                     } else {
2487                                       push( @{$self -> {'prepared_models'}[$ident]{'own'}},
2488                                             model -> new( filename => "./NM_run" . ($ident+1) ."/psn.mod",
2489                                                           target   => 'disk',
2490                                                           ignore_missing_files => 1 ) );
2491                                       $self -> {'ntries'} -> [$ident] = $tries;
2492                                       # $self -> {'evals'} -> [$ident] = \@evals;
2493
2494                                       my $modulus = (($#models+1) <= 10) ? 1 : (($#models+1) / 10)+1;
2495
2496                                       if ( $ident % $modulus == 0
2497                                            or
2498                                            $ident == 0
2499                                            or
2500                                            $ident == $#models ) {
2501                                         if( $fork_finished ){
2502                                           ui -> print( category => 'all',
2503                                                        message  => 'F:'.($ident+1).' .. ',
2504                                                        wrap => 0,
2505                                                        newline => 0)
2506                                               unless ( $self -> {'parent_threads'} > 1 or $self -> {'verbose'} );
2507                                         }
2508                                       }
2509
2510 #                                     if( $self -> {'adaptive'} ){
2511 #                                       my $request = scalar @models - $run - $self -> {'threads'};
2512 #                                       unless( $request + $self -> {'threads'} eq $threads or $request <= 0 ){
2513 #                                         print "Exit Requesting\n"; sleep( rand() % 3 );
2514 #                                         $moshog_client -> request( request => $request );
2515 #                                         $threads = $self -> {'threads'} + $moshog_client -> granted();
2516 #                                         print "Child finished, got new threads: $threads.";
2517
2518 #                                       }
2519
2520 #                                       if( $request <= 0 ){
2521 #                                         $threads = $self -> {'threads'};
2522 #                                       }
2523
2524 #                                       print  "We had ", $pm -> set_max_procs($threads) , " threads \n";
2525 #                                     }
2526
2527                                     }
2528                                   } );
2529         $pm -> run_on_wait(
2530                             sub {
2531
2532                               if( $self -> {'adaptive'} and $self -> {'threads'} > 1 ){
2533                                 my $request = scalar @models - $run - $self -> {'threads'};
2534                                 unless( $request + $self -> {'threads'} eq $threads or $request <= 0 ){
2535                                   $moshog_client -> request( request => $request );
2536                                   $threads = $self -> {'threads'} + $moshog_client -> granted();
2537                                                                 }
2538
2539                                 if( $request <= 0 ){
2540                                   $threads = $self -> {'threads'};
2541                                 }
2542
2543                                 my $old_threads = $pm -> set_max_procs($threads);
2544
2545                                 if( $self -> {'verbose'} and $old_threads != $threads ){
2546                                   ui -> print( category => 'all',
2547                                                message => "Thread: $threads ($old_threads).",
2548                                                newline => 0 );
2549                                 }
2550                               }
2551                             }
2552                             , 5 );
2553
2554       }
2555
2556       # }}}
2557
2558       # {{{ Local execution
2559
2560       my $actual_runs = -1; # Keep track of how many runs we start (only important for the
2561                             # sleep statement below )
2562
2563       # $run is defined above so that we can access it from
2564       # forkmanager callbacks.
2565       for ( $run = 0; $run <= $#models; $run ++ ) {
2566         $actual_runs++;
2567         if ( $threads > 6 and $actual_runs != 0 and not (($actual_runs)%6) ) {
2568           sleep(1);
2569         }
2570         my $done = 0;
2571         my ( @seed, $tries, $final_model, @evals, $evals_ref );
2572         if ( $keepforking ) {
2573
2574           if ( -e $self -> {'models'} -> [$run] -> outputs -> [0] -> full_name ){
2575             mkdir( "./NM_run" . ($run+1) );
2576             open( DONE, ">./NM_run". ($run+1) ."/done" );
2577             print DONE "This is faked\n 0 attepmts\nevals: 0\nseed: 1 1\n" ;
2578             close( DONE );
2579
2580             # Must copy tablefiles if they exist.
2581
2582             cp( $self -> {'models'} -> [$run] -> outputs -> [0] -> full_name, './NM_run' . ($run+1) .'/psn.lst' );
2583           }
2584
2585           if ( -e "./NM_run" . ($run+1) . "/done" and not $self -> {'rerun'} > 0 ) {
2586
2587             # Should check for tablefiles.
2588
2589             $done = 1;
2590             $actual_runs--;
2591             open( DONE, "./NM_run". ($run+1) ."/done" );
2592             my @rows = <DONE>;
2593             close( DONE );
2594             my $junk;
2595             ( $tries, $junk ) = split(' ',$rows[1],3);
2596             @evals = split(' ',$rows[2]);
2597             @seed = split(' ',$rows[3]);
2598             shift( @evals );    # get rid of 'evals'-word
2599             shift( @seed );     # get rid of 'seed'-word
2600             random_set_seed( @seed );
2601
2602             my $modulus = (($#models+1) <= 10) ? 1 : (($#models+1) / 10)+1;
2603
2604             if ( $run % $modulus == 0
2605                  or
2606                  $run == 0
2607                  or
2608                  $run == $#models ) {
2609
2610               ui -> print( category => 'all',
2611                            message  => 'D:'.( $run + 1 ).' .. ',
2612                            wrap => 0,
2613                            newline => 0)
2614                   unless( $self -> {'parent_threads'} > 1 or $self -> {'verbose'} );
2615             }
2616
2617             $final_model = model -> new( filename             => "./NM_run".($run+1)."/psn.mod",
2618                                          target               => 'disk',
2619                                          ignore_missing_files => 1,
2620                                          quick_reload         => 1 );
2621
2622             $self -> print_finish_message( final_model => $final_model,
2623                                            run => $run );
2624           }
2625           $children[$run] = $pm -> start ( $run ) if ( $threads > 1  and not $done);
2626
2627           # We wish to print progress if $threads are 1(that is, no
2628           # forkmanager object exists) or we have multiple threads and
2629           # are in the parent.
2630
2631           if( $threads <= 1 or $children[$run] ) {
2632
2633             # We don't want to print all starting models if they are
2634             # more than ten. But we allways want to print the first
2635             # and last
2636
2637             my $modulus = (($#models+1) <= 10) ? 1 : (($#models+1) / 10);
2638
2639             if ( $run % $modulus == 0
2640                  or
2641                  $run == 0
2642                  or
2643                  $run == $#models ) {
2644
2645               # The unless checks if tool owning the modelfit is
2646               # running more modelfits, in wich case we should be
2647               # silent to avoid confusion. The $done check is made do
2648               # diffrentiate between allready run processes.
2649
2650               ui -> print( category => 'all',
2651                            message  => 'S:'.($run+1).' .. ',
2652                            wrap => 0,
2653                            newline => 0)
2654                   unless ( $self -> {'parent_threads'} > 1 or $done or $self -> {'verbose'} );
2655             }
2656
2657             # Here we skip to the next model if we are running in
2658             # parallel.
2659
2660             next if $threads > 1;
2661           }
2662         } else {
2663           last;
2664         }
2665
2666         if ( not $done ) {
2667           my $cutoff_thetas;
2668           # Thetas that may be fixed to zero if the go below cutoff
2669           if ((ref( $self -> {'cutoff_thetas'} ) eq 'ARRAY')) {
2670             if ((ref($self -> {'cutoff_thetas'}->[$run]) eq 'ARRAY')) {
2671               $cutoff_thetas = $self->{'cutoff_thetas'}->[$run];
2672             } else {
2673               $cutoff_thetas = $self->{'cutoff_thetas'};
2674             }
2675           } else {
2676             if (defined ($self->{'cutoff_thetas'})) {
2677               die("cutoff_thetas: ". $self->{'cutoff_thetas'}. " is not of the correct type in modelfit->run\n");
2678             }
2679           }
2680
2681           my %options_hash = ( 'handle_rounding_errors' => undef,         # <- Handle rounding errors a bit more intelligently
2682                                'handle_hessian_npd' => undef,             # <- Handle hessian not postiv definite a bit more intelligently
2683                                'handle_maxevals' => undef,                # <- Restart after maxeval
2684                                'tweak_inits' => undef,
2685                                'retries' => undef,
2686                                'picky' => undef,
2687                                'significant_digits_rerun' => undef,
2688                                'nm_version' => undef );
2689
2690           foreach my $option ( keys %options_hash ) {
2691             $options_hash{$option} = (ref( $self -> {$option} ) eq 'ARRAY') ?
2692                 $self -> {$option} -> [$run] : $self -> {$option};
2693           }
2694
2695           # This will stop nasty prints from model, output and data
2696           # which are set to print for the scm.
2697           my $old_category = ui -> category();
2698
2699           ui -> category( 'modelfit' );
2700
2701           my $base_dir = getcwd;
2702           my $work_dir = $self -> create_sub_dir( subDir => '/NM_run'.($run+1) );
2703           chdir( $work_dir );
2704
2705           ( $tries, $final_model, $evals_ref ) =
2706               $self -> run_nonmem ( %options_hash,
2707                                     cutoff_thetas               => $cutoff_thetas,
2708                                     model                       => $models[$run],
2709                                     run_no                      => $run );
2710
2711           chdir( $base_dir );
2712
2713           ui -> category( $old_category );
2714
2715           open( DONE, ">>./NM_run". ($run+1) ."/done" );
2716           @seed = random_get_seed;
2717           @evals = @{$evals_ref} if ( defined $evals_ref );
2718           print DONE "This file indicates that PsN has run nonmem",
2719                      "and tried its best to get it through using\n",
2720                      ($tries+1), " number of attempts\n";
2721           print DONE "evals: @evals\n";
2722           print DONE "seed: @seed\n";
2723           close( DONE );
2724
2725         }
2726
2727         if ( not $done ) {
2728           $self -> print_finish_message( final_model => $final_model,
2729                                          run => $run );
2730         }
2731
2732         if ( $threads > 1 and not $done) {
2733           $pm -> finish ( $tries );
2734         } else {
2735           push( @{$self -> {'prepared_models'}[$run]{'own'}}, $final_model );
2736           push( @{$self -> {'ntries'}}, $tries );
2737           # push( @{$self -> {'evals'}}, \@evals );
2738         }
2739       }
2740
2741       ui -> print( category => 'all',
2742                    message  => " done\n") if( $self -> {'parent_threads'} <= 1 and not $self -> {'verbose'} );
2743       ui -> print( category => 'all',
2744                    message  => "Waiting for all NONMEM runs to finish:\n" ) if( $self -> {'parent_threads'} <= 1 and $threads > 1 and not $self -> {'verbose'} );
2745
2746       $fork_finished = 1;
2747       $pm -> wait_all_children if ( $threads > 1 ) ;
2748
2749       ui -> print( category => 'all',
2750                    message  => " done\n" ) if( $self -> {'parent_threads'} <= 1 and $threads > 1 and not $self -> {'verbose'});
2751
2752       # }}}
2753
2754       # {{{ Database code
2755
2756 #       my $dbh;
2757 #       if ( $PsN::config -> {'_'} -> {'use_database'} ) {
2758 #       $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.";databse=psn",
2759 #                             "psn", "psn_test",
2760 #                             {'RaiseError' => 1});
2761
2762 #       for( my $i = 1; $i <= scalar @{$self -> {'models'}}; $i++ ) {
2763 #         my $sth;
2764 #         foreach my $model ( @{$self -> {'prepared_models'}[$i-1]{'own'}} ) {
2765 #           my $model_id = $model -> model_id;
2766 #           if ( defined $model_id ) {
2767 #             $sth = $dbh -> prepare("INSERT INTO psn.tool_model (tool_id,".
2768 #                                    "model_id,prepared_model) ".
2769 #                                    "VALUES ('".$self -> {'tool_id'}."', ".
2770 #                                    $model_id.", 1 )");
2771 #             $sth -> execute;
2772 #           }
2773 #         }
2774 #         $sth -> finish if ( defined $sth );
2775 #       }
2776 #       $dbh -> disconnect;
2777 #       }
2778
2779       # }}}
2780
2781       $self -> prepare_raw_results();
2782
2783       $self -> print_raw_results();
2784
2785        # {{{ Result section
2786
2787 #       my $models_harvest = $self -> harvest_output( accessors => ['filename'],
2788 #                                                   search_models => 1,
2789 #                                                   search_original_models => 0 );
2790
2791 #       # file name result section
2792 #       my @filenames = @{$models_harvest -> {'filename'}};
2793 #       for (my $i = 0; $i <= $#filenames; $i++ ) {
2794 #       my %return_section;
2795 #       $return_section{'name'} = 'filename';
2796 #       $return_section{'labels'} = [];
2797 #       $return_section{'values'} = [[$filenames[$i]{'own'}]];
2798 #       push( @{$self -> {'results'}[$i]{'own'}},\%return_section );
2799 #       }
2800
2801 #       # ntries result section
2802 #       for (my $i = 0; $i <= $#{$self -> {'ntries'}}; $i++ ) {
2803 #       my %return_section;
2804 #       $return_section{'name'} = 'ntries';
2805 #       $return_section{'values'} = [];
2806 #       $return_section{'labels'} = [];
2807 #       $return_section{'values'} = [[ $self -> {'ntries'} -> [$i] ]];
2808 #       push( @{$self -> {'results'}[$i]{'own'}},\%return_section );
2809 #       }
2810
2811 #       # parameter.estimates and minimization result section
2812
2813 #       my $output_harvest = $self -> harvest_output( accessors => ['minimization_successful','ofv','thetas','omegas','sigmas'],
2814 #                                                   search_output => 1 );
2815
2816 #       foreach my $accessor( 'minimization_successful','ofv' ){
2817 #       my @result = @{$output_harvest -> {$accessor}};
2818 #       for ( my $i = 0; $i <= $#result; $i++ ) {
2819 #         my %return_section;
2820 #         my $accessor_renamed = $accessor;
2821 #         $accessor_renamed =~ s/_/./;
2822 #         $return_section{'name'} = $accessor;
2823 #         $return_section{'values'} = [];
2824 #         $return_section{'labels'} = [];
2825 #         $return_section{'values'} = $result[$i]{'own'}[0];
2826 #         push( @{$self -> {'results'}[$i]{'own'}},\%return_section );
2827 #       }
2828 #       }
2829
2830
2831 #       my @thetas = defined $output_harvest -> {'thetas'} ? @{$output_harvest -> {'thetas'}} : ();
2832 #       my @omegas = defined $output_harvest -> {'omegas'} ? @{$output_harvest -> {'omegas'}} : ();
2833 #       my @sigmas = defined $output_harvest -> {'sigmas'} ? @{$output_harvest -> {'sigmas'}} : ();
2834
2835 #       $models_harvest = $self -> harvest_output( accessors => ['labels'],
2836 #                                                accessor_parameters => {parameter_type => 'theta'},
2837 #                                                search_models => 1 );
2838 #       my @theta_labels = @{$models_harvest -> {'labels'}};
2839 #       $models_harvest = $self -> harvest_output( accessors => ['labels'],
2840 #                                                accessor_parameters => {parameter_type => 'omega'},
2841 #                                                search_models => 1 );
2842 #       my @omega_labels = @{$models_harvest -> {'labels'}};
2843 #       $models_harvest = $self -> harvest_output( accessors => ['labels'],
2844 #                                                accessor_parameters => {parameter_type => 'sigma'},
2845 #                                                search_models => 1 );
2846 #       my @sigma_labels = @{$models_harvest -> {'labels'}};
2847
2848 #       my @ests = ();
2849 #       my @labels = ();
2850 #       # Loop over the models
2851 #       for ( my $i = 0; $i <= $#thetas; $i++ ) {
2852 #       my %return_section = ();
2853 #       $return_section{'name'} = 'parameter.estimates';
2854 #       # Loop over the prepared models (only one here (or the same as the orig))
2855 #       my @prep_ests = ();
2856 #       my @prep_labels = ();
2857 #       for ( my $j = 0; $j < scalar @{$thetas[$i]{'own'}}; $j++ ) {
2858 #         # Loop the problems
2859 #         my @prob_ests = ();
2860 #         my @prob_labels = ();
2861 #         for ( my $k = 0; $k < @{$thetas[$i]{'own'}[$j]}; $k++ ) {
2862 #           # Loop the sub problems
2863 #           my @sub_ests = ();
2864 #           for ( my $l = 0; $l < @{$thetas[$i]{'own'}[$j][$k]}; $l++ ) {
2865 #             # Push the parameter estimates
2866 #             my @array;
2867 #             if( defined $thetas[$i]{'own'}[$j][$k][$l] ){
2868 #               push( @array, @{$thetas[$i]{'own'}[$j][$k][$l]} );
2869 #             }
2870 #             if( defined $omegas[$i]{'own'}[$j][$k][$l] ){
2871 #               push( @array, @{$omegas[$i]{'own'}[$j][$k][$l]} );
2872 #             }
2873 #             if( defined $sigmas[$i]{'own'}[$j][$k][$l] ){
2874 #               push( @array, @{$sigmas[$i]{'own'}[$j][$k][$l]} );
2875 #             }
2876 #             push( @sub_ests, \@array );
2877 #           }
2878 #           push( @prob_ests, \@sub_ests );
2879
2880 #           my @array;
2881 #           if( defined $theta_labels[$i]{'own'}[$j][$k] ){
2882 #             push( @array, @{$theta_labels[$i]{'own'}[$j][$k]} );
2883 #           }
2884 #           if( defined $omega_labels[$i]{'own'}[$j][$k] ){
2885 #             push( @array, @{$omega_labels[$i]{'own'}[$j][$k]} );
2886 #           }
2887 #           if( defined $sigma_labels[$i]{'own'}[$j][$k] ){
2888 #             push( @array, @{$sigma_labels[$i]{'own'}[$j][$k]} );
2889 #           }
2890 #           push( @prob_labels, \@array );
2891
2892 #         }
2893 #         push( @prep_ests, \@prob_ests );
2894 #         push( @prep_labels, \@prob_labels );
2895 #       }
2896 #       $return_section{'values'} = \@prep_ests;
2897 #       $return_section{'labels'} = \@prep_labels;
2898 #       push( @{$self -> {'results'}[$i]{'own'}},\%return_section );
2899 #       }
2900
2901 #       @results     = @{$self -> {'results'}};
2902
2903 #       $self -> print_results;
2904        # }}}
2905
2906       chdir( $cwd );
2907     }
2908 end run
2909
2910 # }}} run
2911
2912 # {{{ summarize
2913
2914 start summarize
2915     {
2916
2917       my %raw_header_map;
2918
2919       for( my $i = 0 ; $i <= scalar @{$self -> {'raw_results_header'}}; $i++ ){
2920         $raw_header_map{$self -> {'raw_results_header'} -> [$i]} = $i;
2921       }
2922
2923       my %model_problem_structure;
2924
2925       foreach my $raw_row ( @{$self -> {'raw_results'}} ){
2926         my $model = $raw_row -> [$raw_header_map{'model'}];
2927         my $problem = $raw_row -> [$raw_header_map{'problem'}];
2928         my $subproblem = $raw_row -> [$raw_header_map{'subproblem'}];
2929
2930         $model_problem_structure{ $model }{ $problem }{ $subproblem } = 1;
2931       }
2932
2933       # ------------------------------------------------------------------#
2934
2935       # Model evaluation. Part one. Version 0.1
2936
2937       # This method checks a nonmem output file for flaws.  The parts
2938       # concerning parts of the covariance step obviously needs the
2939       # $COV # record. In addition, the conditioning number test needs
2940       # the PRINT=E # option of the $COV record. The limits for the
2941       # relative tests may be # changed below:
2942
2943       my $correlation_limit = $self -> {'correlation_limit'};      # All correlations above this number
2944                                                                    # will be listed.
2945       my $condition_number_limit = $self -> {'condition_number_limit'};    # An error will be rised if the
2946                                                                            # condition number is greater
2947                                                                            # than this number.
2948       my $near_bound_sign_digits = $self -> {'near_bound_sign_digits'};   # If a parameter estimate is equal
2949                                                                           # to a bound with this many
2950                                                                           # significant digits, a warning
2951                                                                           # will be printed.
2952       my $near_zero_boundary_limit = $self -> {'near_zero_boundary_limit'};  # When the bound is zero, the
2953                                                                              # check above is not valid. Use
2954                                                                              # this limit instead.
2955       my $sign_digits_off_diagonals =  $self -> {'sign_digits_off_diagonals'};   # The off-diagonal elements are
2956                                                                                  # checked against +-1 with this
2957                                                                                  # many significant digits.
2958       my $large_theta_cv_limit = $self -> {'large_theta_cv_limit'};   # Coefficients of variation larger
2959                                                                       # than these numbers will produce
2960                                                                       # warnings.
2961       my $large_omega_cv_limit = $self -> {'large_omega_cv_limit'};
2962       my $large_sigma_cv_limit = $self -> {'large_omega_cv_limit'};
2963
2964       my $confidence_level = $self -> {'confidence_level'};   # Percentage value;
2965       my $precision = $self -> {'precision'};   # Precision in output.
2966
2967       sub acknowledge {
2968         my $name    = shift;
2969         my $outcome = shift;
2970         my $file    = shift;
2971         my $l = (7 - length( $outcome ))/2;
2972         my $c_num = '00';
2973         $c_num    = '07' if ( $outcome eq 'OK' );
2974         $c_num    = '13' if ( $outcome eq 'WARNING' );
2975         $c_num    = '05' if ( $outcome eq 'ERROR' );
2976         # my $text = sprintf( "%-66s%2s%7s%-5s\n\n", $name, '[ ', $outcome. ' ' x $l, ' ]' );
2977
2978         my $text = sprintf( "%-66s%2s%7s%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
2979
2980         #  cprintf( "%-66s%2s\x03$c_num%7s\x030%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
2981         #  my $text = cprintf( "%-66s%2s\x03$c_num%7s\x030%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
2982         ui -> print( category => 'all', message => $text, wrap => 0 );
2983         print $file $text if defined $file;
2984       }
2985
2986       my $raw_line_structure = ext::Config::Tiny -> read( $self -> {'directory'} . "/raw_file_structure" );
2987       my $raw_res_index = 0;
2988
2989       for( my $model = 1; $model <= scalar keys %model_problem_structure; $model ++ ){
2990
2991 #       my $output = $self -> {'models'} -> [ $model - 1 ] -> outputs -> [0];
2992
2993 #       $output -> _read_problems;
2994
2995         open( my $log, ">test.log" );
2996
2997         my $eigen_run     = 1;
2998         my $corr_matr_run = 1;
2999
3000         my $raw_res = $self -> {'raw_results'};
3001
3002         for ( my $problem = 1; $problem <= scalar keys %{$model_problem_structure{ $model }} ; $problem++ ) {
3003           print "PROBLEM ",$problem,"\n" if scalar keys %{$model_problem_structure{ $model }} > 1;
3004           for ( my $subproblem = 1; $subproblem <= scalar keys %{$model_problem_structure{ $model }{ $problem }}; $subproblem++ ) {
3005             print "SUBPROBLEM ",$subproblem,"\n" if scalar keys %{$model_problem_structure{ $model }{ $problem }} > 1;
3006
3007             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'minimization_successful'}] == 1 ) { # Is in raw
3008               acknowledge( 'Successful minimization', 'OK', $log );
3009             } else {
3010               acknowledge( 'Termination problems', 'ERROR', $log );
3011             }
3012             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'rounding_errors'}] eq '0' ) { # Is in raw
3013               acknowledge( 'No rounding errors', 'OK', $log );
3014             } else {
3015               acknowledge( 'Rounding errors', 'ERROR', $log );
3016             }
3017
3018             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'zero_gradients'}] eq '0' ) { # Is in raw
3019               acknowledge( 'No zero gradients', 'OK', $log );
3020             } else {
3021               acknowledge( 'Zero gradients found '.
3022                            $raw_res -> [ $raw_res_index ][$raw_header_map{'zero_gradients'}].
3023                            ' times', 'ERROR', $log );
3024             }
3025
3026             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'final_zero_gradients'}] eq '0' ) { # Is in raw
3027               acknowledge( 'No final zero gradients', 'OK', $log );
3028             } else {
3029               acknowledge( 'Final zero gradients', 'ERROR', $log );
3030             }
3031
3032             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'hessian_reset'}] eq '0' ) { # Is in raw
3033               acknowledge( 'Hessian not reset', 'OK', $log );
3034             } else {
3035               acknowledge( 'Hessian reset '.
3036                            $raw_res -> [ $raw_res_index ][$raw_header_map{'hessian_reset'}].
3037                            ' times', 'WARNING', $log );
3038             }
3039
3040             if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'estimate_near_boundary'}] eq '0' ){
3041               acknowledge( 'Estimate not near boundary', 'OK', $log );
3042             } else {
3043               acknowledge( 'Estimate near boundary','WARNING',$log);
3044             }
3045
3046             my ( $n_b, $f_b, $f_e );# =
3047 #               $output -> near_bounds( zero_limit         => $near_zero_boundary_limit,
3048 #                                       significant_digits => $near_bound_sign_digits,
3049 #                                       off_diagonal_sign_digits => $sign_digits_off_diagonals ); # Is in raw
3050
3051             if ( defined $n_b -> [$problem] and
3052                  defined $n_b -> [$problem][$subproblem] ) {
3053               my @near_bounds = @{$n_b -> [$problem][$subproblem]};
3054               my @found_bounds = @{$f_b -> [$problem][$subproblem]};
3055               my @found_estimates = @{$f_e -> [$problem][$subproblem]};
3056               if ( $#near_bounds < 0 ) {
3057                 acknowledge( 'No parameter near boundary', 'OK', $log );
3058               } else {
3059                 acknowledge( 'Parameter(s) near boundary', 'WARNING', $log );
3060                 for ( my $problem = 0; $problem <= $#near_bounds; $problem++ ) {
3061                   printf( "\t%-20s%10g   %10g\n", $near_bounds[$problem],
3062                           $found_estimates[$problem], $found_bounds[$problem] );
3063                   print $log sprintf( "\t%-20s%10g   %10g\n", $near_bounds[$problem],
3064                                       $found_estimates[$problem], $found_bounds[$problem] );
3065                   print "\n" if ( $problem == $#near_bounds );
3066                   print $log "\n" if ( $problem == $#near_bounds );
3067                 }
3068               }
3069             }
3070
3071             if ( $raw_res -> [ $raw_res_index ][ $raw_header_map{'covariance_step_run'} ] ) { # Is in raw
3072               if (  $raw_res -> [ $raw_res_index ][$raw_header_map{'covariance_step_successful'}] eq '0' ) {
3073                 acknowledge( 'Covariance step', 'ERROR', $log );
3074               } else {
3075                 if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'covariance_step_successful'}] eq '1' ) {
3076                   acknowledge( 'Covariance step ', 'OK', $log );
3077                 }else {
3078                   acknowledge( 'Covariance step', 'WARNING', $log );
3079                 }
3080
3081                 my ( $l_se, $f_cv );# =
3082 #                   $output -> large_standard_errors( theta_cv_limit => $large_theta_cv_limit,
3083 #                                                     omega_cv_limit => $large_omega_cv_limit,
3084 #                                                     sigma_cv_limit => $large_sigma_cv_limit ); # Is NOT in raw
3085                 if ( defined $l_se -> [$problem] and
3086                      defined $l_se -> [$problem][$subproblem] ) {
3087                   my @large_standard_errors = @{$l_se -> [$problem][$subproblem]};
3088                   my @found_cv = @{$f_cv -> [$problem][$subproblem]};
3089                   if ( $#large_standard_errors < 0 ) {
3090                     acknowledge( 'No large standard errors found', 'OK', $log );
3091                   } else {
3092                     acknowledge( 'Large standard errors found', 'WARNING', $log );
3093                     for ( my $problem = 0; $problem <= $#large_standard_errors; $problem++ ) {
3094                       printf( "\t%-20s%10g\n", $large_standard_errors[$problem], $found_cv[$problem] );
3095                       print $log sprintf( "\t%-20s%10g\n", $large_standard_errors[$problem], $found_cv[$problem] );
3096                       print "\n" if ( $problem == $#large_standard_errors );
3097                       print $log "\n" if ( $problem == $#large_standard_errors );
3098                     }
3099                   }
3100                 }
3101
3102                 my $eigens = $raw_res -> [ $raw_res_index ][$raw_header_map{'eigens'}]; # Is NOT in raw
3103                 if ( $eigens ){# defined $eigens and
3104                     # defined $eigens -> [$problem][$subproblem] and
3105                     # scalar @{$eigens -> [$problem][$subproblem]} > 0 ) {
3106                   if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'condition_number'}] < $condition_number_limit ) {
3107                     acknowledge( 'Condition number ', 'OK', $log );
3108                   } else {
3109                     acknowledge( 'Large condition number ('.
3110                                  sprintf("%6.0f",
3111                                          $raw_res -> [ $raw_res_index ][$raw_header_map{'condition_number'}]).
3112                                  ')', 'WARNING', $log );
3113                   }
3114                 } else {
3115                   $eigen_run = 0;
3116                 }
3117
3118                 if ( $raw_res -> [ $raw_res_index ][$raw_header_map{'s_matrix_singular'}] ne '1' ) { # Is in raw
3119                   acknowledge( 'S-matrix ', 'OK', $log );
3120                 } else {
3121                   acknowledge( 'S-matrix is singular', 'WARNING', $log );
3122                 }
3123
3124                 my ( $h_c, $f_c );# = $output high_correlations( limit => $correlation_limit ); # Is NOT in raw
3125                 if ( defined $h_c -> [$problem] and
3126                      defined $h_c -> [$problem][$subproblem] ) {
3127                   my @high_correlations = @{$h_c -> [$problem][$subproblem]};
3128                   my @found_correlations = @{$f_c -> [$problem][$subproblem]};
3129                   if ( $#high_correlations < 0 ) {
3130                     acknowledge( 'Correlations', 'OK', $log );
3131                   } else {
3132                     acknowledge( 'Large correlations between parameter estimates found', 'WARNING', $log );
3133                     for ( my $problem = 0; $problem <= $#high_correlations; $problem++ ) {
3134                       printf( "\t%-20s%10g\n", $high_correlations[$problem],
3135                               $found_correlations[$problem] );
3136                       print $log sprintf( "\t%-20s%10g\n", $high_correlations[$problem],
3137                                           $found_correlations[$problem] );
3138                       print "\n" if ( $problem == $#high_correlations );
3139                       print $log "\n" if ( $problem == $#high_correlations );
3140                     }
3141                   }
3142                 }
3143               }
3144             } else {
3145               print "No Covariance step run\n";
3146               print $log "No Covariance step run\n";
3147             }
3148
3149             # sumo script
3150
3151             my $confidence_interval = 0;
3152             my $csv = 0;
3153             my $form = '%.' . $precision . 'g';
3154
3155             my @output_matrix;
3156             my @output_matrix_sizes;
3157
3158             my %c_levels = ( '90'   => 1.6449,
3159                              '95'   => 1.96,
3160                              '99'   => 2.5758,
3161                              '99.9' => 3.2905 );
3162
3163             if( $confidence_interval ) {
3164               if( not defined $c_levels{$confidence_level} ) {
3165                 debug -> die( message =>  "Sorry, confidence intervals for level ".$confidence_level.
3166                               " can not be output. Valid levels are: ".join(',', keys %c_levels),
3167                               level => 1 );
3168
3169               }
3170             }
3171
3172             my ( %nam, %est, %cest, %ses );
3173
3174
3175             foreach my $estimate ( 'theta','omega','sigma' ){
3176               my @pos_length = split( /,/ , $raw_line_structure -> {($problem-1)}{$estimate});
3177               my @estimates = @{$raw_res -> [ $raw_res_index ]}[ $pos_length[0] .. $pos_length[0] + $pos_length[1] ];
3178               my @names;
3179               foreach my $num ( 0 .. $pos_length[1] - 1){
3180                 my $template = uc(substr($estimate,0,2));
3181                 push( @names, $template.($num+1));
3182               }
3183
3184               my @pos_length = split( /,/ , $raw_line_structure -> {($problem-1)}{'se'.$estimate});
3185               my @se_estimates = @{$raw_res -> [ $raw_res_index ]}[ $pos_length[0] .. $pos_length[0] + $pos_length[1] ];
3186
3187               my @pos_length = split( /,/ , $raw_line_structure -> {($problem-1)}{'cvse'.$estimate});
3188               my @cvse_estimates = @{$raw_res -> [ $raw_res_index ]}[ $pos_length[0] .. $pos_length[0] + $pos_length[1] ];
3189
3190               $nam{$estimate} = \@names;
3191               $est{$estimate} = \@estimates;
3192               $est{'cvse'.$estimate} = \@cvse_estimates;
3193               $ses{$estimate} = \@se_estimates;
3194             }
3195
3196             my $ofv = $raw_res -> [ $raw_res_index ][ $raw_header_map{'ofv'} ];
3197
3198             if ( defined $ofv ) {
3199               print "Objective function value: ",$ofv,"\n\n";
3200             } else {
3201               print "Objective function value: UNDEFINED\n\n";
3202             }
3203
3204             push( @output_matrix, ["","THETA","","","OMEGA","","","SIGMA", ""] );
3205             for( my $i = 0; $i <= $#{$output_matrix[0]}; $i++ ){
3206               if( $output_matrix_sizes[$i] < length( $output_matrix[0][$i] ) ){
3207                 $output_matrix_sizes[$i] = length( $output_matrix[0][$i] );
3208               }
3209             }
3210
3211             my $max_par = $#{$est{'theta'}};
3212             $max_par = $#{$est{'omega'}} if ( $#{$est{'omega'}} > $max_par );
3213             $max_par = $#{$est{'sigma'}} if ( $#{$est{'sigma'}} > $max_par );
3214
3215             for ( my $max = 0; $max <= $max_par; $max++ ) {
3216               my ( @row, %cis );
3217               if( $confidence_interval ) {
3218                 foreach my $param ( 'theta', 'omega', 'sigma' ) {
3219                   if ( defined $est{$param}[$max] ) {
3220                     my $diff = $c_levels{$confidence_level}*$ses{$param}[$max];
3221                     my ( $lo, $up );
3222                     if( defined $diff ) {
3223                       $lo = $est{$param}[$max]-$diff;
3224                       $up = $est{$param}[$max]+$diff;
3225                     }
3226                     my $cis = sprintf( "($form - $form)", $lo, $up );
3227                     push( @row, $nam{$param}[$max],
3228                           sprintf( $form, $est{$param}[$max] ),
3229                           $cis );
3230                   } else {
3231                     push( @row, '','','' );
3232                   }
3233                 }
3234               } else {
3235                 foreach my $estimate ( 'theta', 'omega','sigma' ){ # Loop over comegas instead
3236                   if ( defined $nam{$estimate}->[$max] ) {
3237                     push( @row, $nam{$estimate}->[$max], defined $est{$estimate}->[$max] ? sprintf( $form, $est{$estimate}->[$max] ) : '........',
3238                           $est{'cvse'.$estimate}->[$max] ? sprintf( "($form)", $est{'cvse'.$estimate}->[$max] ) : '(........)' );
3239                   } else {
3240                     push( @row, '','','' );
3241                   }
3242                 }
3243               }
3244
3245               push(@output_matrix, \@row);
3246               for( my $i = 0; $i <= $#row; $i++ ){
3247                 if( $output_matrix_sizes[$i] < length( $row[$i] ) ){
3248                   $output_matrix_sizes[$i] = length( $row[$i] );
3249                 }
3250               }
3251             }
3252
3253             foreach my $row ( @output_matrix ){
3254               for( my $i = 0; $i <= $#{$row}; $i++ ){
3255                 my $spaces = $output_matrix_sizes[$i] - length($row -> [$i]);
3256                 if( $csv ){
3257                   print $row -> [$i], ",";
3258                 } else {
3259                   print " " x $spaces, $row -> [$i], "  ";
3260                 }
3261               }
3262               print "\n";
3263             }
3264
3265             #$output -> flush;
3266             $raw_res_index ++;
3267           }
3268         }
3269         close( $log );
3270
3271       }
3272     }
3273 end summarize
3274
3275 # }}}