// Copyright (C) 2016, 2017 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#ifndef BASEREGRESSION_H
#define BASEREGRESSION_H
#include <memory>
#include <vector>
#include <iostream>
#include <Eigen/Dense>
#include <Eigen/SVD>
#include "StOpt/core/grids/InterpolatorSpectral.h"

/** \file BaseRegression.h
 *  \brief Base class to define regressor for stochastic optimization by Monte Carlo
 *  \author Xavier Warin
 */
namespace StOpt
{
/// \class BaseRegression BaseRegression.h
/// Base class for regression
class BaseRegression
{
protected :

    bool m_bZeroDate ;                    ///< Is the regression date zero ?
    bool m_bRotationAndRescale ; ///< do we rescale particles and do a rotation with SVD on data
    Eigen::ArrayXd  m_meanX ; ///< store  scaled factor in each direction (average of particles values in each direction)
    Eigen::ArrayXd  m_etypX ; ///< store  scaled factor in each direction (standard deviation of particles in each direction)
    Eigen::MatrixXd m_svdMatrix ; ///< svd matrix transposed  used to transform particles
    Eigen::ArrayXd  m_sing      ; ///< singular values associated to SVD
    Eigen::ArrayXXd m_particles;  ///< Particles used to regress: first dimension  : dimension of the problem , second dimension : the  number of particles. These particles are rescaled and a rotation with SVD is achieved to avoid degeneracy in case of high correlations

    // rotation for data and rescaling
    void preProcessData();

public :

    /// \brief Default constructor
    BaseRegression();

    /// \brief Default destructor
    virtual ~BaseRegression() {}

    /// \brief Default constructor
    BaseRegression(const bool &p_bRotationAndRescale);

    /// \brief Constructor storing the particles
    /// \param  p_bZeroDate        first date is 0?
    /// \param  p_particles        particles used for the meshes.
    ///                            First dimension  : dimension of the problem,
    ///                            second dimension : the  number of particles
    /// \param  p_bRotationAndRescale do we rescale particle
    BaseRegression(const bool &p_bZeroDate, const Eigen::ArrayXXd &p_particles, const bool &p_bRotationAndRescale);

    /// \brief Constructor used in simulation, no rotation
    /// \param  p_bZeroDate    first date is 0?
    /// \param  p_bRotationAndRescale do we rescale particle
    BaseRegression(const bool &p_bZeroDate, const bool &p_bRotationAndRescale);


    /// \brief Last constructor used in simulation
    /// \param  p_bZeroDate    first date is 0?
    /// \param  p_meanX            scaled factor in each direction (average of particles values in each direction)
    /// \param  p_etypX            scaled factor in each direction (standard deviation of particles in each direction)
    /// \param  p_svdMatrix        svd matrix transposed  used to transform particles
    /// \param  p_bRotationAndRescale do we rescale particle

    BaseRegression(const bool &p_bZeroDate, const   Eigen::ArrayXd &p_meanX, const   Eigen::ArrayXd   &p_etypX, const   Eigen::MatrixXd   &p_svdMatrix, const bool &p_bRotationAndRescale);

    /// \brief Copy constructor
    /// \param p_object  object to copy
    BaseRegression(const BaseRegression &p_object);

    /// \brief update the particles used in regression  and construct the matrices
    /// \param  p_bZeroDate    first date is 0?
    /// \param  p_particles    particles used for the meshes.
    ///                        Firs dimension  : dimension of the problem,
    ///                        second dimension : the  number of particles
    void updateSimulationsBase(const bool &p_bZeroDate, const Eigen::ArrayXXd &p_particles);

    /// \brief Get some local accessors
    ///@{
    virtual inline Eigen::ArrayXXd  getParticles() const
    {
        return m_particles ;
    }

    /// \brief Get bRotationAndRescale
    virtual inline bool getBRotationAndRescale() const
    {
        return m_bRotationAndRescale ;
    }

    /// \brief Get average of simulation per dimension
    virtual inline Eigen::ArrayXd getMeanX() const
    {
        return m_meanX;
    }

    /// \brief get standard deviation per dimension
    virtual inline Eigen::ArrayXd getEtypX() const
    {
        return m_etypX;
    }

    /// \brief get back the SVD matrix used for rescaling particles
    virtual inline Eigen::MatrixXd getSvdMatrix() const
    {
        return m_svdMatrix;
    }

    /// \brief get back singular values
    virtual inline Eigen::ArrayXd getSing() const
    {
        return m_sing;
    }

    /// \brief Get dimension of the problem
    virtual inline int getDimension() const
    {
        return m_particles.rows();
    }

    /// \brief Get the number of simulations
    virtual inline int  getNbSimul()const
    {
        return m_particles.cols() ;
    }

    /// \brief get back particle by its number
    /// \param p_iPart   particle number
    /// \return the particle (if no particle, send back an empty array)
    virtual Eigen::ArrayXd getParticle(const int &p_iPart) const;

