# include "psphotInternal.h"
# define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY

// XXX add more test information?
bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe) {

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

    // bit-masks to test for good/bad pixels
    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    assert (maskVal);

    // bit-mask to mark pixels not used in analysis
    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    assert (markVal);

    // maskVal is used to test for rejected pixels, and must include markVal
    maskVal |= markVal;

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

    psTimerStart ("modelTest");

    // find the currently selected readout
    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    PS_ASSERT_PTR_NON_NULL (readout, false);

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

    // find the various fitting parameters (try test values first)
    float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
    if (!status || !isfinite(INNER)) {
        INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    }
    float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
    if (!status || !isfinite(OUTER)) {
        OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    }
    float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
    if (!status || !isfinite(RADIUS)) {
        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    }
    float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
    if (!status || !isfinite(mRADIUS)) {
        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");
    if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");

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

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

    // in fitMode, psf sets the model type
    if (fitMode == PM_SOURCE_FIT_PSF) {
        psf = psphotLoadPSF (config, view, recipe);
        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
        modelType = psf->type;
        source->type = PM_SOURCE_TYPE_STAR;
    }
    if (fitMode == PM_SOURCE_FIT_EXT) {
        // find the model: supplied by user or first in the PSF_MODEL list
        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
        if (!status || !strcasecmp (modelName, "DEFAULT")) {
            // get the list pointers for the PSF_MODEL entries

            psList *list = NULL;
            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
            if (mdi == NULL) psAbort("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("missing PSF_MODEL selection");
                list = psMemIncrRefCounter(mdi->data.list);
            }

            // take the first list element
            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
            modelName = item->data.V;
        }
        modelType = pmModelClassGetType (modelName);
        if (modelType < 0) psAbort("unknown model %s", modelName);
        source->type = PM_SOURCE_TYPE_EXTENDED;
    }
    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
        // we need to load BOTH a psf and an ext model
        psf = psphotLoadPSF (config, view, recipe);
        if (!psf) psAbort("PSF_INPUT_FILE not supplied");

        // find the model: supplied by user or first in the PSF_MODEL list
        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
        if (!status || !strcasecmp (modelName, "DEFAULT")) {
            // get the list pointers for the PSF_MODEL entries

            psList *list = NULL;
            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
            if (mdi == NULL) psAbort("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("missing PSF_MODEL selection");
                list = psMemIncrRefCounter(mdi->data.list);
            }

            // take the first list element
            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
            modelName = item->data.V;
        }
        modelType = pmModelClassGetType (modelName);
        if (modelType < 0) psAbort("unknown model %s", modelName);
        source->type = PM_SOURCE_TYPE_EXTENDED;
    }

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

    // get the source moments
    status = pmSourceMoments (source, mRADIUS, 0.0, 1.0, maskVal);
    if (!status) psAbort("psSourceMoments error");
    source->peak->value = source->moments->Peak;

    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->Mx, source->moments->My);
    fprintf (stderr, "moments: %f, %f - %f\n", source->moments->Mxx, source->moments->Myy, source->moments->Mxy);

    psEllipseMoments moments;
    moments.x2 = source->moments->Mxx;
    moments.y2 = source->moments->Myy;
    moments.xy = source->moments->Mxy;
    psEllipseAxes axes = psEllipseMomentsToAxes (moments, 20.0);

    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);
    source->modelEXT = model;

    // if any parameters are defined by the user, take those values
    int nParams = pmModelClassParameterCount (modelType);
    psF32 *params = model->params->data.F32;
    params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value,
    params[PM_PAR_YPOS] = yObj; // XXX or use the centroid
    for (int i = 0; i < nParams; i++) {
        if (i == PM_PAR_XPOS) continue;
        if (i == PM_PAR_YPOS) continue;

        sprintf (name, "TEST_FIT_PAR%d", i);
        value = psMetadataLookupF32 (&status, recipe, name);
        if (status && isfinite (value)) {
            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);

    // for PSF fitting, set the shape parameters based on the PSF & source position
    if (fitMode == PM_SOURCE_FIT_PSF) {
        source->modelPSF = pmModelFromPSF (model, psf);
        psFree (model);
        model = source->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, 20.0);

    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->maskObj, xObj, yObj, RADIUS, "OR", markVal);
    psphotSaveImage (NULL, source->maskObj, "mask1.fits");

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

    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
        // build the psf for the object
        source->modelPSF = pmModelFromPSF (model, psf);
        source->modelEXT = model;

	// what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
	int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
	assert (status);

        model = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
        params = model->params->data.F32;
    } else {
        status = pmSourceFitModel (source, model, fitMode, maskVal);
    }

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

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

    // subtract object, leave local sky
    // pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);

    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->maskObj, "mask.fits");

    psLogMsg ("psphot", PS_LOG_INFO, "model test : %f sec\n", psTimerMark ("modelTest"));

    exit (0);
}
