Commits

Ruben Martinez-Cantin committed cb6fa65

Gaussian process ml working!

Comments (0)

Files changed (20)

   ./src/nonparametricprocess.cpp
   ./src/hierarchical_gaussian_process.cpp
   ./src/gaussian_process.cpp
-  ./src/gaussian_process_ign.cpp
-  ./src/student_t_process.cpp
+  ./src/gaussian_process_ml.cpp
+#  ./src/gaussian_process_ign.cpp
+#  ./src/student_t_process.cpp
   ./src/parameters.cpp
   ./src/kernel_functors.cpp
   ./src/criteria_functors.cpp
   par.kernel.n_theta = 2;
   par.mean.mu[0] = 1.0;
   par.mean.mu[1] = 1.0;
-  par.mean.s_mu[0] = PRIOR_DELTA_SQ;
-  par.mean.s_mu[1] = PRIOR_DELTA_SQ;
+  par.mean.s_mu[0] = MEAN_SIGMA;
+  par.mean.s_mu[1] = MEAN_SIGMA;
   par.mean.n_mu = 2;
   par.alpha = PRIOR_ALPHA;
   par.beta = PRIOR_BETA;
   par.noise = DEFAULT_NOISE;
-  par.surr_name = S_STUDENT_T_PROCESS_JEFFREYS;
+  par.surr_name = S_GAUSSIAN_PROCESS_ML;//S_STUDENT_T_PROCESS_JEFFREYS;
   par.kernel.name = "kSum(kSEISO,kConst)";
   par.mean.name = "mSum(mConst,mConst)";
   par.l_type = L_ML;
   par.alpha = PRIOR_ALPHA;
   par.beta = PRIOR_BETA;
   par.mean.mu[0] = 0.0;
-  par.mean.s_mu[0] = PRIOR_DELTA_SQ;
+  par.mean.s_mu[0] = MEAN_SIGMA;
   par.mean.n_mu = 1;
   par.noise = DEFAULT_NOISE;
   par.surr_name = S_GAUSSIAN_PROCESS_INV_GAMMA_NORMAL;

