/** @file  pmPSFtry.c
 *
 *  XXX: need description of file purpose
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.69 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-01-27 06:39:38 $
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

# include <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAMaskWeight.h"
#include "pmSpan.h"
#include "pmFootprint.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmTrend2D.h"
#include "pmPSF.h"
#include "pmModel.h"
#include "pmSource.h"
#include "pmSourceUtils.h"
#include "pmPSFtry.h"
#include "pmModelClass.h"
#include "pmModelUtils.h"
#include "pmSourceFitModel.h"
#include "pmSourcePhotometry.h"
#include "pmSourceVisual.h"

bool printTrendMap (pmTrend2D *trend) {

    if (!trend->map) return false;
    if (!trend->map->map) return false;

    for (int j = 0; j < trend->map->map->numRows; j++) {
	for (int i = 0; i < trend->map->map->numCols; i++) {
	    fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
	}
	fprintf (stderr, "\t\t\t");
	for (int i = 0; i < trend->map->map->numCols; i++) {
	    fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
	}
	fprintf (stderr, "\n");
    }
    return true;
}

bool psImageMapCleanup (psImageMap *map) {

    if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true;

    // find the weighted average of all pixels
    float Sum = 0.0;
    float Wt = 0.0;
    for (int j = 0; j < map->map->numRows; j++) {
	for (int i = 0; i < map->map->numCols; i++) {
	    if (!isfinite(map->map->data.F32[j][i])) continue;
	    Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
	    Wt += map->error->data.F32[j][i];
	}
    }

    float Mean = Sum / Wt;

    // do any of the pixels in the map need to be repaired?
    // XXX for now, we are just replacing bad pixels with the Mean
    for (int j = 0; j < map->map->numRows; j++) {
	for (int i = 0; i < map->map->numCols; i++) {
	    if (isfinite(map->map->data.F32[j][i]) && 
		(map->error->data.F32[j][i] > 0.0)) continue;
	    map->map->data.F32[j][i] = Mean;
	}
    }
    return true;
}

// ********  pmPSFtry functions  **************************************************
// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
// * fitted version of the model, the pmPSF fit to the fitted model parameters,
// * and the PSF fits to the source. It also includes the statistics from the
// * fits, both the individual sources, and the collection

// free a pmPSFtry structure
static void pmPSFtryFree (pmPSFtry *test)
{
    if (test == NULL) return;

    psFree (test->psf);
    psFree (test->sources);
    psFree (test->metric);
    psFree (test->metricErr);
    psFree (test->fitMag);
    psFree (test->mask);
    return;
}

// allocate a pmPSFtry based on the desired sources and the model (identified by name)
pmPSFtry *pmPSFtryAlloc (const psArray *sources, const pmPSFOptions *options)
{
    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);

    test->psf       = pmPSFAlloc (options);
    test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
    test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
    test->fitMag    = psVectorAlloc (sources->n, PS_TYPE_F32);
    test->mask      = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK);

    psVectorInit (test->mask,        0);
    psVectorInit (test->metric,    0.0);
    psVectorInit (test->metricErr, 0.0);
    psVectorInit (test->fitMag,    0.0);

    test->sources   = psArrayAlloc (sources->n);

    for (int i = 0; i < sources->n; i++) {
        test->sources->data[i] = pmSourceCopy (sources->data[i]);
    }

    return (test);
}

bool psMemCheckPSFtry(psPtr ptr)
{
    PS_ASSERT_PTR(ptr, false);
    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmPSFtryFree);
}


// build a pmPSFtry for the given model:
// - fit each source with the free-floating model
// - construct the pmPSF from the collection of models
// - fit each source with the PSF-parameter models
// - measure the pmPSF quality metric (dApResid)

// sources used in for pmPSFtry may be masked by the analysis
// mask values indicate the reason the source was rejected:

// generate a pmPSFtry with a copy of the test PSF sources
pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
{
    bool status;
    int Next = 0;
    int Npsf = 0;

    // validate the requested model name
    options->type = pmModelClassGetType (modelName);
    if (options->type == -1) {
        psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
        return NULL;
    }

    pmPSFtry *psfTry = pmPSFtryAlloc (sources, options);
    if (psfTry == NULL) {
        psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
        return NULL;
    }

    // maskVal is used to test for rejected pixels, and must include markVal
    maskVal |= markVal;

    // stage 1:  fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF
    psTimerStart ("psf.fit");
    for (int i = 0; i < psfTry->sources->n; i++) {

        pmSource *source = psfTry->sources->data[i];
	if (!source->moments) {
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
	    continue;
	}
	if (!source->moments->nPixels) {
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
	    continue;
	}

        source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
        if (source->modelEXT == NULL) {
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
            continue;
        }

        // set object mask to define valid pixels
        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);

        // fit model as EXT, not PSF
        status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal);

        // clear object mask to define valid pixels
        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));

        // exclude the poor fits
        if (!status) {
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
            continue;
        }
        Next ++;
    }
    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n);
    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);

    if (Next == 0) {
        psError(PS_ERR_UNKNOWN, true, "No sources with good extended fits from which to determine PSF.");
        psFree(psfTry);
        return NULL;
    }

    // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
    if (!pmPSFFromPSFtry (psfTry)) {
        psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
        psFree(psfTry);
        return NULL;
    }

    // stage 3: refit with fixed shape parameters
    psTimerStart ("psf.fit");
    for (int i = 0; i < psfTry->sources->n; i++) {

        pmSource *source = psfTry->sources->data[i];

        // masked for: bad model fit, outlier in parameters
        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) {
            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y);
            continue;
        }

        // set shape for this model based on PSF
        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
        if (source->modelPSF == NULL) {
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
            abort();
            continue;
        }
        source->modelPSF->radiusFit = options->radius;

        // set object mask to define valid pixels
        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);

        // fit the PSF model to the source
        status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal);

        // skip poor fits
        if (!status) {
            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
            continue;
        }

        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal);
        if (!status || isnan(source->apMag)) {
            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
            continue;
        }

        // clear object mask to define valid pixels
        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));

        psfTry->fitMag->data.F32[i] = source->psfMag;
        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
        psfTry->metricErr->data.F32[i] = source->errMag;

        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
        Npsf ++;
    }
    psfTry->psf->nPSFstars = Npsf;

    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n);
    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);

    if (Npsf == 0) {
        psError(PS_ERR_UNKNOWN, true, "No sources with good PSF fits after model is built.");
        psFree(psfTry);
        return NULL;
    }

    // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0])
    // this should be linear for Poisson errors and quadratic for constant sky errors
    psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK);

    // generate the x and y vectors, and mask missing models
    for (int i = 0; i < psfTry->sources->n; i++) {
        pmSource *source = psfTry->sources->data[i];
        if (source->modelPSF == NULL) {
            flux->data.F32[i] = 0.0;
            chisq->data.F32[i] = 0.0;
            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
        } else {
            flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
            chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
        }
    }

    // use 3hi/3lo sigma clipping on the chisq fit
    psStats *stats = options->stats;

    // linear clipped fit of chisq trend vs flux
    bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, mask, 0xff, chisq, NULL, flux);
    psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
    psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev

    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
              psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));

    psFree(flux);
    psFree(mask);
    psFree(chisq);

    if (!result) {
        psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
        psFree(psfTry);
        return NULL;
    }

    for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i,
                  psfTry->psf->ChiTrend->coeff[i]*pow(10000, i),
                  psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i));
    }

    // XXX this function wants aperture radius for pmSourcePhotometry
    if (!pmPSFtryMetric (psfTry, options)) {
        psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName);
        psFree (psfTry);
        return NULL;
    }

    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);

    return (psfTry);
}

bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options)
{
    PS_ASSERT_PTR_NON_NULL(psfTry, false);
    PS_ASSERT_PTR_NON_NULL(options, false);
    PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);

    float RADIUS = options->radius;

    // the measured (aperture - fit) magnitudes (dA == psfTry->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

    // r2rflux = radius^2 * ten(0.4*fitMag);
    psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);

    for (int i = 0; i < psfTry->sources->n; i++) {
        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
            continue;
        r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]);
    }

    // XXX test dump of aperture residual data
    if (psTraceGetLevel("psModules.objects") >= 5) {
        FILE *f = fopen ("apresid.dat", "w");
        for (int i = 0; i < psfTry->sources->n; i++) {
            int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);

            pmSource *source = psfTry->sources->data[i];
            float x = source->peak->x;
            float y = source->peak->y;

            fprintf (f, "%d  %d, %f %f %f  %f %f %f \n",
                     i, keep, x, y,
                     psfTry->fitMag->data.F32[i],
                     r2rflux->data.F32[i],
                     psfTry->metric->data.F32[i],
                     psfTry->metricErr->data.F32[i]);
        }
        fclose (f);
    }

    // This analysis of the apResid statistics is only approximate.  The fitted magnitudes
    // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a
    // function of magnitude.  We re-measure the apResid statistics later in psphot using the
    // linear, constant-error fitting.  Do not reject outliers with excessive vigor here.

    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
    // linear clipped fit of ApResid to r2rflux
    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS)

    // XXX replace this with a psVectorStats call?  since we are not fitting the trend
    bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL,
                                              psfTry->metric, psfTry->metricErr, r2rflux);
    if (!result) {
        psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly");

        psFree(poly);
        psFree(r2rflux);

        return false;
    }
    psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0],
              psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);



    // XXX test dump of fitted model (dump when tracing?)
    if (psTraceGetLevel("psModules.objects") >= 4) {
        FILE *f = fopen ("resid.dat", "w");
        psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux);
        for (int i = 0; i < psfTry->sources->n; i++) {
            int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);

            pmSource *source = psfTry->sources->data[i];
            float x = source->peak->x;
            float y = source->peak->y;

            fprintf (f, "%d  %d, %f %f %f  %f %f %f  %f\n",
                     i, keep, x, y,
                     psfTry->fitMag->data.F32[i],
                     r2rflux->data.F32[i],
                     psfTry->metric->data.F32[i],
                     psfTry->metricErr->data.F32[i],
                     apfit->data.F32[i]);
        }
        fclose (f);
        psFree (apfit);
    }

    // XXX drop the skyBias value, or include above??
    psfTry->psf->skyBias  = poly->coeff[1];
    psfTry->psf->ApResid  = poly->coeff[0];
    psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat);

    psFree (r2rflux);
    psFree (poly);

    return true;
}

/*
  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
*/

