IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9457


Ignore:
Timestamp:
Oct 10, 2006, 8:52:00 AM (20 years ago)
Author:
eugene
Message:

updating bg_stdev, bg_mean_stdev tests

Location:
trunk/ippScripts/scripts
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippScripts/scripts/detrend_reject_exp.pl

    r9453 r9457  
    7878my @expTags;                    # Array of exposure IDs
    7979my @means;                      # Array of means
    80 my @stdevs;                     # Array of stdevs
     80my @variances;                  # Array of variances
    8181my @meanStdevs;                 # Array of mean stdevs
    8282my @accept;                     # Array of accept flags
    … …  
    9191    push @expTags, $exposure->{exp_tag};
    9292    push @means, $exposure->{bg};
    93     push @stdevs, $exposure->{bg_stdev};
    94     push @meanStdevs, $exposure->{bg_mean_stdev};
     93    push @variances, ($exposure->{bg_stdev}*$exposure->{bg_stdev});
     94    push @meanStdevs, $exposure->{bg_mean_stdev}; ### XXX are we keeping this or stdev(mean)?
    9595    push @accept, $exposure->{accept};
    9696    push @include, $exposure->{include};
    … …  
    9898my $meanStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
    9999$meanStats->add_data(@means);
    100 my $stdevStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
    101 $stdevStats->add_data(@stdevs);
     100my $variancestats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
     101$variancestats->add_data(@variances);
    102102my $meanStdevStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
    103103$meanStdevStats->add_data(@meanStdevs);
    … …  
    111111    if not exists REJECT_MEAN_STDEV->{$det_type};
    112112
     113# rejections based on comparison with ensemble statistics
     114my $mean = $meanStats->mean();
     115my $meanStdev = $meanStats->standard_deviation();
     116if (not defined $meanStdev) { $meanStdev = 0; }
     117
     118my $var = $variances->mean();
     119my $varStdev = $variances->standard_deviation();
     120
    113121# Go through again to do rejection, and update the database for each exposure
    114122my $numChanges = 0;             # Number of exposures with changed status
    … …  
    117125    my $command = "$dettool -updateresidexp -det_id $det_id -iteration $iter -exp_tag $expTag"; # Command to run
    118126    my $reject = 0;             # Reject this exposure?
    119     my $not_reject = 0;         # Reject this exposure?
     127
    120128    if (not $accept[$i]) {
    121129        # Rejected this at an earlier stage
    122130        print "Rejecting $expTag based on earlier determination.\n";
    123131        $reject = 1;
    124     } elsif (defined REJECT_MEAN->{$det_type} and
    125              defined $meanStats->standard_deviation() and
    126              $meanStats->standard_deviation() > 0 and
    127              ($means[$i] - $meanStats->mean()) / $meanStats->standard_deviation() > REJECT_MEAN->{$det_type}) {
    128         print "Rejecting $expTag based on bad mean: " .
    129             (($means[$i] - $meanStats->mean()) / $meanStats->standard_deviation()) .
    130             " vs " . REJECT_MEAN->{$det_type} . "\n";
    131         $not_reject = 1;
    132     } elsif (defined REJECT_STDEV->{$det_type} and
    133              defined $stdevStats->standard_deviation() and
    134              $stdevStats->standard_deviation() > 0 and
    135              ($stdevs[$i] - $stdevStats->mean()) / $stdevStats->standard_deviation() >
    136              REJECT_STDEV->{$det_type}) {
    137         print "Rejecting $expTag based on bad stdev: " .
    138             (($stdevs[$i] - $stdevStats->mean()) / $stdevStats->standard_deviation()) .
    139             " vs " . REJECT_STDEV->{$det_type} . "\n";
    140         $not_reject = 1;
    141     } elsif (defined REJECT_MEAN_STDEV->{$det_type} and
    142              defined $meanStdevStats->standard_deviation() and
    143              $meanStdevStats->standard_deviation() > 0 and
    144              ($meanStdevs[$i] - $meanStdevStats->mean()) / $meanStdevStats->standard_deviation() >
    145              REJECT_MEAN_STDEV->{$det_type}) {
    146         print "Rejecting $expTag based on bad stdev: " .
    147             (($meanStdevs[$i] - $meanStdevStats->mean()) / $meanStdevStats->standard_deviation()) .
    148             " vs " . REJECT_MEAN_STDEV->{$det_type} . "\n";
    149         $not_reject = 1;
    150     }
    151    
     132        goto UPDATE;
     133    }
     134    if (REJECT_MEAN->{$det_type} && ($meanStdev > 0)) {
     135        my $nSigma = abs($means[$i] - $mean) / $meanStdev;
     136        if ($nSigma > REJECT_MEAN->{$det_type}) {
     137            print "Rejecting $expTag based on outlier mean value: " .
     138                "$means[$i] is $nSigma vs " . REJECT_MEAN->{$det_type} . "\n";
     139            $reject = 1;
     140            goto UPDATE;
     141        }
     142    } else {
     143        print "no rejection for exposure mean\n";
     144    }
     145    if (defined REJECT_STDEV->{$det_type} && ($varStdev >0)) {
     146        my $nSigma = abs($var[$i] - $var) / $varStdev;
     147        if ($nSigma > REJECT_STDEV->{$det_type}) {
     148            print "Rejecting $expTag based on outlier stdev: " .
     149                sqrt($variances[$i]) " is $nSigma vs " . REJECT_STDEV->{$det_type} . "\n";
     150            $reject = 1;
     151            goto UPDATE;
     152        }
     153    } else {
     154        print "no rejection for exposure stdev\n";
     155    }
     156   
     157  UPDATE:
    152158    if ($reject) {
    153159        $command .= ' -reject';
    … …  
    159165        $numChanges++;
    160166    }
    161    
     167
    162168    unless ($no_update) {
    163169        # Update
    … …  
    182188# Put the result into the database
    183189unless ($no_update) {
    184     my $command = "$dettool -adddetrunsummary -det_id $det_id -iteration $iter " .
    185         "-bg " . $meanStats->mean() . " -bg_stdev " . $stdevStats->mean() .
    186         " -bg_mean_stdev " . $meanStdevStats->mean();
     190    my $command = "$dettool -adddetrunsummary -det_id $det_id -iteration $iter" .
     191        " -bg " . $mean . " -bg_stdev " . sqrt($var) .
     192        " -bg_mean_stdev " . $meanStdev;
    187193    $command .= " -accept" if $master;
    188194   

  • trunk/ippScripts/scripts/detrend_reject_imfile.pl

    r9454 r9457  
    2121my ($det_id, $iter, $exp_tag, $det_type, $no_update, $reject);
    2222GetOptions(
    23     'det_id|d=s'        => \$det_id,
    24     'iteration=s'       => \$iter,
    25     'exp_tag|e=s'       => \$exp_tag,
    26     'det_type|t=s'      => \$det_type,
    27     'no-update'         => \$no_update,
    28     'reject'            => \$reject
    29 ) or pod2usage( 2 );
     23           'det_id|d=s'        => \$det_id,
     24           'iteration=s'       => \$iter,
     25           'exp_tag|e=s'       => \$exp_tag,
     26           'det_type|t=s'      => \$det_type,
     27           'no-update'         => \$no_update,
     28           'reject'            => \$reject
     29           ) or pod2usage( 2 );
    3030
    3131pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
    3232pod2usage(
    33     -msg => "Required options: --det_id --iteration --exp_tag --det_type",
    34     -exitval => 3,
    35 ) unless defined $det_id
     33          -msg => "Required options: --det_id --iteration --exp_tag --det_type",
     34          -exitval => 3,
     35          ) unless defined $det_id
    3636    and defined $iter
    3737    and defined $exp_tag
    … …  
    4242
    4343#### XXXXX these values must come from the config system, and may depend on filter!!!
     44# XXX it is valid to reject on more than one criterion
     45
    4446# The expected mean, as a function of detrend type
    4547use constant EXPECTED_MEAN => {
    … …  
    5052    };
    5153
    52 # Rejection threshold for the mean, in terms of the standard deviation
     54# Rejection threshold for the mean
    5355# This measures how close it is to what's expected
    54 # XXX it is valid to reject on more than one criterion
    55 use constant REJECT_COMPONENT_MEAN => {
     56use constant REJECT_IMFILE_MEAN => {
    5657    'bias' => 0,                # Should be fairly flat; some CRs
    5758    'dark' => 0,                # Lots of CRs
    … …  
    6162
    6263# Rejection threshold for the standard deviation, in ADUs
    63 # This measures how much variation there is in each imfile component
    64 use constant REJECT_COMPONENT_STDEV => {
     64# This measures how much variation there is in each imfile
     65use constant REJECT_IMFILE_STDEV => {
    6566    'bias' => 0,                # Should be fairly flat; some CRs
    6667    'dark' => 0,                # Lots of CRs
    … …  
    6970    };
    7071
    71 # Rejection threshold for the mean of the sample, in terms of the standard deviation of the sample
     72# Rejection threshold for the mean of the exposure, in terms of the standard deviation of the exposure
    7273# This measures how close it is to what's expected
    73 use constant REJECT_SAMPLE_MEAN => {
     74use constant REJECT_EXPOSURE_MEAN => {
    7475    'bias' => 0,                # Should be little variation between chips
    7576    'dark' => 0,                # Could be some glow on some chips
    … …  
    7879    };
    7980
    80 # Rejection threshold for the stdev of the sample, in ADUs
    81 # This measures how much variation there is across the imfile components
    82 use constant REJECT_SAMPLE_STDEV => {
     81# Rejection threshold for the stdev of the exposure, in ADUs
     82# This measures how much variation there is across the imfiles
     83use constant REJECT_EXPOSURE_STDEV => {
    8384    'bias' => 0,                # Should be little variation between chips
    8485    'dark' => 0,                # Could be some glow on some chips
    … …  
    8788    };
    8889
    89 # Rejection threshold for the stdev of the sample, in ADUs
    90 # This measures how much variation there is across the imfile components
    91 use constant REJECT_SAMPLE_MEAN_STDEV => {
     90# Rejection threshold for the stdev of the exposure, in ADUs
     91# This measures how much variation there is across the imfiles
     92use constant REJECT_EXPOSURE_MEAN_STDEV => {
    9293    'bias' => 0,                # Should be little variation between chips
    9394    'dark' => 0,                # Could be some glow on some chips
    … …  
    105106my $mdcParser = PS::IPP::Metadata::Config->new; # Parser for metadata config files
    106107
    107 # Get list of component files
    108 my $files;                      # Array of component files
     108# Get list of imfile files
     109my $files;                      # Array of imfile files
    109110{
    110111    my $command = "$dettool -residimfile -det_id $det_id -iteration $iter -exp_tag $exp_tag"; # Command to run
    … …  
    126127my $list2Name = $outputRoot . '.b2.list'; # Name for the input file list for binning 2
    127128my @means;                      # Array of means
    128 my @stdevs;                     # Array of stdevs
     129my @variances;                  # Array of variances
    129130open my $list1File, '>' . $list1Name;
    130131open my $list2File, '>' . $list2Name;
    … …  
    134135    push @means, $file->{bg};
    135136    ## calculate the root-mean-square of the bd_stdevs
    136     push @stdevs, $file->{bg_mean_stdev}*$file->{bg_mean_stdev};
     137    push @variances, $file->{bg_stdev}*$file->{bg_stdev};
    137138}
    138139close $list1File;
    … …  
    162163die "Unknown expected mean for detrend type $det_type\n"
    163164    if not exists EXPECTED_MEAN->{$det_type};
    164 die "Unknown component mean rejection level for detrend type $det_type\n"
    165     if not exists REJECT_COMPONENT_MEAN->{$det_type};
    166 die "Unknown component stdev rejection level for detrend type $det_type\n"
    167     if not exists REJECT_COMPONENT_STDEV->{$det_type};
    168 die "Unknown sample mean rejection level for detrend type $det_type\n"
    169     if not exists REJECT_SAMPLE_MEAN->{$det_type};
    170 die "Unknown sample stdev rejection level for detrend type $det_type\n"
    171     if not exists REJECT_SAMPLE_STDEV->{$det_type};
    172 
    173 # Reject based on the stats of the component imfile
    174 die "Number of means and number of stdevs differ!\n" if scalar @means != scalar @stdevs;
     165die "Unknown imfile mean rejection level for detrend type $det_type\n"
     166    if not exists REJECT_IMFILE_MEAN->{$det_type};
     167die "Unknown imfile stdev rejection level for detrend type $det_type\n"
     168    if not exists REJECT_IMFILE_STDEV->{$det_type};
     169die "Unknown exposure mean rejection level for detrend type $det_type\n"
     170    if not exists REJECT_EXPOSURE_MEAN->{$det_type};
     171die "Unknown exposure stdev rejection level for detrend type $det_type\n"
     172    if not exists REJECT_EXPOSURE_STDEV->{$det_type};
     173die "Unknown exposure mean stdev rejection level for detrend type $det_type\n"
     174    if not exists REJECT_EXPOSURE_MEAN_STDEV->{$det_type};
     175
     176# Reject based on the stats of the imfiles
     177# it is VALID to reject on more than one criterion
     178die "Number of means and number of variances differ!\n" if scalar @means != scalar @variances;
    175179for (my $i = 0; $i < scalar @means; $i++) {
    176     my $mean = $means[$i];      # Mean for this component
     180    my $mean = $means[$i];      # Mean for this imfile
    177181    $mean -= EXPECTED_MEAN->{$det_type} if defined EXPECTED_MEAN->{$det_type};
    178     my $stdev = sqrt($stdevs[$i]);      # Stdev for this component
    179 
    180     ## IT IS VALID to reject on more than one criterion
    181     if (REJECT_COMPONENT_MEAN->{$det_type} && ($stdev > 0) ) {
    182         if ($mean / $stdev > REJECT_COMPONENT_MEAN->{$det_type}) {
    183             print "Rejecting exposure based on bad component mean for component $i: " .
    184                 ($mean / $stdev) . " vs " . REJECT_COMPONENT_MEAN->{$det_type} . "\n";
     182    my $stdev = sqrt($variances[$i]);   # Stdev for this imfile
     183
     184    if (REJECT_IMFILE_MEAN->{$det_type}) {
     185        if (abs($mean) > REJECT_IMFILE_MEAN->{$det_type}) {
     186            print "Rejecting exposure based on bad imfile mean for imfile $i: " .
     187                $mean . " vs " . REJECT_IMFILE_MEAN->{$det_type} . "\n";
    185188            $reject = 1;
    186189            last;
    187190        }
    188     }
    189     if (REJECT_COMPONENT_STDEV->{$det_type}) {
    190         if ($stdev > REJECT_COMPONENT_STDEV->{$det_type}) {
    191             print "Rejecting exposure based on bad component stdev for component $i: " .
    192                 $stdev . " vs " . REJECT_COMPONENT_STDEV->{$det_type} . "\n";
     191    }  else {
     192        print "no rejection for imfile mean\n";
     193    }
     194    if (REJECT_IMFILE_STDEV->{$det_type}) {
     195        if ($stdev > REJECT_IMFILE_STDEV->{$det_type}) {
     196            print "Rejecting exposure based on bad imfile stdev for imfile $i: " .
     197                $stdev . " vs " . REJECT_IMFILE_STDEV->{$det_type} . "\n";
    193198            $reject = 1;
    194199            last;
    195200        }
     201    } else {
     202        print "no rejection for imfile stdev\n";
    196203    }
    197204}
    198205
    199 # calculate the imfile ensemble statistics
     206# calculate the exposure ensemble statistics
    200207my $meanStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for means
    201208$meanStats->add_data(@means);
    202 my $stdevStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for stdevs
    203 $stdevStats->add_data(@stdevs);
    204 my $mean = $meanStats->mean();  # Mean of the sample of means
    205 my $meanStdev = $meanStats->standard_deviation(); # Stdev of the sample of means
    206 if (not defined $stdev) {
    207     $stdev = 0;
    208 }
    209 my $stdev = sqrt($stdevStats->mean()); # Root-Mean-Square of the sample of stdevs
    210 print "calculating root-mean-square of the sample stdevs: $meanStdev\n";
    211 
    212 ## Reject based on the imfile ensemble stats
    213 # reject if the imfile ensemble
    214 if (REJECT_SAMPLE_MEAN->{$det_type} && ($stdev > 0)) {
    215     if ($mean / $stdev > REJECT_SAMPLE_MEAN->{$det_type}) {
     209my $varianceStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for variances
     210$varianceStats->add_data(@variances);
     211
     212my $mean = $meanStats->mean();  # Mean of the imfile means
     213my $meanStdev = $meanStats->standard_deviation(); # Stdev of the imfile means
     214if (not defined $meanStdev) {
     215    # this is the case for Nimfile == 1
     216    $meanStdev = 0;
     217}
     218my $stdev = sqrt($varianceStats->mean()); # Root-Mean-Square of the imfile stdevs (root mean of variances)
     219print "exposure mean $mean, stdev $stdev, mean stdev $meanStdev\n";
     220
     221## Reject based on the exposure ensemble stats
     222# reject if the exposure ensemble mean is deviant
     223if (REJECT_EXPOSURE_MEAN->{$det_type}) {
     224    if (abs($mean) > REJECT_EXPOSURE_MEAN->{$det_type}) {
    216225        print "Rejecting exposure based on bad mean: " . ($mean / $stdev) . " vs " .
    217             REJECT_SAMPLE_MEAN->{$det_type} . "\n";
     226            REJECT_EXPOSURE_MEAN->{$det_type} . "\n";
    218227        $reject = 1;
    219228    }
    220 }
    221 if (REJECT_SAMPLE_STDEV->{$det_type}) {
    222     if ($stdev > REJECT_SAMPLE_STDEV->{$det_type}) {
     229} else {
     230    print "no rejection for imfile mean\n";
     231}
     232# reject if the exposure ensemble stdev is deviant
     233if (REJECT_EXPOSURE_STDEV->{$det_type}) {
     234    if ($stdev > REJECT_EXPOSURE_STDEV->{$det_type}) {
    223235        print "Rejecting exposure based on bad mean stdev: " . $stdev . " vs " .
    224             REJECT_SAMPLE_STDEV->{$det_type} . "\n";
     236            REJECT_EXPOSURE_STDEV->{$det_type} . "\n";
    225237        $reject = 1;
    226238    }
    227 }
    228 if (REJECT_SAMPLE_MEAN_STDEV->{$det_type}) {
    229     if ($meanStdev > REJECT_SAMPLE_MEAN_STDEV->{$det_type}) {
     239} else {
     240    print "no rejection for imfile stdev\n";
     241}
     242# reject if the exposure ensemble mean stdev is deviant
     243if (REJECT_EXPOSURE_MEAN_STDEV->{$det_type}) {
     244    if ($meanStdev > REJECT_EXPOSURE_MEAN_STDEV->{$det_type}) {
    230245        print "Rejecting exposure based on bad mean stdev: " . $meanStdev . " vs " .
    231             REJECT_SAMPLE_MEAN_STDEV->{$det_type} . "\n";
     246            REJECT_EXPOSURE_MEAN_STDEV->{$det_type} . "\n";
    232247        $reject = 1;
    233248    }
     249} else {
     250    print "no rejection for imfile mean stdev\n";
    234251}
    235252

  • trunk/ippScripts/scripts/detrend_resid.pl

    r9446 r9457  
    115115
    116116# Add the processed file to the database
     117# XXX I think this has the names "bg_stdev" and "bg_mean_stdev" exchanged
     118#     bg_stdev : standard deviation of the background
     119#     bg_mean_stdev : standard deviation of the background means
    117120$outputName = File::Spec->abs2rel ($outputName, $ipprc->workdir() );
    118121$bin1Name = File::Spec->abs2rel( $bin1Name, $ipprc->workdir() );
    … …  
    123126        "-b2_uri $bin2Name"; # Command to run dettool
    124127    $command .= " -bg " . $stats->bg_mean();
     128
     129    # XXX note bg_stdev <--> bg_mean_stdev
    125130    if (defined($stats->bg_stdev())) {
    126         $command .= " -bg_stdev " . $stats->bg_stdev();
     131        $command .= " -bg_mean_stdev " . $stats->bg_stdev();
    127132    } else {
    128133        # May be undefined if there is only a single imfile
    129134        $command .= ' -bg_stdev 0';
    130135    }
    131     $command .= " -bg_mean_stdev " . $stats->bg_mean_stdev();
     136
     137    # XXX note bg_stdev <--> bg_mean_stdev
     138    $command .= " -bg_stdev " . $stats->bg_mean_stdev();
    132139
    133140    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
Note: See TracChangeset for help on using the changeset viewer.