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

// XXX consider this function : add more test information?
bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {

    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, 1.0, 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("PSF_INPUT_FILE not supplied");
        psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, psfFile, FALSE);
        psf = pmPSFfromMetadata (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("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
            item = psListGet (list, PS_LIST_HEAD);
            modelName = item->data.V;
        }
        modelType = pmModelSetType (modelName);
        if (modelType < 0) psAbort("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, maskVal, mark);
    if (!status) psAbort("pmSourceLocalSky error");

    // get the source moments
    status = pmSourceMoments (source, mRADIUS);
    if (!status) psAbort("pmSourceLocalSky error");
    source->peak->value = 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, 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);
    // 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, 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", mark);

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

    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\n", obsMag, fitMag, obsMag - fitMag);

    // 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);

    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");

    exit (0);
}
