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

#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <assert.h>
#include <pslib.h>
#include <psmodules.h>
#include <ppStats.h>

#include "ppMerge.h"
#include "ppMergeData.h"
#include "ppMergeCombine.h"
#include "ppMergeVersion.h"

#if 0
static psMemId memId = 0;

static FILE *dumpFile = NULL;

static psMemId mbAlloc(psMemBlock *mb)
{
    if (!dumpFile) {
        dumpFile = fopen("memBlocks.dat", "w");
    }
    fprintf(dumpFile, "Alloc: %12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
            mb->file, mb->lineno);
    return 1;
}

static void dumpDone(void)
{
    fclose(dumpFile);
    exit(EXIT_FAILURE);
}


static void memDump(void)
{
    psMemBlock **leaks = NULL;
    int numLeaks = psMemCheckLeaks(memId, &leaks, NULL, true);
    FILE *memFile = fopen("mem.dat", "w");
    fprintf(memFile, "# MemBlock Size Source\n");
    for (int i = 0; i < numLeaks; i++) {
        psMemBlock *mb = leaks[i];
        fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
                mb->file, mb->lineno);
    }
    fclose(memFile);
    psFree(leaks);
}

static void memCheck(void)
{
    return;
    if (psTraceGetLevel("ppMerge") > 9) {
        psMemBlock **leaks = NULL;
        int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
        size_t largestSize = 0;
        psMemId largest = 0;
        size_t totalSize = 0;
        for (int i = 0; i < numLeaks; i++) {
            psMemBlock *mb = leaks[i];
            totalSize += mb->userMemorySize;
            if (mb->userMemorySize > largestSize) {
                largestSize = mb->userMemorySize;
                largest = mb->id;
            }
        }
        psFree(leaks);
        psTrace("ppMerge", 0, "Memory in use: %zd\n", totalSize);
        psTrace("ppMerge", 0, "Largest block: %ld\n", largest);
        psTrace("ppMerge", 0, "sbrk(): %zd\n", sbrk(0));
    }
    return;
}
#endif


