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

#include "pmHDU.h"
#include "pmFPA.h"

#include "pmSubtraction.h"
#include "pmSubtractionKernels.h"
#include "pmSubtractionMatch.h"

#define ARRAY_BUFFER 16                 // Number to add to array at a time

// Names of FITS table columns
#define NAME_XMIN "XMIN"                // Region of applicability: minimum x value
#define NAME_XMAX "XMAX"                // Region of applicability: maximum x value
#define NAME_YMIN "YMIN"                // Region of applicability: minimum y value
#define NAME_YMAX "YMAX"                // Region of applicability: maximum y value
#define NAME_KERNEL "KERNEL"            // Kernel description
#define NAME_TYPE "TYPE"                // Kernel type
#define NAME_SIZE "SIZE"                // Kernel half-size
#define NAME_INNER "INNER"              // Size of inner region (only applicable for some kernel types)
#define NAME_SPATIAL "SPATIAL_ORDER"    // Order of spatial polynomial
#define NAME_BG "BG_ORDER"              // Order of background polynomial
#define NAME_MODE "MODE"                // Matching mode
#define NAME_COLS "COLUMNS"             // Number of columns
#define NAME_ROWS "ROWS"                // Number of rows
#define NAME_SOL1 "SOLUTION_1"          // Solution for convolving image 1
#define NAME_SOL2 "SOLUTION_2"          // Solution for convolving image 2


// Generate an extension name for the subtraction kernel
static psString subtractionKernelExtname(const pmReadout *ro // Readout containing kernel
                                         )
{
    psAssert(ro, "Required parameter");

    psString extname = NULL;            // Name of FITS extension
    pmCell *cell = ro->parent;     // Parent cell
    if (cell && cell->parent) {
        pmChip *chip = cell->parent;        // Parent chip
        const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
        const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
        psAssert(chipName && cellName, "Chip and cell names should be defined in concepts.");
        int roNum = -1;                 // Index of readout
        for (int i = 0; i < cell->readouts->n && roNum == -1; i++) {
            if (cell->readouts->data[i] == ro) {
                roNum = i;
            }
        }
        if (roNum != -1) {
            psStringAppend(&extname, "%s_%s_%d", chipName, cellName, roNum);
        } else {
            psWarning("Unable to find readout number");
        }
    }

    if (!extname) {
        extname = psStringCopy("SUBTRACTION_KERNEL");
    }

    return extname;
}

bool pmReadoutWriteSubtractionKernels(pmReadout *ro, psFits *fits)
{
    PM_ASSERT_READOUT_NON_NULL(ro, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);


    // Extract the regions and solutions used in the image matching
    psArray *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
    {
        psString regex = NULL;          // Regular expression
        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
        psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD, regex); // Iterator
        psFree(regex);
        psMetadataItem *item = NULL;// Item from iteration
        while ((item = psMetadataGetAndIncrement(iter))) {
            assert(item->type == PS_DATA_REGION);
            regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
        }
        psFree(iter);
    }
    if (regions->n == 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No subtraction regions found.");
        psFree(regions);
        return false;
    }

    psArray *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
    {
        psString regex = NULL;          // Regular expression
        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
        psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD, regex); // Iterator
        psFree(regex);
        psMetadataItem *item = NULL;// Item from iteration
        while ((item = psMetadataGetAndIncrement(iter))) {
            assert(item->type == PS_DATA_UNKNOWN);
            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
            // images by scans.
            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
            kernel->numCols = ro->image->numCols;
            kernel->numRows = ro->image->numRows;

            kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
        }
        psFree(iter);
    }

    if (regions->n != kernels->n) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of regions (%ld) and kernels (%ld) don't match.\n",
                regions->n, kernels->n);
        psFree(regions);
        psFree(kernels);
        return false;
    }

    // Format for writing a table
    int num = regions->n;              // Number of regions and kernels
    psArray *rows = psArrayAlloc(num); // Array of FITS table rows
    for (int i = 0; i < num; i++) {
        psMetadata *row = psMetadataAlloc(); // Row of interest
        rows->data[i] = psMemIncrRefCounter(row);

        psRegion *region = regions->data[i]; // Region of interest
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_XMIN,  0, "Applicability minimum x", region->x0);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_XMAX,  0, "Applicability maximum x", region->x1);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_YMIN,  0, "Applicability minimum y", region->y0);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_YMAX,  0, "Applicability maximum y", region->y1);

        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel
        psMetadataAddStr(row, PS_LIST_TAIL, NAME_KERNEL, 0, "Kernel description", kernel->description);

#if 0
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_TYPE,  0, "Kernel type (enum)", kernel->type);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_SIZE,  0, "Kernel half-size", kernel->size);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_INNER, 0, "Size of inner region", kernel->inner);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_SPATIAL, 0, "Polynomial order for spatial variations",
                         kernel->spatialOrder);
