// Copyright (C) 2016  Jerome Lelong <jerome.lelong@imag.fr>
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#ifndef _MULTIVARIATEBASIS_H
#define _MULTIVARIATEBASIS_H

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <cassert>


/**
 * \file MultiVariateBasis.h
 * \brief Multivariate function bases for solving regression problems
 * \author Jérôme Lelong
 */

namespace StOpt
{
/**
 * \defgroup MultiVariateBasis Multivariate bases
 * \brief It stores a family of multivariate functions defined as a tensor
 * product of single variate functions. This class provides all the necessary
 * methods for solving regression problems.
 * @{
 */

/**\brief Integer Matrix in row major order */
typedef Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> RowMatrixXi;
/**\brief Sparse integer Matrix in row major order */
typedef Eigen::SparseMatrix<int, Eigen::RowMajor> RowSparseMatrix;

/**
 * \class ComputeDegree
 * \brief Abstract class to compute the p_total degree
 */
class ComputeDegree
{
public:
    virtual ~ComputeDegree() {}
    /**
     * \brief Compute the tensor representation of the given basis
     * \param p_nVariates the number of variates
     * \param p_degree the maximum p_total degree
     */
    virtual RowMatrixXi computeTensor(int p_nVariates, int p_degree) const;

protected:
    /**
     * \brief Compute the p_total degree of a row
     * \param p_tensor a Tensor
     * \param p_row a line number
     */
    virtual int count(const RowMatrixXi &p_tensor, int p_row) const = 0;
    /**
     * \brief Compute the maximum p_partial degree which can be added in order
     * not to go over the specified p_total degree given the running p_total
     * degree
     * \param p_total maximum p_total degree
     * \param p_partial p_partial degree sofar
     */
    virtual int freedom(int p_total, int p_partial) const = 0;
    /**
     * \brief Compute the number of elements with p_total degree less or equal than degree
     * in the basis with (p_tensor->n + 1) variates
     *
     * \param p_tensor the tensor matrix of the basis with n-1 variates
     * \param p_degree the maximum p_total degree requested
     *
     * \return the number of elements with p_total degree less or equal than degree in
     * the basis with n variates
     */
    int computeNumberOfElements(const RowMatrixXi &p_tensor, int p_degree) const;
    /**
     * \brief Copy T_prev(T_previ, :) into p_tensor(Ti,:)
     *
     * \param p_tensor an integer matrix with more columns than that T_prev
     * \param p_previousTensor an integer matrix containing the last tensor computed
     * \param p_tensorRow the index of the line to consider in p_tensor
     * \param p_previousTensorRow the index of the line to consider in T_prev
     *
     */
    void copyPreviousTensor(RowMatrixXi &p_tensor, const RowMatrixXi &p_previousTensor, int p_tensorRow, int p_previousTensorRow) const;
};

/**
 * \class ComputeDegreeSum
 * \brief The p_total degree is the sum of the p_partial degrees
 */
class ComputeDegreeSum : public ComputeDegree
{
protected:
    /**
     * \brief Compute the p_total degree of a row
     * \param p_tensor a Tensor
     * \param p_row a line number
     */
    virtual int count(const RowMatrixXi &p_tensor, int p_row) const;
    /**
     * \brief Compute the maximum p_partial degree which can be added in order
     * not to go over the specified p_total degree given the running p_total
     * degree
     * \param p_total maximum p_total degree
     * \param p_partial p_partial degree sofar
     */
    virtual int freedom(int p_total, int p_partial) const;
};

/**
 * \class ComputeDegreeProd
 * \brief The p_total degree is the product of the p_partial degrees. More
 * precisely,
 * - if not all the p_partial degrees are zero, total_degree= Prod(max(partial_degree, 1))
 * - if all the p_partial degrees are zero, total_degree = 0
 */
class ComputeDegreeProd : public ComputeDegree
{
protected:
    /**
     * \brief Compute the p_total degree of a row
     * \param p_tensor a Tensor
     * \param p_row a line number
     */
    virtual int count(const RowMatrixXi &p_tensor, int p_row) const;
    /**
     * \brief Compute the maximum p_partial degree which can be added in order
     * not to go over the specified p_total degree given the running p_total
     * degree
     * \param p_total maximum p_total degree
     * \param p_partial p_partial degree sofar
     */
    virtual int freedom(int p_total, int p_partial) const;
};

/**
 * \class ComputeDegreeHyperbolic
 * \brief The p_total degree is the result of an hyperbolic semi-norm
 * total_degree = (sum(partial_degree^q)^(1/q)
 */
class ComputeDegreeHyperbolic : public ComputeDegreeSum
{
private:
    double m_q;
    /**
     * \brief Compute the p_total degree of a row
     * \param p_tensor a Tensor
     * \param p_row a line number
     */
    virtual double hyperbolicCount(const RowMatrixXi &p_tensor, int p_row) const;
public:
    ComputeDegreeHyperbolic();
    /**
     * \brief Create an instance of the DegreeHyperbolic class
     * \param p_q the hyperbolic order must be between 0 and 1 to have a sparse
     * representation
     */
    explicit ComputeDegreeHyperbolic(double p_q);
    /**
     * \brief Compute the tensor representation of the given basis
     * \param p_nVariates the number of variates
     * \param p_degree the maximum p_total degree
     */
    virtual RowMatrixXi computeTensor(int p_nVariates, int p_degree) const;
};

/**
 * \class MultiVariateBasis
 * \brief Implement a family of multi-variate functions, each of them being a
 * tensor product of single variate functions.
 *
 * \tparam Func1D A class of 1D functions, which must provide the three
 * following methods
 * - F(double, int) to evaluate the 1d function
 * - DF(double, int) to evaluate the first derivative of the 1d function
 * - D2F(double, int) to evaluate the second derivative of the 1d function
 */
template <class ClassFunc1D> class MultiVariateBasis
{
private:
    ClassFunc1D m_func1D; /*!< Base 1D family */
    int m_numberOfVariates; /*!< number of variates */
    int m_numberOfFunctions; /*!< number of elements in the basis */
    RowMatrixXi m_tensorFull; /*<! Full Basis representation with size m_numberOfFunctions x m_numberOfVariates */
    RowSparseMatrix m_tensorSparse;/*<! Sparse Basis representation with size m_numberOfFunctions x m_numberOfVariates */
    bool m_isReduced; /* true if the basis is reduced */
    Eigen::ArrayXd m_center; /*!< center of the domain */
    Eigen::ArrayXd m_scale;  /*<! inverse of the scaling factor to map the domain to [-1, 1]^nb_variates */


public:
    MultiVariateBasis(): m_numberOfVariates(0), m_numberOfFunctions(0), m_isReduced(false) { }

    /**
     * \brief Create a multi-variate basis
     * \param p_degreeFunction a derived class from ComputeDegree
     * \param p_dimension the number of variates
     * \param p_degree the maximum p_total degree
     */
    MultiVariateBasis(const ComputeDegree &p_degreeFunction, int p_dimension, int p_degree)
    {
        m_numberOfVariates = p_dimension;
        m_isReduced = false;
        m_tensorFull = p_degreeFunction.computeTensor(p_dimension, p_degree);
        m_numberOfFunctions = m_tensorFull.rows();
        m_tensorSparse = m_tensorFull.sparseView();
        m_center = Eigen::ArrayXd(p_dimension);
        m_center.setZero();
        m_scale = Eigen::ArrayXd(p_dimension);
        m_scale.setOnes();
    }

    /**
    * \brief Create a multi-variate basis
    * \param  p_numberOfVariates    number of variates
    * \param  p_numberOfFunctions   number of elements in the basis
    * \param  p_tensorFull          Full Basis representation with size p_numberOfFunctions x p_numberOfVariates
    * \param  p_tensorSparse        Sparse Basis representation with size p_numberOfFunctions x p_numberOfVariates
    * \param  p_isReduced           true if the basis is reduced
    * \param  p_center              center of the domain
    * \param  p_scale               inverse of the scaling factor to map the domain to [-1, 1]^nb_variates
    */
    MultiVariateBasis(const int   &p_numberOfVariates, const int &p_numberOfFunctions, const RowMatrixXi &p_tensorFull, const  RowSparseMatrix &p_tensorSparse,
                      const bool &p_isReduced, const  Eigen::ArrayXd &p_center, const  Eigen::ArrayXd   &p_scale): m_numberOfVariates(p_numberOfVariates),
        m_numberOfFunctions(p_numberOfFunctions), m_tensorFull(p_tensorFull), m_tensorSparse(p_tensorSparse), m_isReduced(p_isReduced), m_center(p_center), m_scale(p_scale)
    {  }
    /**

    * \brief Print a Basis to a stream
    *
    * \param p_os an output stream
    * \param p_basis the basis to be printed out
    */
    friend std::ostream &operator<<(std::ostream &p_os, const MultiVariateBasis &p_basis)
    {
        p_os << "Number of Variates " << p_basis.getNumberOfVariates() << std::endl;
        p_os << "Number of Functions " << p_basis.getNumberOfFunctions() << std::endl;
        p_os << "Tensor Full" << std::endl << p_basis.m_tensorFull << std::endl;
        p_os << "Tensor Sparse" << std::endl << p_basis.m_tensorSparse << std::endl;
        return p_os;
    }

    /**
     * \brief Return the number of functions
     */
    inline int getNumberOfFunctions() const
    {
        return m_numberOfFunctions;
    }

    /**
     * \brief Return the number of variates
     */
    inline int getNumberOfVariates() const
    {
        return m_numberOfVariates;
    }

    /**
     * \brief Return the full tensor representation
     */
    inline const RowMatrixXi &getTensorFull() const
    {
        return m_tensorFull;
    }

    /**
     * \brief Return the sparse tensor representation
     */
    inline const RowSparseMatrix &getTensorSparse() const
    {
        return m_tensorSparse;
    }

    /**
     * \brief Return if the base is reduced
     */
    inline const bool &getIsReduced() const
    {
        return m_isReduced;
    }

    /**
     * \brief Return the center
     */
    inline const Eigen::ArrayXd &getCenter() const
    {
        return m_center;
    }

    /**
     * \brief Return the scale
     */
    inline const Eigen::ArrayXd  &getScale() const
    {
        return m_scale;
    }


    /**
     * \brief Make the basis use a reduced domain using the mapping
     * \f$ p_x \in D \mapsto \left(\frac{x_i - center_i}{scale_i}\right)_i \f$
     *
     * \param p_center center of the domain
     * \param p_scale width of the domain
     */
    void setReduced(const Eigen::ArrayXd &p_center, const Eigen::ArrayXd &p_scale)
    {
        assert(p_center.rows() == p_scale.rows());
        assert(p_center.rows() == m_numberOfVariates);
        m_center = p_center;
        m_scale = 1. / p_scale;
        m_isReduced = true;
    }

    /**
     * \brief Make the basis use a reduced domain using the mapping
     * \f$ p_x \in D \mapsto \left(\frac{x_i - (a_i + b_i)/2}{(b_i - a_i)/2}\right)_i \f$
     * The domain [a, b] is mapped to [-1,1]^d
     *
     * \param p_lower lower bound a of the domain
     * \param p_upper upper bound b of the domain
     */
    void setDomain(const Eigen::ArrayXd &p_lower, const Eigen::ArrayXd &p_upper)
    {
        assert(p_lower.rows() == p_upper.rows());
        assert(p_lower.rows() == m_numberOfVariates);
        m_center = (p_lower + p_upper) / 2;
        m_scale = 2. / (p_upper - p_lower);
        m_isReduced = true;
    }

    /**
     * \brief Evalute the i-th basis function at the point p_x
     *
     * \param p_x evaluation point
     * \param p_i index of the function to evaluate between 0 and m_numberOfFunctions - 1
     */
    double operator()(const Eigen::ArrayXd &p_x, int p_i) const
    {
        assert(p_i >= 0 && p_i < getNumberOfFunctions());
        double aux = 1.;
        for (RowSparseMatrix::InnerIterator it(m_tensorSparse, p_i); it; ++it)
        {
            const int k = it.col();
            const int Tik = it.value();
            aux *= m_func1D.F((p_x(k) - m_center(k)) * m_scale(k), Tik);
        }
        return aux;
    }

    /**
     * \brief An element of a basis writes as a product
     *      p_1(x_1) p2(x_2) .... p_n(x_n)
     * for a polynomial with n variates, where each p_k is a polynomial
     * with only one variate.
     *
     * This functions evaluates the term p_k of the i-th element of the
     * basis at the point p_x
     *
     * \param p_x evaluation point
     * \param p_i index of the function to evaluate between 0 and m_numberOfFunctions - 1
     * \param p_k index of the variate between 0 and m_numberOfVariates - 1
     */
    double operator()(const Eigen::ArrayXd &p_x, int p_i, int p_k) const
    {
        assert(p_i >= 0 && p_i < getNumberOfFunctions());
        assert(p_k >= 0 && p_k < getNumberOfVariates());
        const int Tik = m_tensorFull(p_i, p_k);
        if (Tik == 0) return 1.;
        return m_func1D.F((p_x(p_k) - m_center(p_k)) * m_scale(p_k), Tik);
    }

    /**
      * \brief An element of a basis writes as a product
      *      p_1(x_1) p2(x_2) .... p_n(x_n)
      * for a polynomial with n variates, where each p_k is a polynomial with only one variate.
      *
      * This functions evaluates the term p_k of the i-th element of the basis at the point p_x
      *
      * \param p_x evaluation point
      * \param p_i index of the function to evaluate between 0 and m_numberOfFunctions - 1
      * \param p_k index of the variate between 0 and m_numberOfVariates - 1
      *
      * See #operator()(const Eigen::ArrayXd &p_x, int p_i, int p_k) const
      */
    inline double eval(const Eigen::ArrayXd &p_x, int p_i, int p_k) const
    {
        return operator()(p_x, p_i, p_k);
    }

    /**
      * \brief Evalute the i-th basis function at the point p_x
      *
      * \param p_x evaluation point
      * \param p_i index of the function to evaluate between 0 and m_numberOfFunctions - 1
      *
      * See #operator()(const Eigen::ArrayXd &p_x, int p_i) const
      */
    inline double eval(const Eigen::ArrayXd &p_x, int p_i) const
    {
        return operator()(p_x, p_i);
    }

    /**
     * \brief Evaluate the linear combination of the basis functions given by
     * p_basisCoordinates at the point p_x
     *
     * \param p_x evaluation point
     * \param p_basisCoordinates coefficients of the linear combination
     *
     */
    double eval(const Eigen::ArrayXd &p_x, const Eigen::ArrayXd &p_basisCoordinates) const
    {
        assert(p_basisCoordinates.size() == getNumberOfFunctions());
        assert(p_x.size() == getNumberOfVariates());
        double res = 0.;
        for (int i = 0; i < p_basisCoordinates.size(); i++)
        {
            const double alpha_i = p_basisCoordinates(i);
            res += operator()(p_x, i) * alpha_i;
        }
        return res;
    }

    /**
      * \brief Evaluate all the basis functions at the point p_x and return
      * the values in a vector. The returned vector has size NumberOfFunctions.
      *
      * \param p_x evaluation point
      */
    Eigen::ArrayXd eval(const Eigen::ArrayXd &p_x) const
    {
        assert(p_x.size() == getNumberOfVariates());
        Eigen::ArrayXd res(getNumberOfFunctions());
        for (int i = 0; i < getNumberOfFunctions(); i++)
        {
            res(i) = operator()(p_x, i);
        }
        return res;
    }

    /**
     * \brief Compute the first order partial derivative of the i-th function
     * with respect to the j-th variable
     *
     * \param p_x evaluation point
     * \param p_i index of the function between 0 and m_numberOfFunctions - 1
     * \param p_j index of the partial derivative between 0 and m_numberOfVariates - 1
     *
     * \return (D(f_i)/Dj)(p_x)
     */
    double evalD(const Eigen::ArrayXd &p_x, int p_i, int p_j) const
    {
        assert(p_x.size() == getNumberOfVariates());
        assert(p_i >= 0 && p_i < getNumberOfFunctions());
        assert(p_j >= 0 && p_j < getNumberOfVariates());

        if (m_tensorFull(p_i, p_j) == 0) return 0.;
        double aux = 1.;
        for (RowSparseMatrix::InnerIterator it(m_tensorSparse, p_i); it; ++it)
        {
            const int k = it.col();
            const int Tik = it.value();
            if (k == p_j)
                aux *= m_scale(k) * m_func1D.DF((p_x(k) - m_center(k)) * m_scale(k), Tik);
            else
                aux *= m_func1D.F((p_x(k) - m_center(k)) * m_scale(k), Tik);
        }
        return aux;
    }

    /**
     * \brief Compute the second order partial derivative of the i-th function
     * with respect to the j1-th and j2-th variables
     *
     * \param p_x evaluation point
     * \param p_i index of the function between 0 and m_numberOfFunctions - 1
     * \param p_j1 index of the first partial derivative between 0 and m_numberOfVariates - 1
     * \param p_j2 index of the second partial derivative between 0 and m_numberOfVariates - 1
     *
     * \return (D(f_i)/(Dj1 Dj2))(p_x)
     */
    double evalD2(const Eigen::ArrayXd &p_x, int p_i, int p_j1, int p_j2) const
    {
        assert(p_x.size() == getNumberOfVariates());
        assert(p_i >= 0 && p_i < getNumberOfFunctions());
        assert(p_j1 >= 0 && p_j1 < getNumberOfVariates());
        assert(p_j2 >= 0 && p_j2 < getNumberOfVariates());

        if ((m_tensorFull(p_i, p_j1) == 0) || (m_tensorFull(p_i, p_j2) == 0)) return 0.;
        double res = 1.;

        for (RowSparseMatrix::InnerIterator it(m_tensorSparse, p_i); it; ++it)
        {
            const int k = it.col();
            const int Tik = it.value();
            if ((p_j1 == p_j2) && (k == p_j1))
                res *= m_scale(k) * m_scale(k) * m_func1D.D2F((p_x(k) - m_center(k)) * m_scale(k), Tik);
            else if ((p_j1 != p_j2) && ((k == p_j1) || (k == p_j2)))
                res *= m_scale(k) * m_func1D.DF((p_x(k) - m_center(k)) * m_scale(k), Tik);
            else
                res *= m_func1D.F((p_x(k) - m_center(k)) * m_scale(k), Tik);
        }
        return res;
    }

    /**
     * \brief Compute the first order partial derivative with respect to the j-th variable of the
     * linear combination given by p_basisCoordinates.
     *
     * \param p_x evaluation point
     * \param p_basisCoordinates weights of the decomposition
     * \param p_j index of the partial derivative between 0 and m_numberOfVariates - 1
     *
     */
    double evalD(const Eigen::ArrayXd &p_x, const Eigen::ArrayXd &p_basisCoordinates, int p_j) const
    {
        assert(p_x.size() == getNumberOfVariates());
        assert(p_basisCoordinates.size() == getNumberOfFunctions());
        assert(p_j >= 0 && p_j < getNumberOfVariates());
        double res = 0.;
        for (int i = 0; i < p_basisCoordinates.size(); i++)
        {
            const double alpha_i = p_basisCoordinates(i);
            res += evalD(p_x, i, p_j) * alpha_i;
        }
        return res;
    }

    /**
     * \brief Compute the gradient of the linear combination given by p_basisCoordinates.
     *
     * \param p_x evaluation point
     * \param p_basisCoordinates weights of the decomposition
     * \param[out] p_dfx the gradient vector on output
     */
    void evalGradient(const Eigen::ArrayXd &p_x, const Eigen::ArrayXd &p_basisCoordinates, Eigen::ArrayXd &p_dfx) const
    {
        assert(p_x.size() == getNumberOfVariates());
        assert(p_basisCoordinates.size() == getNumberOfFunctions());
        p_dfx.resize(getNumberOfVariates());
        p_dfx.setZero();
        for (int i = 0; i < p_basisCoordinates.size(); i++)
        {
            const double alpha_i = p_basisCoordinates(i);
            for (int j = 0; j < getNumberOfVariates(); j++)
            {
                p_dfx(j) += evalD(p_x, i, j) * alpha_i;
            }
        }
    }


    /**
     * \brief Compute the second order partial derivative with respect to the j1-th and j2-th
     * variables of the linear combination given by p_basisCoordinates.
     *
     * \param p_x evaluation point
     * \param p_basisCoordinates weights of the decomposition
     * \param p_j1 index of the first partial derivative between 0 and m_numberOfVariates - 1
     * \param p_j2 index of the second partial derivative between 0 and m_numberOfVariates - 1
     *
     */
    double evalD2(const Eigen::ArrayXd &p_x, const Eigen::ArrayXd &p_basisCoordinates, int p_j1, int p_j2) const
    {
        assert(p_basisCoordinates.size() == getNumberOfFunctions());
        double res = 0.;
        for (int i = 0; i < p_basisCoordinates.size(); i++)
        {
            const double alpha_i = p_basisCoordinates(i);
            res += evalD2(p_x, i, p_j1, p_j2) * alpha_i;
        }
        return res;
    }

    /**
     * \brief Compute the gradient of the linear combination given by p_basisCoordinates.
     *
     * \param p_x evaluation point
     * \param p_basisCoordinates weights of the decomposition
     * \param[out] p_d2fx the Hessian matrix on output
     */
    void evalHessian(const Eigen::ArrayXd &p_x, const Eigen::ArrayXd &p_basisCoordinates, Eigen::ArrayXXd &p_d2fx) const
    {
        assert(p_basisCoordinates.size() == getNumberOfFunctions());
        p_d2fx.resize(getNumberOfVariates(), getNumberOfVariates());
        p_d2fx.setZero();
        for (int i = 0; i < p_basisCoordinates.size(); i++)
        {
            const double alpha_i = p_basisCoordinates(i);
            for (int j = 0; j < getNumberOfVariates(); j++)
            {
                for (int l = 0; l <= j; l++)
                {
                    p_d2fx(j, l) += evalD2(p_x, i, j, l) * alpha_i;
                    p_d2fx(l, j) = p_d2fx(j, l);
                }
            }
        }
        // Make d2fx symmetric by copying the lower triangle into the upper one.
        p_d2fx.matrix().triangularView<Eigen::StrictlyUpper>() = p_d2fx.matrix().triangularView<Eigen::StrictlyLower>().transpose();
    }

    /**
     * \brief Compute the gradient and the Hessian matrix of the linear
     * combination given by p_basisCoordinates
     *
     * \param p_x the evaluation point
     * \param p_basisCoordinates the coefficients of the linear combination
     * \param[out] p_fx contains the function value on output
     * \param[out] p_dfx the gradient vector on output
     * \param[out] p_d2fx the Hessian matrix on output
     */
    void evalDerivatives(const Eigen::ArrayXd &p_x, const Eigen::ArrayXd &p_basisCoordinates, double &p_fx, Eigen::ArrayXd &p_dfx, Eigen::ArrayXXd &p_d2fx) const
    {
        assert(p_basisCoordinates.size() == getNumberOfFunctions());

        int nVariates = getNumberOfVariates();
        p_fx = 0.;
        p_dfx.resize(nVariates);
        p_dfx.setZero();
        p_d2fx.resize(nVariates, nVariates);
        p_d2fx.setZero();
        Eigen::ArrayXd f(nVariates);
        Eigen::ArrayXd Df(nVariates);

        for (int i = 0; i < p_basisCoordinates.size(); i++)
        {
            double auxf = 1.;
            const double alpha_i = p_basisCoordinates(i);
            f.setZero();
            Df.setZero();
            for (RowSparseMatrix::InnerIterator it(m_tensorSparse, i); it; ++it)
            {
                const int k = it.col();
                const int Tik = it.value();
                f(k) = m_func1D.F((p_x(k) - m_center(k)) * m_scale(k), Tik);
                auxf *= f(k);
            }
            p_fx += alpha_i * auxf;

            for (int j = 0; j < nVariates; j++)
            {
                auxf = 1.;
                for (RowSparseMatrix::InnerIterator it(m_tensorSparse, i); it; ++it)
                {
                    const int k = it.col();
                    if (k != j) auxf *= f(k);
                }
                const int Tij = m_tensorFull(i, j);
                // Gradient
                Df(j) = m_scale(j) * m_func1D.DF((p_x(j) - m_center(j)) * m_scale(j), Tij);
                p_dfx(j) += alpha_i * Df(j) * auxf;
                // Diagonal terms of the Hessian matrix
                p_d2fx(j, j) += alpha_i * m_scale(j) * m_scale(j) * m_func1D.D2F((p_x(j) - m_center(j)) * m_scale(j), Tij) * auxf;

                for (int l = 0; l < j ; l++)
                {
                    auxf = 1.;
                    for (RowSparseMatrix::InnerIterator it(m_tensorSparse, i); it; ++it)
                    {
                        const int k = it.col();
                        if ((k != j) && (k != l)) auxf *= f(k);
                    }
                    p_d2fx(j, l) += alpha_i * Df(j) * Df(l) * auxf;
                }
            }
        }
        p_d2fx.matrix().triangularView<Eigen::StrictlyUpper>() = p_d2fx.matrix().triangularView<Eigen::StrictlyLower>().transpose();
    }


    /**
     * \brief Find the best least square approximation of the function defined
     * by f(p_x(i,:)) = fx(i) on the MultiViariateBasis
     *
     * \param[out] p_basisCoordinates contains on exit the coefficients of the regression
     * \param p_x the matrix of points at which we know the value of the
     * function. One column of the matrix corresponds to the coordinates of
     * one point. The matrix p_x has size dim p_x NumberOfPoints
     * \param p_fx the values of the function f at the points defined by p_x
     */
    void fitLeastSquare(Eigen::ArrayXd &p_basisCoordinates, const Eigen::ArrayXXd &p_x, const Eigen::ArrayXd &p_fx) const
    {
        assert(p_x.cols() == p_fx.size());
        assert(p_x.rows() == getNumberOfVariates());
        int nFunctions = getNumberOfFunctions();
        p_basisCoordinates.resize(nFunctions);
        p_basisCoordinates.setZero();

#ifdef _OPENMP
        int nbThread = omp_get_max_threads();
#else
        int nbThread = 1 ;
#endif
        Eigen::MatrixXd phi_k(nFunctions, nbThread);
        Eigen::MatrixXd b(nFunctions, nbThread);
        int i;
#ifdef _OPENMP
        #pragma omp parallel for  private(i)
#endif
        for (i = 0 ; i < nbThread; ++i)
        {
            phi_k.col(i).setZero();
            b.col(i).setZero();
        }
        std::vector<Eigen::MatrixXd>  A(nbThread);
        for (int it = 0; it < nbThread; ++it)
            A[it].resize(nFunctions, nFunctions);
#ifdef _OPENMP
        #pragma omp parallel for  private(i)
#endif
        for (i = 0 ; i < nbThread; ++i)
            A[i].setZero();

        /* construct A and b*/
#ifdef _OPENMP
        #pragma omp parallel for  private(i)
#endif
        for (i = 0; i < p_fx.size(); i++)
        {
#ifdef _OPENMP
            int ithread = omp_get_thread_num();
#else
            int ithread = 0;
#endif
            for (int k = 0; k < nFunctions; k++)
            {
                const double tmp = eval(p_x.col(i), k);
                b(k, ithread) += tmp * p_fx(i);
                phi_k(k, ithread) = tmp;
            }
            /* A += phi_k' * phi_k */
            A[ithread] += phi_k.col(ithread) * phi_k.col(ithread).transpose();
        }
        // collapse AThread
        for (int it = 1 ;  it < nbThread; ++it)
        {
            A[0] += A[it];
            b.col(0) += b.col(it);
        }
        /* Because A often comes from simulation, A is not >0. So we use a
         * least-square approach
         */
        p_basisCoordinates = A[0].colPivHouseholderQr().solve(b.col(0));
    }

};

/**@}*/
}
#endif /* _MULTIVARIATEBASIS_H */