/*****************************************************************************
pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
source->modelEXT entries.  The PSF ignores the first 4 (independent) model
parameters and constructs a polynomial fit to the remaining as a function of
image coordinate.
    input: psfTry with fitted source->modelEXT collection, pre-allocated psf
Note: some of the array entries may be NULL (failed fits); ignore them.
 *****************************************************************************/
bool pmPSFFromPSFtry (pmPSFtry *psfTry)
{
    PS_ASSERT_PTR_NON_NULL(psfTry, false);
    PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);

    pmPSF *psf = psfTry->psf;

    // construct the fit vectors from the collection of objects
    psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);

    // construct the x,y terms
    for (int i = 0; i < psfTry->sources->n; i++) {
        pmSource *source = psfTry->sources->data[i];
        if (source->modelEXT == NULL)
            continue;

        x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
        y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
    }

    if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) {
        psFree(x);
        psFree(y);
        psFree(z);
        return false;
    }

    // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
    // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters
    for (int i = 0; i < psf->params->n; i++) {
        switch (i) {
          case PM_PAR_SKY:
          case PM_PAR_I0:
          case PM_PAR_XPOS:
          case PM_PAR_YPOS:
          case PM_PAR_SXX:
          case PM_PAR_SYY:
          case PM_PAR_SXY:
            continue;
          default:
            break;
        }

        // select the per-object fitted data for this parameter
        for (int j = 0; j < psfTry->sources->n; j++) {
            pmSource *source = psfTry->sources->data[j];
            if (source->modelEXT == NULL) continue;
            z->data.F32[j] = source->modelEXT->params->data.F32[i];
        }

        psImageBinning *binning = psImageBinningAlloc();
        binning->nXruff = psf->trendNx;
        binning->nYruff = psf->trendNy;
        binning->nXfine = psf->fieldNx;
        binning->nYfine = psf->fieldNy;

        if (psf->psfTrendMode == PM_TREND_MAP) {
            psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
            psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
        }

        // free existing trend, re-alloc
        psFree (psf->params->data[i]);
        psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
        psFree (binning);

        // fit the collection of measured parameters to the PSF 2D model
        // the mask is carried from previous steps and updated with this operation
        // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
        if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) {
            psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
            psFree(x);
            psFree(y);
            psFree(z);
            return false;
        }
    }

    // test dump of star parameters vs position (compare with fitted values)
    if (psTraceGetLevel("psModules.objects") >= 4) {
        FILE *f = fopen ("params.dat", "w");

        for (int j = 0; j < psfTry->sources->n; j++) {
            pmSource *source = psfTry->sources->data[j];
            if (source == NULL) continue;
            if (source->modelEXT == NULL) continue;

            pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf);

            fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]);

            for (int i = 0; i < psf->params->n; i++) {
                if (psf->params->data[i] == NULL) continue;
                fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]);
            }
            fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter);

            psFree(modelPSF);
        }
        fclose (f);
    }

    psFree (x);
    psFree (y);
    psFree (z);
    return true;
}


bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {

    // we are doing a robust fit.  after each pass, we drop points which are more deviant than
    // three sigma.  mask is currently updated for each pass.

    // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
    // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
    // each source and fit this set of parameters.  These values are less tightly coupled, but
    // are still inter-related.  The fitted values do a good job of constraining the major axis
    // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
    // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
    // parameters, with the constraint that the minor axis must be greater than a minimum
    // threshold.

    // convert the measured source shape paramters to polarization terms
    psVector *e0   = psVectorAlloc (sources->n, PS_TYPE_F32);
    psVector *e1   = psVectorAlloc (sources->n, PS_TYPE_F32);
    psVector *e2   = psVectorAlloc (sources->n, PS_TYPE_F32);
    psVector *mag  = psVectorAlloc (sources->n, PS_TYPE_F32);

    for (int i = 0; i < sources->n; i++) {
	// skip any masked sources (failed to fit one of the model steps or get a magnitude)
	if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
	
        pmSource *source = sources->data[i];
        assert (source->modelEXT); // all unmasked sources should have modelEXT

        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);

        e0->data.F32[i] = pol.e0;
        e1->data.F32[i] = pol.e1;
        e2->data.F32[i] = pol.e2;

        float flux = source->modelEXT->params->data.F32[PM_PAR_I0];
        mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0;
    }

    if (psf->psfTrendMode == PM_TREND_MAP) {
        float scatterTotal = 0.0;
        float scatterTotalMin = FLT_MAX;
        int entryMin = -1;

        psVector *dz = NULL;
        psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK);

        // check the fit residuals and increase Nx,Ny until the error is minimized
        // pmPSFParamTrend increases the number along the longer of x or y
        for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {

	    // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
	    for (int i = 0; i < mask->n; i++) {
		mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
	    }
            if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
                break;
            }

            // store the resulting scatterTotal values and the scales, redo the best
            if (scatterTotal < scatterTotalMin) {
                scatterTotalMin = scatterTotal;
                entryMin = i;
            }
        }
        if (entryMin == -1) {
            psError (PS_ERR_UNKNOWN, true, "failed to find image map for shape params");
            return false;
        }

	// copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
	for (int i = 0; i < mask->n; i++) {
	    mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
	}
        if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
            psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
        }

        pmTrend2D *trend = psf->params->data[PM_PAR_E0];
        psf->trendNx = trend->map->map->numCols;
        psf->trendNy = trend->map->map->numRows;

	// copy mask back to srcMask
	for (int i = 0; i < mask->n; i++) {
	    srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
	}

        psFree (mask);
        psFree (dz);
    } else {

        // XXX iterate Nx, Ny based on scatter?
        // XXX we force the x & y order to be the same
        // modify the order to correspond to the actual number of matched stars:
        int order = PS_MAX (psf->trendNx, psf->trendNy);
        if ((sources->n < 15) && (order >= 3)) order = 2;
        if ((sources->n < 11) && (order >= 2)) order = 1;
        if ((sources->n <  8) && (order >= 1)) order = 0;
        if ((sources->n <  3)) {
            psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
            return false;
        }
        psf->trendNx = order;
        psf->trendNy = order;

        // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
        // This way, the parameters masked by one of the fits will be applied to the others
        for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {

	    psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
	    psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);

	    pmTrend2D *trend = NULL;
	    float mean, stdev;

            // XXX we are using the same stats structure on each pass: do we need to re-init it?
            bool status = true;

	    trend = psf->params->data[PM_PAR_E0];
            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
	    mean = psStatsGetValue (trend->stats, meanOption);
	    stdev = psStatsGetValue (trend->stats, stdevOption);
            psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
	    pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);

	    trend = psf->params->data[PM_PAR_E1];
            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
	    mean = psStatsGetValue (trend->stats, meanOption);
	    stdev = psStatsGetValue (trend->stats, stdevOption);
            psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
	    pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);

	    trend = psf->params->data[PM_PAR_E2];
            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
	    mean = psStatsGetValue (trend->stats, meanOption);
	    stdev = psStatsGetValue (trend->stats, stdevOption);
            psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
	    pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);

            if (!status) {
                psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
                return false;
            }
        }
    }

    // test dump of the psf parameters
    if (psTraceGetLevel("psModules.objects") >= 4) {
        FILE *f = fopen ("pol.dat", "w");
        fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
        for (int i = 0; i < e0->n; i++) {
            fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
                     x->data.F32[i], y->data.F32[i],
                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
                     pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
                     pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
                     srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
        }
        fclose (f);
    }

    psFree (e0);
    psFree (e1);
    psFree (e2);
    psFree (mag);
    return true;
}

