Commits

Ruben Martinez-Cantin committed d0c84f6

Refactoring more elements. Tricks on prediction does speed up things. InnerOptimization encapsulation vs inheritance has a strong penalty.

  • Participants
  • Parent commits c2bd173

Comments (0)

Files changed (15)

examples/bo_cont.cpp

   par.l_type = L_ML;
   par.n_iterations = 200;       // Number of iterations
   par.n_init_samples = 50;
-  par.verbose_level = 2;
+  par.verbose_level = 0;
   /*******************************************/
 
   clock_t start, end;

include/bayesoptcont.hpp

      * @param query point to evaluate in [0,1] hypercube
      * @return actual return value of the target function
      */
-    inline double evaluateSampleInternal( const vectord &query )
-    { 
-      vectord unnormalizedQuery = mBB->unnormalizeVector(query);
-      return evaluateSample(unnormalizedQuery);
-    }; // evaluateSampleInternal
-
+    double evaluateSampleInternal( const vectord &query );
+    
     /** 
      * \brief Wrapper of the innerOptimization class to find the optimal
      * point acording to the criteria.
      * @param xOpt optimal point
      * @return error code
      */
-    inline int findOptimal(vectord &xOpt)
-    { return cOptimizer->run(xOpt); };
-
+    int findOptimal(vectord &xOpt);
 
   private:
     utils::BoundingBox<vectord> *mBB;      ///< Bounding Box (input space limits)
   
   /**@}*/
 
+
+  inline double ContinuousModel::evaluateSampleInternal( const vectord &query )
+  { 
+    vectord unnormalizedQuery = mBB->unnormalizeVector(query);
+    return evaluateSample(unnormalizedQuery);
+  }; // evaluateSampleInternal
+
+  inline int ContinuousModel::findOptimal(vectord &xOpt)
+  { return cOptimizer->run(xOpt); };
+
+
+
 }  //namespace bayesopt
 
 

include/kernel_functors.hpp

     size_t n_inputs;
   };
 
+
+
   template <typename KernelType> Kernel * create_func()
   {
     return new KernelType();
   }
 
