# include "psphotInternal.h"

// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful

static int NfitPSF = 0;
static int NfitBlend = 0;
static int NfitDBL = 0;
static int NfitEXT = 0;
static int NfitPCM = 0;

bool psphotFitInit (int nThreads) {
    psTimerStart ("psphot.fits");
    pmSourceFitSetInit (nThreads);
    return true;
}

bool psphotFitSummary () {

    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
	      NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
    return true;
}

bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {

    float x, y, dR;

    pmSourceFitOptions options = *fitOptions;

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

    // if this source is not a possible blend, just fit as PSF
    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
        bool status = psphotFitPSF (readout, source, psf, fitOptions, maskVal, markVal);
        return status;
    }

    // save the PSF model from the Ensemble fit
    pmModel *PSF = pmModelCopy (source->modelPSF);

    if (isnan(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in blend fit primary");

    x = PSF->params->data.F32[PM_PAR_XPOS];
    y = PSF->params->data.F32[PM_PAR_YPOS];

    psArray *modelSet = psArrayAllocEmpty (source->blends->n + 1);
    psArrayAdd (modelSet, 16, PSF);

    psArray *sourceSet = psArrayAllocEmpty (source->blends->n + 1);
    psArrayAdd (sourceSet, 16, source);

    psTrace ("psphot", 4, "fitting blended source at %f %f : %f\n", source->peak->xf, source->peak->yf, source->peak->rawFlux);

    // we need to include all blends in the fit (unless primary is saturated?)
    dR = 0;
    for (int i = 0; i < source->blends->n; i++) {
        pmSource *blend = source->blends->data[i];

        // find the blend which is furthest from source
        dR = PS_MAX (dR, hypot (blend->peak->xf - x, blend->peak->yf - y));

        // create the model and guess parameters for this blend
        pmModel *model = pmModelAlloc (PSF->type);
        for (int j = 0; j < model->params->n; j++) {
            model->params->data.F32[j] = PSF->params->data.F32[j];
            model->dparams->data.F32[j] = PSF->dparams->data.F32[j];
        }

        // XXX assume local sky is 0.0?
        model->params->data.F32[PM_PAR_I0] = blend->peak->rawFlux;
        model->params->data.F32[PM_PAR_XPOS] = blend->peak->xf;
        model->params->data.F32[PM_PAR_YPOS] = blend->peak->yf;

        // these should never be invalid values
        // XXX drop these tests eventually
        if (isnan(model->params->data.F32[PM_PAR_I0]))   psAbort("nan in blend fit");
        if (isnan(model->params->data.F32[PM_PAR_XPOS])) psAbort("nan in blend fit");
        if (isnan(model->params->data.F32[PM_PAR_YPOS])) psAbort("nan in blend fit");

        // add this blend to the list
        psArrayAdd (modelSet, 16, model);
        psArrayAdd (sourceSet, 16, blend);

	psTrace ("psphot", 5, "adding source at %f %f : %f\n", blend->peak->xf, blend->peak->yf, blend->peak->rawFlux);

        // free to avoid double counting model
        psFree (model);
    }

    // extend source radius as needed
    psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, dR);

    // fit PSF model
    options.mode = PM_SOURCE_FIT_PSF;
    pmSourceFitSet (source, modelSet, &options, maskVal);

    // clear the circular mask
    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 

    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");

    // evaluate the blend objects, subtract if good, free otherwise
    for (int i = 1; i < modelSet->n; i++) {
        pmSource *blend = sourceSet->data[i];
        pmModel *model  = modelSet->data[i];

	if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");

        // if this one failed, skip it
        if (!psphotEvalPSF (blend, model)) {
            psTrace ("psphot", 4, "failed on blend %d of %ld\n", i, modelSet->n);
            continue;
        }

        // otherwise, supply the resulting model to the corresponding blend
        psFree(blend->modelPSF);
        blend->modelPSF = psMemIncrRefCounter (model);
        psTrace ("psphot", 5, "fitted blend as PSF\n");

        // build cached model and subtract
        pmSourceCacheModel (blend, maskVal);
        pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal);
        blend->mode |=  PM_SOURCE_MODE_BLEND_FIT;
    }
    NfitBlend += modelSet->n;

    // evaluate the primary object
    if (!psphotEvalPSF (source, PSF)) {
        psTrace ("psphot", 4, "failed on blend 0 of %ld\n", modelSet->n);
        psFree (PSF);
        psFree (modelSet);
        psFree (sourceSet);
        return false;
    }
    psFree (modelSet);
    psFree (sourceSet);

    // save the new, successful model
    psFree (source->modelPSF);
    source->modelPSF = PSF;
    psTrace ("psphot", 5, "fitted primary as PSF\n");

    // build cached model and subtract
    pmSourceCacheModel (source, maskVal);
    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    source->mode |=  PM_SOURCE_MODE_BLEND_FIT;
    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    return true;
}

bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {

    pmSourceFitOptions options = *fitOptions;

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

    NfitPSF ++;

    // save the PSF model from the Ensemble fit
    pmModel *PSF = pmModelCopy (source->modelPSF);
    if (isnan(PSF->params->data.F32[1])) psAbort("nan in psf fit");

    // extend source radius as needed
    psphotCheckRadiusPSF (readout, source, PSF, markVal);

    // fit PSF model (set/unset the pixel mask)
    options.mode = PM_SOURCE_FIT_PSF;
    pmSourceFitModel (source, PSF, &options, maskVal);

    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");

    // clear the circular mask
    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 

    // does the PSF model succeed?
    if (!psphotEvalPSF (source, PSF)) {
        psFree (PSF);
        return false;
    }

    // free old model, save new model
    psFree (source->modelPSF);
    source->modelPSF = PSF;
    psTrace ("psphot", 5, "fitted as PSF\n");

    // build cached model and subtract
    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    pmSourceCacheModel (source, maskVal);
    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    return true;
}

// save a local, static copy of the EXT model type so we don't have to lookup for each object
static pmModelType modelTypeEXT;

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

    bool status;

    // extended source model descriptions
    char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL");
    modelTypeEXT = pmModelClassGetType (modelNameEXT);

    psphotInitRadiusEXT (recipe, readout);

    return true;
}

bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {

    float fitRadius, windowRadius;
    bool okEXT, okDBL;
    pmModel *ONE = NULL;
    pmSource *tmpSrc = NULL;
    pmModel *EXT = NULL;
    psArray *DBL = NULL;
    // pmMoments psfMoments;

    // skip the source if we don't think it is extended
    // XXX are these robust, or do we have better info from the source-size analysis??
    if (source->type == PM_SOURCE_TYPE_UNKNOWN) return false;
    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;

# define TEST_X -420.0
# define TEST_Y 300.0
   
    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
	fprintf (stderr, "test galaxy\n");
    }

# undef TEST_X
# undef TEST_Y

    // set the radius based on the footprint (also sets the mask pixels)
    if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false;
    // fprintf (stderr, "rad: %6.1f %6.1f  | %5.2f %5.2f %5.2f  ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius);

    psTrace ("psphot", 5, "trying blob...\n");

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

    float chiEXT = NAN;
    float chiDBL = NAN;

    // this temporary source is used as a place-holder by the psphotEval functions below
    tmpSrc = pmSourceAlloc ();
    {
	// DBL will always be defined, but DBL->data[n] might not
	DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal);
	if (!DBL) goto escape;
	if (!DBL->n) goto escape;

	okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
	okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
	okDBL = false; // XXX this is failing badly...
	// XXX should I keep / save the flags set in the eval functions?

	// correct first model chisqs for flux trend
	ONE = DBL->data[0];
	if (ONE) {
	    psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
	    chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
	    ONE->fitRadius = fitRadius;
	}

	// correct second model chisqs for flux trend
	ONE = DBL->data[1];
	if (ONE) {
	    psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
	    ONE->fitRadius = fitRadius;
	}
    }

    { 
	// XXX need to handle failures better here
	EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
	if (!EXT) goto escape;
	if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;

	okEXT = psphotEvalEXT (tmpSrc, EXT);
	chiEXT = EXT ? EXT->chisqNorm : NAN;
	EXT->fitRadius = fitRadius;
    }

    // clear the circular mask
    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 

    psFree (tmpSrc);

    // (void) psTraceSetLevel("psModules.objects.pmSourceFitSet", 0);

    if (okEXT && okDBL) {
        // XXX EAM : a bogus bias: need to examine this better
        if (3*chiEXT > chiDBL) goto keepDBL;
        goto keepEXT;
    }

    if (okEXT && !okDBL) goto keepEXT;
    if (!okEXT && okDBL) goto keepDBL;

    psTrace ("psphot", 4, "both failed: blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);

    // both models failed; reject them both
    // XXX -- change type flags to psf in this case, and make sure we subtract it?
    // reset the psf moments
    // XXX *source->moments = psfMoments;

    psFree (EXT);
    psFree (DBL);
    return false;

keepEXT:
    psTrace ("psphot", 4, "goto EXT : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
    // sub EXT
    psFree (DBL);

    // save new model
    // XXX save the correct radius...
    source->modelEXT = EXT;
    source->type = PM_SOURCE_TYPE_EXTENDED;
    source->mode |= PM_SOURCE_MODE_EXTMODEL;
    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;

    // build cached model and subtract
    pmSourceCacheModel (source, maskVal);
    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);

# if (PS_TRACE_ON)   
    psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
    if (psTraceGetLevel("psphot") >= 6) {
	psLogMsg ("psphot", 1, "source 2:\n");
	for (int i = 0; i < source->modelEXT->params->n; i++) {
	    psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelEXT->params->data.F32[i], source->modelEXT->dparams->data.F32[i]);
	}
    }