include/gaussian_process.hpp

   class GaussianProcess: public NonParametricProcess
   {
   public:
-    GaussianProcess(size_t dim, double noise);
+    GaussianProcess(size_t dim, bopt_params params);
     virtual ~GaussianProcess();
 
     /** 
 
   private:
     vectord mAlphaV;              ///< Precomputed L\y
+    double mSigma;                ///< Signal variance
     GaussianDistribution* d_;     ///< Pointer to distribution function
   };
 

include/gaussian_process_ign.hpp

   class GaussianProcessIGN: public HierarchicalGaussianProcess
   {
   public:
-    GaussianProcessIGN(size_t dim, double noise, double alpha,
-		       double beta, double delta);
+    GaussianProcessIGN(size_t dim, bopt_params params);
 
     virtual ~GaussianProcessIGN();
 

include/gaussian_process_ml.hpp

 #ifndef __GAUSSIAN_PROCESS_ML_HPP__
 #define __GAUSSIAN_PROCESS_ML_HPP__
 
-#include "nonparametricprocess.hpp"
+#include "gauss_distribution.hpp"
+#include "hierarchical_gaussian_process.hpp"
 
 
 namespace bayesopt
    */
   class GaussianProcessML: public HierarchicalGaussianProcess
   {
-    GaussianProcessML(size_t dim, double noise);
+  public:
+    GaussianProcessML(size_t dim, bopt_params params);
     virtual ~GaussianProcessML();
 
     /** 
      * negativeTotalLogLikelihood
      * @return value negative log likelihood
      */
-    double negativeLogLikelihood()
-    { return negativeTotalLogLikelihood(); };
+    double negativeLogLikelihood();
 
     /** 
      * \brief Precompute some values of the prediction that do not depends on
     int precomputePrediction();
 
   private:
-    double mWML, mSigML;           //!< GP ML parameters
-
+    vectord mWML;           //!< GP ML parameters
+    double mSigML;           //!< GP ML parameters
+    
     /// Precomputed GP prediction operations
+    vectord mAlphaF;
+    matrixd mKF, mL2;
 
     GaussianDistribution* d_;      //!< Predictive distributions
   };

include/hierarchical_gaussian_process.hpp

   class HierarchicalGaussianProcess: public NonParametricProcess
   {
   public:
-    HierarchicalGaussianProcess(size_t dim, double noise);
+    HierarchicalGaussianProcess(size_t dim, bopt_params params);
     virtual ~HierarchicalGaussianProcess() {};
 
   protected:

include/inneroptimization.hpp

 #include "dll_stuff.h"
 #include "specialtypes.hpp"
 
-// We plan to add more in the future since nlopt actually support many of them
-enum innerOptAlgorithms {
-  DIRECT,    ///< Global optimization
-  LBFGS,     ///< Local, derivative based
-  BOBYQA,    ///< Local, derivative free
-  COMBINED   ///< Global exploration, local refinement
-};
+namespace bayesopt {
 
+  // We plan to add more in the future since nlopt actually support many of them
+  typedef enum {
+    DIRECT,    ///< Global optimization
+    LBFGS,     ///< Local, derivative based
+    BOBYQA,    ///< Local, derivative free
+    COMBINED   ///< Global exploration, local refinement
+  } innerOptAlgorithms;
 
-class BAYESOPT_API InnerOptimization
-{
-public:
-  InnerOptimization()
-  { 
-    alg = DIRECT;    mDown = 0.;    mUp = 1.;
+
+  class BAYESOPT_API InnerOptimization
+  {
+  public:
+    InnerOptimization()
+      { 
+	alg = DIRECT;    mDown = 0.;    mUp = 1.;
+      };
+
+    virtual ~InnerOptimization(){};
+
+    /** 
+     * Set the optimization algorithm
+     * 
+     * @param newAlg 
+     */
+    void setAlgorithm(innerOptAlgorithms newAlg)
+    { alg = newAlg; }
+
+    /** 
+     * Limits of the hypercube. 
+     * Currently, it assumes that all dimensions have the same limits.
+     * 
+     * @param down 
+     * @param up 
+     */
+    void setLimits(double down, double up)
+    {
+      mDown = down;   mUp = up;
+    }
+
+    /** 
+     * Compute the inner optimization algorithm
+     * 
+     * @param Xnext input: initial guess, output: result
+     * 
+     * @return error_code
+     */
+    int innerOptimize(vectord &Xnext);
+    int innerOptimize(double* x, int n, void* objPointer);	
+
+
+    /** 
+     * Virtual function to be overriden by the actual function to be evaluated
+     * 
+     * @param query input point
+     * 
+     * @return function value at query point
+     */
+    virtual double innerEvaluate(const vectord& query) 
+    {return 0.0;};
+
+
+    /** 
+     * Virtual function to be overriden by the actual function to be evaluated
+     * 
+     * @param query input point
+     * @param grad output gradient at query point
+     * 
+     * @return function value at query point
+     */
+    virtual double innerEvaluate(const vectord& query, 
+				 vectord& grad) 
+    {return 0.0;};
+
+
+  protected:
+
+    innerOptAlgorithms alg;
+    double mDown, mUp;
   };
 
-  virtual ~InnerOptimization(){};
-
-  /** 
-   * Set the optimization algorithm
-   * 
-   * @param newAlg 
-   */
-  void setAlgorithm(innerOptAlgorithms newAlg)
-  { alg = newAlg; }
-
-  /** 
-   * Limits of the hypercube. 
-   * Currently, it assumes that all dimensions have the same limits.
-   * 
-   * @param down 
-   * @param up 
-   */
-  void setLimits(double down, double up)
-  {
-    mDown = down;   mUp = up;
-  }
-
-  /** 
-   * Compute the inner optimization algorithm
-   * 
-   * @param Xnext input: initial guess, output: result
-   * 
-   * @return error_code
-   */
-  int innerOptimize(vectord &Xnext);
-  int innerOptimize(double* x, int n, void* objPointer);	
-
-
-  /** 
-   * Virtual function to be overriden by the actual function to be evaluated
-   * 
-   * @param query input point
-   * 
-   * @return function value at query point
-   */
-  virtual double innerEvaluate(const vectord& query) 
-  {return 0.0;};
-
-
-  /** 
-   * Virtual function to be overriden by the actual function to be evaluated
-   * 
-   * @param query input point
-   * @param grad output gradient at query point
-   * 
-   * @return function value at query point
-   */
-  virtual double innerEvaluate(const vectord& query, 
-			       vectord& grad) 
-  {return 0.0;};
-
-
-protected:
-
-  innerOptAlgorithms alg;
-  double mDown, mUp;
-};
+}//namespace bayesopt
 
 #endif

