Changeset 26259 for trunk/psastro/src/psastroZeroPoint.c

Timestamp:

Nov 22, 2009, 2:57:41 PM (16 years ago)

Author:

eugene

Message:

various fixes to psastro:

1) added bootstrap resampling to zero point error analysis
2) added iterative clump removal from refstars and rawstars
3) added unique reference match option
4) some improved visualizations
5) improved mosaic iterations

File:

• trunk/psastro/src/psastroZeroPoint.c (modified) (14 diffs)

trunk/psastro/src/psastroZeroPoint.c

@@ -30 +30 @@
     // recipe options
     bool byExposure = psMetadataLookupBool (&status, recipe, "ZERO.POINT.BY.EXPOSURE");
-    bool useMean    = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN");
-    float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION");
 
     // select the input data sources
@@ -86 +84 @@
                     // calculate dMag for the matched stars just for this readout (well, chip)
                     psMetadata *header = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
-                    psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean);
+                    psastroZeroPointAnalysis (header, dMag, zeropt, recipe);
                     psFree (dMag);
                     dMag = NULL;
@@ -96 +94 @@
     if (byExposure) {
         psMetadata *header = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
-        psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean);
+        if (!header) {
+          header = psMetadataAlloc ();
+          psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
+          psFree (header);
+        }
+        psastroZeroPointAnalysis (header, dMag, zeropt, recipe);
         psFree (dMag);
         dMag = NULL;
@@ -113 +116 @@
 
     // select the raw objects for this readout
-    psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
+    psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
     if (rawstars == NULL) return dMag;
 
     // select the raw objects for this readout
-    psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
+    psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
     if (refstars == NULL) return dMag;
 
@@ -141 +144 @@
 }
 
-bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean) {
+bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, psMetadata *recipe) {
 
     // XXX make this depend on the mode?
@@ -148 +151 @@
       return false;
     }
+
+    bool status;
+    bool useMean = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN");
 
     // the zero point analysis depends on the type of desired statistic.  For comparisons