# endif

    // reset the psf moments
    // XXX *source->moments = psfMoments;
    return true;

keepDBL:
    psTrace ("psphot", 4, "goto DBL : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
    // sub DLB
    psFree (EXT);

    // drop old model, save new second model...
    psFree (source->modelPSF);
    source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
    source->mode     |= PM_SOURCE_MODE_PAIR;
    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;

    // copy most data from the primary source (modelEXT, blends stay NULL)
    pmSource *newSrc = pmSourceCopy (source);
    newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
    newSrc->peak->footprint = source->peak->footprint; // just a reference; the peak does not own the footprint
    psArrayAdd(newSrc->peak->footprint->peaks, 1, newSrc->peak); // the footprint owns the peak

    // build cached models and subtract
    pmSourceCacheModel (source, maskVal);
    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    pmSourceCacheModel (newSrc, maskVal);
    pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal);

# if (PS_TRACE_ON)   
    psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
    if (psTraceGetLevel("psphot") >= 6) {
	psLogMsg ("psphot", 1, "source 1:\n");
	for (int i = 0; i < newSrc->modelPSF->params->n; i++) {
	    psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, newSrc->modelPSF->params->data.F32[i], newSrc->modelPSF->dparams->data.F32[i]);
	}
	psLogMsg ("psphot", 1, "source 2:\n");
	for (int i = 0; i < source->modelPSF->params->n; i++) {
	    psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
	}
    }
# endif

    psArrayAdd (newSources, 100, newSrc);
    psFree (newSrc);
    psFree (DBL);
    return true;

escape:
    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
    psFree (tmpSrc);
    psFree (EXT);
    psFree (DBL);
    return false;
}

// fit a double PSF source to an extended blob
psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {

    float dx, dy;
    pmModel *DBL;
    pmModel *PSF;
    psEllipseAxes axes;
    psEllipseMoments moments;
    psArray *modelSet;

    pmSourceFitOptions options = *fitOptions;

    NfitDBL ++;

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

    // XXX this is really poor: if we don't have moments for the source, we have no guess.
    // force the measurement?
    psAssert (source->moments, "moments are re-calculated for any extended source");
	
    // make a guess at the position of the two sources
    moments.x2 = source->moments->Mxx;
    moments.xy = source->moments->Mxy;
    moments.y2 = source->moments->Myy;
    axes = psEllipseMomentsToAxes (moments, 20.0);

    if (isnan(axes.major)) return NULL;
    if (isnan(axes.minor)) return NULL;
    if (isnan(axes.theta)) return NULL;

    // XXX this is really arbitrary: 4 pixel separation?
    dx = 2 * cos (axes.theta);
    dy = 2 * sin (axes.theta);

    // save the PSF model from the Ensemble fit
    PSF = source->modelPSF;
    if (isnan(PSF->params->data.F32[1])) psAbort("nan in dbl fit");

    modelSet = psArrayAlloc (2);

    DBL = pmModelCopy (PSF);
    DBL->params->data.F32[PM_PAR_I0]  *= 0.5;
    DBL->params->data.F32[PM_PAR_XPOS] = source->peak->xf + dx;
    DBL->params->data.F32[PM_PAR_YPOS] = source->peak->yf + dy;
    modelSet->data[0] = DBL;

    DBL = pmModelCopy (PSF);
    DBL->params->data.F32[PM_PAR_I0]  *= 0.5;
    DBL->params->data.F32[PM_PAR_XPOS] = source->peak->xf - dx;
    DBL->params->data.F32[PM_PAR_YPOS] = source->peak->yf - dy;
    modelSet->data[1] = DBL;

    // fit PSF model (set/unset the pixel mask)
    options.mode = PM_SOURCE_FIT_PSF;
    pmSourceFitSet (source, modelSet, &options, maskVal);
    return (modelSet);
}

pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {

    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    }

    pmSourceFitOptions options = *fitOptions;

    NfitEXT ++;

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

    // use the source moments, etc to guess basic model parameters
    pmModel *model = pmSourceModelGuess (source, modelType);
    if (!model) {
	psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
	return NULL;
    }

    // for sersic models, use a grid search to choose an index, then float the params there
    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    	psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
    }

    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
	options.mode = PM_SOURCE_FIT_NO_INDEX;
    } else {
	options.mode = PM_SOURCE_FIT_EXT;
    }

    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    pmSourceFitModel (source, model, &options, maskVal);
    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);

    return (model);
}

