/* $Id: montecarlo.c,v 1.5 2013-04-24 14:14:13 cgarcia Exp $
 *
 * This file is part of the FORS Data Reduction Pipeline
 * Copyright (C) 2002-2010 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  02110-1301  USA
 */

/*
 * $Author: cgarcia $
 * $Date: 2013-04-24 14:14:13 $
 * $Revision: 1.5 $
 * $Name: not supported by cvs2svn $
 */

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

#include <string.h>
#include <cpl.h>
#include <moses.h>
#include <fors_dfs.h>

static int montecarlo_create(cpl_plugin *);
static int montecarlo_exec(cpl_plugin *);
static int montecarlo_destroy(cpl_plugin *);
static int montecarlo(cpl_parameterlist *, cpl_frameset *);

static char montecarlo_description[] =
"This recipe is used to test the mos_montecarlo_polyfit() function.\n"
"It accepts a table with columns x, y, y_err, derives the best polynomial\n"
"fit y = p(x), and produces a table with the polynomial 1-sigma accuracy\n"
"on the given set of x coordinates.\n"
"Input files:\n\n"
"  DO category:               Type:       Explanation:         Required:\n"
"  TABLE                      Raw         Table to evaluate       Y\n\n"
"Output files:\n\n"
"  DO category:               Data type:  Explanation:\n"
"  MODEL_ERROR                FITS image  Model error at different x\n\n";

#define montecarlo_exit(message)               \
{                                             \
if ((const char *)message != NULL) cpl_msg_error(recipe, message);  \
cpl_table_delete(table);                  \
cpl_table_delete(points);                  \
cpl_table_delete(model_error);                  \
cpl_polynomial_delete(p);                    \
cpl_msg_indent_less();                        \
return -1;                                    \
}

#define montecarlo_exit_memcheck(message)       \
{                                               \
if ((const char *)message != NULL) cpl_msg_info(recipe, message);     \
cpl_table_delete(table);                    \
cpl_table_delete(points);                  \
cpl_table_delete(model_error);                  \
cpl_polynomial_delete(p);                    \
cpl_msg_indent_less();                          \
return 0;                                       \
}


/**
 * @brief    Build the list of available plugins, for this module. 
 *
 * @param    list    The plugin list
 *
 * @return   0 if everything is ok, -1 otherwise
 *
 * Create the recipe instance and make it available to the application 
 * using the interface. This function is exported.
 */

int cpl_plugin_get_info(cpl_pluginlist *list)
{
    cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe );
    cpl_plugin *plugin = &recipe->interface;

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    FORS_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "montecarlo",
                    "Test function mos_montecarlo_polyfit()",
                    montecarlo_description,
                    "Carlo Izzo",
                    PACKAGE_BUGREPORT,
    "This file is currently part of the FORS Instrument Pipeline\n"
    "Copyright (C) 2002-2010 European Southern Observatory\n\n"
    "This program is free software; you can redistribute it and/or modify\n"
    "it under the terms of the GNU General Public License as published by\n"
    "the Free Software Foundation; either version 2 of the License, or\n"
    "(at your option) any later version.\n\n"
    "This program is distributed in the hope that it will be useful,\n"
    "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
    "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
    "GNU General Public License for more details.\n\n"
    "You should have received a copy of the GNU General Public License\n"
    "along with this program; if not, write to the Free Software Foundation,\n"
    "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
                    montecarlo_create,
                    montecarlo_exec,
                    montecarlo_destroy);

    cpl_pluginlist_append(list, plugin);
    
    return 0;
}


/**
 * @brief    Setup the recipe options    
 *
 * @param    plugin  The plugin
 *
 * @return   0 if everything is ok
 *
 * Defining the command-line/configuration parameters for the recipe.
 */

