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

# include "psphot.h"
static char DEFAULT_MODE[] = "EXT";

bool psphotModelTest (pmReadout *readout, psMetadata *recipe) {

    bool status;
    int modelType;
    unsigned int Nfail;
    float obsMag, fitMag, value;
    char name[64];
    pmPSF *psf = NULL;
    pmSourceFitMode fitMode;

    psMetadataItem *item  = NULL;

    // use poissonian errors or local-sky errors
    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
    if (!status) POISSON_ERRORS = true;
    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);

    // run model fitting tests on a single source
    if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;

    // what fitting mode to use?
    char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
    if (!status) {
	fitModeWord = DEFAULT_MODE;
    }
    fitMode = PM_SOURCE_FIT_EXT;
    if (!strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;

    // in fitMode, psf sets the model type
    if (fitMode == PM_SOURCE_FIT_PSF) {
	char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE");
	if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
	psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, psfFile, FALSE);
	psf = pmPSFfromMD (psfData);
	modelType = psf->type;
    } else {
	// find the model: supplied by user or first in the PSF_MODEL list
	char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
	if (modelName == NULL) {
	    // get the list pointers for the PSF_MODEL entries

	    psList *list = NULL;
	    psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
	    if (mdi == NULL) psAbort ("psphotModelTest", "missing PSF_MODEL selection");
	    if (mdi->type == PS_DATA_STRING) {
		list = psListAlloc(NULL);
		psListAdd (list, PS_LIST_HEAD, mdi);
	    } else {
		if (mdi->type != PS_DATA_METADATA_MULTI) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
		list = psMemIncrRefCounter(mdi->data.list);
	    }

	    // take the first list element
	    item = psListGet (list, PS_LIST_HEAD);
	    modelName = item->data.V;
	}
	modelType = pmModelSetType (modelName);
	if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName);
    }

    // find the fitting parameters (try test values first)
    float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
    if (!status) {
	INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    }

    float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
    if (!status) {
	OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    }

    float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
    if (!status) {
	RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    }

    float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
    if (!status) {
	mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    }

    // define the source of interest
    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
    float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");

    // construct the source structures
    pmSource *source = pmSourceAlloc();
    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
    pmSourceDefinePixels (source, readout, xObj, yObj, OUTER);

    // find the local sky
    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");

    // get the source moments
    status = pmSourceMoments (source, mRADIUS);
    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
    source->peak->counts = source->moments->Peak;

    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->x, source->moments->y);

    psEllipseMoments moments;
    moments.x2 = source->moments->Sx;
    moments.y2 = source->moments->Sy;
    moments.xy = source->moments->Sxy;
    psEllipseAxes axes = psEllipseMomentsToAxes (moments);

    fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);

    // get the initial model parameter guess
    pmModel *model = pmSourceModelGuess (source, modelType);
    // if any parameters are defined, use those values
    int nParams = pmModelParameterCount (modelType);
    psF32 *params = model->params->data.F32;
    for (int i = 0; i < nParams; i++) {
	if (i == 2) {
	    params[i] = xObj;
	    continue;
	}
	if (i == 3) {
	    params[i] = yObj;
	    continue;
	}
	sprintf (name, "TEST_FIT_PAR%d", i);
	value = psMetadataLookupF32 (&status, recipe, name);
	if (status) {
	    params[i] = value;
	}
    }

    float area = params[4]*params[5];
    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);

    if (fitMode == PM_SOURCE_FIT_PSF) {
	pmModel *modelPSF = pmModelFromPSF (model, psf);
	psFree (model);
	model = modelPSF;
	params = model->params->data.F32;
    }

    // list model input shape
    psEllipseShape shape;
    shape.sx  = 1.4 / model->params->data.F32[4];
    shape.sy  = 1.4 / model->params->data.F32[5];
    shape.sxy = model->params->data.F32[6];
    axes = psEllipseShapeToAxes (shape);

    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);


    fprintf (stderr, "input parameters: \n");
    for (int i = 0; i < nParams; i++) {
	fprintf (stderr, "%d : %f\n", i, params[i]);
    }

    // define the pixels used for the fit
    psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PM_MASK_MARK);

    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    if (status) {
	status = psphotFitSet (source, model, fitset, fitMode);
	exit (0);
    }

    status = pmSourceFitModel (source, model, fitMode);

    // measure the source mags
    pmSourcePhotometryModel (&fitMag, model);
    pmSourcePhotometryAper  (&obsMag, model, source->pixels, source->mask);
    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);

    // write out positive object
    psphotSaveImage (NULL, source->pixels, "object.fits");

    // subtract object, leave local sky
    pmModelSub (source->pixels, source->mask, model, false, false);
    
    fprintf (stderr, "output parameters: \n");
    for (int i = 0; i < nParams; i++) {
	fprintf (stderr, "%d : %f\n", i, params[i]);
    }

    // write out 
    psphotSaveImage (NULL, source->pixels, "resid.fits");
    psphotSaveImage (NULL, source->mask, "mask.fits");

    exit (0);
}
