#include <math.h>
#include <Eigen/Dense>
#include <iostream>
#include <iterator>
#include <vector>
#include <algorithm>
#include <functional>
#include <Rcpp.h>
#include "zeros.h"
#include "density_estimation.h"
#include "classData.h"

void
dataManager::readData
(const Eigen::Block<Eigen::Map<Eigen::Matrix<double, -1, -1>,0, Eigen::Stride<0, 0> >, 1, -1, false> & row,
       PRIOR prior, const int & cancel)
{
  numbers.clear();

  if( (row.array()!=0.0).any() )
    BM(numbers, row, prior);
  else
  {
    for(unsigned int i = 0, n = row.size(); i < n; i++)
      numbers.push_back(row(i));
  }

  if(cancel!=-1) numbers.erase(numbers.begin() + cancel);

  howmanyclasses = numbers.size();
};

void
dataManager::transfData
()
{
  // clr transformation of prop_data in transf_data
  double a = 0.0;

  // computing geometric mean
  for (const auto& y:numbers)
    a += log(y);

  // clr transformation
  for (auto& y:numbers)
    y = log(y) - a/howmanyclasses;
};

std::vector<double>
dataManager::getNumbers
()
{
  return numbers;
}

void
dataManager::pacs
(densityEstimator & dens, Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> bspline)
{
  dens.solve(bspline,numbers);
};

void
dataManager::antitData
(Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> x)
{
  // Using trapezoidal integration in continuous setting
  double len = (grid.back() - grid.front())/(grid.size()-1);
  unsigned int x_len = x.size();
  double den = 0.5*exp(x(0))*len + 0.5*exp(x(x_len-1))*len;
  for(int i=1; i<x_len-1;i++){
    den += exp(x(i))*len;
  }

  for(int i=0; i<x.size();i++){
    x(i) = exp(x(i))/den;
  }
};

void
dataManager::fillGrid
(double start, double end, unsigned int numPoints)
{
  double step = (end - start)/numPoints;
  grid.resize(numPoints);
  grid[0] = start;
  for (unsigned int i = 1; i < numPoints-1 ; ++i) {
    grid[i] = grid[i-1] + step;
  }
  grid[numPoints-1] = end;
};

void
dataManager::plotData
(const densityEstimator & dens, unsigned long int numPoints,
  Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> bspline,
  Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> yplot)
{

  double start = dens.get_u();
  double end = dens.get_v();
  unsigned int degree = dens.get_k();
  unsigned int G = dens.get_G();
  const std::vector<double> knots = dens.get_lambda();

  fillGrid(start, end, numPoints);

  Eigen::ArrayXd N;

  for (int i = 0; i < grid.size(); ++i)
  {
    int j = bspline::findspan(degree,grid[i],knots);
    N = Eigen::ArrayXd::Constant(G, 0.0);
    bspline::basisfun(j, grid[i], degree, knots, N);
    long double fvalue = compute_fvalue(bspline, N);
    yplot(i)=fvalue;
  }

  antitData(yplot);
  return;
};

void
dataManager::plotData_Clr
(const densityEstimator & dens, unsigned long int numPoints,
  Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> bspline,
  Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> yplot)
{

  double start = dens.get_u();
  double end = dens.get_v();
  unsigned int degree = dens.get_k();
  unsigned int G = dens.get_G();
  const std::vector<double> knots = dens.get_lambda();

  fillGrid(start, end, numPoints);

  Eigen::ArrayXd N;

  for (int i = 0; i < grid.size(); ++i) {
    int j = bspline::findspan(degree,grid[i],knots);
    N = Eigen::ArrayXd::Constant(G, 0.0);
    bspline::basisfun(j, grid[i], degree, knots, N);
    long double fvalue = compute_fvalue(bspline, N);
    yplot(i)=fvalue;
  }

  return;
};

long double
dataManager::compute_fvalue
(Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> vec1, Eigen::ArrayXd vec2 )
{
  long double res = 0.0;
  unsigned int n = vec1.size();
  if(vec2.size() != n){
    Rcpp::Rcerr << "Error in compute_fvalue function. Check dimensions of the vectors.."
              << std::endl;
    Rcpp::stop("Error in the C++ execution");
    // exit(1);
  }
  for (int i = 0; i < n; ++i) {
    res += vec1[i]*vec2[i];
  }
  return res;
};