From 00195678e68de53b6866328655d49a9f75ccd521 Mon Sep 17 00:00:00 2001 From: pontus_pih Date: Wed, 2 Apr 2008 13:59:10 +0000 Subject: [PATCH] removed mc_subs --- lib/tool/mc_subs.pm | 938 ---------------------------------------------------- 1 file changed, 938 deletions(-) delete mode 100644 lib/tool/mc_subs.pm diff --git a/lib/tool/mc_subs.pm b/lib/tool/mc_subs.pm deleted file mode 100644 index 78badf7..0000000 --- a/lib/tool/mc_subs.pm +++ /dev/null @@ -1,938 +0,0 @@ -# {{{ include statements - -start include statements -use tool::cdd; -use tool::modelfit; -use Math::Random; -use Data::Dumper; -end include statements - -# }}} include statements - -# {{{ new - -start new -foreach my $attribute ( 'logfile', 'raw_results_file' ) { - if ( not( ref($this -> {$attribute}) eq 'ARRAY' or - ref($this -> {$attribute}) eq 'HASH' ) ) { - my $tmp = $this -> {$attribute}; - if ( not defined $tmp and $attribute eq 'logfile' ) { - $tmp = 'mcs.log'; - } - $this -> {$attribute} = []; - for ( my $i = 1; $i <= scalar @{$this -> {'models'}}; $i++ ) { - my $name = $tmp; - if ( $name =~ /\./ ) { - $name =~ s/\./$i\./; - } else { - $name = $name.$i; - } - my $ldir; - ( $ldir, $name ) = - OSspecific::absolute_path( $this -> {'directory'}, $name ); - push ( @{$this -> {$attribute}}, $ldir.$name ) ; - } - } -} - -if ( scalar (@{$this -> {'models'}->[0]-> problems}) != 1 ){ - 'debug' -> die(message => 'The original model must contain exactly one problem.'); -} - -foreach my $alt (@{$this -> {'alternative_models'}}){ - if ( scalar (@{$alt -> problems}) != 1 ){ - 'debug' -> die(message => 'The alternative models must each contain exactly one problem.'); - } -} - -end new - -# }}} - -# {{{ cdd_setup - -start cdd_setup - -my @subm_threads; -if (ref( $self -> {'threads'} ) eq 'ARRAY') { - @subm_threads = @{$self -> {'threads'}}; - shift(@subm_threads); -} else { - @subm_threads = ($self -> {'threads'}); -} -$self -> general_setup( model_number => $model_number, - class => 'tool::cdd', - subm_threads => \@subm_threads ); - -end cdd_setup - -# }}} cdd_setup - -# {{{ modelfit_setup - -start modelfit_setup - -my @subm_threads; -if (ref( $self -> {'threads'} ) eq 'ARRAY') { - @subm_threads = @{$self -> {'threads'}}; - shift(@subm_threads); -} else { - @subm_threads = ($self -> {'threads'}); -} -$self -> general_setup( model_number => $model_number, - class => 'tool::modelfit', - subm_threads => $subm_threads[$model_number-1] ); - -end modelfit_setup - -# }}} modelfit_setup - -# {{{ general_setup - -start general_setup - -my $model = $self -> {'models'} -> [$model_number-1]; -#my $alternative = $self -> {'alternative_models'} -> [$model_number-1]; -my @alternatives = @{$self -> {'alternative_models'}}; #may be empty -my $subm_threads = $parm{'subm_threads'}; -my $own_threads = ref( $self -> {'threads'} ) eq 'ARRAY' ? - $self -> {'threads'} -> [0]:$self -> {'threads'}; -my $done = ( -e $self -> {'directory'}."/m$model_number/done" ) ? 1 : 0; -my ( @seed, $new_datas, $skip_ids, $skip_keys, $skip_values ); -my ( @orig_est_models, @alt_est_models ); -my ( $sim_model, $est_original, $est_alternative ); -my $alternative; - -if ( not $done ) { - my @sim_models; - my @table_header; - my @all_simulated_files; - my ( @orig_table_names, @alt_table_names ); - -# Copy the model to a new simulation models - for( my $sim_no = 1; $sim_no <= $self -> {'samples'}; $sim_no++ ) { - my $sim_name = "mc-$sim_no.sim"; - my $orig_name = "mc-original-$sim_no.mod"; - my $out_name = "mc-$sim_no.lst"; - my $orig_out = "mc-original-$sim_no.lst"; - if( $sim_no == 1 ) { - $sim_model = $model -> - copy( filename => $self -> {'directory'}.'m'.$model_number.'/'.$sim_name, - target => 'disk', - copy_data => 1, - copy_output => 0); - $sim_model -> drop_dropped; - $est_original = $model -> - copy( filename => $self -> {'directory'}.'m'.$model_number.'/'.$orig_name, - target => 'disk', - copy_data => 0, - copy_output => 0); - $est_original -> drop_dropped; - $est_original -> remove_records( type => 'simulation' ); - - - #ignore @ since simdata contains header rows. can skip old ignores since filtered - $est_original -> set_option( record_name => 'data', - option_name => 'IGNORE', - option_value => '@'); - - my $tbl_nm_ref = $model -> table_names(); - @orig_table_names = @{$tbl_nm_ref} if( defined $tbl_nm_ref ); - } else { - $sim_model = $sim_models[0] -> - copy( filename => $self -> {'directory'}.'m'.$model_number.'/'.$sim_name, - target => 'disk', - copy_data => 0, - copy_output => 0); - $est_original = $orig_est_models[0] -> - copy( filename => $self -> {'directory'}.'m'.$model_number.'/'.$orig_name, - target => 'disk', - copy_data => 0, - copy_output => 0); - } - $sim_model -> ignore_missing_files( 1 ); - $sim_model -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$out_name ); - $sim_model -> ignore_missing_files( 0 ); - $est_original -> ignore_missing_files( 1 ); - $est_original -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$orig_out ); - $est_original -> ignore_missing_files( 0 ); - my @new_orig_table_names; - for( my $pr = 0; $pr <= $#orig_table_names; $pr++ ) { - for( my $tbl = 0; $tbl < scalar @{$orig_table_names[$pr]}; $tbl++ ) { - $new_orig_table_names[$pr][$tbl] = $orig_table_names[$pr][$tbl].'-'.$sim_no; - } - } - if( $#new_orig_table_names >= 0 ) { - $est_original -> table_names( new_names => \@new_orig_table_names ); - } - if( $self -> shrinkage() ) { - $est_original -> shrinkage_stats( enabled => 1 ); - $est_alternative -> shrinkage_stats( enabled => 1 ); - my @problems = @{$est_original -> problems}; - for( my $i = 1; $i <= scalar @problems; $i++ ) { - $problems[ $i-1 ] -> shrinkage_module -> model( $est_original ); - } - } - - my @simulated_files; - for( my $j = 0; $j < scalar @{$sim_model -> problems}; $j++ ) { - my $prob = $sim_model -> problems -> [$j]; - - # set $SIMULATION record - - my $sim_record = $sim_model -> record( problem_number => $j, - record_name => 'simulation' ); - - if( scalar(@{$sim_record}) > 0 ){ - - my @new_record; - foreach my $sim_line ( @{$sim_record -> [0]} ){ - my $new_line; - while( $sim_line =~ /([^()]*)(\([^()]+\))(.*)/g ){ - my $seed = random_uniform_integer( 1, 0, 1000000 ); # Upper limit is from nmhelp - - my $head = $1; - my $old_seed = $2; - $sim_line = $3; - - if( $old_seed =~ /\(\d+\s+(NORMAL|UNIFORM|NONPARAMETRIC)\)/ ){ - $new_line .= "$head($seed $1)"; - } else { - $new_line .= "$head($seed)"; - } - - } - - push( @new_record, $new_line.$sim_line ); - } - - $prob -> set_records( type => 'simulation', - record_strings => \@new_record ); - - } else { - - my $seed = random_uniform_integer( 1, 0, 1000000 ); # Upper limit is from nmhelp - $prob -> set_records( type => 'simulation', - record_strings => [ '(' . $seed . - ') ONLYSIMULATION' ] ); - } - - if( $sim_no == 1 ) { - # remove $EST and $COV - - if( $sim_model -> is_option_set( problem_number => $j, - record => 'estimation', - name => 'LIKELIHOOD', - fuzzy_match => 1 ) - or - $sim_model -> is_option_set( problem_number => $j, - record => 'estimation', - name => '-2LOGLIKELIHOOD', - fuzzy_match => 1 ) - or - $sim_model -> is_option_set( problem_number => $j, - record => 'estimation', - name => '-2LLIKELIHOOD', - fuzzy_match => 1 ) - ){ - - $sim_model -> set_option( problem_numbers => [$j], - record_name => 'estimation', - option_name => 'MAXEVALS', - fuzzy_match => 1, - option_value => 0 ); - - } else { - - $prob -> remove_records(type => 'estimation'); - - } - - $prob -> remove_records(type => 'covariance'); - - # set $TABLE record - if( defined $prob -> inputs and defined $prob -> inputs -> [0] -> options ) { - foreach my $option ( @{$prob -> inputs -> [0] -> options} ) { - push( @table_header, $option -> name ); - } - } else { - debug -> die( message => "Trying to construct table for monte-carlo simulation". - " but no headers were found in \$model_number-1NPUT" ); - } - } - my $simulated_file = "mc-sim-".($j+1)."-$sim_no.dat"; - $prob -> set_records( type => 'table', - record_strings => [ join( ' ', @table_header ). - ' NOPRINT NOAPPEND ONEHEADER FILE='. - $simulated_file ] ); - push( @simulated_files, $self -> {'directory'}.'m'.$model_number.'/'. - $simulated_file ); - } - $sim_model -> _write( write_data => $sim_no == 1 ? 1 : 0 ); - push( @sim_models, $sim_model ); - push( @orig_est_models, $est_original ); - push( @all_simulated_files, \@simulated_files ); - } - my $mod_sim = tool::modelfit -> new( reference_object => $self, - models => \@sim_models, - base_directory => $self -> {'directory'}, - directory => $self -> {'directory'}. - 'simulation_dir'.$model_number, - parent_tool_id => $self -> {'tool_id'}, - retries => 1, - logfile => undef, - raw_results => undef, - prepared_models => undef, - threads => 1 ); - - $mod_sim -> run; - for( my $j = 0; $j <= $#orig_est_models; $j++ ) { - $orig_est_models[$j] -> datafiles( new_names => $all_simulated_files[$j] ); - $orig_est_models[$j] -> _write; - $orig_est_models[$j] -> flush_data(); - } - @{$self -> {'mc_models'}} = ( @orig_est_models ); - - foreach $alternative (@alternatives){ - @alt_table_names = (); - @alt_est_models =(); - my $filestem = 'mc-alt-' . $alternative -> filename(); - $filestem =~ s/\.[^.]*$//; - for( my $sim_no = 1; $sim_no <= $self -> {'samples'}; $sim_no++ ) { - my $alt_name = "$filestem-$sim_no.mod"; - my $alt_out = "$filestem-$sim_no.lst"; - - if( $sim_no == 1 ) { - $est_alternative = $alternative -> - copy( filename => $self -> {'directory'}.'m'.$model_number.'/'.$alt_name, - target => 'disk', - copy_data => 0, - copy_output => 0); - $est_alternative -> drop_dropped; - $est_alternative -> remove_records( type => 'simulation' ); - - #ignore @ since simdata contains header rows. cannot skip other old ignores!! - - my $ig_list = $est_alternative -> get_option_value( record_name => 'data', - option_name => 'IGNORE', - option_index => 'all'); - $est_alternative -> remove_option( record_name => 'data', - option_name => 'IGNORE'); - - if (scalar (@{$ig_list})>0){ - foreach my $val (@{$ig_list}){ - unless ($val =~ /^.$/){ - $est_alternative -> add_option( record_name => 'data', - option_name => 'IGNORE', - option_value => $val); - } - } - } - $est_alternative -> add_option( record_name => 'data', - option_name => 'IGNORE', - option_value => '@'); - ##done fixing ignore - - my $tbl_nm_ref = $alternative -> table_names(); - @alt_table_names = @{$tbl_nm_ref} if( defined $tbl_nm_ref ); - } else { - $est_alternative = $alt_est_models[0] -> - copy( filename => $self -> {'directory'}.'m'.$model_number.'/'.$alt_name, - target => 'disk', - copy_data => 0, - copy_output => 0); - } - $est_alternative -> ignore_missing_files( 1 ); - $est_alternative -> outputfile( $self -> {'directory'}.'m'.$model_number.'/'.$alt_out ); - $est_alternative -> ignore_missing_files( 0 ); - my @new_alt_table_names; - for( my $pr = 0; $pr <= $#alt_table_names; $pr++ ) { - for( my $tbl = 0; $tbl < scalar @{$alt_table_names[$pr]}; $tbl++ ) { - $new_alt_table_names[$pr][$tbl] = $alt_table_names[$pr][$tbl].'-'.$sim_no; - } - } - if( $#new_alt_table_names >= 0 ) { - $est_alternative -> table_names( new_names => \@new_alt_table_names ); - } - if( $self -> shrinkage() ) { - my @problems = @{$est_alternative -> problems}; - for( my $i = 1; $i <= scalar @problems; $i++ ) { - $problems[ $i-1 ] -> shrinkage_module -> model( $est_alternative ); - } - } - push( @alt_est_models, $est_alternative ); - } #end loop over sim_no - #is all_simulated_files still in scope? should be, declared inside if not $done - for( my $j = 0; $j <= $#alt_est_models; $j++ ) { - $alt_est_models[$j] -> datafiles( new_names => $all_simulated_files[$j] ); - $alt_est_models[$j] -> _write; - $alt_est_models[$j] -> flush_data(); - } - push (@{$self -> {'mc_models'}}, @alt_est_models); - - - }#end loop over alternatives - - # Create a checkpoint. - open( DONE, ">".$self -> {'directory'}."/m$model_number/done" ) ; - print DONE "Simulation from ",$sim_models[0] -> filename," through ", - $sim_models[$#sim_models] -> filename," performed\n"; - print DONE $self -> {'samples'}." samples\n"; - @seed = random_get_seed; - print DONE "seed: @seed\n"; - close( DONE ); -} else { - #$done=true - # Recreate the datasets and models from a checkpoint - ui -> print( category => 'mc', - message => "Recreating models from a previous run" ); - open( DONE, $self -> {'directory'}."/m$model_number/done" ); - my @rows = ; - close( DONE ); - my ( $junk, $junk, $stored_filename1, $junk, - $stored_filename2, $junk ) = split(' ',$rows[0],4); - my ( $stored_samples, $junk ) = split(' ',$rows[1],2); - @seed = split(' ',$rows[2]); - shift( @seed ); # get rid of 'seed'-word - # Reinitiate the model objects - for ( my $j = 1; $j <= $stored_samples; $j++ ) { - my ($model_dir, $orig_name) = OSspecific::absolute_path( $self -> {'directory'}.'/m'. - $model_number, - "mc-original-$j.mod" ); - my ($out_dir, $orig_out) = OSspecific::absolute_path( $self -> {'directory'}.'/m'. - $model_number, - "mc-original-$j.lst" ); - my $est_original = model -> - new( directory => $model_dir, - filename => $orig_name, - outputfile => $orig_out, - extra_files => $model -> extra_files, - target => 'disk', - ignore_missing_files => 1, - extra_output => $model -> extra_output(), - cwres => $model -> cwres(), - quick_reload => 1 ); - push( @orig_est_models, $est_original ); - my $nl = $j == $stored_samples ? "" : "\r"; - ui -> print( category => 'mc', - message => ui -> status_bar( sofar => $j+1, - goal => $stored_samples+1 ).$nl, - wrap => 0, - newline => 0 ); - } #end loop over samples - @{$self -> {'mc_models'}} = ( @orig_est_models); - - #start alternatives section - foreach $alternative (@alternatives){ - my $filestem = 'mc-alt-' . $alternative -> filename(); - $filestem =~ s/\.[^.]*$//; - @alt_est_models =(); - for ( my $j = 1; $j <= $stored_samples; $j++ ) { - my ($model_dir, $alt_name) = OSspecific::absolute_path( $self -> {'directory'}.'/m'. - $model_number, - "$filestem-$j.mod" ); - my ($out_dir, $alt_out) = OSspecific::absolute_path( $self -> {'directory'}.'/m'. - $model_number, - "$filestem-$j.lst" ); - my $est_alternative = model -> - new( directory => $model_dir, - filename => $alt_name, - outputfile => $alt_out, - extra_files => $alternative -> extra_files, - target => 'disk', - ignore_missing_files => 1, - extra_output => $model -> extra_output(), - cwres => $model -> cwres(), - quick_reload => 1 ); - push( @alt_est_models, $est_alternative ); - my $nl = $j == $stored_samples ? "" : "\r"; - ui -> print( category => 'mc', - message => ui -> status_bar( sofar => $j+1, - goal => $stored_samples+1 ).$nl, - wrap => 0, - newline => 0 ); - } #end loop over samples - push (@{$self -> {'mc_models'}}, @alt_est_models ); - - }#end loop over alternatives - - ui -> print( category => 'mc', message => " ... done." ); - random_set_seed( @seed ); - ui -> print( category => 'mc', - message => "Using $stored_samples data sets, previously simulated ". - "from $stored_filename1 through $stored_filename2" ) - unless $self -> {'parent_threads'} > 1; -} #end $done -my $subdir = $class; -$subdir =~ s/tool:://; -my @subtools = @{$self -> {'subtools'}}; -shift( @subtools ); -my %subargs = (); -if ( defined $self -> {'subtool_arguments'} ) { - %subargs = %{$self -> {'subtool_arguments'}}; -} - -push( @{$self -> {'tools'}}, - $class -> new - ( reference_object => $self, - logfile => undef, - raw_results => undef, - prepared_models => undef, - models => $self -> {'mc_models'}, - nm_version => $self -> {'nm_version'}, - parent_tool_id => $self -> {'tool_id'}, - threads => $subm_threads, - base_directory => $self -> {'directory'}, - directory => $self -> {'directory'}.'/'.$subdir.'_dir'.$model_number, - subtools => $#subtools >= 0 ? \@subtools : undef, - parent_threads => $own_threads, - shrinkage => $self -> {'shrinkage'}, - _raw_results_callback => $self -> - _modelfit_raw_results_callback( model_number => $model_number ), - %subargs ) ); - -end general_setup - -# }}} general_setup - -# {{{ cdd_analyze -start cdd_analyze - $self -> {'tools'}[0] -> print_results; -end cdd_analyze -# }}} cdd_analyze - -# {{{ modelfit_analyze - -start modelfit_analyze -$self -> {'tools'}[0] -> print_results; -end modelfit_analyze - -# }}} - -# {{{ _modelfit_raw_results_callback - -start _modelfit_raw_results_callback - -# Use the mc's raw_results file. -my ($dir,$file) = - OSspecific::absolute_path( $self -> {'directory'}, - $self -> {'raw_results_file'}[$model_number-1] ); -my $orig_mod = $self -> {'models'}[$model_number-1]; -$subroutine = sub { - my $modelfit = shift; - my $mh_ref = shift; - my %max_hash = %{$mh_ref}; - $modelfit -> raw_results_file( $dir.$file ); - - # {{{ New header - - if ( defined $modelfit -> raw_results() ) { - my @rows = @{$modelfit -> raw_results()}; - - #the first 'samples' (a scalar) number of rows is null models - for( my $i = 0; $i < $self -> {'samples'}; $i++ ) { - unshift( @{$rows[$i]}, 'null' ); - } - #remaining rows, if any, belong to alternatives model(s) - my $offset = $self -> {'samples'}; - foreach my $alternative (@{$self -> {'alternative_models'}}){ - if (($offset + $self -> {'samples'}-1)> $#rows){ - debug -> die( message => "Too few rows in raw results file." ); - } - my $altname = 'alt-' . $alternative -> filename(); - $altname =~ s/\.[^.]*$//; - for( my $i = 0; $i < $self -> {'samples'}; $i++ ) { - unshift( @{$rows[($offset+$i)]}, $altname ); - } - $offset += $self -> {'samples'}; - } - if (($offset -1)< $#rows){ - debug -> die( message => "Too many rows in raw results file." ); - } - - unshift( @{$modelfit -> {'raw_results_header'}}, 'hypothesis' ); - - } - @{$self -> {'raw_results_header'}} = @{$modelfit -> {'raw_results_header'}}; - # }}} New header - -}; -return $subroutine; - -end _modelfit_raw_results_callback - -# }}} _modelfit_raw_results_callback - -# {{{ prepare_results - -start prepare_results -{ - # Mean - # Median - # SD - - $self -> {'raw_results'} = $self -> {'tools'}[0] -> {'raw_results'}; - - my $model = $self -> {'models'} -> [0]; - my @alternatives = @{$self -> {'alternative_models'}}; - - my $samples = $self -> {'samples'}; - - my $raw_line_struct = $self -> {'tools'}[0] -> {'raw_line_structure'}; - - my %orig_results_section; - my @alt_results_sections; - - ### Retrieve initial parameter estimates from original model. - - my %initial_values; - my %n_initials; - my (@thetalist,@omegalist,@sigmalist); - my $tmp; - - $n_initials{'theta'} = $model->nthetas(problem_number => 1); - $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'theta'); - @thetalist=@{$tmp->[0]}; - 'debug'->die(message=>'Error number of initial thetas.') unless (scalar(@thetalist) == $n_initials{'theta'}); - $initial_values{'theta'}=\@thetalist; - - $tmp = $model->nomegas(problem_numbers => [1]); - $n_initials{'omega'} = $tmp->[0]; - $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'omega'); - @omegalist=@{$tmp->[0]}; - 'debug'->die(message=>'Error number of initial omegas.') unless (scalar(@omegalist) == $n_initials{'omega'}); - $initial_values{'omega'}=\@omegalist; - - $tmp = $model->nsigmas(problem_numbers => [1]); - $n_initials{'sigma'} = $tmp->[0]; - $tmp = $model->initial_values(problem_numbers => [1], parameter_type => 'sigma'); - @sigmalist=@{$tmp->[0]}; - 'debug'->die(message=>'Error number of initial sigmas.') unless (scalar(@sigmalist) == $n_initials{'sigma'}); - $initial_values{'sigma'}=\@sigmalist; - - $n_initials{'ofv'}=0; - - - ### Loop over the original model (the zero) and over all alternative - ### models. - - for my $model_index ( 0..scalar(@alternatives) ){ - - ### Foreach such model, create a %results_section - - my %results_section; - - ### Each results_section contains statistics of ofv and - ### parameters. - - my @values; - my @labels; - - @{$labels[0]} = ('mean','median','sd','skewness','kurtosis','rmse','bias'); - - if( $model_index == 0 ){ - $results_section{'name'} = 'Original model'; - } else { - $results_section{'name'} = 'Alternative model '.$model_index; - } - - $results_section{'labels'} = \@labels; - $results_section{'values'} = \@values; - - push( @{$self -> {'results'}[0]{'own'}},\%results_section ); - - foreach my $measure ( 'ofv','theta','omega','sigma' ){ - - # Create a header, we might use max_hash from modelfit for this - # later on. - - my $start=0; - my $length=0; - - for( 0..($samples-1) ){ - my ($s, $l) = split(/,/, $raw_line_struct -> {$model_index*$samples + $_} -> {$measure}); - - if( defined $s and $start != 0 and $s != $start ){ - 'debug' -> die( message => "Internal structure corrupted at $model_index, $samples, $_, $measure, this is a bug" ); - } else { - $start = $s; - } - - if( defined $l and $l > $length ){ - $length = $l; - } - - } - - if( $start == 0 ){ - for( 0..6 ){ - push(@{$values[$_]},'Model not terminated'); - } - push( @{$labels[1]},$measure ); - next; - } - - if( $measure eq 'ofv' ){ - push( @{$labels[1]}, $measure ); - } else { - for( 1..$length ){ - push( @{$labels[1]}, uc(substr($measure, 0,2))."_$_" ); - } - } - - my $initial_index=0; - foreach my $col ( $start..($start + $length-1) ){ - my ($skew, $kurt, $mean, $stdev) = $self -> skewness_and_kurtosis( column_index => $col+1, - start_row_index => $model_index*$samples, - end_row_index => $model_index*$samples+$samples-1 ); - - push( @{$values[0]}, $mean ); - push( @{$values[2]}, $stdev ); - push( @{$values[3]}, $skew ); - push( @{$values[4]}, $kurt ); - - my $median = $self -> median( column_index => $col+1, - start_row_index => $model_index*$samples, - end_row_index => $model_index*$samples+$samples-1 ); - - push( @{$values[1]}, $median ); - - if( $initial_index < $n_initials{$measure}){ - - my $initial_value= $initial_values{$measure}->[$initial_index]; - - $initial_index++; - - my $rmse = $self -> rmse_percent( column_index => $col+1, - start_row_index => $model_index*$samples, - end_row_index => $model_index*$samples+$samples-1, - initial_value => $initial_value ); - - push( @{$values[5]}, $rmse ); - - my $bias = $self -> bias_percent( column_index => $col+1, - start_row_index => $model_index*$samples, - end_row_index => $model_index*$samples+$samples-1, - initial_value => $initial_value ); - - - push( @{$values[6]}, $bias ); - }else { - #push placeholders so that values will stay in correct column - push( @{$values[5]}, ' ' ); - push( @{$values[6]}, ' ' ); - } - - } - } - } - - ## Calculate OFV mean calculations - - my ($s, $l) = split(/,/, $raw_line_struct -> {0} -> {'ofv'}); - - my @ofv_limits = (-3.84,-5.99,-7.81,-9.49); - - my %results_section; - my (@labels, @values); - - @{$labels[1]} = ('mean delta_OFV'); - - foreach my $limit( @ofv_limits ){ - push( @{$labels[1]}, "delta_OFV < $limit" ); - } - - $results_section{'name'} = 'OFV Statistics'; - $results_section{'labels'} = \@labels; - $results_section{'values'} = \@values; - - push( @{$self -> {'results'}[0]{'own'}},\%results_section ); - - foreach my $model_index( 1..scalar(@alternatives) ){ - - push( @{$labels[0]}, "original - alternative$model_index" ); - - my $start=0; - - my $delta_sum=0; - - my %nr_below; - - foreach my $sample( 0..($samples-1) ){ - my ($s, $l) = split(/,/, $raw_line_struct -> {$model_index*$samples + $sample} -> {'ofv'}); - - if( $s ne '' and $start != 0 and $s != $start ){ - 'debug' -> die( message => "OFV: Internal structure corrupted at $model_index, $samples, $sample, this is a bug" ); - } else { - $start = $s; - } - - if( $start eq '' ){ - # This model probably did not terminate. - next; - } - - my $alt_ofv = $self -> {'raw_results'} -> [$model_index*$samples + $sample][$s+1]; - my $org_ofv = $self -> {'raw_results'} -> [$sample][$s+1]; - my $delta = $org_ofv - $alt_ofv; - $delta_sum += $delta; - - foreach my $limit ( @ofv_limits ){ - - if( $delta < $limit ){ - $nr_below{$limit}++; - } - - } - - } - - push( @{$values[$model_index-1]}, ($delta_sum / ($samples * scalar(@alternatives))) . "%" ); - foreach my $limit( @ofv_limits ){ - push( @{$values[$model_index-1]}, ($nr_below{$limit} / ($samples * scalar(@alternatives)))*100 . "%" ); - } - } - - -} -end prepare_results - -# }}} - -# {{{ rmse_percent - -start rmse_percent -{ - #input is integers $column_index, $start_row_index, $end_row_index and scalar float $initial_value - #output is scalar $rmse_percent - - unless( $end_row_index ){ - $end_row_index = $#{$self -> {'raw_results'}}; - } - - debug -> die( message => "Bad row index input") if ($start_row_index > $end_row_index); - my $row_count = $end_row_index - $start_row_index +1; - my $sum_squared_errors=0; - for (my $i=$start_row_index; $i<=$end_row_index; $i++){ - $sum_squared_errors += ($self->{'raw_results'}->[$i][$column_index] - $initial_value)**2; - } - - if ($initial_value == 0){ - debug -> warn( level => 2, - message => "Initial value 0, returning absolute rmse instead of relative."); - $rmse_percent= sqrt($sum_squared_errors/$row_count); - }else{ - $rmse_percent= (sqrt($sum_squared_errors/$row_count))*100/abs($initial_value); - } -} -end rmse_percent - -# }}} rmse_percent - -# {{{ bias_percent - -start bias_percent -{ - #input is integers $column_index, $start_row_index, $end_row_index and scalar float $initial_value - #output is scalar $bias_percent - - unless( $end_row_index ){ - $end_row_index = $#{$self -> {'raw_results'}}; - } - - debug -> die( message => "Bad row index input") if ($start_row_index > $end_row_index); - - my $row_count = $end_row_index - $start_row_index +1; - my $sum_errors=0; - for (my $i=$start_row_index; $i<=$end_row_index; $i++){ - $sum_errors += ($self->{'raw_results'}->[$i][$column_index] - $initial_value); - } - - if ($initial_value == 0){ - debug -> warn( level => 2, - message => "Initial value 0, returning absolute bias instead of relative."); - $bias_percent= ($sum_errors/$row_count); - }else{ - $bias_percent= ($sum_errors/$row_count)*100/abs($initial_value); - } -} -end bias_percent - -# }}} bias_percent - -# {{{ skewness_and_kurtosis - -start skewness_and_kurtosis -{ - #input is integers $column_index, $start_row_index, $end_row_index - - unless( $end_row_index ){ - $end_row_index = $#{$self -> {'raw_results'}}; - } - - debug -> die( message => "Bad row index input") if ($start_row_index >= $end_row_index); - - my $row_count = $end_row_index - $start_row_index +1; - my $sum_values=0; - for (my $i=$start_row_index; $i<=$end_row_index; $i++){ - $sum_values += $self->{'raw_results'}->[$i][$column_index]; - } - - $mean=$sum_values/$row_count; - - my $error=0; - my $sum_errors_pow2=0; - my $sum_errors_pow3=0; - my $sum_errors_pow4=0; - - for (my $i=$start_row_index; $i<=$end_row_index; $i++){ - $error = ($self->{'raw_results'}->[$i][$column_index] - $mean); - $sum_errors_pow2 += $error**2; - $sum_errors_pow3 += $error**3; - $sum_errors_pow4 += $error**4; - } - - ## TODO fråga om missing values. och om SD - - $stdev=0; - unless( $sum_errors_pow2 == 0 ){ - - $stdev= sqrt ($sum_errors_pow2/($row_count-1)); - $skewness = $sum_errors_pow3/($row_count*($stdev**3)); - $kurtosis = -3 + $sum_errors_pow4/($row_count*($stdev**4)); - - } -} -end skewness_and_kurtosis - -# }}} skewness_and_kurtosis - -# {{{ median - -start median -{ - #input is integers $column_index, $start_row_index, $end_row_index - - unless( $end_row_index ){ - $end_row_index = $#{$self -> {'raw_results'}}; - } - - debug -> die( message => "Bad row index input") if ($start_row_index >= $end_row_index); - - my @temp; - - for (my $i=$start_row_index; $i<=$end_row_index; $i++){ - push( @temp, $self->{'raw_results'}->[$i][$column_index] ); - } - - @temp = sort({$a <=> $b} @temp); - if( scalar( @temp ) % 2 ){ - $median = $temp[$#temp/2]; - } else { - $median = ($temp[@temp/2]+$temp[(@temp-2)/2]) / 2; - } - -} -end median - -# }}} median -- 2.11.4.GIT