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