    /// \brief get the number of basis functions
    virtual  int getNumberOfFunction() const = 0 ;

    ///@}
    /// \brief Constructor storing the particles
    /// \brief update the particles used in regression  and construct the matrices
    /// \param  p_bZeroDate    first date is 0?
    /// \param  p_particles    particles used for the meshes.
    ///                        First dimension  : dimension of the problem,
    ///                        second dimension : the  number of particles
    virtual void updateSimulations(const bool &p_bZeroDate, const Eigen::ArrayXXd  &p_particles) = 0 ;

    /// \brief conditional expectation basis function coefficient calculation
    /// \param  p_fToRegress  function to regress associated to each simulation used in optimization
    /// \return regression coordinates on the basis  (size : number of meshes multiplied by the dimension plus one)
    /// @{
    virtual Eigen::ArrayXd getCoordBasisFunction(const Eigen::ArrayXd &p_fToRegress) const = 0;
    ///@}
    /// \brief conditional expectation basis function coefficient calculation for multiple functions to regress
    /// \param  p_fToRegress  function to regress associated to each simulation used in optimization (size : number of functions to regress \times the number of Monte Carlo simulations)
    /// \return regression coordinates on the basis  (size :  number of function to regress  \times number of meshes multiplied by the dimension plus one)
    /// @{
    virtual Eigen::ArrayXXd getCoordBasisFunctionMultiple(const Eigen::ArrayXXd &p_fToRegress) const = 0 ;
    ///@}

    /// \brief conditional expectation calculation
    /// \param  p_fToRegress  simulations  to regress used in optimization
    /// \return regressed value function
    /// @{
    virtual Eigen::ArrayXd getAllSimulations(const Eigen::ArrayXd &p_fToRegress) const = 0;
    virtual Eigen::ArrayXXd getAllSimulationsMultiple(const Eigen::ArrayXXd &p_fToRegress) const = 0;
    ///@}

    /// \brief Use basis functions to reconstruct the solution
    /// \param p_basisCoefficients basis coefficients
    ///@{
    virtual Eigen::ArrayXd reconstruction(const Eigen::ArrayXd   &p_basisCoefficients) const = 0 ;
    virtual Eigen::ArrayXXd reconstructionMultiple(const Eigen::ArrayXXd   &p_basisCoefficients) const = 0;
    /// @}

    /// \brief use basis function to reconstruct a given simulation
    /// \param p_isim               simulation number
    /// \param p_basisCoefficients  basis coefficients to reconstruct a given conditional expectation
    virtual double reconstructionASim(const int &p_isim, const Eigen::ArrayXd   &p_basisCoefficients) const = 0 ;

    /// \brief conditional expectation reconstruction
    /// \param  p_coordinates        coordinates to interpolate (uncertainty sample)
    /// \param  p_coordBasisFunction regression coordinates on the basis  (size: number of meshes multiplied by the dimension plus one)
    /// \return regressed value function reconstructed for each simulation
    virtual double getValue(const Eigen::ArrayXd   &p_coordinates,
                            const Eigen::ArrayXd   &p_coordBasisFunction) const = 0;

    /// \brief conditional expectation reconstruction for a lot of simulations
    /// \param  p_coordinates        coordinates to interpolate (uncertainty sample) size uncertainty dimension by number of samples
    /// \param  p_coordBasisFunction regression coordinates on the basis  (size: number of meshes multiplied by the dimension plus one)
    /// \return regressed value function reconstructed for each simulation
    Eigen::ArrayXd  getValues(const Eigen::ArrayXXd   &p_coordinates,
                              const Eigen::ArrayXd   &p_coordBasisFunction) const
    {
        Eigen::ArrayXd valRet(p_coordinates.cols());
        for (int is  = 0; is < p_coordinates.cols(); ++is)
            valRet(is) = getValue(p_coordinates.col(is), p_coordBasisFunction);
        return valRet;
    }

    /// \brief permits to reconstruct a function with basis functions coefficients values given on a grid
    /// \param  p_coordinates          coordinates  (uncertainty sample)
    /// \param  p_ptOfStock            grid point
    /// \param  p_interpFuncBasis      spectral interpolator to interpolate the basis functions  coefficients used in regression on the grid (given for each basis function)
    virtual double getAValue(const Eigen::ArrayXd &p_coordinates,  const Eigen::ArrayXd &p_ptOfStock,
                             const std::vector< std::shared_ptr<InterpolatorSpectral> > &p_interpFuncBasis) const = 0;

    /// \brief is the regression date zero
    inline bool getBZeroDate() const
    {
        return m_bZeroDate;
    }

    /// \brief Clone the regressor
    virtual std::shared_ptr<BaseRegression> clone() const = 0 ;


};

}

#endif