-
   /** 
    * \brief Factory model for kernel functions
    * This factory is based on the libgp library by Manuel Blum
 
     Kernel* getKernel();
     
-    inline int setHyperParameters(const vectord &theta)
-    { return mKernel->setHyperParameters(theta); };
-    
-    inline vectord getHyperParameters(){return mKernel->getHyperParameters();};
-    inline size_t nHyperParameters(){return mKernel->nHyperParameters();};
+    int setHyperParameters(const vectord &theta);
+    vectord getHyperParameters();
+    size_t nHyperParameters();
+
 
     /** 
      * \brief Select kernel (covariance function) for the surrogate process.
     int computeCorrMatrix(const vecOfvec& XX, matrixd& corrMatrix, double nugget);
     int computeDerivativeCorrMatrix(const vecOfvec& XX, matrixd& corrMatrix, int dth_index);
     vectord computeCrossCorrelation(const vecOfvec& XX, const vectord &query);
+    void computeCrossCorrelation(const vecOfvec& XX, 
+				 const vectord &query,
+				 vectord& knx);
     double computeSelfCorrelation(const vectord& query);
     double kernelLogPrior();
 
 
     boost::scoped_ptr<Kernel> mKernel;            ///< Pointer to kernel function
     std::vector<boost::math::normal> priorKernel; ///< Prior of kernel parameters
-    KernelFactory mKFactory;
+  };
 
-  };
+  inline Kernel* KernelModel::getKernel()
+  { return mKernel.get();  }
+
+  inline int KernelModel::setHyperParameters(const vectord &theta)
+  { return mKernel->setHyperParameters(theta); };
+    
+  inline vectord KernelModel::getHyperParameters()
+  {return mKernel->getHyperParameters();};
+  
+  inline size_t KernelModel::nHyperParameters()
+  {return mKernel->nHyperParameters();};
+
+  inline vectord KernelModel::computeCrossCorrelation(const vecOfvec& XX, 
+						      const vectord &query)
+  {
+    vectord knx(XX.size());
+    computeCrossCorrelation(XX,query,knx);
+    return knx;
+  }
+
+  inline void KernelModel::computeCrossCorrelation(const vecOfvec& XX, 
+						   const vectord &query,
+						   vectord& knx)
+  {
+    std::vector<vectord>::const_iterator x_it  = XX.begin();
+    vectord::iterator k_it = knx.begin();
+    while(x_it != XX.end())
+      {	*k_it++ = (*mKernel)(*x_it++, query); }
+  }
+
+
+  inline double KernelModel::computeSelfCorrelation(const vectord& query)
+  { return (*mKernel)(query,query); }
 
   //@}
 

include/mean_functors.hpp

 
 
     virtual size_t nFeatures() = 0;
-    virtual vectord getFeatures(const vectord& x) = 0;  
+    virtual vectord getFeatures(const vectord& x) = 0;
     virtual matrixd getAllFeatures(const vecOfvec& x)
     {
       size_t nf = nFeatures();
     virtual ~MeanModel() {};
 
     ParametricFunction* getMeanFunc();
-    
-    inline int setParameters(const vectord &theta)
-    { return mMean->setParameters(theta); };
-    
-    inline vectord getParameters(){return mMean->getParameters();};
-    inline size_t nParameters(){return mMean->nParameters();};
+   
+    int setParameters(const vectord &theta);
+    vectord getParameters();
+    size_t nParameters();
 
-    inline void setPoints(const vecOfvec &x)
-    { mFeatM = mMean->getAllFeatures(x); };
+    vectord getFeatures(const vectord& x);  
+    void getFeatures(const vectord& x, vectord& kx);  
+    size_t nFeatures();
 
-    inline void addNewPoint(const vectord &x)
-    { 
-      using boost::numeric::ublas::column;
+    void setPoints(const vecOfvec &x);
+    void addNewPoint(const vectord &x);
 
-      vectord feat = mMean->getFeatures(x);
-      mFeatM.resize(feat.size(),mFeatM.size2()+1);  
-      column(mFeatM,mFeatM.size2()-1) = feat;
-    }
-
-    inline vectord muTimesFeat()
-    {  return boost::numeric::ublas::prod(mMu,mFeatM); }
-    
-    inline double muTimesFeat(const vectord& x)
-    { return boost::numeric::ublas::inner_prod(mMu,mMean->getFeatures(x));}
-
+    vectord muTimesFeat();
+    double muTimesFeat(const vectord& x);
 
     /** 
      * \brief Select the parametric part of the surrogate process.
     vectord mS_Mu;    ///< Variance of the params of the mean function W=mS_Mu*I
 
     boost::scoped_ptr<ParametricFunction> mMean;    ///< Pointer to mean function   
-    MeanFactory mPFactory;
+  };
 
-  };
+
+
+  inline ParametricFunction* MeanModel::getMeanFunc()
+  { return mMean.get(); }
+
+  inline int MeanModel::setParameters(const vectord &theta)
+  { return mMean->setParameters(theta); };
+    
+  inline vectord MeanModel::getParameters(){return mMean->getParameters();};
+  inline size_t MeanModel::nParameters(){return mMean->nParameters();};
+
+  inline vectord MeanModel::getFeatures(const vectord& x) 
+  { return mMean->getFeatures(x); }
+
+  inline void MeanModel::getFeatures(const vectord& x, vectord& kx)
+  { kx = mMean->getFeatures(x); }  
+
+  inline size_t MeanModel::nFeatures()
+  { return mMean->nFeatures(); }  
+
+  inline void MeanModel::setPoints(const vecOfvec &x)
+  { mFeatM = mMean->getAllFeatures(x); };
+
+  inline void MeanModel::addNewPoint(const vectord &x)
+  { 
+    using boost::numeric::ublas::column;
+    
+    mFeatM.resize(mFeatM.size1(),mFeatM.size2()+1);  
+    column(mFeatM,mFeatM.size2()-1) = mMean->getFeatures(x);
+  }
+
+  inline vectord MeanModel::muTimesFeat()
+  {  return boost::numeric::ublas::prod(mMu,mFeatM); }
+    
+  inline double MeanModel::muTimesFeat(const vectord& x)
+  { return boost::numeric::ublas::inner_prod(mMu,mMean->getFeatures(x));}
+
+
 
   //@}
 

include/nonparametricprocess.hpp

   class Dataset
   {
   public:
-    Dataset():  
-      mMinIndex(0), mMaxIndex(0) 
-    {};
+    Dataset();
+    Dataset(const matrixd& x, const vectord& y);
+    virtual ~Dataset();
 
+    void setSamples(const matrixd &x, const vectord &y);
+    void addSample(const vectord &x, double y);
+    double getSample(size_t index, vectord &x);
+    double getLastSample(vectord &x);
 
-    Dataset(const matrixd& x, const vectord& y):
-      mMinIndex(0), mMaxIndex(0)
-    {
-      setSamples(x,y);
-    };
-
-    virtual ~Dataset(){};
-
-    void setSamples(const matrixd &x, const vectord &y)
-    {
-      mY = y;
-      for (size_t i=0; i<x.size1(); ++i)
-	{
-	  mX.push_back(row(x,i));
-	  checkBoundsY(i);
-	} 
-    };
-  
-    void addSample(const vectord &x, double y)
-    {
-      mX.push_back(x);
-      mY.resize(mY.size()+1);  mY(mY.size()-1) = y;
-      checkBoundsY(mY.size()-1);
-    }
-
-    inline double getSample(size_t index, vectord &x)
-    { x = mX[index];  return mY(index);  }
-
-    inline double getLastSample(vectord &x)
-    { size_t last = mY.size()-1; return getSample(last, x); }
-
-    inline vectord getPointAtMinimum() { return mX[mMinIndex]; };
-    inline double getValueAtMinimum() { return mY(mMinIndex); };
-    inline size_t getNSamples() { return mY.size(); };
-    inline void checkBoundsY( size_t i )
-    {
-      if ( mY(mMinIndex) > mY(i) )       mMinIndex = i;
-      else if ( mY(mMaxIndex) < mY(i) )  mMaxIndex = i;
-    };
+    vectord getPointAtMinimum();
+    double getValueAtMinimum();
+    size_t getNSamples();
+    void checkBoundsY( size_t i );
 
     vecOfvec mX;                                         ///< Data inputs
     vectord mY;                                          ///< Data values
+
+  private:
     size_t mMinIndex, mMaxIndex;	
   };
 
-  /** \brief Nonparametric process model conditional to the existence of a value for kernel parameters */
-  // class ConditionalProcess
-  // {
-  // public:
-  //   ConditionalProcess(Dataset* data, bopt_params parameters){};
-  //   virtual ~ConditionalProcess(){};
 