@@ -166 +172 @@
         }
     } else {
-        stats = psastroStatsPercentile (dMag, edgeFraction);
+        stats = psastroStatsPercentile (dMag, recipe);
         if (!stats) {
             psError(PS_ERR_UNKNOWN, false, "failure to measure zero point by edge");
@@ -183 +189 @@
 }
 
-#define MAG_RESOLUTION 0.001
+#define MAG_RESOLUTION 0.0002
 
 // set the bin closest to the corresponding value.  if USE_END is +/- 1,
@@ -233 +239 @@
 // return results in stats->sampleMean, sampleStdev
 // XXX this is a misuse of psStats -- make our own structure?
-psStats *psastroStatsPercentile (psVector *myVector, float flimit) {
+psStats *psastroStatsPercentile (psVector *myVector, psMetadata *recipe) {
 
     // search for the 'blue' edge of the dMag distribution:
     // the distribution is not a normal population, but instead has a broad range with fairly hard edges.
     // construct a histogram and look for the 
+
+    bool status;
+    float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION");
+    int edgeSample = psMetadataLookupS32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE");
+    float edgeSampleFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE.FRACTION");
 
     // stats is first used to find the data range
@@ -245 +256 @@
     float min = NAN, max = NAN;         // Mimimum and maximum values
 
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);  
+
     // Get the minimum and maximum values
     if (!psVectorStats(stats, myVector, NULL, NULL, 0)) goto escape;
@@ -250 +263 @@
     max = stats->max;
     if (isnan(min) || isnan(max)) goto escape;
-    psTrace("psastro", 6, "Data min/max is (%.2f, %.2f)\n", min, max);
+    psTrace("psastro", 5, "Data min/max is (%.2f, %.2f)\n", min, max);
 
     // If all data points have the same value, then we set the appropriate members of stats and return.
@@ -264 +277 @@
     float binSize = MAG_RESOLUTION;
     long numBins = (max - min) / binSize; // Number of bins
-    psTrace("psastro", 6, "Numbins is %ld\n", numBins);
-    psTrace("psastro", 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
-
-    // Generate the histogram
+    psTrace("psastro", 5, "Numbins is %ld\n", numBins);
+    psTrace("psastro", 5, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
+
+    // allocate the histogram containers
     histogram = psHistogramAlloc(min, max, numBins);
-    if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
-      // if psVectorHistogram returns false, we have a programming error
-      psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for psastroZeroPointAnalysis\n");
-      psFree(histogram);
-      psFree(stats);
-      return NULL;
-    }
-    if (psTraceGetLevel("psastro") >= 8) {
-      PS_VECTOR_PRINT_F32(histogram->bounds);
-      PS_VECTOR_PRINT_F32(histogram->nums);
-    }
-
-    // Convert the specific histogram to a cumulative histogram
-    // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
     cumulative = psHistogramAlloc(min, max, numBins);
-    cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
-    for (long i = 1; i < histogram->nums->n; i++) {
-      cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
-      cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
-    }
-    if (psTraceGetLevel("psastro") >= 8) {
-      PS_VECTOR_PRINT_F32(cumulative->bounds);
-      PS_VECTOR_PRINT_F32(cumulative->nums);
-    }
-
-    // Find the bin which contains the first data point above the limit
-    long totalDataPoints = cumulative->nums->data.F32[numBins - 1];
-    psTrace("psastro", 6, "Total data points is %ld\n", totalDataPoints);
-
-    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
-    long bin;
-    PS_BIN_FOR_VALUE(bin, cumulative->nums, flimit * totalDataPoints, 0);
-    psTrace("psastro", 6, "The bin is %ld (%.2f to %.2f)\n", bin, cumulative->bounds->data.F32[bin], cumulative->bounds->data.F32[bin+1]);
-
-    // Linear interpolation to the limit value in bin units
-    float value;
-    PS_BIN_INTERPOLATE (value, cumulative->nums, cumulative->bounds, bin, totalDataPoints * flimit);
-    psTrace("psastro", 6, "limit value is %f\n", value);
-
-    stats->clippedMean = value;
-    stats->clippedStdev = 0.0; // XXX derive correct error value
+
+    // find the mean value:
+    stats->clippedMean = psastroStatsPercentileValue (histogram, cumulative, myVector, edgeFraction);
+
+    int nSubset = myVector->n * edgeSampleFraction;
+    psVector *subset = psVectorAlloc (nSubset, PS_TYPE_F32);
+
+    float Sum = 0.0;
+    float S2 = 0.0;
+    for (int i = 0; i < edgeSample; i++) {
+        
+        // generate the subset vector
+        for (long i = 0; i < nSubset; i++) {
+            double frnd = psRandomUniform(rng);
+            int entry = PS_MIN(myVector->n - 1, PS_MAX(0, myVector->n * frnd));
+
+            subset->data.F32[i] = myVector->data.F32[entry];
+        }
+
+        float value = psastroStatsPercentileValue (histogram, cumulative, subset, edgeFraction);
+
+        Sum += value;
+        S2 += value*value;
+    }
+    psTrace("psastro", 6, "subset stats: Sum: %f, S2: %f, Npts: %d\n", Sum, S2, edgeSample);
+
+    stats->clippedStdev = PS_MAX (sqrt(S2 / edgeSample - PS_SQR(Sum/edgeSample)), MAG_RESOLUTION);
+    psTrace("psastro", 5, "percentile stats %f +/- %f\n", stats->clippedMean, stats->clippedStdev);
+
     stats->results |= PS_STAT_CLIPPED_MEAN;
     stats->results |= PS_STAT_CLIPPED_STDEV;
@@ -317 +319 @@
     psFree(histogram);
     psFree(cumulative);
+    psFree(subset);
+    psFree(rng);
     return stats;
 
@@ -332 +336 @@
 }
 
+
+// measure the edge of the sample at flimit
+// return results in stats->sampleMean, sampleStdev
+// XXX this is a misuse of psStats -- make our own structure?
+float psastroStatsPercentileValue (psHistogram *histogram, psHistogram *cumulative, psVector *myVector, float flimit) {
+
+    // need to initialize the histogram on each pass
+    psVectorInit (histogram->nums, 0);
+    if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
+        // if psVectorHistogram returns false, we have a programming error
+        psAbort ("Unable to generate histogram for psastroZeroPointAnalysis");
+    }
+    if (psTraceGetLevel("psastro") >= 8) {
+        PS_VECTOR_PRINT_F32(histogram->bounds);
+        PS_VECTOR_PRINT_F32(histogram->nums);
+    }
+
+    // Convert the specific histogram to a cumulative histogram
+    // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
+    cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
+    for (long i = 1; i < histogram->nums->n; i++) {
+        cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
+        cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
+    }
+    if (psTraceGetLevel("psastro") >= 8) {
+        PS_VECTOR_PRINT_F32(cumulative->bounds);
+        PS_VECTOR_PRINT_F32(cumulative->nums);
+    }
+
+    // Find the bin which contains the first data point above the limit
+    long numBins = cumulative->nums->n;
+    long totalDataPoints = cumulative->nums->data.F32[numBins - 1];
+    psTrace("psastro", 6, "Total data points is %ld\n", totalDataPoints);
+
+    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
+    long bin;
+    PS_BIN_FOR_VALUE(bin, cumulative->nums, flimit * totalDataPoints, 0);
+    psTrace("psastro", 6, "The bin is %ld (%.4f to %.4f)\n", bin, cumulative->bounds->data.F32[bin], cumulative->bounds->data.F32[bin+1]);
+
+    // Linear interpolation to the limit value in bin units
+    float value;
+    PS_BIN_INTERPOLATE (value, cumulative->nums, cumulative->bounds, bin, totalDataPoints * flimit);
+    psTrace("psastro", 6, "limit value is %f\n", value);
+
+    return value;
+}
 
 # define ESCAPE(MSG) { \