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

# include "psphot.h"
static char DEFAULT_OPTION[] = "SKYBIAS";

// measure the aperture residual statistics
bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 

    int Nfail = 0;
    int Nskip = 0;
    int Npsf;
    bool status;
    pmModel *model;
    pmSource *source;

    psTimerStart ("psphot");

    // S/N limit to perform full non-linear fits
    float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET");

    // set limits on the aperture magnitudes 
    pmSourceMagnitudesInit (recipe);

    // measure the aperture loss as a function of radius for PSF
    float REF_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_REF_RADIUS");
    psf->growth = pmGrowthCurveAlloc (3.0, 100.0, REF_RADIUS);
    psphotGrowthCurve (readout, psf);
    
    psVector *mask    = psVectorAlloc (300, PS_TYPE_U8);
    psVector *xPos    = psVectorAlloc (300, PS_TYPE_F64);
    psVector *yPos    = psVectorAlloc (300, PS_TYPE_F64);
    psVector *flux    = psVectorAlloc (300, PS_TYPE_F64);
    psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64);
    psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
    psVector *dMag    = psVectorAlloc (300, PS_TYPE_F64);
    mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0;
    Npsf = 0;

    // select all good PM_SOURCE_TYPE_STAR entries
    for (int i = 0; i < sources->n; i++) {
	source = sources->data[i];
	model = source->modelPSF;

	if (source->type != PM_SOURCE_TYPE_STAR) continue; 
	if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
	if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
	if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
	if (source->mode &  PM_SOURCE_MODE_POOR) continue;

	// get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
	// will fail if below S/N threshold or model is missing
	if (!pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH)) {
	    Nskip ++;
	    continue;
	}

	apResid->data.F64[Npsf] = source->apMag - source->psfMag;

	xPos->data.F64[Npsf]    = model->params->data.F32[2];
	yPos->data.F64[Npsf]    = model->params->data.F32[3];

	flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
	r2rflux->data.F64[Npsf] = PS_SQR(model->radiusTMP) / flux->data.F64[Npsf];

	mask->data.U8[Npsf] = 0;

	// XXX sanity clip?
	// XXX need to see if all data were tossed?
	// XXX need to subtract median?
	// XXX need to put this in the config data...
	// if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;
	if (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET) {
	    Nfail ++;
	    continue;
	}

	dMag->data.F64[Npsf] = model->dparams->data.F32[1] / model->params->data.F32[1];

	psVectorExtend (mask, 	 100, 1);
	psVectorExtend (xPos, 	 100, 1);
	psVectorExtend (yPos, 	 100, 1);
	psVectorExtend (flux,    100, 1);
	psVectorExtend (r2rflux, 100, 1);
	psVectorExtend (dMag,    100, 1);
	psVectorExtend (apResid, 100, 1);
	Npsf ++;
    }
    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n", 
	      psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);

    // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL 
    char *ApTrendOption = psMetadataLookupPtr (&status, recipe, "APTREND");
    if (!status) ApTrendOption = DEFAULT_OPTION;

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

    // no correction
    if (!strcasecmp (ApTrendOption, "NONE")) { 
	// remove ApTrend fit from pmPSFtry
	psf->ApTrend->coeff[0][0][0][0] = 0;
    } else if (!strcasecmp (ApTrendOption, "CONSTANT")) { // constant only
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
	    return false;
	}
    } else if (!strcasecmp (ApTrendOption, "SKYBIAS")) { // constant and skybias only
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_SKYBIAS);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
	    return false;
	}
    } else if (!strcasecmp (ApTrendOption, "SKYSAT")) {
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    return false;
	}

	// apply the fit
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_SKYBIAS);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_SKYSAT);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
	    return false;
	}
    } else if (!strcasecmp (ApTrendOption, "XY_LIN")) {	// constant and linear X,Y only
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_XY_LIN);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
	    return false;
	}
    } else if (!strcasecmp (ApTrendOption, "XY_QUAD")) { // constant and quadratic X,Y only
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_XY_QUAD);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
	    return false;
	}
    } else if (!strcasecmp (ApTrendOption, "SKY_XY_LIN")) { // constant and sky, linear X,Y only
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_SKY_XY_LIN);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
	    return false;
	}
    } else if (!strcasecmp (ApTrendOption, "SKYSAT_XY_LIN")) { // constant and sky, linear X,Y only
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_SKYBIAS);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
	    return false;
	}

	// apply the fit
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_SKYSAT_XY_LIN);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
	    return false;
	}
    } else if (!strcasecmp (ApTrendOption, "ALL")) {
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_CONSTANT);
	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
	    return false;
	}

	// fit just SkyBias and clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_SKYBIAS);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
	    return false;
	}

	// finally, fit x, y, SkyBias and clip out objects which are too far from the median 
	stats->clipIter = 3;
	pmPSF_MaskApTrend (psf->ApTrend, PM_PSF_ALL);
	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
	if (psf->ApTrend == NULL) {
	    psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
	    return false;
	}
    } else {
	psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", ApTrendOption);
	return false;
    }

    // construct the fitted values and the residuals
    psVector *fit   = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux);
    psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit);

    // measure scatter for sources with dMag < 0.01 (S/N = 100)
    int Nkeep = 0;
    psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
    for (int i = 0; i < dMag->n; i++) {
	if (dMag->data.F64[i] > 0.01) {
	    mask->data.U8[i] |= 0x02;
	} 
	if (! mask->data.U8[i]) Nkeep ++;
    }
    residStats  = psVectorStats (residStats, resid, NULL, mask, 0x03);

    // apply ApTrend results
    psf->skyBias  = psf->ApTrend->coeff[0][0][1][0];
    psf->skySat   = psf->ApTrend->coeff[0][0][0][1];
    psf->ApResid  = psf->ApTrend->coeff[0][0][0][0];
    psf->dApResid = residStats->sampleStdev;
    psf->ApTrend->coeff[0][0][1][0] = 0;
    psf->ApTrend->coeff[0][0][0][1] = 0;
    psf->nApResid = Npsf;

    // save results for later output
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->growth->apLoss);
    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);

    psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
    psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n", 
	      psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf);
    psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", 
	      1e3*psf->ApTrend->coeff[1][0][0][0],
	      1e6*psf->ApTrend->coeff[2][0][0][0],
	      1e6*psf->ApTrend->coeff[1][1][0][0],
	      1e3*psf->ApTrend->coeff[0][1][0][0],
	      1e6*psf->ApTrend->coeff[0][2][0][0]);

    psFree (mask);
    psFree (xPos);
    psFree (yPos);
    psFree (flux);
    psFree (r2rflux);
    psFree (apResid);
    psFree (dMag);

    psFree (fit);
    psFree (resid);
    psFree (stats);
    psFree (residStats);
    return true;
}

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