include/nonparametricprocess.hpp

   class NonParametricProcess: public InnerOptimization
   {
   public:
-    NonParametricProcess(size_t dim, double noise);
+    NonParametricProcess(size_t dim, bopt_params parameters);
     virtual ~NonParametricProcess();
 
     /** 
     const double mRegularizer;   ///< Std of the obs. model (also used as nugget)
     vecOfvec mGPXX;                                              ///< Data inputs
     vectord mGPY;                                                ///< Data values
+    
     vectord mMeanV;                           ///< Mean value at the input points
     matrixd mFeatM;           ///< Value of the mean features at the input points
+    vectord mMu;                 ///< Mean of the parameters of the mean function
+    vectord mS_Mu;    ///< Variance of the params of the mean function W=mS_Mu*I
 
     std::vector<boost::math::normal> priorKernel; ///< Prior of kernel parameters
     boost::scoped_ptr<Kernel> mKernel;            ///< Pointer to kernel function

include/parameters.h

 
   typedef enum {  
     S_GAUSSIAN_PROCESS,
+    S_GAUSSIAN_PROCESS_ML,
     S_GAUSSIAN_PROCESS_INV_GAMMA_NORMAL,
     S_STUDENT_T_PROCESS_JEFFREYS,
     S_ERROR = -1
     char* log_filename;          /**< Log file path (if applicable) */
 
     surrogate_name surr_name;    /**< Name of the surrogate function */
+    double sigma_s;              /**< Signal variance (if known) */
+    double noise;                /**< Observation noise (and nugget) */
     double alpha;                /**< Inverse Gamma prior for signal var */
     double beta;                 /**< Inverse Gamma prior for signal var*/
-    double noise;                /**< Observation noise (and nugget) */
     learning_type l_type;        /**< Type of learning for the kernel params*/
 
     kernel_parameters kernel;    /**< Kernel parameters */
     mean_parameters mean;        /**< Mean (parametric function) parameters */
 
     char* crit_name;             /**< Name of the criterion */
-
+    double crit_params[128];     /**< Criterion hyperparameters (if needed) */
+    size_t n_crit_params;        /**< Number of criterion hyperparameters */
   } bopt_params;
 
 
   const double MEAN_SIGMA      = 1000.0;
   const double PRIOR_ALPHA     = 1.0;
   const double PRIOR_BETA      = 1.0;
-  const double PRIOR_DELTA_SQ  = 1000.0;
+  const double DEFAULT_SIGMA   = 1.0;
   const double DEFAULT_NOISE   = 1e-4;
 
   /* Algorithm parameters */

