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

# include "psphot.h"
bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);

// 2006.02.07 : no leaks!
// fit all reasonable sources with the linear PSF model
bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { 

    bool  status;
    float x;
    float y;
    float f;
    float r;

    // psRegion allArray = psRegionSet (0, 0, 0, 0);

    psTimerStart ("psphot");

    // source analysis is done in spatial order
    sources = psArraySort (sources, psphotSortByY);

    // radius for source moments if star is saturated
    // XXX we already do this in psphotSourceStats, is it still needed?
    float INNER_RADIUS     = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");

    // setup the PSF fit details
    psphotInitRadiusPSF (recipe, psf->type);

    // pre-calculate all model pixels
    psArray  *models = psArrayAlloc (sources->n);
    psVector *index  = psVectorAlloc (sources->n, PS_TYPE_U32);
    // DROP models->n = index->n = 0;

    // option to limit analysis to a specific region
    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    psRegion AnalysisRegion = psRegionFromString (region);
    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
    // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");

    for (int i = 0; i < sources->n; i++) {
	pmSource *inSource = sources->data[i];

	// skip non-astronomical objects (very likely defects)
	// XXX EAM : should we try these anyway?
	if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; 
	if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue;
	if (final) {
	    if (inSource->mode &  PM_SOURCE_MODE_SUBTRACTED) continue;
	} else {
	    if (inSource->mode &  PM_SOURCE_MODE_BLEND) continue;
	}

	if (inSource->moments->x < AnalysisRegion.x0) continue;
	if (inSource->moments->y < AnalysisRegion.y0) continue;
	if (inSource->moments->x > AnalysisRegion.x1) continue;
	if (inSource->moments->y > AnalysisRegion.y1) continue;

	pmSource *otSource = pmSourceAlloc ();

	// really saturated stars should be re-measured for a better centroid
	// XXX EAM : move this to a 'clear satstar function'
	if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
	    status = pmSourceMoments (inSource, INNER_RADIUS);
	}

	// XXX EAM : add option to use EXT or PSF form
	// use the source moments, etc to guess basic model parameters
	pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 
	if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
	    modelEXT->params->data.F32[2] = inSource->moments->x;
	    modelEXT->params->data.F32[3] = inSource->moments->y;
	} else {
	    // peak-up on peak (for non-sat objects)

	    // ix,iy must land on inSource->pixels 
	    int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
	    int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));

	    psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
	    psPlane min = psImageBicubeMin (bicube);

	    psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", 
		     modelEXT->params->data.F32[2], modelEXT->params->data.F32[3], min.x + ix, min.y + iy);
	    
	    // if min point is too deviant, keep the old value
	    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
		modelEXT->params->data.F32[2] = min.x + ix;
		modelEXT->params->data.F32[3] = min.y + iy;
	    }
	    psFree (bicube);
	}

	// set PSF parameters for this model
	pmModel *model = pmModelFromPSF (modelEXT, psf);
	psFree (modelEXT);

	// save the original coords
	x = model->params->data.F32[2];
	y = model->params->data.F32[3];

	// set the fit radius based on the object flux limit and the model
	psphotCheckRadiusPSF (readout, inSource, model);

	// make temporary copies of the image pixels and mask
	otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
	otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
	otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);

	// build the model image 
	psImage *flux = otSource->pixels;
	psImage *mask = otSource->mask;

	// set model to unit peak, zero sky (we assume sky is subtracted)
	model->params->data.F32[0] = 0.0;
	model->params->data.F32[1] = 1.0;

	// fill in the model pixel values
	psImageInit (flux, 0.0);
	psImageKeepCircle (mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);
	pmModelAdd (flux, mask, model, false, false);

	// calculate nDOF (nPix - 1)
	// int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT);
	// model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;

	// save source in list
	otSource->modelPSF = model;
	index->data.U32[models->n] = i;
	psArrayAdd (models, 100, otSource);
	psFree (otSource);
    }
    psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n);
    
    // fill out the sparse matrix
    psSparse *sparse = psSparseAlloc (models->n, 100);
    psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
    weight->n = models->n;

    for (int i = 0; i < models->n; i++) {
	int N = index->data.U32[i];
	pmSource *Fi = sources->data[N];
	pmSource *Mi = models->data[i];

	// scale by diagonal element (auto-cross-product)
	r = pmSourceCrossProduct (Mi, Mi);
	weight->data.F32[i] = r;

	psSparseMatrixElement (sparse, i, i, 1.0);

	// find the image x model value
	f = pmSourceCrossProduct (Fi, Mi);
	psSparseVectorElement (sparse, i, f / r);

	// loop over all other stars following this one
	for (int j = i + 1; j < models->n; j++) {
	    pmSource *Mj = models->data[j];

	    // skip over disjoint source images, break after last possible overlap
	    if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break;
	    if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
	    if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue;
	    if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;
	    
	    // got an overlap; calculate cross-product and add to output array
	    f = pmSourceCrossProduct (Mi, Mj);
	    psSparseMatrixElement (sparse, j, i, f / r); 
	}
    }
    psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);

    psSparseConstraint constraint;
    constraint.paramMin   = 0.0;
    constraint.paramMax   = 1e8;
    constraint.paramDelta = 1e8;

    // solve for normalization terms (need include local sky?)
    psSparseResort (sparse);
    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);

    // adjust models, set sources and subtract
    for (int i = 0; i < models->n; i++) {
	int N = index->data.U32[i];
	pmSource *Fi = sources->data[N];
	pmSource *Mi = models->data[i];

	// if we already have a PSF model, free it.
	psFree (Fi->modelPSF);

	// need to increment counter so we can free models here and sources above
	Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF);

	// assign linearly-fitted normalization
	if (isnan(norm->data.F32[i])) {
	    psAbort ("psphot", "ensemble source is nan");
	}
	Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
	Fi->modelPSF->dparams->data.F32[1] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
	// XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it

	// subtract object
	pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
	Fi->mode |= PM_SOURCE_MODE_SUBTRACTED;
	if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB;
    }

    // measure chisq for each source
    for (int i = 0; final && (i < models->n); i++) {
	int N = index->data.U32[i];
	pmSource *Fi = sources->data[N];
	pmModel *model = Fi->modelPSF;

	x = model->params->data.F32[2];
	y = model->params->data.F32[3];

	psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);
	pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
	psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
    }

    psFree (index);
    psFree (sparse);
    psFree (models);
    psFree (norm);
    psFree (weight);

    psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
    return true;
}

bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight) {

    double dC = 0.0;
    int Npix = 0;
    for (int j = 0; j < image->numRows; j++) {
	for (int i = 0; i < image->numCols; i++) {
	    if (mask->data.U8[j][i]) continue;
	    if (weight->data.F32[j][i] <= 0) continue;
	    dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
	    Npix ++;
	}
    }
    model->nDOF = Npix - 1;
    model->chisq = dC;

    return (true);
}
