Commits

Ruben Martinez-Cantin  committed b81889c

Adding full Bayes model

  • Participants
  • Parent commits ad0c283

Comments (0)

Files changed (12)

File CMakeLists.txt

   ./src/gaussian_process.cpp
   ./src/gaussian_process_ml.cpp
   ./src/student_t_process_jef.cpp
+  ./src/student_t_process_nig.cpp
 #  ./src/gaussian_process_ign.cpp
 #  ./src/student_t_process.cpp
   ./src/parameters.cpp

File app/bo_display.cpp

   parameters.n_init_samples = 7;
   parameters.n_iter_relearn = 0;
   parameters.n_iterations = 150;
-  parameters.surr_name = S_STUDENT_T_PROCESS_JEFFREYS;
+  parameters.surr_name = S_STUDENT_T_PROCESS_NORMAL_INV_GAMMA;
   parameters.kernel.hp_mean[0] = 1;
   parameters.kernel.hp_std[0] = 0.1;
   parameters.kernel.n_hp = 1;
   // parameters.mean.name = "mZero";
   //  parameters.crit_name = "cHedge(cEI,cLCB,cExpReturn,cOptimisticSampling)";
   // parameters.epsilon = 0.0;
-  parameters.verbose_level = 0;
+  parameters.verbose_level = 2;
 
   state_ii = 0;
 

File include/parameters.h

   typedef enum {  
     S_GAUSSIAN_PROCESS,
     S_GAUSSIAN_PROCESS_ML,
-    S_GAUSSIAN_PROCESS_INV_GAMMA_NORMAL,
     S_STUDENT_T_PROCESS_JEFFREYS,
+    S_STUDENT_T_PROCESS_NORMAL_INV_GAMMA,
     S_ERROR = -1
   } surrogate_name;
 
   const double KERNEL_THETA    = 1.0;
   const double KERNEL_SIGMA    = 10.0;
   const double MEAN_MU         = 1.0;
-  const double MEAN_SIGMA      = 30.0;
+  const double MEAN_SIGMA      = 1000.0;
   const double PRIOR_ALPHA     = 1.0;
   const double PRIOR_BETA      = 1.0;
   const double DEFAULT_SIGMA   = 1.0;

File include/specialtypes.hpp

 typedef boost::numeric::ublas::zero_vector<double>             zvectord;
 typedef boost::numeric::ublas::scalar_vector<double>           svectord;
 typedef boost::numeric::ublas::matrix<double>                   matrixd;
+typedef boost::numeric::ublas::zero_matrix<double>             zmatrixd;
 
 typedef std::vector<vectord>                                   vecOfvec;
 typedef std::vector<vectord>::iterator                 vecOfvecIterator;

File include/student_t_process_ign.hpp

-/** \file student_t_process_ign.hpp
-    \brief Student's t process with Inverse-Gamma-Normal hyperprior 
-    on mean and signal variance parameters. */
-/*
--------------------------------------------------------------------------
-   This file is part of BayesOpt, an efficient C++ library for 
-   Bayesian optimization.
-
-   Copyright (C) 2011-2013 Ruben Martinez-Cantin <rmcantin@unizar.es>
- 
-   BayesOpt is free software: you can redistribute it and/or modify it 
-   under the terms of the GNU General Public License as published by
-   the Free Software Foundation, either version 3 of the License, or
-   (at your option) any later version.
-
-   BayesOpt is distributed in the hope that it will be useful, but 
-   WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-   GNU General Public License for more details.
-
-   You should have received a copy of the GNU General Public License
-   along with BayesOpt.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------------
-*/
-
-
-#ifndef  _STUDENT_T_PROCESS_IGN_HPP_
-#define  _STUDENT_T_PROCESS_IGN_HPP_
-
-#include "nonparametricprocess.hpp"
- 
-/** \addtogroup  NonParametricProcesses */
-/**@{*/
-
-
-/**
- * \brief Student's t process with Jeffreys hyperprior 
- *        on mean and signal variance parameters.
- */
-class StudentTProcessIGN: public NonParametricProcess
-{
-public:
-  StudentTProcessIGN(double noise, double alpha, double beta, double delta);
-  virtual ~StudentTProcessIGN();
-
-  /** 
-   * \brief Function that returns the prediction of the GP for a query point
-   * in the hypercube [0,1].
-   * 
-   * @param query in the hypercube [0,1] to evaluate the Gaussian process
-   * @return pointer to the probability distribution.
-   */	
-  ProbabilityDistribution* prediction(const vectord &query);
-
-
-private:
-
- /** 
-   * \brief Computes the negative log likelihood and its gradient of the data.
-   * 
-   *
-   * @return value negative log likelihood
-   */
-  double negativeLogLikelihood();
-
-  /** 
-   * \brief Precompute some values of the prediction that do not depends on
-   * the query
-   * @return error code
-   */
-  int precomputePrediction();
-
-private:
-  double mMu, mSig;                   //!< GP posterior parameters
-  double mAlpha, mBeta, mDelta2;
-  //! Precomputed GP prediction operations
-  vectord mUInvR;     
-  vectord mInvRy;         
-  double mUInvRUDelta;
-};
-
-
-/**@}*/
-
-
-#endif