-  //   /** 
-  //    * \brief Factory model generator for surrogate models
-  //    * @param parameters (process name, noise, priors, etc.)
-  //    * @return pointer to the corresponding derivate class (surrogate model)
-  //    */
-  //   static ConditionalProcess* create(Dataset* data, bopt_params parameters);
+  //// Inline methods
+  inline double Dataset::getSample(size_t index, vectord &x)
+  { 
+    x = mX[index];  
+    return mY(index);  
+  }
 
-  //   /** 
-  //    * \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.
-  //    */	
-  //   virtual ProbabilityDistribution* prediction(const vectord &query) = 0;
+  inline double Dataset::getLastSample(vectord &x)
+  { 
+    size_t last = mY.size()-1; 
+    return getSample(last, x); 
+  }
 
-  //   /** 
-  //    * \brief Computes the initial surrogate model. 
-  //    * It also updates the kernel parameters estimation. This
-  //    * function is hightly inefficient.  Use it only at the begining.
-  //    * @return error code
-  //    */
-  //   int fitSurrogateModel();
-  //   int precomputeSurrogate();
+  inline vectord Dataset::getPointAtMinimum() { return mX[mMinIndex]; };
+  inline double Dataset::getValueAtMinimum() { return mY(mMinIndex); };
+  inline size_t Dataset::getNSamples() { return mY.size(); };
+  inline void Dataset::checkBoundsY( size_t i )
+  {
+    if ( mY(mMinIndex) > mY(i) )       mMinIndex = i;
+    else if ( mY(mMaxIndex) < mY(i) )  mMaxIndex = i;
+  };
 
