# include "psphotInternal.h"

psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {

    bool status = false;
    psArray *sources = NULL;

    psTimerStart ("psphot.stats");

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

    // determine properties (sky, moments) of initial sources
    float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    if (!status) return NULL;
    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    if (!status) return NULL;

    psArray *peaks = detections->peaks;
    if (!peaks) {
        psError(PS_ERR_UNEXPECTED_NULL, false, "No peaks found!");
        return NULL;
    }

    // generate the array of sources, define the associated pixel
    sources = psArrayAllocEmpty (peaks->n);

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

        pmPeak *peak = peaks->data[i];
        if (peak->assigned) continue;

        // create a new source
        pmSource *source = pmSourceAlloc();

	// add the peak
        source->peak = psMemIncrRefCounter(peak);

	// allocate space for moments
        source->moments = pmMomentsAlloc();

        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);

	peak->assigned = true;
	psArrayAdd (sources, 100, source);
	psFree (source);
    }

    if (!strcasecmp (breakPt, "PEAKS")) {
	psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
	psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
	psphotVisualShowMoments (sources);
	return sources;
    }

    // threaded measurement of the source magnitudes
    int Nfail = 0;
    int Nmoments = 0;

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

	    psArrayAdd(job->args, 1, cells->data[j]); // sources
	    psArrayAdd(job->args, 1, recipe);
	    PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
	    PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail

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

# if (0)
		int nfail = 0;
		int nmoments = 0;
		if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) {
		    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
		    return NULL;
		}
		Nfail += nfail;
		Nmoments += nmoments;
# endif
	}

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

	// 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 = NULL;
		scalar = job->args->data[2];
		Nmoments += scalar->data.S32;
		scalar = job->args->data[3];
		Nfail += scalar->data.S32;
	    }
	    psFree(job);
	}
    }

    psFree (cellGroups);

    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot.stats"));

    psphotVisualShowMoments (sources);

    return (sources);
}

bool psphotSourceStats_Threaded (psThreadJob *job) {

    bool status = false;
    float BIG_RADIUS;
    psScalar *scalar = NULL;

    psArray *sources                = job->args->data[0];
    psMetadata *recipe              = job->args->data[1];

    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    if (!status) return false;
    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    if (!status) return false;
    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    if (!status) return false;

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

    // threaded measurement of the sources moments
    int Nfail = 0;
    int Nmoments = 0;
    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i];

        // skip faint sources for moments measurement
        if (source->peak->SN < MIN_SN) {
	    source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN;
            continue;
        }

        // measure a local sky value
        // the local sky is now ignored; kept here for reference only
        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
        if (!status) {
	    source->mode |= PM_SOURCE_MODE_SKY_FAILURE;
	    psErrorClear(); // XXX re-consider the errors raised here
	    Nfail ++;
	    continue;
        }

        // measure the local sky variance (needed if noise is not sqrt(signal))
        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
        if (!status) {
	    source->mode |= PM_SOURCE_MODE_SKYVAR_FAILURE;
	    Nfail ++;
	    psErrorClear();
	    continue;
        }

        // measure basic source moments
        status = pmSourceMoments (source, RADIUS);
        if (status) {
            Nmoments ++;
            continue;
        }

        // if no valid pixels, or massive swing, likely saturated source,
        // try a much larger box
        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
        status = pmSourceMoments (source, BIG_RADIUS);
        if (status) {
	    source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
            Nmoments ++;
            continue;
        }

	source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
        Nfail ++;
        psErrorClear();
        continue;
    }

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

    scalar = job->args->data[3];
    scalar->data.S32 = Nfail;
    
    return true;
}

# if (0)
bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) {

    bool status = false;
    float BIG_RADIUS;

    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    if (!status) return false;
    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    if (!status) return false;
    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    if (!status) return false;

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

    // threaded measurement of the sources moments
    int Nfail = 0;
    int Nmoments = 0;
    for (int i = 0; i < sources->n; i++) {
        pmSource *source = sources->data[i];

        // skip faint sources for moments measurement
        if (source->peak->SN < MIN_SN) {
            continue;
        }

        // measure a local sky value
        // the local sky is now ignored; kept here for reference only
        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
        if (!status) {
	    psErrorClear(); // XXX re-consider the errors raised here
	    Nfail ++;
	    continue;
        }

        // measure the local sky variance (needed if noise is not sqrt(signal))
        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
        if (!status) {
	    Nfail ++;
	    psErrorClear();
	    continue;
        }

        // measure basic source moments
        status = pmSourceMoments (source, RADIUS);
        if (status) {
            Nmoments ++;
            continue;
        }

        // if no valid pixels, or massive swing, likely saturated source,
        // try a much larger box
        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
        status = pmSourceMoments (source, BIG_RADIUS);
        if (status) {
            Nmoments ++;
            continue;
        }

        Nfail ++;
        psErrorClear();
        continue;
    }

    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    *nmoments = Nmoments;
    *nfail = Nfail;
    
    return true;
}
# endif 
