Index: trunk/psphot/src/pspsf.c
===================================================================
--- trunk/psphot/src/pspsf.c	(revision 4251)
+++ trunk/psphot/src/pspsf.c	(revision 4375)
@@ -34,25 +34,36 @@
     test->modelFLT    = psArrayAlloc (sources->n);
     test->modelPSF    = psArrayAlloc (sources->n);
-    test->metricStats = NULL;
     test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->fitMag      = psVectorAlloc (sources->n, PS_TYPE_F64);
     test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
+    test->ApResid     = 0;
+    test->dApResid    = 0;
+    test->skyBias     = 0;
 
     for (int i = 0; i < test->modelFLT->n; i++) {
+	test->mask->data.U8[i]  = 0;
 	test->modelFLT->data[i] = NULL;
 	test->modelPSF->data[i] = NULL;
-	test->mask->data.U8[i]  = 0;
+	test->metric->data.F64[i] = 0;
+	test->fitMag->data.F64[i] = 0;
     }	
     return (test);
 }
 
+// test->mask values indicate reason source was rejected:
+// 1: outlier in psf polynomial fit
+// 2: flt model failed to converge
+// 3: psf model failed to converge
+// 4: invalid source photometry
 pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) 
 {
-  
+    bool status;
+    float obsMag;
+    float fitMag;
     float x;
     float y;
-    bool status;
     int Nflt = 0;
     int Npsf = 0;
-    psVector *metricMask = psVectorAlloc (sources->n, PS_TYPE_U8);
+
     pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
 
@@ -66,11 +77,12 @@
 
 	// set temporary object mask and fit object
-	// fit model as FLT, not PSF (mask & skip poor fits)
+	// fit model as FLT, not PSF
 	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
 	status = pmSourceFitModel (source, model, false);
 	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
 
+	// exclude the poor fits
 	if (!status) {
-	  test->mask->data.U8[i] = 1;
+	  test->mask->data.U8[i] = 2;
 	  continue;
 	}
@@ -84,26 +96,16 @@
     // stage 2: construct a psf (pmPSF) from this collection of model fits
     pmPSFFromModels (test->psf, test->modelFLT, test->mask);
-    
-    // count valid sources
-    int Nkeep = 0;
-    for (int i = 0; i < sources->n; i++) {
-      if (test->mask->data.U8[i]) continue;
-      Nkeep++;
-    }
-    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates\n", Nkeep, sources->n);
 
     // stage 3: refit with fixed shape parameters
     psTimerStart ("fit");
     for (int i = 0; i < test->sources->n; i++) {
+	// masked for: bad model fit, outlier in parameters
+	if (test->mask->data.U8[i]) continue;
+
 	psSource *source = test->sources->data[i];
 	psModel  *modelFLT = test->modelFLT->data[i];
 
-	// masked for: bad model fit, outlier in parameters
-	if (test->mask->data.U8[i]) continue;
-	if (modelFLT == NULL) {
-	  psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
-	  continue;
-	}
-	psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
+	// set shape for this model based on PSF
+	psModel *modelPSF = psModelFromPSF (modelFLT, test->psf); 
 	x = source->peak->x;
 	y = source->peak->y;
@@ -111,13 +113,29 @@
 	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
 	status = pmSourceFitModel (source, modelPSF, true);
-	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
 
 	// skip poor fits
 	if (!status) {
-	  test->mask->data.U8[i] = 1;
-	  continue;
+	    test->mask->data.U8[i] = 3;
+	    goto next_source;
 	}
+
+	// otherwise, save the resulting model
 	test->modelPSF->data[i] = modelPSF;
+
+	// XXX : use a different aperture radius from the fit radius?
+	// XXX : use a different estimator for the local sky?
+	// XXX : pass 'source' as input?
+	if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
+	    test->mask->data.U8[i] = 4;
+	    goto next_source;
+	}	    
+
+	test->metric->data.F64[i] = obsMag - fitMag;
+	test->fitMag->data.F64[i] = fitMag;
 	Npsf ++;
+
+    next_source:
+	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+
     }
     psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
@@ -125,51 +143,8 @@
     DumpModelFits (test->modelPSF, "modelsPSF.dat");
 
