Changeset 25352 for branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c

Timestamp:

Sep 13, 2009, 4:50:10 PM (17 years ago)

Author:

eugene

Message:

add function to report on systematic errors; measure sys errors for ap resid, adjust comments a bit

File:

• branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c (modified) (5 diffs)

branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c

@@ -245 +245 @@
         source->modelPSF->radiusFit = options->radius;
 
+        // XXXX use a different radius for the aperture magnitude than for the PSF fit?
+
         // set object mask to define valid pixels
         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
@@ -422 +424 @@
               psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);
 
-
+    float dSys = psVectorSystematicError (psfTry->metric, psfTry->metricErr, 0.1);
+    fprintf (stderr, "systematic error: %f\n", dSys);
+
+    int n = 0;
+    psVector *bright = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32);
+    psVector *brightErr = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32);
+    for (int i = 0; i < psfTry->metric->n; i++) {
+        if (!isfinite(psfTry->metric->data.F32[i])) continue;
+        if (!isfinite(psfTry->metricErr->data.F32[i])) continue;
+        if (psfTry->metricErr->data.F32[i] <= 0.0) continue;
+        if (psfTry->metricErr->data.F32[i] > 0.005) continue;
+        bright->data.F32[n] = psfTry->metric->data.F32[i];
+        brightErr->data.F32[n] = psfTry->metricErr->data.F32[i];
+        n++;
+    }
+    bright->n = brightErr->n = n;
+
+    float dSysBright = psVectorSystematicError (bright, brightErr, 0.1);
+    fprintf (stderr, "bright systematic error: %f\n", dSysBright);
 
     // XXX test dump of fitted model (dump when tracing?)
@@ -582 +602 @@
 }
 
-
 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
 
@@ -596 +615 @@
     // parameters, with the constraint that the minor axis must be greater than a minimum
     // threshold.
+
+    // XXX re-read the sextractor manual on handling 'infinitely thin' sources...
 
     // convert the measured source shape paramters to polarization terms
@@ -1069 +1090 @@
     return true;
 }
+
+float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) {
+
+    psAssert(residuals, "residuals cannot be NULL");
+    psAssert(errors, "errors cannot be NULL");
+    psAssert(residuals->n == errors->n, "residuals and errors must be the same length");
+
+    // given a vector of residuals and their formal errors, calculated the necessary systematic
+    // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction
+    // highest chi-square contributors (allowed outliers)
+
+    psVector *mask  = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK);
+    psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32);
+
+    // calculate the chisq vector:
+    int Ngood = 0;
+    for (int i = 0; i < residuals->n; i++) {
+        chisq->data.F32[i] = PS_MAX_F32;
+        if (!isfinite(residuals->data.F32[i])) continue;
+        if (!isfinite(errors->data.F32[i])) continue;
+        if (errors->data.F32[i] <= 0.0) continue;
+        chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]);
+        Ngood ++;
+    }
+
+    psVector *index = psVectorSortIndex(NULL, chisq);
+
+    // toss out the clipFraction highest chisq values
+    for (int i = 0; i < residuals->n; i++) {
+        int n = index->data.S32[i];
+        if (i < (1.0 - clipFraction)*Ngood) {
+            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0;
+        } else {
+            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1;
+        }
+    }
+
+    // Ndof ~= Ngood
+    // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof
+    // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0
+    
+    // use Newton-Raphson to solve for S2:
+
+    // use median sigma to calculate the initial guess for S2:
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, errors, NULL, mask, 1);
+    float errorMedian = stats->sampleMedian;
+    
+    float nPts = 0.0;
+    float res2mean = 0.0;
+    float ChiSq = 0.0;
+    for (int i = 0; i < residuals->n; i++) {
+        int n = index->data.S32[i];
+        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
+        res2mean += PS_SQR(residuals->data.F32[n]);
+        ChiSq += PS_SQR(residuals->data.F32[n] / errors->data.F32[n]);
+        nPts += 1.0;
+    }
+    res2mean /= nPts;
+    ChiSq /= nPts;
+    
+    float S2guess = res2mean - PS_SQR(errorMedian);
+
+    psLogMsg ("psModules", 3, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n", 
+              ChiSq, residuals->n, Ngood, nPts, S2guess);
+
+    for (int iter = 0; iter < 10; iter++) {
+
+        ChiSq = 0.0;
+        float dRdS = 0.0;
+        for (int i = 0; i < residuals->n; i++) {
+            int n = index->data.S32[i];
+            if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
+            float error2 = PS_SQR(errors->data.F32[n]) + S2guess;
+            ChiSq += PS_SQR(residuals->data.F32[n]) / error2;
+            dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2);
+        }
+        ChiSq /= nPts;
+        dRdS /= nPts;
+
+        // Note the sign on dS: dRdS above is -1 * dR/dS formally
+        float dS = (ChiSq - 1.0) / dRdS;
+        S2guess += dS;
+
+        psLogMsg ("psModules", 3, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess);
+    }
+
+    // free local allocations
+    return (sqrt(S2guess));
+}
+