File include/student_t_process_nig.hpp

+/** \file student_t_process_nig.hpp
+    \brief Student's t process with Normal-Inverse-Gamma hyperprior 
+    on mean and signal variance parameters. */
+/*
+-------------------------------------------------------------------------
+   This file is part of BayesOpt, an efficient C++ library for 
+   Bayesian optimization.
+
+   Copyright (C) 2011-2013 Ruben Martinez-Cantin <rmcantin@unizar.es>
+ 
+   BayesOpt is free software: you can redistribute it and/or modify it 
+   under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   BayesOpt is distributed in the hope that it will be useful, but 
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with BayesOpt.  If not, see <http://www.gnu.org/licenses/>.
+------------------------------------------------------------------------
+*/
+
+
+#ifndef  _STUDENT_T_PROCESS_NIG_HPP_
+#define  _STUDENT_T_PROCESS_NIG_HPP_
+
+#include "student_t_distribution.hpp"
+#include "hierarchical_gaussian_process.hpp"
+
+namespace bayesopt
+{
+
+  /** \addtogroup  NonParametricProcesses */
+  /**@{*/
+
+
+  /**
+   * \brief Student's t process with Normal Inverse-Gamma hyperprior 
+   *        on mean and snigal variance parameters.
+   */
+  class StudentTProcessNIG: public HierarchicalGaussianProcess
+  {
+  public:
+    StudentTProcessNIG(size_t dim, bopt_params params);
+    virtual ~StudentTProcessNIG();
+
+    /** 
+     * \brief Function that returns the prediction of the GP for a query point
+     * in the hypercube [0,1].
+     * 
+     * @param query in the hypercube [0,1] to evaluate the Gaussian process
+     * @return pointer to the probability distribution.
+     */	
+    ProbabilityDistribution* prediction(const vectord &query);
+
+
+  private:
+
+    /** 
+     * \brief Computes the negative log likelihood and its gradient of the data.
+     * @return value negative log likelihood
+     */
+    double negativeLogLikelihood();
+
+    /** 
+     * \brief Precompute some values of the prediction that do not depends on
+     * the query
+     * @return error code
+     */
+    int precomputePrediction();
+
+  private:
+    vectord mWMap;                      //!< GP posterior parameters
+    double mSigmaMap;                   //!< GP posterior parameters
+    double mAlpha, mBeta;
+    vectord mW0;
+    vectord mInvVarW;
+    //! Precomputed GP prediction operations
+    vectord mVf;
+    matrixd mKF, mD;     
+
+    StudentTDistribution* d_;      //!< Predictive distributions
+  };
+
+
+  /**@}*/
+
+} //namespace bayesopt
+
+
+#endif

File src/nonparametricprocess.cpp

 #include "gaussian_process_ml.hpp"
 //#include "gaussian_process_ign.hpp"
 #include "student_t_process_jef.hpp"
+#include "student_t_process_nig.hpp"
 
 
 namespace bayesopt
       case S_STUDENT_T_PROCESS_JEFFREYS: 
       	s_ptr = new StudentTProcessJeffreys(dim,parameters); break;
 
+      case S_STUDENT_T_PROCESS_NORMAL_INV_GAMMA: 
+      	s_ptr = new StudentTProcessNIG(dim,parameters); break;
+
       default:
 	FILE_LOG(logERROR) << "Error: surrogate function not supported.";
 	return NULL;

File src/parameters.cpp

     return S_GAUSSIAN_PROCESS;
   if      (!strcmp(name,  "GAUSSIAN_PROCESS_ML"))
     return S_GAUSSIAN_PROCESS_ML;
