/* $Id: naco_spc.c,v 1.42 2009-05-12 12:43:43 llundin Exp $
 *
 * This file is part of the NACO Pipeline
 * Copyright (C) 2002,2003 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */

/*
 * $Author: llundin $
 * $Date: 2009-05-12 12:43:43 $
 * $Revision: 1.42 $
 * $Name: not supported by cvs2svn $
 */

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

/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/

#include "naco_spc.h"
#include "naco_pfits.h"

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

#ifdef HAVE_GSL
#include <gsl/gsl_multimin.h>
#endif


/*-----------------------------------------------------------------------------
                               Private types
 -----------------------------------------------------------------------------*/
typedef struct {
    double x;
    int    index;
} vector_index;


/*-----------------------------------------------------------------------------
                                   Private functions
 -----------------------------------------------------------------------------*/

static int naco_vector_get_maxpos_window(const cpl_vector *, int, int);

#if 0
static int vector_index_compare(const void *, const void *);
#endif

static cpl_error_code naco_framelist_get_state(const irplib_framelist *,
                                                    int, cpl_boolean,
                                                    cpl_boolean *,
                                                    const cpl_propertylist *);

/*----------------------------------------------------------------------------*/
/**
 * @defgroup naco_spc     Spectro functionality
 */
/*----------------------------------------------------------------------------*/

/**@{*/

/*----------------------------------------------------------------------------*/
/**
  @brief    Use 1D cross-correlation to recompute the offsets in one direction
  @param    offset   The recomputed offsets, XC-changes are rounded off to int
  @param    goffset  The 1st guess of the offsets
  @param    do_wave  CPL_TRUE to recompute Y-offsets, else X-offsets
  @param    self     The imagelist to correlate
  @return   CPL_ERROR_NONE on success, otherwise the relevant CPL error code
  @note     The Y-direction is assumed to be the wavelength direction

  The first element(s) are not updated, the remainder are relative to the 1st.

  The remainder have their signs inverted, so they work with 
    cpl_geom_img_offset_saa()...

 */
