IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 8, 2007, 12:52:10 PM (19 years ago)
Author:
eugene
Message:

rejection is now a recipe; added rejection by flux

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippScripts/scripts/detrend_reject_imfile.pl

    r13698 r13722  
    11#!/usr/bin/env perl
     2
     3# this program has two jobs:
     4# 1: for the given exp_tag, generate the binned & mosaiced JPEG images
     5# 2: examine the collection of per-imfile statistics and reject.  At the moment,
     6#    this program (and the database) only allows rejection at the exposure level,
     7#    not at the imfile level.
    28
    39use warnings;
    … …  
    1319$VERSION = '0.01';
    1420
    15 use IPC::Cmd 0.36 qw( can_run run );
    16 use PS::IPP::Metadata::Config;
    17 use PS::IPP::Metadata::List qw( parse_md_list );
    18 use Statistics::Descriptive;
    19 use File::Temp qw( tempfile );
    20 
    21 use PS::IPP::Config qw($PS_EXIT_SUCCESS
     21use IPC::Cmd 0.36 qw( can_run run );             # tools to run UNIX programs with control over I/O
     22use PS::IPP::Metadata::Config;                   # tools to parse the psMetadataConfig files
     23use PS::IPP::Metadata::List qw( parse_md_list ); # ?
     24use Statistics::Descriptive;                     # tools for calculating basic statistical quantities
     25use File::Temp qw( tempfile );                   # tools to construct temp files
     26use PS::IPP::Config qw($PS_EXIT_SUCCESS         
    2227                       $PS_EXIT_UNKNOWN_ERROR
    2328                       $PS_EXIT_SYS_ERROR
    … …  
    2732                       $PS_EXIT_TIMEOUT_ERROR
    2833                       caturi
    29                        );
     34                       );                        # tools to parse the IPP configuration information
    3035my $ipprc = PS::IPP::Config->new(); # IPP configuration
    3136
    32 use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt );
    33 use Pod::Usage qw( pod2usage );
    34 
     37use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt ); # option parsing
     38use Pod::Usage qw( pod2usage );
     39
     40# parse the command-line options
    3541my ($det_id, $iter, $exp_tag, $det_type, $camera, $filter, $reject, $dbname, $workdir, $no_update, $no_op);
    3642GetOptions(
    … …  
    5864    and defined $camera;
    5965
     66# load IPP config information for the specified camera
    6067$ipprc->define_camera($camera);
    6168
    … …  
    9198}
    9299
    93 my $mdcParser = PS::IPP::Metadata::Config->new; # Parser for metadata config files
     100# Parser for psMetadataConfig files (used for output from dettool)
     101my $mdcParser = PS::IPP::Metadata::Config->new;
    94102
    95103# Get list of imfile files
    96104my $files;                      # Array of imfile files
    97105{
    98     my $command = "$dettool -residimfile -det_id $det_id -iteration $iter -exp_tag $exp_tag"; # Command to run
     106    # dettool command to select imfile data for this exp_tag
     107    my $command = "$dettool -residimfile";
     108    $command .= " -det_id $det_id";
     109    $command .= " -iteration $iter";
     110    $command .= " -exp_tag $exp_tag";
    99111    $command .= " -dbname $dbname" if defined $dbname;
    100112    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    … …  
    105117        exit($error_code);
    106118    }
    107 
    108119    # XXX report an error message if stdout_buf is empty
    109120
    … …  
    114125        &my_die("Unable to parse metadata list", $det_id, $iter, $exp_tag, $PS_EXIT_PROG_ERROR);
    115126}
    116 
    117 # XXX in debug mode, unlink = 0
    118 my ($list1File, $list1Name) = tempfile( "$exp_tag.detresid.$det_id.$iter.b1.list.XXXX", UNLINK => 1 );
    119 my ($list2File, $list2Name) = tempfile( "$exp_tag.detresid.$det_id.$iter.b2.list.XXXX", UNLINK => 1 );
    120 my @means;                      # Array of means
    121 my @variances;                  # Array of variances
    122 my @binVariances;               # Array of binned variances
    123 my @meanStdevs;                 # Array of mean stdevs
    124 my @names;                      # Array of names (class_id)
    125 foreach my $file (@$files) {
    126     print $list1File ($ipprc->filename( "PPIMAGE.BIN1", $file->{path_base}, $file->{class_id} ) . "\n");
    127     print $list2File ($ipprc->filename( "PPIMAGE.BIN2", $file->{path_base}, $file->{class_id} ) . "\n");
    128     push @means, $file->{bg};
    129     push @meanStdevs, $file->{bg_mean_stdev};
    130     ## calculate the root-mean-square of the bd_stdevs
    131     push @variances, $file->{bg_stdev}*$file->{bg_stdev};
    132     push @binVariances, $file->{bin_stdev}*$file->{bin_stdev};
    133     push @names, $file->{class_id};
    134 }
    135 close $list1File;
    136 close $list2File;
    137127
    138128# Output products
    … …  
    146136unless ($no_op) {
    147137    # XXX this will fail if the file exists (because of Nebulous rules)
    148     $logFile = $ipprc->file_create_open( $logName );
    149 }
     138    $logFile = $ipprc->file_create_append( $logName );
     139}
     140
     141#### build the JPEG images ####
     142
     143# XXX in debug mode, unlink = 0
     144my ($list1File, $list1Name) = tempfile( "$exp_tag.detresid.$det_id.$iter.b1.list.XXXX", UNLINK => 1 );
     145my ($list2File, $list2Name) = tempfile( "$exp_tag.detresid.$det_id.$iter.b2.list.XXXX", UNLINK => 1 );
     146foreach my $file (@$files) {
     147    print $list1File ($ipprc->filename( "PPIMAGE.BIN1", $file->{path_base}, $file->{class_id} ) . "\n");
     148    print $list2File ($ipprc->filename( "PPIMAGE.BIN2", $file->{path_base}, $file->{class_id} ) . "\n");
     149}
     150close $list1File;
     151close $list2File;
    150152
    151153# Recipes to use in processing
    … …  
    183185}
    184186
    185 my $expected = rejection_limit( 'EXPECTED', $det_type, $filter ); # Expected mean
    186 # Rejection thresholds
     187#### measure stats and reject if needed ####
     188
     189# Rejection thresholds & related data
     190my $expected                = rejection_limit( 'EXPECTED',         $det_type, $filter );
    187191my $reject_imfile_mean      = rejection_limit( 'IMFILE.MEAN',      $det_type, $filter );
     192my $reject_imfile_flux      = rejection_limit( 'IMFILE.FLUX',      $det_type, $filter );
    188193my $reject_imfile_stdev     = rejection_limit( 'IMFILE.STDEV',     $det_type, $filter );
    189194my $reject_imfile_meanstdev = rejection_limit( 'IMFILE.MEANSTDEV', $det_type, $filter );
    190195my $reject_imfile_sn        = rejection_limit( 'IMFILE.SN',        $det_type, $filter );
     196my $reject_imfile_bin_stdev = rejection_limit( 'IMFILE.BIN.STDEV', $det_type, $filter );
    191197my $reject_imfile_bin_sn    = rejection_limit( 'IMFILE.BIN.SN',    $det_type, $filter );
    192198my $reject_exp_mean         = rejection_limit( 'EXP.MEAN',         $det_type, $filter );
     199my $reject_exp_flux         = rejection_limit( 'EXP.FLUX',         $det_type, $filter );
    193200my $reject_exp_stdev        = rejection_limit( 'EXP.STDEV',        $det_type, $filter );
    194201my $reject_exp_meanstdev    = rejection_limit( 'EXP.MEANSTDEV',    $det_type, $filter );
    195202my $reject_exp_sn           = rejection_limit( 'EXP.SN',           $det_type, $filter );
     203my $reject_exp_bin_stdev    = rejection_limit( 'EXP.BIN.STDEV',    $det_type, $filter );
    196204my $reject_exp_bin_sn       = rejection_limit( 'EXP.BIN.SN',       $det_type, $filter );
     205
     206# storage variables
     207my @means;                      # Array of means
     208my @fluxes;                     # Array of means / exptimes
     209my @variances;                  # Array of variances
     210my @binVariances;               # Array of binned variances
     211my @meanStdevs;                 # Array of mean stdevs
     212my @names;                      # Array of names (class_id)
     213
     214# load the arrays from the imfile output table
     215foreach my $file (@$files) {
     216    push @means,        $file->{bg};                           # mean background counts
     217    if ($file->{exp_time} > 0.0) {
     218        push @fluxes,           $file->{bg} / $file->{exp_time};       # mean background counts / sec
     219    } else {
     220        push @fluxes,           $file->{bg};
     221    }
     222    push @meanStdevs,   $file->{bg_mean_stdev};                # stdev of the mean counts (for imfile components)
     223    push @variances,    $file->{bg_stdev}*$file->{bg_stdev};   # total variance from all imfile components
     224    push @binVariances, $file->{bin_stdev}*$file->{bin_stdev}; # total variance of binned images over all imfile components
     225    push @names,        $file->{class_id};                     # name of the component
     226}
    197227
    198228# Reject based on the stats of the imfiles
    … …  
    200230&my_die("Number of means and number of variances differ!", $det_id, $iter, $exp_tag, $PS_EXIT_PROG_ERROR) unless scalar @means == scalar @variances;
    201231for (my $i = 0; $i < scalar @means; $i++) {
     232    my $flux = $fluxes[$i];     # Flux for this imfile
    202233    my $mean = $means[$i];      # Mean for this imfile
    203234    $mean -= $expected;
    … …  
    208239    last if $no_op;
    209240
     241    # reject exposure if, for any imfiles, the mean residual counts
     242    # deviate from the expected value by more than N sigma, (sigma =
     243    # total pixel variance).  this test is sensible for images which
     244    # should have a predictable nominal residual mean count value (eg,
     245    # 0.0 for a bias).   in general, use with ADDITIVE components
    210246    if ($reject_imfile_mean > 0) {
    211247        if (abs($mean) > $reject_imfile_mean * $stdev) {
    … …  
    215251            print $logFile "Imfile mean for $name meets requirements: ";
    216252        }
    217         print $logFile "$mean vs $reject_imfile_mean\n";
     253        print $logFile "$mean vs $reject_imfile_mean * $stdev\n";
    218254    }  else {
    219255        print $logFile "No rejection on imfile mean for $name\n";
    220256    }
     257
     258    # reject exposure if, for any imfiles, the mean residual flux
     259    # deviates from the expected value by more than N sigma, (sigma =
     260    # total pixel variance).  this test is sensible for images which
     261    # should have a predictable nominal residual flux value (eg, 0.0
     262    # for a bias).  in general, use with ADDITIVE components
     263    if ($reject_imfile_flux > 0) {
     264        if (abs($flux) > $reject_imfile_flux) {
     265            print $logFile "Rejecting exposure based on bad imfile flux for $name: ";
     266            $reject = 1;
     267        } else {
     268            print $logFile "Imfile flux for $name meets requirements: ";
     269        }
     270        print $logFile "$flux vs $reject_imfile_flux\n";
     271    }  else {
     272        print $logFile "No rejection on imfile flux for $name\n";
     273    }
     274
     275    # reject exposure if, for any imfiles, the total pixel variance is
     276    # larger than the limit
    221277    if ($reject_imfile_stdev > 0) {
    222278        if ($stdev > $reject_imfile_stdev) {
    … …  
    231287        print $logFile "No rejection on imfile stdev for $name\n";
    232288    }
     289
     290    # reject exposure if, for any imfiles, the variance of the imfile
     291    # component means is larger than the limit
    233292    if ($reject_imfile_meanstdev > 0) {
    234293        if ($meanStdevs[$i] > $reject_imfile_meanstdev) {
    … …  
    242301        print $logFile "No rejection on imfile mean stdev for $name\n";
    243302    }
     303
     304    # reject exposure if, for any imfiles, the signal-to-noise (ie,
     305    # the mean counts / total pixel variance) of the imfile component
     306    # means are less than the limit.  this test is sensible for images
     307    # which have finite residual flux such as a flat-field image. 
    244308    if ($reject_imfile_sn > 0) {
    245309        if ($mean < $stdev * $reject_imfile_sn) {
    … …  
    253317        print $logFile "No rejection on imfile S/N for $name\n";
    254318    }
     319
     320    # reject exposure if, for any imfiles, the signal-to-noise of the
     321    # imfile component means, based on the stdev of the binned image
     322    # is less than the limit.  this test is sensible for images which
     323    # have finite residual flux such as a flat-field image.
     324    if ($reject_imfile_bin_stdev > 0) {
     325        if ($binStdev > $reject_imfile_bin_stdev) {
     326            print $logFile "Rejecting exposure based on bad imfile binned stdev for $name: ";
     327            $reject = 1;
     328        } else {
     329            print $logFile "Imfile binned stdev for $name meets requirements: ";
     330        }
     331        print $logFile "$binStdev vs $reject_imfile_bin_sn\n";
     332    } else {
     333        print $logFile "No rejection on imfile binned stdev for $name\n";
     334    }
    255335    if ($reject_imfile_bin_sn > 0) {
    256336        if ($mean < $binStdev * $reject_imfile_bin_sn) {
    … …  
    269349my $meanStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for means
    270350$meanStats->add_data(@means);
     351my $fluxStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for means
     352$fluxStats->add_data(@fluxes);
    271353my $varianceStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator for variances
    272354$varianceStats->add_data(@variances);
    … …  
    275357
    276358my $mean = $meanStats->mean();  # Mean of the imfile means
     359my $flux = $fluxStats->mean();  # Mean of the imfile means
    277360my $meanStdev = $meanStats->standard_deviation(); # Stdev of the imfile means
    278361if (not defined $meanStdev) {
    … …  
    290373    print $logFile "Exposure mean $mean, stdev $stdev, mean stdev $meanStdev, exp s/n: $exp_sn\n";
    291374
     375    # reject exposure if the ensemble mean residual counts deviate
     376    # from the expected value by more than N sigma, (sigma = total
     377    # pixel variance).  this test is sensible for images which should
     378    # have a predictable nominal residual mean count value (eg, 0.0
     379    # for a bias).  in general, use with ADDITIVE components
    292380    if ($reject_exp_mean > 0) {
    293381        if (abs($mean) > $reject_exp_mean * $stdev) {
    294             print $logFile "Rejecting exposure based on bad mean: ";
    295             $reject = 1;
    296         } else {
    297             print $logFile "Exposure mean meets requirements: ";
    298         }
    299         print $logFile "mean: $mean, stdev: $stdev (limit is: % $reject_exp_mean\n";
     382            print $logFile "Rejecting exposure based on bad mean counts: ";
     383            $reject = 1;
     384        } else {
     385            print $logFile "Exposure mean counts meets requirements: ";
     386        }
     387        print $logFile "mean: $mean, stdev: $stdev (limit is: $reject_exp_mean\n";
    300388    } else {
    301389        print $logFile "No rejection for exp mean\n";
    302390    }
    303     # reject if the exposure ensemble stdev is deviant
     391
     392    # reject exposure if, for any imfiles, the mean residual flux
     393    # deviates from the expected value by more than N sigma, (sigma =
     394    # total pixel variance).  this test is sensible for images which
     395    # should have a predictable nominal residual flux value (eg, 0.0
     396    # for a bias).  in general, use with ADDITIVE components
     397    if ($reject_exp_flux > 0) {
     398        if (abs($flux) > $reject_exp_flux * $stdev) {
     399            print $logFile "Rejecting exposure based on bad mean flux: ";
     400            $reject = 1;
     401        } else {
     402            print $logFile "Exposure mean flux meets requirements: ";
     403        }
     404        print $logFile "flux: $flux, stdev: $stdev (limit is: $reject_exp_flux\n";
     405    } else {
     406        print $logFile "No rejection for exp mean\n";
     407    }
     408
     409    # reject exposure if the total pixel variance is larger than the
     410    # limit
    304411    if ($reject_exp_stdev > 0) {
    305412        if ($stdev > $reject_exp_stdev) {
    … …  
    313420        print $logFile "No rejection for exp stdev\n";
    314421    }
    315     # reject if the exposure ensemble mean stdev is deviant
     422
     423    # reject exposure if the variance of the imfile means is larger
     424    # than the limit
    316425    if ($reject_exp_meanstdev > 0) {
    317426        if ($meanStdev > $reject_exp_meanstdev) {
    … …  
    325434        print $logFile "No rejection for exp mean stdev\n";
    326435    }
     436
    327437    # reject if the signal-to-noise is insufficient
    328438    if ($reject_exp_sn > 0) {
    … …  
    336446    } else {
    337447        print $logFile "No rejection for exp S/N\n";
     448    }
     449    # reject if the exposure ensemble stdev is deviant
     450    if ($reject_exp_bin_stdev > 0) {
     451        if ($binStdev > $reject_exp_bin_stdev) {
     452            print $logFile "Rejecting exposure based on bad binned stdev: ";
     453            $reject = 1;
     454        } else {
     455            print $logFile "Exposure binned stdev meets requirements: ";
     456        }
     457        print $logFile "$stdev vs $reject_exp_bin_stdev\n";
     458    } else {
     459        print $logFile "No rejection for exp binned stdev\n";
    338460    }
    339461    # reject if the signal-to-noise is insufficient
Note: See TracChangeset for help on using the changeset viewer.