# define TIMING 0

pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {

    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    }

    pmSourceFitOptions options = *fitOptions;

    NfitPCM ++;

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

    // allocate the model
    pmModel *model = pmModelAlloc(modelType);
    if (!model) {
	return NULL;
    }

# define TEST_X -540.0
# define TEST_Y 540.0
   
    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    }

    float t1, t2, t4, t5;
    if (TIMING) { psTimerStart ("psphotFitPCM"); }

    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
    if (!pcm) {
	psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
        model->flags |= PM_MODEL_STATUS_BADARGS; // XXX this is probably already set in pmPCMinit
	return model;
    }
    if (TIMING) { t1 = psTimerMark ("psphotFitPCM"); }

    // NOTE : 65 allocs to here
    // get the guess for sersic models 
    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
	// use the source moments, etc to guess basic model parameters
	if (!psphotSersicModelClassGuessPCM (pcm, source)) {
	    psFree (pcm);
	    model->flags |= PM_MODEL_STATUS_BADARGS;
	    return model;
	}
    } else {
	// use the source moments, etc to guess basic model parameters
	if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
	    psFree (pcm);
	    model->flags |= PM_MODEL_STATUS_BADARGS;
	    return model;
	}
    }

    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }

    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
	options.mode = PM_SOURCE_FIT_NO_INDEX;
    } else {
	options.mode = PM_SOURCE_FIT_EXT;
    }
    // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
    pmPCMupdate(pcm, source, &options, model);
    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }

    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);  // NOTE : 1687 allocs in here
    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }

    if (TIMING) {
	int nPixBig = source->pixels->numCols * source->pixels->numRows;
    	fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
    }

    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    }

    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    psFree (pcm);

    return model;
}

# undef TEST_X
# undef TEST_Y

// note that these should be 1/2n of the standard sersic index
// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0};
// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
float indexGuess[] = {1.0, 2.0, 3.0, 4.0};
# define N_INDEX_GUESS 4

// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {

    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));

    pmSourceFitOptions options = *fitOptions;
    
    // fit EXT (not PSF) model (set/unset the pixel mask)
    options.mode = PM_SOURCE_FIT_NO_INDEX;
    options.nIter = 4;

    int iMin = -1;
    float xMin = NAN;
    float chiSquare[N_INDEX_GUESS];

# define TEST_X -540.0
# define TEST_Y 540.0
   
    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    }

    for (int i = 0; i < N_INDEX_GUESS; i++) {
	model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];

	if (!model->modelGuess(model, source)) {
	    model->flags |= PM_MODEL_STATUS_BADARGS;
	    return false;
	}

	pmSourceFitModel (source, model, &options, maskVal);
	// fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);

	chiSquare[i] = model->chisqNorm;
	if (i == 0) {
	    xMin = chiSquare[i];
	    iMin = i;
	} else {
	    if (chiSquare[i] < xMin) {
		xMin = chiSquare[i];
		iMin = i;
	    }
	}
    }
    assert (iMin >= 0);

    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    }

    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
    model->modelGuess(model, source);

    return true;
}

// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {

    pmModel *model = pcm->modelConv;

    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));

    pmSourceFitOptions options = *fitOptions;
    
    // fit EXT (not PSF) model (set/unset the pixel mask)
    options.mode = PM_SOURCE_FIT_NO_INDEX;
    options.nIter = 4;

    // update the pcm elements if we have changed the circumstance (here, options.mode)
    pmPCMupdate(pcm, source, &options, model);

    int iMin = -1;
    float xMin = NAN;
    float chiSquare[N_INDEX_GUESS];

    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
    }

    for (int i = 0; i < N_INDEX_GUESS; i++) {
	model->params->data.F32[PM_PAR_7] = indexGuess[i];
	
	if (!model->modelGuess(model, source)) {
	    model->flags |= PM_MODEL_STATUS_BADARGS;
	    return false;
	}

# if (0)
	// this block is to test the relative speed of straight and PCM fits
	pmSourceFitModel (source, model, &options, maskVal);
# else
	pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
	pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
# endif
	fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
	// fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);

	chiSquare[i] = model->chisq;
	if (i == 0) {
	    xMin = chiSquare[i];
	    iMin = i;
	} else {
	    if (chiSquare[i] < xMin) {
		xMin = chiSquare[i];
		iMin = i;
	    }
	}
    }
    assert (iMin >= 0);
    
    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
    }

    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];

    pmSourceModelGuessPCM(pcm, source, maskVal, markVal);

    return true;
}
