IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 4, 2007, 1:53:22 PM (19 years ago)
Author:
eugene
Message:

extensive changes to use new stats calculation tools; general cleanups

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippScripts/scripts/detrend_reject_exp.pl

    r13962 r14009  
    4747
    4848pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
    49 pod2usage(
    50     -msg => "Required options: --det_id --iteration --det_type --camera",
    51     -exitval => 3,
    52 ) unless defined $det_id
     49pod2usage( -msg => "Required options: --det_id --iteration --det_type --camera",
     50           -exitval => 3)
     51    unless defined $det_id
    5352    and defined $iter
    5453    and defined $det_type
    5554    and defined $camera;
    5655
     56# values to extract from output metadata and the stats to calculate
     57# XXX -bg_mean_stdev should take rms of bg_mean_stdev if bg_mean_stdev != 0 (A)
     58# XXX -bg_mean_stdev should take stdev of bg_mean if bg_mean_stdev == 0     (B)
     59# XXX  (A) if imfile.Ncomp > 1, (B) if imfile.Ncomp == 1
     60my $STATS =
     61   [   
     62       #          KEYWORD                 STATISTIC          CHIPTOOL FLAG
     63       { name => "bg",             type => "mean",  flag => "-bg" },
     64       { name => "bg_mean_stdev",  type => "stdev", flag => "-bg_mean_stdev" },
     65       { name => "bg_stdev",       type => "rms",   flag => "-bg_stdev" },
     66   ];
     67
     68my $REJSTATS =
     69   [   
     70       #          KEYWORD                 STATISTIC          CHIPTOOL FLAG
     71       { name => "bg",             type => "mean",   flag => "ensMeanMean" },
     72       { name => "bg",             type => "stdev",  flag => "ensMeanStdev" },
     73       { name => "bg_mean_stdev",  type => "mean",   flag => "ensMeanStdevMean" },
     74       { name => "bg_mean_stdev",  type => "stdev",  flag => "ensMeanStdevStdev" },
     75       { name => "bg_stdev",       type => "mean",   flag => "ensStdevMean" },
     76       { name => "bg_stdev",       type => "stdev",  flag => "ensStdevStdev" },
     77   ];
    5778
    5879# Look for programs we need
    … …  
    6485}
    6586
    66 my $mdcParser = PS::IPP::Metadata::Config->new; # Parser for metadata config files
    67 
    6887# Get list of component files
    6988my $exposures;                  # Array of exposures
    7089{
    71     my $command = "$dettool -residexp -det_id $det_id -iteration $iter"; # Command to run
     90    # dettool command to select exp data for this det_run
     91    my $command = "$dettool -residexp";
     92    $command .= " -det_id $det_id";
     93    $command .= " -iteration $iter";
    7294    $command .= " -dbname $dbname" if defined $dbname;
    7395    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    … …  
    78100    }
    79101
     102    # Parse the stdout buffer into a metadata
     103    my $mdcParser = PS::IPP::Metadata::Config->new;     # Parser for metadata config files
    80104    my $metadata = $mdcParser->parse(join "", @$stdout_buf) or
    81105        &my_die("Unable to parse metadata config doc", $det_id, $iter, $PS_EXIT_PROG_ERROR);
    82     # XXX need to test that the result list is 0 length and return with a different error
     106
     107    # parse the file info in the metadata
    83108    $exposures = parse_md_list($metadata) or
    84109        &my_die("Unable to parse metadata list", $det_id, $iter, $PS_EXIT_PROG_ERROR);
     110
     111    # Parse the statistics on the residual image
     112    my $stats = PS::IPP::Metadata::Stats->new($STATS); # Stats parser
     113    $stats->parse($metadata) or &my_die("Unable to find all values in statistics output.", $det_id, $iter, $PS_EXIT_PROG_ERROR);
     114
     115    # Parse the statistics for rejections
     116    my $rejstats = PS::IPP::Metadata::Stats->new($REJSTATS); # Stats parser
     117    $rejstats->parse($metadata) or &my_die("Unable to find all values in statistics output.", $det_id, $iter, $PS_EXIT_PROG_ERROR);
    85118}
    86119
    87120my @expTags;                    # Array of exposure IDs
    88 my @means;                      # Array of means
    89 my @stdevs;                     # Array of standard deviations
    90 my @meanStdevs;                 # Array of standard deviations of the mean (normalised by the mean)
    91 my @variances;                  # Array of variances
    92121my @accept;                     # Array of accept flags
    93122my @include;                    # Array of include flags
    94123foreach my $exposure (@$exposures) {
    95124    &my_die("Unable to find exposure id.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{exp_tag};
    96     &my_die("Unable to find mean.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{bg};
    97     &my_die("Unable to find stdev.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{bg_stdev};
    98     &my_die("Unable to find mean stdev.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{bg_mean_stdev};
    99125    &my_die("Unable to find accept.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{accept};
    100126    &my_die("Unable to find include.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{include};
    101     push @expTags, $exposure->{exp_tag};
    102     push @means, $exposure->{bg};
    103     push @stdevs, $exposure->{bg_stdev};
    104     # XXX why are we using mean_stdev / bd?
    105     # push @meanStdevs, $exposure->{bg_mean_stdev} / $exposure->{bg};
    106     push @meanStdevs, $exposure->{bg_mean_stdev};
    107     push @variances, $exposure->{bg_stdev}**2;
    108     push @accept, $exposure->{accept};
    109     push @include, $exposure->{include};
    110 }
    111 my $meanStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
    112 $meanStats->add_data(@means);
    113 my $stdevStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
    114 $stdevStats->add_data(@stdevs);
    115 my $meanStdevStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
    116 $meanStdevStats->add_data(@meanStdevs);
     127}
     128
     129my $ensMeanMean       = &STATS_value_for_flag ($REJSTATS, "ensMeanMean");     
     130my $ensMeanStdev      = &STATS_value_for_flag ($REJSTATS, "ensMeanStdev");     
     131my $ensMeanStdevMean  = &STATS_value_for_flag ($REJSTATS, "ensMeanStdevMean");
     132my $ensMeanStdevStdev = &STATS_value_for_flag ($REJSTATS, "ensMeanStdevStdev");
     133my $ensStdevMean      = &STATS_value_for_flag ($REJSTATS, "ensStdevMean");     
     134my $ensStdevStdev     = &STATS_value_for_flag ($REJSTATS, "ensStdevStdev");   
    117135
    118136$ipprc->define_camera($camera);
    119137# Rejection thresholds
    120 my $reject_mean = rejection_limit( 'ENSEMBLE.MEAN', $det_type, $filter );
    121 my $reject_stdev = rejection_limit( 'ENSEMBLE.STDEV', $det_type, $filter );
     138my $reject_mean      = rejection_limit( 'ENSEMBLE.MEAN',      $det_type, $filter );
     139my $reject_stdev     = rejection_limit( 'ENSEMBLE.STDEV',    $det_type, $filter );
    122140my $reject_meanstdev = rejection_limit( 'ENSEMBLE.MEANSTDEV', $det_type, $filter );
    123141
    … …  
    129147unless ($no_op) {
    130148    $logFile = $ipprc->file_create_append( $logName );
    131     print $logFile "Ensemble mean " . $meanStats->mean() . " +/- " . $meanStats->standard_deviation .
    132         ", stdev " . $stdevStats->mean() . " +/- " . $stdevStats->standard_deviation() . "\n\n";
     149    print $logFile "Ensemble mean $meanEnsemble +/- $stdevEnsemble, stdev $meanStdevEnsemble\n\n";
    133150}
    134151
    … …  
    136153my $numChanges = 0;             # Number of exposures with changed status
    137154my $numReject = 0;              # Number of exposures rejected
    138 my $command;
    139 my $reject;
    140 my $expTag;
    141 for (my $i = 0; $i < scalar @means; $i++) {
    142     $expTag = $expTags[$i];     # Exposure ID
    143     $command = "$dettool -updateresidexp -det_id $det_id -iteration $iter -exp_tag $expTag"; # Command to run
     155
     156for (my $i = 0; $i < scalar @$exposures; $i++) {
     157    my $file      = $files[$i];
     158    my $mean      = $file->{bg};        # Mean for this exposure
     159    my $stdev     = $file->{bg_stdev}; # Stdev for this exposure
     160    my $meanStdev = $file->{bg_mean_stdev}; # Stdev of Means for this exposure
     161
     162    my $expTag  = $exposure->{exp_tag};
     163    my $accept  = $exposure->{accept};
     164    my $include = $exposure->{include};
     165
     166    my $reject = 0;             # Reject this exposure?
     167
     168    my $command = "$dettool -updateresidexp";
     169    $command .= " -det_id $det_id";
     170    $command .= " -iteration $iter";
     171    $command .= " -exp_tag $expTag";
    144172    $command .= " -dbname $dbname" if defined $dbname;
    145173   
    146     $reject = 0;                # Reject this exposure?
    147 
    148     if (not $accept[$i]) {
     174    if (not $accept) {
    149175        # Rejected this at an earlier stage
    150176        unless ($no_op) {
    … …  
    165191    }
    166192   
    167     if ($reject_mean > 0 and defined $meanStats->standard_deviation() ) {
    168         my $dMean = abs($means[$i] - $meanStats->mean()) ;
    169         if ($dMean > ($reject_mean * $meanStats->standard_deviation())) {
     193    if ($reject_mean > 0 and $ensMeanStdev > 0) {
     194        my $delta = abs($mean - $ensMeanMean);
     195        if ($delta > ($reject_mean * $ensMeanStdev)) {
    170196            print $logFile "Rejecting $expTag based on ensemble mean value: ";
    171197            $reject = 1;
    … …  
    174200            print $logFile "$expTag OK against ensemble mean: ";
    175201        }
    176         print $logFile "$means[$i] --> $dMean vs " . $reject_mean * $meanStats->standard_deviation() . "\n";
     202        print $logFile "$mean --> $delta vs " . $reject_mean * $ensMeanStdev . "\n";
    177203    } else {
    178204        print $logFile "No rejection of $expTag for ensemble mean\n";
    179205    }
    180206
    181     if ($reject_stdev > 0 and defined $stdevStats->standard_deviation() > 0) {
    182         my $dMean = abs($stdevs[$i] - $stdevStats->mean());
    183         if ($dMean > ($reject_stdev * $stdevStats->standard_deviation())) {
     207    if ($reject_stdev > 0 and $ensStdevStdev > 0) {
     208        my $delta = abs($stdev - $ensStdevMean);
     209        if ($delta > ($reject_stdev * $ensStdevStdev)) {
    184210            print $logFile "Rejecting $expTag based on ensemble stdev: ";
    185211            $reject = 1;
    … …  
    188214            print $logFile "$expTag OK against ensemble stdev: ";
    189215        }
    190         print $logFile "$stdevs[$i] --> $dMean sigma vs " . $reject_stdev  * $stdevStats->standard_deviation() . "\n";
     216        print $logFile "$stdev --> $delta sigma vs " . $reject_stdev * $ensStdevStdev . "\n";
    191217    } else {
    192218        print $logFile "No rejection of $expTag for ensemble stdev\n";
    193219    }
    194220   
    195     if ($reject_meanstdev > 0 and defined $meanStdevStats->standard_deviation() > 0) {
    196         my $dMean = abs($meanStdevs[$i] - $meanStdevStats->mean());
    197         if ($dMean > ($reject_meanstdev * $meanStdevStats->standard_deviation())) {
     221    if ($reject_meanstdev > 0 and $ensMeanStdevStdev > 0) {
     222        my $delta = abs($meanStdev - $ensMeanStdevMean);
     223        if ($delta > ($reject_meanstdev * $ensMeanStdevStdev)) {
    198224            print $logFile "Rejecting $expTag based on ensemble mean stdev: ";
    199225            $reject = 1;
    … …  
    202228            print $logFile "$expTag OK against ensemble mean stdev: ";
    203229        }
    204         print $logFile "$meanStdevs[$i] --> $dMean sigma vs " . $reject_meanstdev . "\n";
     230        print $logFile "$meanStdev --> $delta sigma vs " . $reject_meanstdev * $ensMeanStdevStdev. "\n";
    205231    } else {
    206232        print $logFile "No rejection of $expTag for ensemble mean stdev\n";
    … …  
    214240   
    215241    # Check for status changes
    216     if ((not $include[$i] and not $reject) or ($include[$i] and $reject)) {
     242    if ((not $include and not $reject) or ($include and $reject)) {
    217243        unless ($no_op) {
    218244            print $logFile "Status of $expTag has changed.\n";
    … …  
    253279}
    254280
     281my $command = "$dettool -adddetrunsummary";
     282$command .= " -det_id $det_id";
     283$command .= " -iteration $iter";
     284$command .= " -accept" if $master;
     285$command .= " -dbname $dbname" if defined $dbname;
     286
     287# add in the elements from the selected stats above
     288foreach my $entry (@$STATS) {
     289    my $value = $entry->{value};
     290    my $flag = $entry->{flag};
     291    $command .= " $flag $value";
     292}
     293
    255294# Put results into the database
    256295unless ($no_update) {
    257 
    258     {
    259         # Add the summary
    260         my $varianceStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
    261         $varianceStats->add_data(@variances);
    262        
    263         my $bg = ($meanStats->mean() or 'NAN');
    264         my $bg_stdev = (sqrt( $varianceStats->mean() ) or 'NAN');
    265         my $bg_mean_stdev = ($meanStats->standard_deviation() or 'NAN');
    266        
    267         my $command = "$dettool -adddetrunsummary -det_id $det_id -iteration $iter";
    268         $command .= " -bg $bg -bg_stdev $bg_stdev -bg_mean_stdev $bg_mean_stdev";
    269         $command .= " -accept" if $master;
    270         $command .= " -dbname $dbname" if defined $dbname;
    271        
    272         my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    273             run(command => $command, verbose => 1);
    274         unless ($success) {
    275             $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
    276             warn("Unable to perform dettool -adddetrunsummary: $error_code");
    277             exit($error_code);
    278         }
    279     }
    280 
    281     # Re-run processing if required
    282     {
    283         my $command = "$dettool -updatedetrun -det_id $det_id";
    284         if ($stop) {
    285             $command .= ' -state stop';
    286         } else {
    287             $command .= ' -again';
    288         }
    289        
    290         $command .= " -dbname $dbname" if defined $dbname;
    291        
    292         my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    293             run(command => $command, verbose => 1);
    294         unless ($success) {
    295             $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
    296             warn("Unable to perform dettool -updatedetrun: $error_code");
    297             exit($error_code);
    298         }
    299     }
     296    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
     297        run(command => $command, verbose => 1);
     298    unless ($success) {
     299        $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
     300        warn("Unable to perform dettool -adddetrunsummary: $error_code");
     301        exit($error_code);
     302    }
     303
     304    $command = "$dettool -updatedetrun -det_id $det_id";
     305    if ($stop) {
     306        $command .= ' -state stop';
     307    } else {
     308        $command .= ' -again';
     309    }
     310    $command .= " -dbname $dbname" if defined $dbname;
     311
     312    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
     313        run(command => $command, verbose => 1);
     314    unless ($success) {
     315        $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
     316        warn("Unable to perform dettool -updatedetrun: $error_code");
     317        exit($error_code);
     318    }
     319} else {
     320    print "skipping command: $command\n";
    300321}
    301322
    … …  
    317338    carp($msg);
    318339    if ($det_id and $iter and not $no_update) {
    319         my $command = "$dettool -adddetrunsummary -det_id $det_id -iteration $iter -code $exit_code";
     340        my $command = "$dettool -adddetrunsummary";
     341        $command .= " -det_id $det_id";
     342        $command .= " -iteration $iter";
     343        $command .= " -code $exit_code";
    320344        $command .= " -dbname $dbname" if defined $dbname;
    321345###        system ($command);
    … …  
    323347    exit $exit_code;
    324348}
    325 
    326 
    327349
    328350# Retrieve the requested rejection limit, dying unless extant
    … …  
    342364}
    343365
     366sub STATS_value_for_flag
     367{
     368    my $STATS = shift;
     369    my $flag  = shift;
     370
     371    foreach my $entry (@$STATS) {
     372        if ($flag eq $entry->{flag}) {
     373            return $entry->{value};
     374        }
     375    }
     376    return 'NAN';
     377}
    344378
    345379__END__
Note: See TracChangeset for help on using the changeset viewer.