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_imfile.pl

    r13989 r14009  
    2222use IPC::Cmd 0.36 qw( can_run run );             # tools to run UNIX programs with control over I/O
    2323use PS::IPP::Metadata::Config;                   # tools to parse the psMetadataConfig files
    24 use PS::IPP::Metadata::List qw( parse_md_list ); # ?
     24use PS::IPP::Metadata::List qw( parse_md_list ); # tools to parse a metadata into a hash list
    2525use Statistics::Descriptive;                     # tools for calculating basic statistical quantities
    2626use File::Temp qw( tempfile );                   # tools to construct temp files
     27
    2728use PS::IPP::Config qw($PS_EXIT_SUCCESS         
    2829                       $PS_EXIT_UNKNOWN_ERROR
    … …  
    5657
    5758pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
    58 pod2usage(
    59           -msg => "Required options: --det_id --iteration --exp_tag --det_type --camera",
    60           -exitval => 3,
    61           ) unless defined $det_id
     59pod2usage( -msg => "Required options: --det_id --iteration --exp_tag --det_type --camera",
     60           -exitval => 3)
     61    unless defined $det_id
    6262    and defined $iter
    6363    and defined $exp_tag
    … …  
    6868$ipprc->define_camera($camera);
    6969
    70 # Recipes to use, as a function of the detrend type
    71 use constant RECIPES => {
    72     'bin1' => {         # We're creating a master --- already processed the input
    73         'bias'     => 'PPIMAGE_J1_RESID_B',     # Bias only
    74         'dark'     => 'PPIMAGE_J1_RESID_B',     # Dark only
    75         'shutter'  => 'PPIMAGE_J1_RESID_F',     # Shutter only
    76         'flat'     => 'PPIMAGE_J1_RESID_F',     # Flat-field only
    77         'domeflat' => 'PPIMAGE_J1_RESID_F',     # Flat-field only
    78         'skyflat'  => 'PPIMAGE_J1_RESID_F',     # Flat-field only
    79         'fringe'   => 'PPIMAGE_J1_RESID_R',     # Fringe only
    80     },
    81     'bin2' => {         # We're checking the master --- input is not already processed
    82         'bias'     => 'PPIMAGE_J2_RESID_B',     # Bias only
    83         'dark'     => 'PPIMAGE_J2_RESID_B',     # Dark only
    84         'shutter'  => 'PPIMAGE_J2_RESID_F',     # Shutter only
    85         'flat'     => 'PPIMAGE_J2_RESID_F',     # Flat-field only
    86         'domeflat' => 'PPIMAGE_J2_RESID_F',     # Flat-field only
    87         'skyflat'  => 'PPIMAGE_J2_RESID_F',     # Flat-field only
    88         'fringe'   => 'PPIMAGE_J2_RESID_R',     # Fringe only
    89     },
    90 };
     70# Recipes to use based on reduction class
     71$reduction = 'DETREND' unless defined $reduction;
     72
     73my $recipe1 = $ipprc->reduction($reduction, 'JPEG_BIN1_RESID' . uc($det_type); # Recipe to use
     74&my_die("Unrecognised detrend type: $det_type", $det_id, $iter, $PS_EXIT_PROG_ERROR) unless defined $recipe1;
     75
     76my $recipe2 = $ipprc->reduction($reduction, 'JPEG_BIN2_RESID' . uc($det_type); # Recipe to use
     77&my_die("Unrecognised detrend type: $det_type", $det_id, $iter, $PS_EXIT_PROG_ERROR) unless defined $recipe2;
     78
     79# values to extract from output metadata and the stats to calculate
     80# XXX -bg_mean_stdev should take rms of bg_mean_stdev if bg_mean_stdev != 0 (A)
     81# XXX -bg_mean_stdev should take stdev of bg_mean if bg_mean_stdev == 0     (B)
     82# XXX  (A) if imfile.Ncomp > 1, (B) if imfile.Ncomp == 1
     83my $STATS =
     84   [   
     85       #          KEYWORD                 STATISTIC          CHIPTOOL FLAG
     86       { name => "bg",             type => "mean",  flag => "-bg" },
     87       { name => "bg_mean_stdev",  type => "stdev", flag => "-bg_mean_stdev" },
     88       { name => "bg_stdev",       type => "rms",   flag => "-bg_stdev" },
     89       { name => "bin_stdev",      type => "rms",   flag => "-bin_stdev" },
     90       { name => "fringe_0",       type => "mean",  flag => "-fringe_0" },
     91       { name => "fringe_1",       type => "rms",   flag => "-fringe_1" },
     92       { name => "fringe_0",       type => "stdev", flag => "-fringe_2" },
     93       { name => "user_1",         type => "mean",  flag => "-user_1" }, # fringe residual
     94       { name => "user_2",         type => "rms",   flag => "-user_2" }, # fringe residual
     95       { name => "user_3",         type => "stdev", flag => "-user_1" }, # fringe residual
     96   ];
    9197
    9298# Look for programs we need
    … …  
    98104    exit($PS_EXIT_CONFIG_ERROR);
    99105}
    100 
    101 # Parser for psMetadataConfig files (used for output from dettool)
    102 my $mdcParser = PS::IPP::Metadata::Config->new;
    103106
    104107# Get list of imfile files
    … …  
    120123    # XXX report an error message if stdout_buf is empty
    121124
     125    # Parse the stdout buffer into a metadata
     126    my $mdcParser = PS::IPP::Metadata::Config->new;     
    122127    my $metadata = $mdcParser->parse(join "", @$stdout_buf) or
    123128        &my_die("Unable to parse metadata config doc", $det_id, $iter, $exp_tag, $PS_EXIT_PROG_ERROR);
    124129
     130    # parse the file info in the metadata
    125131    $files = parse_md_list($metadata) or
    126132        &my_die("Unable to parse metadata list", $det_id, $iter, $exp_tag, $PS_EXIT_PROG_ERROR);
     133
     134    # Parse the statistics on the residual image
     135    my $stats = PS::IPP::Metadata::Stats->new($STATS); # Stats parser
     136    $stats->parse($metadata) or &my_die("Unable to find all values in statistics output.", $det_id, $iter, $exp_tag, $class_id, $PS_EXIT_PROG_ERROR);
    127137}
    128138
    … …  
    136146my $logFile;
    137147unless ($no_op) {
    138     # XXX this will fail if the file exists (because of Nebulous rules)
    139148    $logFile = $ipprc->file_create_append( $logName );
    140149}
    141 
    142 #### build the JPEG images ####
    143150
    144151# XXX in debug mode, unlink = 0
    … …  
    152159close $list2File;
    153160
    154 # Recipes to use in processing
    155 my $recipe1 = RECIPES->{"bin1"}->{lc($det_type)};
    156 &my_die("Unrecognised detrend type: $det_type", $det_id, $iter, $exp_tag, $PS_EXIT_PROG_ERROR) unless defined $recipe1;
    157 
    158 my $recipe2 = RECIPES->{"bin2"}->{lc($det_type)};
    159 &my_die("Unrecognised detrend type: $det_type", $det_id, $iter, $exp_tag, $PS_EXIT_PROG_ERROR) unless defined $recipe2;
    160 
     161
     162# build the JPEG images
    161163unless ($no_op) {
    162164
    … …  
    206208
    207209# storage variables
    208 my @means;                      # Array of means
    209 my @fluxes;                     # Array of means / exptimes
    210 my @variances;                  # Array of variances
    211 my @binVariances;               # Array of binned variances
    212 my @meanStdevs;                 # Array of mean stdevs
    213 my @names;                      # Array of names (class_id)
    214 my @fringe_means;
    215 my @fringe_vars;
    216 my @dfringe_means;
    217 my @dfringe_vars;
    218210
    219211# load the arrays from the imfile output table
    220212foreach my $file (@$files) {
    221     push @means,        $file->{bg};                           # mean background counts
    222213    if ($file->{exp_time} > 0.0) {
    223214        push @fluxes,           $file->{bg} / $file->{exp_time};       # mean background counts / sec
    … …  
    225216        push @fluxes,           $file->{bg};
    226217    }
    227     push @meanStdevs,        $file->{bg_mean_stdev};                # stdev of the mean counts (for imfile components)
    228     push @variances,         $file->{bg_stdev}*$file->{bg_stdev};   # total variance from all imfile components
    229     push @binVariances,      $file->{bin_stdev}*$file->{bin_stdev}; # total variance of binned images over all imfile components
    230     push @names,             $file->{class_id};                     # name of the component
    231     push @fringe_means,      $file->{fringe_0};                     # fringe amplitude mean for imfile
    232     push @fringe_vars,       $file->{fringe_1}*$file->{fringe_1};   # fringe variance
    233     push @dfringe_means,     $file->{user_1};                       # fringe residual mean for imfile
    234     push @dfringe_vars,      $file->{user_2}*$file->{user_2};       # fringe residual variance
    235     # push @fringe_mean_stdev, $file->{fringe_2};                     # fringe amplitude stdev for imfile
    236218}
    237219
    … …  
    239221# it is VALID to reject on more than one criterion
    240222&my_die("Number of means and number of variances differ!", $det_id, $iter, $exp_tag, $PS_EXIT_PROG_ERROR) unless scalar @means == scalar @variances;
    241 for (my $i = 0; $i < scalar @means; $i++) {
    242     my $flux = $fluxes[$i];     # Flux for this imfile
    243     my $mean = $means[$i];      # Mean for this imfile
     223
     224# storage array
     225my @fluxes;
     226
     227for (my $i = 0; $i < scalar @$files; $i++) {
     228    my $file      = $files[$i];
     229    my $name      = $file->{class_id};
     230    my $mean      = $file->{bg};        # Mean for this imfile
     231    my $stdev     = $file->{bg_stdev}; # Stdev for this imfile
     232    my $meanStdev = $file->{bg_mean_stdev}; # Stdev of Means for this imfile
     233    my $binStdev  = $file->{bin_stdev}; # Binned Stdev for this imfile
     234
     235    # calculate and save the fluxes
     236    my $flux      = $mean / $file->{exp_time};
     237    if ($file->{exp_time} == 0.0) {
     238        push @fluxes, $mean;
     239    }
     240    push @fluxes, $flux;
     241
    244242    $mean -= $expected;
    245     my $stdev = sqrt($variances[$i]);   # Stdev for this imfile
    246     my $binStdev = sqrt($binVariances[$i]);     # Stdev for this imfile
    247     my $name = $names[$i];
    248243
    249244    last if $no_op;
    … …  
    301296    # component means is larger than the limit
    302297    if ($reject_imfile_meanstdev > 0) {
    303         if ($meanStdevs[$i] > $reject_imfile_meanstdev) {
     298        if ($meanStdevs > $reject_imfile_meanstdev) {
    304299            print $logFile "Rejecting exposure based on bad imfile mean stdev for $name: ";
    305300            $reject = 1;
    … …  
    307302            print $logFile "Imfile mean stdev for $name meets requirements: ";
    308303        }
    309         print $logFile "$meanStdevs[$i] vs $reject_imfile_meanstdev\n";
     304        print $logFile "$meanStdevs vs $reject_imfile_meanstdev\n";
    310305    } else {
    311306        print $logFile "No rejection on imfile mean stdev for $name\n";
    … …  
    356351}
    357352
    358 # calculate the exposure ensemble statistics
    359 my $meanStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for means
    360 $meanStats->add_data(@means);
     353# basic ensemble stats
     354my $mean               = &STATS_value_for_flag ($STATS, "-bg");
     355my $meanStdev          = &STATS_value_for_flag ($STATS, "-bg_mean_stdev");
     356my $stdev              = &STATS_value_for_flag ($STATS, "-bg_stdev");
     357my $binStdev           = &STATS_value_for_flag ($STATS, "-bin_stdev");
     358my $fringe_mean        = &STATS_value_for_flag ($STATS, "-fringe_0");
     359my $fringe_err         = &STATS_value_for_flag ($STATS, "-fringe_1");
     360my $fringe_mean_stdev  = &STATS_value_for_flag ($STATS, "-fringe_2");
     361my $dfringe_mean       = &STATS_value_for_flag ($STATS, "-user_1");
     362my $dfringe_err        = &STATS_value_for_flag ($STATS, "-user_2");
     363my $dfringe_mean_stdev = &STATS_value_for_flag ($STATS, "-user_3");
     364
     365# other stats (flux depends on bg and exp_time)
    361366my $fluxStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for means
    362367$fluxStats->add_data(@fluxes);
    363 my $varianceStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for variances
    364 $varianceStats->add_data(@variances);
    365 my $binVarianceStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for variances
    366 $binVarianceStats->add_data(@binVariances);
    367 
    368 # background stats
    369 my $mean = $meanStats->mean();  # Mean of the imfile means
    370368my $flux = $fluxStats->mean();  # Mean of the imfile means
    371 my $meanStdev = $meanStats->standard_deviation(); # Stdev of the imfile means
    372 if (not defined $meanStdev) {
    373     # this is the case for Nimfile == 1
    374     $meanStdev = 0;
    375 }
    376 my $stdev = sqrt($varianceStats->mean()); # Root-Mean-Square of the imfile stdevs (root mean of variances)
    377 my $binStdev = sqrt($binVarianceStats->mean()); # Root-Mean-Square of the imfile stdevs (root mean of variances)
     369
     370# other stats
    378371my $exp_sn = 0.0;
    379372if ($stdev > 0) { $exp_sn = $mean / $stdev; }
    380 
    381 # prepare fringe amplitude stats
    382 my $fringeMeanStats = Statistics::Descriptive::Sparse->new();   # Statistics calculator for means
    383 $fringeMeanStats->add_data(@fringe_means);
    384 my $fringeVarStats = Statistics::Descriptive::Sparse->new();    # Statistics calculator for means
    385 $fringeVarStats->add_data(@fringe_vars);
    386 
    387 # fringe amplitude stats
    388 my $fringe_mean = $fringeMeanStats->mean();
    389 my $fringe_err = sqrt($fringeVarStats->mean());
    390 my $fringe_mean_stdev = $fringeMeanStats->standard_deviation();
    391 if (not defined $fringe_mean_stdev) {
    392     $fringe_mean_stdev = 0;
    393 }
    394 
    395 # prepare fringe residual stast
    396 my $dfringeMeanStats = Statistics::Descriptive::Sparse->new();  # Statistics calculator for means
    397 $dfringeMeanStats->add_data(@dfringe_means);
    398 my $dfringeVarStats = Statistics::Descriptive::Sparse->new();   # Statistics calculator for means
    399 $dfringeVarStats->add_data(@dfringe_vars);
    400 
    401 # fringe amplitude stats
    402 my $dfringe_mean = $dfringeMeanStats->mean();
    403 my $dfringe_err = sqrt($dfringeVarStats->mean());
    404 my $dfringe_mean_stdev = $dfringeMeanStats->standard_deviation();
    405 if (not defined $dfringe_mean_stdev) {
    406     $dfringe_mean_stdev = 0;
    407 }
    408373
    409374## Reject based on the exposure ensemble stats
    … …  
    514479}
    515480
    516 # Add the result into the database
    517 my $bg            = $mean;
    518 my $bg_stdev      = $stdev;
    519 my $bg_mean_stdev = $meanStdev;
    520 my $bin_stdev     = $binStdev;
     481my $command = "$dettool -addresidexp -det_id $det_id -iteration $iter -exp_tag $exp_tag";
     482$command .= " -recip $recipe1,$recipe2 -path_base $outputRoot ";
     483$command .= ' -reject' if $reject;
     484$command .= " -dbname $dbname" if defined $dbname;
     485
     486# add in the elements from the selected stats above
     487foreach my $entry (@$STATS) {
     488    my $value = $entry->{value};
     489    my $flag = $entry->{flag};
     490    $command .= " $flag $value";
     491}
    521492
    522493unless ($no_update) {
    523     my $command = "$dettool -addresidexp -det_id $det_id -iteration $iter -exp_tag $exp_tag";
    524     $command .= " -recip $recipe1,$recipe2 -path_base $outputRoot ";
    525     $command .= " -bg $bg -bg_stdev $bg_stdev -bg_mean_stdev $bg_mean_stdev -bin_stdev $bin_stdev";
    526     $command .= " -fringe_0 $fringe_mean -fringe_1 $fringe_err -fringe_2 $fringe_mean_stdev";
    527     $command .= " -user_1 $dfringe_mean -user_2 $dfringe_err -user_3 $dfringe_mean_stdev";
    528     $command .= ' -reject' if $reject;
    529     $command .= " -dbname $dbname" if defined $dbname;
    530494    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    531495        run(command => $command, verbose => 1);
    … …  
    535499        exit($error_code);
    536500    }
     501} else {
     502    print "skipping command: $command\n";
    537503}
    538504
    … …  
    555521    carp($msg);
    556522    if ($det_id and $iter and $exp_tag and not $no_update) {
    557         my $command = "$dettool -addresidexp -det_id $det_id -iteration $iter -exp_tag $exp_tag -code $exit_code";
     523        my $command = "$dettool -addresidexp";
     524        $command .= " -det_id $det_id";
     525        $command .= " -iteration $iter";
     526        $command .= " -exp_tag $exp_tag";
     527        $command .= " -code $exit_code";
    558528        $command .= " -dbname $dbname" if defined $dbname;
    559529###        system ($command);
    … …  
    561531    exit $exit_code;
    562532}
    563 
    564533
    565534# Retrieve the requested rejection limit, dying if not extant
    … …  
    579548}
    580549
    581 
    582 
     550sub STATS_value_for_flag
     551{
     552    my $STATS = shift;
     553    my $flag  = shift;
     554
     555    foreach my $entry (@$STATS) {
     556        if ($flag eq $entry->{flag}) {
     557            return $entry->{value};
     558        }
     559    }
     560    return 'NAN';
     561}
    583562
    584563__END__
Note: See TracChangeset for help on using the changeset viewer.