-    // stage 4: measure PSF metric (Ap-Fit)
-    // should I be saving the aperture photometry for all sources?
-    // save the metric for all stars
-    for (int i = 0; i < test->sources->n; i++) {
-
-	float obsSum = 0;
-	float fitSum = 0;
-      
-	// is this metricMask redundant with test->mask ?
-	metricMask->data.U8[i] = 1;
-	test->metric->data.F64[i] = 0;
-
-	if (test->mask->data.U8[i]) continue;
-
-	psModel *model = test->modelPSF->data[i];
-	if (model == NULL) {
-	  psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
-	  continue;
-	}
-
-        psImage *image = ((psSource *)test->sources->data[i])->pixels;
-        psImage *mask  = ((psSource *)test->sources->data[i])->mask;
-
-	// this metric is Ap-Fit
-	float sky = model->params->data.F32[0];
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    for (int iy = 0; iy < image->numRows; iy++) {
-		if (mask->data.U8[iy][ix]) continue;
-		obsSum += image->data.F32[iy][ix] - sky;
-		fitSum += psModelEval (model, image, ix, iy) - sky;
-	    }
-	}
-	// fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
-	// keep the good metrics
-	if ((fitSum > 0) && (obsSum > 0)) {
-	    metricMask->data.U8[i] = 0;
-	    test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
-	} else {
-	  test->mask->data.U8[i] = 1;
-	}
-    }
-
-    // XXX use robust stats, not clipped stats
-    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    psVectorStats (stats, test->metric, NULL, metricMask, 1);
-    test->metricStats = stats;
-    psLogMsg ("psphot.pspsf", 3, "test model %s, metric: %f +/- %f\n", modelName, stats->clippedMean, stats->clippedStdev);
+    // XXX this function wants aperture radius from pmSourcePhotometry
+    pmPSFMetricModel (test, RADIUS);
+    psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
+	      modelName, test->ApResid, test->dApResid, test->skyBias);
     return (test);
 }
@@ -181,7 +156,7 @@
     // construct the fit vectors from the collection of objects
     // use the mask to ignore missing fits 
-    psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
     psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
 
@@ -208,17 +183,8 @@
 	    // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
 	}
-	// if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
 
 	psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
+	// psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
 	// psPolynomial2DDump (psf->params->data[i]);
-
-	// count valid sources
-	int Nkeep = 0;
-	for (int j = 0; j < mask->n; j++) {
-	  if (mask->data.U8[j]) continue;
-	  Nkeep++;
-	}
-	psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
-
     }
     return (true);
@@ -241,2 +207,104 @@
 }
 
+bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
+
+  float dBin;
+  int   nKeep, nSkip;
+  
+
+  // the measured (aperture - fit) magnitudes (dA == test->metric)
+  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+  //     dA = dAo + dsky/flux
+  //   where flux is the flux of the star
+  // we fit this trend to find the infinite flux aperture correction (dAo),
+  //   the nominal sky bias (dsky), and the error on dAo
+  // the values of dA are contaminated by stars with close neighbors in the aperture
+  //   we use an outlier rejection to avoid this bias
+
+  // rflux = ten(0.4*fitMag);
+  psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+  for (int i = 0; i < test->sources->n; i++) {
+    if (test->mask->data.U8[i]) continue;
+    rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
+  }
+  // psScalar *t1 = psScalar(0.4);
+  // psVector *v1 = psBinaryOp (NULL, test->fitMag, "*", t1);
+  // psVector *rflux = psUnaryOp (NULL, v1, "ten");
+  // psFree (t1);
+  // psFree (v1);
+
+  // find min and max of (1/flux):
+  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+  psVectorStats (stats, rflux, NULL, test->mask, 0xff);
+  
+  // build binned versions of rflux, metric
+  dBin = (stats->max - stats->min) / 10.0;
+  psVector *rfBin = psVectorCreate (stats->min, stats->max, dBin, PS_TYPE_F64);
+  psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
+  psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
+  psFree (stats);
+
+  psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
+
+  for (int i = 0; i < daBin->n; i++) {
+
+    psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+    tmp->n = 0;
+
+    // accumulate data within bin range
+    for (int j = 0; j < test->sources->n; j++) {
+      // masked for: bad model fit, outlier in parameters
+      if (test->mask->data.U8[j]) continue;
+    
+      // skip points with extreme dA values
+      if (fabs(test->metric->data.F64[j]) > 0.5) continue;
+
+      // skip points outside of this bin
+      if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue;
+      if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
+
+      tmp->data.F64[tmp->n] = test->metric->data.F64[j];
+      tmp->n ++;
+    }
+
+    // is this a valid point?
+    maskB->data.U8[i] = 0;
+    if (tmp->n < 2) {
+      maskB->data.U8[i] = 1;
+      continue;
+    } 
+
+    // dA values are contaminated with low outliers 
+    // measure statistics only on upper 50% of points
+    psVectorSort (tmp, tmp);
+    nKeep = 0.5*tmp->n;
+    nSkip = tmp->n - nKeep;
+
+    psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
+    for (int j = 0; j < tmp2->n; j++) {
+      tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
+    }
+
+    stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, tmp2, NULL, NULL, 0);
+    psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
+
+    daBin->data.F64[i] = stats->sampleMedian;
+  }
+
+  // linear fit to rfBin, daBin
+  psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+  poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
+
+  psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
+  psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
+
+  stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
+  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
+
+  test->ApResid = poly->coeff[0];
+  test->dApResid = stats->clippedStdev;
+  test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+  return true;
+}