// Combine the inputs
bool ppMergeCombine(psImage *scales,    // Scales for each cell of each integration, or NULL
                    psImage *zeros,     // Zeros for each cell of each integration, or NULL
                    ppMergeData *data,  // Data
                    ppMergeOptions *options, // Options
                    pmConfig *config    // Configuration
    )
{
    psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
    assert(filenames);                  // It should be here --- it's put here in ppMergeConfig

    // Sanity checks
    assert(!options->scale || scales);
    assert(!scales || (scales->type.type == PS_TYPE_F32 &&
                       scales->numCols == data->numCells &&
                       scales->numRows == filenames->n));
    assert(!options->zero || zeros);
    assert(!zeros || (zeros->type.type == PS_TYPE_F32 &&
                      zeros->numCols == data->numCells &&
                      zeros->numRows == filenames->n));

    // Iterate over the FPA
    pmFPA *fpa = data->out;             // Output FPA
    pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration
    int cellNum = -1;                   // Cell number in the whole FPA
    pmFPAWrite(data->out, data->outFile, config->database, true, false); // Write header only
    pmChip *chip;                       // Chip of interest
    psRandom *rng = NULL;               // Random number generator; required for building a mask
    pmHDU *lastHDU = NULL;              // Last HDU to be updated
    if (options->mask) {
        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    }
    while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
        pmChipWrite(chip, data->outFile, config->database, true, false); // Write header only
        pmCell *cell;                   // Cell of interest
        while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
            cellNum++;
            pmCellWrite(cell, data->outFile, config->database, true); // Write header only
            pmReadout *readout = pmReadoutAlloc(cell); // Output readout of interest
            psArray *stack = psArrayAlloc(filenames->n); // Stack of readouts to combine
            psVector *cellScales = NULL; // Scales for this cell
            if (scales) {
                cellScales = psImageCol(NULL, scales, cellNum);
            }
            psVector *cellZeros = NULL;  // Zeros for this cell
            if (zeros) {
                cellZeros = psImageCol(NULL, zeros, cellNum);
            }

            // Read bit by bit
            int numRead;  // Number of inputs read
            int numScan = 0;

            // Put version metadata into header
            pmHDU *hdu = pmHDUFromCell(cell);
            if (hdu && hdu != lastHDU) {
                if (!hdu->header) {
                    hdu->header = psMetadataAlloc();
                }
                ppMergeVersionMetadata(hdu->header);
                lastHDU = hdu;
            }

            do {
                numRead = 0;
                for (int i = 0; i < filenames->n; i++) {
                    if (! filenames->data[i] || strlen(filenames->data[i]) == 0) {
                        continue;
                    }
                    psFits *fits = data->files->data[i]; // FITS file handle
                    if (!fits) {
                        continue;
                    }

                    if (!stack->data[i]) {
                        pmFPA *fpaIn = data->in->data[i]; // Input FPA
                        pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip
                        pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell
                        stack->data[i] = pmReadoutAlloc(cellIn); // Input readout
                    }

                    // Only reading and writing the first readout in each cell (plane 0)
                    if (pmReadoutReadNext(stack->data[i], fits, 0, options->rows)) {
                        // If the mask value is set, we want to generate a mask
                        if (options->combine->maskVal) {
                            pmReadoutSetMask(stack->data[i]);
                        }

                        // If we're combining with weights, we want to generate weights.
                        if (options->combine->weights) {

                            // If it's a bias or dark, set the gain to zero: noise only contributed by read
                            if ((!options->zero && !options->scale) || options->darktime) {
                                pmReadoutSetWeight(stack->data[i], false);
                            } else {
                                pmReadoutSetWeight(stack->data[i], true);
                            }
                        }

                        numRead++;
                    } else {
                        psTrace("ppMerge", 3, "Unable to read from file %d for chip %d, "
                                "cell %d, scan %d\n", i, view->chip, view->cell, numScan);
                    }

                }

                psTrace("ppMerge", 5, "Chip %d, cell %d, scan %d\n", view->chip, view->cell, numScan);
                if (numRead > 0) {

                    if (options->shutter) {
                        pmShutterCorrectionMeasure(readout, stack, options->shutterSize,
                                                   options->mean, options->stdev,
                                                   options->shutterIter,
                                                   options->shutterRej,
                                                   options->combine->maskVal);
                    } else {
                        pmReadoutCombine(readout, stack, cellZeros, cellScales, options->combine);
                    }
                }
                numScan++;

            } while (numRead > 0);

            psFree(stack);

            // Set the dark time for the output image, since we normalised
            if (options->darktime) {
                psMetadataItem *darkItem = psMetadataLookup(cell->concepts, "CELL.DARKTIME");
                darkItem->data.F32 = 1.0;
                psMetadataItem *expItem = psMetadataLookup(cell->concepts, "CELL.EXPOSURE");
                expItem->data.F32 = 1.0;
            }

            // Measure the fringes for this cell
            //
            // XXX Need to deal with multiple components: we will do this by building up the components
            // Read the existing fringe measurements
            // Use existing regions to measure fringe statistics
            // Add the new fringe measurements to the existing fringe measurements.
            // Write the appended fringe measurements.
            // Read in the "output" file to get the existing components.
            // Put the new readout into the cell after the existing readouts.
            if (options->fringe && readout->image) {
                pmFringeRegions *regions = pmFringeRegionsAlloc(options->fringeNum, options->fringeSize,
                                                                options->fringeSize, options->fringeSmoothX,
                                                                options->fringeSmoothY); // Fringe regions
                pmFringeStats *fringe = pmFringeStatsMeasure(regions, readout, options->combine->maskVal);
                psFree(regions);
                if (!fringe) {
                    psError(PS_ERR_UNKNOWN, false, "Unable to measure fringe statistics.\n");
                    psFree(readout);
                    return false;
                }

                psArray *fringes = psArrayAlloc(1); // Array of fringes
                fringes->data[0] = fringe;

                pmFringesFormat(cell, NULL, fringes);
                psFree(fringes);        // Drop reference
            }

            // Add MD5 information for cell
            if (readout->image) {
                pmHDU *hdu = pmHDUFromCell(cell); // HDU that owns the cell
                const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
                const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");

                psString headerName = NULL; // Header name for MD5
                psStringAppend(&headerName, "MD5_%s_%s", chipName, cellName);

                psVector *md5 = psImageMD5(readout->image); // md5 hash
                psString md5string = psMD5toString(md5); // String
                psFree(md5);
                psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
                                 "Image MD5", md5string);
                psFree(md5string);
                psFree(headerName);
            }

            psFree(readout);            // Drop reference

            // Statistics on the merged cell
            if (data->statsFile) {
                data->stats = ppStats(data->stats, data->out, view, config); // Statistics for output FPA
            }

            // We threw away the bias sections --- record this
            psMetadataItem *biassecItem = psMetadataLookup(cell->concepts, "CELL.BIASSEC"); // Item of BIASSEC
            psList *biassecList = biassecItem->data.V; // List of BIASSECs
            while (psListRemove(biassecList, PS_LIST_TAIL)); // Removing all entries

            // Blow away the cell data
            for (int i = 0; i < filenames->n; i++) {
                pmFPA *fpaIn = data->in->data[i]; // Input FPA
                pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip
                pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell
                pmCellFreeData(cellIn);
            }

            // Write the pixels
            if (cell->hdu && !cell->hdu->blankPHU) {
                psTrace("ppMerge", 5, "Writing out cell HDU.\n");
                pmCellWrite(cell, data->outFile, config->database, false);
                if (options->fringe) {
                    pmCellWriteTable(data->outFile, cell, "FRINGE");
                }

                pmCellFreeData(cell);
            }
        }

        // Blow away the chip data
        for (int i = 0; i < filenames->n; i++) {
            pmFPA *fpaIn = data->in->data[i]; // Input FPA
            pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip
            pmChipFreeData(chipIn);
        }

        // Write the pixels
        if (chip->hdu && !chip->hdu->blankPHU) {
            psTrace("ppMerge", 5, "Writing out chip HDU.\n");
            pmChipWrite(chip, data->outFile, config->database, false, false);
            if (options->fringe) {
                pmChipWriteTable(data->outFile, chip, "FRINGE");
            }

            pmChipFreeData(chip);
        }
    }

    // Blow away the FPA data
    for (int i = 0; i < filenames->n; i++) {
        pmFPA *fpaIn = data->in->data[i]; // Input FPA
        pmFPAFreeData(fpaIn);
    }

    if (data->out->hdu && !data->out->hdu->blankPHU) {
        // Write the pixels
        psTrace("ppMerge", 5, "Writing out FPA HDU.\n");
        pmFPAWrite(data->out, data->outFile, config->database, false, false);
        if (options->fringe) {
            pmFPAWriteTable(data->outFile, fpa, "FRINGE");
        }
    }

    pmFPAFreeData(data->out);

    psFree(view);
    psFree(rng);

    return true;

}