include/student_t_process.hpp

   class StudentTProcess: public HierarchicalGaussianProcess
   {
   public:
-    StudentTProcess(size_t dim, double noise);
+    StudentTProcess(size_t dim, bopt_params params);
     virtual ~StudentTProcess();
 
     /** 

src/gaussian_process.cpp

 #include "gaussian_process.hpp"
 #include "cholesky.hpp"
 #include "trace_ublas.hpp"
-#include "gauss_distribution.hpp"
 
 namespace bayesopt
 {
 
   namespace ublas = boost::numeric::ublas;
 
-  GaussianProcess::GaussianProcess(size_t dim, double noise):
-    NonParametricProcess(dim, noise)
+  GaussianProcess::GaussianProcess(size_t dim, bopt_params params):
+    NonParametricProcess(dim, params), mSigma(params.sigma_s)
   {
     d_ = new GaussianDistribution();
   }  // Constructor
     matrixd L(n,n);
     utils::cholesky_decompose(K,L);
 
-    vectord alpha(mGPY-mMeanV);
+    vectord alpha(mGPY-prod(mMu,mFeatM));
     inplace_solve(L,alpha,ublas::lower_tag());
-    double loglik = .5*ublas::inner_prod(alpha,alpha) + utils::log_trace(L);
+    double loglik = ublas::inner_prod(alpha,alpha)/(2*mSigma) + 
+      utils::log_trace(L);
     return loglik;
   }
 
   ProbabilityDistribution* GaussianProcess::prediction(const vectord &query)
   {
-    const double kn = (*mKernel)(query, query);
-    const vectord kStar = computeCrossCorrelation(query);
-    double yPred = ublas::inner_prod(kStar,mAlphaV);
+    const double kq = (*mKernel)(query, query);
+    const vectord kn = computeCrossCorrelation(query);
 
-    vectord vd(kStar);
+    vectord vd(kn);
     ublas::inplace_solve(mL,vd,ublas::lower_tag());
-    double sPred = sqrt(kn - ublas::inner_prod(vd,vd));
+    double yPred = ublas::inner_prod(vd,mAlphaV);
+    double sPred = sqrt(mSigma*(kq - ublas::inner_prod(vd,vd)));
 
     d_->setMeanAndStd(yPred,sPred);
     return d_;
     const size_t n = mGPY.size();
   
     mAlphaV.resize(n,false);
-    mAlphaV = mGPY-mMeanV;
+    mAlphaV = mGPY-prod(mMu,mFeatM);
     utils::cholesky_solve(mL,mAlphaV,ublas::lower());
 
     return 1; 

src/gaussian_process_ign.cpp

   using utils::cholesky_solve;
   using utils::cholesky_decompose;
 
-  GaussianProcessIGN::GaussianProcessIGN(size_t dim, double noise, double alpha, 
-					 double beta, double delta):
-    HierarchicalGaussianProcess(dim,noise),
-    mAlpha(alpha), mBeta (beta), mDelta2(delta)
+  GaussianProcessIGN::GaussianProcessIGN(size_t dim, bopt_params params):
+    GaussianProcess(dim,params), mAlpha(params.alpha), 
+    mBeta (params.beta), mDelta2(parameters.mean.s_mu[0])
   {
     d_ = new GaussianDistribution();
   }  // Constructor

src/gaussian_process_ml.cpp

 -----------------------------------------------------------------------------
 */
 
-#include "gaussian_process_ml.hpp"
+
 #include "cholesky.hpp"
 #include "trace_ublas.hpp"
-#include "gauss_distribution.hpp"
+#include "gaussian_process_ml.hpp"
 
-using boost::numeric::ublas::inplace_solve;
-using boost::numeric::ublas::lower_tag;
-using boost::numeric::ublas::lower;
+namespace bayesopt
+{
+
+  using boost::numeric::ublas::inplace_solve;
+  using boost::numeric::ublas::inner_prod;
+  using boost::numeric::ublas::lower_tag;
+  using boost::numeric::ublas::lower;
+  using boost::numeric::ublas::trans;
   
 
-GaussianProcessML::GaussianProcessIGN(double noise):
-  NonParametricProcess(noise)
-{
-  d_.reset(new GaussianDistribution());
-}  // Constructor
+  GaussianProcessML::GaussianProcessML(size_t dim, bopt_params params):
+    HierarchicalGaussianProcess(dim, params)
+  {
+    d_ = new GaussianDistribution();
+  }  // Constructor
 
 
 
-GaussianProcessIGN::~GaussianProcessIGN()
-{} // Default destructor
+  GaussianProcessML::~GaussianProcessML()
+  {
+    delete d_;
+  } // Default destructor
 
 
 
 
-double GaussianProcessIGN::negativeLogLikelihood()
-{
-  matrixd K = computeCorrMatrix();
-  size_t n = K.size1();
-  matrixd L(n,n);
-  utils::cholesky_decompose(K,L);
+  double GaussianProcessML::negativeLogLikelihood()
+  {
+    /* In this case, they are equivalent */
+    return negativeTotalLogLikelihood();
+  }
 
-  matrixd featAll = mMean->getFeatures(mGPXX);
 
-  vectord alphU(mMeanV);
-  inplace_solve(L,alphU,lower_tag());
-  double eta = inner_prod(alphU,alphU) + 1/mDelta2;
+  ProbabilityDistribution* GaussianProcessML::prediction( const vectord &query )
+  {
+    double kq = (*mKernel)(query, query);;
+    vectord kn = computeCrossCorrelation(query);
+    vectord phi = mMean->getFeatures(query);
   
-  vectord alphY(mGPY);
-  inplace_solve(L,alphY,lower_tag());
-  double mu     = inner_prod(alphU,alphY) / eta;
-  double YInvRY = inner_prod(alphY,alphY);
+    vectord v(kn);
+    inplace_solve(mL,v,lower_tag());
+
+    vectord rq = phi - prod(v,mKF);
+
+    vectord rho(rq);
+    inplace_solve(mL2,rho,lower_tag());
     
-  double sigma = (mBeta + YInvRY - mu*mu*eta) / (mAlpha + (n+1) + 2);
+    double yPred = inner_prod(phi,mWML) + inner_prod(v,mAlphaF);
+    double sPred = sqrt( mSigML * (kq - inner_prod(v,v) 
+				   + inner_prod(rho,rho)));
+
+    d_->setMeanAndStd(yPred,sPred);
+    return d_;
+  }
+
+  int GaussianProcessML::precomputePrediction()
+  {
+    size_t n = mGPXX.size();
+    size_t p = mMean->nFeatures();
+
+    mKF = trans(mFeatM);
+    inplace_solve(mL,mKF,lower_tag());
+
+    matrixd FKF = prod(trans(mKF),mKF);
+    mL2 = FKF;
+    utils::cholesky_decompose(FKF,mL2);
+
+    vectord Ky(mGPY);
+    inplace_solve(mL,Ky,lower_tag());
+
+    mWML = prod(Ky,mKF);
+    utils::cholesky_solve(mL2,mWML,lower());
+
+    mAlphaF = mGPY - prod(mWML,mFeatM);
+    inplace_solve(mL,mAlphaF,lower_tag());
+    mSigML = inner_prod(mAlphaF,mAlphaF)/(n-p);
   
-  vectord yumu = mGPY - mMeanV*mu;
-  alphY = yumu;
-  inplace_solve(L,alphY,lower_tag());
+    return 1;
+  }
 
-  double lik1 = inner_prod(alphY,alphY) / (2*sigma); 
-  double lik2 = log_trace(L) + 0.5*n*log(sigma) + n*0.91893853320467; //log(2*pi)/2
-  
-  //TODO: This must be wrong.
-  size_t index = 1;
-  vectord th = mKernel->getScale(index);
-
-  return lik1 + lik2 + mBeta/2 * th(0) -
-    (mAlpha+1) * log(th(0));
-}
-
-
-ProbabilityDistribution* GaussianProcessIGN::prediction( const vectord &query )
-{
-  double kn;
-  double uInvRr, rInvRr, rInvRy;
-  double meanf = mMean->getMean(query);
-
-  vectord colR = computeCrossCorrelation(query);
-  kn = (*mKernel)(query, query);
-
-#if USE_CHOL
-  vectord invRr(colR);
-  inplace_solve(mL,invRr,lower_tag());
-  rInvRr = inner_prod(invRr,invRr);
-  uInvRr = inner_prod(mUInvR, invRr);  
-  rInvRy = inner_prod(invRr, mInvRy );
-#else
-  vectord rInvR = prod(colR,mInvR);
-	
-  rInvRr = inner_prod(rInvR,colR);
-  uInvRr = inner_prod(mUInvR,colR);
-  rInvRy = inner_prod(colR,mInvRy);
-#endif
-  
-  double yPred = meanf*mMu + rInvRy;
-  double sPred = sqrt( mSig * (kn - rInvRr + (meanf - uInvRr) * (meanf - uInvRr) 
-			       / mUInvRUDelta ) );
-
-  d_->setMeanAndStd(yPred,sPred);
-  return d_.get();
-}
-
-int GaussianProcessIGN::precomputePrediction()
-{
-  size_t ns = mGPXX.size();
-  size_t nf = mMean->nFeatures();
-
-  mFKF =
-
-  mUInvR.resize(nSamples,false);
-  mInvRy.resize(nSamples,false);
-
-#if USE_CHOL
-  mAlphaV.resize(nSamples,false);
-  mAlphaV = mGPY;
-  utils::cholesky_solve(mL,mAlphaV,lower());
-
-  //  vectord alphaMean(mMeanV);
-  mUInvR = mMeanV;
-  inplace_solve(mL,mUInvR,lower_tag());
-  mUInvRUDelta = inner_prod(mUInvR,mUInvR) + 1/mDelta2;
-
-  mMu = inner_prod(mMeanV,mAlphaV) / mUInvRUDelta;
-  double YInvRY = inner_prod(mGPY,mAlphaV);
-
-  mInvRy = mGPY - mMeanV*mMu;
-  inplace_solve(mL,mInvRy,lower_tag());
-#else
-  // mInvR.resize(nSamples,nSamples);
-  // mInvR.assign(boost::numeric::ublas::identity_matrix<double>(nSamples));
-  // cholesky_solve(mL,mInvR,lower());
-  mUInvR = prod(mMeanV,mInvR);
-  mUInvRUDelta = inner_prod(mUInvR,mMeanV) + 1/mDelta2;
-  
-  vectord YInvR(nSamples);
-  double YInvRY;
-  
-  mMu =  inner_prod(mUInvR,mGPY) / mUInvRUDelta;
-  
-  noalias(YInvR) = prod(mGPY,mInvR);
-  YInvRY = inner_prod(YInvR,mGPY);
-  vectord ymu = mGPY - mMeanV*mMu;
-  mInvRy = prod(mInvR,ymu);
-#endif  
-
-  mSig = (mBeta + YInvRY - mMu*mMu*mUInvRUDelta) / (mAlpha + (nSamples+1) + 2);
-  
-  return 1;
-}
+} //namespace bayesopt

src/hierarchical_gaussian_process.cpp

   namespace ublas = boost::numeric::ublas;
 
   HierarchicalGaussianProcess::HierarchicalGaussianProcess(size_t dim, 
-							   double noise):
-    NonParametricProcess(dim, noise) {};
+							   bopt_params params):
+    NonParametricProcess(dim, params) {};
 
   double HierarchicalGaussianProcess::negativeTotalLogLikelihood()
   {
+    /*This is the restricted version. For the unrestricted, make p=0
+      and remove the last term of loglik*/
+
     const matrixd K = computeCorrMatrix();
     const size_t n = K.size1();
     const size_t p = mFeatM.size1(); 
     matrixd L(n,n);
     utils::cholesky_decompose(K,L);
 
-    matrixd KF(n,p);
+    matrixd KF(ublas::trans(mFeatM));
     inplace_solve(L,KF,ublas::lower_tag());
     
     matrixd FKF = prod(ublas::trans(KF),KF);
 
     vectord Ky(mGPY);
     inplace_solve(L,Ky,ublas::lower_tag());
+
     vectord wML = prod(Ky,KF);
-
     utils::cholesky_solve(L2,wML,ublas::lower());
 
     vectord alpha = mGPY - prod(wML,mFeatM);
     inplace_solve(L,alpha,ublas::lower_tag());
+    double sigma = ublas::inner_prod(alpha,alpha)/(n-p);
 
-    double loglik = .5*n*log(ublas::inner_prod(alpha,alpha)) 
-                    + utils::log_trace(L);
+    double loglik = .5*(n-p)*log(ublas::inner_prod(alpha,alpha)) 
+      + utils::log_trace(L) + utils::log_trace(L2);
     return loglik;
   }
 }