-  else if (!strcmp(name,  "GAUSSIAN_PROCESS_INV_GAMMA_NORMAL"))
-    return S_GAUSSIAN_PROCESS_INV_GAMMA_NORMAL;
+  else if (!strcmp(name,  "STUDENT_T_PROCESS_NORMAL_INV_GAMMA"))
+    return S_STUDENT_T_PROCESS_NORMAL_INV_GAMMA;
   else if (!strcmp(name,  "STUDENT_T_PROCESS_JEFFREYS"))
     return S_STUDENT_T_PROCESS_JEFFREYS;
   else return S_ERROR;
     {
     case S_GAUSSIAN_PROCESS: return "GAUSSIAN_PROCESS"; 
     case S_GAUSSIAN_PROCESS_ML: return "GAUSSIAN_PROCESS_ML"; 
-    case S_GAUSSIAN_PROCESS_INV_GAMMA_NORMAL: return "GAUSSIAN_PROCESS_INV_GAMMA_NORMAL"; 
-    case S_STUDENT_T_PROCESS_JEFFREYS: return "S_STUDENT_T_PROCESS_JEFFREYS"; 
+    case S_STUDENT_T_PROCESS_NORMAL_INV_GAMMA: return "STUDENT_T_PROCESS_NORMAL_INV_GAMMA";
+    case S_STUDENT_T_PROCESS_JEFFREYS: return "STUDENT_T_PROCESS_JEFFREYS"; 
     case S_ERROR:
     default: return "ERROR!";
     }

File src/student_t_process_ign.cpp

-/*
------------------------------------------------------------------------------
-   Copyright (C) 2011 Ruben Martinez-Cantin <rmcantin@unizar.es>
- 
-   This program is free software: you can redistribute it and/or modify
-   it under the terms of the GNU General Public License as published by
-   the Free Software Foundation, either version 3 of the License, or
-   (at your option) any later version.
-
-   This program is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-   GNU General Public License for more details.
-
-   You should have received a copy of the GNU General Public License
-   along with this program.  If not, see <http://www.gnu.org/licenses/>.
------------------------------------------------------------------------------
-*/
-
-#include "student_t_process.hpp"
-#include "cholesky.hpp"
-#include "trace_ublas.hpp"
-#include "student_t_distribution.hpp"
-
-using boost::numeric::ublas::inplace_solve;
-using boost::numeric::ublas::lower_tag;
-using boost::numeric::ublas::lower;
-
-  
-StudentTProcessIGN::StudentTProcessIGN(double noise, double alpha, 
-					double beta, double delta):
-  NonParametricProcess(noise),
-  mAlpha(alpha), mBeta (beta), mDelta2(delta)
-{}  // Constructor
-
-
-
-StudentTProcessIGN::~StudentTProcessIGN()
-{} // Default destructor
-
-
-double StudentTProcessIGN::negativeLogLikelihoodIGN()
-{
-  
-  vectord h = mMean->getFeatureVector()

File src/student_t_process_jef.cpp

     mSigma = inner_prod(mAlphaF,mAlphaF)/(n-p);
     
     d_->setDof(n-p);  
-    return 1;
+    return 0;
   }
 
 } //namespace bayesopt

File src/student_t_process_nig.cpp

