/***************************************************************************
 *   Copyright (C) 2007 by Faubet Pierre   *
 *   pierre.faubet@e.ujf-grenoble.fr   *
 *                                                                         *
 *   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.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/


#include <QReadLocker>

#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>

#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>

#include "BIMrErrors.h"
#include "rawdata.h"

QReadWriteLock RawData::lock ;
int RawData::N = 512;

RawData::RawData(string filename,int column,int size)
{
  /* Sample size */
  n = size;

  /* Create array */
  iter = new double [n];
  if (!iter)
    {
      EXCEPTION_INFOS();
      throw MemError();
    }

  obs = new double [n];
  if (!obs)
    {
      EXCEPTION_INFOS();
      throw MemError();
    }

  x_kde = new double [N];
  if (!x_kde)
    {
      EXCEPTION_INFOS();
      throw MemError();
    }

  y_kde = new double [N];
  if (!y_kde)
    {
      EXCEPTION_INFOS();
      throw MemError();
    }

  QReadLocker locker(&lock );

  /* Read data from file */
  clog << "Get raw data from file " << filename.data() << endl;

  ifstream datafile(filename.data());
  string line;
  getline(datafile,line);
  istringstream colname_str(line);
  /* Get column title */
  colname_str >> colname;
  for (int col=1;col<=column;col++)
    colname_str >> colname;
  /* Get sample */
  for (int i=0;i<n;i++)
    {
      getline(datafile,line);
      istringstream data(line);
      data >> iter[i];
      obs[i] = iter[i];
      for (int col=1;col<=column;col++)
        data >> obs[i];
    }
    clog << "Get raw data from file " << filename.data() << " ... (done)" << endl;
		    
  /* Compute statistics */
  getStatistics();
  getKde(Epanechnikov);
  getHpdi();
}

RawData::RawData(string filename,int filter,int column,int size)
{
  /* Sample size */
  n = size;

  iter = new double [n];
  if (!iter)
    cerr << "RawData error: can not allocate memory" << endl;

  obs = new double [n];
  if (!obs)
    cerr << "RawData error: can not allocate memory" << endl;

  x_kde = new double [N];
  if (!x_kde)
    cerr << "RawData error: can not allocate memory" << endl;

  y_kde = new double [N];
  if (!y_kde)
    cerr << "RawData error: can not allocate memory" << endl;

  QReadLocker locker(&lock );

  /* Read data from file */
  clog << "Get raw data from file " << filename.data() << endl;

  ifstream datafile(filename.data());
  string line;
  getline(datafile,line);
  istringstream colname_str(line);
  /* Get column title */
  colname_str >> colname;
  for (int col=1;col<=column;col++)
    colname_str >> colname;
  /* Get sample */
  int index=0;
  for (int i=0;i<n;i++)
    {
      getline(datafile,line);
      istringstream data(line);
      /* Get iteration number */
      data >> iter[index];
      obs[index] = iter[index];

      /* Get model */
      int model;
      data >> model;
      if (model == filter)
        {
          data >> obs[index];
          for (int col=3;col<=column;col++)
            data >> obs[index];
          index++;
        }
    }

  clog << "Compute statistics for parameter " << colname.data() << " under model " << filter << endl;
  /* Sample size under the model */
  n = index;
  if (index != 0)
    {
      /* Compute statistics */
      getStatistics();
      getKde(Epanechnikov);
      getHpdi();
    }
}

