# include "psphotInternal.h"
// XXXX this code fails if there are too few sources to measure the aperture residual
// the larger problem is that the rules for accepting more polynomial terms are weak.

static pmPSFApTrendOptions DEFAULT_OPTION = PM_PSF_APTREND_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;

    PS_ASSERT_PTR_NON_NULL(psf, false);
    PS_ASSERT_PTR_NON_NULL(psf->ApTrend, false);
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_PTR_NON_NULL(sources, false);
    PS_ASSERT_PTR_NON_NULL(recipe, false);

    psTimerStart ("psphot");

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

    // this is the smallest radius allowed: need to at least extend growth curve down to this...
    float PSF_FIT_PAD   = psMetadataLookupF32 (&status, recipe, "PSF_FIT_PADDING");
    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
    int NSTAR_APERTURE_CORRECTION_MIN =
	psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN");
    if (!status) {
	NSTAR_APERTURE_CORRECTION_MIN = 5;
    }

    pmSourcePhotometryMode photMode = 0;
    if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH;
    if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP;

    // 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 (PSF_FIT_PAD, 100.0, REF_RADIUS);

    if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH)) {
	psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
	psFree(psf->growth); psf->growth = NULL;
	return false;
    }

    psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_U8);
    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
    psVector *flux    = psVectorAllocEmpty (300, PS_TYPE_F64);
    psVector *r2rflux = psVectorAllocEmpty (300, PS_TYPE_F64);
    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64);
    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F64);
    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, photMode)) {
            Nskip ++;
            continue;
        }

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

        xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
        yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];

        flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
        r2rflux->data.F64[Npsf] = PS_SQR(model->radiusFit) / 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 ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {
            Nfail ++;
            continue;
        }

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

        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", PS_LOG_DETAIL, "measure aperture residuals for %d objects (%d skipped, %d failed, %ld invalid)\n",
              Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);

    // XXX choose a better value here?
    if (Npsf < NSTAR_APERTURE_CORRECTION_MIN) {
        psError(PSPHOT_ERR_APERTURE, true, "Only %d valid aperture residual sources (need %d), giving up",
		Npsf, NSTAR_APERTURE_CORRECTION_MIN);
        return false;
    }

    // APTREND options : NONE CONSTANT SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL
    // APTREND options are used in the switch block below
    pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION;
    char *optionName = psMetadataLookupPtr (&status, recipe, "APTREND");
    if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName);
    if (ApTrendOption == PM_PSF_APTREND_ERROR) {
        psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
        return false;
    }

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

#define P_APTREND_SWITCH_CLEANUP	/* Cleanup memory on error in ApTrendOption switch */ \
    psFree(psf->growth); psf->growth = NULL; \
    psFree(mask); \
    psFree(xPos); \
    psFree(yPos); \
    psFree(flux); \
    psFree(r2rflux); \
    psFree(apResid); \
    psFree(dMag); \
    psFree(stats)

    // no correction
    switch (ApTrendOption) {
      case PM_PSF_APTREND_NONE:
        // remove ApTrend fit from pmPSFtry
        psf->ApTrend->coeff[0][0][0][0] = 0;
        break;
      case PM_PSF_APTREND_CONSTANT:
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      case PM_PSF_APTREND_SKYBIAS:
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      case PM_PSF_APTREND_SKYSAT:
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      case PM_PSF_APTREND_XY_LIN:
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      case PM_PSF_APTREND_XY_QUAD:
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      case PM_PSF_APTREND_SKY_XY_LIN:
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      case PM_PSF_APTREND_SKYSAT_XY_LIN:
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// apply the fit
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      case PM_PSF_APTREND_ALL:
	// first clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
	if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// fit just SkyBias and clip out objects which are too far from the median 
	stats->clipIter = 2;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	// finally, fit x, y, SkyBias and clip out objects which are too far from the median 
	stats->clipIter = 3;
	pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
	if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
	    psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
	    P_APTREND_SWITCH_CLEANUP;
	    return false;
	}
	break;
      default:
        psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
        return false;
    }
#undef P_APTREND_SWITCH_CLEANUP
    
    // 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 ++;
    }
    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", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n", 
	      Nkeep, Npsf, psTimerMark ("psphot"));
    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f : %f bias, %f skysat\n",
              psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat);
    psLogMsg ("psphot.apresid", PS_LOG_MINUTIA, "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)
*/