static int montecarlo_create(cpl_plugin *plugin)
{
    cpl_recipe    *recipe;
    cpl_parameter *p;


    /* 
     * Check that the plugin is part of a valid recipe 
     */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    /* 
     * Create the parameters list in the cpl_recipe object 
     */

    recipe->parameters = cpl_parameterlist_new(); 


    /*
     * Name of input table columns
     */

    p = cpl_parameter_new_value("fors.montecarlo.x",
                                CPL_TYPE_STRING,
                                "Name of independent variable column",
                                "fors.montecarlo",
                                "x");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "x");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    p = cpl_parameter_new_value("fors.montecarlo.y",
                                CPL_TYPE_STRING,
                                "Name of dependent variable column",
                                "fors.montecarlo",
                                "y");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "y");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    p = cpl_parameter_new_value("fors.montecarlo.sigma",
                                CPL_TYPE_STRING,
                                "Name of error column on dependent variable",
                                "fors.montecarlo",
                                "");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sigma");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Order of fitting polynomial
     */

    p = cpl_parameter_new_value("fors.montecarlo.order",
                                CPL_TYPE_INT,
                                "Order of fitting polynomial",
                                "fors.montecarlo",
                                1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "order");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Origin of x for fit
     */

    p = cpl_parameter_new_value("fors.montecarlo.zero",
                                CPL_TYPE_DOUBLE,
                                "Origin of x for fit",
                                "fors.montecarlo",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "zero");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Start x for evaluation
     */

    p = cpl_parameter_new_value("fors.montecarlo.start",
                                CPL_TYPE_DOUBLE,
                                "Start x for evaluation",
                                "fors.montecarlo",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "start");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * End x for evaluation
     */

    p = cpl_parameter_new_value("fors.montecarlo.end",
                                CPL_TYPE_DOUBLE,
                                "End x for evaluation",
                                "fors.montecarlo",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "end");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Evaluation sampling
     */
    
    p = cpl_parameter_new_value("fors.montecarlo.step",
                                CPL_TYPE_DOUBLE,
                                "x sampling interval",
                                "fors.montecarlo",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "step");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Statistical sample
     */

    p = cpl_parameter_new_value("fors.montecarlo.trials",
                                CPL_TYPE_INT,
                                "Size of statistical sample",
                                "fors.montecarlo",
                                100);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "trials");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    return 0;
}


/**
 * @brief    Execute the plugin instance given by the interface
 *
 * @param    plugin  the plugin
 *
 * @return   0 if everything is ok
 */

static int montecarlo_exec(cpl_plugin *plugin)
{
    cpl_recipe *recipe;
    
    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    return montecarlo(recipe->parameters, recipe->frames);
}


/**
 * @brief    Destroy what has been created by the 'create' function
 *
 * @param    plugin  The plugin
 *
 * @return   0 if everything is ok
 */

static int montecarlo_destroy(cpl_plugin *plugin)
{
    cpl_recipe *recipe;
    
    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    cpl_parameterlist_delete(recipe->parameters); 

    return 0;
}


/**
 * @brief    Interpret the command line options and execute the data processing
 *
 * @param    parlist     The parameters list
 * @param    frameset    The set-of-frames
 *
 * @return   0 if everything is ok
 */