src/inneroptimization.cpp

 #include "log.hpp"
 #include "inneroptimization.hpp"
 
-void checkNLOPTerror(nlopt_result errortype)
+namespace bayesopt
 {
-  switch(errortype)
+  void checkNLOPTerror(nlopt_result errortype)
+  {
+    switch(errortype)
       {
       case -1: FILE_LOG(logERROR) << "NLOPT: General failure"; break;
       case -2: FILE_LOG(logERROR) << "NLOPT: Invalid arguments. Check bounds."; break;
       case -5: FILE_LOG(logERROR) << "NLOPT: Force stop."; break;
       default: ;
       }
-}
+  }
 
 
-int InnerOptimization::innerOptimize(vectord &Xnext)
-{   
+  int InnerOptimization::innerOptimize(vectord &Xnext)
+  {   
     void *objPointer = static_cast<void *>(this);
     int n = static_cast<int>(Xnext.size());
     int error;
     error = innerOptimize(&Xnext(0), n, objPointer);
 
     return error;
-} // innerOptimize (uBlas)
+  } // innerOptimize (uBlas)
 
-int InnerOptimization::innerOptimize(double* x, int n, void* objPointer)
-{
+  int InnerOptimization::innerOptimize(double* x, int n, void* objPointer)
+  {
     double u[128], l[128];
     double fmin = 1;
     int maxf = MAX_INNER_EVALUATIONS*n;    
     int ierror;
 
     for (int i = 0; i < n; ++i) {
-	l[i] = mDown;	u[i] = mUp;
-	// What if x is undefined?
-	if (x[i] < l[i] || x[i] > u[i])
-	  x[i]=(l[i]+u[i])/2.0;
+      l[i] = mDown;	u[i] = mUp;
+      // What if x is undefined?
+      if (x[i] < l[i] || x[i] > u[i])
+	x[i]=(l[i]+u[i])/2.0;
     }
 
     nlopt_opt opt;
     ierror = static_cast<int>(errortype);
     return ierror;
 
-} // innerOptimize (C array)
+  } // innerOptimize (C array)
 
 
+}// namespace bayesopt
+

