#include "density_estimation.h"
#include "classData.h"
#include "zeros.h"
#include <vector>
#include <string>
#include <cmath>
#include <iostream>
#include <sstream>
#include <iterator>
#include <RcppEigen.h>
#include <Rcpp.h>
#include <ctime>
#ifdef _OPENMP
  #include <omp.h>
#else
  #define omp_get_thread_num() 0
#endif
#include "timer.h"
#include <chrono>



using namespace Rcpp;
/*! \file
 @brief Main functions that read data from R input and retun the density estimation function
 */

extern "C"{
/*!
 @brief	Main function that read data from R input, transform it and process it to estimate the density function

 @param k_ Spline degree
 @param l_ Derivative degree
 @param alpha_ Penalization parameter
 @param data_ Data for density estimation
 @param Xcp_ Control points
 @param knots_ Knot sequence
 @param weights_ Weights sequence
 @param numPoints_ Number of points of the grid for plotting the density
 @param prior_ 1:Perks, 2: Jeffreys, 3: Bayes-Laplace, 4: Square root
 @see PRIOR
 @param nCPU_ Number of core to use for parallelization
 @param fast_ Unable/disable progress bar

 @return List with B-spline coefficients, points to plot density (also clr-transformed version) and the number of points.

 */
SEXP smoothingSplines_(SEXP k_, SEXP l_, SEXP alpha_, SEXP data_, SEXP Xcp_, SEXP knots_, SEXP weights_, SEXP numPoints_, SEXP prior_, SEXP nCPU_, SEXP fast_)
{
  cns::timer<> t1;
  t1.start();

  #ifdef _OPENMP
    omp_set_dynamic(0);         // Explicitly disable dynamic teams
    omp_set_num_threads(INTEGER(nCPU_)[0]);
  #endif

  bool furious = INTEGER(fast_)[0];

  // Read parameters
  unsigned int k = INTEGER(k_)[0];
  unsigned int l = INTEGER(l_)[0];
  double alpha = REAL(alpha_)[0];
  unsigned int numPoints = INTEGER(numPoints_)[0];
  unsigned int prior_num = INTEGER(prior_)[0];

  PRIOR prior = PRIOR::DEFAULT;

  switch(prior_num)
  {
    case 1: // "perks"
      prior = PRIOR::PERKS;
      break;
    case 2: // "jeffreys"
      prior = PRIOR::JEFFREYS;
      break;
    case 3: // "bayes_laplace"
      prior = PRIOR::BAYES_LAPLACE;
      break;
    case 4:  // "sq"
      prior = PRIOR::SQ;
      break;
    default: {};
  }

  dataManager obj;
  densityEstimator dens(parametersManager(k,l,alpha));

  // Read xcp
  double *Xcp = REAL(Xcp_);
  unsigned int Xcpsize = LENGTH(Xcp_);
  dens.readXcp(Xcp,Xcpsize);

  // Read knots
  double *knots = REAL(knots_);
  unsigned int knotsSize = LENGTH(knots_);
  dens.readKnots(knots,knotsSize);

  // Read data
  Eigen::Map<Eigen::MatrixXd> data(as<Eigen::Map<Eigen::MatrixXd>> (data_));
  unsigned int nrow = data.rows();
  furious = furious || (nrow < 100); // if not useful, progress bar will not be shown
  dens.set_matrix();

  // Read weights
  Eigen::Map<Eigen::MatrixXd> weights(as<Eigen::Map<Eigen::MatrixXd>> (weights_));

  Eigen::MatrixXd bsplineMat(nrow,dens.get_G());
  Eigen::MatrixXd yvalueMat(nrow,numPoints);
  Eigen::MatrixXd yvalueMatClr(nrow,numPoints);

  int barWidth = 70; // for progress bar plot
  int count = 0;
  int pos = 0;

  #pragma omp parallel private(obj) firstprivate(dens)// default(shared)
  {
  #pragma omp for
    for(std::size_t i = 0; i < nrow; i++)
    {
      obj.readData(data.row(i), prior);
      obj.transfData();
      dens.set_weights(weights.row(i));
      obj.pacs(dens, bsplineMat.row(i));
      obj.plotData(dens, numPoints, bsplineMat.row(i), yvalueMat.row(i));
      obj.plotData_Clr(dens, numPoints, bsplineMat.row(i), yvalueMatClr.row(i));

      // Progress bar
      if(!furious)
      {
        #pragma omp critical
        {
          Rcpp::Rcout << "[";
          pos = barWidth * (double)count/nrow;
          for (int j = 0; j < barWidth; ++j)
          {
            if (j < pos) Rcpp::Rcout << "=";
            else if (j == pos) Rcpp::Rcout << ">";
            else Rcpp::Rcout << " ";
          }
          Rcpp::Rcout << "] " << int((double)count/(nrow-1) * 100.0) << "%\r";
          Rcpp::Rcout.flush();
          count++;
        }
      }
    }
  }

  List result = List::create(Named("bspline") = bsplineMat,
                             Named("Y") = yvalueMat,
                             Named("Y_clr") = yvalueMatClr,
                             Named("Xcp") = Xcp_,
                             Named("NumPoints") = numPoints_);
  t1.stop();
  Rcpp::Rcout << "\nIt took "<< t1.elapsed<std::chrono::milliseconds>() <<" milliseconds. " << std::endl;


  return wrap(result);
};


/*!
 @brief	Cross-validation function to choose the best \f$ \alpha \f$

 @param k_ Spline degree
 @param l_ Derivative degree
 @param alpha_ Penalization parameter vector
 @param data_ Data for density estimation
 @param Xcp_ Control points
 @param knots_ Knot sequence
 @param weights_ Weights sequence
 @param numPoints_ Number of points of the grid for plotting the density
 @param prior_ 1:Perks, 2: Jeffreys, 3: Bayes-Laplace, 4: Square root
 @param nCPU_ Number of core to use for parallelization

 @return List with the best \f$ \alpha \f$, the related value of the functional and error

 */
SEXP smoothingSplinesValidation_(SEXP k_, SEXP l_, SEXP alpha_, SEXP data_, SEXP Xcp_, SEXP knots_, SEXP weights_, SEXP prior_, SEXP nCPU_)
{
  cns::timer<> t;
  t.start();

  #ifdef _OPENMP
    omp_set_dynamic(0);         // Explicitly disable dynamic teams
    omp_set_num_threads(INTEGER(nCPU_)[0]);
  #endif

  // Read parameters
  unsigned int k = INTEGER(k_)[0];     // Spline degree
  unsigned int l = INTEGER(l_)[0];
  double * alpha = REAL(alpha_);  // penalization parameters vector
  unsigned int alpha_size= LENGTH(alpha_);
  unsigned int prior_num = INTEGER(prior_)[0];

  PRIOR prior = PRIOR::DEFAULT;

  switch(prior_num)
  {
    case 1: // "perks"
      prior = PRIOR::PERKS;
      break;
    case 2: // "jeffreys"
      prior = PRIOR::JEFFREYS;
      break;
    case 3: // "bayes_laplace"
      prior = PRIOR::BAYES_LAPLACE;
      break;
    case 4:  // "sq"
      prior = PRIOR::SQ;
      break;
    default: {};
  }

  dataManager obj;
  densityEstimator dens(parametersManager(k,l,alpha[0]));

  // Read xcp
  double *Xcp = REAL(Xcp_);
  unsigned int Xcpsize = LENGTH(Xcp_);

  // Read knots
  double *knots = REAL(knots_);
  unsigned int knotsSize = LENGTH(knots_);
  dens.readKnots(knots,knotsSize);

  // Read data
  Eigen::Map<Eigen::MatrixXd> data(as<Eigen::Map<Eigen::MatrixXd>> (data_));
  // Read weights
  Eigen::Map<Eigen::MatrixXd> weights(as<Eigen::Map<Eigen::MatrixXd>> (weights_));

  unsigned int nrow = data.rows();
  unsigned int ncol = data.cols();

  int barWidth = 70; // for bar progress plot
  int count = 0;
  int pos = 0;
  unsigned int min_index = 0;
  Eigen::VectorXd Jvalues = Eigen::VectorXd::Zero(alpha_size);
  Eigen::VectorXd CVerror = Eigen::VectorXd::Zero(alpha_size,0.0);
  double CVopt = -1.0;
  #pragma omp parallel private(obj) firstprivate(dens)// default(shared)
  {
    Eigen::MatrixXd threadBsplineMat(nrow,dens.get_G());
    Eigen::ArrayXd N;
    int span;
    long double fvalue;

    #pragma omp for
    for(std::size_t z = 0; z < alpha_size; z++)
    {
      dens.set_alpha(alpha[z]);
      for(std::size_t i = 0; i < nrow; i++)
      {
        for(std::size_t j = 1; j < ncol-1; j++)    // we do not leave out extremals, otherwise we have to change the knots
        {
          // LOO-CV, leave out column j, fit spline and compute error
          dens.readXcp(Xcp,Xcpsize,j);
          dens.set_matrix();

          obj.readData(data.row(i),prior,j);
          obj.transfData();
          dens.set_weights(weights.row(i));
          obj.pacs(dens, threadBsplineMat.row(i));
          Jvalues(z)+=dens.eval_J(obj.getNumbers())/nrow;

          // computing error using spline coefficients, fvalue is the predicted value
          span = bspline::findspan(k,Xcp[j],dens.get_lambda());
          N = Eigen::ArrayXd::Constant(dens.get_G(), 0.0);
          bspline::basisfun(span,Xcp[j],k,dens.get_lambda(), N);
          fvalue = obj.compute_fvalue(threadBsplineMat.row(i), N);
          CVerror(z) += (fvalue - obj.getNumbers()[j])*(fvalue - obj.getNumbers()[j])/nrow;
        }
      }
    }
  }
  for(std::size_t z = 0; z < nrow; z++)
  {
    if(CVerror(z)<CVopt || CVopt<0)
    {
      CVopt = CVerror(z);
      min_index = z;
    }
  }

  // Printing out J[alpha]
  for(std::size_t k = 0; k < alpha_size; k++)
  {
    Rcout << std::string((k==min_index)?" ***":" ") << "\talpha = "<< alpha[k] << "\tJ = " << Jvalues(k)
          << "\t CV-error: " << CVerror(k) << '\n';
  }



  List result = List::create(Named("alpha") = alpha_,
                             Named("Jvalues") = Jvalues,
                             Named("CVerror") = CVerror);


  t.stop();
  Rcpp::Rcout << "\nIt took "<< t.elapsed<std::chrono::milliseconds>() <<" milliseconds. " << std::endl;

  return wrap(result);
};
}