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