#ifndef __COSMO_DINTERPOLATE_HPP
#define __COSMO_DINTERPOLATE_HPP

#include "config.hpp"
#include "mykdtree.hpp"
#include <gsl/gsl_eigen.h>

namespace CosmoTool {
 
  template<typename PType, typename IType, int N>
  class DelaunayInterpolate {
  public:
    struct SimplexAccess {
      int32_t *simplex_list;      
    };
 
    typedef KDTree<N, SimplexAccess, PType> QuickTree;
    typedef typename QuickTree::Cell QuickCell;
    typedef PType CoordType[N];
    
    QuickTree *quickAccess;
    QuickCell *cells;
    PType *all_preweight;
    int32_t *point_to_simplex_list_base;    
    IType *values;
    CoordType *positions;    
    uint32_t numPoints;
    uint32_t numSimplex;
    uint32_t *simplex_list;
    gsl_eigen_symmv_workspace *eigen_work;
    
    /**
     * This construct the interpolator. The construction is time consuming so
     * please do it the less possible number of times, especially if you have
     * a large number of points.
     *
     * @param positions list of the positions
     * @param values list of the values taken at each position
     * @param simplex_list list of points for each simplex. The packing
     * is the following:
     * [t(0,1),t(0,2),...,t(0,n+1),t(1,0),t(1,1),...,t(1,n+1),..],
     * with t(i,j) the i-th simplex and j-th point of the simplex. The indexes
     * refer to the previous list of points.
     * @param numPoints the number of points
     */
    DelaunayInterpolate(CoordType *positions, IType *values, uint32_t *simplex_list,
			uint32_t numPoints, uint32_t numSimplex)
    {
      this->positions = positions;
      this->values = values;
      this->simplex_list = simplex_list;
      this->numPoints = numPoints;
      this->numSimplex = numSimplex;
      
      buildPreweight();
      buildQuickAccess();

      eigen_work = gsl_eigen_symmv_alloc(N);

    }
    
    ~DelaunayInterpolate()
    {
      delete[] cells;
      delete quickAccess;
      delete[] point_to_simplex_list_base;      
      delete[] all_preweight;

      gsl_eigen_symmv_free(eigen_work);
    }

    void buildPreweight();
    void buildQuickAccess();
    void buildHyperplane(const PType *v, CoordType& hyper);

    bool checkPointInSimplex(const CoordType& pos, uint32_t simplex);
    
    uint32_t findSimplex(const CoordType& pos)
      throw (InvalidArgumentException);

    IType computeValue(const CoordType& pos)
      throw (InvalidArgumentException);
  };
};

#include "dinterpolate.tcc"

#endif