?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
*** empty log message ***
[PsN.git] / lib / tool / cdd_subs.pm
blobfb7d383244add59b3d254f9fa19d1483808bb4cd
1 # {{{ include statements
2
3 start include statements
4 use Cwd;
5 use File::Copy 'cp';
6 use tool::llp;
7 use tool::modelfit;
8 use Math::Random;
9 use ext::Math::MatrixReal;
10 use Data::Dumper;
11 if ( $PsN::config -> {'_'} -> {'use_database'} ) {
12 # Testing DBD::mysql:
13 require DBI;
14 }
15 end include statements
16
17 # }}} include statements
18
19 # {{{ new
20
21 start new
22 {
23 foreach my $attribute ( 'logfile', 'raw_results_file' ) {
24 if ( not( ref($this -> {$attribute}) eq 'ARRAY' or
25 ref($this -> {$attribute}) eq 'HASH' ) ) {
26 my $tmp = $this -> {$attribute};
27 if ( not defined $tmp and $attribute eq 'logfile' ) {
28 $tmp = 'cdd.log';
29 }
30 $this -> {$attribute} = [];
31 for ( my $i = 1; $i <= scalar @{$this -> {'models'}}; $i++ ) {
32 my $name = $tmp;
33 if ( $name =~ /\./ ) {
34 $name =~ s/\./$i\./;
35 } else {
36 $name = $name.$i;
37 }
38 my $ldir;
39 ( $ldir, $name ) =
40 OSspecific::absolute_path( $this -> {'directory'}, $name );
41 push ( @{$this -> {$attribute}}, $ldir.$name ) ;
42 }
43 }
44 }
45 if ( $PsN::config -> {'_'} -> {'use_database'} ) {
46 my( $found_log, $found_cdd_id ) = $this -> read_cdd_log;
47
48 $this -> register_cdd_in_database unless ( $found_cdd_id );
49
50 $this -> log_object unless ( $found_log and $found_cdd_id );
51 print "Found ",$this -> {'cdd_id'},"\n";
52 }
53 }
54 end new
55
56 # }}} new
57
58 # {{{ register_cdd_in_database
59
60 start register_cdd_in_database
61 {
62 if ( $PsN::config -> {'_'} -> {'use_database'} ) {
63 my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
64 ";databse=".$PsN::config -> {'_'} -> {'project'},
65 $PsN::config -> {'_'} -> {'user'},
66 $PsN::config -> {'_'} -> {'password'},
67 {'RaiseError' => 1});
68 my $sth;
69 # bins and case_column can be defined for more than one model. Skip
70 # registration of these for now.
71 # $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
72 # ".cdd ( tool_id, bins, case_column ) ".
73 # "VALUES (".$self -> {'tool_id'}.", '".$self -> {'bins'}.
74 # "', '".$self -> {'case_column'}."' )");
75 $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
76 ".cdd ( tool_id ) ".
77 "VALUES (".$self -> {'tool_id'}." )");
78 $sth -> execute;
79 $self -> {'cdd_id'} = $sth->{'mysql_insertid'};
80 $sth -> finish;
81 $dbh -> disconnect;
82 }
83 }
84 end register_cdd_in_database
85
86 # }}} register_cdd_in_database
87
88 # {{{ register_mfit_results
89
90 start register_mfit_results
91 {
92 if ( $PsN::config -> {'_'} -> {'use_database'} ) {
93 my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
94 ";databse=".$PsN::config -> {'_'} -> {'project'},
95 $PsN::config -> {'_'} -> {'user'},
96 $PsN::config -> {'_'} -> {'password'},
97 {'RaiseError' => 1});
98 $dbh -> do( "LOCK TABLES ".$PsN::config -> {'_'} -> {'project'}.
99 ".cdd_modelfit_results WRITE" );
100 my $sth = $dbh -> prepare( "SELECT MAX(cdd_modelfit_results_id)".
101 " FROM ".$PsN::config -> {'_'} -> {'project'}.
102 ".cdd_modelfit_results" );
103 $sth -> execute or debug -> die( message => $sth->errstr ) ;
104 my $select_arr = $sth -> fetchall_arrayref;
105 $first_res_id = defined $select_arr -> [0][0] ? ($select_arr -> [0][0] + 1) : 0;
106 $last_res_id = $first_res_id + $#cook_score;
107 $sth -> finish;
108
109 my $insert_values;
110 for( my $i = 0; $i <= $#cook_score; $i++ ) {
111 $insert_values = $insert_values."," if ( defined $insert_values );
112 $insert_values = $insert_values.
113 "('".$self -> {'cdd_id'}."', '".$self -> {'model_ids'}[$model_number-1].
114 "', '".$self -> {'prepared_model_ids'}[($model_number-1)*($#cook_score+1)+$i].
115 "', '".($i+1).
116 "','$cook_score[$i]', '$covariance_ratio[$i]', '$projections[$i][0]', '$projections[$i][01]', '$outside_n_sd[$i]' )";
117 }
118 $dbh -> do("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
119 ".cdd_modelfit_results ".
120 "( cdd_id, orig_model_id, model_id, ".
121 "bin, cook_score, covariance_ratio, ".
122 "pca_component_1, pca_component_2, outside_n_sd ) ".
123 "VALUES $insert_values");
124 $dbh -> do( "UNLOCK TABLES" );
125 $dbh -> disconnect;
126 }
127 }
128 end register_mfit_results
129
130 # }}} register_mfit_results
131
132 # {{{ read_cdd_log
133 start read_cdd_log
134 {
135 if( -e $self -> {'directory'}.'object.txt' ) {
136 $found_log = 1;
137 open( OLOG, '<'.$self -> {'directory'}.'object.txt' );
138 my @olog = <OLOG>;
139 my $str = "(";
140 for ( my $i = 1; $i < $#olog; $i++ ) {
141 $str = $str.$olog[$i];
142 }
143 $str = $str.")";
144 my %tmp = eval( $str );
145
146 if( exists $tmp{'cdd_id'} ) {
147 $self -> {'cdd_id'} = $tmp{'cdd_id'};
148 $found_cdd_id = 1;
149 }
150 close( OLOG );
151 }
152 }
153 end read_cdd_log
154 # }}} read_cdd_log
155
156 # {{{ llp_pre_fork_setup
157
158 start llp_pre_fork_setup
159 {
160 $self -> modelfit_pre_fork_setup;
161 }
162 end llp_pre_fork_setup
163
164 # }}} llp_pre_fork_setup
165
166 # {{{ modelfit_pre_fork_setup
167
168 start modelfit_pre_fork_setup
169 {
170 # These attributes can be given as a
171 # 1. A scalar : used for all models and problems
172 # 2. A 1-dim. array : specified per problem but same for all models
173 # 3. A 2-dim. array : specified per problem and model
174 my $bins = $self -> {'bins'};
175 # my $idxs = $self -> {'grouping_indexes'};
176 my $case_columns = $self -> {'case_columns'};
177
178 if ( defined $bins ) {
179 unless ( ref( \$bins ) eq 'SCALAR' or
180 ( ref( $bins ) eq 'ARRAY' and scalar @{$bins} > 0 ) ) {
181 debug -> die( message => "Attribute bins is ",
182 "defined as a ",ref( $bins ),
183 "and is neither a scalar or a non-zero size array." );
184 } elsif ( ref( \$bins ) eq 'SCALAR' ) {
185 my @mo_bins = ();
186 foreach my $model ( @{$self -> {'models'}} ) {
187 my @pr_bins = ();
188 foreach my $problem ( @{$model -> problems} ) {
189 push( @pr_bins, $bins );
190 }
191 push( @mo_bins, \@pr_bins );
192 }
193 $self -> {'bins'} = \@mo_bins;
194 } elsif ( ref( $bins ) eq 'ARRAY' ) {
195 unless ( ref( \$bins -> [0] ) eq 'SCALAR' or
196 ( ref( $bins -> [0] ) eq 'ARRAY' and scalar @{$bins -> [0]} > 0 ) ) {
197 debug -> die( message => "Attribute bins is ",
198 "defined as a ",ref( $bins -> [0] ),
199 "and is neither a scalar or a non-zero size array." );
200 } elsif ( ref(\$bins -> [0]) eq 'SCALAR' ) {
201 my @mo_bins = ();
202 foreach my $model ( @{$self -> {'models'}} ) {
203 push( @mo_bins, $bins );
204 }
205 $self -> {'bins'} = \@mo_bins;
206 }
207 }
208 } else {
209 my @mo_bins = ();
210 foreach my $model ( @{$self -> {'models'}} ) {
211 my @pr_bins = ();
212 foreach my $data ( @{$model -> datas} ) {
213 push( @pr_bins, $data -> count_ind );
214 }
215 push( @mo_bins, \@pr_bins );
216 }
217 $self -> {'bins'} = \@mo_bins;
218 }
219
220 if ( defined $case_columns ) {
221 if ( defined $case_columns ) {
222 unless ( ref( \$case_columns ) eq 'SCALAR' or
223 ( ref( $case_columns ) eq 'ARRAY' and scalar @{$case_columns} > 0 ) ) {
224 debug -> die( message => "Attribute case_columns is ",
225 "defined as a ",ref( $case_columns ),
226 "and is neither a scalar or a non-zero size array." );
227 } elsif ( ref( \$case_columns ) eq 'SCALAR' ) {
228 # SCALAR!
229 my @mo_case_columns = ();
230 foreach my $model ( @{$self -> {'models'}} ) {
231 my @pr_case_columns = ();
232 for( my $i = 1; $i <= scalar @{$model -> problems}; $i++ ) {
233 if ( not $case_columns =~ /^\d/ ) {
234 # STRING
235 my ( $junk, $column_position ) = $model ->
236 _get_option_val_pos( name => $case_columns,
237 record_name => 'input',
238 problem_numbers => [$i] );
239 # We assume that there is no duplicate column names
240 push ( @pr_case_columns, $column_position->[0][0] );
241 } else {
242 # NUMBER
243 push ( @pr_case_columns, $case_columns );
244 }
245 }
246 push( @mo_case_columns, \@pr_case_columns );
247 }
248 $self -> {'case_columns'} = \@mo_case_columns;
249 } elsif ( ref( $case_columns ) eq 'ARRAY' ) {
250 # ARRAY!
251 unless ( ref( \$case_columns -> [0] ) eq 'SCALAR' or
252 ( ref( $case_columns -> [0] ) eq 'ARRAY' and
253 scalar @{$case_columns -> [0]} > 0 ) ) {
254 debug -> die( message => "Second dimension of attribute case_columns is ",
255 "defined as a ",ref( $case_columns -> [0]),
256 "and is neither a scalar or a non-zero size array." );
257 } elsif ( ref(\$case_columns -> [0]) eq 'SCALAR' ) {
258 # ARRAY -> SCALAR!
259 my @mo_case_columns = ();
260 foreach my $model ( @{$self -> {'models'}} ) {
261 my @pr_case_columns = ();
262 for( my $i = 1; $i <= scalar @{$model -> problems}; $i++ ) {
263 if ( not $case_columns =~ /^\d/ ) {
264 # STRING
265 my ( $junk, $column_position ) = $model ->
266 _get_option_val_pos( name => $case_columns->[$i-1],
267 record_name => 'input',
268 problem_numbers => [$i] );
269 # We assume that there is no duplicate column names
270 push ( @pr_case_columns, $column_position->[0][0] );
271 } else {
272 # NUMBER
273 push ( @pr_case_columns, $case_columns -> [$i-1] );
274 }
275 }
276 push( @mo_case_columns, \@pr_case_columns );
277 }
278 $self -> {'case_columns'} = \@mo_case_columns;
279 } elsif ( ref($case_columns -> [0]) eq 'ARRAY' ) {
280 # ARRAY -> ARRAY!
281 my @mo_case_columns = ();
282 for( my $j = 0; $j < scalar @{$self -> {'models'}}; $j++ ) {
283 my @pr_case_columns = ();
284 for( my $i = 1; $i <= scalar @{$self -> {'models'} -> problems}; $i++ ) {
285 if ( not $case_columns =~ /^\d/ ) {
286 # STRING
287 my ( $junk, $column_position ) = $self -> {'models'} -> [$j] ->
288 _get_option_val_pos( name => $case_columns->[$j]->[$i-1],
289 record_name => 'input',
290 problem_numbers => [$i] );
291 # We assume that there is no duplicate column names
292 push ( @pr_case_columns, $column_position->[0][0] );
293 } else {
294 # NUMBER
295 push ( @pr_case_columns, $case_columns -> [$j]->[$i-1] );
296 }
297 }
298 push( @mo_case_columns, \@pr_case_columns );
299 }
300 $self -> {'case_columns'} = \@mo_case_columns;
301 }
302 }
303 } else {
304 debug -> die( message => "case_columns is not specified for model $model_number" );
305 }
306 }
307 }
308 end modelfit_pre_fork_setup
309
310 # }}} modelfit_pre_fork_setup
311
312 # {{{ modelfit_setup
313
314 start modelfit_setup
315 {
316 my $subm_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
317 $self -> {'threads'} -> [1]:$self -> {'threads'};
318 $self -> general_setup( model_number => $model_number,
319 class => 'tool::modelfit',
320 subm_threads => $subm_threads );
321 }
322 end modelfit_setup
323
324 # }}} modelfit_setup
325
326 # {{{ llp_setup
327 start llp_setup
328 {
329 my @subm_threads;
330 if (ref( $self -> {'threads'} ) eq 'ARRAY') {
331 @subm_threads = @{$self -> {'threads'}};
332 unshift(@subm_threads);
333 } else {
334 @subm_threads = ($self -> {'threads'});
335 }
336 $self -> general_setup( model_number => $model_number,
337 class => 'tool::llp',
338 subm_threads => \@subm_threads );
339 }
340 end llp_setup
341 # }}} llp_setup
342
343 # {{{ general_setup
344
345 start general_setup
346 {
347 # Sub tool threads can be given as scalar or reference to an array?
348 my $subm_threads = $parm{'subm_threads'};
349 my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
350 $self -> {'threads'} -> [0]:$self -> {'threads'};
351 # case_column names are matched in the model, not the data!
352
353 my $model = $self -> {'models'} -> [$model_number-1];
354 my @bins = @{$self -> {'bins'} -> [$model_number-1]};
355
356 # Check which models that hasn't been run and run them
357 # This will be performed each step but will only result in running
358 # models at the first step, if at all.
359
360 # If more than one process is used, there is a VERY high risk of interaction
361 # between the processes when creating directories for model fits. Therefore
362 # the {'directory'} attribute is given explicitly below.
363
364
365 unless ( $model -> is_run ) {
366
367 # ----------------------- Run original run ------------------------------
368
369 # {{{ orig run
370
371 my %subargs = ();
372 if ( defined $self -> {'subtool_arguments'} ) {
373 %subargs = %{$self -> {'subtool_arguments'}};
374 }
375 if( $self -> {'nonparametric_etas'} or
376 $self -> {'nonparametric_marginals'} ) {
377 $model -> add_nonparametric_code;
378 }
379
380 my $orig_fit = tool::modelfit ->
381 new( reference_object => $self,
382 models => [$model],
383 threads => $subm_threads,
384 directory => $self -> {'directory'}.'/orig_modelfit_dir'.$model_number,
385 subtools => [],
386 parent_threads => $own_threads,
387 parent_tool_id => $self -> {'tool_id'},
388 logfile => undef,
389 raw_results => undef,
390 prepared_models => undef,
391 %subargs );
392
393 # $Data::Dumper::Maxdepth=1;
394 # die Dumper $orig_fit;
395 # ( models => [$model],
396
397 # run_on_lsf => $self -> {'run_on_lsf'},
398 # no_remote_execution => $self -> {'no_remote_execution'},
399 # no_remote_compile => $self -> {'no_remote_compile'},
400 # lsf_queue => $self -> {'lsf_queue'},
401 # lsf_options => $self -> {'lsf_options'},
402 # lsf_job_name => $self -> {'lsf_job_name'},
403 # lsf_project_name => $self -> {'lsf_project_name'},
404
405 # run_on_nordugrid => $self -> {'run_on_nordugrid'},
406 # cpu_time => $self -> {'cpu_time'},
407
408 # parent_tool_id => $self -> {'tool_id'},
409
410 # subtools => [],
411 # nm_version => $self -> {'nm_version'},
412 # picky => $self -> {'picky'},
413 # compress => $self -> {'compress'},
414 # threads => $subm_threads,
415 # retries => $self -> {'retries'},
416 # remove_temp_files => $self -> {'remove_temp_files'},
417 # base_directory => $self -> {'directory'},
418 # directory => $self -> {'directory'}.'/orig_modelfit_dir'.$model_number,
419 # parent_threads => $own_threads );
420
421 ui -> print( category => 'cdd',
422 message => 'Executing base model.' );
423
424
425 $orig_fit -> run;
426
427 # }}} orig run
428
429 }
430
431 # ------------------------ Print a log-header -----------------------------
432
433 # {{{ log header
434
435 open( LOG, ">>".$self -> {'logfile'}[$model_number-1] );
436 my $ui_text = sprintf("%-5s",'RUN').','.sprintf("%20s",'FILENAME ').',';
437 print LOG sprintf("%-5s",'RUN'),',',sprintf("%20s",'FILENAME '),',';
438 foreach my $param ( 'ofv', 'theta', 'omega', 'sigma' ) {
439 my $accessor = $param eq 'ofv' ? $param : $param.'s';
440 my $orig_ests = $model -> outputs -> [0] -> $accessor;
441 # Loop the problems
442 if( defined $orig_ests ){
443 for ( my $j = 0; $j < scalar @{$orig_ests}; $j++ ) {
444 if ( $param eq 'ofv' ) {
445 my $label = uc($param)."_".($j+1);
446 $ui_text = $ui_text.sprintf("%12s",$label).',';
447 print LOG sprintf("%12s",$label),',';
448 } else {
449 # Loop the parameter numbers (skip sub problem level)
450 if( defined $orig_ests -> [$j] and
451 defined $orig_ests -> [$j][0] ){
452 for ( my $num = 1; $num <= scalar @{$orig_ests -> [$j][0]}; $num++ ) {
453 my $label = uc($param).$num."_".($j+1);
454 $ui_text = $ui_text.sprintf("%12s",$label).',';
455 print LOG sprintf("%12s",$label),',';
456 }
457 }
458 }
459 }
460 }
461 }
462
463 print LOG "\n";
464
465 # }}} log header
466
467 # ------------------------ Log original run -------------------------------
468
469 # {{{ log orig
470
471 open( LOG, ">>".$self -> {'logfile'}[$model_number-1] );
472 $ui_text = sprintf("%5s",'0').','.sprintf("%20s",$model -> filename).',';
473 print LOG sprintf("%5s",'0'),',',sprintf("%20s",$model -> filename),',';
474 foreach my $param ( 'ofv', 'theta', 'omega', 'sigma' ) {
475 my $accessor = $param eq 'ofv' ? $param : $param.'s';
476 my $orig_ests = $model -> outputs -> [0] -> $accessor;
477 # Loop the problems
478 if( defined $orig_ests ) {
479 for ( my $j = 0; $j < scalar @{$orig_ests}; $j++ ) {
480 if ( $param eq 'ofv' ) {
481 $ui_text = $ui_text.sprintf("%12f",$orig_ests -> [$j][0]).',';
482 print LOG sprintf("%12f",$orig_ests -> [$j][0]),',';
483 } else {
484 # Loop the parameter numbers (skip sub problem level)
485 if( defined $orig_ests -> [$j] and
486 defined $orig_ests -> [$j][0] ){
487 for ( my $num = 0; $num < scalar @{$orig_ests -> [$j][0]}; $num++ ) {
488 $ui_text = $ui_text.sprintf("%12f",$orig_ests -> [$j][0][$num]).',';
489 print LOG sprintf("%12f",$orig_ests -> [$j][0][$num]),',';
490 }
491 }
492 }
493 }
494 }
495 }
496
497 print LOG "\n";
498
499 # }}} log orig
500
501 # --------------------- Initiate some variables ----------------------------
502
503 # {{{ inits
504
505 my $case_column = $self -> {'case_columns'} -> [$model_number-1] -> [0];
506
507 # Case-deletion Diagnostics will only work for models with one problem.
508 my $orig_data = $model -> datas -> [0];
509
510 if ( not defined $orig_data ) {
511 debug -> die( message => "No data file to resample from found in ".$model -> full_name );
512 }
513
514 my @problems = @{$model -> problems};
515 my @new_models = ();
516
517 my ( @skipped_ids, @skipped_keys, @skipped_values );
518
519 my %orig_factors = %{$orig_data -> factors( column => $case_column )};
520 my $maxbins = scalar keys %orig_factors;
521 my $pr_bins = ( defined $bins[0] and $bins[0] <= $maxbins ) ?
522 $bins[0] : $maxbins;
523
524 my $done = ( -e $self -> {'directory'}."/m$model_number/done" ) ? 1 : 0;
525
526 my ( @seed, $new_datas, $skip_ids, $skip_keys, $skip_values, $remainders );
527
528 # }}} inits
529
530 if ( not $done ) {
531
532 # -------------- Create new case-deleted data sets ----------------------
533
534 # {{{ create new
535
536 # Keep the new sample data objects i memory and write them to disk later
537 # with a proper name.
538
539 ( $new_datas, $skip_ids, $skip_keys, $skip_values, $remainders )
540 = $orig_data -> case_deletion( case_column => $case_column,
541 selection => $self -> {'selection_method'},
542 bins => $pr_bins,
543 target => 'mem' );
544 my $ndatas = scalar @{$new_datas};
545 open( DB, ">".$self -> {'directory'}."m$model_number/done.database.models" );
546 my @model_ids;
547 for ( my $j = 1; $j <= $ndatas; $j++ ) {
548 my @names = ( 'cdd_'.$j, 'rem_'.$j );
549 my @datasets = ( $new_datas -> [$j-1], $remainders -> [$j-1] );
550 foreach my $i ( 0, 1 ) {
551 my $dataset = $datasets[$i];
552 my ($data_dir, $data_file) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.$model_number,
553 $names[$i].'.dta' );
554 $dataset -> directory( $data_dir );
555 $dataset -> filename( $data_file );
556 $dataset -> flush;
557 my $newmodel = $model -> copy( filename => $data_dir.$names[$i].'.mod',
558 copy_data => 0,
559 copy_output => 0);
560 $newmodel -> ignore_missing_files(1);
561 $newmodel -> datafiles( new_names => [$names[$i].'.dta'] );
562 $newmodel -> outputfile( $data_dir.$names[$i].".lst" );
563 $newmodel -> datas -> [0] = $dataset;
564 if( $i == 1 ) {
565 # set MAXEVAL=0. Again, CDD will only work for one $PROBLEM
566 $newmodel -> maxeval( new_values => [[0]] );
567 }
568
569 if( $self -> {'nonparametric_etas'} or
570 $self -> {'nonparametric_marginals'} ) {
571 $newmodel -> add_nonparametric_code;
572 }
573
574 $newmodel -> _write;
575 push( @{$new_models[$i]}, $newmodel );
576 $model_ids[$j*$i+$j-1] = $newmodel -> register_in_database( force => 1 ) ;
577 print DB $model_ids[$j*$i+$j-1],"\n";
578 $self -> {'prepared_model_ids'}[($model_number-1)*$ndatas*2+$j*$i+$j-1] =
579 $model_ids[$j*$i+$j-1];
580 }
581 }
582
583 # my $new_data = $new_datas -> [$j-1];
584 # my ($data_dir, $data_file) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.$model_number,
585 # 'cdd_'.$j.'.dta' );
586 # $new_data -> directory( $data_dir );
587 # $new_data -> filename( $data_file );
588 # $new_data -> flush;
589 # my $newmodel = $model -> copy( filename => $data_dir."cdd_$j.mod",
590 # copy_data => 0,
591 # copy_output => 0);
592 # $newmodel -> ignore_missing_files(1);
593 # $newmodel -> datafiles( new_names => ["cdd_$j.dta"] );
594 # $newmodel -> outputfile( $data_dir."cdd_$j.lst" );
595 # $newmodel -> datas -> [0] = $new_data;
596 # $newmodel -> _write;
597 # push( @new_models, $newmodel );
598 # $model_ids[$j-1] = $newmodel -> register_in_database( force => 1 );
599 # print DB $model_ids[$j-1],"\n";
600 # $self -> {'prepared_model_ids'}[($model_number-1)*$ndatas+$j-1] =
601 # $model_ids[$j-1];
602 # }
603 close( DB );
604 if ( not -e $self -> {'directory'}."m$model_number/done.database.tool_models" ) {
605 $self -> register_tm_relation( model_ids => \@model_ids,
606 prepared_models => 1 );
607 open( DB, ">".$self -> {'directory'}."m$model_number/done.database.tool_models" );
608 print DB "";
609 close( DB );
610 }
611 # Create a checkpoint. Log the samples and individuals.
612 open( DONE, ">".$self -> {'directory'}."/m$model_number/done" ) ;
613 print DONE "Sampling from ",$orig_data -> filename, " performed\n";
614 print DONE "$pr_bins bins\n";
615 print DONE "Skipped individuals:\n";
616 for( my $k = 0; $k < scalar @{$skip_ids}; $k++ ) {
617 print DONE join(',',@{$skip_ids -> [$k]}),"\n";
618 }
619 print DONE "Skipped keys:\n";
620 for( my $k = 0; $k < scalar @{$skip_keys}; $k++ ) {
621 print DONE join(',',@{$skip_keys -> [$k]}),"\n";
622 }
623 print DONE "Skipped values:\n";
624 for( my $k = 0; $k < scalar @{$skip_values}; $k++ ) {
625 print DONE join(',',@{$skip_values -> [$k]}),"\n";
626 }
627 @seed = random_get_seed;
628 print DONE "seed: @seed\n";
629 close( DONE );
630
631 open( SKIP, ">".$self -> {'directory'}."skipped_individuals".$model_number.".csv" ) ;
632 for( my $k = 0; $k < scalar @{$skip_ids}; $k++ ) {
633 print SKIP join(',',@{$skip_ids -> [$k]}),"\n";
634 }
635 close( SKIP );
636 open( SKIP, ">".$self -> {'directory'}."skipped_keys".$model_number.".csv" ) ;
637 for( my $k = 0; $k < scalar @{$skip_keys}; $k++ ) {
638 print SKIP join(',',@{$skip_keys -> [$k]}),"\n";
639 }
640 close( SKIP );
641
642 # }}} create new
643
644 } else {
645
646 # --------- Recreate the datasets and models from a checkpoint ----------
647
648 # {{{ resume
649
650 ui -> print( category => 'cdd',
651 message => "Recreating models from a previous run" );
652 open( DONE, $self -> {'directory'}."/m$model_number/done" );
653 my @rows = <DONE>;
654 close( DONE );
655 my ( $junk, $junk, $stored_filename, $junk ) = split(' ',$rows[0],4);
656 my ( $stored_bins, $junk ) = split(' ',$rows[1],2);
657 my ( @stored_ids, @stored_keys, @stored_values );
658 for ( my $k = 3; $k < 3+$stored_bins; $k++ ) {
659 chomp($rows[$k]);
660 my @bin_ids = split(',', $rows[$k] );
661 push( @stored_ids, \@bin_ids );
662 }
663 for ( my $k = 4+$stored_bins; $k < 4+2*$stored_bins; $k++ ) {
664 chomp($rows[$k]);
665 my @bin_keys = split(',', $rows[$k] );
666 push( @stored_keys, \@bin_keys );
667 }
668 for ( my $k = 5+2*$stored_bins; $k < 5+3*$stored_bins; $k++ ) {
669 chomp($rows[$k]);
670 my @bin_values = split(',', $rows[$k] );
671 push( @stored_values, \@bin_values );
672 }
673 @seed = split(' ',$rows[5+3*$stored_bins]);
674 $skip_ids = \@stored_ids;
675 $skip_keys = \@stored_keys;
676 $skip_values = \@stored_values;
677 shift( @seed ); # get rid of 'seed'-word
678
679 # Reinitiate the model objects
680 my @model_ids;
681 my $reg_relations = 0;
682 if ( -e $self -> {'directory'}."m$model_number/done.database.models" ) {
683 open( DB, $self -> {'directory'}."m$model_number/done.database.models" );
684 @model_ids = <DB>;
685 chomp( @model_ids );
686 } else {
687 open( DB, ">".$self -> {'directory'}."m$model_number/done.database.models" );
688 $reg_relations = 1;
689 }
690 for ( my $j = 1; $j <= $stored_bins; $j++ ) {
691 my @names = ( 'cdd_'.$j, 'rem_'.$j );
692 foreach my $i ( 0, 1 ) {
693 my ($model_dir, $filename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
694 $model_number,
695 $names[$i].'.mod' );
696 my ($out_dir, $outfilename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
697 $model_number,
698 $names[$i].'.lst' );
699 my $new_mod = model ->
700 new( directory => $model_dir,
701 filename => $filename,
702 outputfile => $outfilename,
703 extra_files => $model -> extra_files,
704 target => 'disk',
705 ignore_missing_files => 1,
706 quick_reload => 1,
707 model_id => $model_ids[$j*$i+$j-1] );
708 push( @{$new_models[$i]}, $new_mod );
709 my $model_id;
710 if( not defined $model_ids[$j*$i+$j-1] ) {
711 $model_ids[$j*$i+$j-1] = $new_mod -> register_in_database;
712 print DB $model_ids[$j-1],"\n";
713 }
714 $self -> {'prepared_model_ids'}[($model_number-1)*
715 $stored_bins+$j*$i+$j-1] =
716 $model_ids[$j*$i+$j-1];
717 }
718
719 # my ($model_dir, $filename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
720 # $model_number,
721 # 'cdd_'.$j.'.mod' );
722 # my ($out_dir, $outfilename) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.
723 # $model_number,
724 # '/cdd_'.$j.'.lst' );
725 # my $new_mod = model ->
726 # new( directory => $model_dir,
727 # filename => $filename,
728 # outputfile => $outfilename,
729 # extra_files => $model -> extra_files,
730 # target => 'disk',
731 # ignore_missing_files => 1,
732 # quick_reload => 1,
733 # model_id => $model_ids[$j-1] );
734 # push( @new_models, $new_mod );
735 # my $model_id;
736 # if( not defined $model_ids[$j-1] ) {
737 # $model_ids[$j-1] = $new_mod -> register_in_database;
738 # print DB $model_ids[$j-1],"\n";
739 # }
740 # $self -> {'prepared_model_ids'}[($model_number-1)*$stored_bins+$j-1] =
741 # $model_ids[$j-1];
742 my $nl = $j == $stored_bins ? "" : "\r";
743 ui -> print( category => 'cdd',
744 message => ui -> status_bar( sofar => $j+1,
745 goal => $stored_bins+1 ).$nl,
746 wrap => 0,
747 newline => 0 );
748 }
749 close( DB );
750 if ( not -e $self -> {'directory'}."m$model_number/done.database.tool_models" ) {
751 $self -> register_tm_relation( model_ids => \@model_ids,
752 prepared_models => 1 ) if ( $reg_relations );
753 open( DB, ">".$self -> {'directory'}."m$model_number/done.database.tool_models" );
754 print DB "";
755 close( DB );
756 }
757 ui -> print( category => 'cdd',
758 message => " ... done." );
759 random_set_seed( @seed );
760 ui -> print( category => 'cdd',
761 message => "Using $stored_bins previously sampled case-deletion sets ".
762 "from $stored_filename" )
763 unless $self -> {'parent_threads'} > 1;
764
765 # }}} resume
766
767 }
768 push( @skipped_ids, $skip_ids );
769 push( @skipped_keys, $skip_keys );
770 push( @skipped_values, $skip_values );
771
772 # Use only the first half (the case-deleted) of the data sets.
773 $self -> {'prepared_models'}[$model_number-1]{'own'} = $new_models[0];
774
775 # The remainders are left until the analyze step
776 $self -> {'prediction_models'}[$model_number-1]{'own'} = $new_models[1];
777
778 # --------------------- Create the sub tools ------------------------------
779
780 # {{{ sub tools
781
782 my $subdir = $class;
783 $subdir =~ s/tool:://;
784 my @subtools = @{$self -> {'subtools'}};
785 shift( @subtools );
786 my %subargs = ();
787 if ( defined $self -> {'subtool_arguments'} ) {
788 %subargs = %{$self -> {'subtool_arguments'}};
789 }
790 push( @{$self -> {'tools'}},
791 $class ->
792 new( reference_object => $self,
793 models => $new_models[0],
794 threads => $subm_threads,
795 directory => $self -> {'directory'}.'/'.$subdir.'_dir'.$model_number,
796 _raw_results_callback => $self ->
797 _modelfit_raw_results_callback( model_number => $model_number ),
798 subtools => \@subtools,
799 parent_threads => $own_threads,
800 parent_tool_id => $self -> {'tool_id'},
801 logfile => undef,
802 raw_results => undef,
803 prepared_models => undef,
804 %subargs ) );
805
806
807 # ( models => $new_models[0],
808 # nm_version => $self -> {'nm_version'},
809 # run_on_lsf => $self -> {'run_on_lsf'},
810 # no_remote_execution => $self -> {'no_remote_execution'},
811 # no_remote_compile => $self -> {'no_remote_compile'},
812 # lsf_queue => $self -> {'lsf_queue'},
813 # lsf_options => $self -> {'lsf_options'},
814 # lsf_job_name => $self -> {'lsf_job_name'},
815 # lsf_project_name => $self -> {'lsf_project_name'},
816
817 # run_on_nordugrid => $self -> {'run_on_nordugrid'},
818 # cpu_time => $self -> {'cpu_time'},
819
820 # parent_tool_id => $self -> {'tool_id'},
821
822 # picky => $self -> {'picky'},
823 # compress => $self -> {'compress'},
824 # threads => $subm_threads,
825 # retries => $self -> {'retries'},
826 # remove_temp_files => $self -> {'remove_temp_files'},
827 # base_directory => $self -> {'directory'},
828 # directory => $self -> {'directory'}.'/'.$subdir.'_dir'.$model_number,
829 # _raw_results_callback => $self ->
830 # _modelfit_raw_results_callback( model_number => $model_number ),
831 # subtools => \@subtools,
832 # parent_threads => $own_threads,
833 # %subargs ) );
834
835 # }}} sub tools
836
837 open( SKIP, ">".$self -> {'directory'}."skipped_values".$model_number.".csv" ) ;
838 for( my $k = 0; $k < scalar @{$skip_values}; $k++ ) {
839 print SKIP join(',',@{$skip_values -> [$k]}),"\n";
840 }
841 close( SKIP );
842
843 }
844 end general_setup
845
846 # }}} general_setup
847
848 # {{{ llp_analyze
849
850 start llp_analyze
851 {
852 print "POSTFORK\n";
853 my %proc_results;
854 # $proc_results{'skipped.section.identifiers'} = $self -> {'skipped.section.identifiers'};
855 # $proc_results{'skipped_ids'} = $self -> {'skipped_ids'};
856 # $proc_results{'skipped_keys'} = $self -> {'skipped_keys'};
857 # $proc_results{'skipped_values'} = $self -> {'skipped_values'};
858
859 push( @{$self -> {'results'} -> {'own'}}, \%proc_results );
860 }
861 end llp_analyze
862
863 # }}} llp_analyze
864
865 # {{{ _modelfit_raw_results_callback
866
867 start _modelfit_raw_results_callback
868 {
869
870 # Use the cdd's raw_results file.
871 # The cdd and the bootstrap's callback methods are identical
872 # in the beginning, then the cdd callback adds cook.scores and
873 # cov.ratios.
874
875 my ($dir,$file) =
876 OSspecific::absolute_path( $self -> {'directory'},
877 $self -> {'raw_results_file'}[$model_number-1] );
878 my $orig_mod = $self -> {'models'}[$model_number-1];
879 my $xv = $self -> {'cross_validate'};
880 $subroutine = sub {
881 my $modelfit = shift;
882 my $mh_ref = shift;
883 my %max_hash = %{$mh_ref};
884 $modelfit -> raw_results_file( $dir.$file );
885 if( $cross_validation_set ) {
886 $modelfit -> raw_results_append( 1 ) if( not $self -> {'bca_mode'} ); # overwrite when doing a jackknife
887 my ( @new_header, %param_names );
888 foreach my $row ( @{$modelfit -> {'raw_results'}} ) {
889 unshift( @{$row}, 'cross_validation' );
890 }
891 $modelfit -> {'raw_results_header'} = undef; # May be a bit silly to do...
892 } else {
893 # my @orig_res = ( 0, 1, 1 ); # Model zero: original run. Problem 1 and subproblem 1
894 # foreach my $param ( @{$modelfit -> raw_results_header} ){
895 # next if( $param eq 'model' or $param eq 'problem' or $param eq 'subproblem' );
896 # my ( $accessor, $res );
897 # if ( $param eq 'theta' or $param eq 'omega' or $param eq 'sigma' or
898 # $param eq 'setheta' or $param eq 'seomega' or $param eq 'sesigma' or
899 # $param eq 'npomega' ) {
900
901 # $accessor = $param.'s';
902 # $res = $orig_mod -> {'outputs'} -> [0] -> $accessor;
903
904 # } elsif ( $param eq 'shrinkage_etas' ) {
905
906 # # Shrinkage does not work for subproblems right now.
907 # $res = $orig_mod -> eta_shrinkage;
908
909 # } elsif ( $param eq 'shrinkage_wres' ) {
910
911 # $res = $orig_mod -> wres_shrinkage;
912
913 # } else {
914
915 # $res = $orig_mod -> {'outputs'} -> [0] -> $param;
916
917 # }
918
919 # # {{{ Push results
920
921 # $res = defined $res ? $res -> [0][0] : undef;
922 # if( defined $res ) {
923 # if ( ref $res eq 'ARRAY' ) {
924 # push( @orig_res, @{$res} );
925 # push( @orig_res, (undef) x ($max_hash{$param}- scalar @{$res}) );
926 # } else {
927 # push( @orig_res, $res );
928 # }
929 # } else {
930 # push( @orig_res, (undef) x $max_hash{$param} );
931 # }
932
933 # # }}} Push results
934
935 # }
936
937 my ($raw_results_row,$row_structure) = $self -> create_raw_results_rows( max_hash => $mh_ref,
938 model => $orig_mod );
939
940 unshift( @{$modelfit -> {'raw_results'}}, @{$raw_results_row} );
941
942 &{$self -> {'_raw_results_callback'}}( $self, $modelfit )
943 if ( defined $self -> {'_raw_results_callback'} );
944
945 my ( @new_header, %param_names );
946 my @params = ( 'theta', 'omega', 'sigma' );
947 foreach my $param ( @params ) {
948 my $labels = $orig_mod -> labels( parameter_type => $param );
949 $param_names{$param} = $labels -> [0] if ( defined $labels );
950 }
951
952 foreach my $name ( @{$modelfit -> {'raw_results_header'}} ) {
953 my $success;
954 foreach my $param ( @params ) {
955 if ( $name eq $param ){
956 if( defined $param_names{$name} ) {
957 push( @new_header , @{$param_names{$name}} );
958 } else {
959 for ( my $i = 1; $i <= $max_hash{ $name }; $i++ ) {
960 push ( @new_header, uc(substr($name,0,2)).$i );
961 }
962 }
963 $success = 1;
964 last;
965 } elsif ( $name eq 'se'.$param ) {
966 my @new_names = ();
967 for ( my $i = 1; $i <= $max_hash{ $name }; $i++ ) {
968 push ( @new_names, uc(substr($param,0,2)).$i );
969 }
970 map ( $_ = 'se'.$_, @new_names );
971
972 push( @new_header, @new_names );
973 $success = 1;
974 last;
975 }
976 }
977 unless( $success ){
978 push( @new_header, $name );
979 }
980 }
981
982 $modelfit -> {'raw_results_header'} = \@new_header;
983
984 if( $xv and not $self -> {'bca_mode'} ) {
985 foreach my $row ( @{$modelfit -> {'raw_results'}} ) {
986 unshift( @{$row}, 'cdd' );
987 }
988 unshift( @{$modelfit -> {'raw_results_header'}}, 'method' );
989 }
990 }
991 };
992 return $subroutine;
993 }
994 end _modelfit_raw_results_callback
995
996 # }}} _modelfit_raw_results_callback
997
998 # {{{ modelfit_analyze
999
1000 start modelfit_analyze
1001 {
1002 # Only valid for one problem and one sub problem.
1003
1004 if ( $self -> {'cross_validate'} ) {
1005
1006 # --- Evaluate the models on the remainder data sets ----
1007
1008 # {{{ eval models
1009
1010 for ( my $i = 0;
1011 $i < scalar @{$self -> {'prediction_models'}[$model_number-1]{'own'}};
1012 $i++ ) {
1013 $self -> {'prediction_models'}[$model_number-1]{'own'}[$i] ->
1014 update_inits( from_model => $self ->
1015 {'prepared_models'}[$model_number-1]{'own'}[$i]);
1016 $self -> {'prediction_models'}[$model_number-1]{'own'}[$i] ->
1017 remove_records( type => 'covariance' );
1018 $self -> {'prediction_models'}[$model_number-1]{'own'}[$i] -> _write;
1019 }
1020 my ($dir,$file) =
1021 OSspecific::absolute_path( $self -> {'directory'},
1022 $self -> {'raw_results_file'}[$model_number-1] );
1023 my $mod_eval = tool::modelfit ->
1024 new( reference_object => $self,
1025 models => $self ->
1026 {'prediction_models'}[$model_number-1]{'own'},
1027 base_directory => $self -> {'directory'},
1028 directory => $self -> {'directory'}.
1029 'evaluation_dir'.$model_number,
1030 _raw_results_callback => $self ->
1031 _modelfit_raw_results_callback( model_number => $model_number,
1032 cross_validation_set => 1 ),
1033 parent_tool_id => $self -> {'tool_id'},
1034 logfile => undef,
1035 raw_results => undef,
1036 prepared_models => undef,
1037 retries => 1 );
1038 $Data::Dumper::Maxdepth = 2;
1039 # die Dumper $mod_eval;
1040 print "Running xv runs\n";
1041 $mod_eval -> run;
1042
1043 # }}} eval models
1044
1045 }
1046
1047 # ---------------------- Cook-score -----------------------
1048
1049 # {{{ Cook-score
1050
1051 my ( @cook_score, @cov_ratio );
1052 if( $self -> models -> [$model_number-1] ->
1053 outputs -> [0] -> covariance_step_successful -> [0][0]) {
1054
1055 ui -> print( category => 'cdd',
1056 message => "Calculating the cook-scores" );
1057 my @orig_ests;
1058 my @cdd_ests;
1059 my @changes;
1060 my @rel_changes;
1061
1062 my $output_harvest = $self -> harvest_output( accessors => ['thetas','raw_omegas','raw_sigmas'],
1063 search_output => 1 );
1064
1065 foreach my $param ( 'thetas', 'raw_omegas', 'raw_sigmas' ) {
1066 my $orig_est = $self -> models -> [$model_number-1] -> outputs -> [0] -> $param;
1067 my @mod_ests;
1068 my $est = defined $output_harvest -> {$param} ?
1069 $output_harvest -> {$param} -> [$model_number-1]{'own'} : [];
1070 if( defined $est ) {
1071 for ( my $i = 0; $i < scalar @{$est}; $i++ ) {
1072 my $n_par = defined $est->[$i][0][0] ?
1073 scalar @{$est->[$i][0][0]} : 0;
1074 if( defined $est->[$i][0][0] ) {
1075 push( @{$cdd_ests[$i]}, @{$est->[$i][0][0]} );
1076 for ( my $j = 0; $j < $n_par; $j++ ) {
1077 push( @{$changes[$i]}, $orig_est->[0][0][$j]-$est->[$i][0][0][$j]);
1078 }
1079 } else {
1080 push( @{$cdd_ests[$i]}, (undef) x $n_par );
1081 for ( my $j = 0; $j < $n_par; $j++ ) {
1082 push( @{$changes[$i]}, (undef) x $n_par);
1083 }
1084 }
1085 }
1086 }
1087 }
1088 my $inv_cov_linear = $self -> models -> [$model_number-1] ->
1089 outputs -> [0] -> raw_invcovmatrix -> [0][0];
1090
1091 # Equation: sqrt((orig_est-est(i))'*inv_cov_matrix*(orig_est-est(i)))
1092 for ( my $i = 0; $i <= $#changes; $i++ ) {
1093 # Rightmost multiplication:
1094 my @prod;
1095 my $m = 1;
1096 my $n = 1;
1097 foreach my $element ( @{$inv_cov_linear} ) {
1098 $element = 0 if ( $element eq '.........' );
1099 if ( $m == $n ) {
1100 # Diagonal
1101 $prod[$n-1] = $prod[$n-1] + $element * $changes[$i][$n-1];
1102 } else {
1103 # Upper triangle
1104 $prod[$n-1] = $prod[$n-1] + $element * $changes[$i][$m-1];
1105 # Lower triangle
1106 $prod[$m-1] = $prod[$m-1] + $element * $changes[$i][$n-1];
1107 }
1108 if ( $m == $n ) {
1109 $m++;
1110 $n = 1;
1111 } else {
1112 $n++;
1113 }
1114 }
1115 # Leftmost multiplication:
1116 for( my $j = 0; $j <= $#prod; $j++ ) {
1117 $cook_score[$i] = $cook_score[$i] + $changes[$i][$j] * $prod[$j];
1118 }
1119 my $nl = $i == $#changes ? "" : "\r";
1120 ui -> print( category => 'cdd',
1121 message => ui -> status_bar( sofar => $i+1,
1122 goal => $#changes+1 ).$nl,
1123 wrap => 0,
1124 newline => 0 );
1125 }
1126 # Calculate the square root
1127 for ( my $i = 0; $i <= $#cook_score; $i++ ) {
1128 if( $cook_score[$i] >= 0 ) {
1129 $cook_score[$i] = sqrt($cook_score[$i]);
1130 } else {
1131 open( LOG, ">>".$self -> {'logfile'}[$model_number-1] );
1132 print LOG "Negative squared cook-score ",$cook_score[$i],
1133 "; can't take the square root.\n",
1134 "The cook-score for model $model_number and cdd bin $i was set to zero\n";
1135 # "The cook-score for model $model_number and cdd bin $i was set to undef\n";
1136 close( LOG );
1137 $cook_score[$i] = 0;
1138 # $cook_score[$i] = undef;
1139 }
1140 }
1141 }
1142
1143 $self -> {'cook_scores'} = \@cook_score;
1144
1145 ui -> print( category => 'cdd',
1146 message => " ... done." );
1147
1148 # }}} Cook-score
1149
1150 # ------------------- Covariance Ratio --------------------
1151
1152 # {{{ Covariance Ratio
1153
1154 # {{{ sub clear dots
1155
1156 sub clear_dots {
1157 my $m_ref = shift;
1158 my @matrix = @{$m_ref};
1159 # get rid of '........'
1160 my @clear;
1161 foreach ( @matrix ) {
1162 push( @clear, $_ ) unless ( $_ eq '.........' );
1163 }
1164 # print Dumper \@clear;
1165 return \@clear;
1166 }
1167
1168 # }}}
1169
1170 # {{{ sub make square
1171
1172 sub make_square {
1173 my $m_ref = shift;
1174 my @matrix = @{$m_ref};
1175 # Make the matrix square:
1176 my $elements = scalar @matrix; # = M*(M+1)/2
1177 my $M = -0.5 + sqrt( 0.25 + 2 * $elements );
1178 my @square;
1179 for ( my $m = 1; $m <= $M; $m++ ) {
1180 for ( my $n = 1; $n <= $m; $n++ ) {
1181 push( @{$square[$m-1]}, $matrix[($m-1)*$m/2 + $n - 1] );
1182 unless ( $m == $n ) {
1183 push( @{$square[$n-1]}, $matrix[($m-1)*$m/2 + $n - 1] );
1184 }
1185 }
1186 }
1187 return \@square;
1188 }
1189
1190 # }}}
1191
1192 ui -> print( category => 'cdd',
1193 message => "Calculating the covariance-ratios" );
1194
1195 # Equation: sqrt(det(cov_matrix(i))/det(cov_matrix(orig)))
1196 my $cov_linear = $self -> models -> [$model_number-1] ->
1197 outputs -> [0] -> raw_covmatrix -> [0][0];
1198 my $orig_det;
1199 if( defined $cov_linear ) {
1200 my $orig_cov = Math::MatrixReal ->
1201 new_from_cols( make_square( clear_dots( $cov_linear ) ) );
1202 $orig_det = $orig_cov -> det();
1203 }
1204 # AUTOLOAD: raw_covmatrix
1205
1206 my $output_harvest = $self -> harvest_output( accessors => ['raw_covmatrix'],
1207 search_output => 1 );
1208
1209 my $est_cov = defined $output_harvest -> {'raw_covmatrix'} ? $output_harvest -> {'raw_covmatrix'} -> [$model_number-1]{'own'} : [];
1210
1211 my $mods = scalar @{$est_cov};
1212 for ( my $i = 0; $i < scalar @{$est_cov}; $i++ ) {
1213 if ( $orig_det != 0 and defined $est_cov->[$i][0][0] ) {
1214 my $cov = Math::MatrixReal ->
1215 new_from_cols( make_square( clear_dots( $est_cov->[$i][0][0] ) ) );
1216 my $ratio = $cov -> det() / $orig_det;
1217 if( $ratio > 0 ) {
1218 push( @cov_ratio, sqrt( $ratio ) );
1219 } else {
1220 open( LOG, ">>".$self -> {'logfile'}[$model_number-1] );
1221 print LOG "Negative covariance ratio ",$ratio,
1222 "; can't take the square root.\n",
1223 "The covariance ratio for model $model_number and cdd bin $i was set to one (1)\n";
1224 # "The covariance ratio for model $model_number and cdd bin $i was set to undef\n";
1225 close( LOG );
1226 push( @cov_ratio, 1 );
1227 # push( @cov_ratio, undef );
1228 }
1229 } else {
1230 open( LOG, ">>".$self -> {'logfile'}[$model_number-1] );
1231 print LOG "The determinant of the cov-matrix of the original run was zero\n",
1232 "or the determinant of cdd bin $i was undefined\n",
1233 "The covariance ratio for model $model_number and cdd bin $i was set to one (1)\n";
1234 # "The covariance ratio for model $model_number and cdd bin $i was set to undef\n";
1235 close( LOG );
1236 push( @cov_ratio, 1 );
1237 # push( @cov_ratio, undef );
1238 }
1239
1240 my $nl = $i == $mods-1 ? "" : "\r";
1241 ui -> print( category => 'cdd',
1242 message => ui -> status_bar( sofar => $i+1,
1243 goal => $mods ).$nl,
1244 wrap => 0,
1245 newline => 0 );
1246 }
1247
1248 $self -> {'covariance_ratios'} = \@cov_ratio;
1249
1250 ui -> print( category => 'cdd',
1251 message => " ... done." );
1252
1253 # }}} Covariance Ratio
1254
1255 # - Perform a PCA on the cook-score:covariance-ratio data --
1256
1257 my ( $eig_ref, $eig_vec_ref, $proj_ref, $std_ref ) =
1258 $self -> pca( data_matrix => [\@cook_score,\@cov_ratio] );
1259 my @projections = @{$proj_ref};
1260 my @standard_deviation = @{$std_ref};
1261
1262 # ---- Mark the runs with CS-CR outside N standard deviations of the PCA ----
1263
1264 # {{{ mark runs
1265
1266 my @outside_n_sd;
1267 for( my $i = 0; $i <= $#projections; $i++ ) {
1268 my $vector_length = 0;
1269 for( my $j = 0; $j <= 1; $j++ ) {
1270 $vector_length += $projections[$i][$j]**2;
1271 }
1272 $vector_length = sqrt( $vector_length );
1273 my $n_sd = 0;
1274 for( my $j = 0; $j <= 1; $j++ ) {
1275 $n_sd += (($projections[$i][$j]/$vector_length)*$standard_deviation[$j])**2;
1276 }
1277 $n_sd = ( $self -> {'outside_n_sd_check'} * sqrt( $n_sd ) );
1278 $outside_n_sd[$i] = $vector_length > $n_sd ? 1 : 0;
1279 }
1280
1281 $self -> {'outside_n_sd'} = \@outside_n_sd;
1282
1283 # }}} mark runs
1284
1285 my %covariance_return_section;
1286 $covariance_return_section{'name'} = 'Diagnostics';
1287 $covariance_return_section{'labels'} = [[],['cook.scores','covariance.ratios','outside.n.sd']];
1288
1289 my @res_array;
1290 for( my $i = 0; $i <= $#cov_ratio; $i ++ ){
1291 push( @res_array , [$cook_score[$i],$cov_ratio[$i],$outside_n_sd[$i]] );
1292 }
1293
1294 $covariance_return_section{'values'} = \@res_array;
1295
1296 push( @{$self -> {'results'}[$model_number-1]{'own'}},\%covariance_return_section );
1297
1298 # --------------- Relative estimate change ----------------
1299
1300 # {{{ Relative change of the parameter estimates
1301
1302 my $output_harvest = $self -> harvest_output( accessors => ['ofv', 'thetas', 'omegas', 'sigmas','sethetas', 'seomegas', 'sesigmas'],
1303 search_output => 1 );
1304
1305
1306 my %return_section;
1307 $return_section{'name'} = 'relative.changes';
1308 $return_section{'labels'} = [[],[]];
1309
1310 my @rel_ests;
1311
1312 for ( my $i = 0; $i < scalar @{$output_harvest -> {'ofv'} -> [$model_number-1]{'own'}}; $i++ ) {
1313 my @values;
1314 foreach my $param ( 'ofv', 'thetas', 'omegas', 'sigmas',
1315 'sethetas', 'seomegas', 'sesigmas',) {
1316
1317 my $orig_est = $self -> {'models'} -> [$model_number-1] -> outputs -> [0] -> $param;
1318 my $est = defined $output_harvest -> {$param} ? $output_harvest -> {$param} -> [$model_number-1]{'own'} : [];
1319
1320 if ( $param eq 'ofv' ) {
1321 if ( defined $orig_est->[0][0] and $orig_est->[0][0] != 0 ) {
1322 push( @values, ($orig_est->[0][0]-$est->[$i][0][0])/$orig_est->[0][0] );
1323 } else {
1324 push( @values, 'INF' );
1325 }
1326 if( $i == 0 ){
1327 push( @{$return_section{'labels'} -> [1]}, $param );
1328 }
1329 } else {
1330 my @in_rel_ests;
1331 if( defined $est->[$i][0][0] ){
1332 for ( my $j = 0; $j < scalar @{$est->[$i][0][0]}; $j++ ) {
1333 if ( defined $orig_est->[0][0][$j] and $orig_est->[0][0][$j] != 0 ) {
1334 push( @values, ($orig_est->[0][0][$j]-$est->[$i][0][0][$j])/$orig_est->[0][0][$j]*100);
1335 } else {
1336 push( @values, 'INF' );
1337 }
1338 if( $i == 0 ){
1339 push( @{$return_section{'labels'} -> [1]}, $param . '_' . $j );
1340 }
1341 }
1342 }
1343 }
1344
1345 }
1346 push( @rel_ests, \@values );
1347 }
1348
1349 $return_section{'values'} = \@rel_ests ;
1350 push( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
1351
1352 # }}} Relative change of the parameter estimates
1353
1354 $self -> update_raw_results(model_number => $model_number);
1355
1356 # ------------- Register the results in a Database ----------------
1357
1358 if( not -e $self -> {'directory'}."m$model_number/done.database.results" ) {
1359 open( DB, ">".$self -> {'directory'}."m$model_number/done.database.results" );
1360 my ( $start_id, $last_id ) = $self ->
1361 register_mfit_results( model_number => $model_number,
1362 cook_score => \@cook_score,
1363 covariance_ratio => \@cov_ratio,
1364 projections => $proj_ref,
1365 outside_n_sd => \@outside_n_sd );
1366 print DB "$start_id-$last_id\n";
1367 close( DB );
1368 }
1369
1370 # experimental: to save memory
1371 $self -> {'prepared_models'}[$model_number-1]{'own'} = undef;
1372 }
1373 end modelfit_analyze
1374
1375 # }}} modelfit_analyze
1376
1377 # {{{ pca
1378 start pca
1379
1380 my $D = Math::MatrixReal ->
1381 new_from_rows( \@data_matrix );
1382 my @n_dim = @{$data_matrix[0]};
1383 my @d_dim = @data_matrix;
1384 my $n = scalar @n_dim;
1385 my $d = scalar @d_dim;
1386 map( $_=(1/$n), @n_dim );
1387 my $frac_vec_n = Math::MatrixReal ->
1388 new_from_cols( [\@n_dim] );
1389 map( $_=1, @n_dim );
1390 map( $_=1, @d_dim );
1391 my $one_vec_n = Math::MatrixReal ->
1392 new_from_cols( [\@n_dim] );
1393 my $one_vec_d = Math::MatrixReal ->
1394 new_from_cols( [\@d_dim] );
1395 my $one_vec_d_n = $one_vec_d * ~$one_vec_n;
1396 my $M = $D*$frac_vec_n;
1397 my $M_matrix = $M * ~$one_vec_n;
1398
1399 # Calculate the mean-subtracted data
1400 my $S = $D-$M_matrix;
1401
1402 # compue the empirical covariance matrix
1403 my $C = $S * ~$S;
1404
1405 # compute the eigenvalues and vectors
1406 my ($l, $V) = $C -> sym_diagonalize();
1407
1408 # Project the original data on the eigenvectors
1409 my $P = ~$V * $S;
1410
1411
1412 # l, V and projections are all MatrixReal objects.
1413 # We need to return the normal perl equivalents.
1414 @eigenvalues = @{$l->[0]};
1415 @eigenvectors = @{$V->[0]};
1416 @std = @{$self -> std( data_matrix => $P -> [0] )};
1417 # Make $P a n * d matrix
1418 $P = ~$P;
1419 @projections = @{$P->[0]};
1420
1421 end pca
1422 # }}} pca
1423
1424 # {{{ std
1425 start std
1426
1427 my ( @sum, @pow_2_sum );
1428 if ( defined $data_matrix[0] ) {
1429 my $n = scalar @{$data_matrix[0]};
1430 for( my $i = 0; $i <= $#data_matrix; $i++ ) {
1431 for( my $j = 0; $j < $n; $j++ ) {
1432 $sum[$i] = $sum[$i]+$data_matrix[$i][$j];
1433 $pow_2_sum[$i] += $data_matrix[$i][$j]*$data_matrix[$i][$j];
1434 }
1435 $std[$i] = sqrt( ( $n*$pow_2_sum[$i] - $sum[$i]*$sum[$i] ) / ($n*$n) );
1436 }
1437 }
1438
1439 end std
1440 # }}} std
1441
1442 # {{{ modelfit_post_fork_analyze
1443
1444 start modelfit_post_fork_analyze
1445 {
1446 # my @modelfit_results = @{ $self -> {'results'} -> {'subtools'} };
1447 my @modelfit_results = @{ $self -> {'results'} };
1448
1449 ui -> print( category => 'cdd',
1450 message => "Soon done" );
1451 }
1452 end modelfit_post_fork_analyze
1453
1454 # }}} modelfit_post_fork_analyze
1455
1456 # {{{ modelfit_results
1457
1458 start modelfit_results
1459 {
1460 my @orig_models = @{$self -> {'models'}};
1461 my @orig_raw_results = ();
1462 foreach my $orig_model ( @orig_models ) {
1463 my $orig_output = $orig_model -> outputs -> [0];
1464 push( @orig_raw_results, $orig_output -> $accessor );
1465 }
1466 # my @models = @{$self -> {'prepared_models'}};
1467 my @outputs = @{$self -> {'results'}};
1468
1469 my @raw_results = ();
1470
1471 foreach my $mod ( @outputs ) {
1472 my @raw_inner = ();
1473 foreach my $output ( @{$mod -> {'subset_outputs'}} ) {
1474 push( @raw_inner, $output -> $accessor );
1475 }
1476 push( @raw_results, \@raw_inner );
1477 }
1478 if ( $format eq 'relative' or $format eq 'relative_percent' ) {
1479 @results = ();
1480 for ( my $i = 0; $i <= $#orig_raw_results; $i++ ) {
1481 print "Model\t$i\n";
1482 my @rel_subset = ();
1483 for ( my $i2 = 0; $i2 < scalar @{$raw_results[$i]}; $i2++ ) {
1484 print "Subset Model\t$i2\n";
1485 my @rel_prob = ();
1486 for ( my $j = 0; $j < scalar @{$orig_raw_results[$i]}; $j++ ) {
1487 print "Problem\t$j\n";
1488 if( ref( $orig_raw_results[$i][$j] ) eq 'ARRAY' ) {
1489 my @rel_subprob = ();
1490 for ( my $k = 0; $k < scalar @{$orig_raw_results[$i][$j]}; $k++ ) {
1491 print "Subprob\t$k\n";
1492 if( ref( $orig_raw_results[$i][$j][$k] ) eq 'ARRAY' ) {
1493 my @rel_instance = ();
1494 for ( my $l = 0; $l < scalar @{$orig_raw_results[$i][$j][$k]}; $l++ ) {
1495 print "Instance\t$l\n";
1496 my $orig = $orig_raw_results[$i][$j][$k][$l];
1497 my $res = $raw_results[$i][$i2][$j][$k][$l];
1498 if( defined $orig and ! $orig == 0 ) {
1499 print "ORIGINAL $orig\n";
1500 print "SUBSET $res\n";
1501 print "RELATIVE ",$res/$orig,"\n";
1502 if ( $format eq 'relative_percent' ) {
1503 push( @rel_instance, ($res/$orig-1)*100 );
1504 } else {
1505 push( @rel_instance, $res/$orig );
1506 }
1507 } else {
1508 push( @rel_instance, 'NA' );
1509 }
1510 push( @rel_subprob,\@rel_instance );
1511 }
1512 } elsif( ref( $orig_raw_results[$i][$j][$k] ) eq 'SCALAR' ) {
1513 print "One instance per problem\n";
1514 my $orig = $orig_raw_results[$i][$j][$k];
1515 my $res = $raw_results[$i][$i2][$j][$k];
1516 if( defined $orig and ! $orig == 0 ) {
1517 print "ORIGINAL $orig\n";
1518 print "SUBSET $res\n";
1519 print "RELATIVE ",$res/$orig,"\n";
1520 if ( $format eq 'relative_percent' ) {
1521 push( @rel_subprob, ($res/$orig-1)*100 );
1522 } else {
1523 push( @rel_subprob, $res/$orig );
1524 }
1525 } else {
1526 push( @rel_subprob, 'NA' );
1527 }
1528 } else {
1529 print "WARNING: tool::cdd -> modelfit_results: neither\n\t".
1530 "array or scalar reference found at layer 4 in result data\n\t".
1531 "structure (found ",ref( $orig_raw_results[$i][$j][$k] ),")\n";
1532 }
1533 }
1534 push( @rel_prob, \@rel_subprob );
1535 } elsif( ref( $orig_raw_results[$i][$j] ) eq 'SCALAR' ) {
1536 print "One instance per problem\n";
1537 my $orig = $orig_raw_results[$i][$j];
1538 my $res = $raw_results[$i][$i2][$j];
1539 if( defined $orig and ! $orig == 0 ) {
1540 print "ORIGINAL $orig\n";
1541 print "SUBSET $res\n";
1542 print "RELATIVE ",$res/$orig,"\n";
1543 if ( $format eq 'relative_percent' ) {
1544 push( @rel_prob, ($res/$orig-1)*100 );
1545 } else {
1546 push( @rel_prob, $res/$orig );
1547 }
1548 } else {
1549 push( @rel_prob, 'NA' );
1550 }
1551 } else {
1552 print "WARNING: tool::cdd -> modelfit_results: neither\n\t".
1553 "array or scalar reference found at layer 3 in result data\n\t".
1554 "structure (found ",ref( $orig_raw_results[$i][$j] ),")\n";
1555 }
1556 }
1557 push( @rel_subset, \@rel_prob );
1558 }
1559 push( @results, \@rel_subset );
1560 }
1561 } else {
1562 @results = @raw_results;
1563 }
1564 }
1565 end modelfit_results
1566
1567 # }}} modelfit_results
1568
1569 # {{{ relative_estimates
1570
1571 start relative_estimates
1572 {
1573 my $accessor = $parameter.'s';
1574 my @params = $self -> $accessor;
1575
1576 # print "Parameter: $parameter\n";
1577 # sub process_inner_results {
1578 # my $res_ref = shift;
1579 # my $pad = shift;
1580 # $pad++;
1581 # foreach my $res ( @{$res_ref} ) {
1582 # if ( ref ( $res ) eq 'ARRAY' ) {
1583 # process_inner_results( $res, $pad );
1584 # } else {
1585 # print "RELEST $pad\t$res\n";
1586 # }
1587 # }
1588 # }
1589 # process_inner_results( \@params, 0 );
1590
1591 my @orig_params = $self -> $accessor( original_models => 1 );
1592 # [?][model][prob][subp][#]
1593 # print "ORIG TH1: ",$orig_params[0][0][0][0][0],"\n";
1594 for ( my $i = 0; $i < scalar @params; $i++ ) {
1595 # Loop over models
1596 my @mod = ();
1597 for ( my $j = 0; $j < scalar @{$params[$i]}; $j++ ) {
1598 # Loop over data sets
1599 my @prep = ();
1600 for ( my $k = 1; $k < scalar @{$params[$i]->[$j]}; $k++ ) {
1601 # Loop over problems (sort of, at least)
1602 my @prob = ();
1603 for ( my $l = 0; $l < scalar @{$params[$i]->[$j]->[$k]}; $l++ ) {
1604 # Loop over sub problems (sort of, at least)
1605 my @sub = ();
1606 for ( my $m = 0; $m < scalar @{$params[$i]->[$j]->[$k]->[$l]}; $m++ ) {
1607 # Loop over the params
1608 my @par = ();
1609 for ( my $n = 0; $n < scalar @{$params[$i][$j][$k][$l][$m]}; $n++ ) {
1610 my $orig = $orig_params[$i][$j][$l][$m][$n];
1611 # my $orig = $params[$i][$j][0][$l][$m][$n];
1612 my $prep = $params[$i][$j][$k][$l][$m][$n];
1613 if ( $orig != 0 ) {
1614 if ( $percentage ) {
1615 push( @par, ($prep/$orig*100)-100 );
1616 } else {
1617 push( @par, $prep/$orig );
1618 }
1619 } else {
1620 push( @par, $PsN::out_miss_data );
1621 }
1622 }
1623 push( @sub, \@par );
1624 }
1625 push( @prob, \@sub );
1626 }
1627 push( @prep, \@prob );
1628 }
1629 push( @mod, \@prep );
1630 }
1631 push( @relative_estimates, \@mod );
1632 }
1633 }
1634 end relative_estimates
1635
1636 # }}} relative_estimates
1637
1638 # {{{ relative_confidence_limits
1639
1640 start relative_confidence_limits
1641 {
1642 my @params = @{$self -> confidence_limits( class => 'tool::llp',
1643 parameter => $parameter )};
1644 for ( my $i = 0; $i < scalar @params; $i++ ) {
1645 # Loop over models
1646 my @mod = ();
1647 for ( my $j = 1; $j < scalar @{$params[$i]}; $j++ ) {
1648 # Loop over data sets
1649 my %num_lim;
1650 my @nums = sort {$a <=> $b} keys %{$params[$i][$j]};
1651 foreach my $num ( @nums ) {
1652 my @prob_lim = ();
1653 for ( my $n = 0; $n < scalar @{$params[$i][$j]->{$num}}; $n++ ) {
1654 my @side_lim = ();
1655 for ( my $o = 0; $o < scalar @{$params[$i][$j]->{$num}->[$n]}; $o++ ) {
1656 # OBS: the [0] in the $j position points at the first element i.e
1657 # the results of the tool run on the original model
1658 my $orig = $params[$i][0]->{$num}->[$n][$o];
1659 my $prep = $params[$i][$j]->{$num}->[$n][$o];
1660 print "ORIG: $orig, PREP: $prep\n";
1661 if ( $orig != 0 ) {
1662 if ( $percentage ) {
1663 push( @side_lim, ($prep/$orig*100)-100 );
1664 } else {
1665 push( @side_lim, $prep/$orig );
1666 }
1667 } else {
1668 push( @side_lim, $PsN::out_miss_data );
1669 }
1670 }
1671 push( @prob_lim, \@side_lim );
1672 }
1673 $num_lim{$num} = \@prob_lim;
1674 }
1675 push( @mod, \%num_lim );
1676 }
1677 push( @relative_limits, \@mod );
1678 }
1679 }
1680 end relative_confidence_limits
1681
1682 # }}} relative_confidence_limits
1683
1684 # {{{ llp_print_results
1685
1686 start llp_print_results
1687 {
1688 # NOTE! Only valid for models with one problem and one sub problem!
1689
1690 my %relative_values;
1691 $relative_values{'theta_cis'} = $self ->
1692 relative_confidence_limits( parameter => 'theta',
1693 percentage => 1 );
1694 $relative_values{'omega_cis'} = $self ->
1695 relative_confidence_limits( parameter => 'omega',
1696 percentage => 1 );
1697 $relative_values{'sigma_cis'} = $self ->
1698 relative_confidence_limits( parameter => 'sigma',
1699 percentage => 1 );
1700 $relative_values{'thetas'} = $self ->
1701 relative_estimates( parameter => 'theta',
1702 percentage => 1 );
1703 $relative_values{'omegas'} = $self ->
1704 relative_estimates( parameter => 'omega',
1705 percentage => 1 );
1706 $relative_values{'sigmas'} = $self ->
1707 relative_estimates( parameter => 'sigma',
1708 percentage => 1 );
1709 $relative_values{'sethetas'} = $self ->
1710 relative_estimates( parameter => 'setheta',
1711 percentage => 1 );
1712 $relative_values{'seomegas'} = $self ->
1713 relative_estimates( parameter => 'seomega',
1714 percentage => 1 );
1715 $relative_values{'sesigmas'} = $self ->
1716 relative_estimates( parameter => 'sesigma',
1717 percentage => 1 );
1718
1719 my %prep_values;
1720 $prep_values{'theta_cis'} = $self -> confidence_limits( class => 'tool::llp',
1721 parameter => 'theta' );;
1722 $prep_values{'omega_cis'} = $self -> confidence_limits( class => 'tool::llp',
1723 parameter => 'omega' );;
1724 $prep_values{'sigma_cis'} = $self -> confidence_limits( class => 'tool::llp',
1725 parameter => 'sigma' );;
1726 $prep_values{'thetas'} = $self -> thetas;
1727 $prep_values{'omegas'} = $self -> omegas;
1728 $prep_values{'sigmas'} = $self -> sigmas;
1729 $prep_values{'sethetas'} = $self -> sethetas;
1730 $prep_values{'seomegas'} = $self -> seomegas;
1731 $prep_values{'sesigmas'} = $self -> sesigmas;
1732
1733
1734
1735 open( RES, ">".$self -> {'results_file'} );
1736 print RES "Case-Deletion Diagnostic with Log-Likelihood Profiling\n";
1737 # Loop over models
1738 for ( my $i = 0; $i < scalar @{$relative_values{'theta_cis'}}; $i++ ) {
1739 print RES "MODEL:;",$i+1,"\n";
1740 foreach my $param ( 'theta_cis', 'omega_cis', 'sigma_cis' ) {
1741 print RES "\n",uc($param),":\n";
1742 # Loop over data sets
1743
1744 my @nums = sort {$a <=> $b} keys %{$relative_values{$param}[$i][0]};
1745 print RES ";";
1746 foreach my $num ( @nums ) {
1747 printf RES "$num;;;;";
1748 }
1749 print RES "\n";
1750 foreach my $num ( @nums ) {
1751 for ( my $o = 0; $o < scalar @{$relative_values{$param}[$i][0]->{$num}[0]}; $o++ ) {
1752 my $side = $o == 0 ? 'lower' : 'upper';
1753 printf RES ";$side;rel diff (%)";
1754 }
1755 }
1756 print RES "\n";
1757 print RES "orig";
1758 foreach my $num ( @nums ) {
1759 for ( my $o = 0; $o < scalar @{$relative_values{$param}[$i][0]->{$num}[0]}; $o++ ) {
1760 printf RES ";%7.5f",$prep_values{$param}[$i][0]->{$num}[0][$o];
1761 print RES ";0";
1762 }
1763 }
1764 print RES "\n";
1765 for ( my $j = 0; $j < scalar @{$relative_values{$param}[$i]}; $j++ ) {
1766 printf RES "%-7d",$j+1;
1767 my @nums = sort {$a <=> $b} keys %{$relative_values{$param}[$i][$j]};
1768 foreach my $num ( @nums ) {
1769 for ( my $n = 0; $n < scalar @{$relative_values{$param}[$i][$j]->{$num}}; $n++ ) {
1770 for ( my $o = 0; $o < scalar @{$relative_values{$param}[$i][$j]->{$num}[$n]}; $o++ ) {
1771 my $rel = $relative_values{$param}[$i][$j]->{$num}[$n][$o];
1772 my $prep = $prep_values{$param}[$i][$j+1]->{$num}[$n][$o];
1773 printf RES ";%7.5f",$prep;
1774 printf RES ";%3.0f",$rel;
1775 }
1776 }
1777 }
1778 print RES "\n";
1779 }
1780 }
1781 }
1782 # Skipped id's, keys and values:
1783 # Loop over models
1784
1785
1786 # sub process_inner_results {
1787 # my $res_ref = shift;
1788 # my $pad = shift;
1789 # $pad++;
1790 # foreach my $res ( @{$res_ref} ) {
1791 # if ( ref ( $res ) eq 'ARRAY' ) {
1792 # print "$pad ARRAY size ",scalar @{$res},"\n";
1793 # process_inner_results( $res, $pad );
1794 # } elsif ( ref ( $res ) eq 'HASH' ) {
1795 # print "$pad HASH keys ",keys %{$res},"\n";
1796 # } else {
1797 # print "$pad OTHER\n";
1798 # }
1799 # }
1800 # }
1801 # process_inner_results( $self -> {'results'}, 0 );
1802
1803 my $i = 1;
1804 foreach my $own ( @{$self -> {'results'} -> {'own'}} ) {
1805 # print "REF1: ",ref($mod),"\n";
1806 # foreach my $prob ( @{$mod} ) {
1807 # print "REF2: ",ref($prob),"\n";
1808 # foreach my $subprob ( @{$prob} ) {
1809 # print "REF3: ",ref($subprob),"\n";
1810 # print "KEYS: ",keys %{$subprob},"\n";
1811 # }
1812 # }
1813 print RES "MODEL $i\n";
1814 foreach my $param ( 'skipped_ids', 'skipped_keys', 'skipped_values' ) {
1815 print RES uc($param),"\n";
1816 my $j = 1;
1817 foreach my $prep ( @{$own -> {$param}} ) {
1818 print RES "Bin no;$j;";
1819 foreach my $val ( @{$prep} ) {
1820 print RES ";$val";
1821 }
1822 print RES "\n";
1823 $j++;
1824 }
1825 }
1826 $i++;
1827 }
1828
1829
1830 # for ( my $j = 0; $j < scalar @{$relative_values{'thetas_cis'}->[$i]}; $j++ ) {
1831 # print RES "MODEL:;",$j+1,"\n";
1832 # # Loop over problems (sort of, at least)
1833 # for ( my $l = 0; $l < scalar @{$relative_values{'thetas_cis'}->[$i]->[$j]->[0]}; $l++ ) {
1834 # # Loop over sub problems (sort of, at least)
1835 # for ( my $m = 0; $m < scalar @{$relative_values{'thetas_cis'}->[$i]->[$j]->[0]->[$l]}; $m++ ) {
1836 # # foreach my $param ( 'thetas', 'omegas', 'sigmas',
1837 # # 'sethetas', 'seomegas', 'sesigmas' ) {
1838 # foreach my $param ( 'theta_cis' ) {
1839 # print RES uc($param),":\n\n";
1840 # # Here one could add printing of parameter names, i.e. 'CL V...' or 'TH1 TH2...'
1841 # # Loop over data sets
1842 # for ( my $k = 0; $k < scalar @{$relative_values{$param}->[$i]->[$j]}; $k++ ) {
1843 # printf RES "%-7d",$k+1;
1844 # for ( my $n = 0; $n < scalar @{$relative_values{$param}[$i][$j][$k][$l][$m]}; $n++ ) {
1845 # for ( my $o = 0; $o < scalar @{$relative_values{$param}[$i][$j][$k][$l][$m][$n]}; $o++ ) {
1846 # my $rel = $relative_values{$param}->[$i][$j][$k][$l][$m][$n];
1847 # my $prep = $prep_values{$param}->[$j][$k][$l][$m][$n];
1848 # printf RES ";%7.5f",$prep;
1849 # printf RES ";%3.0f",$rel;
1850 # }
1851 # print RES "\n";
1852 # }
1853 # print RES "\n";
1854 # }
1855 # print RES "\n";
1856 # }
1857 # }
1858 # print "\n";
1859 # }
1860 # print "\n\n";
1861 # }
1862 # }
1863 close( RES );
1864 }
1865 end llp_print_results
1866
1867 # }}} llp_print_results
1868
1869 # {{{ general_print_results
1870
1871 start general_print_results
1872 {
1873 unless ( defined $self -> {'results'} ) {
1874 print "WARNING: cdd->general_print_results: no return values defined;\n"
1875 ."cannot print results\n";
1876 return;
1877 }
1878 my %results = %{$self -> {'results'}};
1879
1880 open( RES, ">".$self -> {'results_file'} );
1881 print RES "Case-Deletion Diagnostic\n";
1882
1883 # Print meta data
1884
1885 unless ( defined $results{'own'} ) {
1886 print "WARNING: cdd->general_print_results: no own return values defined;\n"
1887 ."cannot print results\n";
1888 return;
1889 }
1890
1891
1892 my @own_results = @{$results{'own'}};
1893 foreach my $result_unit ( @own_results ) {
1894 print RES $result_unit -> {'name'},"\n";
1895 print RES $result_unit -> {'comment'},"\n";
1896 my @values = defined $result_unit{'values'} ? @{$result_unit{'values'}} : ();
1897 my @labels = defined $result_unit{'labels'} ? @{$result_unit{'labels'}} : ();
1898 # Loop the models
1899 for ( my $i = 0; $i <= $#values; $i++ ) {
1900 # Loop the problems
1901 for ( my $j = 0; $j <= $#values[$i]; $j++ ) {
1902 # Loop the sub problems
1903 for ( my $k = 0; $k <= $#values[$i][$j]; $k++ ) {
1904 # Loop the first result dimension
1905 for (my $l = 0; $l <= $#values[$i][$j][$k]; $l++ ) {
1906 # Loop the second result dimension
1907 for ( my $m = 0; $m <= $#values[$i][$j][$k][$l]; $m++ ) {
1908 # Loop the second result dimension
1909 for ( my $m = 0; $m <= $#values[$i][$j]$k][$l]; $m++ ) {
1910 foreach my $model_res ( @values ) {
1911 foreach my $prob_res ( @{$model_unit} ) {
1912 foreach my $subprob_res ( @{$prob_unit} ) {
1913 foreach my $subprob_res ( @{$prob_unit} ) {
1914
1915 close ( RES );
1916 }
1917 end general_print_results
1918
1919 # }}} general_print_results
1920
1921 # {{{ modelfit_print_results
1922
1923 start modelfit_print_results
1924 {
1925 my @parameters = ( 'theta', 'omega', 'sigma',
1926 'setheta', 'seomega', 'sesigma' );
1927 my %relative_values;
1928 my %prep_values;
1929 my %orig_values;
1930 foreach my $parameter ( @parameters ) {
1931 my $accessor = $parameter.'s';
1932 $relative_values{$parameter} = $self ->
1933 relative_estimates( parameter => $parameter,
1934 percentage => 1 );
1935 $prep_values{$parameter} = $self -> $accessor;
1936 $orig_values{$parameter} = $self -> $accessor( original_models => 1 );
1937 }
1938 # sub process_results {
1939 # my $res_ref = shift;
1940 # my $pad = shift;
1941 # $pad++;
1942 # foreach my $res ( @{$res_ref} ) {
1943 # if ( ref ( $res ) eq 'ARRAY' ) {
1944 # process_results( $res, $pad );
1945 # } else {
1946 # print "final $pad\t$res\n";
1947 # }
1948 # }
1949 # }
1950 # process_results( $relative_values{'thetas'}, 0 );
1951
1952
1953 my %nparam;
1954 print "Calling nthetas\n";
1955 $nparam{'thetas'} = $self -> nthetas( original_models => 1 );
1956 print "Done that\n";
1957 open( RES, ">".$self -> {'results_file'} );
1958 print RES "Case-Deletion Diagnostic\n";
1959 # Date information to be added
1960 # print RES "Date:;;;;",$self -> {'date'},"\n";
1961
1962 print RES "Modelfiles:";
1963 foreach my $model ( @{$self -> {'models'}} ) {
1964 print RES ";;;;",$model -> filename,"\n";
1965 }
1966
1967 # Based on columns and number of datasets might better be shown if split by
1968 # model and problem:
1969 print RES "Based on columns:";
1970 my $vars = $self -> {'case_columns'};
1971 if ( ref( $vars ) eq 'ARRAY' ) {
1972 foreach my $vars2 ( @{$vars} ) {
1973 if ( ref( $vars2 ) eq 'ARRAY' ) {
1974 foreach my $vars3 ( @{$vars2} ) {
1975 print RES ";;;;$vars3\n";
1976 }
1977 } else {
1978 print RES ";;;;$vars2\n";
1979 }
1980 }
1981 } else {
1982 print RES ";;;;$vars\n";
1983 }
1984
1985 print RES "Number of data sets:";
1986 my $bins = $self -> {'bins'};
1987 if ( ref( $bins ) eq 'ARRAY' ) {
1988 foreach my $bins2 ( @{$bins} ) {
1989 if ( ref( $bins2 ) eq 'ARRAY' ) {
1990 foreach my $bins3 ( @{$bins2} ) {
1991 print RES ";;;;$bins3\n";
1992 }
1993 } else {
1994 print RES ";;;;$bins2\n";
1995 }
1996 }
1997 } else {
1998 print RES ";;;;$bins\n";
1999 }
2000
2001 print RES "Selection:;;;;",$self-> {'selection_method'},"\n";
2002 if ( defined $self -> {'seed'} ) {
2003 print RES "Seed number:;;;;",$self -> {'seed'},"\n";
2004 } else {
2005 print RES "No seed number specified\n";
2006 }
2007
2008 # TODO: $skip_keys etc from data->case_deletion must be transferred back to
2009 # the main process and appropriate attributes set.
2010
2011 print RES "\n\n\n\n";
2012
2013 # process_results( $relative_values{'thetas'}->[0]->[0]->[0], 0 );
2014
2015 for ( my $i = 0; $i < scalar @{$relative_values{'theta'}}; $i++ ) {
2016 # Loop over models
2017 for ( my $j = 0; $j < scalar @{$relative_values{'theta'}->[$i]}; $j++ ) {
2018 print RES "MODEL:;",$j+1,"\n";
2019 # Loop over problems (sort of, at least)
2020 for ( my $l = 0; $l < scalar @{$relative_values{'theta'}->[$i]->[$j]->[0]}; $l++ ) {
2021 # Loop over sub problems (sort of, at least)
2022 for ( my $m = 0; $m < scalar @{$relative_values{'theta'}->[$i]->[$j]->[0]->[$l]}; $m++ ) {
2023 foreach my $param ( @parameters ) {
2024 print RES uc($param),":\n\n";
2025 # Here one could add printing of parameter names, i.e. 'CL V...' or 'TH1 TH2...'
2026 # Loop over data sets
2027 print RES ";";
2028 for ( my $n = 1; $n <=scalar @{$relative_values{$param}[$i][$j][0][$l][$m]}; $n++ ) {
2029 printf RES "estimate;rel diff (%);";
2030 }
2031 print RES "\n;";
2032 for ( my $n = 1; $n <= scalar @{$relative_values{$param}[$i][$j][0][$l][$m]}; $n++ ) {
2033 printf RES "$n;;";
2034 }
2035 print RES "\n";
2036 print RES "orig";
2037 for ( my $n = 0; $n < scalar @{$relative_values{$param}[$i][$j][0][$l][$m]}; $n++ ) {
2038 printf RES ";%7.5f",$orig_values{$param}[$j][$l][$m][$n];
2039 print RES ";0";
2040 }
2041 print RES "\n";
2042 for ( my $k = 0; $k < scalar @{$relative_values{$param}->[$i]->[$j]}; $k++ ) {
2043 printf RES "%-7d",$k+1;
2044 for ( my $n = 0; $n < scalar @{$relative_values{$param}[$i][$j][$k][$l][$m]}; $n++ ) {
2045 my $rel = $relative_values{$param}->[$i][$j][$k][$l][$m][$n];
2046 my $prep = $prep_values{$param}->[$j][$k+1][$l][$m][$n];
2047 printf RES ";%7.5f",$prep;
2048 printf RES ";%3.0f",$rel;
2049 }
2050 print RES "\n";
2051 }
2052 print RES "\n";
2053 }
2054 }
2055 print "\n";
2056 }
2057 print "\n\n";
2058 }
2059 }
2060
2061 # sub process_inner_results {
2062 # my $res_ref = shift;
2063 # my $pad = shift;
2064 # $pad++;
2065 # foreach my $res ( @{$res_ref} ) {
2066 # if ( ref ( $res ) eq 'ARRAY' ) {
2067 # print "$pad ARRAY size ",scalar @{$res},"\n";
2068 # process_inner_results( $res, $pad );
2069 # } elsif ( ref ( $res ) eq 'HASH' ) {
2070 # print "$pad HASH keys ",keys %{$res},"\n";
2071 # } else {
2072 # print "$pad OTHER\n";
2073 # }
2074 # }
2075 # }
2076 # process_inner_results( $self -> {'results'}, 0 );
2077 # die;
2078 # Skipped id's, keys and values:
2079 # Loop over models
2080 my $i = 1;
2081 foreach my $own ( @{$self -> {'results'} -> {'own'}} ) {
2082 print RES "MODEL $i\n";
2083 foreach my $param ( 'skipped_ids', 'skipped_keys', 'skipped_values' ) {
2084 print RES uc($param),"\n";
2085 my $j = 1;
2086 foreach my $prep ( @{$own -> {$param}} ) {
2087 print RES "Bin no;$j;";
2088 foreach my $val ( @{$prep} ) {
2089 print RES ";$val";
2090 }
2091 print RES "\n";
2092 $j++;
2093 }
2094 }
2095 $i++;
2096 }
2097
2098 close ( RES );
2099 }
2100 end modelfit_print_results
2101
2102 # }}} modelfit_print_results
2103
2104 # {{{ prepare_results
2105
2106 start prepare_results
2107 {
2108 if ( not defined $self -> {'raw_results'} ) {
2109 $self -> read_raw_results();
2110 }
2111 }
2112 end prepare_results
2113
2114 # }}}
2115
2116 # {{{ print_summary
2117 start print_summary
2118 {
2119 my ($outside_n_sd );
2120 for( my $i = 0; $i < scalar @{$self -> {'raw_results_header'} -> [0]} ; $i++) {
2121 if( $self -> {'raw_results_header'} -> [0][$i] eq 'outside.n.sd' ){
2122 $outside_n_sd = $i;
2123 }
2124 }
2125 for( my $model_i; $model_i <= $#{$self -> {'raw_results'}}; $model_i++ ){
2126 for( my $problem_i; $problem_i <= $#{$self -> {'raw_results'} -> [$model_i]}; $problem_i++ ){
2127 my $test_val = $self -> {'raw_results'} -> [$model_i] -> [$problem_i] -> [$outside_n_sd];
2128 if( $test_val ne '' ){
2129 print "Outside ". $self -> {'outside_n_sd_check'};
2130 if( $test_val == 1 ){
2131 print " \t[ WARNING ]\n";
2132 } else {
2133 print " \t[ OK ]\n";
2134 }
2135 }
2136 }
2137 }
2138 }
2139 end print_summary
2140 # }}} print_summary
2141
2142 # {{{ update_raw_results
2143 start update_raw_results
2144 {
2145 my $cook_scores;
2146 my $cov_ratios;
2147 my $outside_n_sd;
2148
2149 # foreach my $section( @{$self -> {'results'}[0] -> {'own'}} ){
2150 # if( $section -> {'name'} eq 'cook.scores' ){
2151 # $cook_scores = $section -> {'values'};
2152 # }
2153 # if( $section -> {'name'} eq 'cov.ratio' ){
2154 # $cov_ratios = $section -> {'values'};
2155 # }
2156 # if( $section -> {'name'} eq 'outside.n.sd' ){
2157 # $outside_n_sd = $section -> {'values'};
2158 # }
2159 # }
2160
2161 my ($dir,$file) =
2162 OSspecific::absolute_path( $self -> {'directory'},
2163 $self -> {'raw_results_file'}[$model_number-1] );
2164 open( RRES, $dir.$file );
2165 my @rres = <RRES>;
2166 close( RRES );
2167 open( RRES, '>',$dir.$file );
2168
2169 chomp( $rres[0] );
2170 print RRES $rres[0] . ",cook.scores,cov.ratios,outside.n.sd\n";
2171 chomp( $rres[1] );
2172 print RRES $rres[1] . ",0,1,0\n";
2173
2174 my @new_rres;
2175 for( my $i = 2 ; $i <= $#rres; $i ++ ) {
2176 my $row_str = $rres[$i];
2177 chomp( $row_str );
2178 $row_str .= sprintf( ",%.5f,%.5f,%1f\n" ,
2179 $self -> {'cook_scores'} -> [$i-2],
2180 $self -> {'covariance_ratios'} -> [$i-2],
2181 $self -> {'outside_n_sd'} -> [$i-2] );
2182 print RRES $row_str;
2183 }
2184 close( RRES );
2185 }
2186 end update_raw_results
2187
2188 # }}} update_raw_results
2189
2190 # {{{ create_R_scripts
2191 start create_R_scripts
2192 {
2193 unless( -e $PsN::lib_dir . '/R-scripts/cdd.R' ){
2194 'debug' -> die( message => 'CDD R-script are not installed, no matlab scripts will be generated.' );
2195 return;
2196 }
2197 cp ( $PsN::lib_dir . '/R-scripts/cdd.R', $self -> {'directory'} );
2198 }
2199 end create_R_scripts
2200 # }}}
Perl-speaks-NONMEM CVS import