/** @file ippToPspsBatchDetection.c
 *
 *  @ingroup ippToPsps
 *
 *  @author IfA
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

#include <limits.h>
#include "ippToPsps.h"
#include "ippToPspsConfig.h"
#include "ippToPspsTestEnums.h"
#include "fitsio.h"

// Gets flux from magnitude
static __inline bool ippToPsps_getFlux(const float zeroPoint, const float exposureTime, const float magnitude, float* flux) {

    // use uncalibrated instrumental flux
    *flux = powf(10.0, -0.4*magnitude) / exposureTime;

    // use calibrated flux in Janskys, where 3631 Jy is the zero point flux for the filter: constant over all filters because we're using AB magnitudes
    //float flux = 3631 * powf(10.0, -0.4*(magnitude + zeroPoint)) / exposureTime;

    //    printf("mag=%f\texpTime=%f\tzeroPoint=%f\tflux=%f\n", magnitude, exposureTime, zeroPoint, flux);
    return (!isfinite(*flux) || *flux < 0.000001) ? false : true;
}

/**
  \brief test data for checking PSPS values correspond with IPP 

  */
int ippToPsps_batchTest(IppToPsps *this) {

    if (this->numOfInputFiles < 1) return PS_EXIT_DATA_ERROR;

    int status = 0;
    fitsfile *fitsIn;         
    int nKeys;

    if (fits_open_file(&fitsIn, this->inputFiles[0], READONLY, &status)) {

        fits_report_error(stderr, status);
        return PS_EXIT_SYS_ERROR;
    }

    // get primary header and pull stuff out for later
    fits_get_hdrspace(fitsIn, &nKeys, NULL, &status);

    float zptObs, zeroPoint, exposureTime;
    char filterType[20];
    status=0; fits_read_key(fitsIn, TFLOAT, "ZPT_OBS", &zptObs, NULL, &status);
    status=0; fits_read_key(fitsIn, TFLOAT, "ZPT_REF", &zeroPoint, NULL, &status);
    status=0; fits_read_key(fitsIn, TFLOAT, "EXPREQ", &exposureTime, NULL, &status);
    status=0; fits_read_key(fitsIn, TSTRING, "FILTERID", filterType, NULL, &status);
    
    ippToPspsConfig_writeTable(this->config, fitsIn, this->fitsOut, 1, "FrameMeta", true);

    // FrameMeta values
    fits_write_col(this->fitsOut, TLONG, FRAMEMETA_FRAMEID, 1, 1, 1, &this->expId, &status);

    int8_t filterID = -1;
    ippToPspsConfig_getFilterId(this->config, filterType, &filterID);
    fits_write_col(this->fitsOut, TBYTE, FRAMEMETA_FILTERID, 1, 1, 1, &filterID, &status);

    // stuff to keep from psf.hdr header
    double obsTime = 0.0;
    int sourceId = -1;
    int imageId = -1;
    float fwhmMaj;
    float fwhmMin;
    float momentMaj;
    float momentMin;
    float psfFwhm;
    float momentFwhm;
    long pspsImageId = -1;

    // DVO variables
    SkyList *skylist = NULL;
    dvoDetection *dvoDetections = NULL;
    int maxDvoDetId = -1;
    int numDvoDetections = -1;
    Image *pImage = NULL;

    // stuff for detections table
    time_t now;
    struct tm *ts;
    char timeStr[12];
    now = time(NULL);
    ts = gmtime(&now); // TODO should be UTC
    strftime(timeStr, sizeof(timeStr), "%Y-%m-%d", ts);

    uint32_t totalDetections = 0;
    uint32_t s,d;
    uint32_t invalidDvoRows;
    uint32_t smfJumps;
    uint32_t unmatched;

    long longnull = -999;
    float floatnull = -999.0;
    int anynull = 0;

    char ccdNumber[3];
    char extensionName[15];

    long ippIDet[MAXDETECT];
    float instMag[MAXDETECT];
    float instMagErr[MAXDETECT];
    float peakMag[MAXDETECT];
    long objID[MAXDETECT];
    long detectID[MAXDETECT];
    long ippObjID[MAXDETECT];
    long ippDetectID[MAXDETECT];
    long imageID[MAXDETECT];
    float obsTimes[MAXDETECT];
    float instFlux[MAXDETECT];
    float instFluxErr[MAXDETECT];
    float peakFlux[MAXDETECT];
    int8_t filterIDs[MAXDETECT];

    char** assocDate = (char**)calloc(MAXDETECT, sizeof(char**));
    for (uint32_t i=0; i<MAXDETECT;i++) assocDate[i] = (char*)calloc(20,sizeof(char));

    // some stuff is the same for all detections so we can populate here
    for (long s = 0; s<MAXDETECT; s++) {

        filterIDs[s] = filterID;
        strcpy(assocDate[s], timeStr);
    }

    long maxObjID = LONG_MIN; 
    long minObjID = LONG_MAX;
    short nOta = 0;
    long i;
    bool isDuplicate;
    uint32_t numOfDuplicates;
    uint32_t numInvalidFlux;
    long numDetectionsOut;
    long removeList[MAXDETECT];
    long thisObjId;

    // loop round all 60 chips
    for (int x=3; x<4; x++) {
        for (int y=3; y<4; y++) {

            // dodge the corners
            if (x==0 && y==0) continue;
            if (x==0 && y==7) continue;
            if (x==7 && y==0) continue;
            if (x==7 && y==7) continue;

            sprintf(ccdNumber, "%d%d", x, y);

            // check we can move to detections table in smf
            sprintf(extensionName, "XY%s.psf", ccdNumber);
            status=0;
            if (fits_movnam_hdu(fitsIn, BINARY_TBL, extensionName, 0, &status)) {
                psError(PS_ERR_IO, false, "Can't move to extension: %s skipping this chip\n", extensionName);
                continue;
            }

            // move to header extension
            sprintf(extensionName, "XY%s.hdr", ccdNumber);
            status=0;
            if (fits_movnam_hdu(fitsIn, IMAGE_HDU, extensionName, 0, &status)) {
                psError(PS_ERR_IO, false, "Can't move to extension: %s skipping this chip\n", extensionName);
                continue;
            }

            // stuff to save from psf.hdr
            status=0; fits_read_key(fitsIn, TFLOAT, "FWHM_MAJ", &fwhmMaj, NULL, &status);
            status=0; fits_read_key(fitsIn, TFLOAT, "FWHM_MIN", &fwhmMin, NULL, &status);
            status=0; fits_read_key(fitsIn, TFLOAT, "IQ_FW1", &momentMaj, NULL, &status);
            status=0; fits_read_key(fitsIn, TFLOAT, "IQ_FW2", &momentMin, NULL, &status);
            status=0; fits_read_key(fitsIn, TDOUBLE, "MJD-OBS", &obsTime, NULL, &status);

            status=0; fits_read_key(fitsIn, TLONG, "IMAGEID", &imageId, NULL, &status);
            status=0; fits_read_key(fitsIn, TLONG, "SOURCEID", &sourceId, NULL, &status);

            // access DVO database
            skylist = dvoSkyListByExternID(this->dvoConfig, sourceId, imageId, &pImage);
            if (skylist == NULL) {
                psError(PS_ERR_IO, false, "DVO: can't find SkyList for sourceId='%d' imageId='%d' (CCD = XY%s): skipping\n", sourceId, imageId, ccdNumber);
                continue;
            }

            // create unique int from 'frameID' (aka exposure ID) and ccd number
            pspsImageId = (this->expId*100) + pImage->ccdnum;

            // now get DVO detections
            dvoDetections = NULL;
            numDvoDetections = dvoGetDetections(skylist, pImage->imageID, &dvoDetections, &maxDvoDetId);

            // TODO check nDet < MAXDETECT

            // create ImageMeta
            ippToPspsConfig_writeTable(this->config, fitsIn, this->fitsOut, 1, "ImageMeta", true);
            psfFwhm = (fwhmMaj+fwhmMin)/2;
            momentFwhm = (momentMaj+momentMin)/2;
            fits_write_col(this->fitsOut, TLONG, IMAGEMETA_PHOTOCALID, 1, 1, 1, &pImage->photcode, &status); 
            fits_write_col(this->fitsOut, TBYTE, IMAGEMETA_FILTERID, 1, 1, 1, &filterID, &status);
            fits_write_col(this->fitsOut, TFLOAT, IMAGEMETA_PHOTOSCAT, 1, 1, 1, &zptObs, &status);
            fits_write_col(this->fitsOut, TFLOAT, IMAGEMETA_PSFFWHM, 1, 1, 1, &psfFwhm, &status);
            fits_write_col(this->fitsOut, TFLOAT, IMAGEMETA_PHOTOZERO, 1, 1, 1, &zptObs, &status);

            // now move BACK to detections table in smf
            sprintf(extensionName, "XY%s.psf", ccdNumber);
            status=0;
            if (fits_movnam_hdu(fitsIn, BINARY_TBL, extensionName, 0, &status)) {
                psError(PS_ERR_IO, false, "Can't move to extension: %s skipping\n", extensionName);
                continue;
            }

            // keep a running count of 'images' we find in order to write total to FrameMeta at the end
            nOta++;

            long nDet = 0;
            if (fits_get_num_rows(fitsIn, &nDet, &status)) {
                fits_report_error(stderr, status);
            }

            // loop round detections to populate some extra stuff
            s = d = totalDetections = invalidDvoRows = smfJumps = unmatched = 0;

            anynull = 0;
            fits_read_col(fitsIn, TLONG, 1, 1, 1, nDet, &longnull, ippIDet, &anynull, &status);
            fits_read_col(fitsIn, TFLOAT, 8, 1, 1, nDet, &floatnull, instMag, &anynull, &status);
            fits_read_col(fitsIn, TFLOAT, 9, 1, 1, nDet, &floatnull, instMagErr, &anynull, &status);
            fits_read_col(fitsIn, TFLOAT, 12, 1, 1, nDet, &floatnull, peakMag, &anynull, &status);

            // DVO detections are ordered by IPP_IDET, which increments from 0 in SMF table
            // both DVO detections and smf detections should increment. however, there may be gaps in smf and 'invalid' rows in dvo
            numOfDuplicates = 0;
            numInvalidFlux = 0;
            for (long s = 0; s<nDet; s++) {

                isDuplicate = false;

                // count jumps in smf file
                if (s>0 && (ippIDet[s] != (ippIDet[s-1] + 1))) smfJumps++;

                while (ippIDet[s] > dvoDetections[d].meas.detID && dvoDetections[d].meas.detID <= maxDvoDetId) {

                    d++;
                    if (!dvoDetections[d].valid) invalidDvoRows++;
                }

                if (ippIDet[s] == dvoDetections[d].meas.detID ) {

                    // TODO HACK duplicates bug - don't include any duplicates
                    thisObjId = dvoDetections[d].ave.extID;
                    for (i=0; i<totalDetections; i++) {
                        if (thisObjId == objID[i]) {
                            removeList[numOfDuplicates+numInvalidFlux] = s+1;
                            numOfDuplicates++;
                            isDuplicate=true;
                            break;
                        }
                    }
                    // populate some arrays to be shoved in table in bulk later
                    objID[s] = thisObjId;
                    detectID[s] = dvoDetections[d].meas.extID;
                    ippObjID[s] = (uint64_t)dvoDetections[d].ave.catID*1000000000 + (uint64_t)dvoDetections[d].ave.objID;
                    ippDetectID[s] = dvoDetections[d].meas.detID;
                    imageID[s] = pspsImageId;
                    obsTimes[s] = obsTime;

                    ippToPsps_getFlux(zeroPoint, exposureTime, instMagErr[s], &instFluxErr[s]);

                    // check for invalid flux values (if not already labelled as a duplicate)
                    if ((!ippToPsps_getFlux(zeroPoint, exposureTime, peakMag[s], &peakFlux[s]) ||
                                !ippToPsps_getFlux(zeroPoint, exposureTime, instMag[s], &instFlux[s])) && !isDuplicate) {
                        removeList[numOfDuplicates+numInvalidFlux] = s+1;
                        numInvalidFlux++;
                    }

                    // store max/min objID
                    if (objID[s] > maxObjID) maxObjID = objID[s];
                    if (objID[s] < minObjID) minObjID = objID[s];

                    totalDetections++;
                }
                // if we get here then we cannot find this smf detection in DVO 
                else {

                    unmatched++;
                    continue;
                }

                if (dvoDetections[d].meas.detID == maxDvoDetId) break;
            }

            // write number of rows (detections) to ImageMeta
            numDetectionsOut = totalDetections - numOfDuplicates - numInvalidFlux;
            fits_write_col(this->fitsOut, TLONG, IMAGEMETA_NDETECT, 1, 1, 1, &numDetectionsOut, &status);

            // detections
            ippToPspsConfig_writeTable(this->config, fitsIn, this->fitsOut, nDet, "Detection", false);
            fits_write_col(this->fitsOut, TLONGLONG, DETECTION_OBJID, 1, 1, nDet, objID, &status);
            fits_write_col(this->fitsOut, TLONGLONG, DETECTION_DETECTID, 1, 1, nDet, detectID, &status);
            fits_write_col(this->fitsOut, TLONGLONG, DETECTION_IPPOBJID, 1, 1, nDet, ippObjID, &status);
            fits_write_col(this->fitsOut, TLONGLONG, DETECTION_IPPDETECTID, 1, 1, nDet, ippDetectID, &status);
            fits_write_col(this->fitsOut, TBYTE, DETECTION_FILTERID, 1, 1, nDet, filterIDs, &status);
            fits_write_col(this->fitsOut, TFLOAT, DETECTION_INSTFLUX, 1, 1, nDet, instFlux, &status);
            fits_write_col(this->fitsOut, TFLOAT, DETECTION_INSTFLUXERR, 1, 1, nDet, instFluxErr, &status);
            fits_write_col(this->fitsOut, TFLOAT, DETECTION_PEAKADU, 1, 1, nDet, peakFlux, &status);
            if (numOfDuplicates||numInvalidFlux) fits_delete_rowlist(this->fitsOut, removeList, numOfDuplicates+numInvalidFlux, &status);


            psLogMsg("ippToPsps", PS_LOG_INFO,
                    "+-----------+---------+----------+------------------+---------------+--------------+\n"
                    "| Extension | Rows in | Rows out | Missing from DVO | Duplicate IDs | Invalid Flux |\n"
                    "|  %s |  %5ld  |   %5ld  |      %5d       |    %5d      |    %5d     |\n",
                    extensionName, nDet, numDetectionsOut, unmatched, numOfDuplicates, numInvalidFlux);

            if (dvoDetections!= NULL) dvoFree(dvoDetections);
            SkyListFree(skylist);
        }
    }

    for (uint32_t i=0; i<MAXDETECT;i++) free(assocDate[i]);
    free(assocDate);

    // write number of images we have found into FrameMeta table
    status=0; 
    if (fits_movnam_hdu(this->fitsOut, BINARY_TBL, "FrameMeta", 0, &status))
        psError(PS_ERR_IO, false, "Can't move back to FrameMeta extension to write number of OTAs\n");
    else 
        fits_write_col(this->fitsOut, TSHORT, FRAMEMETA_NOTA, 1, 1, 1, &nOta, &status);

    status=0;
    if (fits_close_file(fitsIn, &status)) fits_report_error(stderr, status);

    // write results
    if (this->resultsFile) {

        if (fprintf(this->resultsFile, "%ld\n", minObjID) < 1 ) 
            psError(PS_ERR_IO, false, "Unable to write min Object ID to'%s'\n", this->resultsPath);
        if (fprintf(this->resultsFile, "%ld\n", maxObjID) < 1 ) 
            psError(PS_ERR_IO, false, "Unable to write max Object ID to'%s'\n", this->resultsPath);
    }

    psLogMsg("ippToPsps", PS_LOG_INFO, "Data written for a total of %d chips/OTAs", nOta);

    return PS_EXIT_SUCCESS;
}


