?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Lots of fixes for 2.2, see README.txt for details
[PsN.git] / lib / tool / scm_subs.pm
blobf458ce8fe31e57e4365b23629b278c5d91065af5
1 # TODO, possible improvements.
2 # * Interface for gof-function needs an overhaul
3 # * General readability improvments in setup would be nice, but not
4 # * necessary.
5 #
6 # 2004-06-10:
7 # - The header in the data file is still needed. Make the scm
8 # calculate the correct columns from the $INPUT record.
9 # - Make it possible to supply statistics (min, max, median) of
10 # covariates as arguments to new().
11
12 # 2004-08-26
13 # Verkar som om kovariatnamnen måste vara unika. (Se substr i
14 # _option_val_pos i problem_subs) Annars kan kovariater som används
15 # droppas när en annan kovariat med kortare namn som börjar lika dant
16 # skall droppas (CANT och CA tex)
17
18 # Attributes:
19 #
20 # Name: Type: Default Value:
21 # 'label_type' Hash of strings
22 #
23 # label_type has keys named $parm.$cov, e.g CLWGT or VAGE. Values are
24 # either 'con' for continuous or 'cat' for categorical.
25 #
26 # Name: Type: Default Value:
27 # 'work_queue' Array of step objects
28 #
29 # Name: Type: Default Value:
30 # 'spec_code' Array of step objects
31 #
32
33 # {{{ include
34
35 start include statements
36 use strict;
37 use Cwd;
38 use tool::modelfit;
39 use OSspecific;
40 use Data::Dumper;
41 use File::Copy 'cp';
42 use status_bar;
43 end include statements
44
45 # }}}
46
47 # {{{ new
48
49 start new
50 {
51 # <I>test_relations, fix and p_value</I> can be specified
52 # as either a reference to a hash or as a reference to an
53 # array (for different settings for each model).
54
55 if( defined $this -> {'config_file_name'} or defined $this -> {'config_file'} ){
56 $this -> read_config_file;
57 }
58
59 # This block is a duplicate of the settign of the 'directory' attribute in tool.pm
60 # It needs to be here too to make sure that the correct directory is set when
61 # resuming an scm.
62 if ( defined $parm{'directory'} ) {
63 my $dummy;
64 ( $this -> {'directory'}, $dummy ) = OSspecific::absolute_path( $parm{'directory'}, '');
65 }
66
67 debug -> die( message => "You need to specify \'models'\ either as argument or in the config file.")
68 unless ( defined $this -> {'models'} and scalar @{$this -> {'models'}} > 0 );
69
70 # debug -> die( message => "tool::scm->new: You need to specify \'parameters\' either".
71 # " as argument or in the config file.")
72 # unless ( defined $this -> {'parameters'} and scalar @{$this -> {'parameters'}} > 0 );
73
74 if ( not( ref($this -> {'p_value'}) eq 'ARRAY' or
75 ref($this -> {'p_value'}) eq 'HASH' ) ) {
76 if ( not defined $this -> {'p_value'} or
77 $this -> {'p_value'} eq '' ) {
78 $this -> {'p_value'} = 0.05;
79 }
80 my $p_value = $this -> {'p_value'};
81 $this -> {'p_value'} = [];
82 foreach my $model ( @{$this -> {'models'}} ) {
83 push ( @{$this -> {'p_value'}}, $p_value ) ;
84 }
85 }
86
87 if ( not( ref($this -> {'fix'}) eq 'ARRAY' or
88 ref($this -> {'fix'}) eq 'HASH' ) ) {
89 my $fix = $this -> {'fix'};
90 $this -> {'fix'} = [];
91 foreach my $model ( @{$this -> {'models'}} ) {
92 push ( @{$this -> {'fix'}}, $fix ) ;
93 }
94 }
95
96 foreach my $file ( 'logfile', 'short_logfile', 'raw_results_file',
97 'final_model_name', 'covariate_statistics_file',
98 'relations_file' ) {
99 if ( not( ref($this -> {$file}) eq 'ARRAY' or
100 ref($this -> {$file}) eq 'HASH' ) ) {
101 my ($ldir, $filename) = OSspecific::absolute_path( $this -> {'directory'},
102 $this -> {$file} );
103 $this -> {$file} = [];
104 for ( my $i = 1; $i <= scalar @{$this -> {'models'}}; $i++ ) {
105 my $tmp = $filename;
106 $tmp =~ s/\./$i\./;
107 push ( @{$this -> {$file}}, $ldir.$tmp ) ;
108 }
109 }
110 }
111
112 # if ( not( ref($this -> {'final_model_name'}) eq 'ARRAY' or
113 # ref($this -> {'final_model_name'}) eq 'HASH' ) ) {
114 # my ($ldir, $final_model_name) = OSspecific::absolute_path( $this -> {'directory'},
115 # $this -> {'final_model_name'} );
116 # $this -> {'final_model_name'} = [];
117 # for ( my $i = 1; $i <= scalar @{$this -> {'models'}}; $i++ ) {
118 # my $tmp = $final_model_name;
119 # $tmp =~ s/\./$i\./;
120 # push ( @{$this -> {'final_model_name'}}, $ldir.$tmp ) ;
121 # }
122 # }
123
124 if ( defined $this -> {'test_relations'} ) {
125 if ( ref ( $this -> {'test_relations'} ) eq 'HASH' ) {
126 my $relations = $this-> {'test_relations'};
127 $this -> {'test_relations'} = [];
128 foreach my $model ( @{$this -> {'models'}} ) {
129 push ( @{$this -> {'test_relations'}}, $relations ) ;
130 }
131 }
132 } else {
133 debug -> die( message => "You need to specify \'test_relations'\ either as argument or in the config file." );
134 }
135
136 # Check Errors and init
137 unless ( defined( $this -> {'categorical_covariates'} )
138 or defined( $this -> {'continuous_covariates'} )) {
139 debug -> die( message => "You must specify either " .
140 "categorical and/or continuous covariates either as argument or in the config file" );
141 }
142 unless ( $this -> {'search_direction'} eq 'forward' or
143 $this -> {'search_direction'} eq 'backward' or
144 $this -> {'search_direction'} eq 'both' ) {
145 debug -> die( message => "You must specify the search direction ",
146 "either as \"forward\" or \"backward\". Default is \"forward\"." );
147 }
148
149 if ( ref(\$this -> {'gof'}) eq 'SCALAR' ) {
150 if ( lc($this -> {'gof'}) eq 'ofv' ) {
151 $this -> {'gof'} = \&gof_ofv;
152 } elsif ( lc($this -> {'gof'}) eq 'crc') {
153 $this -> {'gof'} = \&gof_ofv;
154 }
155 }
156
157 # check the validity of the covariates and the relations to be tested
158 for ( my $i = 0; $i < scalar @{$this -> {'models'}}; $i++ ) {
159 my $extra_data_headers = $this -> {'models'}[$i] -> extra_data_headers;
160
161 if ( defined $this -> {'continuous_covariates'} or defined $this -> {'categorical_covariates'} ) {
162
163 my @continuous = defined $this -> {'continuous_covariates'} ? @{$this -> {'continuous_covariates'}} : ();
164 my @categorical = defined $this -> {'categorical_covariates'} ? @{$this -> {'categorical_covariates'}} : ();
165
166 my @not_found = ();
167 foreach my $cov ( @continuous, @categorical ) {
168 my $extra_data_test = 0;
169 if ( defined $extra_data_headers and defined $extra_data_headers -> [0] ) {
170 foreach my $head ( @{$extra_data_headers->[0]} ) {
171 $extra_data_test = 1 and last if ( $cov eq $head );
172 }
173 }
174 push( @not_found, $cov )
175 unless ( $extra_data_test or
176 $this -> {'models'}[$i] -> is_option_set( record => 'input',
177 name => $cov ) );
178 }
179 if ( scalar @not_found and
180 ( not defined $this -> {'models'}[$i] -> extra_files or
181 scalar @{$this -> {'models'}[$i] -> extra_files} == 0 ) ) {
182 debug -> die( message => "Covariate(s) [ ". join(',', @not_found). " ] specified is not defined in " .
183 $this -> {'models'}[$i] -> filename );
184 }
185 }
186
187 if ( defined $this -> {'test_relations'} ) {
188 foreach my $par ( sort keys %{$this -> {'test_relations'}[$i]} ){
189 my @not_found = ();
190 foreach my $cov ( @{$this -> {'test_relations'}[$i]{$par}} ){
191 my @continuous = defined $this -> {'continuous_covariates'} ? @{$this -> {'continuous_covariates'}} : ();
192 my @categorical = defined $this -> {'categorical_covariates'} ? @{$this -> {'categorical_covariates'}} : ();
193
194 my $covariate_test = 0;
195 foreach my $specified_cov ( @continuous, @categorical ) {
196 $covariate_test = 1 and last if( $cov eq $specified_cov );
197 }
198 push( @not_found, $cov ) unless ( $covariate_test );
199 }
200 if ( scalar @not_found and
201 ( not defined $this -> {'models'}[$i] -> extra_files or
202 scalar @{$this -> {'models'}[$i] -> extra_files} == 0 ) ) {
203 debug -> die( message => "Covariate(s) [ " . join( ',', @not_found ). " ] specified for parameter $par " .
204 "in test_relations is not defined as a covariate" );
205 }
206 }
207 }
208 }
209
210 # If no previous information on the statistics of the
211 # covariates is available, initiate this.
212 # First; the continuous covariates:
213 for ( my $i = 0; $i < scalar @{$this -> {'models'}}; $i++ ) {
214 if( -e $this -> {'covariate_statistics_file'}[$i] ) {
215 open( STAT, '<'.$this -> {'covariate_statistics_file'}[$i] );
216 my $tmp;
217 for ( <STAT> ) {
218 $tmp = $tmp.$_;
219 }
220 close( STAT );
221 my $VAR1;
222 eval( $tmp );
223 $this -> {'covariate_statistics'}[$i] = $VAR1;
224 } else {
225
226 my $model = $this -> {'models'} -> [$i];
227 # Assume one $PROBLEM
228
229 my $data = $model -> datas -> [0];
230 my $extra_data;
231 if ( defined $model -> problems -> [0] -> extra_data ) {
232 $extra_data = $model -> problems -> [0] -> extra_data;
233 }
234 unless ( defined $this -> {'covariate_statistics'}[$i] ) {
235 $data -> target('mem');
236 $extra_data -> target('mem') if defined $extra_data;
237 ui -> print( category => 'scm',
238 message => "Calculating continuous covariate statistics" );
239 my $ncov = scalar @{$this -> {'continuous_covariates'}};
240 my $status_bar = status_bar -> new( steps => $ncov );
241 ui -> print( category => 'scm',
242 message => $status_bar -> print_step(),
243 newline => 0);
244
245 for ( my $j = 1; $j <= $ncov; $j++ ) {
246 my $cov = $this -> {'continuous_covariates'} -> [$j-1];
247 # Factors
248 $this -> {'covariate_statistics'}[$i]{$cov}{'factors'} =
249 $model -> factors( problem_number => 1,
250 return_occurences => 1,
251 unique_in_individual => 0,
252 column_head => $cov );
253 # Statistics
254 $this -> {'covariate_statistics'}[$i]{$cov}{'have_missing_data'} =
255 $model -> have_missing_data( column_head => $cov );
256
257 ( $this -> {'covariate_statistics'}[$i]{$cov}{'median'},
258 $this -> {'covariate_statistics'}[$i]{$cov}{'min'},
259 $this -> {'covariate_statistics'}[$i]{$cov}{'max'} ) =
260 $this -> calculate_continuous_statistics( model => $model,
261 covariate => $cov );
262 #my $nl = $j == $ncov ? "" : "\r";
263 if( $status_bar -> tick () ){
264 ui -> print( category => 'scm',
265 message => $status_bar -> print_step(),
266 wrap => 0,
267 newline => 0 );
268 }
269 }
270 ui -> print( category => 'scm',
271 message => " ... done\nCalculating categorical covariate statistics" );
272 $ncov = defined $this -> {'categorical_covariates'} ? scalar @{$this -> {'categorical_covariates'}} : 0;
273
274 my $status_bar = status_bar -> new( steps => $ncov );
275 ui -> print( category => 'scm',
276 message => $status_bar -> print_step(),
277 newline => 0);
278
279 for ( my $j = 1; $j <= $ncov; $j++ ) {
280 my $cov = $this -> {'categorical_covariates'} -> [$j-1];
281 # Factors
282 $this -> {'covariate_statistics'}[$i]{$cov}{'factors'} =
283 $model -> factors( problem_number => 1,
284 return_occurences => 1,
285 unique_in_individual => 0,
286 column_head => $cov );
287 # Fractions
288 $this -> {'covariate_statistics'}[$i]{$cov}{'fractions'} =
289 $model -> fractions( unique_in_individual => 1,
290 column_head => $cov,
291 ignore_missing => 1 );
292 # Statistics
293 $this -> {'covariate_statistics'}[$i]{$cov}{'have_missing_data'} =
294 $model -> have_missing_data( column_head => $cov );
295
296 ( $this -> {'covariate_statistics'}[$i]{$cov}{'median'},
297 $this -> {'covariate_statistics'}[$i]{$cov}{'min'},
298 $this -> {'covariate_statistics'}[$i]{$cov}{'max'} ) =
299 $this -> calculate_categorical_statistics
300 ( factors => $this -> {'covariate_statistics'}[$i]{$cov}{'factors'},
301 have_missing_data => $this -> {'covariate_statistics'}[$i]{$cov}{'have_missing_data'},
302 fix => $this -> {'fix'}[$i] );
303
304 if( $status_bar -> tick () ){
305 ui -> print( category => 'scm',
306 message => $status_bar -> print_step(),
307 wrap => 0,
308 newline => 0 );
309 }
310 }
311 ui -> print( category => 'scm',
312 message => " ... done" );
313 $data -> target('disk');
314 $extra_data -> target('disk') if defined $extra_data;
315 }
316 open( STAT, '>'.$this -> {'covariate_statistics_file'}[$i] );
317 $Data::Dumper::Purity = 1;
318 print STAT Dumper $this -> {'covariate_statistics'}[$i];
319 $Data::Dumper::Purity = 0;
320 close( STAT );
321 }
322 }
323
324 # Default ofv drops at desired p-values (assuming chi-squared
325 # distribution of hirerchical models)
326 my %p_values;
327 #For unlimited stepping: p = 100% /JR
328 $p_values{'1'} = {1=>0,
329 2=>0,
330 3=>0,
331 4=>0,
332 5=>0,
333 6=>0,
334 7=>0,
335 8=>0,
336 9=>0,
337 10=>0};
338 ## p= 0.05
339 $p_values{'0.05'} = {1=>3.84,
340 2=>5.99,
341 3=>7.81,
342 4=>9.49,
343 5=>11.07,
344 6=>12.59,
345 7=>14.07,
346 8=>15.51,
347 9=>16.92,
348 10=>18.31};
349 ## p=0.01
350 $p_values{'0.01'} = {1=>6.63,
351 2=>9.21,
352 3=>11.34,
353 4=>13.28,
354 5=>15.09,
355 6=>16.81,
356 7=>18.48,
357 8=>20.09,
358 9=>21.67,
359 10=>23.21};
360 ## p=0.005
361 $p_values{'0.005'} = {1=>7.88,
362 2=>10.60,
363 3=>12.84,
364 4=>14.86,
365 5=>16.75,
366 6=>18.55,
367 7=>20.28,
368 8=>21.95,
369 9=>23.59,
370 10=>25.19};
371 ## p=0.001
372 $p_values{'0.001'} = {1=>10.83,
373 2=>13.82,
374 3=>16.27,
375 4=>18.47,
376 5=>20.52,
377 6=>22.46,
378 7=>24.32,
379 8=>26.12,
380 9=>27.88,
381 10=>29.59};
382
383
384 # If no previous information on the relations is available,
385 # create the basic parameter-covariate relation data structure
386 # including the information about states (1=not included,
387 # 2=linear relation, 3=hockey-stick relation).
388 # unless ( defined $this -> {'relations'} ) {
389 for ( my $i = 0; $i < scalar @{$this -> {'models'}}; $i++ ) {
390 my $first = 1;
391 if( -e $this -> {'relations_file'}[$i] ) {
392 open( RELATIONS, '<'.$this -> {'relations_file'}[$i] );
393 my $tmp;
394 for ( <RELATIONS> ) {
395 $tmp = $tmp.$_;
396 }
397 close( RELATIONS );
398 my $VAR1;
399 eval( $tmp );
400 $this -> {'relations'}[$i] = $VAR1;
401 foreach my $par ( sort keys %{$this -> {'test_relations'}[$i]} ) {
402 foreach my $cov ( @{$this -> {'test_relations'}[$i]{$par}} ){
403 # Is this covariate continuous or not?
404 my $continuous = 1;
405 foreach my $cat ( @{$this -> {'categorical_covariates'}} ) {
406 $continuous = 0 if ( $cov eq $cat );
407 }
408 my @valid_states;
409 if ( $continuous ) {
410 @valid_states = @{$this -> {'valid_states'}{'continuous'}};
411 } else {
412 @valid_states = @{$this -> {'valid_states'}{'categorical'}};
413 }
414 $this -> {'relations'}[$i]{$par}{$cov}{'state'} = 1;
415 foreach my $state ( @valid_states ) {
416 if ( defined $this -> {'included_relations'}[$i]{$par}{$cov} and
417 $this -> {'included_relations'}[$i]{$par}{$cov}{'state'} == $state ) {
418 $this -> {'included_relations'}[$i]{$par}{$cov}{'code'} =
419 $this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state};
420 $this -> {'included_relations'}[$i]{$par}{$cov}{'nthetas'} =
421 $this -> {'relations'}[$i]{$par}{$cov}{'nthetas'}{$state};
422 $this -> {'included_relations'}[$i]{$par}{$cov}{'inits'} =
423 $this -> {'relations'}[$i]{$par}{$cov}{'inits'}{$state};
424 $this -> {'included_relations'}[$i]{$par}{$cov}{'bounds'} =
425 $this -> {'relations'}[$i]{$par}{$cov}{'bounds'}{$state};
426 }
427 }
428 }
429 }
430 } else {
431 foreach my $par ( sort keys %{$this -> {'test_relations'}[$i]} ) {
432 foreach my $cov ( @{$this -> {'test_relations'}[$i]{$par}} ){
433 # Here the ofv-drops should be defined
434 # I've started 2004-11-09
435 my $ofv_changes = $p_values{$this -> {'p_value'}[$i]};
436 if ( defined $this -> {'ofv_change'} ) {
437 if ( ref($this -> {'ofv_change'}) eq 'HASH' ) {
438 # If only one ofv_drop given for all models and all relations
439 while ( my ( $df, $ofv ) = each %{$this -> {'ofv_change'}} ) {
440 $ofv_changes -> {$df} = $ofv;
441 }
442 if ( $first ) {
443 open( LOG, ">>".$this -> {'logfile'} -> [$i] );
444 print LOG "Using user-defined ofv change criteria\n";
445 print LOG "Degree of freedom | Required ofv change\n";
446 my @dfs = sort {$a <=> $b} keys %{$ofv_changes};
447 foreach my $df ( @dfs ) {
448 print LOG " $df - ",
449 $ofv_changes -> {$df},"\n";
450 }
451 close( LOG );
452 }
453 } elsif ( ref($this -> {'ofv_change'}) eq 'ARRAY' and
454 scalar @{$this -> {'ofv_change'}} > 0 ) {
455 # If ofv_drops given for all relations
456 if ( defined $this -> {'ofv_change'} -> [$i] and
457 ref( $this -> {'ofv_change'} -> [$i]) eq 'HASH' and
458 defined $this -> {'ofv_change'} -> [$i]{$par} and
459 ref( $this -> {'ofv_change'} -> [$i]{$par}) eq 'ARRAY' and
460 defined $this -> {'ofv_change'} -> [$i]{$par}{$cov} and
461 ref( $this -> {'ofv_change'} -> [$i]{$par}{$cov}) eq 'HASH' ) {
462 while ( my ( $df, $ofv ) = each
463 %{$this -> {'ofv_change'} -> [$i]{$par}{$cov}} ) {
464 $ofv_changes -> {$df} = $ofv;
465 }
466 }
467 } else {
468 # Attempt to give ofv drops for all relations, but wrong structure
469 debug -> warn( level => 2,
470 message => "Incorrect structure used for ofv_change, using default" );
471 }
472 }
473 $this -> {'relations'}[$i]{$par}{$cov}{'ofv_changes'} = $ofv_changes;
474 # Is this covariate continuous or not?
475 my $continuous = 1;
476 foreach my $cat ( @{$this -> {'categorical_covariates'}} ) {
477 $continuous = 0 if ( $cov eq $cat );
478 }
479 $this -> {'relations'}[$i]{$par}{$cov}{'continuous'} = $continuous;
480 my @valid_states;
481 if ( $continuous ) {
482 @valid_states = @{$this -> {'valid_states'}{'continuous'}};
483 } else {
484 @valid_states = @{$this -> {'valid_states'}{'categorical'}};
485 }
486 $this -> {'relations'}[$i]{$par}{$cov}{'state'} = 1;
487 foreach my $state ( @valid_states ) {
488 if ( defined $this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state} ) {
489 if ( not ref $this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state} eq 'ARRAY' ) {
490 debug -> die( message => "The code specified for $par $cov $state is not ".
491 "an array\n" );
492 } else {
493 for ( @{$this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state}} ) {
494 s/PARCOV/$par$cov/g;
495 s/PAR/$par/g;
496 s/COV/$cov/g;
497 }
498 }
499 } else {
500 $this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state} = [];
501 }
502 $this -> {'relations'}[$i]{$par}{$cov}{'inits'}{$state} = [] unless
503 ( defined $this -> {'relations'}[$i]{$par}{$cov}{'inits'}{$state} );
504 $this -> {'relations'}[$i]{$par}{$cov}{'bounds'}{$state} = {} unless
505 ( defined $this -> {'relations'}[$i]{$par}{$cov}{'bounds'}{$state} );
506 ( $this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state},
507 $this -> {'relations'}[$i]{$par}{$cov}{'nthetas'}{$state},
508 $this -> {'relations'}[$i]{$par}{$cov}{'inits'}{$state},
509 $this -> {'relations'}[$i]{$par}{$cov}{'bounds'}{$state} ) =
510 $this -> create_code( start_theta => 1,
511 parameter => $par,
512 covariate => $cov,
513 continuous => $continuous,
514 state => $state,
515 code => $this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state},
516 inits => $this -> {'relations'}[$i]{$par}{$cov}{'inits'}{$state},
517 bounds => $this -> {'relations'}[$i]{$par}{$cov}{'bounds'}{$state},
518 statistics => $this -> {'covariate_statistics'}[$i]{$cov},
519 missing_data_token => $this -> {'missing_data_token'},
520 fix => $this -> {'fix'}[$i] );
521 if ( defined $this -> {'included_relations'}[$i]{$par}{$cov} and
522 $this -> {'included_relations'}[$i]{$par}{$cov}{'state'} == $state ) {
523 $this -> {'included_relations'}[$i]{$par}{$cov}{'code'} =
524 $this -> {'relations'}[$i]{$par}{$cov}{'code'}{$state};
525 $this -> {'included_relations'}[$i]{$par}{$cov}{'nthetas'} =
526 $this -> {'relations'}[$i]{$par}{$cov}{'nthetas'}{$state};
527 $this -> {'included_relations'}[$i]{$par}{$cov}{'inits'} =
528 $this -> {'relations'}[$i]{$par}{$cov}{'inits'}{$state};
529 $this -> {'included_relations'}[$i]{$par}{$cov}{'bounds'} =
530 $this -> {'relations'}[$i]{$par}{$cov}{'bounds'}{$state};
531 }
532 }
533 $first = 0;
534 }
535 }
536 open( RELATIONS, '>'.$this -> {'relations_file'}[$i] );
537 $Data::Dumper::Purity = 1;
538 print RELATIONS Dumper $this -> {'relations'}[$i];
539 $Data::Dumper::Purity = 0;
540 close( RELATIONS );
541 }
542 }
543 # }
544 }
545 end new
546
547 # }}}
548
549 # {{{ register_in_database
550 start register_in_database
551 {
552 if ( $PsN::config -> {'_'} -> {'use_database'} ) {
553 tool::register_in_database( $self );
554 my $dbh = DBI -> connect("DBI:mysql:host=".$PsN::config -> {'_'} -> {'database_server'}.
555 ";databse=".$PsN::config -> {'_'} -> {'project'},
556 $PsN::config -> {'_'} -> {'user'},
557 $PsN::config -> {'_'} -> {'password'},
558 {'RaiseError' => 1});
559 my $sth;
560 $sth = $dbh -> prepare("INSERT INTO ".$PsN::config -> {'_'} -> {'project'}.
561 ".scm (tool_id,direction,step_number) ".
562 "VALUES (".$self -> {'tool_id'}.", '".
563 $self -> {'search_direction'}."', '".
564 $self -> {'step_number'}."' )");
565 $sth -> execute;
566 $self -> {'scm_id'} = $sth->{'mysql_insertid'};
567 $sth -> finish;
568 $dbh -> disconnect;
569 }
570 }
571 end register_in_database
572 # }}} register_in_database
573
574 # {{{ read_config_file
575
576 start read_config_file
577 {
578 unless( defined $self -> {'config_file'} ){
579 $self -> {'config_file'} = config_file -> new( file -> new( path => './', name => $self -> {'config_file_name'} ) );
580 }
581
582 my $config_file = $self -> {'config_file'};
583
584 if( defined( $config_file -> relations ) ){
585 $self -> {'relations'} = [$config_file -> relations];
586 }
587
588 foreach my $config_option ( keys %{$config_file -> valid_scalar_options},
589 keys %{$config_file -> valid_array_options},
590 keys %{$config_file -> valid_hash_options},
591 keys %{$config_file -> valid_code_options} ){
592
593 # These are options passed to the modelfile in bin/scm.pl
594 next if( $config_option eq 'extra_data_files' );
595 next if( $config_option eq 'extra_files' );
596
597 # Handle some special cases (where the option is per model)
598 if( $config_option eq 'base_criteria_values' and defined $config_file -> base_criteria_values ){
599 @{$self -> {'base_criteria_values'}}[0] = $config_file -> base_criteria_values;
600 } elsif( $config_option eq 'test_relations' and defined $config_file -> test_relations ){
601 @{$self -> {'test_relations'}}[0] = $config_file -> test_relations;
602 } elsif( $config_option eq 'included_relations' and defined $config_file -> included_relations ){
603 @{$self -> {'included_relations'}}[0] = $config_file -> included_relations;
604 } elsif( defined $config_file -> $config_option ){
605
606 # This is the general case where we just copy the option.
607 $self -> { $config_option } = $config_file -> $config_option;
608 }
609 }
610 }
611 end read_config_file
612
613 # }}}
614
615 # {{{ _read_scm_file
616
617 start _read_scm_file
618 {
619 # Do we need to set an absolute bath of this file somewhere?
620 open(OLDFILE, "<" . $self -> {'scm_file'}) ||
621 die "Unable to open oldstyle file".$self -> {'scm_file'}."for parsing\n";
622 my $have_pval = 0;
623 while ( <OLDFILE> ) {
624 chomp;
625 if ( /^\s*;;;\s*DECL:\s*(.*)\s*$/ ) { # Find DECL lines, $1 is string of covs
626 foreach (split(/\s+/, $1)) {
627 $self -> {'labels'}{$_} = '';
628 }
629 # ...else Find BOUNDS lines, $1 is bounds, $2 is cov
630 } elsif ( /^\s*;;;\s*BOUNDS:\s*(.*)\s*;\s*(.*)\s*$/ ) {
631 my $label = $2;
632 $1 =~ /\s*\((.*),(.*),(.*)\)/;
633 $self -> {'bounds'}{$label} = [$1,$3];
634 if ( length($2) > 0 ) {
635 $self -> {'inits'}{$label} = $2;
636 }
637 } elsif (/^\s*;;;\s*GLOBAL_INIT:\s*(.*)\s*$/) {
638 if ( defined $self -> {'global_init'} ) {
639 debug -> warn( level => 1,
640 message => "Global init already specified, ".
641 "either as parameter to scm -> new, or multiple".
642 " times in scm file, using last" );
643 }
644 $self -> {'global_init'} = $1;
645 } elsif ( /^\s*;;;\s*DYNAMIC_INIT\s*$/) {
646 $self -> {'dynamic_init'} = 1;
647 # ... else Find MISS lines, $1 is missing data token
648 } elsif ( /^\s*;;;\s*MISS:\s*(.*)\s*/) {
649 debug -> warn( level => 2,
650 message => "Found MISS: $1" );
651 } elsif ( /^\s*;;;\s*GOF:\s*(.*)\s*/) {
652 if ( not lc($1) eq 'extern' ) {
653 $self -> {'gof'} = $1;
654 }
655 if ( $have_pval and not lc($self -> {'gof'}) eq 'ofv' ) {
656 debug -> warn( level => 1,
657 message => "SCM Warning: PFORW and/or PBACK given in modelfile, but gof-function is not specified to OFV" );
658 }
659 } elsif ( /^\s*;;;\s*TASK:\s*(.*)\s*/) {
660 $self -> {'search_direction'} = $1;
661 } elsif ( /^\s*;;;\s*FIX:\s*(.*)\s*/) {
662 # Must be updated for array format of fix
663 $self -> {'fix'} = $1;
664 } elsif ( /^\s*;;;\s*(\w+)(\d+)\s*START/ ) {
665 my $label = $1;
666 my $level = $2;
667 my @code_array;
668 while (<OLDFILE>) {
669 last if /\s*;;;\s*.*\s*END\s*/;
670 push( @code_array, $_ );
671 }
672 $self -> {'spec_code'}{$label}{$level} = \@code_array;
673 } elsif ( /^\s*;;;\s*PFORW:\s*(\d+)/ ) {
674 if ( defined $self -> {'gof_data'} ) {
675 if ( $have_pval ) {
676 $self -> {'gof_data'} -> {'pforw'} = $1;
677 } else {
678 debug -> warn( level => 1,
679 message => "SCM Warning: both PFORW and \"gof_data\" ".
680 "is specified, PFORW value will be ignored\n" );
681 }
682 } else {
683 $self -> {'gof_data'} -> {'pforw'} = $1;
684 $have_pval = 1;
685 }
686 } elsif ( /^\s*;;;\s*PBACK:\s*(\d+)/ ) {
687 if ( defined $self -> {'gof_data'} ) {
688 if ( $have_pval ) {
689 $self -> {'gof_data'} -> {'pback'} = $1;
690 } else {
691 debug -> warn( level => 1,
692 message => "SCM Warning: both PBACK and \"gof_data\" ".
693 "is specified, PBACK value will be ignored\n" );
694 }
695 } else {
696 $self -> {'gof_data'} -> {'pback'} = $1;
697 $have_pval = 1;
698 }
699 } elsif ( /^\s*;;;\s*LST:\s*(.*)\s*$/ ) {
700 $self -> {'listfile'} = $1;
701 }
702 }
703 if ( defined $self -> {'models'} and
704 scalar @{$self -> {'models'}} > 0 ) {
705 die "Error in scm -> _read_scm_file: You can't ",
706 "specify both a modelobject and a scm modelfile\n";
707 } else {
708 my $scm_model = model -> new( 'filename' => $self -> {'scm_file'},
709 'outputfile' => $self -> {'listfile'},
710 'target' => 'disk',
711 'ignore_missing_files' => 1 );
712
713 my @code_block = @{$scm_model -> pk}; my @new_code_block;
714 my $prev_code_row = '';
715 my $in_block = 0;
716 foreach my $code_row ( @code_block ) {
717 if ( $code_row =~ /^\s*;;;\s*(\w+)(\d+)\s*START/ ) {
718 $in_block = 1;
719 }
720 if ( $in_block and $code_row =~ /\s*;;;\s*.*\s*END\s*/ ) {
721 $in_block = 0;
722 }
723 unless( ($prev_code_row =~ /^$/ and $code_row =~ /^$/) ){
724 unless( $code_row =~ /^\s*;;;/ or $in_block ){
725 push(@new_code_block, $code_row);
726 $prev_code_row = $code_row;
727 }
728 }
729 }
730 $scm_model -> pk( 'new_pk' => \@new_code_block);
731 push( @{$self -> {'models'}}, $scm_model);
732 }
733 }
734 end _read_scm_file
735
736 # }}} _read_scm_file
737
738 # {{{ _raw_results_callback
739
740 start _raw_results_callback
741 {
742 # The goal is to transfer the default modelfit raw_results format of
743 # diagnostic_values-all_thetas-omegas-sigmas
744 # to the format:
745 # diagnostic_values-orig_thetas-omegas-sigmas-scm_thetas
746 # where the scm_thetas are formatted to only hold values where there are active
747 # relations in this step. Do a "print Dumper $self -> {'raw_results'}" before and
748 # after to see the transformation
749
750 my $orig_mod = $self -> {'models'}[$model_number-1];
751 my ( %param_names, %npar_orig );
752 my @params = ( 'theta', 'omega', 'sigma' );
753 my $cols_orig = 0;
754 foreach my $param ( @params ) {
755 my $labels = $orig_mod -> labels( parameter_type => $param );
756 if ( defined $labels ) {
757 $param_names{$param} = $labels -> [0];
758 $npar_orig{$param} = scalar @{$param_names{$param}};
759 $cols_orig += $npar_orig{$param};
760 }
761 }
762
763 my ( @rel_header, @rel_flag, @step_rel_names, %npars );
764
765 # In this loop we create a mesh of all (allowed and) possible parameter-covariate
766 # relations. The active relations of each model [$i] is stored in $rel_flag[$i] as
767 # a 1. All inactive relations are indicated by a 0. A header for the raw_results
768 # file is stored in @rel_header. %npars is a bit superfluous since this
769 # information may be reach through relations-{par}{cov}{'nthetas'}{state} later on.
770 foreach my $parameter ( sort keys %{$self -> {'relations'}[$model_number-1]} ) {
771 foreach my $covariate ( sort keys %{$self -> {'relations'}[$model_number-1]{$parameter}} ) {
772 my ( $in_step, $in_step_state ) = ( 0, undef );
773 my $type = $self -> {'relations'}[$model_number-1]{$parameter}{$covariate}{'continuous'} == 1 ?
774 'continuous' : 'categorical';
775 for ( my $j = 0; $j < scalar @{$self -> {'valid_states'}{$type}}; $j++ ) {
776 my $state = $self -> {'valid_states'}{$type}[$j];
777 my $npar =
778 $self -> {'relations'}[$model_number-1]{$parameter}{$covariate}{'nthetas'}{$state};
779 if ( defined $npar ) { # Skip states without parameters
780 $npars{$parameter}{$covariate}{$state} = $npar;
781 for ( my $k = 1; $k <= $npar; $k++ ) {
782 push( @rel_header, $parameter.$covariate.'-'.$state.'-'.$k );
783 for ( my $i = 0; $i < scalar @{$self -> {'step_relations'}}; $i++ ) {
784 if ( $parameter eq $self -> {'step_relations'} -> [$i]{'parameter'} and
785 $covariate eq $self -> {'step_relations'} -> [$i]{'covariate'} and
786 $state eq $self -> {'step_relations'} -> [$i]{'state'} or
787 ( defined $self -> {'included_relations'}[$model_number-1]{$parameter} and
788 defined $self -> {'included_relations'}[$model_number-1]{$parameter}{$covariate} and
789 $state eq $self -> {'included_relations'}[$model_number-1]{$parameter}{$covariate}{'state'} and not
790 ( $parameter eq $self -> {'step_relations'} -> [$i]{'parameter'} and
791 $covariate eq $self -> {'step_relations'} -> [$i]{'covariate'} ) ) ) {
792 push( @{$rel_flag[$i]}, 1);
793 } else {
794 push( @{$rel_flag[$i]}, 0);
795 }
796 }
797 }
798 }
799 }
800 }
801 }
802
803 # Add up all relations and save the maximum number of thetas used. This number is
804 # used to count backwards from the the end of each raw_results row produced by the
805 # model fit.
806 # IMPORTANT: This "counting backwards" from the end of the raw_results of each
807 # model fit requires that the raw_results from the modelfit object is allways
808 # well-formatted, no matter how the results from each individual model fit are.
809 my $nmax = 0;
810 for ( my $i = 0; $i <= $#rel_flag; $i++ ) {
811 my $sum = 0;
812 for ( my $j = 0; $j < scalar @{$rel_flag[$i]}; $j++ ) {
813 $sum += $rel_flag[$i][$j];
814 }
815 $nmax = $nmax > $sum ? $nmax : $sum;
816 }
817
818 # Use the scm's raw_results file.
819 my ($dir,$file) =
820 OSspecific::absolute_path( $self -> {'directory'},
821 $self -> {'raw_results_file'}[$model_number-1] );
822 my $step_number = $self -> {'step_number'};
823 for ( my $i = 0; $i < scalar @{$self -> {'step_relations'}}; $i++ ) {
824 push( @step_rel_names, $self -> {'step_relations'} -> [$i]{'parameter'}.
825 $self -> {'step_relations'} -> [$i]{'covariate'}.'-'.
826 $self -> {'step_relations'} -> [$i]{'state'} );
827 }
828 my %included_relations = %{$self -> {'included_relations'}[$model_number-1]};
829 my %relations = %{$self -> {'relations'}[$model_number-1]};
830 my @step_relations = @{$self -> {'step_relations'}};
831 my %valid_states = %{$self -> {'valid_states'}};
832 my $action = $self -> {'search_direction'} eq 'forward' ? 'added' : 'removed';
833 $subroutine = sub {
834 my $modelfit = shift;
835 my $mh_ref = shift;
836 my %max_hash = %{$mh_ref};
837 my @new_header = ('step.number', $action.'.relation' );
838 $modelfit -> raw_results_file( $dir.$file );
839 $modelfit -> raw_results_append( 1 ) if ( $step_number > 1 );
840
841 my $raw_results_header = $modelfit -> raw_results_header;
842 my $raw_results = $modelfit -> raw_results;
843 my $cols = scalar @{$modelfit -> {'raw_results'} -> [0]}; # first non-header row
844 for ( my $i = 0; $i < scalar @{$modelfit -> {'raw_results'}}; $i++ ) {
845 my @new_raw_results = ( $step_number, $step_rel_names[$i] );
846
847 my ( @diagnostics, @thetas, @omsi, @sethetas, @seomsi );
848
849
850 # {{{ Get diagnostic results:
851
852 for ( my $j = 0; $j < 16; $j++ ) {
853 push( @diagnostics, $modelfit -> {'raw_results'} -> [$i][$j] );
854 }
855
856 # }}}
857
858 # {{{ Get the thetas that were present in the original model
859
860 for ( my $j = 16;
861 $j < (16+$npar_orig{'theta'}); $j++ ) {
862 push( @thetas, $modelfit -> {'raw_results'} -> [$i][$j] );
863 }
864
865 # }}}
866
867 # {{{ Get the results for all par-cov-relation
868
869 # Initiate $j as starting position for the relation thetas
870 my %res;
871
872 foreach my $kind ( 'estimate', 'se' ) {
873 my $j = $kind eq 'estimate' ?
874 0 : ($max_hash{'theta'} + $npar_orig{'omega'} + $npar_orig{'sigma'});
875 $j += (16+$npar_orig{'theta'});
876 # Important to loop over the sorted hash
877 # Add all included relations estimates
878 foreach my $incl_par ( sort keys %included_relations ) {
879 foreach my $incl_cov ( sort keys %{$included_relations{$incl_par}} ) {
880 next if ( $incl_par eq $step_relations[$i]->{'parameter'} and
881 $incl_cov eq $step_relations[$i]->{'covariate'} );
882 my $npar = $included_relations{$incl_par}{$incl_cov}{'nthetas'};
883 my $state = $included_relations{$incl_par}{$incl_cov}{'state'};
884 for ( my $l = 1; $l <= $npar; $l++ ) {
885 push( @{$res{$incl_par}{$incl_cov}{$state}{$kind}},
886 $modelfit -> {'raw_results'} -> [$i][$j++] );
887 }
888 }
889 }
890
891 # Add the estimates of the relation unique to the model [$i]
892 for ( my $l = 1; $l <= $npars{$step_relations[$i]->{'parameter'}}
893 {$step_relations[$i]->{'covariate'}}
894 {$step_relations[$i]->{'state'}}; $l++ ) {
895 push( @{$res{$step_relations[$i]->{'parameter'}}
896 {$step_relations[$i]->{'covariate'}}
897 {$step_relations[$i]->{'state'}}{$kind}},
898 $modelfit -> {'raw_results'} -> [$i][$j++] );
899 }
900
901 # Sort the results accorrding to the order they appear in (a sorted) $self ->
902 # {'relations'}
903 foreach my $par ( sort keys %npars ) {
904 foreach my $cov ( sort keys %{$npars{$par}} ) {
905 my $type = $relations{$par}{$cov}{'continuous'} == 1 ?
906 'continuous' : 'categorical';
907 foreach my $state ( @{$valid_states{$type}} ) {
908 my $val = ( defined $res{$par} and defined $res{$par}{$cov} and
909 defined $res{$par}{$cov}{$state} and
910 defined $res{$par}{$cov}{$state}{$kind} ) ?
911 $res{$par}{$cov}{$state}{$kind} :
912 [(undef) x $npars{$par}{$cov}{$state}];
913 push( @thetas, @{$val} ) if( defined $val and $kind eq 'estimate' );
914 push( @sethetas, @{$val} ) if( defined $val and $kind eq 'se' );
915 }
916 }
917 }
918 }
919
920 # }}}
921
922 # {{{ Get all the omegas and sigmas
923
924 for ( my $j = (16+$max_hash{'theta'}); $j <
925 ( 16+$max_hash{'theta'} + $npar_orig{'theta'} +
926 $npar_orig{'omega'} + $npar_orig{'sigma'} ); $j++ ) {
927 push( @omsi, $modelfit -> {'raw_results'} -> [$i][$j] );
928 }
929
930 # }}}
931
932 # {{{ Get all the standard errors of the omegas and sigmas
933
934 for ( my $j = (16+($max_hash{'theta'}*2 + $npar_orig{'omega'} + $npar_orig{'sigma'}));
935 $j < (16+($max_hash{'theta'} + $npar_orig{'omega'} + $npar_orig{'sigma'})*2); $j++ ) {
936 push( @seomsi, $modelfit -> {'raw_results'} -> [$i][$j] );
937 }
938
939 # }}}
940
941 push( @new_raw_results, ( @diagnostics, @thetas, @omsi,
942 @sethetas, @seomsi ) );
943
944 $modelfit -> {'raw_results'} -> [$i] = \@new_raw_results;
945 }
946 if ( $step_number == 1 ) {
947 # Loop through the unchanged header and use the header names as accessors
948 # for the original model raw results.
949 my @orig_res = ( 0, undef, 1, 1, 1 );
950 foreach my $param ( @{$raw_results_header} ){
951 next if( $param eq 'model' or $param eq 'problem' or $param eq 'subproblem' );
952 my ( $accessor, $res );
953 if ( $param eq 'theta' or $param eq 'omega' or $param eq 'sigma' or
954 $param eq 'setheta' or $param eq 'seomega' or $param eq 'sesigma' or
955 $param eq 'npomega' or $param eq 'eigen' ) {
956 $accessor = $param.'s';
957 $res = $orig_mod -> {'outputs'} -> [0] -> $accessor;
958
959 } elsif ( $param eq 'shrinkage_etas' ) {
960
961 # Shrinkage does not work for subproblems right now.
962 $res = $orig_mod -> eta_shrinkage;
963
964 } elsif ( $param eq 'shrinkage_wres' ) {
965
966 $res = $orig_mod -> wres_shrinkage;
967
968 } else {
969
970 $res = $orig_mod -> {'outputs'} -> [0] -> $param;
971
972 }
973 # To handle settings on problem level.
974
975 if( defined $res and ref $res -> [0] eq 'ARRAY' ){
976 $res = $res -> [0][0];
977 } else {
978 $res = $res -> [0];
979 }
980
981 if( $max_hash{$param} > 1 or ref $res eq 'ARRAY' ) {
982 if( defined $res ) {
983 push( @orig_res, @{$res} );
984 push( @orig_res, (undef) x ($max_hash{$param} - scalar @{$res}) )
985 unless ( $param eq 'theta' or $param eq 'setheta' );
986 } else {
987 push( @orig_res, (undef) x $max_hash{$param} )
988 unless ( $param eq 'theta' or $param eq 'setheta' );
989 }
990 } else {
991 push( @orig_res, $res );
992 }
993
994 if( $param eq 'theta' or $param eq 'setheta' ) {
995 # Push undef's for all possible relations
996 foreach my $par ( sort keys %npars ) {
997 foreach my $cov ( sort keys %{$npars{$par}} ) {
998 my $type = $relations{$par}{$cov}{'continuous'} == 1 ?
999 'continuous' : 'categorical';
1000 foreach my $state ( @{$valid_states{$type}} ) {
1001 my $val = [(undef) x $npars{$par}{$cov}{$state}];
1002 push( @orig_res, @{$val} );
1003 }
1004 }
1005 }
1006 }
1007
1008 }
1009 unshift( @{$raw_results}, \@orig_res );
1010
1011 my @new_header = ('step.number', $action.'.relation' );
1012 foreach my $name ( @{$raw_results_header} ) {
1013 my @new_names = ();
1014 foreach my $param ( @params ) {
1015 if ( $name eq $param ) {
1016 @new_names = @{$param_names{$param}};
1017 if( $param eq 'theta' ) {
1018 push( @new_names, @rel_header );
1019 }
1020 last;
1021 }
1022 if ( $name eq 'se'.$param ) {
1023 if( $param eq 'theta' ) {
1024 foreach my $head_str ( @{$param_names{$param}}, @rel_header ) {
1025 push( @new_names, 'se'.$head_str );
1026 }
1027 }
1028 last;
1029 }
1030 }
1031 if ( $#new_names >= 0 ) {
1032 push( @new_header, @new_names );
1033 } else {
1034 push( @new_header, $name );
1035 }
1036 }
1037 $modelfit -> {'raw_results_header'} = \@new_header;
1038 }
1039
1040 # $modelfit -> {'raw_results_header'} = [@new_header, @rel_header];
1041 };
1042 return $subroutine;
1043 }
1044 end _raw_results_callback
1045
1046 # }}} _raw_results_callback
1047
1048 # {{{ modelfit_setup
1049
1050 start modelfit_setup
1051 {
1052 my $model = $self -> {'models'} -> [$model_number-1];
1053 # If the number of threads are given per tool, e.g. [2,5] meaning 2 threads for
1054 # scm and 5 for the modelfit.
1055 my $mfit_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
1056 $self -> {'threads'} -> [1]:$self -> {'threads'};
1057 my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
1058 $self -> {'threads'} -> [0]:$self -> {'threads'};
1059 # More threads than models?
1060 my $num = scalar @{$self -> {'models'}};
1061 $own_threads = $num if ( $own_threads > $num );
1062
1063 # Check which models that hasn't been run and run them
1064 # This will be performed each step but will only result in running
1065 # models at the first step, if at all.
1066
1067 # If more than one process is used, there is a VERY high risk of interaction
1068 # between the processes when creating directories for model fits. Therefore
1069 # the {'directory'} attribute is given explicitly below.
1070
1071 unless ( $model -> is_run ) {
1072 my $orig_fit = tool::modelfit -> new
1073 ( reference_object => $self,
1074 prepared_models => undef,
1075 subtools => undef,
1076 results => undef,
1077 logfile => undef,
1078 raw_results => undef,
1079 base_directory => $self -> {'directory'},
1080 directory => $self -> {'directory'}.'/orig_modelfit_dir'.$model_number.'/',
1081 models => [$model],
1082 subtools => undef,
1083 top_tool => 0,
1084 parent_tool_id => $self -> {'tool_id'},
1085 threads => $mfit_threads,
1086 parent_threads => $own_threads );
1087 ui -> print( category => 'scm',
1088 message => "Evaluating basic model" ) unless ( $self -> {'parent_threads'} > 1 );
1089 $orig_fit -> run;
1090 if( $self -> {'scm_file'} and not -e $self -> {'scm_file'} ){
1091 $self -> _write( model_number => $model_number,
1092 filename => $self -> {'scm_file'} );
1093 }
1094 }
1095 ( $self -> {'prepared_models'}[$model_number-1]{'own'}, $self -> {'step_relations'} ) =
1096 $self -> _create_models( model_number => $model_number,
1097 orig_model => $self -> {'models'} -> [$model_number-1],
1098 base_model => $self -> {'base_models'} -> [$model_number-1],
1099 relations => $self -> {'relations'}[$model_number-1],
1100 included_relations =>
1101 $self -> {'included_relations'}[$model_number-1]);
1102 # Create a modelfit tool for all the models of this step.
1103 # This is the last setup part before running the step.
1104 my %subargs = ();
1105 if ( defined $self -> {'subtool_arguments'} ) {
1106 %subargs = %{$self -> {'subtool_arguments'}};
1107 }
1108
1109 push( @{$self -> {'tools'}},
1110 tool::modelfit -> new
1111 ( reference_object => $self,
1112 prepared_models => undef,
1113 subtools => undef,
1114 results => undef,
1115 raw_results => undef,
1116 _raw_results_callback => $self ->
1117 _raw_results_callback( model_number => $model_number ),
1118 models => $self -> {'prepared_models'}[$model_number-1]{'own'},
1119 threads => $mfit_threads,
1120 logfile => $self -> {'directory'}."/modelfit".$model_number.".log",
1121 base_directory => $self -> {'directory'},
1122 directory => $self -> {'directory'}.'/modelfit_dir'.$model_number,
1123 parent_threads => $own_threads,
1124 parent_tool_id => $self -> {'tool_id'},
1125 top_tool => 0,
1126 %subargs ) );
1127 }
1128 end modelfit_setup
1129
1130 # }}}
1131
1132 # {{{ modelfit_analyze
1133
1134 start modelfit_analyze
1135 {
1136 # Own_threads is used to set parent_threads for child tools
1137 my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
1138 $self -> {'threads'} -> [0]:$self -> {'threads'};
1139 # More threads than models?
1140 my $num = scalar @{$self -> {'models'}};
1141 $own_threads = $num if ( $own_threads > $num );
1142
1143 my @results = @{$self -> {'results'}};
1144
1145 my @minimizations;
1146 for ( my $i = 0; $i < scalar @{$self -> {'prepared_models'}[$model_number-1]{'own'}}; $i++ ) {
1147 my $term = $self -> {'prepared_models'}[$model_number-1]{'own'}[$i] ->
1148 outputs -> [0] -> minimization_successful;
1149 push( @minimizations, $term->[0][0] );
1150 }
1151 # This loop checks the termination of all model fits
1152 foreach my $mod_return_type ( @{$results[$model_number-1]{'subtools'}[0]} ){
1153 my $crash = 0;
1154 my $crash_name;
1155 foreach my $type ( @{$mod_return_type -> {'own'}} ) {
1156 if ( $type -> {'name'} eq 'minimization_successful' ){
1157 for ( my $i = 0; $i < scalar @{$type -> {'values'}}; $i++ ) {
1158 for ( my $j = 0; $j < scalar @{$type -> {'values'}[$i]}; $j++ ) {
1159 if ( $self -> {'picky'} ) {
1160 # Did minimization just loose one dimension?
1161 # if ( not defined $type -> {'values'}[$i][$j][0] or
1162 # $type -> {'values'}[$i][$j][0] != 1 ) {
1163 if ( not defined $type -> {'values'}[$i][$j] or
1164 $type -> {'values'}[$i][$j] != 1 ) {
1165 $crash = 2;
1166 }
1167 } else {
1168 # if ( not defined $type -> {'values'}[$i][$j][0] or
1169 # $type -> {'values'}[$i][$j][0] <= 0 ) {
1170 if ( not defined $type -> {'values'}[$i][$j] or
1171 $type -> {'values'}[$i][$j] <= 0 ) {
1172 $crash = 1;
1173 }
1174 }
1175 }
1176 }
1177 }
1178 $crash_name = $type -> {'values'}[0] if ( $type -> {'name'} eq 'filename' );
1179 }
1180 if ( $crash == 2 and $self -> {'abort_on_fail'} ) {
1181 die "The execution of the modelfile ",
1182 $crash_name," failed (picky==1), aborting SCM, sorry\n";
1183 } elsif ( $crash == 1 and $self -> {'abort_on_fail'} ) {
1184 die "The execution of the modelfile ",
1185 $crash_name," failed (picky==0), aborting SCM, sorry\n";
1186 }
1187 }
1188
1189 my %return_section;
1190 $return_section{'name'} = 'ofv.drop';
1191 $return_section{'values'} = [];
1192 $return_section{'labels'} = [];
1193
1194 # Collect the drops. The $i is needed since the prepared
1195 # models are not indexed by parameters and covariate
1196 # combination
1197 my $base_ofv;
1198 if ( defined $self -> {'base_criteria_values'} and
1199 defined $self -> {'base_criteria_values'} ->
1200 [$model_number-1] -> {'ofv'}) {
1201 $base_ofv = $self -> {'base_criteria_values'} ->
1202 [$model_number-1] -> {'ofv'}
1203 } else {
1204 $base_ofv = $self -> {'models'}[$model_number-1] ->
1205 outputs -> [0] -> ofv -> [0][0];
1206 }
1207 my $i = 0;
1208
1209 foreach my $parameter ( sort keys %{$self -> {'relations'}[$model_number-1]} ) {
1210 foreach my $covariate ( sort keys %{$self -> {'relations'}[$model_number-1]{$parameter}} ) {
1211 # Is this covariate continuous or not?
1212 my $continuous = 1;
1213 foreach my $cat ( @{$self -> {'categorical_covariates'}} ) {
1214 $continuous = 0 if ( $covariate eq $cat );
1215 }
1216 my @valid_states;
1217 if ( $continuous ) {
1218 @valid_states = @{$self -> {'valid_states'}{'continuous'}};
1219 } else {
1220 @valid_states = @{$self -> {'valid_states'}{'categorical'}};
1221 }
1222
1223 my $state;
1224 # Important: just setting $state to $self->incl_rel....{'state'} initiates
1225 # included_relations for this parameter and covariate. Avoid this.
1226 if ( defined $self -> {'included_relations'}[$model_number-1]{$parameter}{$covariate} ) {
1227 $state = $self -> {'included_relations'}[$model_number-1]{$parameter}{$covariate}{'state'};
1228 }
1229 $state = defined $state ? $state : $valid_states[0];
1230
1231 next if ( ( $self -> {'search_direction'} eq 'forward' and
1232 $state == $valid_states[$#valid_states] ) or
1233 # $state == $self -> {'valid_states'} ->
1234 # [scalar @{$self -> {'valid_states'}}-1] ) or
1235 ( $self -> {'search_direction'} eq 'backward' and
1236 $state == $valid_states[0] ) );
1237 # $state == $self -> {'valid_states'} -> [0] ) );
1238 my $new_ofv = $self -> {'prepared_models'}[$model_number-1]{'own'}[$i] -> outputs -> [0] -> ofv -> [0][0];
1239 $new_ofv = $base_ofv unless ( defined $new_ofv );
1240 # Only one problem and one sub problem
1241 push( @{$return_section{'values'}[0][0]}, $base_ofv - $new_ofv );
1242 $i++;
1243 }
1244 }
1245 push( @{$self -> {'results'}[$model_number-1]{'own'}},
1246 \%return_section );
1247 my @ofv_changes;
1248 foreach my $par ( sort keys %{$self -> {'test_relations'}[$model_number-1]} ) {
1249 foreach my $cov ( @{$self -> {'test_relations'}[$model_number-1]{$par}} ){
1250 push( @ofv_changes,
1251 $self -> {'relations'}[$model_number-1]{$par}{$cov}{'ofv_changes'} );
1252 }
1253 }
1254 my ( $chosen_parameter,
1255 $chosen_covariate,
1256 $sign_models_ref,
1257 $test_vals_ref,
1258 $criterion,
1259 $test_log_ref,
1260 $new_base_crit_val_ref,
1261 $short_log_text,
1262 $short_log_ref ) = ( undef, undef, undef, undef,
1263 undef, undef, undef, undef, undef );
1264 ( $self -> {'resulting_model'},
1265 $chosen_parameter,
1266 $chosen_covariate,
1267 $sign_models_ref,
1268 $test_vals_ref,
1269 $criterion,
1270 $test_log_ref,
1271 $new_base_crit_val_ref,
1272 $short_log_text,
1273 $short_log_ref )
1274 = $self -> {'gof'} -> ( $self,
1275 $self -> {'search_direction'},
1276 $self -> {'models'} -> [$model_number-1],
1277 $model_number,
1278 \@ofv_changes);
1279
1280 # Print a short summary of the step (All logging should be kept in a log-method in the future)
1281 open( LOG, ">>".$self -> {'short_logfile'} -> [$model_number-1] );
1282 print LOG $short_log_text;
1283 close( LOG );
1284
1285 my %return_section;
1286 $return_section{'name'} = 'relation.chosen.in.step';
1287 $return_section{'values'} = [];
1288 $return_section{'labels'} = [];
1289 if ( defined $chosen_parameter and defined $chosen_covariate ) {
1290 $return_section{'values'}[0][0] = $chosen_parameter.$chosen_covariate;
1291 my $task = $self -> {'search_direction'} eq 'forward' ? 'Adding' : 'Removing';
1292 ui -> print( category => 'scm',
1293 message => "$task $chosen_covariate on $chosen_parameter" )
1294 unless $self -> {'parent_threads'} > 1;
1295 }
1296 push( @{$self -> {'results'}[$model_number-1]{'own'}},
1297 \%return_section );
1298
1299 # my $i = 0;
1300 # if ( defined $chosen_parameter and defined $chosen_covariate ) {
1301 # while ( my( $parameter, $cov_ref ) = each %{$self -> {'relations'}[$model_number-1]} ) {
1302 # while ( my( $covariate, $state_ref ) = each %{$self -> {'relations'}[$model_number-1]{$parameter}} ) {
1303 # if ( $parameter eq $chosen_parameter and
1304 # $covariate eq $chosen_covariate ) {
1305 # $self -> {'resulting_model'} = $self -> {'prepared_models'}[$model_number-1]{'own'}[$i];
1306 # }
1307 # $i++;
1308 # }
1309 # }
1310 # }
1311
1312
1313
1314 my $final_model;
1315 if ( defined $self -> {'resulting_model'} ) {
1316 # Promote and log the included relation:
1317 # Is this covariate continuous or not?
1318 my $continuous = 1;
1319 foreach my $cat ( @{$self -> {'categorical_covariates'}} ) {
1320 $continuous = 0 if ( $chosen_covariate eq $cat );
1321 }
1322 my @valid_states;
1323 if ( $continuous ) {
1324 @valid_states = @{$self -> {'valid_states'}{'continuous'}};
1325 } else {
1326 @valid_states = @{$self -> {'valid_states'}{'categorical'}};
1327 }
1328
1329 my $state;
1330 # Important: just setting $state to $self->incl_rel....{'state'} initiates
1331 # included_relations for this parameter and covariate. Avoid this.
1332 if ( defined $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}
1333 {$chosen_covariate} ) {
1334 $state = $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}
1335 {$chosen_covariate}{'state'};
1336 }
1337 $state = defined $state ? $state : $valid_states[0];
1338 # my $state = $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}
1339 # {$chosen_covariate}{'state'};
1340 if ( $self -> {'search_direction'} eq 'forward' ) {
1341 my $flag = 0;
1342 for( my $s_idx = 0; $s_idx <= $#valid_states; $s_idx++ ) {
1343 if ( $flag ) {
1344 $state = $valid_states[$s_idx];
1345 last;
1346 }
1347 $flag = 1 if( $state == $valid_states[$s_idx] );
1348 }
1349 } elsif ( $self -> {'search_direction'} eq 'backward' ) {
1350 my $flag = 0;
1351 for( my $s_idx = $#valid_states; $s_idx >= 0; $s_idx-- ) {
1352 if ( $flag ) {
1353 $state = $valid_states[$s_idx];
1354 last;
1355 }
1356 $flag = 1 if( $state == $valid_states[$s_idx] );
1357 }
1358 }
1359 # If the state is 1 (not included); remove the relation.
1360 if ( $state == 1 ) {
1361 # Check if this relation is the last for the parameter
1362 if ( scalar keys %{$self -> {'included_relations'}[$model_number-1]
1363 {$chosen_parameter}} == 1 ) {
1364 delete $self -> {'included_relations'}[$model_number-1]{$chosen_parameter};
1365 } else {
1366 delete $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}
1367 {$chosen_covariate};
1368 }
1369 } else {
1370 $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'state'} = $state;
1371 $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'theta_estimates'} =
1372 $self -> {'resulting_model'} -> outputs -> [0] -> thetas -> [0][0];
1373 $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'code'} =
1374 $self -> {'relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'code'}{$state};
1375 $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'nthetas'} =
1376 $self -> {'relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'nthetas'}{$state};
1377 $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'inits'} =
1378 $self -> {'relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'inits'}{$state};
1379 $self -> {'included_relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'bounds'} =
1380 $self -> {'relations'}[$model_number-1]{$chosen_parameter}{$chosen_covariate}{'bounds'}{$state};
1381 }
1382
1383 $self -> write_log
1384 ( direction => $self -> {'search_direction'},
1385 logfile => $self -> {'logfile'}[$model_number-1],
1386 included_relations => $self -> {'included_relations'}[$model_number-1],
1387 chosen_parameter => $chosen_parameter,
1388 chosen_covariate => $chosen_covariate,
1389 results => $self -> {'results'}[$model_number-1]{'own'},
1390 criterion => $criterion,
1391 test_log => $test_log_ref,
1392 append => $self -> {'append_log'} );
1393
1394 # Check if there still are some states to test
1395 # To do, maybe fixed 2005-01-12: this is only valid for the forward search, fix the backward
1396 my $still_one_left = 0;
1397 foreach my $par ( sort keys %{$self -> {'test_relations'}
1398 [$model_number-1]} ) {
1399 foreach my $cov ( @{$self -> {'test_relations'}[$model_number-1]{$par}} ){
1400 # bad coding, I know. /Lasse
1401 my $kind = 'continuous';
1402 foreach my $cat ( @{$self -> {'categorical_covariates'}} ) {
1403 $kind = 'categorical' if ( $cov eq $cat );
1404 }
1405 if ( defined $self -> {'included_relations'}[$model_number-1]{$par}{$cov} ) {
1406 $still_one_left = 1 unless
1407 ( $self -> {'included_relations'}[$model_number-1]{$par}{$cov}{'state'} ==
1408 $self -> {'valid_states'}{$kind} ->
1409 [ scalar @{$self -> {'valid_states'}{$kind}} - 1] );
1410 } else {
1411 $still_one_left = 1;
1412 }
1413 }
1414 }
1415
1416 my ( $returns, $prep_models );
1417 if ( $still_one_left ) {
1418 my $cpu_time = defined $self -> {'cpu_time'} ? int($self -> {'cpu_time'}*1.2) : undef;
1419 my $internal_scm =
1420 tool::scm ->
1421 new( reference_object => $self,
1422 cpu_time => $cpu_time,
1423 covariate_statistics => [$self -> {'covariate_statistics'}[$model_number-1]],
1424 directory => $self -> {'directory'}.'/scm_dir'.$model_number,
1425 models => [$self -> {'models'}[$model_number-1]],
1426 relations => [$self -> {'relations'}[$model_number-1]],
1427 base_models => [$self -> {'resulting_model'}],
1428 included_relations => [$self -> {'included_relations'}[$model_number-1]],
1429 fix => [$self -> {'fix'}[$model_number-1]],
1430 append_log => 1,
1431 step_number => ($self -> {'step_number'} + 1),
1432 raw_results_file => [$self -> {'raw_results_file'}[$model_number-1]],
1433 logfile => [$self -> {'logfile'}[$model_number-1]],
1434 base_criteria_values => [$new_base_crit_val_ref],
1435 parent_tool_id => $self -> {'tool_id'},
1436 parent_threads => $own_threads ,
1437 top_tool => 0,
1438 config_file => undef,
1439 resulting_model => undef,
1440 tools => undef,
1441 prepared_models => undef,
1442 results => undef );
1443
1444 ui -> print( category => 'scm',
1445 message => "Taking a step " . $self -> {'search_direction'} )
1446 unless $self -> {'parent_threads'} > 1;
1447 $internal_scm -> run;
1448 $returns = $internal_scm -> {'results'};
1449 $prep_models = $internal_scm -> {'prepared_models'};
1450 ui -> print( category => 'scm',
1451 message => $self -> {'search_direction'} . " step done." )
1452 unless $self -> {'parent_threads'} > 1;
1453
1454 foreach my $return ( @{$returns ->[0]{'own'}} ) {
1455 if ( $return -> {'name'} eq 'base.criteria.values' ){
1456 $self -> base_criteria_values( $return -> {'values'} );
1457 }
1458 }
1459 } else {
1460 my @tmp_ret;
1461 $tmp_ret[0]{'own'}[0]{'name'} = 'final.model';
1462 $tmp_ret[0]{'own'}[0]{'values'}[0][0] = 'basic.model';
1463 $returns = \@tmp_ret;
1464
1465 $self -> base_criteria_values( [$new_base_crit_val_ref] );
1466 }
1467
1468 if ( defined $prep_models ) {
1469 debug -> warn( level => 2,
1470 message => " have called internal scm " .
1471 scalar @{$prep_models} );
1472
1473 # Enclose $prep_models in array ref to reflect the
1474 # per-tool level, even though a modelfit does not
1475 # prepare any models itself
1476 #push ( @{$self -> {'prepared_models'}[$model_number-1]{'subtools'}}, $prep_models -> [0]);
1477
1478 # Push the results of the internal scm on the results attribute:
1479 #push( @{$self -> {'results'}[$model_number-1]{'subtools'}}, $returns );
1480 } else {
1481 debug -> warn( level => 2,
1482 message => " no prep_models defined from internal scm " );
1483 }
1484
1485 # set final model to this steps' best model if the internal scm returned 'basic_model'.
1486 foreach my $return ( @{$returns ->[0]{'own'}} ) {
1487 $final_model = $return -> {'values'}[0][0] if ( $return -> {'name'} eq 'final.model' );
1488 }
1489
1490 if ( not defined $final_model or $final_model eq 'basic_model' ) {
1491 $final_model = $self -> {'resulting_model'};
1492 my ( $fdir, $fname ) = OSspecific::absolute_path( '.',$self -> {'final_model_name'}[$model_number-1] );
1493 $final_model -> filename( $fname );
1494 $final_model -> directory( $fdir );
1495 $fname =~ s/\.mod/\.lst/;
1496 cp( $self -> {'resulting_model'} -> outputfile, "$fdir$fname" );
1497 $final_model -> outputfile( "$fdir$fname" );
1498 $final_model -> _write;
1499 $self -> write_log
1500 ( direction => $self -> {'search_direction'},
1501 logfile => $self -> {'short_logfile'}[$model_number-1],
1502 included_relations => $self -> {'included_relations'}[$model_number-1],
1503 chosen_parameter => $chosen_parameter,
1504 chosen_covariate => $chosen_covariate,
1505 results => $self -> {'results'}[$model_number-1]{'own'},
1506 criterion => $criterion,
1507 test_log => $test_log_ref,
1508 append => 1 );
1509 }
1510 } else {
1511 # This is the last step.
1512 $self -> write_log ( direction => $self -> {'search_direction'},
1513 logfile => $self -> {'logfile'}[$model_number-1],
1514 results => $self -> {'results'}[$model_number-1]{'own'},
1515 append => $self -> {'append_log'} );
1516
1517 # Leave base_criteria_values as they are
1518 }
1519
1520 my %return_section;
1521 $return_section{'name'} = 'base.criteria.values';
1522 $return_section{'values'} = $self -> base_criteria_values;
1523 $return_section{'labels'} = undef;
1524 push( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
1525
1526 my %return_section;
1527 $return_section{'name'} = 'included.relations';
1528 $return_section{'values'} = $self -> {'included_relations'}[$model_number-1];
1529 $return_section{'labels'} = undef;
1530 push( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
1531
1532
1533 # This loop tries to minimize the data written to disc.
1534 for ( my $i = 0; $i < scalar @{$self -> {'prepared_models'}[$model_number-1]{'own'}}; $i++ ) {
1535 $self -> {'prepared_models'}[$model_number-1]{'own'}[$i] -> {'datas'} = undef;
1536 $self -> {'prepared_models'}[$model_number-1]{'own'}[$i] -> {'outputs'} = undef;
1537 }
1538
1539 my %return_section;
1540 $return_section{'name'} = 'final.model';
1541 $return_section{'values'}[0][0] = $final_model;
1542 $return_section{'labels'} = undef;
1543 push( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
1544
1545 debug -> warn( level => 2,
1546 message => "Finished in modelfit_analyze!" );
1547 }
1548 end modelfit_analyze
1549
1550 # }}} modelfit_analyze
1551
1552 # {{{ modelfit_post_fork_analyze
1553
1554 start modelfit_post_fork_analyze
1555 {
1556 # Check the thread number of this tool level:
1557 my $threads = ref( $self -> {'threads'} ) eq 'ARRAY' ?
1558 $self -> {'threads'} -> [0]:$self -> {'threads'};
1559 # More threads than models?
1560 my @models = @{$self -> {'models'}};
1561 $threads = $#models + 1 if ( $threads > $#models + 1);
1562
1563 # It is not necessary to collect the included relations from
1564 # the runs if no parallelism has been used or if only one run
1565 # has been performed. Now there is allways parallelism! see
1566 # tool::run and Readme.txt
1567
1568 my @included_relations = ();
1569 foreach my $return ( @{$self -> {'results'}->[0]{'own'}} ) {
1570 if ( $return -> {'name'} eq 'included.relations' ){
1571 $self -> {'included_relations'} = [$return -> {'values'}];
1572 }
1573 if ( $return -> {'name'} eq 'base.criteria.values' ){
1574 $self -> base_criteria_values( $return -> {'values'} );
1575 }
1576 }
1577 }
1578 end modelfit_post_fork_analyze
1579
1580 # }}}
1581
1582 # {{{ gof_crc
1583
1584 start gof_crc
1585 {
1586 my ($direction, $basic_model, $model_number) = @_;
1587 my ( $chosen_parameter, $chosen_covariate );
1588 my %p_values;
1589
1590 #For unlimited stepping: p = 100% /JR
1591 $p_values{'1'} = {1=>0,
1592 2=>0,
1593 3=>0,
1594 4=>0,
1595 5=>0,
1596 6=>0,
1597 7=>0,
1598 8=>0,
1599 9=>0,
1600 10=>0};
1601 ## Set the critical values to use
1602 ## p= 0.05
1603 $p_values{'0.05'} = {1=>3.84,
1604 2=>5.99,
1605 3=>7.81,
1606 4=>9.49,
1607 5=>11.07,
1608 6=>12.59,
1609 7=>14.07,
1610 8=>15.51,
1611 9=>16.92,
1612 10=>18.31};
1613 ## p=0.01
1614 $p_values{'0.01'} = {1=>6.63,
1615 2=>9.21,
1616 3=>11.34,
1617 4=>13.28,
1618 5=>15.09,
1619 6=>16.81,
1620 7=>18.48,
1621 8=>20.09,
1622 9=>21.67,
1623 10=>23.21};
1624 ## p=0.005
1625 $p_values{'0.005'} = {1=>7.88,
1626 2=>10.60,
1627 3=>12.84,
1628 4=>14.86,
1629 5=>16.75,
1630 6=>18.55,
1631 7=>20.28,
1632 8=>21.95,
1633 9=>23.59,
1634 10=>25.19};
1635 ## p=0.001
1636 $p_values{'0.001'} = {1=>10.83,
1637 2=>13.82,
1638 3=>16.27,
1639 4=>18.47,
1640 5=>20.52,
1641 6=>22.46,
1642 7=>24.32,
1643 8=>26.12,
1644 9=>27.88,
1645 10=>29.59};
1646
1647 my $resulting_model;
1648 my @ofv_drops;
1649 foreach my $return ( @{$self -> {'results'}[$model_number-1]{'own'}} ) {
1650 @ofv_drops = @{$return -> {'values'}[0][0]} if ( $return -> {'name'} eq 'ofv.drop' );
1651 }
1652 my $base_n_param = $self ->
1653 {'models'}[$model_number-1] -> nthetas( problem_number => 1 );
1654 print "Orig model NPARAM: $base_n_param\n";
1655 if ( defined $self -> {'included_relations'}[$model_number-1] ) {
1656 my %included_relations = %{$self -> {'included_relations'}[$model_number-1]};
1657 foreach my $incl_par ( sort keys %included_relations ) {
1658 foreach my $incl_cov ( sort keys %{$included_relations{$incl_par}} ) {
1659 $base_n_param += $included_relations{$incl_par}{$incl_cov}{'nthetas'};
1660 print "Adding $included_relations{$incl_par}{$incl_cov}{'nthetas'}\n";
1661 }
1662 }
1663 }
1664 print "STEP model NPARAM: $base_n_param\n";
1665 my %sign_drops;
1666 my @step_relations = @{$self -> {'step_relations'}};
1667 for ( my $i = 0; $i <= $#step_relations; $i++ ) {
1668 my $n_param_diff =
1669 $self -> {'prepared_models'}[$model_number-1]{'own'}[$i] ->
1670 nthetas( problem_number => 1 ) - $base_n_param;
1671 print "Diff in parameters for model $i: $n_param_diff\n";
1672 if( $self -> {'search_direction'} eq 'forward'){
1673 if( $ofv_drops[$i] > $p_values{'0.05'}{$n_param_diff} ){ # Significant ?
1674 print "FORW SIGNIFICANT! OFV_DROP=$ofv_drops[$i] ,NPARM_DIFF= $n_param_diff ",
1675 $step_relations[$i]{'parameter'},$step_relations[$i]{'covariate'},"\n";
1676 $sign_drops{$i} = $ofv_drops[$i];
1677 }
1678 } else {
1679 print "Testing against p-value: $p_values{'0.01'}{-$n_param_diff}\n";
1680 if( $ofv_drops[$i] > -$p_values{'0.01'}{-$n_param_diff} ){ # Significant ?
1681 print "BACKW SIGNIFICANT! OFV_DROP=$ofv_drops[$i] ,NPARM_DIFF= $n_param_diff ",
1682 $step_relations[$i]{'parameter'},$step_relations[$i]{'covariate'},"\n";
1683 $sign_drops{$i} = $ofv_drops[$i];
1684 }
1685 }
1686 }
1687
1688 my %omega_number;
1689 $omega_number{'iiv_CL'} = 1;
1690 $omega_number{'iiv_V'} = 2;
1691 my %crc_changes;
1692 foreach my $crit ( 'iiv_CL', 'iiv_V' ) {
1693 foreach my $model ( @{$self -> {'prepared_models'}[$model_number-1]{'own'}} ) {
1694 $crc_changes{$crit} = $model -> outputs -> [0] -> omegas -> [0][0][$omega_number{$crit}];
1695 }
1696 }
1697
1698 my %clin_sign;
1699 # Clinical relevance
1700 for ( my $i = 0; $i <= $#step_relations; $i++ ) {
1701 next unless abs($sign_drops{$i});
1702 foreach my $crit ( 'iiv_CL', 'iiv_V' ) {
1703 # if ( $self -> {'search_direction'} eq 'forward'){
1704 # if (
1705 # $clin_sign{$crit}{$i}
1706 }
1707 }
1708 return ( $resulting_model, $chosen_parameter, $chosen_covariate,
1709 'CRC', {'BASE_MODEL_TVCL' => $base_cl[0][0],
1710 'CHOSEN_MODEL_TVCL' => $chosen_cl } );
1711 }
1712 end gof_crc
1713
1714 # }}}
1715
1716 # {{{ gof_ofv
1717
1718 start gof_ofv
1719 {
1720 my ( $direction, $basic_model, $model_number, $ofv_ch_ref ) = @_;
1721 my @ofv_changes = @{$ofv_ch_ref};
1722 my $base_ofv;
1723 if ( defined $self -> {'base_criteria_values'} and
1724 defined $self -> {'base_criteria_values'} ->
1725 [$model_number-1] -> {'ofv'} ) {
1726 $base_ofv = $self -> {'base_criteria_values'} ->
1727 [$model_number-1] -> {'ofv'};
1728 } elsif ( $direction eq 'backward' ) {
1729 'debug' -> die( message => "Backward search needs a 'base' OFV estimate" );
1730 } else {
1731 if ( defined $self -> {'models'} -> [$model_number-1] -> outputs -> [0] ->
1732 ofv -> [0][0] ) {
1733 $base_ofv = $self -> {'models'} -> [$model_number-1] -> outputs -> [0] ->
1734 ofv -> [0][0];
1735 } else {
1736 'debug' -> die( message => "OFV estimates not available from model" .
1737 $self -> {'models'} -> [$model_number-1] -> full_name );
1738 }
1739 }
1740
1741 my @models = @{$self -> {'prepared_models'} -> [$model_number-1]{'own'}};
1742 my @ofvs;
1743 foreach my $model ( @models ) {
1744 push( @ofvs, $model -> outputs -> [0] -> ofv ->[0][0] );
1745 }
1746
1747 my ( @ofv_drops, @log_texts );
1748 my $base_n_param = $self ->
1749 {'models'}[$model_number-1] -> nthetas( problem_number => 1 );
1750 if ( defined $self -> {'included_relations'}[$model_number-1] ) {
1751 my %included_relations = %{$self -> {'included_relations'}[$model_number-1]};
1752 # while ( my( $incl_par, $cov_ref ) = each %included_relations ) {
1753 # while ( my( $incl_cov, $state_ref ) = each %{$included_relations{$incl_par}} ) {
1754 foreach my $incl_par ( sort keys %included_relations ) {
1755 foreach my $incl_cov ( sort keys %{$included_relations{$incl_par}} ) {
1756 $base_n_param += $included_relations{$incl_par}{$incl_cov}{'nthetas'};
1757 }
1758 }
1759 }
1760
1761 open( LOG, ">>".$self -> {'logfile'} -> [$model_number-1] );
1762 my $un = $direction eq 'backward' ? '(UN)' : '';
1763 print LOG sprintf("%-8s",'MODEL'),
1764 sprintf("%12s",'TEST NAME'),
1765 sprintf("%12s",'BASE VAL'),
1766 sprintf("%12s",'NEW VAL'),
1767 sprintf("%50s",'TEST VAL (DROP)'),
1768 sprintf("%10s","GOAL"),
1769 sprintf("%14s"," $un".'SIGNIFICANT'),"\n";
1770 my ( %drop_sign, @n_param_diffs );
1771 my @step_relations = @{$self -> {'step_relations'}};
1772 for ( my $i = 0; $i <= $#step_relations; $i++ ) {
1773 my $n_param_diff =
1774 $self -> {'prepared_models'}[$model_number-1]{'own'}[$i] ->
1775 nthetas( problem_number => 1 ) - $base_n_param;
1776 push( @n_param_diffs, $n_param_diff );
1777 my $change = $direction eq 'forward' ?
1778 $ofv_changes[$model_number-1]{$n_param_diff} :
1779 -$ofv_changes[$model_number-1]{-$n_param_diff};
1780 my $test_val;
1781 my $ofv;
1782 unless( defined( $ofvs[$i] ) ){
1783 $test_val = ' ' x 43 . 'FAILED';
1784 $ofv = ' ' x 4 . 'FAILED';
1785 } else {
1786 $test_val = $base_ofv - $ofvs[$i];
1787 $test_val = sprintf("%47.5f",$test_val);
1788 $ofv = sprintf("%12.5f",$ofvs[$i])
1789 }
1790
1791 push ( @ofv_drops, $test_val );
1792 my $log_text = sprintf("%-8s",$step_relations[$i]{'parameter'}.
1793 $step_relations[$i]{'covariate'}.'-'.
1794 $step_relations[$i]{'state'}).
1795 sprintf("%12s",'OFV ').
1796 sprintf("%12.5f",$base_ofv).
1797 $ofv.
1798 $test_val. ' >'.
1799 sprintf("%10.5f",$change);
1800 print LOG $log_text;
1801 # Significant ?
1802 if( defined $ofvs[$i] and $test_val > $change ){
1803 my $yes_text = sprintf("%12s",'YES! ');
1804 $log_text = $log_text.$yes_text;
1805 print LOG $yes_text;
1806 $drop_sign{$i} = 1;
1807 }
1808 print LOG "\n";
1809 push( @log_texts, $log_text."\n" );
1810 }
1811 print LOG "\n";
1812 close ( LOG );
1813
1814 my ( %sign, %l_text );
1815 for ( my $i = 0; $i <= $#models; $i++ ) {
1816 my $od = defined $drop_sign{$i} ? 1 : 0;
1817 $sign{$i} = $ofv_drops[$i] if ($od);
1818 $l_text{$i} = $log_texts[$i] if ($od);
1819 }
1820
1821 my $chosen_ofv;
1822 my $resulting_model;
1823 my ( $chosen_parameter, $chosen_covariate, $chosen_log_text );
1824 if ( scalar keys %sign > 0 ) {
1825 my @sorted_ids = sort { $sign{$b} <=> $sign{$a} } keys %sign;
1826 $resulting_model = $self ->
1827 {'prepared_models'}[$model_number-1]{'own'}[$sorted_ids[0]];
1828 $chosen_ofv = $self ->
1829 {'prepared_models'}[$model_number-1]{'own'}[$sorted_ids[0]] ->
1830 outputs -> [0] -> ofv -> [0][0];
1831 $chosen_log_text = $l_text{$sorted_ids[0]};
1832 $chosen_parameter = $step_relations[$sorted_ids[0]]{'parameter'};
1833 $chosen_covariate = $step_relations[$sorted_ids[0]]{'covariate'};;
1834 }
1835
1836 return ( $resulting_model,
1837 $chosen_parameter,
1838 $chosen_covariate,
1839 \%drop_sign,
1840 \@ofv_drops,
1841 'OFV',
1842 {'BASE_MODEL_OFV' => $base_ofv,
1843 'CHOSEN_MODEL_OFV' => $chosen_ofv },
1844 {'ofv' => $chosen_ofv},
1845 # \@n_param_diffs,
1846 $chosen_log_text,
1847 \%l_text );
1848 }
1849 end gof_ofv
1850
1851 # }}}
1852
1853 # {{{ cov_code_print
1854
1855 start cov_code_print
1856 {
1857 my $cov = $self -> {'labels'}{$label};
1858 my $basic_model = @{$self -> {'models'}}[0];
1859 my $data = @{$basic_model -> datas}[0];
1860 my $basic_theta_num;
1861 my $real_theta_num;
1862 # Print the start line
1863 push @TX, "\n;;; $label".$level." START\n";
1864 my $comment;
1865 if( $print_oldstyle ){
1866 $comment = ";;; ";
1867 $basic_theta_num = 0;
1868 } else {
1869 $comment = "";
1870
1871 # This assumes one problem per file:
1872 $basic_theta_num = $basic_model -> nthetas;
1873 }
1874 if( $self -> {'scm_file'} and -e $self -> {'scm_file'}){
1875 if( defined $self -> {'spec_code'}{$label}{$level} ){
1876 foreach my $code_row( @{$self -> {'spec_code'}{$label}{$level}} ){
1877 $code_row =~ s/^\s*;;;\s*//;
1878 if( $code_row =~ /THETA/ ){
1879 $theta_numbering++;
1880 $real_theta_num = $basic_theta_num + $theta_numbering;
1881 $code_row =~ s/THETA\(\d+\)/THETA($real_theta_num)/;
1882 }
1883 push @TX, $code_row;
1884 }
1885 } else {
1886 if( $level == 1 ){
1887 push @TX, "$label = 0\n";
1888 }
1889 }
1890 } else {
1891 # IF at first level
1892 if( $level == 1 ){
1893 push @TX, "$comment$label = 0\n";
1894 } else {
1895 if($self -> {'label_type'}{$label} eq "con") {
1896 # {{{ Code for continuous covariates
1897 # {{{ If at second level
1898
1899 if($level == 2) {
1900 my $med = $self -> {'median'}{$cov};
1901 $med = sprintf "%6.2f", $med;
1902 $med =~ s/\s*//;
1903 if($data -> have_missing_data('column_head' => $cov)) {
1904 my $miss_data_lbl = $self -> {'missing_data_token'};
1905 push @TX, $comment."IF($cov.EQ.$miss_data_lbl) THEN\n";
1906 push @TX, "$comment $label = 0\n";
1907 push @TX, $comment."ELSE\n";
1908 $theta_numbering++;
1909 $real_theta_num = $basic_theta_num + $theta_numbering;
1910 my $sign = '-';
1911 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
1912 push @TX, "$comment $label = THETA($real_theta_num)*($cov $sign $med)\n";
1913 push @TX, $comment."ENDIF\n";
1914 } else {
1915 $theta_numbering++;
1916 $real_theta_num = $basic_theta_num + $theta_numbering;
1917 my $sign = '-';
1918 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
1919 push @TX, "$comment$label = THETA($real_theta_num)*($cov $sign $med)\n";
1920 }
1921 }
1922
1923 # }}}
1924 # {{{ If at third level
1925
1926 if($level == 3) {
1927 my $med = $self -> {'median'}{$cov};
1928 $med = sprintf "%6.2f", $med;
1929 $med =~ s/\s*//;
1930 $theta_numbering++;
1931 $real_theta_num = $basic_theta_num + $theta_numbering;
1932 my $sign = '-';
1933 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
1934 push @TX, $comment."IF($cov.LE.$med) $label = THETA($real_theta_num)*($cov $sign $med)\n";
1935 $theta_numbering++;
1936 $real_theta_num = $basic_theta_num + $theta_numbering;
1937 my $sign = '-';
1938 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
1939 push @TX, $comment."IF($cov.GT.$med) $label = THETA($real_theta_num)*($cov $sign $med)\n";
1940 if ($data -> have_missing_data('column_head' => $cov)) {
1941 push @TX, $comment."IF($cov.EQ.$self -> {'missing_data_token'}) $label = 0\n";
1942 }
1943 }
1944
1945 # }}}
1946 # }}}
1947 } else {
1948 # {{{ Code for categorical variables
1949
1950 my %levels;
1951
1952 $levels{$cov} = $data -> factors( 'return_occurences' => 1,
1953 'unique_in_individual' => 0,
1954 'column_head' => $cov );
1955 my @levels = sort {$a<=>$b} keys %{$levels{$cov}};
1956 my $numlvs = scalar @levels;
1957 $numlvs = $numlvs -1 if $data -> have_missing_data('column_head' => $cov);
1958 # If there are two levels and -f center it
1959 if ($self -> {'fix'}[$model_number] and $numlvs == 2) {
1960 # {{{ If fixing and at second level
1961
1962 my ($lev1, $lev2);
1963 my %fractions = $data -> fractions( 'unique_in_individual' => 1,
1964 'column_head' => $cov,
1965 'ignore_missing' => 1);
1966 my $fraction;
1967 if ($fractions{$levels[0]} > $fractions{$levels[1]}) {
1968 $lev1 = 0;
1969 $lev2 = 1;
1970 $fraction = $fractions{$levels[0]};
1971 } else {
1972 $lev1 = 1;
1973 $lev2 = 0;
1974 $fraction = $fractions{$levels[0]};
1975 }
1976 $fraction = sprintf("%4.2f", $fraction);
1977 if ($data -> have_missing_data('column_head' => $cov)) {
1978 push @TX, $comment."IF($cov.EQ.$levels[$lev1].OR.$cov.EQ.$self -> {'missing_data_token'})";
1979 $theta_numbering++;
1980 $real_theta_num = $basic_theta_num + $theta_numbering;
1981 push @TX, "$comment$label = -1*THETA($real_theta_num)*(1-$fraction)\n";
1982 push @TX, $comment."IF($cov.EQ.$levels[$lev2])";
1983 $theta_numbering++;
1984 $real_theta_num = $basic_theta_num + $theta_numbering;
1985 push @TX, "$comment$label = $fraction*THETA($real_theta_num)\n";
1986 } else {
1987 push @TX, $comment."IF($cov.EQ.$levels[$lev1])";
1988 $theta_numbering++;
1989 $real_theta_num = $basic_theta_num + $theta_numbering;
1990 push @TX, "$comment$label = -1*THETA($real_theta_num)*(1-$fraction)\n";
1991 push @TX, $comment."IF($cov.EQ.$levels[$lev2])";
1992 $theta_numbering++;
1993 $real_theta_num = $basic_theta_num + $theta_numbering;
1994 push @TX, "$comment$label = $fraction*THETA($real_theta_num)\n";
1995 }
1996
1997 # }}}
1998 } elsif ($self -> {'fix'}[$model_number] and $numlvs > 2) {
1999 # {{{ If at fixing and third level
2000
2001 # Here we have two cases:
2002 # 1) The median is the same as the highest or lowest level
2003 # 2) The median is in the middle
2004 my $med = sprintf "%6.2f",$self -> {'median'}{$cov};
2005 $med =~ s/\s*//;
2006 # Second level
2007 if ($med == $levels[0] or $med == $levels[$#levels]) {
2008 if ($data -> have_missing_data('column_head' => $cov)) {
2009 push @TX, $comment."IF($cov.EQ.$self -> {'missing_data_token'}) THEN\n";
2010 push @TX, "$comment $label = 0\n";
2011 push @TX, $comment."ELSE\n";
2012 $theta_numbering++;
2013 $real_theta_num = $basic_theta_num + $theta_numbering;
2014 my $sign = '-';
2015 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
2016 push @TX, "$comment $label = THETA($real_theta_num)*($cov $sign $med)\n";
2017 push @TX, $comment."ENDIF\n";
2018 } else {
2019 $theta_numbering++;
2020 $real_theta_num = $basic_theta_num + $theta_numbering;
2021 my $sign = '-';
2022 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
2023 push @TX, "$comment$label = THETA($real_theta_num)*($cov $sign $med)\n";
2024 }
2025 } else {
2026 $theta_numbering++;
2027 $real_theta_num = $basic_theta_num + $theta_numbering;
2028 my $sign = '-';
2029 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
2030 push @TX, $comment."IF($cov.LE.$med) $label = THETA($real_theta_num)*($cov $sign $med)\n";
2031 $theta_numbering++;
2032 $real_theta_num = $basic_theta_num + $theta_numbering;
2033 my $sign = '-';
2034 if ( $med < 0 ) { $sign = '+'; $med = -$med; }
2035 push @TX, $comment."IF($cov.GT.$med) $label = THETA($real_theta_num)*($cov $sign $med)\n";
2036 if ($data -> have_missing_data('column_head' => $cov)) {
2037 push @TX, $comment."IF($cov.EQ.$self -> {'missing_data_token'}) $label = 0\n";
2038 }
2039 }
2040
2041 # }}}
2042 } else {
2043 # {{{ Otherwise
2044
2045 if ($level == 2) {
2046 my $first = 0;
2047 ## Print out the code. Note that for the missing covariate to work
2048 ## it has to have a value lower than all the other levels.
2049 my $havemiss;
2050 foreach my $lev (@levels) {
2051 if ($first == 0) {
2052 # Register that there are a missing value
2053 if ($lev == $self -> {'missing_data_token'}) {
2054 $havemiss = 1;
2055 next;
2056 }
2057 # OK, there is a missing
2058 if ($havemiss) {
2059 push @TX, $comment."IF($cov.EQ.$lev.OR.$cov.EQ.$self -> {'missing_data_token'}) $label = 0\n";
2060 $havemiss = 0;
2061 } else {
2062 push @TX, $comment."IF($cov.EQ.$lev) $label = 0\n";
2063 }
2064 $first = 1;
2065 } else {
2066 # This if statement is probably not necessary
2067 if ($havemiss) {
2068 $theta_numbering++;
2069 $real_theta_num = $basic_theta_num + $theta_numbering;
2070 push @TX, $comment."IF($cov.EQ.$lev.OR.$cov.EQ.$self -> {'missing_data_token'}) $label = THETA($real_theta_num)\n";
2071 $havemiss = 0;
2072 } else {
2073 $theta_numbering++;
2074 $real_theta_num = $basic_theta_num + $theta_numbering;
2075 push @TX, $comment."IF($cov.EQ.$lev) $label = THETA($real_theta_num)\n";
2076 }
2077 }
2078 }
2079
2080 # }}}
2081 }
2082
2083 # }}}
2084 }
2085 }
2086 }
2087
2088 }
2089
2090 push @TX, ";;; $label".$level." END\n\n";
2091 }
2092
2093 end cov_code_print
2094
2095 # }}}
2096
2097 # {{{ cov_comb_print
2098 start cov_comb_print
2099 {
2100 # Combining of covariates
2101 foreach my $par (@{$self -> {'parameters'}}) {
2102 my $ind = 0;
2103 my $parmind = 0;
2104 my @covmods;
2105
2106 foreach( @{$self -> {'parm_cov_labels'}{$par}} ){
2107 $ind++;
2108 push( @{$covmods[$parmind]}, $_ );
2109 if( $ind == 4 ){
2110 $parmind ++;
2111 $ind = 0;
2112 }
2113 }
2114 foreach my $i (0.. $#covmods) {
2115 if($i == 0) {
2116 push @TX, "$par"."COV=";
2117 } else {
2118 push @TX, "$par"."COV=$par"."COV*";
2119 }
2120 my @tmp = @{$covmods[$i]};
2121 if($self -> {'exp'}) {
2122 @tmp = map( $_ = "EXP($_)", @tmp);
2123 push( @TX, join('*', @tmp) );
2124 } else {
2125 @tmp = map( $_ = "(1+$_)", @tmp);
2126 push( @TX, join('*', @tmp) );
2127 }
2128 push @TX, "\n";
2129 }
2130 push @TX, "\n\n";
2131 }
2132 }
2133 end cov_comb_print
2134 # }}}
2135
2136 # {{{ _write
2137
2138 start _write
2139 {
2140 my @TX;
2141 push @TX, "\n";
2142 # {{{ Print declarations
2143 {
2144 my $ind = 0;
2145 my $tmp;
2146 foreach( keys %{$self -> {'labels'}}) {
2147 $tmp .= " $_";
2148 $ind ++;
2149 if( $ind == 4 ) {
2150 push @TX, ';;; DECL:'.$tmp,"\n";
2151 $ind = 0;
2152 $tmp = '';
2153 }
2154 }
2155 if( $ind > 0 ){
2156 push @TX, ';;; DECL:'. $tmp . "\n";
2157 }
2158 }
2159 # }}}
2160 # Print Bounds
2161 foreach( keys %{$self -> {'labels'}} ){
2162 my ($lower, $upper) = @{$self -> {'bounds'}{$_}};
2163 push @TX, ';;; BOUNDS: '."($lower,,$upper) ; $_\n";
2164 }
2165 # {{{ Print the extra information before the actual models
2166
2167 ## The missing data label
2168 push @TX, ';;; MISS: '. $self -> {'missing_data_token'} ."\n";
2169 ## P-values
2170 if( $self -> {'gof_description'} ){
2171 push @TX, ';;; GOF: '.$self->{'gof_description'}."\n";
2172 } else {
2173 push @TX, ';;; GOF: EXTERN\n';
2174 }
2175 #push @TX, ';;; PFORW: '. $self -> {'pforward'} . "\n";
2176 #push @TX, ';;; PBACK: '. $self -> {'pbackward'} . "\n";
2177 ## List file
2178 # push @TX, ';;; LST: '. @{@{$self -> {'models'}}[0] -> outputs}[0] -> filename . "\n";
2179 ## Log file
2180 # push @TX, ';;; LOG: '. $self -> logfile . "_scm\n";
2181 ## Task
2182 push @TX, ';;; TASK: '. $self -> {'search_direction'} . "\n";
2183 ## If the thetas were to be fixed
2184 push @TX, ';;; FIX: '. $self -> {'fix'}[$model_number] . "\n";
2185
2186 # }}}
2187
2188 foreach( keys %{$self -> {'labels'}} ) {
2189 my ($code, undef) = $self -> cov_code_print( model_number => $model_number,
2190 'label' => $_,
2191 'level' => 1,
2192 'print_oldstyle' => 1,
2193 'theta_numbering' => 0);
2194 push( @TX, @{$code} );
2195 ($code, undef) = $self -> cov_code_print( model_number => $model_number,
2196 'label' => $_,
2197 'level' => 2,
2198 'print_oldstyle' => 1,
2199 'theta_numbering' => 0);
2200 push( @TX, @{$code} );
2201 ($code, undef) = $self -> cov_code_print( model_number => $model_number,
2202 'label' => $_,
2203 'level' => 3,
2204 'print_oldstyle' => 1,
2205 'theta_numbering' => 0);
2206 push( @TX, @{$code} );
2207 }
2208 my $outmodel = @{$self -> {'models'}}[0] -> copy(filename => $filename,
2209 target => 'mem',
2210 ignore_missing_files => 1);
2211 push( @TX, @{$outmodel -> pk()} );
2212 $outmodel -> pk( 'new_pk' => \@TX );
2213 $outmodel -> write;
2214 }
2215 end _write
2216
2217 # }}}
2218
2219 # {{{ _create_models
2220
2221 start _create_models
2222 {
2223 # Set names of tested parameter-covariate-state combinations in results
2224 # Done in loop below after state change
2225 my %return_section;
2226 $return_section{'name'} = 'combinations';
2227 $return_section{'values'} = [];
2228 $return_section{'labels'} = [];
2229
2230 my $done = ( -e $self -> {'directory'}."/m$model_number/done" ) ? 1 : 0;
2231
2232 if ( not $done ) {
2233 open( DONE_LOG, '>'.$self -> {'directory'}."/m$model_number/done.log" );
2234 foreach my $parameter ( sort keys %relations ) {
2235 foreach my $covariate ( sort keys %{$relations{$parameter}} ) {
2236 # Is this covariate continuous or not?
2237 my $continuous = 1;
2238 foreach my $cat ( @{$self -> {'categorical_covariates'}} ) {
2239 $continuous = 0 if ( $covariate eq $cat );
2240 }
2241 my @valid_states;
2242 if ( $continuous ) {
2243 @valid_states = @{$self -> {'valid_states'}{'continuous'}};
2244 } else {
2245 @valid_states = @{$self -> {'valid_states'}{'categorical'}};
2246 }
2247 my $state;
2248 # Important: just setting $state to $self->incl_rel....{'state'} initiates
2249 # included_relations for this parameter and covariate. Avoid this.
2250 if ( defined $included_relations{$parameter}{$covariate} ) {
2251 $state = $included_relations{$parameter}{$covariate}{'state'};
2252 }
2253 $state = defined $state ? $state : $valid_states[0];
2254 #1. Create a new model if the state is not yet included and at highest level
2255 # (forward search) or base level (backward search).
2256 next if ( ( $self -> {'search_direction'} eq 'forward' and
2257 $state == $valid_states[$#valid_states] ) or
2258 ( $self -> {'search_direction'} eq 'backward' and
2259 $state == $valid_states[0] ) );
2260 my $old_state = $state;
2261 # Increment (forward search) or decrement (backward search) the state
2262 if ( $self -> {'search_direction'} eq 'forward' ) {
2263 my $flag = 0;
2264 for( my $s_idx = 0; $s_idx <= $#valid_states; $s_idx++ ) {
2265 if ( $flag ) {
2266 $state = $valid_states[$s_idx];
2267 last;
2268 }
2269 $flag = 1 if( $state == $valid_states[$s_idx] );
2270 }
2271 } elsif ( $self -> {'search_direction'} eq 'backward' ) {
2272 my $flag = 0;
2273 for( my $s_idx = $#valid_states; $s_idx >= 0; $s_idx-- ) {
2274 if ( $flag ) {
2275 $state = $valid_states[$s_idx];
2276 last;
2277 }
2278 $flag = 1 if( $state == $valid_states[$s_idx] );
2279 }
2280 }
2281
2282 # Only one problem and one sub problem
2283 push( @{$return_section{'values'}[0][0]}, "$parameter$covariate-$state" );
2284
2285 my ( $dir, $filename ) =
2286 OSspecific::absolute_path( $self -> {'directory'}.
2287 '/m'.$model_number.'/',
2288 $parameter.$covariate.$state.".mod" );
2289 my ( $odir, $outfilename ) =
2290 OSspecific::absolute_path( $self -> {'directory'}.
2291 '/m'.$model_number.'/',
2292 $parameter.$covariate.$state.".lst" );
2293 my $applicant_model;
2294 $applicant_model = $orig_model ->
2295 copy( filename => $dir.$filename,
2296 copy_data => 0,
2297 copy_output => 0 );
2298
2299 $applicant_model -> ignore_missing_files(1);
2300 $applicant_model -> outputfile( $odir.$outfilename );
2301 my @table_names = @{$applicant_model -> table_names};
2302 for ( my $i = 0; $i <= $#table_names; $i++ ) {
2303 for ( my $j = 0; $j < scalar @{$table_names[$i]}; $j++ ) {
2304 $table_names[$i][$j] = $self -> {'directory'}.
2305 '/m'.$model_number.'/'.
2306 $filename.'.'.
2307 OSspecific::nopath($table_names[$i][$j]);
2308 }
2309 }
2310 $applicant_model -> table_names( new_names => \@table_names,
2311 ignore_missing_files => 1);
2312 # print "$parameter-$covariate: new state: $state\n";
2313 my @used_covariates = ();
2314 # $included_relations is a reference to $self -> included_relations
2315 # and should only be initialized for truly incorporated relations
2316 # see beginning of loop above.
2317 foreach my $incl_par ( sort keys %included_relations ) {
2318 foreach my $incl_cov ( sort keys %{$included_relations{$incl_par}} ) {
2319 next if ( $incl_par eq $parameter and $incl_cov eq $covariate );
2320 $self ->
2321 add_code( definition_code => $included_relations{$incl_par}{$incl_cov}{'code'},
2322 nthetas => $included_relations{$incl_par}{$incl_cov}{'nthetas'},
2323 inits => $included_relations{$incl_par}{$incl_cov}{'inits'},
2324 bounds => $included_relations{$incl_par}{$incl_cov}{'bounds'},
2325 applicant_model => $applicant_model,
2326 parameter => $incl_par,
2327 covariate => $incl_cov );
2328 push( @used_covariates, $incl_cov );
2329 }
2330 }
2331 # If the new state is base level (backward search) don't add this relation, otherwise:
2332 unless ( $state == $valid_states[0] ) {
2333 $self -> add_code( definition_code => $relations{$parameter}{$covariate}{'code'}{$state},
2334 nthetas => $relations{$parameter}{$covariate}{'nthetas'}{$state},
2335 inits => $relations{$parameter}{$covariate}{'inits'}{$state},
2336 bounds => $relations{$parameter}{$covariate}{'bounds'}{$state},
2337 applicant_model => $applicant_model,
2338 parameter => $parameter,
2339 covariate => $covariate );
2340 push( @used_covariates, $covariate );
2341 }
2342 $self -> drop_undrop_covariates( applicant_model => $applicant_model,
2343 used_covariates => \@used_covariates,
2344 all_covariates => [@{$self -> {'categorical_covariates'}},
2345 @{$self -> {'continuous_covariates'}}],
2346 do_not_drop => $self -> {'do_not_drop'});
2347
2348 if ( defined $base_model ) {
2349 $applicant_model -> update_inits( from_model => $base_model,
2350 update_thetas => 1,
2351 update_omegas => 1,
2352 update_sigmas => 1,
2353 ignore_missing_parameters => 1 );
2354 } else {
2355 $applicant_model -> update_inits( from_model => $orig_model,
2356 update_thetas => 1,
2357 update_omegas => 1,
2358 update_sigmas => 1,
2359 ignore_missing_parameters => 1 );
2360 }
2361 $applicant_model -> _write;
2362 push( @new_models, $applicant_model );
2363 my $model_id = $applicant_model -> register_in_database;
2364 $self -> register_tm_relation( model_ids => [$model_id],
2365 prepared_models => 1 );
2366 my %st_rel;
2367 $st_rel{'parameter'} = $parameter;
2368 $st_rel{'covariate'} = $covariate;
2369 $st_rel{'state'} = $state;
2370 $st_rel{'continuous'} = $continuous;
2371 push( @step_relations, \%st_rel );
2372 print DONE_LOG "$parameter $covariate $continuous $old_state $state\n";
2373 }
2374 }
2375 open( TMP, ">".$self -> {'directory'}."/m$model_number/done" );
2376 close( TMP );
2377 close( DONE_LOG );
2378 } else {
2379 ui -> print( category => 'scm',
2380 message => "Recreating models from previously run step" );
2381 debug -> warn( level => 2,
2382 message => "Creating applicant model from file on disk" );
2383 if ( not -e $self -> {'directory'}."/m$model_number/done.log" ) {
2384 debug -> die( message => "No file ".$self -> {'directory'}.
2385 "/m$model_number/done.log seem to exist although the existance".
2386 " of the file ".$self -> {'directory'}.
2387 "/m$model_number/done indicates so.");
2388 }
2389 open( DONE_LOG, $self -> {'directory'}."/m$model_number/done.log" );
2390 my @rows = <DONE_LOG>;
2391 close( DONE_LOG );
2392 for( my $i = 0; $i <= $#rows; $i++ ) { # skip first row
2393 chomp( $rows[$i] );
2394 my ( $parameter, $covariate, $continuous, $old_state, $state ) =
2395 split(' ',$rows[$i],5);
2396 my @valid_states;
2397 if ( $continuous ) {
2398 @valid_states = @{$self -> {'valid_states'}{'continuous'}};
2399 } else {
2400 @valid_states = @{$self -> {'valid_states'}{'categorical'}};
2401 }
2402 #1. Recreate the model if the state is not yet included and at highest level
2403 # (forward search) or base level (backward search).
2404 next if ( ( $self -> {'search_direction'} eq 'forward' and
2405 $old_state == $valid_states[$#valid_states] ) or
2406 ( $self -> {'search_direction'} eq 'backward' and
2407 $old_state == $valid_states[0] ) );
2408
2409 my ( $dir, $filename ) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.$model_number.'/',$parameter.$covariate.$state.".mod");
2410 my ( $odir, $outfilename ) = OSspecific::absolute_path( $self -> {'directory'}.'/m'.$model_number.'/',$parameter.$covariate.$state.".lst");
2411
2412 my $applicant_model = model -> new( reference_object => $orig_model,
2413 outputs => undef,
2414 datas => undef,
2415 synced => undef,
2416 problems => undef,
2417 active_problems => undef,
2418 filename => $dir.$filename,
2419 outputfile => $odir.$outfilename,
2420 target => 'disk',
2421 ignore_missing_files => 1 );
2422 # Set the correct data file for the object
2423 my $moddir = $orig_model -> directory;
2424 my @datafiles = @{$orig_model -> datafiles};
2425 for( my $df = 0; $df <= $#datafiles; $df++ ) {
2426 $datafiles[$df] = $moddir.'/'.$datafiles[$df];
2427 }
2428 $applicant_model -> synchronize;
2429 $applicant_model -> datafiles( new_names => \@datafiles );
2430 push( @new_models, $applicant_model );
2431 my %st_rel;
2432 $st_rel{'parameter'} = $parameter;
2433 $st_rel{'covariate'} = $covariate;
2434 $st_rel{'state'} = $state;
2435 $st_rel{'continuous'} = $continuous;
2436 push( @step_relations, \%st_rel );
2437 my $nl = $i == $#rows ? "" : "\r";
2438 ui -> print( category => 'scm',
2439 message => ui -> status_bar( sofar => $i+1,
2440 goal => $#rows+1 ).$nl,
2441 wrap => 0,
2442 newline => 0 )
2443 }
2444 ui -> print( category => 'scm',
2445 message => " ... done." );
2446 }
2447 push( @{$self -> {'results'}[$model_number-1]{'own'}},\%return_section );
2448 }
2449 end _create_models
2450
2451 # }}} _create_models
2452
2453 # {{{ calculate_categorical_statistics
2454
2455 start calculate_categorical_statistics
2456 {
2457 # Levels = the number of unique factors of a covariate
2458 # my @levels = sort {$a<=>$b} keys %factors;
2459 # Sort by frequency
2460 my @sorted = sort {$factors{$a}<=>$factors{$b}} keys (%factors);
2461 # These lines will set the most common value in $medians{$cov}
2462 $median = $sorted[0]; # First element of the sorted array
2463 # (the factor that most subjects have)
2464
2465 # $min = $sorted[$#sorted]; # The first element of the (also sorted) level array
2466 # $max = $sorted[0]; # The last element of the (also sorted) level array
2467
2468 }
2469 end calculate_categorical_statistics
2470
2471 # }}} calculate_categorical_statistics
2472
2473 # {{{ calculate_continuous_statistics
2474
2475 start calculate_continuous_statistics
2476 {
2477 # Max, median and min
2478 $median = $model -> median( problem_number => 1,
2479 'column_head' => $covariate,
2480 'unique_in_individual' => 1);
2481 $max = $model -> max( problem_number => 1,
2482 'column_head' => $covariate );
2483 $min = $model -> min( problem_number => 1,
2484 'column_head' => $covariate );
2485
2486 $median = sprintf("%6.2f", $median );
2487 }
2488 end calculate_continuous_statistics
2489
2490 # }}} calculate_continuous_statistics
2491
2492 # {{{ create_linear_code
2493
2494 start create_linear_code
2495 {
2496 # Creates code for linear continuous and categorical
2497 # parameter-covariate relations Returns code, the number of
2498 # thetas necessary for the relation plus initial values and
2499 # boundaries for these.
2500 if ( defined $bounds{'upper'} and ref $bounds{'upper'} ne 'ARRAY' ) {
2501 debug -> die( message => "Wrong data type of upper bounds. The upper bounds should be kept in an array" );
2502 }
2503 if ( defined $bounds{'lower'} and ref $bounds{'lower'} ne 'ARRAY' ) {
2504 debug -> die( message => "Wrong data type of lower bounds. The lower bounds should be kept in an array" );
2505 }
2506
2507 my $theta_number = defined $start_theta ? $start_theta : 1;
2508 my $comment = "";
2509 my $median = $statistics{'median'};
2510 my $min = $statistics{'min'};
2511 my $max = $statistics{'max'};
2512 $median = sprintf "%6.2f", $median;
2513 $median =~ s/\s*//;
2514 if ( $continuous ) {
2515 if ( scalar @code < 1 or ( scalar @code == 1 and $code[0] eq '' ) ) {
2516 if( $statistics{'have_missing_data'} ) {
2517 $code[0] = $comment."IF($covariate.EQ.$missing_data_token) THEN\n";
2518 $code[1] = "$comment $parameter$covariate = 0\n";
2519 $code[2] = $comment."ELSE\n";
2520 my $sign = '-';
2521 if ( $median < 0 ) { $sign = '+'; $median = -$median; }
2522 $code[3] = "$comment $parameter$covariate = THETA(".$theta_number++.")*($covariate $sign $median)\n";
2523 $code[4] = $comment."ENDIF\n";
2524 } else {
2525 my $sign = '-';
2526 if ( $median < 0 ) { $sign = '+'; $median = -$median; }
2527 $code[0] = "$comment$parameter$covariate = THETA(".$theta_number++.")*($covariate $sign $median)\n";
2528 }
2529 } else {
2530 # count the thetas. One theta per row is allowed
2531 for ( @code ) {
2532 $theta_number++ if /THETA\((\d+)\)/;
2533 }
2534 }
2535 # Boundaries
2536 if ( $median-$min == 0 ) {
2537 debug -> die( message => "median($median)-min($min)",
2538 " for column $covariate is zero. Cannot calculate an upper boundary." );
2539 }
2540 if ( $median-$max == 0 ) {
2541 debug -> die( message => "median($median)-min($min)",
2542 " for column $covariate is zero. Cannot calculate an lower boundary." );
2543 }
2544
2545 ## Boundaries:
2546 unless ( defined $bounds{'upper'} and defined $bounds{'upper'}[0] ) {
2547 my $upper_bound = 1/($median-$min);
2548 my ($big,$small) = split('\.',$upper_bound);
2549 $small = substr($small,0,3);
2550 $upper_bound = $big.'.'.$small;
2551 $upper_bound = '0' if eval($upper_bound) == 0;
2552 $bounds{'upper'}[0] = $upper_bound;
2553 }
2554 unless ( defined $bounds{'lower'} and defined $bounds{'lower'}[0] ) {
2555 my $lower_bound = 1/($median - $max);
2556 my ($big,$small) = split('\.',$lower_bound);
2557 $small = substr($small,0,3);
2558 $lower_bound = $big.'.'.$small;
2559 $lower_bound = '0' if eval($lower_bound) == 0;
2560 $bounds{'lower'}[0] = $lower_bound;
2561 }
2562
2563 } else {
2564 my %factors = %{$statistics{'factors'}};
2565 my @levels = sort {$a<=>$b} keys %factors;
2566 my @sorted = sort {$factors{$b}<=>$factors{$a}} keys(%factors);
2567 my $numlvs = scalar @levels;
2568 $numlvs = $numlvs -1 if $statistics{'have_missing_data'};
2569
2570 my $first_non_missing = 1;
2571 if ( scalar @code < 1 or ( scalar @code == 1 and $code[0] eq '' ) ) {
2572 @code = ();
2573 for ( my $i = 0; $i <= $#sorted; $i++ ) {
2574 if ( $statistics{'have_missing_data'} and
2575 ( $self -> {'missing_data_token'} eq $sorted[$i] ) ) {
2576 push @code, $comment."IF($covariate.EQ.$sorted[$i]) $parameter$covariate".
2577 " = 0 ; Missing data\n";
2578 } else {
2579 if ( $first_non_missing ) {
2580 push @code, $comment."IF($covariate.EQ.$sorted[$i]) $parameter$covariate".
2581 " = 0 ; Most common\n";
2582 $first_non_missing = 0;
2583 } else {
2584 push @code, $comment."IF($covariate.EQ.$sorted[$i]) $parameter$covariate".
2585 " = THETA(".$theta_number++.")\n";
2586 }
2587 }
2588 }
2589 } else {
2590 # count the thetas. One theta per row is allowed
2591 for ( @code ) {
2592 $theta_number++ if /THETA\((\d+)\)/;
2593 }
2594 }
2595 ## Boundaries:
2596 for ( my $i = 0; $i < $theta_number - $start_theta; $i++ ) {
2597 $bounds{'upper'}[$i] = 5 unless ( defined $bounds{'upper'}[$i] );
2598 $bounds{'lower'}[$i] = -1 unless ( defined $bounds{'lower'}[$i] );
2599 }
2600 }
2601 # Initial values
2602 for ( my $i = 0; $i < $theta_number - $start_theta; $i++ ) {
2603 next if ( defined $inits[$i] );
2604 my $fraction = $self -> {'global_init'};
2605 my $tmp;
2606 if( ( abs($bounds{'upper'}[$i]) >= 1000000 or
2607 not defined $bounds{'upper'}[$i] ) and
2608 ( abs($bounds{'lower'}[$i]) >= 1000000 or
2609 not defined $bounds{'lower'}[$i] ) ) {
2610 $tmp = $fraction*100;
2611 } else {
2612 if ( abs($bounds{'upper'}[$i]) < abs($bounds{'lower'}[$i]) ) {
2613 $tmp = $bounds{'upper'}[$i] == 0 ? $bounds{'lower'}[$i]*$fraction :
2614 $bounds{'upper'}[$i]*$fraction;
2615 } else {
2616 $tmp = $bounds{'lower'}[$i] == 0 ? $bounds{'upper'}[$i]*$fraction :
2617 $bounds{'lower'}[$i]*$fraction;
2618 }
2619 }
2620 $inits[$i] = $tmp;
2621 }
2622 $end_theta = --$theta_number;
2623 }
2624 end create_linear_code
2625
2626 # }}} create_linear_code
2627
2628 # {{{ create_exponential_code
2629
2630 start create_exponential_code
2631 {
2632 my $theta_number = $start_theta;
2633 my $median = $statistics{'median'};
2634 $median = sprintf "%6.2f", $median;
2635 $median =~ s/\s*//;
2636 if ( $continuous ) {
2637 if ( scalar @code < 1 or ( scalar @code == 1 and $code[0] eq '' ) ) {
2638 my $sign = '-';
2639 if ( $median < 0 ) { $sign = '+'; $median = -$median; }
2640 $code[0] = " $parameter$covariate = EXP(THETA(".$theta_number++.")*($covariate $sign $median))\n";
2641 if ( $statistics{'have_missing_data'} ) {
2642 $code[1] = "IF($covariate.EQ.$missing_data_token) ".
2643 "$parameter$covariate = 0\n";
2644 }
2645 } else {
2646 # count the thetas. One theta per row is allowed
2647 for ( @code ) {
2648 $theta_number++ if /THETA\((\d+)\)/;
2649 }
2650 }
2651 ## Boundaries:
2652 unless ( defined $bounds{'upper'} and defined $bounds{'upper'}[0] ) {
2653 $bounds{'upper'}[0] = 1000000;
2654 }
2655 unless ( defined $bounds{'lower'} and defined $bounds{'lower'}[0] ) {
2656 $bounds{'lower'}[0] = -1000000;
2657 }
2658 }
2659 # Initial values
2660 for ( my $i = 0; $i < $theta_number - $start_theta; $i++ ) {
2661 next if ( defined $inits[$i] );
2662 $inits[$i] = $self -> {'global_init'};
2663 }
2664 $end_theta = --$theta_number;
2665 }
2666 end create_exponential_code
2667
2668 # }}} create_exponential_code
2669
2670 # {{{ create_hockey_stick_code
2671
2672 start create_hockey_stick_code
2673 {
2674 my $theta_number = $start_theta;
2675 my $median = $statistics{'median'};
2676 $median = sprintf "%6.2f", $median;
2677 my $min = $statistics{'min'};
2678 my $max = $statistics{'max'};
2679 $median =~ s/\s*//;
2680 if ( $continuous ) {
2681 if ( scalar @code < 1 or ( scalar @code == 1 and $code[0] eq '' ) ) {
2682 my $sign = '-';
2683 if ( $median < 0 ) { $sign = '+'; $median = -$median; }
2684 $code[0] = "IF($covariate.LE.$median) $parameter$covariate = THETA(".
2685 $theta_number++.")*($covariate $sign $median)\n";
2686 $code[1] = "IF($covariate.GT.$median) $parameter$covariate ".
2687 "= THETA(".$theta_number++.")*($covariate $sign $median)\n";
2688 if ( $statistics{'have_missing_data'} ) {
2689 $code[2] = "IF($covariate.EQ.$missing_data_token) ".
2690 "$parameter$covariate = 0\n";
2691 }
2692 } else {
2693 # count the thetas. One theta per row is allowed
2694 for ( @code ) {
2695 $theta_number++ if /THETA\((\d+)\)/;
2696 }
2697 }
2698 ## Boundaries:
2699 if ( not defined $bounds{'upper'} ) {
2700 if ( not defined $bounds{'upper'}[0] ) {
2701 my $upper_bound = 1/($median-$min);
2702 my ($big,$small) = split('\.',$upper_bound);
2703 $small = substr($small,0,3);
2704 $upper_bound = $big.'.'.$small;
2705 $upper_bound = '0' if eval($upper_bound) == 0;
2706 $bounds{'upper'}[0] = $upper_bound;
2707 }
2708 $bounds{'upper'}[1] = 1000000 if ( not defined $bounds{'upper'}[1] );
2709 }
2710 if ( not defined $bounds{'lower'} ) {
2711 $bounds{'lower'}[0] = -1000000 if ( not defined $bounds{'lower'}[0] );
2712 if ( not defined $bounds{'lower'}[1] ) {
2713 my $lower_bound = 1/($median - $max);
2714 my ($big,$small) = split('\.',$lower_bound);
2715 $small = substr($small,0,3);
2716 $lower_bound = $big.'.'.$small;
2717 $lower_bound = '0' if eval($lower_bound) == 0;
2718 $bounds{'lower'}[1] = $lower_bound;
2719 }
2720 }
2721 }
2722 # Initial values
2723 for ( my $i = 0; $i < $theta_number - $start_theta; $i++ ) {
2724 next if ( defined $inits[$i] );
2725 my $fraction = $self -> {'global_init'};
2726 my $tmp;
2727 if( ( abs($bounds{'upper'}[$i]) >= 1000000 or
2728 not defined $bounds{'upper'}[$i] ) and
2729 ( abs($bounds{'lower'}[$i]) >= 1000000 or
2730 not defined $bounds{'lower'}[$i] ) ) {
2731 $tmp = $fraction*100;
2732 } else {
2733 if ( abs($bounds{'upper'}[$i]) < abs($bounds{'lower'}[$i]) ) {
2734 $tmp = $bounds{'upper'}[$i] == 0 ? $bounds{'lower'}[$i]*$fraction :
2735 $bounds{'upper'}[$i]*$fraction;
2736 } else {
2737 $tmp = $bounds{'lower'}[$i] == 0 ? $bounds{'upper'}[$i]*$fraction :
2738 $bounds{'lower'}[$i]*$fraction;
2739 }
2740 }
2741 $inits[$i] = $tmp;
2742 }
2743 $end_theta = --$theta_number;
2744 }
2745 end create_hockey_stick_code
2746
2747 # }}} create_hockey_stick_code
2748
2749 # {{{ create_code
2750
2751 start create_code
2752 {
2753 my $new_code = [];
2754 my $new_inits = [];
2755 my $new_bounds = {};
2756 if ( $state == 1 ) {
2757 # This is the default basic state. No relation included
2758 } elsif ( $state == 2 ) {
2759 # First state
2760 ( $new_code, $end_theta, $new_inits, $new_bounds ) =
2761 $self -> create_linear_code( start_theta => $start_theta,
2762 parameter => $parameter,
2763 covariate => $covariate,
2764 continuous => $continuous,
2765 bounds => \%bounds,
2766 inits => \@inits,
2767 code => \@code,
2768 statistics => \%statistics,
2769 missing_data_token => $missing_data_token,
2770 fix => $fix );
2771 } elsif ( $state == 3 ) {
2772 ( $new_code, $end_theta, $new_inits, $new_bounds ) =
2773 $self -> create_hockey_stick_code( start_theta => $start_theta,
2774 parameter => $parameter,
2775 covariate => $covariate,
2776 continuous => $continuous,
2777 bounds => \%bounds,
2778 inits => \@inits,
2779 code => \@code,
2780 statistics => \%statistics,
2781 missing_data_token => $missing_data_token,
2782 fix => $fix );
2783 } elsif ( $state == 4 ) {
2784 ( $new_code, $end_theta, $new_inits, $new_bounds ) =
2785 $self -> create_exponential_code( start_theta => $start_theta,
2786 parameter => $parameter,
2787 covariate => $covariate,
2788 continuous => $continuous,
2789 bounds => \%bounds,
2790 inits => \@inits,
2791 code => \@code,
2792 statistics => \%statistics,
2793 missing_data_token => $missing_data_token,
2794 fix => $fix );
2795 } else {
2796 debug -> die( message => "No such state $state\n" );
2797 }
2798 @inits = @{$new_inits};
2799 %bounds = %{$new_bounds};
2800 @code = @{$new_code};
2801 }
2802 end create_code
2803
2804 # }}} create_code
2805
2806 # {{{ add_code
2807
2808 start add_code
2809 {
2810 my @labels;
2811 for ( my $i = 1; $i <= $nthetas; $i++ ) {
2812 push( @labels, $parameter.$covariate.$i );
2813 }
2814
2815 my $start_theta = $applicant_model -> nthetas + 1;
2816 my $end_theta = $start_theta + $nthetas - 1;
2817
2818 my $tmp = $start_theta;
2819 for ( @definition_code ) {
2820 if ( /THETA\((\d+)\)/ ) {
2821 s/THETA\((\d+)\)/THETA\($tmp\)/;
2822 $tmp++;
2823 }
2824 }
2825
2826 my $start_theta = $applicant_model -> nthetas + 1;
2827
2828 # Add the definition_code to the PK or PRED block
2829 my @code;
2830 @code = @{$applicant_model -> pk( problem_number => 1 )};
2831 my $use_pred = 0;
2832 unless ( $#code > 0 ) {
2833 @code = @{$applicant_model -> pred( problem_number => 1 )};
2834 $use_pred = 1;
2835 }
2836 if ( $#code <= 0 ) {
2837 debug -> die( message => "Neither PK or PRED defined in " .
2838 $applicant_model -> filename . "\n" );
2839 }
2840
2841 my $found_REL = 0;
2842 my $i = 0;
2843 my $relationarea = 0;
2844 my @row;
2845 my $found_correct_REL = 0;
2846 for ( @code ) {
2847 if ( /^;;; (\w+)-RELATION START/ and $1 eq $parameter ) {
2848 $relationarea = 1;
2849 $i++;
2850 next;
2851 }
2852 if ( /^;;; (\w+)-RELATION END/ and $1 eq $parameter ) {
2853 $relationarea = 0;
2854 last;
2855 }
2856 if ( $relationarea ) {
2857 $found_REL = $i;
2858 debug -> warn( level => 2,
2859 message => $parameter . "COV has already been added to the code" );
2860 $found_correct_REL = 1 and last if ( /$parameter$covariate/ );
2861 @row = split(/\)\*\(/);
2862 }
2863 $i++;
2864 }
2865
2866 # If we have old scm code present.
2867 if ( $found_REL ) {
2868 unless ( $found_correct_REL ) {
2869 if ( $#row > 2 ) {
2870 @code = (@code[0..$found_REL],
2871 "$parameter"."COV=$parameter"."COV*(1+$parameter$covariate)\n",
2872 @code[$found_REL+1..$#code]);
2873 } else {
2874 chomp($code[$found_REL]);
2875 $code[$found_REL] = $code[$found_REL]."*(1+$parameter$covariate)\n";
2876 }
2877 }
2878 } else {
2879 @code = ( ";;; $parameter-RELATION START\n",
2880 "$parameter"."COV=(1+$parameter$covariate)\n",
2881 ";;; $parameter-RELATION END\n\n",
2882 @code );
2883 }
2884
2885 @code = ( "\n;;; $parameter$covariate-DEFINITION START\n",
2886 @definition_code,
2887 ";;; $parameter$covariate-DEFINITION END\n\n",
2888 @code );
2889
2890 # Add to the parameter code
2891 unless ( $found_REL ) {
2892 my $success = 0;
2893 for ( @code ) {
2894 if ( /^\s*TV(\w+)\s*=\s*/ and $1 eq $parameter ) {
2895 my ($line,$comment) = split( ';', $_, 2 );
2896 $_ = $line;
2897 # if ( /^\s*TV(\w+)\s*=\s*THETA\((\d+)\)/ and $1 eq $parameter ) {
2898 # s/^(\s*)TV(\w+)\s*=\s*THETA\((\d+)\)/$1TV$2 = $2COV*THETA($3)/;
2899 chomp;
2900 s/^(\s*)TV(\w+)\s*=\s*(.*)/$1TV$2 = $2COV*\($3\)/;
2901 $success = 1;
2902 $_ = $_.' ;'.$comment."\n";
2903 }
2904 }
2905 unless ( $success ) {
2906 debug -> die( message => "Could not determine a good place to add the covariate relation.\n".
2907 " i.e. No TV$parameter was found\n" );
2908 }
2909 }
2910 if ( $use_pred ) {
2911 $applicant_model -> pred( problem_number => 1,
2912 new_pred => \@code );
2913 } else {
2914 $applicant_model -> pk( problem_number => 1,
2915 new_pk => \@code );
2916 }
2917
2918
2919 $applicant_model -> initial_values( parameter_numbers => [[$start_theta..$end_theta]],
2920 new_values => [\@inits],
2921 add_if_absent => 1,
2922 parameter_type => 'theta',
2923 problem_numbers => [1]);
2924 $applicant_model -> upper_bounds( parameter_type => 'theta',
2925 parameter_numbers => [[$start_theta..$end_theta]],
2926 problem_numbers => [1],
2927 new_values => [$bounds{'upper'}] );
2928 $applicant_model -> lower_bounds( parameter_type => 'theta',
2929 parameter_numbers => [[$start_theta..$end_theta]],
2930 problem_numbers => [1],
2931 new_values => [$bounds{'lower'}] );
2932 $applicant_model -> labels( parameter_type => 'theta',
2933 parameter_numbers => [[$start_theta..$end_theta]],
2934 problem_numbers => [1],
2935 new_values => [\@labels] );
2936 }
2937 end add_code
2938
2939 # }}} add_code
2940
2941 # {{{ write_log
2942
2943 start write_log
2944 {
2945 my %parm = @_;
2946 my $direction = $parm{'direction'};
2947 my $logfile = $parm{'logfile'};
2948 my $included_relations = $parm{'included_relations'};
2949 my $chosen_parameter = $parm{'chosen_parameter'};
2950 my $chosen_covariate = $parm{'chosen_covariate'};
2951 my $results = $parm{'results'};
2952 my $criterion = $parm{'criterion'};
2953 my $test_log_ref = $parm{'test_log'};
2954 my $append = $parm{'append'};
2955 my %test_log = %{$test_log_ref} if defined ( $test_log_ref );
2956 my @names = ();
2957 my @drops = ();
2958 my $chosen;
2959 foreach my $result ( @{$results} ) {
2960 @names = @{$result -> {'values'}} if ($result -> {'name'} eq 'combinations');
2961 @drops = @{$result -> {'values'}} if ($result -> {'name'} eq 'ofv.drop');
2962 $chosen = $result -> {'values'} if ($result -> {'name'} eq 'relation.chosen.in.step');
2963 }
2964
2965 open( LOG, ">>$logfile" );
2966 # for( my $i = 0; $i <= $#names; $i++ ) {
2967 # print LOG 'DROP ',sprintf("%10s",$names[$i]),sprintf("%10.1f",$drops[$i]),"\n";
2968 # }
2969
2970 if ( defined $chosen ) {
2971 print LOG "Parameter-covariate relation chosen in this $direction step: ",
2972 "$chosen_parameter-$chosen_covariate\n";
2973 print LOG sprintf("%-15s",'CRITERION'),uc( $criterion ),"\n";
2974 my @names = sort keys %test_log;
2975 foreach my $name ( @names ) {
2976 my $val = $test_log{$name};
2977 # while ( my ( $name, $val ) = sort {$a <=> $b} each %test_log ) {
2978 if ( ref($val) eq 'HASH' ) {
2979 foreach my $name2 ( sort keys %{$val} ) {
2980 # while ( my ( $name2, $val2 ) = each %{$val} ) {
2981 print LOG sprintf("%-20s",uc($name)),sprintf("%-30s",uc( $name2) ),
2982 sprintf("%12.5f",uc( $val -> {$name2} )),"\n";
2983 }
2984 } else {
2985 print LOG sprintf("%-20s",uc($name)),sprintf("%12.5f",uc( $val )),"\n";
2986 }
2987 }
2988 }
2989 if (defined $included_relations and scalar %{$included_relations} > 0 ) {
2990 print LOG "Relations included after this step:\n";
2991 foreach my $kind ( 'categorical', 'continuous' ) {
2992 print LOG sprintf("%-15s",'KIND'),uc( $kind ),"\n";
2993 print LOG sprintf("%-15s",'VALID_STATES');
2994 foreach my $state ( @{$self -> {'valid_states'}{$kind}} ) {
2995 print LOG sprintf("%-4s",$state)
2996 }
2997 print LOG "\n";
2998 print LOG sprintf("%-15s",'COVARIATES');
2999 foreach my $cov ( @{$self -> {lc($kind).'_covariates'}} ) {
3000 print LOG sprintf("%-8s",$cov)
3001 }
3002 print LOG "\n";
3003 foreach my $par ( sort keys %{$included_relations} ) {
3004 print LOG sprintf("%-15s",'INCL'),sprintf("%-8s",$par);
3005 foreach my $cov ( sort keys %{$included_relations -> {$par}} ) {
3006 # Is this covariate continuous or not?
3007 my $continuous = 1;
3008 foreach my $cat ( @{$self -> {'categorical_covariates'}} ) {
3009 $continuous = 0 if ( $cov eq $cat );
3010 }
3011 if ( ( $continuous and $kind eq 'continuous' ) or
3012 not $continuous and $kind eq 'categorical' ) {
3013 print LOG sprintf("%-8s",$cov.'-'.$included_relations -> {$par}{$cov}{'state'});
3014 }
3015 }
3016 print LOG "\n";
3017 }
3018 }
3019 }
3020 print LOG "--------------------\n\n";
3021 close( LOG );
3022 }
3023 end write_log
3024
3025 # }}}
3026
3027 # This method is not used at the moment and a method with the same name exists in tool.pm
3028 # Keep the code for now but rename the method if it is to be used again.
3029 # # {{{ read_log
3030
3031 # start read_log
3032 # {
3033 # my @log;
3034 # open( LOG, $self -> {'logfile'} );
3035 # @log = <LOG>;
3036 # close( LOG );
3037 # my @step_indexes;
3038 # for ( my $i = 0; $i <= $#log; $i++ ) {
3039 # push( @step_indexes, $i ) if ( $log[$i] =~ /--------------------/ );
3040 # }
3041 # # my %incl_relations;
3042 # # my %valid_states;
3043 # my ( $junk, $kind );
3044 # for ( my $i = $step_indexes[$#step_indexes-1]; $i <= $step_indexes[$#step_indexes]; $i++ ) {
3045 # if ( $log[$i] =~ /^KIND/ ) {
3046 # ( $junk, $kind ) = split( ' ', $log[$i] );
3047 # }
3048 # if ( $log[$i] =~ /^VALID_STATES/ ) {
3049 # my @valid_states;
3050 # @valid_states = split( ' ', $log[$i] );
3051 # shift @valid_states;
3052 # $valid_states{lc($kind)} = \@valid_states;
3053 # }
3054 # if ( $log[$i] =~ /^INCL/ ) {
3055 # my @row = split( ' ', $log[$i] );
3056 # my $par = $row[1];
3057 # my @cov_states = @row[2..$#row];
3058 # foreach my $covstate ( @cov_states ) {
3059 # my ( $cov, $state ) = split('-',$covstate);
3060 # if ( not defined $state or $state eq '' ) {
3061 # $state = $valid_states{$kind}[0];
3062 # }
3063 # $incl_relations{$par}{$cov}{'state'} = $state;
3064 # }
3065 # }
3066 # }
3067 # print;
3068 # }
3069 # end read_log
3070
3071 # # }}}
3072
3073 # {{{ drop_undrop_covariates
3074
3075 start drop_undrop_covariates
3076 {
3077 OUTER: foreach my $cov ( @all_covariates ) {
3078 my $used = 0;
3079 foreach my $do_not_cov ( @do_not_drop ) {
3080 next OUTER if ( $cov eq $do_not_cov );
3081 }
3082 foreach my $used_cov ( @used_covariates ) {
3083 $used = 1 if ( $cov eq $used_cov );
3084 }
3085 if ( $used ) {
3086 debug -> warn( level => 2,
3087 message => "undropping $cov" );
3088 $applicant_model -> _option_val_pos ( problem_numbers => [1],
3089 instance_numbers => [[1]],
3090 name => $cov,
3091 record_name => 'input',
3092 new_values => [['']],
3093 exact_match => 1 );
3094 } else {
3095 debug -> warn( level => 2,
3096 message => "dropping $cov" );
3097 $applicant_model -> _option_val_pos ( problem_numbers => [1],
3098 instance_numbers => [[1]],
3099 name => $cov,
3100 record_name => 'input',
3101 new_values => [['DROP']],
3102 exact_match => 1 );
3103 }
3104 }
3105 }
3106 end drop_undrop_covariates
3107
3108 # }}}
Perl-speaks-NONMEM CVS import