src/nonparametricprocess.cpp

 #include "ublas_extra.hpp"
 
 #include "gaussian_process.hpp"
+#include "gaussian_process_ml.hpp"
 #include "gaussian_process_ign.hpp"
 #include "student_t_process.hpp"
 
 namespace bayesopt
 {
   
-  NonParametricProcess::NonParametricProcess(size_t dim, double noise):
-    InnerOptimization(), mRegularizer(noise), dim_(dim)
+  NonParametricProcess::NonParametricProcess(size_t dim, bopt_params parameters):
+    InnerOptimization(), mRegularizer(parameters.noise), dim_(dim)
   { 
     mMinIndex = 0;     mMaxIndex = 0;   
     setAlgorithm(BOBYQA);
     setLimits(0.,100.);
+    setLearnType(parameters.l_type);
+    setKernel(parameters.kernel,dim);
+    setMean(parameters.mean,dim);
   }
 
   NonParametricProcess::~NonParametricProcess(){}
     switch(parameters.surr_name)
       {
       case S_GAUSSIAN_PROCESS: 
-	s_ptr = new GaussianProcess(dim,parameters.noise); 
-	break;
+	s_ptr = new GaussianProcess(dim,parameters); break;
 
-      case S_GAUSSIAN_PROCESS_INV_GAMMA_NORMAL:
-	s_ptr = new GaussianProcessIGN(dim,parameters.noise, parameters.alpha,
-				       parameters.beta,parameters.mean.s_mu[0]);  
-	break;
+      case S_GAUSSIAN_PROCESS_ML: 
+	s_ptr = new GaussianProcessML(dim,parameters); break;
 
-      case S_STUDENT_T_PROCESS_JEFFREYS: 
-	if (strcmp(parameters.mean.name,"mZero") == 0)
-	  {
-	    FILE_LOG(logWARNING) << "Zero mean incompatible with Student's t "
-				 << "process, using one-mean instead.";
-	    parameters.mean.name = "mOne";
-	  }
-	s_ptr = new StudentTProcess(dim,parameters.noise); 
-	break;
+      // case S_GAUSSIAN_PROCESS_INV_GAMMA_NORMAL:
+      // 	s_ptr = new GaussianProcessIGN(dim,parameters);	break;
+
+      // case S_STUDENT_T_PROCESS_JEFFREYS: 
+      // 	s_ptr = new StudentTProcess(dim,parameters); break;
 
       default:
 	FILE_LOG(logERROR) << "Error: surrogate function not supported.";
 	return NULL;
       }
 
-    s_ptr->setLearnType(parameters.l_type);
-    s_ptr->setKernel(parameters.kernel,dim);
-    s_ptr->setMean(parameters.mean,dim);
     return s_ptr;
   };
 
 				     size_t dim)
   {
     mMean.reset(mPFactory.create(m_name,dim));
+    mMu = muv; mS_Mu = smu;
 
     if (mMean == NULL) 	return -1; 
 

src/parameters.cpp

 {
   if      (!strcmp(name,  "GAUSSIAN_PROCESS"))
     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_JEFFREYS"))
   switch(name)
     {
     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_ERROR:
 static const bopt_params DEFAULT_PARAMS = {
   DEFAULT_ITERATIONS, MAX_INNER_EVALUATIONS, DEFAULT_SAMPLES, 
   DEFAULT_VERBOSE, DEF_LOG_FILE,
-  S_GAUSSIAN_PROCESS,
-  PRIOR_ALPHA, PRIOR_BETA,  DEFAULT_NOISE,
+  S_GAUSSIAN_PROCESS, 
+  DEFAULT_SIGMA, DEFAULT_NOISE,
+  PRIOR_ALPHA, PRIOR_BETA, 
   L_MAP,
   DEFAULT_KERNEL, DEFAULT_MEAN,
-  DEF_CRITERIA_NAME,
+  DEF_CRITERIA_NAME, {}, 0
 };
 
 

src/student_t_process.cpp

   using utils::cholesky_decompose;
   
   
-  StudentTProcess::StudentTProcess(size_t dim, double noise):
-    HierarchicalGaussianProcess(dim, noise)
+  StudentTProcess::StudentTProcess(size_t dim, bopt_params params):
+    HierarchicalGaussianProcess(dim, params)
   {
     d_ = new StudentTDistribution();
+    if (strcmp(parameters.mean.name,"mZero") == 0)
+      {
+	FILE_LOG(logWARNING) << "Zero mean incompatible with Student's t "
+			     << "process, using one-mean instead.";
+	parameters.mean.name = "mOne";
+      }
+    setMean(parameters.mean,dim);
   }  // Constructor
 
 

wrappers/nloptwpr.cpp

 namespace NLOPT_WPR
 {
 
-  using namespace boost::numeric::ublas;	
+  namespace ublas = boost::numeric::ublas;
+  using bayesopt::InnerOptimization;
 
 
   double evaluate_nlopt (unsigned int n, const double *x,
     double xcopy[128];
     for (unsigned int i=0;i<n;i++)
       xcopy[i] = x[i];
-    array_adaptor<double> shared(n, xcopy);
-    vector<double, array_adaptor<double> > sharedN(n, shared); 
+    ublas::array_adaptor<double> shared(n, xcopy);
+    ublas::vector<double, ublas::array_adaptor<double> > sharedN(n, shared); 
     
     // This is not very clever... but works!
     void *objPointer = my_func_data;
     double xcopy[128];
     for (unsigned int i=0;i<n;i++)
       xcopy[i] = x[i];
-    array_adaptor<double> shared(n, xcopy);
-    vector<double, array_adaptor<double> > sharedN(n, shared); 
+    ublas::array_adaptor<double> shared(n, xcopy);
+    ublas::vector<double, ublas::array_adaptor<double> > sharedN(n, shared); 
     
     // This is not very clever... but works!
     void *objPointer = my_func_data;
     InnerOptimization* OPTIMIZER = static_cast<InnerOptimization*>(objPointer);
     
-    vector<double> vgrad = zero_vector<double>(n);
+    ublas::vector<double> vgrad = ublas::zero_vector<double>(n);
 
     double f =  OPTIMIZER->innerEvaluate(sharedN,vgrad);