# include "psphotInternal.h"

bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {

    bool status = false;
    int Nap = 0;

    psTimerStart ("psphot.mags");

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    assert (recipe);

    // determine the number of allowed threads
    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    if (!status) {
	nThreads = 0;
    }

    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
    pmReadout *backModel = psphotSelectBackground (config, view);
    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);

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

    pmSourceMagnitudesInit (recipe);

    // the binning details are saved on the analysis metadata
    psImageBinning *binning = psMetadataLookupPtr(&status, recipe, "PSPHOT.BACKGROUND.BINNING");

    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");

    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_APCORR | PM_SOURCE_PHOT_WEIGHT;
    if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH;
    if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP;

    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    int Cx = 1, Cy = 1;
    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);

    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);

    for (int i = 0; i < cellGroups->n; i++) {

	psArray *cells = cellGroups->data[i];

	for (int j = 0; j < cells->n; j++) {

	    // allocate a job -- if threads are not defined, this just runs the job
	    psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES");

	    psArrayAdd(job->args, 1, cells->data[j]); // sources
	    psArrayAdd(job->args, 1, psf);
	    psArrayAdd(job->args, 1, binning);
	    psArrayAdd(job->args, 1, backModel);
	    psArrayAdd(job->args, 1, backStdev);

	    PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
	    PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
	    PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp

	    if (!psThreadJobAddPending(job)) {
		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
		psFree (job);
		return false;
	    }
	    psFree(job);

# if (0)
		int nap = 0;
		if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) {
		    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
		    return false;
		}
		Nap += nap;
# endif
	}

	// wait for the threads to finish and manage results
	if (!psThreadPoolWait (false)) {
	    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
	    return false;
	}

	// we have only supplied one type of job, so we can assume the types here
	psThreadJob *job = NULL;
	while ((job = psThreadJobGetDone()) != NULL) {
	    if (job->args->n < 1) {
		fprintf (stderr, "error with job\n");
	    } else {
		psScalar *scalar = job->args->data[7];
		Nap += scalar->data.S32;
	    }
	    psFree(job);
	}
    }

    psFree (cellGroups);

    psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot.mags"), sources->n, Nap);
    return true;
}

bool psphotMagnitudes_Threaded (psThreadJob *job) {

    bool status;
    int Nap = 0;

    psArray *sources                = job->args->data[0];
    pmPSF *psf                      = job->args->data[1];
    psImageBinning *binning         = job->args->data[2];
    pmReadout *backModel            = job->args->data[3];
    pmReadout *backStdev            = job->args->data[4];
    pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[5],S32);
    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);

    for (int i = 0; i < sources->n; i++) {
        pmSource *source = (pmSource *) sources->data[i];
        status = pmSourceMagnitudes (source, psf, photMode, maskVal);
        if (status && isfinite(source->apMag)) Nap ++;

	if (backModel) {
	    psAssert (binning, "if backModel is defined, so should binning be");
	    source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
	    if (isnan(source->sky) && false) {
		psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
	    }
	} else {
	    source->sky = NAN;
	}

	if (backStdev) {
	    psAssert (binning, "if backStdev is defined, so should binning be");
	    source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
	    if (isnan(source->skyErr) && false) {
		psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
	    }
	} else {
	    source->skyErr = NAN;
	}
    }

    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    psScalar *scalar = job->args->data[7];
    scalar->data.S32 = Nap;

    return true;
}

# if (0)
bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {

    bool status;
    int Nap = 0;

    for (int i = 0; i < sources->n; i++) {
        pmSource *source = (pmSource *) sources->data[i];
        status = pmSourceMagnitudes (source, psf, photMode, maskVal);
        if (status && isfinite(source->apMag)) Nap ++;

	if (backModel) {
	    psAssert (binning, "if backModel is defined, so should binning be");
	    source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
	    if (isnan(source->sky) && false) {
		psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
	    }
	} else {
	    source->sky = NAN;
	}

	if (backStdev) {
	    psAssert (binning, "if backStdev is defined, so should binning be");
	    source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
	    if (isnan(source->skyErr) && false) {
		psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
	    }
	} else {
	    source->skyErr = NAN;
	}
    }

    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    *nap = Nap;

    return true;
}
# endif

bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources) {

    bool status = false;

    psTimerStart ("psphot.mags");

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    assert (recipe);

    // determine the number of allowed threads
    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    if (!status) {
	nThreads = 0;
    }
    nThreads = 0; // XXX until testing is complete, do not thread this function

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

    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    int Cx = 1, Cy = 1;
    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);

    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);

    for (int i = 0; i < cellGroups->n; i++) {

	psArray *cells = cellGroups->data[i];

	for (int j = 0; j < cells->n; j++) {

	    // allocate a job -- if threads are not defined, this just runs the job
	    psThreadJob *job = psThreadJobAlloc ("PSPHOT_PSF_WEIGHTS");

	    psArrayAdd(job->args, 1, cells->data[j]); // sources
	    PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);

	    if (!psThreadJobAddPending(job)) {
		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
		psFree (job);
		return false;
	    }
	    psFree(job);

	}

	// wait for the threads to finish and manage results
	if (!psThreadPoolWait (false)) {
	    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
	    return false;
	}

	// we have only supplied one type of job, so we can assume the types here
	psThreadJob *job = NULL;
	while ((job = psThreadJobGetDone()) != NULL) {
	    if (job->args->n < 1) {
		fprintf (stderr, "error with job\n");
	    } 
	    psFree(job);
	}
    }

    psFree (cellGroups);

    psLogMsg ("psphot", PS_LOG_DETAIL, "measure psf weights : %f sec for %ld objects\n", psTimerMark ("psphot.mags"), sources->n);
    return true;
}

bool psphotPSFWeights_Threaded (psThreadJob *job) {

    bool status;
    bool isPSF;

    psArray *sources                = job->args->data[0];
    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[1],PS_TYPE_IMAGE_MASK_DATA);

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

	// we must have a valid model
	pmModel *model = pmSourceGetModel (&isPSF, source);
	if (model == NULL) {
	  psTrace ("psphot", 3, "fail mag : no valid model");
	  source->pixWeight = NAN;
	  continue;
	}

        status = pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
	if (!status) {
	  psTrace ("psphot", 3, "fail to measure pixel weight");
	  source->pixWeight = NAN;
	  continue;
	}

    }

    return true;
}
