/* -*- Mode: C; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
/* vim:set sw=2 sts=2 et cin: */
/*
 * This file is part of the MUSE Instrument Pipeline
 * Copyright (C) 2007-2011 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 Street, Fifth Floor, Boston, MA  02110-1301, USA.
 */

#define _DEFAULT_SOURCE
#define _BSD_SOURCE /* get mkdtemp() from stdlib.h */
#include <unistd.h> /* rmdir() */
#include <stdio.h> /* remove() */
#include <sys/stat.h> /* mkdir() */
#include <string.h>

#include <muse.h>
#include <muse_instrument.h>

#define DIR_TEMPLATE "/tmp/muse_sky_master_test_XXXXXX"

static void
test_lsf_params(const char *aDir)
{
  muse_lsf_params **par 
    = cpl_calloc(kMuseNumIFUs+1, sizeof(muse_lsf_params *));

  // Fill with some real-world specs.
  cpl_size i;
  for (i = 0; i < kMuseNumIFUs; i++) {
    par[i] = cpl_malloc(sizeof(muse_lsf_params));
    par[i]->ifu = 1;
    par[i]->slice = i+1;
    par[i]->lambda_ref = 7000.;
    par[i]->sensitivity = cpl_array_new(2, CPL_TYPE_DOUBLE);
    cpl_array_set(par[i]->sensitivity, 0, 1.0);
    cpl_array_set(par[i]->sensitivity, 1, 0.0);
    par[i]->slit_width = 2.51;
    par[i]->bin_width = 1.30;
    par[i]->lsf_width = cpl_array_new(3, CPL_TYPE_DOUBLE);
    cpl_array_set(par[i]->lsf_width, 0, 0.444);
    cpl_array_set(par[i]->lsf_width, 1, -7.77e-6);
    cpl_array_set(par[i]->lsf_width, 2, 6e-8);
    par[i]->hermit[0] = cpl_array_new(1, CPL_TYPE_DOUBLE);
    par[i]->hermit[1] = cpl_array_new(2, CPL_TYPE_DOUBLE);
    par[i]->hermit[2] = cpl_array_new(1, CPL_TYPE_DOUBLE);
    par[i]->hermit[3] = cpl_array_new(2, CPL_TYPE_DOUBLE);
    cpl_array_set(par[i]->hermit[0], 0, -0.47);
    cpl_array_set(par[i]->hermit[1], 0, 1.777);
    cpl_array_set(par[i]->hermit[1], 1, 3.6e-4);
    cpl_array_set(par[i]->hermit[2], 0, -0.3);
    cpl_array_set(par[i]->hermit[3], 0, 2.86);
    cpl_array_set(par[i]->hermit[3], 1, 3.4e-4);
  }

  char *sfilename = cpl_sprintf("%s/sky_par01.fits", aDir);
  char *sfilename2 = cpl_sprintf("%s/sky_par02.fits", aDir);

  cpl_propertylist_save(NULL, sfilename, CPL_IO_CREATE);
  cpl_test(muse_lsf_params_save((const muse_lsf_params **)par,
                                sfilename) == CPL_ERROR_NONE);
  for (i = 0; i < kMuseNumIFUs; i++) {
    par[i]->ifu = 2;
  }
  cpl_propertylist_save(NULL, sfilename2, CPL_IO_CREATE);
  cpl_test(muse_lsf_params_save((const muse_lsf_params **)par,
                                sfilename2) == CPL_ERROR_NONE);

  muse_lsf_params **spar = muse_lsf_params_load(sfilename, NULL, 0);
  cpl_test(spar != NULL);
  if (spar != NULL) {
    cpl_test(spar[0] != NULL);
    if (spar[0] != NULL) {
      cpl_test(spar[0]->sensitivity != NULL);
      cpl_test(spar[0]->lsf_width != NULL);
      cpl_test(spar[0]->hermit[0] != NULL);
    }
    cpl_test(muse_lsf_params_get_size(spar) == kMuseNumIFUs);
  }

  spar = muse_lsf_params_load(sfilename2, spar, 0);
  cpl_test(spar != NULL);
  if (spar != NULL) {
    cpl_test(spar[0] != NULL);
    if (spar[0] != NULL) {
      cpl_test(spar[0]->sensitivity != NULL);
      cpl_test(spar[0]->lsf_width != NULL);
      cpl_test(spar[0]->hermit[0] != NULL);
    }
    cpl_test(muse_lsf_params_get_size(spar) == 2 * kMuseNumIFUs);
  }

  muse_lsf_params_delete_all(spar);

  spar = muse_lsf_params_load(sfilename, NULL, 2);
  cpl_test(spar != NULL);
  if (spar != NULL) {
    cpl_test(muse_lsf_params_get_size(spar) == 0);
  }

  spar = muse_lsf_params_load(sfilename2, spar, 2);
  cpl_test(spar != NULL);
  if (spar != NULL) {
    cpl_test(spar[0] != NULL);
    if (spar[0] != NULL) {
      cpl_test(spar[0]->sensitivity != NULL);
      cpl_test(spar[0]->lsf_width != NULL);
      cpl_test(spar[0]->hermit[0] != NULL);
    }
    cpl_test(muse_lsf_params_get_size(spar) == kMuseNumIFUs);
  }

  // Test that the values are within the some reasonable values. Also test that
  // we can use values within a certain range (+-20).
  cpl_array *lsf = cpl_array_new(5, CPL_TYPE_DOUBLE);
  cpl_array_set_double(lsf, 0, -20);
  cpl_array_set_double(lsf, 1, -2);
  cpl_array_set_double(lsf, 2, 0.1);
  cpl_array_set_double(lsf, 3, 2.7);
  cpl_array_set_double(lsf, 4, 20);

  double lambda = 6000;
  cpl_array *val = cpl_array_duplicate(lsf);
  muse_lsf_params_apply(par[0], val, lambda);
  int res;
  cpl_test_leq(cpl_array_get(val, 0, &res), 1e-3);
  cpl_test_eq(res, 0);
  cpl_test_leq(cpl_array_get(val, 1, &res), 1e-1);
  cpl_test_eq(res, 0);
  cpl_test_leq(1e-1, cpl_array_get(val, 2, &res));
  cpl_test_eq(res, 0);
  cpl_test_leq(cpl_array_get(val, 3, &res), 1e-1);
  cpl_test_eq(res, 0);
  cpl_test_leq(cpl_array_get(val, 4, &res), 1e-3);
  cpl_test_eq(res, 0);

  cpl_array_delete(val);
  cpl_array_delete(lsf);

  // Test that the LSF is well centered and normalized.
  cpl_size n_val = 1000;
  double step = 20./(n_val - 1);
  lsf = cpl_array_new(n_val, CPL_TYPE_DOUBLE);
  for (i = 0; i < n_val; i++) {
    cpl_array_set(lsf, i, -10. + step * i);
  }
  for (lambda = 4500; lambda < 9500; lambda += 122) {
    val = cpl_array_duplicate(lsf);
    muse_lsf_params_apply(par[0], val, lambda);
    cpl_array * v = cpl_array_duplicate(lsf);
    cpl_test_abs(cpl_array_get_mean(val) * n_val * step, 1.0, 1e-3);
    cpl_array_multiply(v, val);
    cpl_test_abs(cpl_array_get_mean(v)/cpl_array_get_mean(val), 0.0, 1e-4);
    cpl_array_delete(val);
    cpl_array_delete(v);
  }
  cpl_array_delete(lsf);

  muse_lsf_params_delete_all(par);
  muse_lsf_params_delete_all(spar);

  remove(sfilename);
  remove(sfilename2);
  cpl_free(sfilename);
  cpl_free(sfilename2);
}