+/*
+-------------------------------------------------------------------------
+   This file is part of BayesOpt, an efficient C++ library for 
+   Bayesian optimization.
+
+   Copyright (C) 2011-2013 Ruben Martinez-Cantin <rmcantin@unizar.es>
+ 
+   BayesOpt is free software: you can redistribute it and/or modify it 
+   under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   BayesOpt is distributed in the hope that it will be useful, but 
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with BayesOpt.  If not, see <http://www.gnu.org/licenses/>.
+------------------------------------------------------------------------
+*/
+
+#include <boost/numeric/ublas/banded.hpp>
+#include "log.hpp"
+#include "student_t_process_nig.hpp"
+#include "cholesky.hpp"
+#include "trace_ublas.hpp"
+#include "elementwise_ublas.hpp"
+#include "student_t_distribution.hpp"
+
+namespace bayesopt
+{
+
+  namespace ublas = boost::numeric::ublas; 
+  
+  StudentTProcessNIG::StudentTProcessNIG(size_t dim, 
+					 bopt_params params):
+    HierarchicalGaussianProcess(dim,params),
+    mAlpha(params.alpha), mBeta (params.beta), 
+    mW0(params.mean.n_coef), mInvVarW(params.mean.n_coef), 
+    mD(params.mean.n_coef,params.mean.n_coef)
+  {  
+    mW0 = utils::array2vector(params.mean.coef_mean,params.mean.n_coef);
+    for (size_t ii = 0; ii < params.mean.n_coef; ++ii)
+      {
+	double varii = params.mean.coef_std[ii] * params.mean.coef_std[ii];
+	mInvVarW(ii) = 1/varii;
+      }
+     d_ = new StudentTDistribution();
+  }  // Constructor
+
+
+
+  StudentTProcessNIG::~StudentTProcessNIG()
+  {
+    delete d_;
+  } // Default destructor
+
+
+  ProbabilityDistribution* StudentTProcessNIG::prediction(const vectord &query)
+  {
+    double kq = (*mKernel)(query, query);;
+    vectord kn = computeCrossCorrelation(query);
+    vectord phi = mMean->getFeatures(query);
+  
+    vectord v(kn);
+    inplace_solve(mL,v,ublas::lower_tag());
+
+    vectord rq = phi - prod(v,mKF);
+
+    vectord rho(rq);
+    inplace_solve(mD,rho,ublas::lower_tag());
+    
+    double yPred = inner_prod(phi,mWMap) + inner_prod(v,mVf);
+    double sPred = sqrt( mSigmaMap * (kq - inner_prod(v,v) 
+				   + inner_prod(rho,rho)));
+
+    d_->setMeanAndStd(yPred,sPred);
+    return d_;
+  }
+
+
+  double StudentTProcessNIG::negativeLogLikelihood()
+  {
+    /* TODO: For testing */
+    return negativeTotalLogLikelihood();
+  }
+
+
+
+  int StudentTProcessNIG::precomputePrediction()
+  {
+    size_t n = mGPXX.size();
+    size_t p = mMean->nFeatures();
+
+    mKF = trans(mFeatM);
+    inplace_solve(mL,mKF,ublas::lower_tag());
+    //TODO: make one line
+    matrixd DD(p,p);
+    DD = prod(trans(mKF),mKF);
+    utils::addToDiagonal(DD,mInvVarW);
+    utils::cholesky_decompose(DD,mD);
+
+    vectord vn = mGPY;
+    inplace_solve(mL,vn,ublas::lower_tag());
+    mWMap = prod(mFeatM,vn) + utils::ublas_elementwise_prod(mInvVarW,mW0);
+    utils::cholesky_solve(mD,mWMap,ublas::lower());
+
+    mVf = mGPY - prod(trans(mFeatM),mWMap);
+    inplace_solve(mL,mVf,ublas::lower_tag());
+
+    vectord v0 = mGPY - prod(trans(mFeatM),mW0);
+    //TODO: check for "cheaper" version
+    matrixd KK = prod(mL,trans(mL));
+    matrixd WW = zmatrixd(p,p);  //TODO: diagonal matrix
+    utils::addToDiagonal(WW,mInvVarW);
+    matrixd FW = prod(trans(mFeatM),WW);
+    KK += prod(FW,mFeatM);
+    matrixd BB(n,n);
+    utils::cholesky_decompose(KK,BB);
+    inplace_solve(BB,v0,ublas::lower_tag());
+    mSigmaMap = (mBeta/mAlpha + inner_prod(v0,v0))/(n+2*mAlpha);
+    
+    int dof = static_cast<int>(n+2*mAlpha);
+    
+    if (dof <= 0)  
+      {
+	FILE_LOG(logERROR) << "ERROR: Incorrect alpha. Dof invalid.";
+	dof = 1;
+      }
+
+    d_->setDof(dof);  
+    return 0;
+  }
+
+} //namespace bayesopt

File utils/ublas_extra.hpp

       return 0;
     };
 
+    template<class M, class V>
+    int addToDiagonal(M& mat, const V& vec)
+    {
+      assert(mat.size1()==mat.size2());
+      assert(mat.size1()==vec.size());
+      const size_t ll = vec.size();
+      for(size_t ii = 0; ii < ll; ++ii)
+	{
+	  mat(ii,ii) += vec(ii);
+	}
+      return 0;
+    }
+
     boost::numeric::ublas::vector<double> array2vector(const double array[], 
 						       const size_t n);