// fit the shape variations as a psImageMap for the given scale factor
bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {

    int Nx, Ny;

    // set the map scale to match the aspect ratio : for a scale of 1, we guarantee
    // that we have a single cell
    if (psf->fieldNx > psf->fieldNy) {
        Nx = scale;
	float AR = psf->fieldNy / (float) psf->fieldNx;
	Ny = (int) (Nx * AR + 0.5);
        Ny = PS_MAX (1, Ny);
    } else {
        Ny = scale;
	float AR = psf->fieldNx / (float) psf->fieldNy;
	Nx = (int) (Ny * AR + 0.5);
        Nx = PS_MAX (1, Nx);
    }

    // do we have enough sources for this fine of a grid?
    if (x->n < 10*Nx*Ny) {
        return false;
    }

    // XXX check this against the exising type
    pmTrend2DMode psfTrendMode = PM_TREND_MAP;

    psImageBinning *binning = psImageBinningAlloc();
    binning->nXruff = Nx;
    binning->nYruff = Ny;
    binning->nXfine = psf->fieldNx;
    binning->nYfine = psf->fieldNy;
    psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);

    psFree (psf->params->data[PM_PAR_E0]);
    psFree (psf->params->data[PM_PAR_E1]);
    psFree (psf->params->data[PM_PAR_E2]);

    int nIter = psf->psfTrendStats->clipIter;
    psf->psfTrendStats->clipIter = 1;
    psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    psFree (binning);

    // if the map is 1x1 (a single value), we measure the resulting ensemble scatter

    // if the map is more finely sampled, divide the values into two sets: measure the fit from
    // one set and the scatter from the other set.
    psVector *x_fit = NULL;
    psVector *y_fit = NULL;
    psVector *x_tst = NULL;
    psVector *y_tst = NULL;

    psVector *e0obs_fit = NULL;
    psVector *e1obs_fit = NULL;
    psVector *e2obs_fit = NULL;
    psVector *e0obs_tst = NULL;
    psVector *e1obs_tst = NULL;
    psVector *e2obs_tst = NULL;

    if (scale == 1) {
	x_fit  = psMemIncrRefCounter (x);
	y_fit  = psMemIncrRefCounter (y);
	x_tst  = psMemIncrRefCounter (x);
	y_tst  = psMemIncrRefCounter (y);
	e0obs_fit = psMemIncrRefCounter (e0obs);
	e1obs_fit = psMemIncrRefCounter (e1obs);
	e2obs_fit = psMemIncrRefCounter (e2obs);
	e0obs_tst = psMemIncrRefCounter (e0obs);
	e1obs_tst = psMemIncrRefCounter (e1obs);
	e2obs_tst = psMemIncrRefCounter (e2obs);
    } else {
	x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
	for (int i = 0; i < e0obs_fit->n; i++) {
	    // e0obs->n ==  8 or 9:
	    // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
	    // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
	    x_fit->data.F32[i] = x->data.F32[2*i+0];  
	    x_tst->data.F32[i] = x->data.F32[2*i+1];  
	    y_fit->data.F32[i] = y->data.F32[2*i+0];  
	    y_tst->data.F32[i] = y->data.F32[2*i+1];  

	    e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0];  
	    e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1];  
	    e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
	    e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
	    e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
	    e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
	}
    }

    // the mask marks the values not used to calculate the ApTrend
    psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK);
    // copy mask values to fitMask as a starting point
    for (int i = 0; i < fitMask->n; i++) {
	fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    }

    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    // This way, the parameters masked by one of the fits will be applied to the others
    for (int i = 0; i < nIter; i++) {
        // XXX we are using the same stats structure on each pass: do we need to re-init it?
	psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
	psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);

	pmTrend2D *trend = NULL;
	float mean, stdev;

	// XXX we are using the same stats structure on each pass: do we need to re-init it?
	bool status = true;

	trend = psf->params->data[PM_PAR_E0];
	status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
	mean = psStatsGetValue (trend->stats, meanOption);
	stdev = psStatsGetValue (trend->stats, stdevOption);
	psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
	// printTrendMap (trend);
	psImageMapCleanup (trend->map);
	// printTrendMap (trend);
	pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);

	trend = psf->params->data[PM_PAR_E1];
	status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
	mean = psStatsGetValue (trend->stats, meanOption);
	stdev = psStatsGetValue (trend->stats, stdevOption);
	psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
	// printTrendMap (trend);
	psImageMapCleanup (trend->map);
	// printTrendMap (trend);
	pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);

	trend = psf->params->data[PM_PAR_E2];
	status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
	mean = psStatsGetValue (trend->stats, meanOption);
	stdev = psStatsGetValue (trend->stats, stdevOption);
	psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
	// printTrendMap (trend);
	psImageMapCleanup (trend->map);
	// printTrendMap (trend);
	pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
    }
    psf->psfTrendStats->clipIter = nIter; // restore default setting

    // construct the fitted values and the residuals
    psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst);
    psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst);
    psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst);

    psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit);
    psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit);
    psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit);

    // measure scatter for the unfitted points
    // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10);
    // psTraceSetLevel ("psLib.math.vectorClippedStats", 10);
    pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options));
    // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0);
    // psTraceSetLevel ("psLib.math.vectorClippedStats", 0);

    psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
    psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal);

    psFree (x_fit);
    psFree (y_fit);
    psFree (x_tst);
    psFree (y_tst);

    psFree (e0obs_fit);
    psFree (e1obs_fit);
    psFree (e2obs_fit);
    psFree (e0obs_tst);
    psFree (e1obs_tst);
    psFree (e2obs_tst);

    psFree (e0fit);
    psFree (e1fit);
    psFree (e2fit);

    psFree (e0res);
    psFree (e1res);
    psFree (e2res);

    // XXX copy fitMask values back to mask
    for (int i = 0; i < fitMask->n; i++) {
	mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    }
    psFree (fitMask);

    return true;
}

