?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
removed "remove_temp_files" option and documented "clean" a bit
[PsN.git] / lib / tool / llp_subs.pm
blob81bff08787b5014c8c48e7b8f435d0222e790d50
1 # Det finns ingen mening med att kunna specificera subproblem eftersom vi bara kan fixera
2 # parametervärden på problemnivå.
3 # Det finns dessutom lite värde i att kunna spec. problem eftersom vi än så länge inte kan
4 # köra problem var för sig via PsN och om man har flera problem i en modellfil så skulle
5 # alla köras även om man bara är intresserad av ett. Så man får helt enkelt editera
6 # modellfilen och ta bort onödiga problem.
7
8 # {{{ include
9
10 start include statements
11 use File::Copy 'cp';
12 use ext::Math::SigFigs;
13 use tool::modelfit;
14 use Data::Dumper;
15 end include
16
17 # }}} include statements
18
19 # {{{ new
20
21 start new
22 {
23 my @various_input_formats = ( 'run_thetas', 'run_omegas', 'run_sigmas',
24 'rse_thetas', 'rse_omegas', 'rse_sigmas' );
25 foreach my $var ( @various_input_formats ) {
26 if ( defined $this -> {$var} ) {
27 if ( ref( $this -> {$var} ) eq 'ARRAY' and ref( $this -> {$var} -> [0] ) eq 'ARRAY' ) {
28 if ( scalar @{$this -> {$var}} != scalar @{$this -> {'models'}} ) {
29 debug -> die( message => 'If you define the thetas per model, the first '.
30 'dimension of run_thetas must match the number of models' );
31 }
32 } else {
33 my $variable = $this -> {$var};
34 $this -> {$var} = [];
35 if ( scalar @{$this -> {'models'}} > 1) {
36 for ( my $j = 1; $j <= scalar @{$variable}; $j++ ) {
37 debug -> warn( level => 2,
38 message => "All models: using the same $var $j" );
39 }
40 }
41 foreach my $model ( @{$this -> {'models'}} ) {
42 push ( @{$this -> {$var}}, $variable ) ;
43 }
44 foreach my $param ( 'theta', 'omega', 'sigma' ) {
45 if ( defined $this -> {$var} and defined $this -> {'rse_'.$param}
46 and $var eq 'rse_'.$param.'s' ) {
47 for ( my $i = 0; $i < scalar @{$this -> {$var}}; $i++ ) {
48 debug -> die( message => "The number of run_".$param."s (".
49 scalar @{$this -> {'run_'.$param}[$i]}.
50 ") does not match the number of supplied ".$var.
51 " (".scalar @{$this -> {$var}[$i]}.")" )
52 if ( not scalar @{$this -> {'run_'.$param}[$i]} ==
53 scalar @{$this -> {$var}[$i]} );
54 }
55 }
56 }
57 }
58 }
59 }
60
61 # First check if parameters are of BLOCK type, then we must adjust numbers.
62 # Then check if parameters are of SAME or FIX type, then we will print an ERROR!!
63 if( 0 ) { #-- PONTUS -- This is work in progress. Don't remove please -- PONTUS --
64 my %run;
65 $run{'thetas'} = (ref( $this -> {'run_thetas'} -> [0] ) eq 'ARRAY') ?
66 $this -> {'run_thetas'} -> [0]:$this -> {'run_thetas'};
67 $run{'omegas'} = (ref( $this -> {'run_omegas'} -> [0] ) eq 'ARRAY') ?
68 $this -> {'run_omegas'} -> [0]:$this -> {'run_omegas'};
69 $run{'sigmas'} = (ref( $this -> {'run_sigmas'} -> [0] ) eq 'ARRAY') ?
70 $this -> {'run_sigmas'} -> [0]:$this -> {'run_sigmas'};
71
72 my $model = $this -> {'models'} -> [0];
73
74 # Loop over the parameters
75 foreach my $parameter_type ( keys %run ) {
76 # jump to next parameter if no parameter of this type should be run
77 next unless ( defined $run{$parameter_type} and
78 scalar @{$run{$parameter_type}} > 0 and
79 $run{$parameter_type}->[0] ne '' );
80
81 my @parameter_numbers = @{$run{$parameter_type}};
82
83 foreach my $parameter_number ( @parameter_numbers ){
84
85 if( defined( $model -> problems -> [0] -> $parameter_type )) {
86 my $global_parameter_number = 0;
87 for( my $i = 0; $i < scalar @{$model -> problems -> [0] -> $parameter_type} ; $i++ ){
88 my $init_record = $model -> problems -> [0] -> $parameter_type -> [$i];
89 unless( defined $init_record ){
90 debug -> die( message => "The LLP tool can not be run on $parameter_type number $parameter_number because it does not exist.");
91 }
92 print "$parameter_type $init_record $parameter_number\n";
93 for( my $j = 0 ; $j < scalar @{$init_record -> options}; $j++ ){
94 if( $parameter_number - 1 == $global_parameter_number ){
95 if( $init_record -> same ){
96 debug -> die( message => "The LLP tool can not be run on $parameter_type number $parameter_number because it is \"SAME\"" );
97 } elsif( $init_record -> fix ){
98 debug -> die( message => "The LLP tool can not be run on $parameter_type number $parameter_number because it is \"FIX\"" );
99 }
100 }
101 $global_parameter_number ++;
102 }
103 }
104 }
105 }
106 }
107 }
108
109 foreach my $file ( 'logfile', 'unstacked_logfile', 'raw_results_file' ) {
110 if ( not( ref($this -> {$file}) eq 'ARRAY' or
111 ref($this -> {$file}) eq 'HASH' ) ) {
112 my ($ldir, $filename) = OSspecific::absolute_path( $this -> {'directory'},
113 $this -> {$file} );
114 $this -> {$file} = [];
115 for ( my $i = 1; $i <= scalar @{$this -> {'models'}}; $i++ ) {
116 my $tmp = $filename;
117 $tmp =~ s/\./$i\./;
118 push ( @{$this -> {$file}}, $ldir.$tmp ) ;
119 }
120 }
121 }
122 }
123 end new
124
125 # }}}
126
127
128 # {{{ modelfit_setup
129
130 start modelfit_setup
131 {
132 my $model = $self -> {'models'} -> [$model_number-1];
133
134 my $mfit_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
135 $self -> {'threads'} -> [1]:$self -> {'threads'};
136 my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
137 $self -> {'threads'} -> [0]:$self -> {'threads'};
138
139 # Check which models that hasn't been run and run them
140 # This will be performed each step but will only result in running
141 # models at the first step, if at all.
142
143 # If more than one process is used, there is a VERY high risk of interaction
144 # between the processes when creating directories for model fits. Therefore
145 # the {'directory'} attribute is given explicitly below.
146
147 unless ( $model -> is_run ) {
148
149 # ----------------------- Run original run ------------------------------
150
151 # {{{ orig run
152
153 my %subargs = ();
154 if ( defined $self -> {'subtool_arguments'} ) {
155 %subargs = %{$self -> {'subtool_arguments'}};
156 }
157 my $orig_fit = tool::modelfit ->
158 new( reference_object => $self,
159 models => [$model],
160 threads => $mfit_threads,
161 base_directory => $self -> {'directory'},
162 directory => $self -> {'directory'}.'/orig_modelfit_dir'.$model_number,
163 subtools => [],
164 parent_threads => $own_threads,
165 parent_tool_id => $self -> {'tool_id'},
166 logfile => undef,
167 raw_results => undef,
168 prepared_models => undef,
169 %subargs );
170 # ( models => [$model],
171 # subtools => [],
172
173 # run_on_nordugrid => $self -> {'run_on_nordugrid'},
174 # cpu_time => $self -> {'cpu_time'},
175
176 # run_on_lsf => $self -> {'run_on_lsf'},
177
178 # lsf_queue => $self -> {'lsf_queue'},
179 # lsf_options => $self -> {'lsf_options'},
180 # lsf_job_name => $self -> {'lsf_job_name'},
181 # lsf_project_name => $self -> {'lsf_project_name'},
182
183 # clean => $self -> {'clean'},
184 # nm_version => $self -> {'nm_version'},
185 # picky => $self -> {'picky'},
186 # compress => $self -> {'compress'},
187 # threads => $mfit_threads,
188 # retries => $self -> {'retries'},
189 # base_directory => $self -> {'directory'},
190 # directory => $self -> {'directory'}.'/orig_modelfit_dir'.$model_number.'/',
191 # parent_tool_id => $self -> {'tool_id'},
192 # parent_threads => $own_threads );
193 ui -> print( category => 'llp',
194 message => "Evaluating basic model" ) unless $self -> {'parent_threads'} > 1;
195 $orig_fit -> run;
196
197 # }}} orig run
198
199 }
200
201 my $first = 0;
202 # Prepare the step
203 unless ( defined $self -> {'theta_log'} or
204 defined $self -> {'omega_log'} or
205 defined $self -> {'sigma_log'} ) {
206 # First step, register the estimates from the original runs
207 debug -> warn( level => 2,
208 message => "First llp step, register estimates from original run" );
209 $self -> _register_estimates( first => 1,
210 model_number => $model_number);
211 $first = 1;
212 }
213
214 # If we are resuming, there should exist a 'done' file
215 my $done = 0;
216 if ( -e $self -> {'directory'}.'/m'.$model_number."/done" ) {
217 # Recreate the datasets and models from a checkpoint
218 $done = 1;
219 }
220
221 # Make new guesses for the parameter values
222 $self -> _make_new_guess( first => $first,
223 model_number => $model_number,
224 done => $done );
225
226 # Construct new models for this step from the observed OFV's from the
227 # previous step
228 # They should be placed in m1, m2 and so on for each model.
229 $self -> {'prepared_models'} =
230 $self -> _create_models( model_number => $model_number,
231 done => $done );
232
233 # --------------------- Create the modelfit -------------------------------
234
235 # {{{ modelfit
236
237 if( defined $self -> {'prepared_models'} and scalar @{$self -> {'prepared_models'}} > 0 ) {
238 # Create a modelfit tool for all the models of this step.
239 # This is the last setup part before running the step.
240 my %subargs = ();
241 if ( defined $self -> {'subtool_arguments'} ) {
242 %subargs = %{$self -> {'subtool_arguments'}};
243 }
244 push( @{$self -> {'tools'}},
245 tool::modelfit ->
246 new( reference_object => $self,
247 models => $self -> {'prepared_models'},
248 threads => $mfit_threads,
249 base_directory => $self -> {'directory'},
250 directory => $self -> {'directory'}.'/modelfit_dir'.$model_number,
251 _raw_results_callback => $self -> _modelfit_raw_results_callback,
252 subtools => [],
253 parent_threads => $own_threads,
254 parent_tool_id => $self -> {'tool_id'},
255 raw_results_file => $self -> {'raw_results_file'}[$model_number-1],
256 logfile => undef,
257 raw_results => undef,
258 prepared_models => undef,
259 raw_results_header => undef,
260 tools => undef,
261 %subargs ) );
262 }
263 # $Data::Dumper::Maxdepth = 3;
264 # die Dumper $self -> {'tools'} if not $first;
265
266 # ( raw_results_file => $self -> {'raw_results_file'}[$model_number-1],
267
268 # run_on_lsf => $self -> {'run_on_lsf'},
269
270 # lsf_queue => $self -> {'lsf_queue'},
271 # lsf_options => $self -> {'lsf_options'},
272 # lsf_job_name => $self -> {'lsf_job_name'},
273 # lsf_project_name => $self -> {'lsf_project_name'},
274
275 # run_on_nordugrid => $self -> {'run_on_nordugrid'},
276 # cpu_time => $self -> {'cpu_time'},
277
278 # parent_tool_id => $self -> {'tool_id'},
279
280 # clean => $self -> {'clean'},
281 # models => $self -> {'prepared_models'},
282 # nm_version => $self -> {'nm_version'},
283 # picky => $self -> {'picky'},
284 # compress => $self -> {'compress'},
285 # threads => $mfit_threads,
286 # retries => $self -> {'retries'},
287 # base_directory => $self -> {'directory'},
288 # directory => $self -> {'directory'}.'/modelfit_dir'.$model_number.'/',
289 # subtools => [],
290 # parent_tool_id => $self -> {'tool_id'},
291 # parent_threads => $own_threads,
292 # _raw_results_callback => $self -> _modelfit_raw_results_callback,
293 # %subargs ) );
294
295 # }}} modelfit
296
297 }
298 end modelfit_setup
299
300 # }}}
301
302 # {{{ _create_models
303
304 start _create_models
305 {
306
307 # -------------------- Initiate the run parameters -------------------------
308
309 # {{{ initiate params
310
311 my %run;
312 $run{'thetas'} = (ref( $self -> {'run_thetas'} -> [0] ) eq 'ARRAY') ?
313 $self -> {'run_thetas'} -> [$model_number-1]:$self -> {'run_thetas'};
314 $run{'omegas'} = (ref( $self -> {'run_omegas'} -> [0] ) eq 'ARRAY') ?
315 $self -> {'run_omegas'} -> [$model_number-1]:$self -> {'run_omegas'};
316 $run{'sigmas'} = (ref( $self -> {'run_sigmas'} -> [0] ) eq 'ARRAY') ?
317 $self -> {'run_sigmas'} -> [$model_number-1]:$self -> {'run_sigmas'};
318
319 my $model = $self -> {'models'} -> [$model_number-1];
320
321 # }}} initiate params
322
323 # ------------------ Create the fixed parameter models ---------------------
324
325 # {{{ create models
326
327 # Loop over the parameters
328 foreach my $param ( 'theta', 'omega', 'sigma' ) {
329 # jump to next parameter if no parameter of this type should be run
330 next unless ( defined $run{$param.'s'} and
331 scalar @{$run{$param.'s'}} > 0 and
332 $run{$param.'s'}->[0] ne '' );
333 my @par_nums = @{$run{$param.'s'}};
334 my %par_log = %{$self -> {$param.'_log'}};
335 # Loop over the parameter numbers of interest
336 foreach my $num ( @par_nums ) {
337 my @active = @{$model -> active_problems};
338 my $skip_model = 1;
339 foreach my $val ( @active ) {
340 $skip_model = 0 if ( $val );
341 }
342 foreach my $side ( 'lower', 'upper' ) {
343 # todo: Maybe not necessary to copy data file as well. This is done by
344 # default in model->copy.
345 my $filename = substr($param,0,2).$num.$side.'.mod';
346 my $model_dir = $self -> {'directory'}.'/m'.$model_number.'/';
347 my ($output_dir, $outputfilename) =
348 OSspecific::absolute_path( $self -> {'directory'}.'/m'.$model_number.'/',
349 substr($param,0,2).$num.
350 $side.'.lst');
351 if ( not $done ) {
352
353 # ---------------------- Create model -----------------------------
354
355 # {{{ create
356
357 my $new_mod = $model -> copy( filename => $filename,
358 directory => $model_dir,
359 copy_data => 0,
360 copy_output => 0 );
361 # According to Lars 2005-01-31 removing tables and covariances is a bug.
362 #foreach my $problem ( @{$new_mod -> problems} ) {
363 # $problem -> tables([]);
364 # $problem -> covariances([]);
365 #}
366 $new_mod -> extra_files( $model -> extra_files ),
367 $new_mod -> ignore_missing_files( 1 );
368 $new_mod -> outputfile( $output_dir . $outputfilename );
369 $new_mod -> ignore_missing_files( 0 );
370
371 $new_mod -> _write;
372
373 $new_mod -> update_inits( from_output => $model -> outputs -> [0] );
374 my $active_flag = 0;
375 # Loop over the problems:
376 for ( my $j = 1; $j <= scalar @{$par_log{$num}}; $j++ ) {
377 # Is this side of the problem finished?
378 debug -> warn( level => 2,
379 message => "This side is finished!" )
380 if ( $self->{$param.'_log'}->{$num}->[$j-1]->[2]->{$side} );
381 next if $self->{$param.'_log'}->{$num}->[$j-1]->[2]->{$side};
382 my $sofar = scalar @{$par_log{$num}->[$j-1]->[0]};
383 my $guess;
384 if ( $side eq 'lower' ) {
385 $guess = $par_log{$num}->[$j-1]->[0]->[0];
386 } else {
387 $guess = $par_log{$num}->[$j-1]->[0]->[$sofar-1];
388 }
389 my @diagnostics =
390 @{$new_mod -> initial_values( parameter_type => $param,
391 parameter_numbers => [[$num]],
392 problem_numbers => [$j],
393 new_values => [[$guess]] )};
394 if ( $side eq 'lower' ) {
395 $par_log{$num}->[$j-1]->[0]->[0] = $diagnostics[0][2];
396 } else {
397 $par_log{$num}->[$j-1]->[0]->[$sofar-1] = $diagnostics[0][2];
398 }
399 debug -> warn( level => 2,
400 message => "Registering used value ".
401 $diagnostics[0][2]." for the $side side" );
402 $new_mod -> fixed( parameter_type => $param,
403 parameter_numbers => [[$num]],
404 problem_numbers => [$j],
405 new_values => [[1]] );
406 $active_flag = 1;
407 }
408 if ( $active_flag ) {
409 $new_mod -> _write;
410 push( @new_models, $new_mod );
411 $self->{$param.'_models'}->{$num}->{$side} = $new_mod;
412 }
413
414 # }}} create
415
416 } else {
417
418 # ------------------------- Resume --------------------------------
419
420 # {{{ resume
421
422 my $active_flag = 0;
423 # Loop over the problems:
424 for ( my $j = 1; $j <= scalar @{$par_log{$num}}; $j++ ) {
425 # Is this side of the problem finished?
426 debug -> warn( level => 2,
427 message => "This side is finished!" )
428 if ( $self->{$param.'_log'}->{$num}->[$j-1]->[2]->{$side} );
429 next if $self->{$param.'_log'}->{$num}->[$j-1]->[2]->{$side};
430 $active_flag = 1;
431 }
432 if ( $active_flag ) {
433 my $new_mod = model -> new( directory => $model_dir,
434 filename => $filename,
435 outputfile => $outputfilename,
436 extra_files => $model -> extra_files,
437 target => 'disk',
438 ignore_missing_files => 1 );
439 # Set the correct data file for the object
440 my $moddir = $model -> directory;
441 my @datafiles = @{$model -> datafiles};
442 for( my $df = 0; $df <= $#datafiles; $df++ ) {
443 $datafiles[$df] = $moddir.'/'.$datafiles[$df];
444 }
445 $new_mod -> datafiles( new_names => \@datafiles );
446 push( @new_models, $new_mod );
447 $self->{$param.'_models'}->{$num}->{$side} = $new_mod;
448 }
449
450 # }}} resume
451
452 }
453 }
454 }
455 }
456
457 # }}} create models
458
459 }
460 end _create_models
461
462 # }}} _create_models
463
464 # {{{ _register_estimates
465
466 start _register_estimates
467 {
468
469 # Förenkla loggen: param_log{parameternr}[prob][[estimates...][ofv...]]
470 # Antag att man spec. paramnr som gäller ? alla modeller och problem,
471 # dvs att run_param = [2,4,3] tex.
472 my %run;
473 $run{'thetas'} = (ref( $self -> {'run_thetas'} -> [0] ) eq 'ARRAY') ?
474 $self -> {'run_thetas'} -> [$model_number-1]:$self -> {'run_thetas'};
475 $run{'omegas'} = (ref( $self -> {'run_omegas'} -> [0] ) eq 'ARRAY') ?
476 $self -> {'run_omegas'} -> [$model_number-1]:$self -> {'run_omegas'};
477 $run{'sigmas'} = (ref( $self -> {'run_sigmas'} -> [0] ) eq 'ARRAY') ?
478 $self -> {'run_sigmas'} -> [$model_number-1]:$self -> {'run_sigmas'};
479
480 my %models;
481 my @mod_variants;
482
483 my @prep_models;
484 my ( $orig_output, $orig_ofvs );
485 my $ui_text;
486
487 # Global activity flag for this model:
488 $active = 0;
489
490 my $orig_model = $self -> {'models'} -> [$model_number-1];
491
492 my $mode = $first ? '>' : '>>';
493 open( LOG, $mode.$self -> {'logfile'}[$model_number-1] );
494 open( LOG2, $mode.$self -> {'unstacked_logfile'}[$model_number-1] );
495 print LOG2 sprintf("%10s",'Parameter'),',',sprintf("%10s",'Side'),
496 ',',sprintf("%10s",'Value'),',',sprintf("%10s",'OFV.diff'),"\n" if( $first );
497
498 if ( $first ) {
499 'debug' -> die( message => "No output object defined through model" )
500 unless ( defined $orig_model -> outputs -> [0] );
501 $orig_output = $orig_model -> outputs -> [0];
502 $orig_ofvs = $orig_output -> ofv;
503 # Save the ofvs in {'orig_ofvs'}
504 $self -> {'orig_ofvs'} = $orig_ofvs;
505
506 # Print a log-header
507 foreach my $param ( 'theta', 'omega', 'sigma' ) {
508 next unless ( defined $run{$param.'s'} and
509 scalar @{$run{$param.'s'}} > 0 and
510 $run{$param.'s'}->[0] ne '' );
511 my @par_nums = @{$run{$param.'s'}};
512 my $accessor = $param.'s';
513 my $orig_ests = $orig_output -> $accessor;
514 # Loop the parameter numbers
515 foreach my $num ( @par_nums ) {
516 # Loop the problems
517 for ( my $j = 0; $j < scalar @{$orig_ests}; $j++ ) {
518 my $label = uc($param).$num."_".($j+1);
519 $ui_text = $ui_text.sprintf("%10s",$label).','.sprintf("%10s",'OFV_DIFF').',';
520 print LOG sprintf("%10s",$label),',',sprintf("%10s",'OFV_DIFF'),',';
521 }
522 }
523 }
524 ui -> print( category => 'llp',
525 message => $ui_text,
526 wrap => 0 );
527 print LOG "\n";
528
529 # It is much more time efficient to store all bounds, estimates etc in
530 # an array and access that when doing the registration than calling the
531 # tools within the loop:
532 # Loop over the parameter names
533 $ui_text = '';
534 foreach my $param ( 'theta', 'omega', 'sigma' ) {
535 # jump to next parameter if no parameter of this type should be run
536 next unless ( defined $run{$param.'s'} and
537 scalar @{$run{$param.'s'}} > 0 and
538 $run{$param.'s'}->[0] ne '' );
539 my $accessor = $param.'s';
540 my @par_nums = @{$run{$param.'s'}};
541 my $orig_ests = $orig_output -> $accessor;
542 my $se_accessor = 'se'.$param.'s';
543 my $orig_se_ests = $orig_output -> $se_accessor;
544 # Loop over the parameter numbers of interest
545 my %all_num_est = ();
546 foreach my $num ( @par_nums ) {
547 my @all_prob_est = ();
548 # Loop over the problems:
549 for ( my $j = 0; $j < scalar @{$orig_ests}; $j++ ) {
550 debug -> die( message => "Subproblems are not allowed for the log-likelihood profiling tool" )
551 if ( scalar @{$orig_ests->[$j]} > 1 );
552 # For the original estimates, we need to use the estimates from the original output file.
553 my $orig = $orig_ests->[$j][0][$num-1];
554 # my $orig = $orig_model -> initial_values( parameter_type => 'theta' ) -> [0][$num - 1];
555 # Store original estimates. The [0] is there for print_results.
556 push(@{$self -> {'original_estimates'}[$model_number-1][$j][0]}, $orig);
557 # all_prob_est: [fixed estimates][ofv_diffs][finished_flag(lower,upper)]
558 my %finished = ( 'lower' => 0, 'upper' => 0 );
559 push( @all_prob_est, [[$orig],[0],\%finished]);
560 # If we end up here, the the llp of this model is still active
561 $active = 1;
562 $ui_text = $ui_text.sprintf("%10f",$orig).','.sprintf("%10f",0).',';
563 print LOG sprintf("%10f",$orig),',',sprintf("%10f",0),',';
564 print LOG2 sprintf("%10s",$param.$num),',',sprintf("%10s",'orig'),
565 ',',sprintf("%10f",$orig),',',sprintf("%10f",0),"\n";
566 }
567 $all_num_est{$num} = \@all_prob_est;
568 }
569 $self -> {$param.'_log'} = \%all_num_est;
570 }
571 ui -> print( category => 'llp',
572 message => $ui_text,
573 wrap => 0 );
574 print LOG "\n";
575 } else {
576 @prep_models = @{$self -> {'prepared_models'}};
577
578 # It is much more time efficient to store all bounds, estimates etc in
579 # an array and access that when doing the registration than calling the
580 # tools within the loop:
581 # Loop over the parameter names
582 foreach my $side ( 'lower', 'upper' ) {
583 # reset the user interface text
584 $ui_text = '';
585 foreach my $param ( 'theta', 'omega', 'sigma' ) {
586 # jump to next parameter if no parameter of this type should be run
587 next unless ( defined $run{$param.'s'} and
588 scalar @{$run{$param.'s'}} > 0 and
589 $run{$param.'s'}->[0] ne '' );
590 my %bounds;
591 $bounds{'lower'} =
592 $orig_model -> lower_bounds( parameter_type => $param );
593 $bounds{'upper'} =
594 $orig_model -> upper_bounds( parameter_type => $param );
595
596 my $accessor = $param.'s';
597 my @par_nums = @{$run{$param.'s'}};
598 # get the stored ofvs
599 $orig_ofvs = $self -> {'orig_ofvs'};
600 # Loop over the parameter numbers of interest
601 foreach my $num ( @par_nums ) {
602 my ( $model, $output, $ofvs, $ests );
603 # Check if there exists any model for this side (if not, it is
604 # probably (hopefully) finished
605 if (defined $self->{$param.'_models'} -> {$num} -> {$side}) {
606 $model = $self->{$param.'_models'} -> {$num} -> {$side};
607 $output = $model -> outputs -> [0];
608 $ofvs = $output -> ofv;
609 my $accessor = $param.'s';
610 $ests = $output -> $accessor;
611 }
612 # Loop over the problems:
613 for ( my $j = 0; $j < scalar @{$orig_ofvs}; $j++ ) {
614 # Is this side of the problem finished?
615 if ( $self->{$param.'_log'}->{$num}->[$j]->[2]->{$side} ) {
616 $ui_text = $ui_text.sprintf("%10s",$PsN::out_miss_data).','.sprintf("%10s",$PsN::out_miss_data).',';
617 print LOG sprintf("%10s",$PsN::out_miss_data),',',sprintf("%10s",$PsN::out_miss_data),',';
618 next;
619 }
620 # Collect the model outputs
621 my $diff = $ofvs->[$j][0]-$orig_ofvs->[$j][0];
622 #my $est = $ests -> [$j][0][$num-1];
623 my $est = $model -> initial_values( parameter_type => $param ) -> [0][$num - 1];
624 $ui_text = $ui_text.sprintf("%10f",$est).','.sprintf("%10f",$diff).',';
625 print LOG sprintf("%10f",$est),',',sprintf("%10f",$diff),',';
626 print LOG2 sprintf("%10s",$param.$num),',',sprintf("%10s",$side),',',
627 sprintf("%10f",$est),',',sprintf("%10f",$diff),"\n";
628 my $bound = $bounds{$side}->[$j][$num-1];
629 debug -> warn( level => 2,
630 message => "Registering diff $diff" );
631 if ( $side eq 'lower' ) {
632 unshift( @{$self->{$param.'_log'}->{$num}->[$j]->[1]}, $diff );
633 $est = $self->{$param.'_log'}->{$num}->[$j]->[0]->[0];
634 } else {
635 push( @{$self->{$param.'_log'}->{$num}->[$j]->[1]}, $diff );
636 my $sofar = scalar @{$self->{$param.'_log'}->{$num}->[$j]->[0]};
637 $est = $self->{$param.'_log'}->{$num}->[$j]->[0]->[$sofar-1];
638 }
639 my $finished = $self -> _test_sigdig( number => $diff,
640 goal => $self -> {'ofv_increase'},
641 sigdig => $self -> {'significant_digits'} );
642 my $print_diff = &FormatSigFigs($diff, $self -> {'significant_digits'} );
643 debug -> warn( level => 2,
644 message => "New OFV diff: $print_diff" );
645 debug -> warn( level => 2,
646 message => " equals goal ".$self->{'ofv_increase'}.
647 ". This side is finished" ) if $finished;
648 if ( $finished ) {
649 $self->{$param.'_log'}->{$num}->[$j]->[2]->{$side} = 1;
650 } elsif ( defined $bound and
651 $self -> _test_sigdig( number => $est,
652 goal => $bound,
653 sigdig => 2 ) ) {
654 debug -> warn( level => 1,
655 message => "Estimate $est too close to $side boundary $bound".
656 " terminating search" );
657 $self->{$param.'_log'}->{$num}->[$j]->[2]->{$side} = 2;
658 } else {
659 $active = 1;
660 }
661 $model -> datas -> [0] -> flush;
662 $model -> outputs -> [0] -> flush;
663 }
664 }
665 }
666 # Print to screen if we are in the right context
667 ui -> print( category => 'llp',
668 message => $ui_text,
669 wrap => 0 );
670 print LOG "\n";
671 }
672 }
673 close ( LOG );
674 close ( LOG2 );
675 }
676 end _register_estimates
677
678 # }}} _register_estimates
679
680 # {{{ modelfit_analyze
681
682 start modelfit_analyze
683 {
684 my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
685 $self -> {'threads'} -> [0]:$self -> {'threads'};
686
687 # 1. Register fixed values and ofvs
688 # and check if ofvs are close enough to target
689 my $active = $self -> _register_estimates( model_number => $model_number );
690
691 my ( $returns, $prep_models );
692 # 3. Return logged values-ofv pairs to calling tool
693 if ( $active ) {
694 if ( $self -> {'max_iterations'} > 1 ) {
695 $self -> update_raw_results( model_number => $model_number );
696 # !!! The results_file attribute should not be set when calling
697 # llp recursively. If set, a result file will be written for each
698 # recurrence level which is redundant.
699 my %run;
700 $run{'thetas'} = (ref( $self -> {'run_thetas'} -> [0] ) eq 'ARRAY') ?
701 $self -> {'run_thetas'} -> [$model_number-1]:$self -> {'run_thetas'};
702 $run{'omegas'} = (ref( $self -> {'run_omegas'} -> [0] ) eq 'ARRAY') ?
703 $self -> {'run_omegas'} -> [$model_number-1]:$self -> {'run_omegas'};
704 $run{'sigmas'} = (ref( $self -> {'run_sigmas'} -> [0] ) eq 'ARRAY') ?
705 $self -> {'run_sigmas'} -> [$model_number-1]:$self -> {'run_sigmas'};
706 my %subargs = ();
707 if ( defined $self -> {'subtool_arguments'} ) {
708 %subargs = %{$self -> {'subtool_arguments'}};
709 }
710 my $internal_llp =
711 tool::llp ->
712 new( reference_object => $self,
713 raw_results_file => [$self -> {'raw_results_file'}[$model_number-1]],
714 parent_tool_id => $self -> {'tool_id'},
715 base_directory => $self -> {'directory'},
716 directory => $self -> {'directory'}.'/llp_dir'.$model_number,
717 logfile => [$self -> {'logfile'}[$model_number-1]],
718 unstacked_logfile => [$self -> {'unstacked_logfile'}[$model_number-1]],
719 max_iterations => ($self -> {'max_iterations'} - 1),
720 iteration => ($self -> {'iteration'} + 1),
721 models => [$self -> {'models'} -> [$model_number-1]],
722 run_thetas => [$self -> {'run_thetas'} -> [$model_number-1]],
723 run_omegas => [$self -> {'run_omegas'} -> [$model_number-1]],
724 run_sigmas => [$self -> {'run_sigmas'} -> [$model_number-1]],
725 parent_tool_id => $self -> {'tool_id'},
726 parent_threads => $own_threads,
727 raw_results => undef,
728 prepared_models => undef,
729 rse_thetas => undef,
730 rse_omegas => undef,
731 rse_sigmas => undef,
732 raw_results_header => undef,
733 tools => undef,
734 prepared_models => undef,
735 results => undef,
736 %subargs );
737
738 # $Data::Dumper::Maxdepth = 2;
739 # die Dumper $internal_llp;
740
741 # new ( raw_results_file => [$self -> {'raw_results_file'}[$model_number-1]],
742 # nm_version => $self -> {'nm_version'},
743
744 # run_on_lsf => $self -> {'run_on_lsf'},
745 # lsf_queue => $self -> {'lsf_queue'},
746 # lsf_options => $self -> {'lsf_options'},
747 # lsf_job_name => $self -> {'lsf_job_name'},
748 # lsf_project_name => $self -> {'lsf_project_name'},
749
750 # run_on_nordugrid => $self -> {'run_on_nordugrid'},
751 # cpu_time => $self -> {'cpu_time'},
752
753 # parent_tool_id => $self -> {'tool_id'},
754
755 # picky => $self -> {'picky'},
756 # compress => $self -> {'compress'},
757 # threads => $self -> {'threads'},
758 # retries => $self -> {'retries'},
759 # base_directory => $self -> {'directory'},
760 # directory => $self -> {'directory'}.'/llp_dir'.$model_number.'/',
761 # logfile => [$self -> {'logfile'}[$model_number-1]],
762 # unstacked_logfile => [$self -> {'unstacked_logfile'}[$model_number-1]],
763 # max_iterations => ($self -> {'max_iterations'} - 1),
764 # iteration => ($self -> {'iteration'} + 1),
765 # significant_digits => $self -> {'significant_digits'},
766 # normq => $self -> {'normq'},
767 # ofv_increase => $self -> {'ofv_increase'},
768 # models => [$self -> {'models'} -> [$model_number-1]],
769 # subtools => $self -> {'subtools'},
770 # run_thetas => [$self -> {'run_thetas'} -> [$model_number-1]],
771 # run_omegas => [$self -> {'run_omegas'} -> [$model_number-1]],
772 # run_sigmas => [$self -> {'run_sigmas'} -> [$model_number-1]],
773 # theta_log => $self -> {'theta_log'},
774 # omega_log => $self -> {'omega_log'},
775 # sigma_log => $self -> {'sigma_log'},
776 # orig_ofvs => $self -> {'orig_ofvs'},
777 # parent_tool_id => $self -> {'tool_id'},
778 # parent_threads => $own_threads );
779
780 ( $returns, $prep_models ) = $internal_llp -> run;
781
782 if ( defined $prep_models ) {
783 debug -> warn( level => 2,
784 message => "Inside ".ref($self).
785 " have called internal llp ".
786 scalar @{$prep_models} );
787 push ( @{$self -> {'prepared_models'}[$model_number-1]{'subtools'}},
788 $prep_models );
789 } else {
790 debug -> warn( level => 1,
791 message => "Inside analyze".ref($self).
792 " but no prep_models defined from internal llp" );
793 }
794
795 #push( @{$self -> {'raw_results'}[$model_number-1]},
796 # @{$internal_llp -> raw_results -> [0]} );
797
798 } else {
799 $self -> {'raw_results'}[$model_number-1] =
800 $self -> {'tools'} -> [0] -> raw_results;
801 $self -> update_raw_results( model_number => $model_number );
802 }
803 } else {
804 $self -> {'raw_results'}[$model_number-1] =
805 $self -> {'tools'} -> [0] -> raw_results if( defined $self -> {'tools'} -> [0] );
806 }
807 if( $self -> iteration() < 2 and
808 defined $PsN::config -> {'_'} -> {'R'} and
809 -e $PsN::lib_dir . '/R-scripts/llp.R' ) {
810 # copy the llp R-script
811 cp ( $PsN::lib_dir . '/R-scripts/llp.R', $self -> {'directory'} );
812 # Execute the script
813 system( $PsN::config -> {'_'} -> {'R'}." CMD BATCH llp.R" );
814 }
815 }
816 end modelfit_analyze
817
818 # }}} modelfit_analyze
819
820 # {{{ prepare_results
821
822 start prepare_results
823 {
824 # if ( not defined $self -> {'raw_results'} ) {
825 $self -> read_raw_results();
826 # }
827
828 for ( my $i = 0; $i < scalar @{$self->{'raw_results'}}; $i++ ) { # All models
829 my $orig_mod = $self -> {'models'}[$i];
830 my @params = ( 'theta', 'omega', 'sigma' );
831 my ( %param_nums, %labels, %orig_estimates );
832 foreach my $param ( @params ) {
833 my $modlabels = $orig_mod -> labels( parameter_type => $param );
834 $labels{$param} = $modlabels -> [0]; # Only one problem
835 $param_nums{$param} = scalar @{$modlabels -> [0]} if ( defined $modlabels );
836 my $accessor = $param.'s';
837 $orig_estimates{$param} = $orig_mod -> outputs -> [0] -> $accessor -> [0][0];
838 }
839 my ( %ci, %near_bound, %max_iterations, %interval_ratio,
840 %skewed, %within_interval );
841
842 # The 9 on the row below is offset for iteration,
843 # parameter.type, parameter.number, side, finish.message,
844 # model, problem, subproblem, ofv
845
846 my $cols = scalar @{$self -> {'diagnostic_parameters'}} + 9;
847
848 # my $cols = scalar @{$self -> {'raw_results'}[$i][0]}; # first non-header row
849 # Skip original run:
850 for ( my $j = 1; $j < scalar @{$self -> {'raw_results'}[$i]}; $j++ ) {
851 my $num;
852 my $row = $self -> {'raw_results'}[$i][$j];
853 if ( $row -> [1] eq 'theta' ) {
854 $num = $cols + ($row -> [2] - 1);
855 } elsif ( $row -> [1] eq 'omega' ) {
856 $num = $cols + $param_nums{'theta'} + ($row -> [2] - 1);
857 } else {
858 $num = $cols + $param_nums{'theta'} + $param_nums{'omega'} + ($row -> [2] - 1);
859 }
860 my $row = $self -> {'raw_results'}[$i][$j];
861 $ci{$row -> [1]}{$row -> [2]}{$row -> [3]} = $row -> [$num];
862 if ( $row -> [4] eq 'near.boundary' ) {
863 $near_bound{$row -> [1]}{$row -> [2]}{$row -> [3]} = 1;
864 } elsif ( $row -> [4] eq 'max.iterations' ) {
865 $max_iterations{$row -> [1]}{$row -> [2]}{$row -> [3]} = 1;
866 }
867 }
868 my ( @ci_labels, @ci_values, @li_values );
869 $ci_labels[1] = [ 'lower', 'maximum.likelihood.estimate',
870 'upper', 'interval.ratio', 'near.bound','max.iterations' ];
871 foreach my $param ( @params ) {
872 next if ( not defined $ci{$param} );
873 my @nums = sort { $a <=> $b } keys %{$ci{$param}};
874 foreach my $num ( @nums ) {
875 push( @{$ci_labels[0]}, $labels{$param}[$num-1] );
876 if ( defined $ci{$param}{$num}{'lower'} and
877 defined $ci{$param}{$num}{'upper'} ) {
878 if( abs( $ci{$param}{$num}{'lower'} - $orig_estimates{$param}[$num-1] ) == 0 ){
879 $interval_ratio{$param}{$num} = 'INF';
880 } else {
881 $interval_ratio{$param}{$num} =
882 abs( $ci{$param}{$num}{'upper'} - $orig_estimates{$param}[$num-1] ) /
883 abs( $ci{$param}{$num}{'lower'} - $orig_estimates{$param}[$num-1] );
884 if ( $interval_ratio{$param}{$num} > $self -> {$param.'_interval_ratio_check'} or
885 $interval_ratio{$param}{$num} < 1/$self -> {$param.'_interval_ratio_check'} ) {
886 $skewed{$param}{$num} = 1;
887 } else {
888 $skewed{$param}{$num} = 0;
889 }
890 if ( $self -> {'within_interval_check'} < $ci{$param}{$num}{'upper'} and
891 $self -> {'within_interval_check'} > $ci{$param}{$num}{'lower'} ) {
892 $within_interval{$param}{$num} = 1;
893 } else {
894 $within_interval{$param}{$num} = 0;
895 }
896 }
897 }
898 my @row;
899 push( @row, $ci{$param}{$num}{'lower'} );
900 push( @row, $orig_estimates{$param}[$num-1] );
901 push( @row, $ci{$param}{$num}{'upper'} );
902 push( @row, $interval_ratio{$param}{$num} );
903 push( @row, $near_bound{$param}{$num}{'upper'} ? 1 : $near_bound{$param}{$num}{'lower'} ? 1 : 0 );
904 push( @row, $max_iterations{$param}{$num}{'upper'} ? 1 : $max_iterations{$param}{$num}{'lower'} ? 1 : 0 );
905 push( @ci_values, \@row );
906 }
907 }
908 #print Dumper \%ci;die;
909 $self -> {'confidence_intervals'}[$i] = \%ci;
910 $self -> {'interval_ratio'}[$i] = \%interval_ratio;
911 $self -> {'skewed_intervals'}[$i] = \%skewed;
912 $self -> {'hit_max_iterations'}[$i] = \%max_iterations;
913 $self -> {'near_boundary'}[$i] = \%near_bound;
914 $self -> {'within_interval'}[$i] = \%within_interval;
915 my %return_section;
916 $return_section{'name'} = 'confidence.intervals';
917 $return_section{'labels'} = \@ci_labels;
918 $return_section{'values'} = \@ci_values;
919 unshift( @{$self -> {'results'}[$i]{'own'}},\%return_section );
920 }
921 # # Lasse 2005-04-29: This section may be reduced when the results array turns into a hash
922 # my ( @ofv_log1, @labels0 );
923 # # print Dumper $self -> {'theta_log'};
924 # foreach my $param ( 'theta', 'omega', 'sigma' ) {
925 # my $accessor = $param.'names';
926 # my @param_names = @{$self -> models -> [$model_number -1] -> outputs -> [0] -> $accessor};
927 # my @par_nums = keys %{$self -> {$param.'_log'}};
928 # if ( defined $par_nums[0] ) {
929 # # Loop the problems
930 # for ( my $i = 0; $i < scalar @{$self -> {$param.'_log'}{$par_nums[0]}}; $i++ ) {
931 # foreach my $num ( @par_nums ) {
932 # for ( my $j = 0; $j < scalar @{$self -> {$param.'_log'}{$num}[$i][0]}; $j++ ) {
933 # push( @{$ofv_log1[$i][0]}, [ $self -> {$param.'_log'}{$num}[$i][0][$j],
934 # $self -> {$param.'_log'}{$num}[$i][1][$j] ] );
935 # push( @{$labels0[$i][0][0]}, $param_names[$i][$num-1] );
936 # }
937 # }
938 # $labels0[$i][0][1] = ['fixed.value','OFV'];
939 # }
940 # }
941 # }
942 # # print Dumper \@ofv_log1;
943 # my %return_section;
944 # $return_section{'name'} = 'ofv.trace';
945 # $return_section{'labels'} = \@labels0;
946 # $return_section{'values'} = \@ofv_log1;
947 # push( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
948
949 # my ( @last_values, @labels1, @labels2, @l_dist_ratio, @labels3, @all, @labels4 );
950 # foreach my $param ( 'theta', 'omega', 'sigma' ) {
951 # my $accessor = $param.'names';
952 # my @param_names = @{$self -> models -> [$model_number -1] -> outputs -> [0] -> $accessor};
953 # $accessor = $param.'s';
954 # my @origs = @{$self -> models -> [$model_number -1] -> outputs -> [0] -> $accessor};
955 # my @par_nums = keys %{$self -> {$param.'_log'}};
956 # if ( defined $par_nums[0] ) {
957 # # Loop the problems
958 # for ( my $i = 0; $i < scalar @{$self -> {$param.'_log'}{$par_nums[0]}}; $i++ ) {
959 # foreach my $num ( @par_nums ) {
960 # my $last = scalar @{$self -> {$param.'_log'}{$num}[$i][0]} - 1;
961 # # the first [0] in last_values is there only to add a subproblem level
962 # # for the general print_results method. The first [0] of the param_log
963 # # indicates the log of the fixed parameter values.
964 # # Lower end:
965 # push( @{$last_values[$i][0][0]}, $self -> {$param.'_log'}{$num}[$i][0][0] );
966 # # Upper end:
967 # push( @{$last_values[$i][0][1]}, $self -> {$param.'_log'}{$num}[$i][0][$last] );
968 # my $ratio = (abs( $self -> {$param.'_log'}{$num}[$i][0][$last] -
969 # $origs[$i][0][$num-1] ) / $origs[$i][0][$num-1] ) /
970 # (abs( $self -> {$param.'_log'}{$num}[$i][0][0] -
971 # $origs[$i][0][$num-1] ) / $origs[$i][0][$num-1] );
972 # push( @{$l_dist_ratio[$i][0][0]}, $ratio );
973 # push( @{$all[$i][0]},[ $origs[$i][0][$num-1],
974 # $self -> {$param.'_log'}{$num}[$i][0][0],
975 # $self -> {$param.'_log'}{$num}[$i][0][$last],
976 # $ratio ] );
977 # push( @{$labels1[$i][0][1]}, $param_names[$i][$num-1] );
978 # push( @{$labels2[$i][0][1]}, $param_names[$i][$num-1] );
979 # push( @{$labels3[$i][0][1]}, $param_names[$i][$num-1] );
980 # push( @{$labels4[$i][0][0]}, $param_names[$i][$num-1] );
981 # }
982 # $labels1[$i][0][0] = ['lower', 'upper'];
983 # $labels2[$i][0][0] = [];
984 # $labels3[$i][0][0] = [];
985 # $labels4[$i][0][1] = ['original', 'lower', 'upper', 'ratio'];
986 # }
987 # }
988 # }
989 # my %return_section;
990 # $return_section{'name'} = '';
991 # $return_section{'labels'} = \@labels4;
992 # $return_section{'values'} = \@all;
993 # unshift( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
994 # my %return_section;
995 # $return_section{'name'} = 'limit.distance.ratio';
996 # $return_section{'labels'} = \@labels3;
997 # $return_section{'values'} = \@l_dist_ratio;
998 # unshift( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
999 # my %return_section;
1000 # $return_section{'name'} = 'confidence.intervals';
1001 # $return_section{'labels'} = \@labels1;
1002 # $return_section{'values'} = \@last_values;
1003 # unshift( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
1004 # my %return_section;
1005 # $return_section{'name'} = 'original.estimates';
1006 # $return_section{'labels'} = \@labels2;
1007 # $return_section{'values'} = $self -> {'original_estimates'};
1008 # unshift( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
1009 # # print Dumper $self -> {'results'}[$model_number-1]{'own'};
1010 # }
1011 }
1012 end prepare_results
1013
1014 # }}} prepare_results
1015
1016 # {{{ _modelfit_raw_results_callback
1017
1018 start _modelfit_raw_results_callback
1019 {
1020
1021 # This functions creates a subrouting which will be called by
1022 # the modelfit subtool. The subroutine will add columns to the
1023 # result rows from the modelfit. The result rows will then be
1024 # printed to the "raw_results" file by the modelfit.
1025
1026 # First we create some variables that the subroutine needs. To
1027 # avoid confuision we avoid using the $self variable in the
1028 # subroutine, instead we create these variables which become
1029 # available to the subroutine.
1030
1031 my $iteration = $self -> {'iteration'};
1032 my %run;
1033 my @models = @{$self -> {'models'}};
1034 my $model_number = $self -> {'model_number'};
1035 my ($dir,$file) =
1036 OSspecific::absolute_path( $self -> {'directory'},
1037 $self -> {'raw_results_file'}[$model_number-1] );
1038
1039 $run{'thetas'} = (ref( $self -> {'run_thetas'} -> [0] ) eq 'ARRAY') ?
1040 $self -> {'run_thetas'} -> [$model_number-1]:$self -> {'run_thetas'};
1041 $run{'omegas'} = (ref( $self -> {'run_omegas'} -> [0] ) eq 'ARRAY') ?
1042 $self -> {'run_omegas'} -> [$model_number-1]:$self -> {'run_omegas'};
1043 $run{'sigmas'} = (ref( $self -> {'run_sigmas'} -> [0] ) eq 'ARRAY') ?
1044 $self -> {'run_sigmas'} -> [$model_number-1]:$self -> {'run_sigmas'};
1045 my %log;
1046 foreach my $param ( 'theta', 'omega', 'sigma' ) {
1047 $log{$param} = $self -> {$param.'_log'};
1048 }
1049
1050 # The raw results of the original model should be put as the first row
1051 my $orig_mod = $self -> {'models'}[$model_number-1];
1052
1053 # And here is the subroutine which is given to the modelfit.
1054 my $callback = sub {
1055 # The modelfit object is given to us. It could have been
1056 # retrieved from $self but this was easier.
1057
1058 my $modelfit = shift;
1059 my $mh_ref = shift;
1060 my $raw_results_header = $modelfit -> raw_results_header;
1061 my $raw_results = $modelfit -> raw_results;
1062
1063 $modelfit -> raw_results_append( 1 ) if( $iteration > 1 );
1064
1065 if ( $iteration == 1 ) {
1066
1067 my %dummy;
1068
1069 my ($raw_results_rows, $nonp_rows) = $self -> create_raw_results_rows( max_hash => $mh_ref,
1070 model => $orig_mod,
1071 raw_line_structure => \%dummy );
1072
1073 $orig_mod -> outputs -> [0] -> flush;
1074
1075 unshift( @{$raw_results_rows -> [0]}, ( 0, undef, undef, undef, undef ));
1076
1077 unshift( @{$raw_results}, @{$raw_results_rows} );
1078
1079 # Set the header once.
1080
1081 unshift( @{$raw_results_header} , ('iteration', 'parameter.type',
1082 'parameter.number', 'side', 'finish.message' ) );
1083 }
1084
1085 # {{{ New header
1086
1087 # First prepend llp specific stuff.
1088
1089 # It is implicitly true that the inner loop will execute once
1090 # for each row in the raw_results array. We keep track of its
1091 # the row index with $result_row.
1092
1093 my $result_row = $iteration == 1 ? 1 : 0; # skip the original results row
1094
1095 foreach my $param ( 'theta', 'omega', 'sigma' ) {
1096 foreach my $num ( @{$run{$param.'s'}} ) {
1097 foreach my $side ( 'lower', 'upper' ) {
1098 next unless ( defined $run{$param.'s'} and
1099 scalar @{$run{$param.'s'}} > 0 );
1100 next if $log{$param}{$num}[0][2]{$side};
1101 if ( defined $raw_results -> [$result_row] ) {
1102 unshift( @{$raw_results -> [$result_row]}, ($iteration, $param, $num, $side, undef ) );
1103 }
1104 $result_row++;
1105 }
1106 }
1107 }
1108
1109 # }}} New header
1110
1111 };
1112 return $callback;
1113 }
1114 end _modelfit_raw_results_callback
1115
1116 # }}} _modelfit_raw_results_callback
1117
1118 # {{{ update_raw_results
1119
1120 start update_raw_results
1121 {
1122 my ($dir,$file) = OSspecific::absolute_path( $self -> {'directory'},
1123 $self -> {'raw_results_file'}[$model_number-1] );
1124 open( RRES, $dir.$file );
1125 my @rres = <RRES>;
1126 close( RRES );
1127 open( RRES, '>',$dir.$file );
1128 my @new_rres;
1129 foreach my $row_str ( @rres ) {
1130 chomp( $row_str );
1131 my @row = split( ',', $row_str );
1132 if ( $row[0] eq $self -> {'iteration'} ) {
1133 # The [0] is the problem level, should be removed
1134 if ( $self -> {$row[1].'_log'}{$row[2]}[0][2]{$row[3]} == 1 ) {
1135 $row[4] = 'limit.found';
1136 } elsif ( $self -> {$row[1].'_log'}{$row[2]}[0][2]{$row[3]} == 2 ) {
1137 $row[4] = 'near.boundary';
1138 } elsif ( $self -> {'max_iterations'} <= 1 ) {
1139 $row[4] = 'max.iterations';
1140 }
1141 }
1142 push( @new_rres, \@row );
1143 print RRES join(',',@row ),"\n";
1144 }
1145 close( RRES );
1146 $self -> {'raw_results'}[$model_number-1] = \@new_rres;
1147 }
1148 end update_raw_results
1149
1150 # }}} update_raw_results
1151
1152 # {{{ print_summary
1153
1154 start print_summary
1155 {
1156 sub acknowledge {
1157 my $name = shift;
1158 my $outcome = shift;
1159 my $file = shift;
1160 my $l = (7 - length( $outcome ))/2;
1161 my $c_num = '00';
1162 $c_num = '07' if ( $outcome eq 'OK' );
1163 $c_num = '13' if ( $outcome eq 'WARNING' );
1164 $c_num = '05' if ( $outcome eq 'ERROR' );
1165 # my $text = sprintf( "%-66s%2s%7s%-5s\n\n", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1166 my $text = sprintf( "%-66s%2s%7s%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1167 # cprintf( "%-66s%2s\x03$c_num%7s\x030%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1168 # my $text = cprintf( "%-66s%2s\x03$c_num%7s\x030%-5s", $name, '[ ', $outcome. ' ' x $l, ' ]' );
1169 print $text, "\n\n";
1170 print $file $text if defined $file;
1171 }
1172
1173 my @params = ( 'theta', 'omega', 'sigma' );
1174
1175 sub sum {
1176 my $arr = shift;
1177 my $sum = 0;
1178 foreach my $param ( @params ) {
1179 next if ( not defined $arr -> {$param} );
1180 my @nums = sort { $a <=> $b } keys %{$arr -> {$param}};
1181 foreach my $num ( @nums ) {
1182 if ( ref($arr -> {$param}{$num}) eq 'HASH' ) {
1183 my @sides = sort { $a <=> $b } keys %{$arr -> {$param}{$num}};
1184 foreach my $side ( @sides ) {
1185 $sum += $arr -> {$param}{$num}{$side};
1186 }
1187 } else {
1188 $sum += $arr -> {$param}{$num};
1189 }
1190 }
1191 }
1192 return $sum;
1193 }
1194 my $diag_number = scalar @{$self -> {'diagnostic_parameters'}} - 1;
1195 my %diag_idxs;
1196 for ( my $i = 0; $i <= $diag_number; $i++ ) {
1197 $diag_idxs{$self -> {'diagnostic_parameters'} -> [$i]} = $i;
1198 }
1199
1200 open( my $log, ">test.log" );
1201 for ( my $i = 0; $i < scalar @{$self -> {'raw_results'}} ; $i++ ) { # All models
1202 print "MODEL ",$i+1,"\n" if ( scalar @{$self -> {'raw_results'}} > 1 );
1203 my $sum = sum( $self -> {'within_interval'}[$i] );
1204 if ( not defined $sum or $sum < 1 ) {
1205 acknowledge( 'No confidence intervals include '.
1206 $self -> {'within_interval_check'}, 'OK', $log );
1207 } else {
1208 acknowledge( "$sum confidence intervals include ".
1209 $self -> {'within_interval_check'}, 'WARNING', $log );
1210 foreach my $param ( @params ) {
1211 next if ( not defined $self -> {'within_interval'}[$i]{$param} );
1212 my @nums = sort { $a <=> $b } keys %{$self -> {'within_interval'}[$i]{$param}};
1213 foreach my $num ( @nums ) {
1214 printf( "\t%-20s%3.3f\t%3.3f\n", uc(substr($param,0,2)).$num,
1215 $self -> {'confidence_intervals'}[$i]{$param}{$num}{'lower'},
1216 $self -> {'confidence_intervals'}[$i]{$param}{$num}{'upper'} )
1217 if ( $self -> {'within_interval'}[$i]{$param}{$num} );
1218 print $log sprintf( "\t%-20s\n",uc($param).$num )
1219 if ( $self -> {'within_interval'}[$i]{$param}{$num} );
1220 }
1221 }
1222 }
1223
1224 my $sum = sum( $self -> {'near_boundary'}[$i] );
1225 if ( not defined $sum or $sum < 1 ) {
1226 acknowledge( 'No confidence interval limit search ended near a boundary '.
1227 $self -> {'near_boundary_check'}, 'OK', $log );
1228 } else {
1229 acknowledge( "$sum searches for confidence intervals limits ended near a boundary ",
1230 'WARNING', $log );
1231 foreach my $param ( @params ) {
1232 next if ( not defined $self -> {'near_boundary'}[$i]{$param} );
1233 my @nums = sort { $a <=> $b } keys %{$self -> {'near_boundary'}[$i]{$param}};
1234 foreach my $num ( @nums ) {
1235 next if ( not defined $self -> {'near_boundary'}[$i]{$param} );
1236 foreach my $side ( 'lower', 'upper' ) {
1237 printf( "\t%-20s%-8s%3.3f\n", uc(substr($param,0,2)).$num, $side,
1238 $self -> {'confidence_intervals'}[$i]{$param}{$num}{$side} )
1239 if ( $self -> {'near_boundary'}[$i]{$param}{$num}{$side} );
1240 print $log sprintf( "\t%-20s%-8s%3.3f\n", uc(substr($param,0,2)).$num, $side,
1241 $self -> {'confidence_intervals'}[$i]{$param}{$num}{$side} )
1242 if ( $self -> {'near_boundary'}[$i]{$param}{$num}{$side} );
1243 }
1244 }
1245 }
1246 }
1247
1248 my $sum = sum( $self -> {'hit_max_iterations'}[$i] );
1249 if ( not defined $sum or $sum < 1 ) {
1250 acknowledge( 'No confidence interval limit search reached the maximum number of iterations '.
1251 $self -> {'hit_max_iterations_check'}, 'OK', $log );
1252 } else {
1253 acknowledge( "$sum searches for confidence intervals limits reached the maximum number of iterations ",
1254 'WARNING', $log );
1255 foreach my $param ( @params ) {
1256 next if ( not defined $self -> {'hit_max_iterations'}[$i]{$param} );
1257 my @nums = sort { $a <=> $b } keys %{$self -> {'hit_max_iterations'}[$i]{$param}};
1258 foreach my $num ( @nums ) {
1259 next if ( not defined $self -> {'hit_max_iterations'}[$i]{$param} );
1260 foreach my $side ( 'lower', 'upper' ) {
1261 printf( "\t%-20s%-8s\n", uc(substr($param,0,2)).$num, $side )
1262 if ( $self -> {'hit_max_iterations'}[$i]{$param}{$num}{$side} );
1263 print $log sprintf( "\t%-20s%-8s\n", uc(substr($param,0,2)).$num, $side )
1264 if ( $self -> {'hit_max_iterations'}[$i]{$param}{$num}{$side} );
1265 }
1266 }
1267 }
1268 }
1269
1270 my $sum = sum( $self -> {'skewed_intervals'}[$i] );
1271 if ( not defined $sum or $sum < 1 ) {
1272 acknowledge( 'No CI limit distance ratio exceeds the given thresholds', 'OK', $log );
1273 } else {
1274 acknowledge( "$sum CI limit distance ratios exceed the given thresholds", 'WARNING', $log );
1275 foreach my $param ( @params ) {
1276 next if ( not defined $self -> {'skewed_intervals'}[$i]{$param} );
1277 my @nums = sort { $a <=> $b } keys %{$self -> {'skewed_intervals'}[$i]{$param}};
1278 foreach my $num ( @nums ) {
1279 printf( "\t%-15s%3.3f is less than %3.3f or larger than %3.3f\n", uc(substr($param,0,2)).$num,
1280 $self -> {'interval_ratio'}[$i]{$param}{$num},
1281 (1/$self -> {$param.'_interval_ratio_check'}),
1282 $self -> {$param.'_interval_ratio_check'} )
1283 if ( $self -> {'skewed_intervals'}[$i]{$param}{$num} );
1284 print $log sprintf( "\t%-20s\n",uc($param).$num )
1285 if ( $self -> {'skewed_intervals'}[$i]{$param}{$num} );
1286 }
1287 }
1288 }
1289 }
1290 }
1291 end print_summary
1292
1293 # }}} print_summary
1294
1295 # {{{ _try_bounds
1296
1297 start _try_bounds
1298 {
1299 # Anv䮤 sort f?tt f堦ram l䧳ta och h?a pr?e v?et
1300 # anv䮤 de f?tt fixa estimat utanf?r䮳erna.
1301
1302 if ( ( $side eq 'lower' and $guess < $bound ) or
1303 ( $side eq 'upper' and $guess > $bound ) ) {
1304 my @s_log = sort { $a <=> $b } @param_log;
1305 if ( $side eq 'lower' and $guess < $bound ) {
1306 $guess = ( $bound - $s_log[0])*0.9+$s_log[0];
1307 debug -> warn( level => 1,
1308 message => "Corrected lower guess to $guess" );
1309 }
1310 if ( $side eq 'upper' and $guess > $bound ) {
1311 $guess = ( $bound - $s_log[$#s_log])*0.9+$s_log[$#s_log];
1312 debug -> warn( level => 1,
1313 message => "Corrected upper guess to $guess" );
1314 }
1315 }
1316
1317 $new_value = $guess;
1318 }
1319 end _try_bounds
1320
1321 # }}} _try_bounds
1322
1323 # {{{ _aag
1324
1325 start _aag
1326 {
1327 my $total = scalar(@x);
1328 debug -> die( message => "No data supplied to the polynomial approximation".
1329 " algorithm in the log-likelihood profiling tool" ) if ( $total < 1 );
1330
1331 my $y=0; my $y2=0; my $x1=0; my $x2=0;
1332 my $x3=0; my $x4=0; my $x1y=0; my $x2y=0;
1333
1334 my $count=0;
1335 while ($count<$total){
1336 $y+=$y[$count];
1337 $y2+=($y[$count])**2;
1338 $x1+=$x[$count];
1339 $x2+=($x[$count])**2;
1340 $x3+=($x[$count])**3;
1341 $x4+=($x[$count])**4;
1342 $x1y+=$x[$count]*$y[$count];
1343 $x2y+=(($x[$count])**2)*$y[$count];
1344 ++$count;
1345 }
1346 my $a = $x1y - $x1*$y/$total;
1347 my $b = $x3 - $x1*$x2/$total;
1348 my $c = $x2 - $x1**2/$total;
1349 my $d = $x2y - $x2*$y/$total;
1350 my $e = $x4 - $x2**2/$total;
1351
1352 # Try to avoid division by zero
1353 $c = $c == 0 ? 0.00001 : $c;
1354 my $tmp = ($b**2/$c - $e);
1355 $tmp = $tmp == 0 ? 0.00001 : $tmp;
1356
1357 my $apr1 = ($a*$b/$c - $d) / $tmp;
1358 my $apr2 = $a/$c - $b*$apr1/$c;
1359 my $apr3 = ($y - $apr2*$x1 - $apr1*$x2)/$total;
1360
1361 @polynomial = ($apr1,$apr2,$apr3);
1362 }
1363 end _aag
1364
1365 # }}} _aag
1366
1367 # {{{ _guess
1368
1369 start _guess
1370 {
1371 my @x = @{$param_log[0]};
1372 my @y = @{$param_log[1]};
1373
1374 debug -> die( message => 'The number of logged parameter values ('.
1375 scalar @{$param_log[0]}.
1376 ') does not match the number of logged ofv-diffs ('.
1377 scalar @{$param_log[1]}.')' ) if ( scalar @{$param_log[0]} !=
1378 scalar @{$param_log[1]} );
1379 my ( @x1, @y1 );
1380
1381 my $points = scalar(@x);
1382
1383 my $zero = 0;
1384
1385 while ($y[$zero] > 0){
1386 $zero++;
1387 }
1388
1389 if ( $side eq 'lower' ) {
1390 @x1 = @x[0..2];
1391 @y1 = @y[0..2];
1392
1393 } else {
1394 @x1 = @x[$points-3..$points-1];
1395 @y1 = @y[$points-3..$points-1];
1396 }
1397
1398 my $goal = $y[$zero]+ $self -> {'ofv_increase'};
1399
1400 my @pol = @{$self -> _aag( x => \@x1,
1401 y => \@y1 ) };
1402
1403 # print "X: @x1\n";
1404 # print "Y: @y1\n";
1405 # print "pol: @pol\n";
1406 # print "goal: $goal\n";
1407 if( $pol[0] == 0 ) {
1408 print "The log-likelihood profile could not be approximated by a second order polynomial\n".
1409 "The output may not be correct\n";
1410 $pol[0] = 0.00001;
1411 }
1412
1413 if ( $side eq 'lower' ){
1414 if ($pol[0] > 0){
1415 $guess = -$pol[1]/2/$pol[0] -
1416 (($pol[1]/2/$pol[0])**2 - ($pol[2]-$goal)/$pol[0])**(0.5);
1417 # print "guess: $guess\n";
1418 } else {
1419 $guess = -$pol[1]/2/$pol[0] +
1420 (($pol[1]/2/$pol[0])**2 - ($pol[2]-$goal)/$pol[0])**(0.5);
1421 }
1422
1423 } else {
1424 if ($pol[0] > 0){
1425 $guess = -$pol[1]/2/$pol[0] +
1426 (($pol[1]/2/$pol[0])**2 - ($pol[2]-$goal)/$pol[0])**(0.5);
1427 # print "guess: $guess\n";
1428 } else {
1429 $guess = -$pol[1]/2/$pol[0] -
1430 (($pol[1]/2/$pol[0])**2 - ($pol[2]-$goal)/$pol[0])**(0.5);
1431 }
1432 }
1433
1434 if ($guess eq 'nan' or $guess eq '-1.#IND' ){ #'-1.#IND' - is that a compiler spec. signal?
1435 if ( ($y[0] - $y[1]) == 0 or ($x[0] - $x[1]) == 0 or
1436 ($y[$points-1] - $y[$points-2]) == 0 or ($x[$points-1] - $x[$points-2]) == 0 ) {
1437 $guess = undef;
1438 } else {
1439 if ( $side eq 'lower' ){
1440 $guess = $x[0] - ($goal - $y[0]) / ( ($y[0] - $y[1])/($x[0] - $x[1]));
1441 } else {
1442 $guess = $x[$points-1] + ($goal - $y[$points-1]) / ( ($y[$points-1] - $y[$points-2])/($x[$points-1] - $x[$points-2]));
1443 }
1444 }
1445 }
1446
1447 # print "Goal: $goal\n";
1448 # print LOG "X1: @x1, Y1: @y1\n";
1449 # print "Polynom: @pol\n";
1450 # print "Guess side $side: $guess\n";
1451 }
1452 end _guess
1453
1454 # }}} _guess
1455
1456 # {{{ _test_sigdig
1457
1458 start _test_sigdig
1459 {
1460 $number = &FormatSigFigs($number, $sigdig );
1461 if ( $goal == 0 ) {
1462 $number = sprintf( "%.4f", $number );
1463 $goal = sprintf( "%.4f", $goal );
1464 } else {
1465 $goal = &FormatSigFigs($goal, $sigdig );
1466 }
1467 $test = $number eq $goal ? 1 : 0;
1468 }
1469 end _test_sigdig
1470
1471 # }}} _test_sigdig
1472
1473 # {{{ _make_new_guess
1474
1475 start _make_new_guess
1476 {
1477 # if ( not $done ) {
1478 my %run;
1479 $run{'thetas'} = (ref( $self -> {'run_thetas'} -> [0] ) eq 'ARRAY') ?
1480 $self -> {'run_thetas'} -> [$model_number-1]:$self -> {'run_thetas'};
1481 $run{'omegas'} = (ref( $self -> {'run_omegas'} -> [0] ) eq 'ARRAY') ?
1482 $self -> {'run_omegas'} -> [$model_number-1]:$self -> {'run_omegas'};
1483 $run{'sigmas'} = (ref( $self -> {'run_sigmas'} -> [0] ) eq 'ARRAY') ?
1484 $self -> {'run_sigmas'} -> [$model_number-1]:$self -> {'run_sigmas'};
1485
1486 my $orig_output;
1487 my $orig_model = $self -> {'models'} -> [$model_number-1];
1488 if ( $first ) {
1489 'debug' -> die( message => "No output object defined through model" )
1490 unless ( defined $orig_model -> outputs -> [0] );
1491 $orig_output = $orig_model -> outputs -> [0];
1492 }
1493 # Loop over the parameter names
1494 foreach my $param ( 'theta', 'omega', 'sigma' ) {
1495 # jump to next parameter if no parameter of this type should be run
1496 next unless ( defined $run{$param.'s'} and
1497 scalar @{$run{$param.'s'}} > 0 and
1498 $run{$param.'s'}->[0] ne '' );
1499 my %bounds;
1500 $bounds{'lower'} =
1501 $orig_model -> lower_bounds( parameter_type => $param );
1502 $bounds{'upper'} =
1503 $orig_model -> upper_bounds( parameter_type => $param );
1504 my $accessor = $param.'s';
1505 my @par_nums = @{$run{$param.'s'}};
1506 if ( $first ) {
1507 my $orig_ests = $orig_output -> $accessor;
1508 my $se_accessor = 'se'.$param.'s';
1509 my $orig_se_ests = $orig_output -> $se_accessor;
1510 # Loop over the parameter numbers of interest
1511 foreach my $num ( @par_nums ) {
1512 # Loop over the problems:
1513 for ( my $j = 0; $j < scalar @{$orig_ests}; $j++ ) {
1514 die "Subproblems are not allowed for the log-likelihood profiling tool\n"
1515 if ( scalar @{$orig_ests->[$j]} > 1 );
1516 my $orig = $orig_ests->[$j][0][$num-1];
1517 my $upbnd = $bounds{'upper'}->[$j][$num-1];
1518 my $lobnd = $bounds{'lower'}->[$j][$num-1];
1519 my $width;
1520 if ( defined $orig_se_ests->[$j][0][$num-1] ) {
1521 $width = abs( $orig_se_ests->[$j][0][$num-1] *
1522 $self -> {'normq'} );
1523 } elsif ( defined $self -> {'rse_'.$param.'s'}->[$model_number-1]{$num} ) {
1524 $width = abs( $self -> {'rse_'.$param.'s'}->[$model_number-1]{$num}/100*abs($orig) *
1525 $self -> {'normq'} );
1526 } else {
1527 die "No estimate of the standard error of $param $num is available from the output file nor from the command line\n";
1528 }
1529 my $upper = $orig + $width;
1530 my $lower = $orig - $width;
1531 # print "Upper: $upper, lower: $lower\n";
1532 $lower = ( defined $lobnd and $lower < $lobnd ) ?
1533 ($lobnd-$orig)*0.9+$orig : $lower;
1534 $upper = ( defined $upbnd and $upper > $upbnd ) ?
1535 ($upbnd-$orig)*0.9+$orig : $upper;
1536 # print "Upbnd: $upbnd, lobnd: $lobnd\n";
1537 # print "Upper: $upper, lower: $lower\n";
1538 unshift( @{$self->{$param.'_log'}->{$num}->[$j]->[0]}, $lower );
1539 push( @{$self->{$param.'_log'}->{$num}->[$j]->[0]}, $upper );
1540 }
1541 }
1542 } else {
1543 # Loop over the parameter numbers of interest
1544 foreach my $num ( @par_nums ) {
1545 # Loop over the problems:
1546 for ( my $j = 0; $j < scalar @{$bounds{'lower'}}; $j++ ) {
1547 my %guesses;
1548 foreach my $side ( 'lower', 'upper' ) {
1549 # Is this side of the problem finished?
1550 next if $self->{$param.'_log'}->{$num}->[$j]->[2]->{$side};
1551 # Collect the model outputs
1552 debug -> warn( level => 2,
1553 message => "Making new guess for $param number $num on the $side side" );
1554 my $bound = $bounds{$side}->[$j][$num-1];
1555 my $guess =
1556 $self -> _guess( param_log => $self->{$param.'_log'}->{$num}->[$j],
1557 side => $side );
1558 # print "G1: $num : $guess\n";
1559 if ( defined $bounds{$side}->[$j][$num-1] ) {
1560 $guess =
1561 $self -> _try_bounds( guess => $guess,
1562 side => $side,
1563 bound => $bounds{$side}->[$j][$num-1],
1564 param_log => $self->{$param.'_log'}->
1565 {
1566 $num}->[$j]->[0]);
1567 # print "G2: $guess\n" if $num == 9;
1568 }
1569 $guesses{$side} = $guess;
1570 if ( not defined $guess ) {
1571 print "Warning: The search for the $side CI-limit for $param $num ".
1572 "could not continue due to numerical difficulties\n";
1573 $self->{$param.'_log'}->{$num}->[$j]->[2]->{$side} = 1;
1574 }
1575 }
1576 unshift( @{$self->{$param.'_log'}->{$num}->[$j]->[0]}, $guesses{'lower'} )
1577 if ( defined $guesses{'lower'} );
1578 push( @{$self->{$param.'_log'}->{$num}->[$j]->[0]}, $guesses{'upper'} )
1579 if ( defined $guesses{'upper'} );
1580 }
1581 }
1582 }
1583 }
1584 # Logging must be done fairly quick, therefore this loop by itself
1585 open( DONE, '>'.$self -> {'directory'}."/m$model_number/done" );
1586 foreach my $param ( 'theta', 'omega', 'sigma' ) {
1587 next unless ( defined $self->{$param.'_log'} );
1588 while ( my ( $num, $probs ) = each %{$self->{$param.'_log'}} ) {
1589 # Loop over the problems:
1590 for ( my $prob = 0; $prob < scalar @{$probs}; $prob++ ) {
1591 foreach my $side ( 'lower', 'upper' ) {
1592 next if $self->{$param.'_log'}->{$num}->[$prob]->[2]->{$side};
1593 my $log_size = scalar @{$probs -> [$prob] -> [0]};
1594 if ( $side eq 'lower' ) {
1595 print DONE "$param $num $prob $side ",
1596 $probs -> [$prob] -> [0] -> [0],"\n";
1597 } elsif ( $side eq 'upper' ) {
1598 print DONE "$param $num $prob $side ",
1599 $probs -> [$prob] -> [0] -> [$log_size-1],"\n";
1600 }
1601 }
1602 }
1603 }
1604 }
1605 close( DONE );
1606 # } else {
1607 # open( DONE, $self -> {'directory'}."/m$model_number/done" );
1608 # my @rows = <DONE>;
1609 # close( DONE );
1610 # for( my $i = 0; $i <= $#rows; $i++ ) { # skip first row
1611 # my ( $param, $num, $prob, $side, $guess ) = split(' ',$rows[$i],5);
1612 # next if $self->{$param.'_log'}->{$num}->[$prob]->[2]->{$side};
1613 # unshift( @{$self->{$param.'_log'}->{$num}->[$prob]->[0]}, $guess )
1614 # if ( $side eq 'lower' );
1615 # push( @{$self->{$param.'_log'}->{$num}->[$prob]->[0]}, $guess )
1616 # if ( $side eq 'upper' );
1617 # }
1618 # }
1619 }
1620 end _make_new_guess
1621
1622 # }}} _make_new_guess
1623
1624 # {{{ confidence_limits
1625
1626 start confidence_limits
1627 {
1628 if ( defined $self -> {'results'} ){
1629 for ( my $i = 0; $i < scalar @{$self -> {'results'} -> {'own'}}; $i++ ) {
1630 my %num_lim;
1631 if ( defined $self->{'results'} -> {'own'}->[$i]->{$parameter.'_log'} ) {
1632 # print "$parameter defined\n";
1633 # print "REF: ",ref($self->{'results'}->
1634 # [$i]->{$parameter.'_log'}),"\n";
1635 # print "KEYS: ",keys %{$self->{'results'}->
1636 # [$i]->{$parameter.'_log'}},"\n";
1637 my @nums = sort {$a <=> $b} keys %{$self->{'results'} -> {'own'} ->
1638 [$i]->{$parameter.'_log'}};
1639 foreach my $num ( @nums ) {
1640 my @prob_lim = ();
1641 # print "REF: ",ref($self->{'results'}->
1642 # [$i]->{$parameter.'_log'}->{$num}),"\n";
1643 for ( my $j = 0; $j < scalar @{$self->{'results'} -> {'own'}->
1644 [$i]->{$parameter.'_log'}->{$num}}; $j++ ) {
1645 my @last_lvl = @{$self->{'results'} -> {'own'}->
1646 [$i]->{$parameter.'_log'}->{$num}->[$j]};
1647 push( @prob_lim, [$last_lvl[0][0],$last_lvl[0][scalar @{$last_lvl[0]}-1]] );
1648 # print "REF: ",ref($self->{'results'}->
1649 # [$i]->{$parameter.'_log'}->{$num}->[$j]),"\n";
1650 # print "UP fini: ",$last_lvl[2]->{'upper'},"\n";
1651 # print "LO limit: ",$last_lvl[0][0],"\n";
1652 # print "UP limit: ",$last_lvl[0][scalar @{$last_lvl[0]}-1],"\n";
1653 }
1654 $num_lim{$num} = \@prob_lim;
1655 }
1656 }
1657 push( @limits, \%num_lim );
1658 }
1659 }
1660 }
1661 end confidence_limits
1662
1663 # }}} confidence_limits
1664
1665 # {{{ print_results
1666
1667 start print_results
1668 {
1669 # Run the print_results specific for the subtool
1670 my $sub_print_results = $self -> {'subtools'} -> [0];
1671 if ( defined $sub_print_results ) {
1672 sub get_dim {
1673 my $arr = shift;
1674 my $dim = shift;
1675 my $size_ref = shift;
1676 $dim++;
1677 if ( defined $arr and ref($arr) eq 'ARRAY' ) {
1678 push( @{$size_ref}, scalar @{$arr} );
1679 ( $dim, $size_ref ) = get_dim( $arr->[0], $dim, $size_ref );
1680 }
1681 return ( $dim, $size_ref );
1682 }
1683 sub format_value {
1684 my $val = shift;
1685 if ( not defined $val or $val eq '' ) {
1686 return sprintf("%10s",$PsN::output_style).',';
1687 } else {
1688 $_ = $val;
1689 my $nodot = /.*\..*/ ? 0 : 1;
1690 $_ =~ s/\.//g;
1691 if ( /.*\D+.*/ or $nodot) {
1692 return sprintf("%10s",$val).',';
1693 } else {
1694 return sprintf("%10.5f",$val).',';
1695 }
1696 }
1697 }
1698 debug -> die( message => "No results_file defined" )
1699 unless ( defined $self -> {'results_file'} );
1700 open ( RES, ">".$self -> {'directory'}.'/'.$self -> {'results_file'} );
1701 if ( defined $self -> {'results'} ) {
1702 my @all_results = @{$self -> {'results'}};
1703 for ( my $i = 0; $i <= $#all_results; $i++ ) {
1704 if ( defined $all_results[$i]{'own'} ) {
1705 my @my_results = @{$all_results[$i]{'own'}};
1706 for ( my $j = 0; $j <= $#my_results; $j++ ) {
1707 # These size estimates include the problem and sub_problem dimensions:
1708 my ( $ldim, $lsize_ref ) = get_dim( $my_results[$j]{'labels'}, -1, [] );
1709 my ( $vdim, $vsize_ref ) = get_dim( $my_results[$j]{'values'}, -1, [] );
1710 print RES $my_results[$j]{'name'},"\n" if ( $vdim > 1 );
1711 if ( defined $my_results[$j]{'values'} and
1712 scalar @{$my_results[$j]{'values'}} >= 0 ) {
1713 my @values = @{$my_results[$j]{'values'}};
1714 my @labels;
1715 if ( defined $my_results[$j]{'labels'} and
1716 scalar @{$my_results[$j]{'labels'}} >= 0 ) {
1717 @labels = @{$my_results[$j]{'labels'}};
1718 }
1719 # Print Header Labels
1720 if ( $ldim == 0 ) {
1721 my $label = \@labels;
1722 print RES ','.format_value($label),"\n";
1723 } elsif ( $ldim == 2 ) {
1724 print RES ',';
1725 for ( my $n = 0; $n < scalar @{$labels[1]}; $n++ ) {
1726 my $label = $labels[1][$n];
1727 print RES format_value($label);
1728 }
1729 print RES "\n" if ( scalar @{$labels[1]} );
1730 }
1731 # Print the values:
1732 if ( $vdim == 0 ) {
1733 print RES ','.format_value(\@values),"\n";
1734 } elsif ( $vdim == 1 ) {
1735 for ( my $m = 0; $m < scalar @{\@values}; $m++ ) {
1736 my $label = $labels[$m];
1737 print RES ','.format_value($label);
1738 my $val = $values[$m];
1739 print RES ','.format_value($val),"\n";
1740 }
1741 } elsif ( $vdim == 2 ) {
1742 for ( my $m = 0; $m < scalar @{\@values}; $m++ ) {
1743 my $label;
1744 if ( $ldim == 1 ) {
1745 $label = $labels[$m];
1746 } elsif ( $ldim == 2 ) {
1747 $label = $labels[0][$m];
1748 }
1749 print RES format_value($label);
1750 if ( defined $values[$m] ) {
1751 for ( my $n = 0; $n < scalar @{$values[$m]}; $n++ ) {
1752 print RES format_value($values[$m][$n]);
1753 }
1754 }
1755 print RES "\n";
1756 }
1757 }
1758 }
1759 }
1760 }
1761 }
1762 }
1763 close( RES );
1764 } else {
1765 debug -> warn( level => 2,
1766 message => "No subtools defined".
1767 ", using default printing routine" );
1768 }
1769 }
1770 end print_results
1771
1772 # }}}
1773
1774 # {{{ create_matlab_scripts
1775 start create_matlab_scripts
1776 {
1777 if( defined $PsN::lib_dir ){
1778 unless( -e $PsN::lib_dir . '/profiles.m') {
1779 'debug' -> die( message => 'LLP matlab template scripts are not installed, no matlab scripts will be generated.' );
1780 return;
1781 }
1782
1783 open( PROF, $PsN::lib_dir . '/profiles.m' );
1784 my @file = <PROF>;
1785 close( PROF );
1786 my $found_code;
1787 my $code_area_start=0;
1788 my $code_area_end=0;
1789
1790 for(my $i = 0;$i < scalar(@file); $i++) {
1791 if( $file[$i] =~ /% ---------- Autogenerated code below ----------/ ){
1792 $found_code = 1;
1793 $code_area_start = $i;
1794 }
1795 if( $file[$i] =~ /% ---------- End autogenerated code ----------/ ){
1796 unless( $found_code ){
1797 'debug' -> die ( message => 'LLP matlab template script is malformated, no matlab scripts will be generated' );
1798 return;
1799 }
1800 $code_area_end = $i;
1801 }
1802 }
1803
1804 my @auto_code;
1805 push( @auto_code, "str_format = '%30s';\n\n" );
1806
1807 my %param_names;
1808
1809 push( @auto_code, "col_names = [ " );
1810
1811 foreach my $param ( 'theta','omega','sigma' ) {
1812 my $labels = $self -> {'models'} -> [0] -> labels( parameter_type => $param );
1813 if ( defined $labels ){
1814 foreach my $label ( @{$labels -> [0]} ){
1815 push( @auto_code, " sprintf(str_format,'",$label,"');\n" );
1816 }
1817 }
1818 }
1819 push( @auto_code, " ];\n\n" );
1820 push( @auto_code, "goal = 3.84;\n\n" );
1821
1822 push( @auto_code, "filename = 'llplog1_no_header.csv';\n" );
1823
1824 splice( @file, $code_area_start, ($code_area_end - $code_area_start), @auto_code );
1825 open( OUTFILE, ">", $self -> {'directory'} . "/profiles.m" );
1826 print OUTFILE "addpath " . $PsN::lib_dir . ";\n";
1827 print OUTFILE @file ;
1828 close OUTFILE;
1829
1830 open( LOGFILE, "<", $self -> {'directory'} . "/llplog1.csv" );
1831 my @log = <LOGFILE>;
1832 close LOGFILE;
1833
1834 open( OUTFILE, ">", $self -> {'directory'} . "/llplog1_no_header.csv" );
1835 for( my $i = 1; $i <= $#log; $i ++ ){ #Skip header
1836 print OUTFILE $log[$i];
1837 }
1838 close OUTFILE;
1839
1840 } else {
1841 'debug' -> die( message => 'matlab_dir not configured, no matlab scripts will be generated.');
1842 return;
1843 }
1844 }
1845 end create_matlab_scripts
1846 # }}}
1847
1848 # {{{ create_R_scripts
1849 start create_R_scripts
1850 {
1851 unless( -e $PsN::lib_dir . '/R-scripts/llp.R' ){
1852 'debug' -> die( message => 'LLP R-script are not installed, no R scripts will be generated.' );
1853 return;
1854 }
1855 cp ( $PsN::lib_dir . '/R-scripts/llp.R', $self -> {'directory'} );
1856 }
1857 end create_R_scripts
1858 # }}}
Perl-speaks-NONMEM CVS import