RawData::~RawData()
{
  delete [] x_kde;
  delete [] y_kde;
  delete [] obs;
  delete [] iter;
}
//
void RawData::getStatistics(int burnin,int burnout)
{
  double *data = new double [n-(burnin+burnout)];

  mean=0;
  sd=0;
  for (int i=0;i<n-(burnin+burnout);i++)
    {
      data[i] = obs[burnin+i];
      mean += data[i];
      sd += data[i]*data[i];
    }
  mean /= n-(burnin+burnout);
  sd = sqrt(sd/(n-(burnin+burnout)) - mean*mean);

  gsl_sort(data,1,n-(burnin+burnout));
  min = data[0];
  max = data[n-(burnin+burnout)-1];
  double q1 = gsl_stats_quantile_from_sorted_data(data,1,n-(burnin+burnout),0.25);
  double q3 = gsl_stats_quantile_from_sorted_data(data,1,n-(burnin+burnout),0.75);
  /* Bandwidth for kernel density estimation */
  h0 = 1.36374*std::min(sd,(q3-q1)/1.34)*pow(n-(burnin+burnout),-1.0/5.0);
  
  delete [] data;
}

void RawData::getKde(kernel_t kernel,double adjust,int burnin,int burnout)
{
  /* Bandwidth for kernel density estimation */
  double bandwidth = h0*adjust;
  switch (kernel)
    {
    case Gaussian:
      bandwidth *= 0.7764;
      break;
    case Epanechnikov:
      bandwidth *= 1.7188;
      break;
    case Triangular:
      bandwidth *= 1.8882;
      break;
    case Boxcar:
      bandwidth *= 1.3510;
      break;
    default:
      cerr << "Kernel error" << endl;
      break;
    }

  double dx = (max - min) /(N-1);
  mode = min;
  double y_max = 0;
  for (int i=0;i<N;i++)
    {
      x_kde[i] = min + i*dx;
      y_kde[i] = kde(x_kde[i],kernel,bandwidth,obs+burnin,n-(burnin+burnout));

      if (y_kde[i] > y_max)
        {
          y_max = y_kde[i];
          mode = x_kde[i];
        }
    }
}

double RawData::kde(double x_i,kernel_t k_t,double bandwidth, double *x, int size)
{
  if (bandwidth == 0)
    return(x_i == min);

  if ((k_t != Gaussian) && ((x_i < min - bandwidth) || (max + bandwidth < x_i)))
    return(0);

  double sum = 0;
  switch (k_t)
    {
    case Epanechnikov:
      for (int i=0;i<size;++i)
        if ((x[i] - bandwidth <= x_i) && (x_i <= x[i] + bandwidth))
          {
            double t = (x_i-x[i])/bandwidth;
            sum += 1 - t*t;
          }
      sum *= 3.0/4.0;
      break;
    case Boxcar:
      for (int i=0; i<size; ++i)
        if ((x[i] - bandwidth <= x_i) && (x_i <= x[i] + bandwidth))
          sum += 0.5;
      break;
    case Triangular:
      for (int i=0; i<size; ++i)
        if ((x[i] - bandwidth <= x_i) && (x_i <= x[i] + bandwidth))
          {
            double t = (x_i-x[i])/bandwidth;
            sum += 1 - fabs(t);
          }
      break;
    case Gaussian:
      for (int i=0; i<size; ++i)
        {
          double t = (x_i-x[i])/bandwidth;
          sum += exp(-t*t/2);
        }
      sum /= sqrt(2*M_PI);
      break;
    default:
      break;
    }
  return sum/(size*bandwidth);
}

void RawData::getHpdi(double alpha)
{
  int m = ((int) std::max(1.0, ceil(alpha*N)));
  double *a = new double [m];
  if (!a)
    {
      cerr << "Hpdi error" << endl;
      exit(0);
    }
  double *b = new double [m];
  if (!b)
    {
      cerr << "Hpdi error" << endl;
      exit(0);
    }

  int ind = 0;
  double width = x_kde[N-1] - x_kde[0];
  for (int i=0;i<m;i++)
    {
      a[i] = x_kde[i];
      b[i] = x_kde[N-m+i] ;
      if (width > b[i]-a[i])
        {
          width = b[i]-a[i];
          ind = i;
        }
    }
  hpdilo = a[ind];
  hpdihi = b[ind];
  
  delete [] a;
  delete [] b;
}