-  //   /** 
-  //    * \brief Sequential update of the surrogate model by adding a new point.
-  //    *  Add new point efficiently using Matrix Decomposition Lemma for the
-  //    *  inversion of the correlation matrix or adding new rows to the
-  //    *  Cholesky decomposition.  If retrain is true, it recomputes the
-  //    *  kernel hyperparameters and computes a full covariance matrix.
-  //    * @return error code
-  //    */   
-  //   int updateSurrogateModel(const vectord &Xnew, double Ynew);
-
-  // };
-
-    
 
   /**
    * \brief Abstract class to implement non-parametric processes
     // Getters and setters
     void setSamples(const matrixd &x, const vectord &y);
     void addSample(const vectord &x, double y);
-    Dataset* getData(){return &mData;}
+    Dataset* getData() {return &mData;}
 
-    inline vectord getPointAtMinimum() { return mData.getPointAtMinimum(); };
-    inline double getValueAtMinimum() { return mData.getValueAtMinimum(); };
-    inline double getSignalVariance() { return mSigma; };
+    vectord getPointAtMinimum() { return mData.getPointAtMinimum(); };
+    double getValueAtMinimum() { return mData.getValueAtMinimum(); };
+    double getSignalVariance() { return mSigma; };
 
     /** Sets the kind of learning methodology for kernel hyperparameters */
-    inline void setLearnType(learning_type l_type) { mLearnType = l_type; };
+    void setLearnType(learning_type l_type) { mLearnType = l_type; };
 
   protected:
     /** 
     const double mRegularizer;   ///< Std of the obs. model (also used as nugget)
   };
 
+  //// Inline methods
+  inline int NonParametricProcess::computeCorrMatrix(matrixd& corrMatrix)
+  {
+    return mKernel.computeCorrMatrix(mData.mX,corrMatrix,mRegularizer);
+  }
 
+  inline  matrixd NonParametricProcess::computeCorrMatrix()
+  {
+    const size_t nSamples = mData.getNSamples();
+    matrixd corrMatrix(nSamples,nSamples);
+    int error = mKernel.computeCorrMatrix(mData.mX,corrMatrix,mRegularizer);
+    return corrMatrix;
+  }
+
+  inline vectord NonParametricProcess::computeCrossCorrelation(const vectord &query)
+  {
+    return mKernel.computeCrossCorrelation(mData.mX,query);
+  }
+
+  inline  double NonParametricProcess::computeSelfCorrelation(const vectord& query)
+  {
+    return mKernel.computeSelfCorrelation(query);
+  }
 
   /**@}*/
   

include/optimization.hpp

     virtual double evaluate(const vectord& query) = 0;
   };
 
+  class DiscreteOptimization: public Optimization
+  {
+  public:
+    DiscreteOptimization(){};
+    DiscreteOptimization(vecOfvec *validSet): 
+      mInputSet(validSet){};
+    virtual ~DiscreteOptimization(){};
+
+    void setValidSet(vecOfvec* input)
+    { mInputSet = input; }
+
+    int run(vectord& result)
+    {
+      double current, min;
+  
+      result = *mInputSet->begin();
+      min = evaluate(result);
+  
+      for(vecOfvecIterator it = mInputSet->begin();
+	  it != mInputSet->end(); ++it)
+	{
+	  current = evaluate(*it);
+	  if (current < min)
+	    {
+	      result = *it;  
+	      min = current;
+	    }
+	}
+      return 0;
+    }
+
+  protected:
+    vecOfvec* mInputSet;               ///< List of input points
+  };
+
 }
 
 

include/optimizecriteria.hpp

     virtual ~OptimizeCriteria(){};
 
     double evaluate(const vectord& query)
-    {
-      return model_->evaluateCriteria(query);
-    }
+    {  return model_->evaluateCriteria(query);  }
     
   private:
     ContinuousModel* model_;

include/student_t_process_jef.hpp

 
   private:
     vectord mWML;           //!< GP ML parameters
-    
+    vectord mKn;
     /// Precomputed GP prediction operations
     vectord mAlphaF;
     matrixd mKF, mL2;
 
     StudentTDistribution* d_;      //!< Predictive distributions
+    clock_t timer;
   };
 
   /**@}*/