/*----------------------------------------------------------------------------*/
cpl_error_code naco_vector_correlate_imagelist_1d(cpl_vector * offset,
                                                  const cpl_vector * goffset,
                                                  cpl_boolean do_wave,
                                                  const cpl_imagelist * self)
{

    cpl_image  * d2d      = NULL;
    cpl_image  * dspat1d  = NULL;
    cpl_vector * vspat1dp = NULL;
    cpl_vector * vspat1di = NULL;
    cpl_vector * vxc      = NULL;
    const int    nsize    = cpl_imagelist_get_size(self);
    const int    maxerr   = 20;
    int          gleni, glenprev = 0; /* Avoid (false) uninit warning */

    /* Wavelenth/Spatial resolution */
    const int    nz = (do_wave ? cpl_image_get_size_y : cpl_image_get_size_x)
        (cpl_imagelist_get_const(self, 0));
    const char   dirch = do_wave ? 'Y' : 'X';
    const int    direc = do_wave ? 1 : 0;
    int          ii, i, iprev = 0; /* Avoid (false) uninit warning */
    double       xcmin    = 1.0;


    bug_if(self     == NULL);

    bug_if(offset  == NULL);
    bug_if(goffset == NULL);
    bug_if(cpl_vector_get_size(offset)  != nsize);
    bug_if(cpl_vector_get_size(goffset) != nsize);

    /* Process frames in order of offset */
    for (ii = 0; ii < nsize; ii++) {
        i = ii;

        cpl_msg_info(cpl_func, "%c-offset(%d:%d): %g", dirch, ii, i,
                     cpl_vector_get(goffset, i));

    }

    vxc = cpl_vector_new(maxerr);

    /* Process frames in order of offset */
    for (ii = 0; ii < nsize; ii++) {

        int maxoff, halfoff;
        int ixc, ixc0;
        int irealfirst, ireallast;
        int ioff;
        double xci;
        double doff, dnewoff;

        i = ii;

        gleni = (int)round(cpl_vector_get(goffset, i));

        /* Need double for the cpl_vector */
        d2d = cpl_image_cast(cpl_imagelist_get_const(self, i), CPL_TYPE_DOUBLE);

        /* Sum over the spatial dimension */
        dspat1d = cpl_image_collapse_create(d2d, direc);
        /* A skip before vspat1di has wrapped will leak the pixel buffer */
        cpl_image_delete(d2d);
        d2d = NULL;

        cpl_vector_delete(vspat1di);
        vspat1di = cpl_vector_wrap(nz, cpl_image_get_data_double(dspat1d));
        (void)cpl_image_unwrap(dspat1d);
        dspat1d = NULL;
        bug_if(0);

#if 0
        cpl_plot_vector("set grid;", "t '1D-Profile i' w linespoints", "",
                        vspat1di);
#endif
        if (ii == 0) {
            vspat1dp = vspat1di;
            vspat1di = NULL;
            glenprev = gleni;
            iprev = i;
            continue;
        }

        halfoff = maxerr + abs(gleni - glenprev);
        if (halfoff > nz-1) halfoff = nz - 1;
        maxoff = 2 * halfoff + 1;

        bug_if(cpl_vector_set_size(vxc, maxoff));

        ixc0 = cpl_vector_correlate(vxc, vspat1di, vspat1dp);
        bug_if(0);

        /* Limit the search range to 1+2*maxerr to reduce risk of false max */
        irealfirst = halfoff - (glenprev - gleni) - maxerr;
        ireallast  = halfoff - (glenprev - gleni) + maxerr;

        ixc = naco_vector_get_maxpos_window(vxc, irealfirst, ireallast);
        bug_if(0);

        if (ixc != ixc0) {
            cpl_msg_warning(cpl_func, "%c-False maximum(%d:%d): %d <+. %d. "
                            "0 <= %d => %d < %d", dirch, ii, i, ixc0, ixc,
                            irealfirst, ireallast, maxoff);
        }

        ioff = ixc - halfoff;
        doff = cpl_vector_get(goffset, i) - cpl_vector_get(goffset, iprev);

        xci = cpl_vector_get(vxc, ixc);

        if (xci < xcmin) xcmin = xci;

        if (ioff != (int)round(doff)) {
            cpl_msg_warning(cpl_func, "%c-XC(%d)=%g changes offset: %g - %g = "
                            "%g => %d = %d - %d", dirch, i, xci,
                            cpl_vector_get(goffset, i),
                            cpl_vector_get(goffset, iprev), doff, ioff,
                            gleni, glenprev);
            dnewoff = (double)-ioff;
        } else {
            cpl_msg_info(cpl_func, "%c-XC(%d)=%g confirms offset: %g - %g = %g "
                         "<=> %d = %d - %d", dirch, i, xci,
                         cpl_vector_get(goffset, i),
                         cpl_vector_get(goffset, iprev), doff, ioff,
                         gleni, glenprev);
            dnewoff = -doff;
        }
        bug_if(0);
        bug_if(cpl_vector_set(offset, i, dnewoff));

#if 0
        cpl_plot_vector("set grid;", "t 'Cross-correlation' w linespoints",
                        "", vxc);

#endif
    }

    cpl_msg_info(cpl_func, "Minimum 1D-spatial XC for %d sets: %g",
                 nsize, xcmin);

    end_skip;

    cpl_image_delete(d2d);
    (void)cpl_image_unwrap(dspat1d);
    cpl_vector_delete(vspat1dp);
    cpl_vector_delete(vspat1di);
    cpl_vector_delete(vxc);

    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Readd all the pairs, preserving the number of images
  @param    self      The imagelist to readd
  @return   CPL_ERROR_NONE on success, otherwise the relevant CPL error code
  @see naco_imagelist_split()

 */
/*----------------------------------------------------------------------------*/
cpl_error_code naco_imagelist_add_split(cpl_imagelist * self)
{

    cpl_image     * copy  = NULL;
    int             n     = cpl_imagelist_get_size(self);
    int             i;

    bug_if(self == NULL);

    bug_if(n&1);


    for (i=0; i < n; i += 2) {
        cpl_image * first  = cpl_imagelist_get(self, i);
        cpl_image * second = cpl_imagelist_get(self, i+1);

        bug_if(cpl_image_subtract(first, second));

        copy = cpl_image_multiply_scalar_create(first, -1.0);

        bug_if(cpl_imagelist_set(self, copy, i+1));
        copy = NULL;

    }

    end_skip;

    cpl_image_delete(copy);

    return cpl_error_get_code();

}



/*----------------------------------------------------------------------------*/
/**
  @brief    Fo each image append also its inverted
  @param    self      The imagelist to extend
  @return   CPL_ERROR_NONE on success, otherwise the relevant CPL error code
  @note The list is modified to D1, D2, ..., Dn, -D1, -D2, ..., -Dn.

 */
/*----------------------------------------------------------------------------*/
cpl_error_code naco_imagelist_append_invert(cpl_imagelist * self)
{

    cpl_image     * BA = NULL;
    const int       n      = cpl_imagelist_get_size(self);
    int             i;


    bug_if(self == NULL);

    for (i=0; i < n; i++) {
        const cpl_image * AB  = cpl_imagelist_get(self, i);

        BA = cpl_image_multiply_scalar_create(AB, -1.0);

        bug_if(cpl_imagelist_set(self, BA, n+i));
        BA = NULL;

    }

    end_skip;

    cpl_image_delete(BA);

    return cpl_error_get_code();

}


/*----------------------------------------------------------------------------*/
/**
  @brief    Split the images in an imagelist into positive and negative images
  @param    self      The imagelist to split
  @return   CPL_ERROR_NONE on success, otherwise the relevant CPL error code
  @note Each input image is replaced by two, a positive and a negative, the
  length of the list thus doubles. The positive image is always placed before
  the negative.

 */
/*----------------------------------------------------------------------------*/
cpl_error_code naco_imagelist_split(cpl_imagelist * self)
{

    cpl_imagelist * copy   = cpl_imagelist_new();
    cpl_image     * image  = NULL;
    cpl_image     * im_neg = NULL;
    int             j      = cpl_imagelist_get_size(self);
    int             i      = 0;


    skip_if(self == NULL);

    /* Move all images to the copy - and revert their order */
    j--;
    for (; j >= 0; j--, i++) {
        bug_if(cpl_imagelist_set(copy, cpl_imagelist_unset(self, j), i));
    }

    j++;
    i--;

    bug_if(j != 0);

    /* Move all images back to self, revert their order again while splitting */
    for (; i >= 0; i--, j += 2) {
        image = cpl_imagelist_unset(copy, i);

        /* The negative flux */
        im_neg = cpl_image_duplicate(image);
        bug_if(cpl_image_threshold(im_neg, -FLT_MAX, 0.0, 0.0, 0.0));
        bug_if(cpl_image_multiply_scalar(im_neg, -1.0));


        /* The positive flux */
        bug_if(cpl_image_threshold(image, 0.0, FLT_MAX, 0.0, 0.0));

        bug_if(cpl_imagelist_set(self, image, j));
        image = NULL;

        bug_if(cpl_imagelist_set(self, im_neg, j+1));
        im_neg = NULL;

    }

    end_skip;

    cpl_image_delete(image);
    cpl_image_delete(im_neg);
    cpl_imagelist_delete(copy);

    return CPL_ERROR_NONE;

}

/*----------------------------------------------------------------------------*/
/**
  @brief    Create a string suitable for frame comparison in spectroscopy
  @param    self   Frame to create the new tag from
  @param    plist  The propertylist of the frame
  @param    dummy  A non-negative number (required in the API, but unused here)
  @return   Comparison string or NULL on error
  @note The comparison string must be deallocated with cpl_free().

  The tag currently concatenates the SPECMODE, OPTI1.ID, CWLEN and DIT.

 */
/*----------------------------------------------------------------------------*/
char * naco_spc_make_tag(const cpl_frame* self, const cpl_propertylist * plist,
			 int dummy)
{

    char       * tag = NULL;
    double       wlen, dit;
    const char * specmode; 
    const char * slitname; 


    bug_if (0);

    bug_if(self  == NULL);
    bug_if(plist == NULL);
    bug_if(dummy < 0); /* Avoid warning of unused variable */

    specmode = irplib_pfits_get_string(plist, NACO_PFITS_STRING_SPECMODE);
    skip_if(cpl_error_get_code());

    slitname = irplib_pfits_get_string(plist, NACO_PFITS_STRING_SLITNAME);
    skip_if(cpl_error_get_code());

    dit  = irplib_pfits_get_double(plist, NACO_PFITS_DOUBLE_DIT);
    skip_if(cpl_error_get_code());

    wlen = irplib_pfits_get_double(plist, NACO_PFITS_DOUBLE_CWLEN);
    skip_if(cpl_error_get_code());

    tag = cpl_sprintf("%s:%s:%.5f:%.5f", specmode, slitname, wlen, dit);
    bug_if(tag == NULL);

    end_skip;

    if (cpl_error_get_code()) {
        cpl_free(tag);
        tag = NULL;
    }

    return tag;

}


/*----------------------------------------------------------------------------*/
/**
  @brief    Fill the list of difference images from on/off frames
  @param    self      The imagelist to fill
  @param    onofflist The list of frames with the on/off-frames
  @param    onoffkeys The on/off keys, positive: on; 0: off; neg.: don't care
  @return   CPL_ERROR_NONE on success, otherwise the relevant CPL error code

  The integer properties in onoffkeys specify which keys in the frames to check
  for the on/off status.
  In the propertylists of the frames
  keys of type boolean are assumed to indicate an on-state iff true,
  keys of numerical type are assumed to indicate an on-state iff positive,
  keys of other types (string) are not supported.

  For an off-frame all on/offs specified in onoffkeys must be off.

  For an on-frame at least one on/off specified in onoffkeys must be on.

  Additionally for the on-frames the requirement of each individual on/off
  must be met. These requirements are defined by the values of the integer
  type properties in onoffkeys:
  1) A positive property value means that the state must be on,
  2) A property value of 0 means that the state must be off,
  3) A negative property value means that state may be either on or off.
  This implies that onoffkeys must have at least one property with a
  non-zero value.

  All frames must have identical NAXIS* and DIT.

 */
/*----------------------------------------------------------------------------*/
cpl_error_code naco_imagelist_load_diff(cpl_imagelist * self,
                                        const irplib_framelist * onofflist,
                                        const cpl_propertylist * onoffkeys)
{
    cpl_image     * img_off = NULL;
    cpl_image     * diff = NULL;
    int             nfiles, ndiff;
    int             ion, ioff;

    skip_if (0);
    skip_if (self      == NULL);
    skip_if (onofflist == NULL);

    skip_if(cpl_imagelist_get_size(self) != 0);

    skip_if (irplib_framelist_contains(onofflist, "NAXIS1", CPL_TYPE_INT,
                                       CPL_TRUE, 0.0));
    skip_if (irplib_framelist_contains(onofflist, "NAXIS2", CPL_TYPE_INT,
                                       CPL_TRUE, 0.0));

    skip_if (irplib_framelist_contains(onofflist, NACO_PFITS_DOUBLE_DIT,
                                       CPL_TYPE_DOUBLE, CPL_TRUE, 1e-5));
   
    nfiles = irplib_framelist_get_size(onofflist);

    for (ioff = 0, ion = 0, ndiff = 0; ioff < nfiles && ion < nfiles;
         ioff++, ion++) {
        /* Each completed iteration produces one on-off image */
        cpl_boolean need_on = CPL_TRUE;

        const cpl_frame * frame;
        const char * name;


        /* Find next off-frame - and perhaps also the next on-frame */
        for (; ioff < nfiles; ioff++) {
            cpl_boolean is_off;

            skip_if(naco_framelist_get_state(onofflist, ioff, CPL_FALSE,
                                                  &is_off, onoffkeys));

            if (is_off) {
                break;
            } else if (need_on && ioff >= ion) {
                ion = ioff;
                need_on = CPL_FALSE;
            }
        }

        bug_if(0);

        if (ioff == nfiles) break;

        if (need_on) {
            /* Find next on-frame */
            if (ion == ioff) ion++; /* No need to check an off-frame */
            for (; ion < nfiles; ion++) {

                cpl_boolean is_on;

                skip_if(naco_framelist_get_state(onofflist, ion, CPL_TRUE,
                                                      &is_on, onoffkeys));

                if (is_on) break;

            }

            bug_if(0);

            if (ion == nfiles) break;
        }

        frame = irplib_framelist_get_const(onofflist, ion);
        name = cpl_frame_get_filename(frame);

        bug_if(0);

        diff = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
        error_if(diff == NULL, cpl_error_get_code(),
                 "Could not load on-image from %s", name);

        frame = irplib_framelist_get_const(onofflist, ioff);
        name = cpl_frame_get_filename(frame);

        bug_if(0);

        img_off = cpl_image_load(name, CPL_TYPE_FLOAT, 0, 0);
        error_if(img_off == NULL, cpl_error_get_code(),
                 "Could not load off-image from %s", name);

        error_if(cpl_image_subtract(diff, img_off), cpl_error_get_code(),
                 "Could not subtract frame %d from %d", ioff, ion);

        cpl_image_delete(img_off);
        img_off = NULL;

        bug_if(cpl_imagelist_set(self, diff, ndiff));
        diff = NULL;

        ndiff++;
    }

    bug_if(cpl_imagelist_get_size(self) != ndiff);

    error_if(ndiff == 0 && ion == nfiles, CPL_ERROR_DATA_NOT_FOUND,
             "The %d frame(s) contain(s) no on-frames", nfiles);
    error_if(ndiff == 0 && ioff == nfiles, CPL_ERROR_DATA_NOT_FOUND,
             "The %d frame(s) contain(s) no off-frames", nfiles);

    if (2 * ndiff < nfiles) {
        cpl_msg_warning(cpl_func, "The %d frames contains only %d pairs of "
                        "on/off-frames, ignoring the rest", nfiles, ndiff);
    }

    bug_if(cpl_imagelist_get_size(self) < 1);

    end_skip;

    cpl_image_delete(img_off);
    cpl_image_delete(diff);

    return cpl_error_get_code();
}


/**@}*/



/*----------------------------------------------------------------------------*/
/**
  @internal
  @brief    Determine if a frame is an on- or off- frame
  @param    self      The list of frames with the on/off-frames
  @param    iframe    The frame to inquire about
  @param    is_on     Which state (on/off) to look for
  @param    pis_ok    The state is as requested by is_on
  @param    onoffkeys The on/off keys, positive: on; 0: off; neg.: don't care
  @return   CPL_TRUE iff the state is as required (undefined on error)
  @see naco_imagelist_load_diff()
  @note The values of onoffkeys are only used with is_on is true.

 */
/*----------------------------------------------------------------------------*/
static
cpl_error_code naco_framelist_get_state(const irplib_framelist * self,
                                        int iframe,
                                        cpl_boolean is_on,
                                        cpl_boolean * pis_ok,
                                        const cpl_propertylist * onoffkeys)
{

    cpl_boolean state = CPL_FALSE;
    const cpl_propertylist * plist
        = irplib_framelist_get_propertylist_const(self, iframe);
    const int nonoffkeys = cpl_propertylist_get_size(onoffkeys);
    int i;
    cpl_boolean onoffkeys_has_one_on = CPL_FALSE;

    bug_if (0);
    bug_if (self     == NULL);
    bug_if (iframe   <  0);
    bug_if (iframe   >= irplib_framelist_get_size(self));
    bug_if (pis_ok   == NULL);
    bug_if (onoffkeys == NULL);
    bug_if (plist    == NULL);

    for (i=0; i < nonoffkeys; i++) {
        const cpl_property * prop = cpl_propertylist_get_const(onoffkeys, i);
        const char         * name = cpl_property_get_name(prop);
        const int            ireq = cpl_property_get_int(prop);
        const cpl_type  statetype = cpl_propertylist_get_type(plist, name);
        cpl_boolean state_is_on;

        error_if(0, cpl_error_get_code(), "On/off property number %d, %s",
                 1+i, name ? name : "<NULL>");

        switch (statetype) {
        case CPL_TYPE_CHAR:
            state_is_on = cpl_propertylist_get_char(plist, name) > 0
                ? CPL_TRUE : CPL_FALSE;
            break;
        case CPL_TYPE_BOOL:
            state_is_on = cpl_propertylist_get_bool(plist, name);
            break;
        case CPL_TYPE_INT:
            state_is_on = cpl_propertylist_get_int(plist, name) > 0
                ? CPL_TRUE : CPL_FALSE;
            break;
        case CPL_TYPE_LONG:
            state_is_on = cpl_propertylist_get_long(plist, name) > 0
                ? CPL_TRUE : CPL_FALSE;
            break;
        case CPL_TYPE_FLOAT:
            state_is_on = cpl_propertylist_get_float(plist, name) > 0.0
                ? CPL_TRUE : CPL_FALSE;
            break;
        case CPL_TYPE_DOUBLE:
            state_is_on = cpl_propertylist_get_double(plist, name) > 0.0
                ? CPL_TRUE : CPL_FALSE;
            break;
        default:
            cpl_ensure_code(0, CPL_ERROR_UNSUPPORTED_MODE);
        }

        if (is_on) {
            if (state_is_on) {
                if (ireq == 0) break; /* This state is on, but shouldn't be */
                onoffkeys_has_one_on = CPL_TRUE; /* A state is corrrectly on */
            } else if (ireq >  0) {
                break; /* This state is off, but shouldn't be */
            }
        } else if (state_is_on) {
            break; /* This state is on, but all states should be off */
        }
    }

    if (i == nonoffkeys) {
        if (is_on) {
            /* All states were on or off, as specified by onoffkeys */
            state = onoffkeys_has_one_on; /* At least one state must be on */
        } else {
            /* All states were off, as they should be */
            state = CPL_TRUE;
        }
    }

    *pis_ok = state;

    end_skip;

    return cpl_error_get_code();
}

#if 0
/*----------------------------------------------------------------------------*/
/**
   @internal
   @brief   Define sorting order of vector_index struct for ANSIs qsort()
   @param   a     1st operand, points to a vector_index struct
   @param   b     2nd operand, points to a vector_index struct
   @note Uses void pointers due to the API of ANSIs qsort().
   @see qsort()

*/
/*----------------------------------------------------------------------------*/
static int vector_index_compare(const void * a, const void * b)
{
    return
        ( (const vector_index*)a)->x < ((const vector_index*)b)->x ? -1 :
        (((const vector_index*)a)->x > ((const vector_index*)b)->x ?  1 : 0);  
}

#endif
/*----------------------------------------------------------------------------*/
/**
  @brief    Get the index of the maximum of the cpl_vector
  @param    v    const cpl_vector
  @param    i0   The first index to check, 0 for first
  @param    i1   The last index to check
  @return   the index of the maximum value of the vector or negative on error
  @see      cpl_vector_get_max()
 */
/*----------------------------------------------------------------------------*/
static int naco_vector_get_maxpos_window(const cpl_vector * v, int i0, int i1)
{
    int    i    = i0;
    int    imax = i0;
    double vmax = cpl_vector_get(v, imax);

    cpl_ensure(v != NULL, CPL_ERROR_NULL_INPUT,    -1);
    cpl_ensure(0  <= i0,  CPL_ERROR_ILLEGAL_INPUT, -2);
    cpl_ensure(i0 <= i1,  CPL_ERROR_ILLEGAL_INPUT, -3);
    cpl_ensure(i1 < cpl_vector_get_size(v), CPL_ERROR_ILLEGAL_INPUT, -4);

    for (i = i0; i <= i1; i++) if (cpl_vector_get(v, i) > vmax) {
        imax = i;
        vmax = cpl_vector_get(v, imax);
    }

    return imax;
}