#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>

#include "ppStack.h"
#include "ppStackLoop.h"

bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
{
    psAssert(stack, "Require stack");
    psAssert(options, "Require options");
    psAssert(config, "Require configuration");

    psTimerStart("PPSTACK_FINAL");

    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    psAssert(recipe, "We've thrown an error on this before.");

    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    int overlap = 2 * psMetadataLookupS32(NULL, ppsub,
                                          "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans

    pmFPAview *view = NULL;             // View to readout

    int row0, col0;                 // Offset for readout
    int numCols, numRows;           // Size of readout
    ppStackThread *thread = stack->threads->data[0]; // Representative thread
    if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
        psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
        return false;
    }

    pmFPAfileActivate(config->files, false, NULL);
    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true);
    view = ppStackFilesIterateDown(config);
    if (!view) {
        return false;
    }

    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
    options->outRO = pmReadoutAlloc(outCell); // Output readout
    psFree(view);

    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
        return false;
    }

    psFree(options->cells); options->cells = NULL;

    bool status;                    // Status of read
    int numChunk;                   // Number of chunks
    for (numChunk = 0; true; numChunk++) {
        ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap);
        if (!status) {
            // Something went wrong
            psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
            return false;
        }
        if (!thread) {
            // Nothing more to read
            break;
        }

        // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
        psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
        psArrayAdd(job->args, 1, thread);
        psArrayAdd(job->args, 1, options);
        psArrayAdd(job->args, 1, config);
        if (!psThreadJobAddPending(job)) {
            psFree(job);
            return false;
        }
        psFree(job);
    }

    if (!psThreadPoolWait(false)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination.");
        return false;
    }

    // Harvest the jobs, gathering the inspection lists
    options->inspect = psArrayAlloc(options->num);
    for (int i = 0; i < options->num; i++) {
        if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
            continue;
        }
        options->inspect->data[i] = psArrayAllocEmpty(numChunk);
    }
    psThreadJob *job;               // Completed job
    while ((job = psThreadJobGetDone())) {
        psAssert(strcmp(job->type, "PPSTACK_INITIAL_COMBINE") == 0,
                 "Job has incorrect type: %s", job->type);
        psArray *results = job->results; // Results of job
        for (int i = 0; i < options->num; i++) {
            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
                continue;
            }
            options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, results->data[i]);
        }
        psFree(job);
    }

    ppStackMemDump("initial");

#ifdef TESTING
    ppStackWriteImage("combined_initial.fits", NULL, outRO->image, config);
    pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits");
#endif

    // Sum covariance matrices
    for (int i = 0; i < options->num; i++) {
        if (options->inputMask->data.U8[i]) {
            psFree(options->covariances->data[i]);
            options->covariances->data[i] = NULL;
        }
    }
    options->outRO->covariance = psImageCovarianceSum(options->covariances);
    psFree(options->covariances); options->covariances = NULL;
    psFree(options->matchChi2); options->matchChi2 = NULL;


    if (options->stats) {
        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_FINAL", 0,
                         "Time to make final stack", psTimerMark("PPSTACK_FINAL"));
    }

    return true;
}