static void
test_lsf_params_fit(const char *aDir)
{
  UNUSED_ARGUMENT(aDir);
  char *sfilename = cpl_sprintf(BASEFILENAME"_pixtable.fits");
  muse_pixtable *pixtable = cpl_calloc(1, sizeof(muse_pixtable));
  pixtable->table = cpl_table_load(sfilename, 1, 0);
  /* header can stay empty for this test */
  cpl_free(sfilename);
  if (pixtable->table == NULL) {
    muse_pixtable_delete(pixtable);
    return;
  }

  cpl_size i_row;
  cpl_size n_rows = cpl_table_get_nrow(pixtable->table);
  cpl_size n_lines = 0;
  double llambda = -1;
  int ifu = 7;
  int slice = 12;
  cpl_table *lines = cpl_table_new(1);
  cpl_table_new_column(lines, "lambda", CPL_TYPE_FLOAT);
  cpl_table_new_column(lines, "flux", CPL_TYPE_FLOAT);
  for (i_row = 0; i_row < n_rows; i_row++) {
    uint32_t origin = (uint32_t)cpl_table_get_int(pixtable->table,
                                                  MUSE_PIXTABLE_ORIGIN,
                                                  i_row, NULL);
    if (muse_pixtable_origin_get_ifu(origin) != ifu) {
      continue;
    }
    if ( muse_pixtable_origin_get_slice(origin) != slice) {
      continue;
    }
    double line_lambda = cpl_table_get(pixtable->table, "line_lambda",
                                       i_row, NULL);
    if (llambda == line_lambda) {
      continue;
    }
    cpl_size i_line;
    for (i_line = 0; i_line < n_lines; i_line++) {
      if (cpl_table_get(lines, "lambda", i_line, NULL) == line_lambda) {
        break;
      }
    }
    if (i_line < n_lines) {
      continue;
    }
    llambda = line_lambda;
    n_lines++;
    cpl_table_set_size(lines, n_lines);
    cpl_table_set(lines, "lambda", n_lines-1, line_lambda);
    // we dont use the calculated flux here since it should be
    // calculated by themself.
    cpl_table_set(lines, "flux", n_lines-1, 1.0);
  }
  cpl_table_erase_column(pixtable->table, "line_lambda");
  cpl_table_erase_column(pixtable->table, "line_flux");

  /* sort list by decreasing flux and cut it to the 20 brightest lines */
  cpl_propertylist *flux_sort = cpl_propertylist_new();
  cpl_propertylist_append_bool(flux_sort, "flux", CPL_TRUE);
  cpl_table_sort(lines, flux_sort);
  cpl_propertylist_delete(flux_sort);
  cpl_size nrows = 20;
  if (nrows < cpl_table_get_nrow(lines)) {
    cpl_table_erase_window(lines, nrows, cpl_table_get_nrow(lines));
  }

  muse_lsf_params *lsfParams = muse_lsf_params_fit(pixtable, lines, 40);
  cpl_test_nonnull(lsfParams);

  cpl_size n_val = 100;
  double step = 12./(n_val - 1);
  cpl_array *lsf = cpl_array_new(n_val, CPL_TYPE_DOUBLE);
  cpl_size i;
  for (i = 0; i < n_val; i++) {
    cpl_array_set(lsf, i, -6. + step * i);
  }
  double lambda;
  for (lambda = 4500; lambda < 9500; lambda += 122) {
    cpl_array *val = cpl_array_duplicate(lsf);
    muse_lsf_params_apply(lsfParams, val, lambda);
    // Check for some specific points.(left border, center, right border)
    cpl_test_lt(cpl_array_get(val, 0, NULL), 1e-2);
    cpl_test_lt(cpl_array_get(val, n_val-1, NULL), 1e-2);
    cpl_test_lt(0.2, cpl_array_get(val, n_val/2, NULL));

    // Test that the function raises on the left half and decreases
    // on the right half.
    double old = 0;
    int up = 1;
    for (i = 0; i < n_val; i++) {
      double z = cpl_array_get(val, i, NULL);
      if (i < 0.45 * n_val) {
        cpl_test_leq(old, z + 1e-2);
      } else if (!up || (i > 0.55 * n_val)) {
        cpl_test_leq(z, old + 1e-2);
      } else if (up && (old > z)) {
        up = 0;
      }
      old = z;
    }

    // Test for normalization
    cpl_test_abs(cpl_array_get_mean(val) * n_val * step, 1.0, 1e-3);

    // Test that the LSF is well centered.
    cpl_array * v = cpl_array_duplicate(lsf);
    cpl_array_multiply(v, val);
    cpl_test_abs(cpl_array_get_mean(v)/cpl_array_get_mean(val), 0.0, 1e-4);
    cpl_array_delete(val);
    cpl_array_delete(v);
  }
  cpl_array_delete(lsf);

  muse_lsf_params_delete(lsfParams);
  cpl_table_delete(lines);
  muse_pixtable_delete(pixtable);
}

/*----------------------------------------------------------------------------*/
/**
  @brief    short test program to check sky subtraction related stuff.
 */
/*----------------------------------------------------------------------------*/
int main(int argc, char **argv)
{
  UNUSED_ARGUMENTS(argc, argv);
  cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_INFO);

  char dirtemplate[FILENAME_MAX] = DIR_TEMPLATE;
  char *dir = mkdtemp(dirtemplate);

  test_lsf_params(dir);
  test_lsf_params_fit(dir);

  cpl_test_zero(rmdir(dir));
  return cpl_test_end(0);
}