static int montecarlo(cpl_parameterlist *parlist, cpl_frameset *frameset)
{

    const char *recipe = "montecarlo";


    /*
     * Input parameters
     */

    const char *x;
    const char *y;
    const char *sigma;
    int         order;
    double      start;
    double      end;
    double      step;
    int         trials;
    double      zero;

    /*
     * CPL objects
     */

    cpl_table      *table = NULL;
    cpl_table      *model_error = NULL;
    cpl_table      *points = NULL;
    cpl_polynomial *p = NULL;

    /*
     * Auxiliary variables
     */

    char   *version = "0.1";
    char   *table_tag = "TABLE";
    char   *model_error_tag = "MODEL_ERROR";
    double *xdata;
    int     ntables;
    int     i, count;


    cpl_msg_set_indentation(2);


    /* 
     * Get configuration parameters
     */

    cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
    cpl_msg_indent_more();

    x = dfs_get_parameter_string(parlist, "fors.montecarlo.x", NULL);
    y = dfs_get_parameter_string(parlist, "fors.montecarlo.y", NULL);
    sigma = dfs_get_parameter_string(parlist, "fors.montecarlo.sigma", NULL);
    order = dfs_get_parameter_int(parlist, "fors.montecarlo.order", NULL);
    if (order < 0)
        montecarlo_exit("Invalid polynomial order");
    start = dfs_get_parameter_double(parlist, "fors.montecarlo.start", NULL);
    end = dfs_get_parameter_double(parlist, "fors.montecarlo.end", NULL);
    if (end <= start)
        montecarlo_exit("Invalid interval");
    step = dfs_get_parameter_double(parlist, "fors.montecarlo.step", NULL);
    if (step >= end - start || step <= 0.0)
        montecarlo_exit("Invalid step");
    trials = dfs_get_parameter_int(parlist, "fors.montecarlo.trials", NULL);
    if (trials < 2)
        montecarlo_exit("At least 2 trials should be performed");
    zero = dfs_get_parameter_double(parlist, "fors.montecarlo.zero", NULL);

    if (cpl_error_get_code())
        montecarlo_exit("Failure getting the configuration parameters");

    ntables = cpl_frameset_count_tags(frameset, table_tag);
    if (ntables == 0) {
        cpl_msg_error(recipe, "Missing required input: %s", table_tag);
        montecarlo_exit(NULL);
    }
    if (ntables > 1) {
        cpl_msg_error(recipe, "Too many in input: %s", table_tag);
        montecarlo_exit(NULL);
    }

    table = dfs_load_table(frameset, table_tag, 1);
    points = cpl_table_new(cpl_table_get_nrow(table));

    if (cpl_table_has_column(table, x)) {
        cpl_table_move_column(points, x, table);
        if (!cpl_table_has_column(points, "x")) {
            cpl_table_name_column(points, x, "x");
        }
    }
    else {
        cpl_msg_error(recipe, "Missing column: %s", x);
        montecarlo_exit(NULL);
    }

    cpl_table_subtract_scalar(points, "x", zero);

    if (cpl_table_has_column(table, y)) {
        cpl_table_move_column(points, y, table);
        if (!cpl_table_has_column(points, "y")) {
            cpl_table_name_column(points, y, "y");
        }
    }
    else {
        cpl_msg_error(recipe, "Missing column: %s", y);
        montecarlo_exit(NULL);
    }

    if (sigma[0] != '\0') {
        if (cpl_table_has_column(table, sigma)) {
            cpl_table_move_column(points, sigma, table);
            if (!cpl_table_has_column(points, "y_err")) {
                cpl_table_name_column(points, sigma, "y_err");
            }
        }
        else {
            cpl_msg_error(recipe, "Missing column: %s", sigma);
            montecarlo_exit(NULL);
        }
    }

    cpl_table_delete(table); table = NULL;
    cpl_table_erase_invalid(points);

    count = 1;
    while (start + step*count <= end)
        count++;

    model_error = cpl_table_new(count);
    cpl_table_new_column(model_error, "x", CPL_TYPE_DOUBLE);
    cpl_table_fill_column_window_double(model_error, "x", 0, count, 0.0);
    xdata = cpl_table_get_data_double(model_error, "x");

    for (i = 0; i < count; i++)
        xdata[i] = start + step*i - zero;

    p = mos_montecarlo_polyfit(points, model_error, trials, order);

    if (cpl_error_get_code() != CPL_ERROR_NONE) {
        cpl_msg_error(recipe, "%s", cpl_error_get_message());
        montecarlo_exit(NULL);
    }

    for (i = 0; i < count; i++)
        xdata[i] += zero;

    cpl_polynomial_delete(p); p = NULL;
    cpl_table_delete(points); points = NULL;

    if (dfs_save_table(frameset, model_error, model_error_tag, NULL,
                       parlist, recipe, version))
        montecarlo_exit(NULL);

    cpl_table_delete(model_error); model_error = NULL;

    return 0;
}