/** @file  pmPSFtry.c
 *
 *  XXX: need description of file purpose
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-06-03 01:02:08 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

# include <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAMaskWeight.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmModel.h"
#include "pmSource.h"
#include "pmGrowthCurve.h"
#include "pmPSF.h"
#include "pmPSFtry.h"
#include "pmModelGroup.h"
#include "pmSourceFitModel.h"
#include "pmSourcePhotometry.h"

// ********  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->modelEXT);
    psFree (test->modelPSF);
    psFree (test->metric);
    psFree (test->fitMag);
    psFree (test->mask);
    return;
}

// allocate a pmPSFtry based on the desired sources and the model (identified by name)
pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors)
{

    pmModelType type;

    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));

    // XXX probably need to increment ref counter
    type           = pmModelSetType (modelName);
    test->psf      = pmPSFAlloc (type, poissonErrors);
    test->sources  = psMemIncrRefCounter(sources);
    test->modelEXT = psArrayAlloc (sources->n);
    test->modelPSF = psArrayAlloc (sources->n);
    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->metric->n = test->fitMag->n = test->mask->n = sources->n;
    test->modelEXT->n = test->modelPSF->n = sources->n;

    for (int i = 0; i < test->modelEXT->n; i++) {
        test->mask->data.U8[i]  = 0;
        test->modelEXT->data[i] = NULL;
        test->modelPSF->data[i] = NULL;
        test->metric->data.F64[i] = 0;
        test->fitMag->data.F64[i] = 0;
    }

    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
    return (test);
}

// 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:

pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors)
{
    bool status;
    float obsMag;
    float fitMag;
    float x;
    float y;
    int Next = 0;
    int Npsf = 0;

    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors);

    // stage 1:  fit an independent model (freeModel) to all sources
    psTimerStart ("fit");
    for (int i = 0; i < psfTry->sources->n; i++) {

        pmSource *source = psfTry->sources->data[i];
        pmModel  *model  = pmSourceModelGuess (source, psfTry->psf->type);
        x = source->peak->x;
        y = source->peak->y;

        // set temporary object mask and fit object
        // fit model as EXT, not PSF

        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
        status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT);
        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", NOT_U8(PM_MASK_MARK));

        // exclude the poor fits
        if (!status) {
            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
            psFree (model);
            continue;
        }
        psfTry->modelEXT->data[i] = model;
        Next ++;
    }
    psLogMsg ("psphot.psftry", 4, "fit ext:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (EXT)\n", Next, sources->n);

    // stage 2: construct a psf (pmPSF) from this collection of model fits
    pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask);

    // stage 3: refit with fixed shape parameters
    psTimerStart ("fit");
    for (int i = 0; i < psfTry->sources->n; i++) {
        // masked for: bad model fit, outlier in parameters
        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
            continue;

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

        // set shape for this model based on PSF
        pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf);
        x = source->peak->x;
        y = source->peak->y;

        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
        status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF);

        // skip poor fits
        if (!status) {
            psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
            psFree (modelPSF);
            goto next_source;
        }

        // otherwise, save the resulting model
        psfTry->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 (!pmSourcePhotometryModel (&fitMag, modelPSF)) {
            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
            goto next_source;
        }
        if (!pmSourcePhotometryAper (&obsMag, modelPSF, source->pixels, source->mask)) {
            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
            goto next_source;
        }

        psfTry->metric->data.F64[i] = obsMag - fitMag;
        psfTry->fitMag->data.F64[i] = fitMag;
        Npsf ++;

next_source:
        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", NOT_U8(PM_MASK_MARK));

    }
    psfTry->psf->nPSFstars = Npsf;

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

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

    // write sources with models first
    for (int i = 0; i < psfTry->sources->n; i++) {
        pmModel *model = psfTry->modelPSF->data[i];
        if (model == NULL) {
            flux->data.F64[i] = 0.0;
            chisq->data.F64[i] = 0.0;
            mask->data.U8[i] = 1;
        } else {
            flux->data.F64[i] = model->params->data.F32[1];
            chisq->data.F64[i] = model[0].chisq/model[0].nDOF;
            mask->data.U8[i] = 0;
        }
    }
    // use 3hi/3lo sigma clipping on the r2rflux vs metric fit
    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    stats->clipSigma = 3.0;
    stats->clipIter = 3;
    // linear clipped fit of ApResid to r2rflux
    psVectorClipFitPolynomial1D (psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux);
    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));
    }
    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    psFree (stats);
    psFree (flux);
    psFree (mask);
    psFree (chisq);

    // XXX this function wants aperture radius for pmSourcePhotometry
    pmPSFtryMetric (psfTry, RADIUS);
    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, float 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_F64);
    r2rflux->n = psfTry->sources->n;

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

    // use 3hi/1lo sigma clipping on the r2rflux vs metric fit
    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    stats->min = 1.0;
    stats->max = 3.0;
    stats->clipIter = 3;

    // 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 = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, r2rflux);
    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);

    pmPSF_MaskApTrend (psfTry->psf->ApTrend, PM_PSF_SKYBIAS);
    psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0];
    psfTry->psf->ApTrend->coeff[0][0][1][0] = 0;

    psfTry->psf->skyBias  = poly->coeff[1];
    psfTry->psf->ApResid  = poly->coeff[0];
    psfTry->psf->dApResid = stats->sampleStdev;

    psFree (r2rflux);
    psFree (poly);
    psFree (stats);

    return true;
}

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