#endif

        psMetadataAddS32(row, PS_LIST_TAIL, NAME_BG,  0, "Polynomial order for background fitting",
                         kernel->bgOrder);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_MODE,  0, "Matching mode (enum)", kernel->mode);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_COLS,  0, "Number of columns", kernel->numCols);
        psMetadataAddS32(row, PS_LIST_TAIL, NAME_ROWS,  0, "Number of rows", kernel->numRows);
        if (kernel->mode == PM_SUBTRACTION_MODE_1 || kernel->mode == PM_SUBTRACTION_MODE_DUAL) {
            psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL1, 0, "Solution vector 1", kernel->solution1);
        }
        if (kernel->mode == PM_SUBTRACTION_MODE_2 || kernel->mode == PM_SUBTRACTION_MODE_DUAL) {
            psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL2, 0, "Solution vector 2", kernel->solution2);
        }
    }
    psFree(regions);
    psFree(kernels);

    psMetadata *header = psMetadataAlloc(); // Header for FITS file

    // CVS tags, used to identify the version of this file (in case incompatibilities are introduced)
    psString cvsFile = psStringCopy("$RCSfile: pmSubtractionIO.c,v $");
    psString cvsRev  = psStringCopy("$Revision: 1.1 $");
    psString cvsDate = psStringCopy("$Date: 2008-06-16 21:59:43 $");
    psStringSubstitute(&cvsFile, NULL, "RCSfile: ");
    psStringSubstitute(&cvsRev,  NULL, "Revision: ");
    psStringSubstitute(&cvsDate, NULL, "Date: ");

    psString version = NULL;            // Version information, for header
    psStringAppend(&version, "%s %s %s", cvsFile, cvsRev, cvsDate);
    psFree(cvsFile);
    psFree(cvsRev);
    psFree(cvsDate);
    psStringSubstitute(&version, NULL, "$");
    psMetadataAddStr(header, PS_LIST_TAIL, "PSVERSION", 0, "S/W version", version);
    psFree(version);

    psString extname = subtractionKernelExtname(ro); // Extension name

    if (!psFitsWriteTable(fits, header, rows, extname)) {
        psError(PS_ERR_IO, false, "Unable to write subtraction kernel to FITS table.");
        psFree(header);
        psFree(rows);
        psFree(extname);
        return false;
    }

    psFree(header);
    psFree(rows);
    psFree(extname);

    return true;
}


bool pmReadoutReadSubtractionKernels(pmReadout *ro, psFits *fits)
{
    PM_ASSERT_READOUT_NON_NULL(ro, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    psString extname = subtractionKernelExtname(ro); // Extension name
    if (!psFitsMoveExtName(fits, extname)) {
        psError(PS_ERR_IO, false, "Unable to move to subtraction kernel table.");
        psFree(extname);
        return false;
    }
    psFree(extname);

    psArray *table = psFitsReadTable(fits); // Table of interest
    if (!table) {
        psError(PS_ERR_IO, false, "Unable to read FITS table");
        return false;
    }

    // Look up a column value for a row
#define TABLE_LOOKUP(TYPE, SUFFIX, TARGET, NAME) \
    TYPE TARGET; \
    { \
        bool mdok; \
        TARGET = psMetadataLookup##SUFFIX(&mdok, row, NAME); \
        if (!mdok) { \
            psError(PS_ERR_UNKNOWN, false, "Unable to find column %s in subtraction kernel table.", NAME); \
            psFree(table); \
            return false; \
        } \
    }

    for (int i = 0; i < table->n; i++) {
        psMetadata *row = table->data[i]; // Table row

        TABLE_LOOKUP(int, S32, xMin, NAME_XMIN);
        TABLE_LOOKUP(int, S32, xMax, NAME_XMAX);
        TABLE_LOOKUP(int, S32, yMin, NAME_YMIN);
        TABLE_LOOKUP(int, S32, yMax, NAME_YMAX);

        psRegion *region = psRegionAlloc(xMin, xMax, yMin, yMax); // Region of applicability
        psMetadataAddPtr(ro->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_DATA_REGION,
                         "Subtraction region", region);
        psFree(region);

        TABLE_LOOKUP(const char *, Str, description, NAME_KERNEL);
        TABLE_LOOKUP(pmSubtractionMode, S32, mode,    NAME_MODE);

#if 0
        TABLE_LOOKUP(int, S32, size,    NAME_SIZE);
        TABLE_LOOKUP(int, S32, inner,   NAME_INNER);
        TABLE_LOOKUP(int, S32, spatial, NAME_SPATIAL);
#endif

        TABLE_LOOKUP(int, S32, bg,      NAME_BG);
        TABLE_LOOKUP(int, S32, numCols, NAME_COLS);
        TABLE_LOOKUP(int, S32, numRows, NAME_ROWS);

        pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode);
        kernels->numCols = numCols;
        kernels->numRows = numRows;

        bool mdok;                      // Status of MD lookup
        if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) {
            kernels->solution1 = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, row, NAME_SOL1));
            if (!mdok) {
                psError(PS_ERR_UNKNOWN, false, "Unable to find column %s in subtraction kernel table.",
                        NAME_SOL1);
                psFree(kernels);
                psFree(table);
                return false;
            }
        }
        if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) {
            kernels->solution2 = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, row, NAME_SOL2));
            if (!mdok) {
                psError(PS_ERR_UNKNOWN, false, "Unable to find column %s in subtraction kernel table.",
                        NAME_SOL2);
                psFree(kernels);
                psFree(table);
                return false;
            }
        }

        psMetadataAddPtr(ro->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, PS_DATA_UNKNOWN,
                         "Subtraction kernels", kernels);
        psFree(kernels);
    }
    psFree(table);
    return true;
}