src/gaussian_process_ml.cpp

   {
     double kq = computeSelfCorrelation(query);
     vectord kn = computeCrossCorrelation(query);
-    vectord phi = mMean.getMeanFunc()->getFeatures(query);
+    vectord phi = mMean.getFeatures(query);
   
     vectord v(kn);
     inplace_solve(mL,v,ublas::lower_tag());

src/gaussian_process_normal.cpp

   {
     double kq = computeSelfCorrelation(query);
     vectord kn = computeCrossCorrelation(query);
-    vectord phi = mMean.getMeanFunc()->getFeatures(query);
+    vectord phi = mMean.getFeatures(query);
   
     vectord v(kn);
     inplace_solve(mL,v,ublas::lower_tag());

src/kernel_functors.cpp

       }
   }
 
-  Kernel* KernelModel::getKernel()
-  {
-    return mKernel.get();
-  }
-
   int KernelModel::setKernel (const vectord &thetav, 
 			      const vectord &stheta,
 			      std::string k_name, 
 			      size_t dim)
   {
+    KernelFactory mKFactory;
+
     mKernel.reset(mKFactory.create(k_name, dim));
     int error = setKernelPrior(thetav,stheta);
     
     return 0;
   }
 
-
-
-  vectord KernelModel::computeCrossCorrelation(const vecOfvec& XX, 
-					       const vectord &query)
-  {
-    vectord knx(XX.size());
-
-    std::vector<vectord>::const_iterator x_it  = XX.begin();
-    std::vector<vectord>::const_iterator x_end = XX.end();
-    vectord::iterator k_it = knx.begin();
-    while(x_it != x_end)
-      {
-	*k_it++ = (*mKernel)(*x_it++, query);
-      }
-    
-    return knx;
-  }
-
-  double KernelModel::computeSelfCorrelation(const vectord& query)
-  {
-    return (*mKernel)(query,query);
-  }
   
   double KernelModel::kernelLogPrior()
   {

src/mean_functors.cpp

       }
   }
 