// calculate the scatter of the parameters
bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt)
{

    // psStats *stats = psStatsAlloc(stdevOpt);
    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV);

    // calculate the root-mean-square of E0, E1, E2
    float dEsquare = 0.0;
    psStatsInit (stats);
    psVectorStats (stats, e0res, NULL, mask, maskValue);
    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));

    psStatsInit (stats);
    psVectorStats (stats, e1res, NULL, mask, maskValue);
    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));

    psStatsInit (stats);
    psVectorStats (stats, e2res, NULL, mask, maskValue);
    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));

    *scatterTotal = sqrtf(dEsquare);

    psFree(stats);
    return true;
}

// calculate the minimum scatter of the parameters
bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res,
                            psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt)
{

    psStats *statsS = psStatsAlloc(stdevOpt);

    // measure the trend in bins with 10 values each; if < 10 total, use them all
    int nBin = PS_MAX (mag->n / nGroup, 1);

    // use mag to group parameters in sequence
    psVector *index = psVectorSortIndex (NULL, mag);

    // subset vectors for mag and dap values within the given range
    psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    psVector *mkSubset  = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);

    int n = 0;
    float min = INFINITY;               // Minimum error
    for (int i = 0; i < nBin; i++) {
        int j;
	int nValid = 0;
        for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
            int N = index->data.U32[n];
            dE0subset->data.F32[j] = e0res->data.F32[N];
            dE1subset->data.F32[j] = e1res->data.F32[N];
            dE2subset->data.F32[j] = e2res->data.F32[N];

            mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j]   = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
	    if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
        }
	if (nValid < 3) continue;

        dE0subset->n = j;
        dE1subset->n = j;
        dE2subset->n = j;
        mkSubset->n = j;

        // calculate the root-mean-square of E0, E1, E2
        float dEsquare = 0.0;
        psStatsInit (statsS);
        psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff);
        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));

        psStatsInit (statsS);
        psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff);
        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));

        psStatsInit (statsS);
        psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff);
        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));

        if (isfinite(dEsquare)) {
            float err = sqrtf(dEsquare);
            if (err < min) {
                min = err;
            }
        }
    }
    *errorFloor = min;

    psFree (dE0subset);
    psFree (dE1subset);
    psFree (dE2subset);
    psFree (mkSubset);

    psFree(index);

    psFree(statsS);

    return true;
}
