/** File wavelet.c
 *  Compute Wavelet transform
 */

#include <stdlib.h>
#include <string.h>

#include "defines.h"
#include "graphs.h"
#include "utils.h"
#include "noxprotos.h"
#include "ssdata.h"

#include "wavelet.h"


static WAVParams wavparams;


void wavelet_set_defaults (void)
{
  wavparams.invflag = 0;
  wavparams.typ     = WAVELET_DAUBECHIES;
  wavparams.stride  = 1;
}

int wavelet_params_get (WAVParams *p)
{
  if (p != NULL) {
    memcpy (p ,&wavparams ,sizeof(WAVParams));
    return RETURN_SUCCESS;
  } else {
    return RETURN_FAILURE;
  }
}

int wavelet_params_set (WAVParams *p)
{
  if (p != NULL) {
    memcpy (&wavparams ,p ,sizeof(WAVParams));
    return RETURN_SUCCESS;
  } else {
    return RETURN_FAILURE;
  }
}

int do_wavelet (int gno1 ,int setno1 ,int gno2 ,int setno2)
{
  int i ,inlen ,outlen ,i2 ,k;
  double *in_x ,*in_y ,*out_x ,*out_y;
  double xspace;
  gsl_wavelet *wav;
  gsl_wavelet_workspace *work;
  gsl_wavelet_direction dir;
  char buf[256];
  /* get input and perform safety check on it */
  inlen = getsetlength (gno1 ,setno1);
  if (inlen < 2) {
    errmsg("Set length < 2");
    return RETURN_FAILURE;
  }
  if ((i2 = ilog2 (inlen)) <= 0) {
    errmsg("Set length not a power of 2");
    return RETURN_FAILURE;
  }
  in_x = getcol (gno1 ,setno1 ,DATA_X);
  in_y = getcol (gno1 ,setno1 ,DATA_Y);
  if (monospaced (in_x ,inlen ,&xspace) != TRUE) {
    errmsg("do_wavelet: abscissas not monospaced.\n\t (you can interpolate the set before doing wavelet transform if you want)");
    return RETURN_FAILURE;
  } else {
    if (xspace == 0.0) {
      errmsg("do_wavelet: The set spacing is 0, can't continue");
      return RETURN_FAILURE;
    }
  }  
  /* Reallocate set(s) */
  outlen = inlen;
  setlength (gno2 ,setno2 ,outlen);
  out_x = getx (gno2 ,setno2);
  out_y = gety (gno2 ,setno2);
  /* Compute abscissas */
  for (i = 0; i < outlen; i++) out_x[i] = (double) i;
  /* Prepare call to gsl */
  k = wavparams.k;
  wav = NULL;
  switch(wavparams.typ) {
  case WAVELET_DAUBECHIES:
    wav = gsl_wavelet_alloc (gsl_wavelet_daubechies ,k);
    break;
  case WAVELET_DAUBECHIES_CENTERED:
    wav = gsl_wavelet_alloc (gsl_wavelet_daubechies_centered ,k);
    sprintf (buf, "centered Daubechies wavelet on G%d.s%d", gno1 ,setno1);
    break;
  case WAVELET_HARR:
    wav = gsl_wavelet_alloc (gsl_wavelet_haar ,2);
    break;
  case WAVELET_HARR_CENTERED:
    wav = gsl_wavelet_alloc (gsl_wavelet_haar_centered ,k);
    break;
  case WAVELET_BSPLINE:
    wav = gsl_wavelet_alloc (gsl_wavelet_bspline ,k);
    break;
  case WAVELET_BSPLINE_CENTERED:
    wav = gsl_wavelet_alloc (gsl_wavelet_bspline_centered ,k);
    break;
  }
  work = gsl_wavelet_workspace_alloc (inlen);
  dir = (wavparams.invflag == TRUE) ? gsl_wavelet_backward : gsl_wavelet_forward;
    
  /* EFFECTIVE COMPUTATION of the Wavelet transform */
  gsl_wavelet_transform (wav ,out_y ,wavparams.stride ,inlen ,dir ,work);
  
    
  /* Write appropriate comment */
  sprintf (buf, "%s coef. for G%d.s%d" ,gsl_wavelet_name (wav) ,gno1 ,setno1);
  comment_and_legend_set (gno2 ,setno2 ,buf);

  gsl_wavelet_free (wav);
  gsl_wavelet_workspace_free (work); 
  return RETURN_SUCCESS;
}

void do_wavelet_command (int gno1 ,int setno1 ,WAVELET_type typ ,int stride ,int k ,int invflag)
{
  wavparams.typ     = typ;
  wavparams.invflag = invflag;
  wavparams.stride  = stride;
  wavparams.k       = k;
  int setno2  = nextset (gno1);
  do_copyset (gno1 ,setno1 ,gno1 ,setno2);
  do_wavelet (gno1 ,setno1 ,gno1 ,setno2);
#ifdef WITH_GTK
  gg_wavelet_update ();
#endif

}