-  ParametricFunction* MeanModel::getMeanFunc()
-  {
-    return mMean.get();
-  }
-
-
   int MeanModel::setMean (const vectord &muv,
 				     const vectord &smu,
 				     std::string m_name,
 				     size_t dim)
   {
+    MeanFactory mPFactory;
+
     mMean.reset(mPFactory.create(m_name,dim));
+
     if ("mZero" == m_name) 
       {
 	mMu = zvectord(1);

src/nonparametricprocess.cpp

 
 namespace bayesopt
 {
+
+  Dataset::Dataset(): mMinIndex(0), mMaxIndex(0) {};
+
+  Dataset::Dataset(const matrixd& x, const vectord& y):
+    mMinIndex(0), mMaxIndex(0)
+  {
+    setSamples(x,y);
+  };
+
+  Dataset::~Dataset(){};
+
+  void Dataset::setSamples(const matrixd &x, const vectord &y)
+  {
+    mY = y;
+    for (size_t i=0; i<x.size1(); ++i)
+      {
+	mX.push_back(row(x,i));
+	checkBoundsY(i);
+      } 
+  };
   
+  void Dataset::addSample(const vectord &x, double y)
+  {
+    mX.push_back(x);
+    mY.resize(mY.size()+1);  mY(mY.size()-1) = y;
+    checkBoundsY(mY.size()-1);
+  }
+
+  
+
+  ///////////////////////////////////////////////////////////////////////////
   NonParametricProcess::NonParametricProcess(size_t dim, bopt_params parameters):
     dim_(dim), mRegularizer(parameters.noise),
     mKernel(dim, parameters), mMean(dim, parameters)
     else  if(!name.compare("sGaussianProcessNormal"))
       s_ptr = new GaussianProcessNormal(dim,parameters);
     else if (!name.compare("sStudentTProcessJef"))
-      s_ptr = new StudentTProcessNIG(dim,parameters); 
+      s_ptr = new StudentTProcessJeffreys(dim,parameters); 
     else if (!name.compare("sStudentTProcessNIG"))
       s_ptr = new StudentTProcessNIG(dim,parameters); 
     else
   int NonParametricProcess::updateSurrogateModel( const vectord &Xnew,
 						  double Ynew, bool retrain)
   {
-    assert( Dataset.mX[1].size() == Xnew.size() );
+    assert( mData.mX[1].size() == Xnew.size() );
 
     if (retrain)
       {
     return utils::cholesky_decompose(K,mL);
   }
 
-
-  int NonParametricProcess::computeCorrMatrix(matrixd& corrMatrix)
-  {
-    return mKernel.computeCorrMatrix(mData.mX,corrMatrix,mRegularizer);
-  }
-
-
-
-  matrixd NonParametricProcess::computeCorrMatrix()
-  {
-    const size_t nSamples = mData.getNSamples();
-    matrixd corrMatrix(nSamples,nSamples);
-    int error = mKernel.computeCorrMatrix(mData.mX,corrMatrix,mRegularizer);
-    return corrMatrix;
-  }
-
   matrixd NonParametricProcess::computeDerivativeCorrMatrix(int dth_index)
   {
     const size_t nSamples = mData.getNSamples();
     return corrMatrix;
   }
 
-
-  vectord NonParametricProcess::computeCrossCorrelation(const vectord &query)
-  {
-    return mKernel.computeCrossCorrelation(mData.mX,query);
-  }
-
-  double NonParametricProcess::computeSelfCorrelation(const vectord& query)
-  {
-    return mKernel.computeSelfCorrelation(query);
-  }
-
 } //namespace bayesopt

src/student_t_process_jef.cpp

 -----------------------------------------------------------------------------
 */
 
-
+#include <ctime>
 #include "cholesky.hpp"
 #include "trace_ublas.hpp"
 #include "student_t_process_jef.hpp"
 						   bopt_params params):
     HierarchicalGaussianProcess(dim, params)
   {
+    timer = 0;
     mSigma = params.sigma_s;
     d_ = new StudentTDistribution();
   }  // Constructor
 
   StudentTProcessJeffreys::~StudentTProcessJeffreys()
   {
+    std::cout << "Time in pred:" << (double)(timer) / (double)CLOCKS_PER_SEC << std::endl;
     delete d_;
   } // Default destructor
 
   ProbabilityDistribution* 
   StudentTProcessJeffreys::prediction(const vectord &query )
   {
+    clock_t start = clock();
     double kq = computeSelfCorrelation(query);
-    vectord kn = computeCrossCorrelation(query);
-    vectord phi = mMean.getMeanFunc()->getFeatures(query);
+    //    vectord kn = computeCrossCorrelation(query);
+    mKernel.computeCrossCorrelation(mData.mX,query,mKn);
+    vectord phi = mMean.getFeatures(query);
   
-    vectord v(kn);
-    inplace_solve(mL,v,ublas::lower_tag());
+    //    vectord v(mKn);
+    inplace_solve(mL,mKn,ublas::lower_tag());
 
-    vectord rq = phi - prod(v,mKF);
+    vectord rho = phi - prod(mKn,mKF);
 
-    vectord rho(rq);
+    //    vectord rho(rq);
     inplace_solve(mL2,rho,ublas::lower_tag());
     
-    double yPred = inner_prod(phi,mWML) + inner_prod(v,mAlphaF);
-    double sPred = sqrt( mSigma * (kq - inner_prod(v,v) 
+    double yPred = inner_prod(phi,mWML) + inner_prod(mKn,mAlphaF);
+    double sPred = sqrt( mSigma * (kq - inner_prod(mKn,mKn) 
 				   + inner_prod(rho,rho)));
 
     d_->setMeanAndStd(yPred,sPred);
+    timer += clock() - start;
     return d_;
   }
 
   int StudentTProcessJeffreys::precomputePrediction()
   {
     size_t n = mData.getNSamples();
-    size_t p = mMean.getMeanFunc()->nFeatures();
+    size_t p = mMean.nFeatures();
+
+    mKn.resize(n);
 
     mKF = trans(mMean.mFeatM);
     inplace_solve(mL,mKF,ublas::lower_tag());

src/student_t_process_nig.cpp

   {
     double kq = computeSelfCorrelation(query);
     vectord kn = computeCrossCorrelation(query);
-    vectord phi = mMean.getMeanFunc()->getFeatures(query);
+    vectord phi = mMean.getFeatures(query);
   
     vectord v(kn);
     inplace_solve(mL,v,ublas::lower_tag());