gp.h

Go to the documentation of this file.

#ifndef GP_H
#define GP_H
 
 
#include <distribution.h>
#include <scaler.h>
#include <covariance.h>
#include <mean++.h>
#include <prior++.h>
#include <optimization.h>
#include <memory>
 
namespace cmp::gp {
 
enum method {
    MLE,     
    LOO,     
    LOO_MSE  
};
 
enum type {
    PRIOR,     
    POSTERIOR  
};
 
class GaussianProcess {
  private:
 
    // Hyperparameters and prior mean and covariance functions
    Eigen::VectorXd par_;                             
    std::shared_ptr<covariance::Covariance> pKernel_; 
    std::shared_ptr<mean::Mean> pMean_;               
    double nugget_;                                   
 
    // Observations (the GP can own the data or point to external data, if it owns the data the pointer will point to it)
    std::optional<Eigen::MatrixXd> xObs_;                      
    std::optional<Eigen::VectorXd> yObs_;                      
    std::optional<Eigen::Ref<const Eigen::MatrixXd>> pXObs_;   
    std::optional<Eigen::Ref<const Eigen::VectorXd>> pYObs_;   
 
    // Internal storage for the covariance matrix, its decomposition and related vectors (for fast query)
    Eigen::LDLT<Eigen::MatrixXd> covDecomposition_;            
    Eigen::VectorXd alpha_;                                    
    Eigen::VectorXd diagCovInverse_;                           
    Eigen::VectorXd residual_;                                 
 
    // Internal flag to check whether we normalize y or not
    bool normalizeY_ = false;               
    cmp::scaler::StandardScaler yScaler_;    
 
// Private members
  private:
    void compute(const Eigen::Ref<const Eigen::VectorXd> &par);
 
  public:
 
    // Default constructor and destructor
    GaussianProcess();
    ~GaussianProcess() = default;
 
    // Constructor with parameters
    GaussianProcess(const std::shared_ptr<covariance::Covariance> &kernel, const std::shared_ptr<mean::Mean> &mean, Eigen::Ref<const Eigen::VectorXd> params, double nugget = 1e-8);
 
    // Copy constructors and operators
    GaussianProcess(const GaussianProcess &other);
    GaussianProcess &operator=(const GaussianProcess &other);
    GaussianProcess(GaussianProcess &&other) noexcept;
    GaussianProcess &operator=(GaussianProcess &&other) noexcept;
 
    void set(const std::shared_ptr<covariance::Covariance> &kernel, const std::shared_ptr<mean::Mean> &mean, Eigen::Ref<const Eigen::VectorXd> params, double nugget = 1e-8);
 
    void condition(const Eigen::Ref<const Eigen::MatrixXd> &xObs, const Eigen::Ref<const Eigen::VectorXd> &yObs, bool copyData = true, bool normalizeY = false);
 
    void fit(const Eigen::Ref<const Eigen::MatrixXd> &xObs, const Eigen::Ref<const Eigen::VectorXd> &yObs, const Eigen::Ref<const Eigen::VectorXd> &lb, const Eigen::Ref<const Eigen::VectorXd> &ub, const method &method, const nlopt::algorithm &alg, const double &tol_rel, bool copyData = true, bool normalizeY = false, const std::shared_ptr<cmp::prior::Prior> &prior = cmp::prior::Uniform::make(), const std::vector<bool> &logScale = {});
 
    Eigen::VectorXd getParameters() const {
        return par_;
    }
 
    std::shared_ptr<cmp::mean::Mean> getMean() const {
        return this->pMean_;
    }
 
    std::shared_ptr<cmp::covariance::Covariance> getKernel() const {
        return this->pKernel_;
    }
 
    double getNugget() const {
        return this->nugget_;
    }
 
    Eigen::MatrixXd covariance(Eigen::Ref<const Eigen::VectorXd> par) const;
 
    Eigen::MatrixXd covarianceGradient(Eigen::Ref<const Eigen::VectorXd> par, const int &i) const;
 
    Eigen::MatrixXd covarianceHessian(Eigen::Ref<const Eigen::VectorXd> par, const size_t &i, const size_t &j) const;
 
    const Eigen::LDLT<Eigen::MatrixXd> &getCovDecomposition() const {
        return covDecomposition_;
    }
 
    const Eigen::VectorXd &getAlpha() const {
        return alpha_;
    }
 
    const Eigen::VectorXd &getDiagCovInverse() const {
        return diagCovInverse_;
    }
 
    const Eigen::VectorXd &getResidualVector() const {
        return residual_;
    }
 
    const Eigen::Ref<const Eigen::MatrixXd> &getXObs() const {
        return pXObs_.value();
    }
 
    const Eigen::Ref<const Eigen::VectorXd> &getYObs() const {
        return pYObs_.value();
    }
 
    size_t nObs() const {
        if(!pXObs_.has_value()) {
            return 0;
        }
        return static_cast<size_t>(pXObs_->rows());
    }
    Eigen::VectorXd priorMean(Eigen::Ref<const Eigen::VectorXd> par) const;
 
    Eigen::VectorXd priorMeanGradient(Eigen::Ref<const Eigen::VectorXd> par, const int &i) const;
 
    Eigen::VectorXd residual(Eigen::Ref<const Eigen::VectorXd> par) const;
 
    std::pair<double, double> predict(const Eigen::Ref<const Eigen::VectorXd> &x, type predictionType = type::POSTERIOR) const;
 
    double predictMean(const Eigen::Ref<const Eigen::VectorXd> &x, type predictionType = type::POSTERIOR) const;
 
    std::pair<Eigen::VectorXd, Eigen::MatrixXd> predictMultiple(const Eigen::Ref<const Eigen::MatrixXd> &x_pts, type predictionType = type::POSTERIOR) const;
 
    Eigen::VectorXd predictMeanMultiple(const Eigen::Ref<const Eigen::MatrixXd> &x_pts, type predictionType = type::POSTERIOR) const;
 
 
    std::pair<double, double> predictLOO(const size_t &i) const;
 
    double logLikelihood() const;
 
    double logLikelihoodLOO(const size_t &i) const;
 
    Eigen::MatrixXd expectedVarianceImprovement(const Eigen::Ref<const Eigen::MatrixXd> &x_pts,
                                                const Eigen::Ref<const Eigen::MatrixXd> &x_pending,
                                                double nu) const;
 
    double objectiveFunction(const Eigen::Ref<const Eigen::VectorXd> &x, Eigen::Ref<Eigen::VectorXd> grad, const std::shared_ptr<cmp::prior::Prior> &prior);
 
    double objectiveFunctionLOO(const Eigen::Ref<const Eigen::VectorXd> &x, Eigen::Ref<Eigen::VectorXd> grad, const std::shared_ptr<cmp::prior::Prior> &prior);
 
    double objectiveFunctionLOOMSE(const Eigen::Ref<const Eigen::VectorXd> &x, Eigen::Ref<Eigen::VectorXd> grad, const std::shared_ptr<cmp::prior::Prior> &prior);
};
}
#endif // MACRO