Changeset 26224 for branches/eam_branches/20091113/psastro/src/psastroZeroPoint.c

Timestamp:

Nov 20, 2009, 3:22:59 PM (17 years ago)

Author:

eugene

Message:

add bootstrap-resampling to get zero point error for edge

File:

• branches/eam_branches/20091113/psastro/src/psastroZeroPoint.c (modified) (12 diffs)

branches/eam_branches/20091113/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);
+        psastroZeroPointAnalysis (header, dMag, zeropt, recipe);
         psFree (dMag);
         dMag = NULL;
@@ -141 +139 @@
 }
 
-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 +146 @@
       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 +167 @@
         }
     } else {
-        stats = psastroStatsPercentile (dMag, edgeFraction);
+        stats = psastroStatsPercentile (dMag, recipe);
         if (!stats) {
             psError(PS_ERR_UNKNOWN, false, "failure to measure zero point by edge");
@@ -233 +234 @@
 // 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 +251 @@
     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 +258 @@
     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 +272 @@
     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 = sqrt(S2 / edgeSample - PS_SQR(Sum/edgeSample));
+    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 +314 @@
     psFree(histogram);
     psFree(cumulative);
+    psFree(subset);
+    psFree(rng);
     return stats;
 
@@ -332 +331 @@
 }
 
+
+// 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) { \