/** @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"

// 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)
{
    // 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
    if (!pmPSFtryFitEXT(psfTry, options, maskVal, markVal)) {
        psError(PS_ERR_UNKNOWN, false, "failed to fit EXT models to sources for psf model");
        psFree(psfTry);
        return NULL;
    }

    // hard limit on minimum number of stars
    if ((sources->n <  3)) {
        psError (PS_ERR_UNKNOWN, true, "failed to determine PSF parameters");
        return NULL;
    }

    // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different 
    // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a 
    // user-specified MAX order, which we should respect, regardless of the mode

    // set the max order (0 = constant) which the number of psf stars can support:
    // rule of thumb: require 3 stars per 'cell' (order+1)^2
    int MaxOrderForStars = 0;
    if (sources->n >= 12) MaxOrderForStars = 1; // 4 cells
    if (sources->n >= 27) MaxOrderForStars = 2; // 9 cells
    if (sources->n >= 48) MaxOrderForStars = 3; // 16 cells
    if (sources->n >  75) MaxOrderForStars = 4; // 25 cells

    int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy);
    int orderMin = 0;
    if (options->psfTrendMode == PM_TREND_MAP) {
	MaxOrderForStars ++;
	orderMin ++;
    }
    orderMax = PS_MIN (orderMax, MaxOrderForStars);

    // save the raw source mask (generated by pmPSFtryFitEXT)
    psVector *srcMask = psVectorCopy (NULL, psfTry->mask, PS_TYPE_VECTOR_MASK);

    // we will save the PSF with the best fit (min systematic error)
    pmPSF *minPSF = NULL;
    psVector *minMask = NULL;

    // as we loop over orders, we need to refer to the initial selection, but we modify the
    // option values to match the current guess: save the max values here:
    int Nx = options->psfTrendNx;
    int Ny = options->psfTrendNy;
    for (int i = orderMin; i <= orderMax; i++) {

        if (Nx > Ny) {
            options->psfTrendNx = i;
            options->psfTrendNy = PS_MAX (orderMin, (int)(i * (Ny / Nx) + 0.5));
        } else {
            options->psfTrendNy = i;
            options->psfTrendNx = PS_MAX (orderMin, (int)(i * (Nx / Ny) + 0.5));
        }

        // free existing data, if any
        psFree(psfTry->psf);
        psFree(psfTry->mask);

        // allocate a mask and a psf model using the current Nx,Ny order values;
        psfTry->psf = pmPSFAlloc (options);
        psfTry->mask = psVectorCopy (NULL, srcMask, PS_TYPE_VECTOR_MASK);

        // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
        if (!pmPSFtryMakePSF (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, measure pmPSFtry->metric
        if (!pmPSFtryFitPSF (psfTry, options, maskVal, markVal)) {
            psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
            psFree(psfTry);
            return NULL;
        }

        // stage 4: measure systematic error in pmPSFtry->metric
        if (!pmPSFtryMetric (psfTry)) {
            psError(PS_ERR_UNKNOWN, false, "failed to measure systematic error of metric");
            psFree(psfTry);
            return NULL;
        }

        if (!minPSF) {
            minPSF = psMemIncrRefCounter(psfTry->psf);
            minMask = psMemIncrRefCounter(psfTry->mask);
        }

        if (psfTry->psf->dApResid < minPSF->dApResid) {
            psFree (minPSF);
            psFree (minMask);
            minPSF = psMemIncrRefCounter(psfTry->psf);
            minMask = psMemIncrRefCounter(psfTry->mask);
        }
    }
    psFree (srcMask);

    // keep the ones matching the min systematic error:
    psFree (psfTry->psf);
    psFree (psfTry->mask);
    psfTry->psf = minPSF;
    psfTry->mask = minMask;

    // XXXXX this is probably not used any more.  Are the chisq of the fits so bad? can we
    // fix them by softening the errors on the brightest pixels?

    // 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
    if (options->chiFluxTrend) {
        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));
    }

    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);
}

