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

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

#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPALevel.h"
#include "pmConcepts.h"
#include "pmConceptsAverage.h"

// Update a metadata entry directly
#define MD_UPDATE(MD, NAME, TYPE, VALUE) \
{ \
    psMetadataItem *item = psMetadataLookup(MD, NAME); \
    item->data.TYPE = VALUE; \
}

// Update a metadata string entry directly
#define MD_UPDATE_STR(MD, NAME, VALUE) \
{ \
    psMetadataItem *item = psMetadataLookup(MD, NAME); \
    psFree(item->data.str); \
    item->data.str = psStringCopy(VALUE); \
}

bool pmConceptsAverageFPAs(pmFPA *target, psList *sources)
{
    PS_ASSERT_PTR_NON_NULL(target, false);
    PS_ASSERT_PTR_NON_NULL(sources, false);
    PS_ASSERT_INT_POSITIVE(sources->n, false);

    double time      = 0.0;             // Time of observation
    psTimeType timeSys = 0;             // Time system
    char *filter     = NULL;            // Filter

    int num = 0;                        // Number of FPAs
    psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources
    pmFPA *fpa = NULL;                  // Source FPA from iteration
    while ((fpa = psListGetAndIncrement(sourcesIter))) {
        if (!fpa) {
            continue;
        }

        num++;

        psTime *fpaTime = psMetadataLookupPtr(NULL, fpa->concepts, "FPA.TIME");
        time       += psTimeToMJD(fpaTime);
        if (num == 1) {
            timeSys = psMetadataLookupS32(NULL, fpa->concepts, "FPA.TIMESYS");
            filter = psMetadataLookupStr(NULL, fpa->concepts, "FPA.FILTER");
        } else {
            if (timeSys != psMetadataLookupS32(NULL, fpa->concepts, "FPA.TIMESYS")) {
                psWarning("Differing FPA.TIMESYS in use: %d vs %d\n",
                          timeSys, psMetadataLookupS32(NULL, fpa->concepts, "FPA.TIMESYS"));
            }
            if (strcmp(filter, psMetadataLookupStr(NULL, fpa->concepts, "FPA.FILTER")) != 0) {
                psWarning("Differing FPA.FILTER in use: %s vs %s\n",
                          filter, psMetadataLookupStr(NULL, fpa->concepts, "FPA.FILTER"));
            }
        }
    }
    psFree(sourcesIter);

    time      /= (double)num;

    MD_UPDATE(target->concepts, "FPA.TIMESYS", S32, timeSys);
    MD_UPDATE_STR(target->concepts, "FPA.FILTER", filter);

    // FPA.TIME needs special care
    {
        psMetadataItem *timeItem = psMetadataLookup(target->concepts, "FPA.TIME");
        psFree(timeItem->data.V);
        timeItem->data.V = psTimeFromMJD(time);
    }

    return true;
}


