/** @file  pmSourceIO.c
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-18 02:44:19 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

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

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

#include "pmConfig.h"
#include "pmDetrendDB.h"

#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPALevel.h"
#include "pmFPAview.h"
#include "pmFPAfile.h"

#include "pmTrend2D.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmSpan.h"
#include "pmFootprintSpans.h"
#include "pmFootprint.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmModelFuncs.h"
#include "pmModel.h"
#include "pmModelUtils.h"
#include "pmModelClass.h"
#include "pmSourceMasks.h"
#include "pmSourceExtendedPars.h"
#include "pmSourceDiffStats.h"
#include "pmSourceSatstar.h"
#include "pmSourceLensing.h"
#include "pmSource.h"
#include "pmSourceFitModel.h"
#include "pmPSF.h"
#include "pmPSFtry.h"

#include "pmSourceIO.h"
#include "pmSourceOutputs.h"

// read in a readout from the fits file
psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header)
{
    PS_ASSERT_PTR_NON_NULL(fits, false);
    PS_ASSERT_PTR_NON_NULL(header, false);

    bool status;
    psEllipseAxes axes;

    // define PSF model type
    // XXX need to carry the extra model parameters
    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");

    // Read lookup table for model classes (if defined)
    pmModelClassReadHeader(header);

    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL");
    if (PSF_NAME != NULL) {
        modelType = pmModelClassGetType (PSF_NAME);
    }
    PS_ASSERT_INT_NONNEGATIVE(modelType, NULL);

    // We get the size of the table, and allocate the array of sources first because the table
    // is large and ephemeral --- when the table gets blown away, whatever is allocated after
    // the table is read blocks the free.  In fact, it's better to read the table row by row.
    long numSources = psFitsTableSize(fits); // Number of sources in table
    psArray *sources = psArrayAlloc(numSources); // Array of sources, to return

    // convert the table to the pmSource entriesa
    for (int i = 0; i < numSources; i++) {
        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
        if (!row) {
            psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i);
            psFree(sources);
            return NULL;
        }

	// here are the things we read from the FITS table (XXX modify names if needed)

	unsigned int ID  = psMetadataLookupU32 (&status, row, "ID");
	float X          = psMetadataLookupF32 (&status, row, "X");
        assert(status);  // poorly constructed cff
	float Y          = psMetadataLookupF32 (&status, row, "Y");
        assert(status);  // poorly constructed cff
        float flux       = psMetadataLookupF32 (&status, row, "FLUX"); // guess at the instrumental flux
        if (!isfinite(X) || !isfinite(Y) || !isfinite(flux)) {
            psError(PS_ERR_IO, true, "source ID %d is invalid x: %f y: %f flux: %f", ID, X, Y, flux);
            psFree(sources);
            return NULL;
        }


        float apRadius   = psMetadataLookupF32 (&status, row, "AP_RADIUS");
        float kronRadius = psMetadataLookupF32 (&status, row, "KRON_RADIUS");
        float petRadius  = psMetadataLookupF32 (&status, row, "PETRO_RADIUS");
        bool fitGalaxy   = psMetadataLookupU8 (&status, row, "FIT_GALAXY");
        bool psfStar     = psMetadataLookupU8 (&status, row, "PSF_STAR");

        float Rmajor     = psMetadataLookupF32 (&status, row, "R_MAJOR");
        float Rminor     = psMetadataLookupF32 (&status, row, "R_MINOR");
        float theta      = psMetadataLookupF32 (&status, row, "THETA");

        int   galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE");
        if (status) {
            galaxyModelType = pmModelClassGetLocalType(galaxyModelType);
        } else {
            galaxyModelType = -1;
        }
        float Sindex     = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index

        pmSource *source = pmSourceAlloc ();
        pmModel *model = pmModelAlloc (modelType);
        source->modelPSF  = model;
//        RoughClass wants source type to be unknown
//        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
        source->type = PM_SOURCE_TYPE_UNKNOWN;

	// XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and
        // stars that are masked
        if (psfStar) {
	    source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
	}

	// NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
	psF32 *PAR = model->params->data.F32;
	psF32 *dPAR = model->dparams->data.F32;

        source->seq       = ID;
        PAR[PM_PAR_XPOS]  = X;
        PAR[PM_PAR_YPOS]  = Y;

	dPAR[PM_PAR_XPOS] = 0.0;
	dPAR[PM_PAR_YPOS] = 0.0;

	PAR[PM_PAR_SKY]   = 0.0;
	dPAR[PM_PAR_SKY]  = 0.0;

	PAR[PM_PAR_I0]    = 1.0;
	dPAR[PM_PAR_I0]   = 0.0;

	source->sky       = PAR[PM_PAR_SKY];
	source->skyErr    = dPAR[PM_PAR_SKY];

	source->psfMag    = 0.0;
	source->psfMagErr = 0.0;
	source->apMag     = 0.0;
        source->apRadius  = apRadius;

	// we generate a somewhat fake PSF model here -- 
	// in most (all?) contexts, we will replace this with a measured psf model
	// elsewhere
	axes.major        = 1.0;
	axes.minor        = 1.0;
	axes.theta        = 0.0;
	pmPSF_AxesToModel (PAR, axes, modelType);

	// peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second
        source->peak      = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE);
        source->peak->xf  = X; // pmPeakAlloc converts X,Y to int, so reset here
        source->peak->yf  = Y; // pmPeakAlloc converts X,Y to int, so reset here
        source->peak->dx  = 0.0;
        source->peak->dy  = 0.0;

        source->moments = pmMomentsAlloc ();
	source->moments->Mx = X;
	source->moments->My = Y;
	source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment

        // Don't mark the moments as measured because that causes many fields to be left blank.
        // The moments code knows not to change the position or the Mrf for external sources
        // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;

	if (isfinite(petRadius)) {
	    source->extpars = pmSourceExtendedParsAlloc ();
	    source->extpars->petrosianRadius = petRadius;
	}

        if (fitGalaxy && galaxyModelType >= 0) {
            source->modelFits = psArrayAllocEmpty (1);
	    pmModel *model = pmModelAlloc(galaxyModelType);
	    psF32 *xPAR = model->params->data.F32;

	    xPAR[PM_PAR_SKY]  = 0.0;
	    xPAR[PM_PAR_I0]   = 1.0;
	    xPAR[PM_PAR_XPOS] = X;
	    xPAR[PM_PAR_YPOS] = Y;
	    
	    psEllipseAxes galaxyAxes;
	    galaxyAxes.major = Rmajor;
	    galaxyAxes.minor = Rminor;
	    galaxyAxes.theta = theta * PS_RAD_DEG;

	    pmPSF_AxesToModel (xPAR, galaxyAxes, galaxyModelType);
	    if (model->params->n > 7) {
                xPAR[PM_PAR_7] = 0.5 / Sindex;
	    }

	    psArrayAdd (source->modelFits, 1, model);

#ifdef notyet
            // XXX: set source->modelEXT to this model and flag as extended so that we have a better
            // shot at subtracting extended sources?

            // This doesn't work right. Need to do some more work to flesh out the model before it can be
            // used in subtaction. Better idea might be to do 2 passes in psphotFullForceReadout to
            // First find the best extended models and then do everything else.
            source->modelEXT = psMemIncrRefCounter(model);
            // is this safe? (no see above)
            source->type = PM_SOURCE_TYPE_EXTENDED;
            source->mode |= PM_SOURCE_MODE_EXTMODEL; 
#endif

	    psFree (model);
        }

        sources->data[i] = source;
        psFree(row);
    }

    return sources;
}

bool pmSourcesWrite_CFF(pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe) {
    
    char *extname = "SkyChip.cff";
    
    bool mdok;
    psF32 exptime = psMetadataLookupF32(&mdok, readout->parent->concepts, "CELL.EXPOSURE");
    PS_ASSERT(mdok, false);

    // write the definition of the model class type values to the header
    psMetadata *outputHeader = psMetadataAlloc();
    pmModelClassWriteHeader(outputHeader);

    psArray *table = psArrayAllocEmpty(sources->n);

    pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC");

    psString modelToChoose = psMetadataLookupStr(&mdok, recipe, "EXT_MODEL_TYPE_FOR_CFF");
    pmModelType selectedModelType = -1;
    if (mdok && modelToChoose != NULL) {
        if (strcmp(modelToChoose, "BEST")) {
            selectedModelType = pmModelClassGetType(modelToChoose);
        }
    }

    for (int i = 0; i < sources->n; i++) {
        pmSource *thisSource = sources->data[i];
        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;

        psF32 xPos, yPos, flux, rMajor, rMinor, theta;
        psS32 modelType = 0;
        bool fitGalaxy = false;
        bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false;
        psF32 sersicIndex = 0;
        // For now only perform galaxy fits on extended objects
        if (source->modelEXT == NULL) {
            pmModel *model = source->modelPSF;
            if (model == NULL) continue;
            psF32 *PAR = model->params->data.F32;
            if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY]) ||
                !isfinite(source->psfFlux)) {
                continue;
            }

            xPos = model->params->data.F32[PM_PAR_XPOS];
            yPos = model->params->data.F32[PM_PAR_YPOS];
            flux = source->psfFlux;
            rMajor = 0;
            rMinor = 0;
            theta = 0;
        } else {
            //   Find the model with the selected type. If selected type is -1 choose the one selected ad
            //   modelEXT which was the best
            int iModel = -1;
            if (source->modelEXT) {
                pmModelType ext_model_type =  selectedModelType != -1  ? selectedModelType : source->modelEXT->type;
                for (int j=0; j<source->modelFits->n; j++) {
                    pmModel *aModel = source->modelFits->data[j];
                    if (aModel->type == ext_model_type) {
                        iModel = j;
                        break;
                    }
                }
            }
            if (iModel == -1) {
                // Can this happen? perhaps if model type is qgauss (extended but below s/n for model fits?
                // XXX: Should this be an assert?
                continue;
            }
            pmModel *model = source->modelFits->data[iModel];
            psF32 *PAR = model->params->data.F32;
            xPos = PAR[PM_PAR_XPOS];
            yPos = PAR[PM_PAR_YPOS];
            if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
                // XXX: do we need to handle this type
                continue;
	    } else {
		if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY]) ||
                    !isfinite(model->mag)) { 
                    // bad model
                    continue;
		}
                psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
                rMajor = axes.major;
                rMinor = axes.minor;
                theta  = axes.theta*PS_DEG_RAD;
                flux = pow(10.0, -0.4*model->mag);
                fitGalaxy = true;
            }
            modelType = model->type;
            if (modelType == sersicModelType) {
                PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false);
                sersicIndex = 0.5 / PAR[PM_PAR_7];
            }
        }
        psMetadata *row = psMetadataAlloc();
        psMetadataAddU32 (row, PS_LIST_TAIL, "ID",         0,   "IPP detection identifier",  source->seq);
        psMetadataAddF32 (row, PS_LIST_TAIL, "X",                0, "x coordinate",          xPos);
        psMetadataAddF32 (row, PS_LIST_TAIL, "Y",                0, "y coordinate",          yPos);
        psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX",             0, "flux per second",       flux/exptime);
        psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS",        0, "aperture radius",       source->apRadius);
        psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS",      0, "Kron radius",           source->moments->Mrf * 2.5);
        psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS",     0, "Petrosian Radius",      source->extpars ? source->extpars->petrosianRadius : NAN);
        psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY",      0, "source has xfit",       fitGalaxy); 
        psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR",        0, "source was psf star",   psfStar);
        psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR",          0, "radius of major axis",  rMajor);
        psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR",          0, "radius of minor axis",  rMinor);
        psMetadataAddF32 (row, PS_LIST_TAIL, "THETA",            0, "theta",                 theta);
        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",            modelType);
        psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex);

        psArrayAdd(table, 100, row);
        psFree(row);
    }

    if (!psFitsWriteTable(fits, outputHeader, table, extname)) {
        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
        psFree(table);
        psFree(outputHeader);
        return false;
    }
    psFree(table);
    psFree(outputHeader);

    return true;
}