// Set a variety of concepts in a cell by averaging over several
// XXX does not properly set XSIZE, YSIZE
bool pmConceptsAverageCells(pmCell *target, psList *sources, psRegion *trimsec, psRegion *biassec, bool same)
{
    PS_ASSERT_PTR_NON_NULL(target, false);
    PS_ASSERT_PTR_NON_NULL(sources, false);
    PS_ASSERT_INT_POSITIVE(sources->n, false);

    float gain       = 0.0;             // Gain
    float readnoise  = 0.0;             // Read noise
    float saturation = INFINITY;        // Saturation level
    float bad        = -INFINITY;       // Bad level
    float exposure   = 0.0;             // Exposure time
    float darktime   = 0.0;             // Dark time
    double time      = 0.0;             // Time of observation
    psTimeType timeSys = 0;             // Time system
    int readdir      = 0;               // Cell read direction
    int xBin = 0, yBin = 0;             // Binning
    int x0 = 0, y0 = 0;                 // Offset
    int xParity = 0, yParity = 0;       // Parity

    int nCells = 0;                     // Number of cells;
    psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources
    pmCell *cell = NULL;                // Source cell from iteration
    while ((cell = psListGetAndIncrement(sourcesIter))) {
        if (!cell) {
            continue;
        }

        nCells++;
        gain       += psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN");
        readnoise  += psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");
        exposure   += psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE");
        darktime   += psMetadataLookupF32(NULL, cell->concepts, "CELL.DARKTIME");
        psTime *cellTime = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME");
        time       += psTimeToMJD(cellTime);
        if (nCells == 1) {
            timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS");
            readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR");
            xBin    = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
            yBin    = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN");

            if (same) {
                // Only makes sense to update these if they are the same cell
                x0 = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
                y0 = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
                xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
                yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
            }
        } else {
            if (timeSys != psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS")) {
                psWarning("Differing CELL.TIMESYS in use: %d vs %d\n",
                          timeSys, psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"));
            }
            if (readdir != psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR")) {
                psWarning("Differing CELL.READDIR in use: %d vs %d\n",
                          readdir, psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"));
            }
            if (xBin != psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN")) {
                psWarning("Differing CELL.XBIN in use: %d vs %d\n",
                          xBin, psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"));
            }
            if (yBin != psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN")) {
                psWarning("Differing CELL.YBIN in use: %d vs %d\n",
                          yBin, psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"));
            }
            if (same) {
                if (x0 != psMetadataLookupS32(NULL, cell->concepts, "CELL.X0")) {
                    psWarning("Differing CELL.X0 in use: %d vs %d\n",
                              x0, psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"));
                }
                if (y0 != psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0")) {
                    psWarning("Differing CELL.Y0 in use: %d vs %d\n",
                              y0, psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"));
                }
                if (xParity != psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY")) {
                    psWarning("Differing CELL.XPARITY in use: %d vs %d\n",
                              xParity, psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"));
                }
                if (yParity != psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY")) {
                    psWarning("Differing CELL.YPARITY in use: %d vs %d\n",
                              yParity, psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"));
                }
            }
        }

        float cellSaturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
        if (cellSaturation < saturation) {
            saturation = cellSaturation;
        }
        float cellBad = psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD");
        if (cellBad > bad) {
            bad = cellBad;
        }
    }
    psFree(sourcesIter);

    gain      /= (float)nCells;
    readnoise /= (float)nCells;
    exposure  /= (float)nCells;
    darktime  /= (float)nCells;
    time      /= (double)nCells;

    MD_UPDATE(target->concepts, "CELL.GAIN", F32, gain);
    MD_UPDATE(target->concepts, "CELL.READNOISE", F32, readnoise);
    MD_UPDATE(target->concepts, "CELL.SATURATION", F32, saturation);
    MD_UPDATE(target->concepts, "CELL.BAD", F32, bad);
    MD_UPDATE(target->concepts, "CELL.EXPOSURE", F32, exposure);
    MD_UPDATE(target->concepts, "CELL.DARKTIME", F32, darktime);
    MD_UPDATE(target->concepts, "CELL.TIMESYS", S32, timeSys);
    MD_UPDATE(target->concepts, "CELL.READDIR", S32, readdir);
    MD_UPDATE(target->concepts, "CELL.XBIN", S32, xBin);
    MD_UPDATE(target->concepts, "CELL.YBIN", S32, yBin);
    if (same) {
        MD_UPDATE(target->concepts, "CELL.X0", S32, x0);
        MD_UPDATE(target->concepts, "CELL.Y0", S32, y0);
        MD_UPDATE(target->concepts, "CELL.XPARITY", S32, xParity);
        MD_UPDATE(target->concepts, "CELL.YPARITY", S32, yParity);
    }

    // CELL.TIME needs special care
    {
        psMetadataItem *timeItem = psMetadataLookup(target->concepts, "CELL.TIME");
        psFree(timeItem->data.V);
        timeItem->data.V = psTimeFromMJD(time);
    }

    // CELL.TRIMSEC needs special care
    if (trimsec) {
        psMetadataItem *trimsecItem = psMetadataLookup(target->concepts, "CELL.TRIMSEC");
        psFree(trimsecItem->data.V);
        trimsecItem->data.V = psMemIncrRefCounter(trimsec);
    }

    // CELL.BIASSEC needs special care
    if (biassec) {
        psMetadataItem *biassecItem = psMetadataLookup(target->concepts, "CELL.BIASSEC");
        psFree(biassecItem->data.V);
        biassecItem->data.V = psMemIncrRefCounter(biassec);
    }

    return true;
}

