Commits

Ruben Martinez-Cantin committed 3cbfe02

Adding a new interface for nonparametric process

Comments (0)

Files changed (30)

 
 INCLUDE(UseDoxygen)
 
+# Sobol sequences are hardcoded tables, so it might take a lot of time
+# to compile.
 IF(BAYESOPT_BUILD_SOBOL)
   ADD_DEFINITIONS(-DUSE_SOBOL)	
   SET(SOBOL_SRC
   ./src/bayesoptdisc.cpp
   ./src/bayesoptbase.cpp
   ./src/inneroptimization.cpp
+  ./src/bayesianregressor.cpp
   ./src/nonparametricprocess.cpp
   ./src/empiricalbayesprocess.cpp
   ./src/hierarchical_gaussian_process.cpp
 LINK_DIRECTORIES( ${CMAKE_SOURCE_DIR}/lib )
 
 IF(NLOPT_BUILD)
-  ADD_SUBDIRECTORY(nlopt2)
-  include_directories(${CMAKE_SOURCE_DIR}/nlopt2/api)
+  ADD_SUBDIRECTORY(nlopt)
+  include_directories(${CMAKE_SOURCE_DIR}/nlopt/api)
   SET(EXT_LIBS nlopt)
 ELSE(NLOPT_BUILD)
   SET(EXT_LIBS ${NLOPT})

include/bayesoptbase.hpp

 #ifndef  _BAYESOPTBASE_HPP_
 #define  _BAYESOPTBASE_HPP_
 
+#include <boost/scoped_ptr.hpp>
 #include "criteria_functors.hpp"
 
 /**
     };
 
     
-    inline NonParametricProcess* getSurrogateModel()
+    inline BayesianRegressor* getSurrogateModel()
     { return mGP.get(); };
 
     int setSurrogateModel();    
 
 
   protected:
-    boost::scoped_ptr<NonParametricProcess> mGP;  ///< Pointer to surrogate model
+    boost::scoped_ptr<BayesianRegressor> mGP;  ///< Pointer to surrogate model
     boost::scoped_ptr<Criteria> mCrit;                    ///< Metacriteria model
     bopt_params mParameters;                        ///< Configuration parameters
     size_t mDims;                                       ///< Number of dimensions

include/criteria_atomic.hpp

   {
   public:
     virtual ~AtomicCriteria(){};
-    virtual int init(NonParametricProcess* proc)
+    virtual int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       return 0;
   {
   public:
     virtual ~ExpectedImprovement(){};
-    int init(NonParametricProcess* proc)
+    int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       mExp = 1;
   {
   public:
     virtual ~BiasedExpectedImprovement(){};
-    int init(NonParametricProcess* proc)
+    int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       mBias = 0.01;
   {
   public:
     virtual ~LowerConfidenceBound(){};
-    int init(NonParametricProcess* proc)
+    int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       mBeta = 1.0;
   {
   public:
     virtual ~ProbabilityOfImprovement(){};
-    int init(NonParametricProcess* proc)
+    int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       mEpsilon = 0.01;
   {
   public:
     virtual ~AnnealedExpectedImprovement(){};
-    int init(NonParametricProcess* proc)
+    int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       reset();
   {
   public:
     virtual ~AnnealedLowerConfindenceBound(){};
-    int init(NonParametricProcess* proc)
+    int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       reset();
   class InputDistance: public AtomicCriteria
   {
   public:
-    int init(NonParametricProcess* proc)
+    int init(BayesianRegressor* proc)
     { 
       mProc = proc;
       mW = 1;

include/criteria_combined.hpp

 	delete mCriteriaList[i];
       }
     };
-    virtual int init(NonParametricProcess *proc, 
+    virtual int init(BayesianRegressor *proc, 
 		     const std::vector<Criteria*>& list) 
     { 
       mProc = proc;
 
 
   protected:
-    NonParametricProcess* mProc;
+    BayesianRegressor* mProc;
     std::vector<Criteria*> mCriteriaList;
   };
 
   public:
     GP_Hedge();
     virtual ~GP_Hedge() {};
-    int init(NonParametricProcess *proc, 
+    int init(BayesianRegressor *proc, 
 	     const std::vector<Criteria*>& list);
 
     bool requireComparison(){ return true; };

include/criteria_functors.hpp

 #ifndef  _CRITERIA_FUNCTORS_HPP_
 #define  _CRITERIA_FUNCTORS_HPP_
 
+#include <map>
 #include <algorithm>
 #include "optimizable.hpp"
-#include "nonparametricprocess.hpp"
+#include "bayesianregressor.hpp"
 
 namespace bayesopt
 {
   {
   public:
     virtual ~Criteria() {};
-    virtual int init(NonParametricProcess *proc) { return 0; };
-    virtual int init(NonParametricProcess *proc, 
+    virtual int init(BayesianRegressor *proc) { return 0; };
+    virtual int init(BayesianRegressor *proc, 
 		     const std::vector<Criteria*>& list) { return 0; };
 
     virtual bool requireComparison() = 0;
 			     std::string& name,
 			     int& error_code) { assert(false); return false; };
   protected:
-    NonParametricProcess *mProc;
+    BayesianRegressor *mProc;
   };
 
 
     CriteriaFactory ();
     virtual ~CriteriaFactory () {};
   
-    //Criteria* create(criterium_name name, NonParametricProcess* proc);
-    Criteria* create(std::string name, NonParametricProcess* proc);
+    //Criteria* create(criterium_name name, BayesianRegressor* proc);
+    Criteria* create(std::string name, BayesianRegressor* proc);
     
   private:
     typedef Criteria* (*create_func_definition)();

include/empiricalbayesprocess.hpp

     /** 
      * \brief Updates the kernel parameters acording with a point
      * estimate (ML, MAP, etc.)
-     * @return error code
      */
-    int updateKernelParameters();
+    void updateKernelParameters();
 
     /** 
-     * \brief Updates and computes the score (eg:likelihood) of the kernel
-     * parameters.
-     * @param query set of parameters.
+     * \brief Computes the score (eg:likelihood) of the kernel
+     * parameters.  
+     * Warning: To evaluate the score, it is necessary to change the parameters
+     * @param x set of parameters.  
      * @return score
      */
-    double evaluateKernelParams(const vectord& query);
-    double evaluate(const vectord &x) {return evaluateKernelParams(x);}
+    double evaluate(const vectord &x);
 
     /** 
      * \brief Computes the score (eg:likelihood) of the current kernel
     NLOPT_Optimization* kOptimizer;
   };
 
+
+
+  inline double EmpiricalBayesProcess::evaluate(const vectord& x)
+  { 
+    mKernel.setHyperParameters(x);
+    return evaluateKernelParams();
+  };
+
+
   /**@}*/
   
 } //namespace bayesopt

include/fullbayesprocess.hpp

     ProbabilityDistribution* prediction(const vectord &query);
 
     /** 
-     * \brief Updates the kernel parameters acording with a point
-     * estimate (ML, MAP, etc.)
-     * @return error code
+     * \brief Updates the kernel parameters acording with a Bayesian
+     * estimate (grid sampling, MCMC, etc.)
      */
-    int updateKernelParameters();
+    void updateKernelParameters();
 
   private:
     std::vector<NonParametricProcess*>   mVProc;

include/gaussian_process.hpp

      */
     double negativeLogLikelihood();
 
-    /** 
-     * \brief Precompute some values of the prediction that do not depends on
-     * the query
-     * @return error code
+    /** Precompute some values of the prediction that do not depends
+     *	on the query
      */
-    int precomputePrediction();
+    void precomputePrediction();
 
   private:
     vectord mAlphaV;              ///< Precomputed L\y

include/gaussian_process_ml.hpp

      */
     double negativeLogLikelihood();
 
-    /** 
-     * \brief Precompute some values of the prediction that do not depends on
+    /** Precompute some values of the prediction that do not depends on
      * the query
-     * @return error code
      */
-    int precomputePrediction();
+    void precomputePrediction();
 
   private:
     vectord mWML;           //!< GP ML parameters

include/gaussian_process_normal.hpp

      */
     double negativeLogLikelihood();
 
-    /** 
-     * \brief Precompute some values of the prediction that do not depends on
-     * the query
-     * @return error code
+    /** Precompute some values of the prediction that do not depends
+     *	on the query
      */
-    int precomputePrediction();
+    void precomputePrediction();
 
   private:
     vectord mWMap;                      //!< GP posterior parameters

include/kernel_atomic.hpp

       n_inputs = input_dim;
       return 0;
     };
-    int setHyperParameters(const vectord &theta) 
+    void setHyperParameters(const vectord &theta) 
     {
       if(theta.size() != n_params)
 	{
 	  FILE_LOG(logERROR) << "Wrong number of kernel hyperparameters"; 
-	  return -1; 
+	  throw std::invalid_argument("Wrong number of kernel hyperparameters");
 	}
       params = theta;
-      return 0;
     };
+
     vectord getHyperParameters() {return params;};
     size_t nHyperParameters() {return n_params;};
 

include/kernel_combined.hpp

       this->right = right;
       return 0;
     };
-    int setHyperParameters(const vectord &theta) 
+    void setHyperParameters(const vectord &theta) 
     {
       using boost::numeric::ublas::subrange;
-      int error = 0;
 
       size_t n_lhs = left->nHyperParameters();
       size_t n_rhs = right->nHyperParameters();
       if (theta.size() != n_lhs + n_rhs)
 	{
 	  FILE_LOG(logERROR) << "Wrong number of kernel hyperparameters"; 
-	  return -1; 
+	  throw std::invalid_argument("Wrong number of kernel hyperparameters");
 	}
-      error += left->setHyperParameters(subrange(theta,0,n_lhs));
-      error += right->setHyperParameters(subrange(theta,n_lhs,n_lhs+n_rhs));
-      return error;
+      left->setHyperParameters(subrange(theta,0,n_lhs));
+      right->setHyperParameters(subrange(theta,n_lhs,n_lhs+n_rhs));
     };
 
     vectord getHyperParameters() 

include/kernel_functors.hpp

     virtual int init(size_t input_dim) {return 0;};
     virtual int init(size_t input_dim, Kernel* left, Kernel* right) {return 0;};
 
-    virtual int setHyperParameters(const vectord &theta) = 0;
+    virtual void setHyperParameters(const vectord &theta) = 0;
     virtual vectord getHyperParameters() = 0;
     virtual size_t nHyperParameters() = 0;
 
 
     Kernel* getKernel();
     
-    int setHyperParameters(const vectord &theta);
+    void setHyperParameters(const vectord &theta);
     vectord getHyperParameters();
     size_t nHyperParameters();
 
      * @param thetav kernel parameters (mean)
      * @param stheta kernel parameters (std)
      * @param k_name kernel name
-     * @return error_code
      */
-    int setKernel (const vectord &thetav, const vectord &stheta, 
+    void setKernel (const vectord &thetav, const vectord &stheta, 
 		   std::string k_name, size_t dim);
 
     /** Wrapper of setKernel for C kernel structure */
-    int setKernel (kernel_parameters kernel, size_t dim);
+    void setKernel (kernel_parameters kernel, size_t dim);
 
-    int computeCorrMatrix(const vecOfvec& XX, matrixd& corrMatrix, double nugget);
-    int computeDerivativeCorrMatrix(const vecOfvec& XX, matrixd& corrMatrix, int dth_index);
+    void computeCorrMatrix(const vecOfvec& XX, matrixd& corrMatrix, double nugget);
+    void 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,
+    void computeCrossCorrelation(const vecOfvec& XX, const vectord &query,
 				 vectord& knx);
     double computeSelfCorrelation(const vectord& query);
     double kernelLogPrior();
 
   private:
     /** Set prior (Gaussian) for kernel hyperparameters */
-    int setKernelPrior (const vectord &theta, const vectord &s_theta);
+    void setKernelPrior (const vectord &theta, const vectord &s_theta);
 
     boost::scoped_ptr<Kernel> mKernel;            ///< Pointer to kernel function
     std::vector<boost::math::normal> priorKernel; ///< Prior of kernel parameters
   inline Kernel* KernelModel::getKernel()
   { return mKernel.get();  }
 
-  inline int KernelModel::setHyperParameters(const vectord &theta)
-  { return mKernel->setHyperParameters(theta); };
+  inline void KernelModel::setHyperParameters(const vectord &theta)
+  { mKernel->setHyperParameters(theta); };
     
   inline vectord KernelModel::getHyperParameters()
   {return mKernel->getHyperParameters();};
   inline double KernelModel::computeSelfCorrelation(const vectord& query)
   { return (*mKernel)(query,query); }
 
+  inline void KernelModel::setKernelPrior (const vectord &theta, 
+					   const vectord &s_theta)
+  {
+    for (size_t i = 0; i<theta.size(); ++i)
+      {
+	boost::math::normal n(theta(i),s_theta(i));
+	priorKernel.push_back(n);
+      }
+  };
+
+
   //@}
 
 } //namespace bayesopt

include/nonparametricprocess.hpp

 #ifndef __NONPARAMETRICPROCESS_HPP__
 #define __NONPARAMETRICPROCESS_HPP__
 
-#include "ublas_extra.hpp"
+#include "cholesky.hpp"
+#include "bayesianregressor.hpp"
 #include "kernel_functors.hpp"
-#include "mean_functors.hpp"
-#include "prob_distribution.hpp"
 
 namespace bayesopt
 {
    */
   /**@{*/
 
-  /** \brief Dataset model to deal with the vector (real) based datasets */
-  class Dataset
-  {
-  public:
-    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);
-
-    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;	
-  };
-
-
-  //// Inline methods
-  inline double Dataset::getSample(size_t index, vectord &x)
-  { 
-    x = mX[index];  
-    return mY(index);  
-  }
-
-  inline double Dataset::getLastSample(vectord &x)
-  { 
-    size_t last = mY.size()-1; 
-    return getSample(last, x); 
-  }
-
-  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 Abstract class to implement non-parametric processes
    */
-  class BAYESOPT_API NonParametricProcess
+  class NonParametricProcess: public BayesianRegressor
   {
   public:
     NonParametricProcess(size_t dim, bopt_params parameters);
     virtual ~NonParametricProcess();
 
     /** 
-     * \brief Factory model generator for surrogate models
-     * @param parameters (process name, noise, priors, etc.)
-     * @return pointer to the corresponding derivate class (surrogate model)
+     * \brief Updates the kernel parameters with a point estimate
+     * (empirical Bayes) or a full Bayesian distribution
      */
-    static NonParametricProcess* create(size_t dim, bopt_params parameters);
+    virtual void updateKernelParameters() = 0;
 
     /** 
-     * \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;
-		 		 
-    /** 
-     * \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
+     * \brief Computes the initial surrogate model and updates the
+     * kernel parameters estimation. 
+     *
+     * This function requires to recompute all covariance matrixes,
+     * inverses, etc.  Use it with precaution.
      */
-    int fitSurrogateModel();
-    virtual int updateKernelParameters() = 0;
-    int precomputeSurrogate();
+    void fitSurrogateModel();
 
     /** 
      * \brief Sequential update of the surrogate model by adding a new point.
      *  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, bool retrain);
+    void updateSurrogateModel(const vectord &Xnew, double Ynew, bool retrain);
 
 
     // Getters and setters
-    void setSamples(const matrixd &x, const vectord &y);
-    void addSample(const vectord &x, double y);
-    Dataset* getData() {return &mData;}
-
-    vectord getPointAtMinimum() { return mData.getPointAtMinimum(); };
-    double getValueAtMinimum() { return mData.getValueAtMinimum(); };
-    double getSignalVariance() { return mSigma; };
+    double getSignalVariance();
 
     /** Sets the kind of learning methodology for kernel hyperparameters */
-    void setLearnType(learning_type l_type) { mLearnType = l_type; };
+    void setLearnType(learning_type l_type);
 
   protected:
-    /** 
-     * \brief Precompute some values of the prediction that do not
-     * depends on the query.
-     * @return error code
-     */
-    virtual int precomputePrediction() = 0;
 
-    //TODO: Refactorize from KernelModel
-    int computeCorrMatrix(matrixd& corrMatrix);
+    /** Precompute some values of the prediction that do not depends
+     * on the query. */
+    virtual void precomputePrediction() = 0;
+
+    /** Computes the Correlation (Kernel or Gram) matrix */
+    void computeCorrMatrix(matrixd& corrMatrix);
+
+    /** Computes the Correlation (Kernel or Gram) matrix */
     matrixd computeCorrMatrix();
+
+    /** Computes the derivative of the correlation matrix with respect
+     *	to the dth hyperparameter */
     matrixd computeDerivativeCorrMatrix(int dth_index);
     vectord computeCrossCorrelation(const vectord &query);
     double computeSelfCorrelation(const vectord& query);
 
+    /** Computes the Cholesky decomposition of the Correlation matrix */
+    void computeCholeskyCorrelation();
+
 
   protected:
-    Dataset mData;
-    
-    double mSigma;                                   //!< GP posterior parameters
     matrixd mL;             ///< Cholesky decomposition of the Correlation matrix
-    size_t dim_;
     learning_type mLearnType;
     KernelModel mKernel;
-    MeanModel mMean;
 
   private:
-    /** 
-     * Computes the Cholesky decomposition of the Correlation (Kernel
-     * or Gram) matrix 
-     * @return error code
-     */
-    int computeCholeskyCorrelation();
-
-    /** 
-     * Adds a new point to the Cholesky decomposition of the Correlation 
-     * matrix.
-     * @return error code
-     */
-    int addNewPointToCholesky(const vectord& correlation,
+    /** Adds a new point to the Cholesky decomposition of the
+     * Correlation matrix. */
+    void addNewPointToCholesky(const vectord& correlation,
 			      double selfcorrelation);
 
 
   };
 
   //// Inline methods
-  inline int NonParametricProcess::computeCorrMatrix(matrixd& corrMatrix)
+  inline void NonParametricProcess::fitSurrogateModel()
   {
-    return mKernel.computeCorrMatrix(mData.mX,corrMatrix,mRegularizer);
-  }
+    updateKernelParameters();
+    computeCholeskyCorrelation();
+    precomputePrediction(); 
+  };
+
+  inline void NonParametricProcess::setLearnType(learning_type l_type) 
+  { mLearnType = l_type; };
+
+  inline void NonParametricProcess::computeCorrMatrix(matrixd& corrMatrix)
+  { 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);
+    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); }
+
+  inline void NonParametricProcess::addNewPointToCholesky(const vectord& correlation,
+							  double selfcorrelation)
   {
-    return mKernel.computeCrossCorrelation(mData.mX,query);
+    vectord newK(correlation);
+    utils::append(newK, selfcorrelation);
+    utils::cholesky_add_row(mL,newK);
   }
 
-  inline  double NonParametricProcess::computeSelfCorrelation(const vectord& query)
-  {
-    return mKernel.computeSelfCorrelation(query);
-  }
 
   /**@}*/
   

include/student_t_process_jef.hpp

      */
     double negativeLogLikelihood();
 
-    /** 
-     * \brief Precompute some values of the prediction that do not depends on
-     * the query
-     * @return error code
+    /** Precompute some values of the prediction that do not depends
+     *	on the query
      */
-    int precomputePrediction();
+    void precomputePrediction();
 
   private:
     vectord mWML;           //!< GP ML parameters

include/student_t_process_nig.hpp

      */
     double negativeLogLikelihood();
 
-    /** 
-     * \brief Precompute some values of the prediction that do not depends on
-     * the query
-     * @return error code
+    /** Precompute some values of the prediction that do not depends
+     *	on the query
      */
-    int precomputePrediction();
+    void precomputePrediction();
 
   private:
     vectord mWMap;                      //!< GP posterior parameters

src/bayesoptbase.cpp

   int BayesOptBase::setSurrogateModel()
   {
     // Configure Surrogate and Criteria Functions
-    mGP.reset(NonParametricProcess::create(mDims,mParameters));
+    mGP.reset(BayesianRegressor::create(mDims,mParameters));
     if (mGP == NULL) 
       {
 	FILE_LOG(logERROR) << "Error setting the surrogate function"; 

src/criteria_combined.cpp

     sampleUniform( mtRandom, realUniformDist(0,1))
   {};
 
-  int GP_Hedge::init(NonParametricProcess *proc, 
+  int GP_Hedge::init(BayesianRegressor *proc, 
 		     const std::vector<Criteria*>& list) 
   { 
     mProc = proc;

src/criteria_functors.cpp

    * @return criteria pointer
    */
   Criteria* CriteriaFactory::create(std::string name,
-				    NonParametricProcess* proc)
+				    BayesianRegressor* proc)
   {
     Criteria *cFunc;
     std::string os;

src/empiricalbayesprocess.cpp

   }
 
 
-  int EmpiricalBayesProcess::updateKernelParameters()
+  void EmpiricalBayesProcess::updateKernelParameters()
   {
     if (mLearnType == L_FIXED)
       {
 	FILE_LOG(logDEBUG) << "Fixed hyperparameters. Not learning";
-	return 0;
       }
     else
       {
-	int error = -1;
 	vectord optimalTheta = mKernel.getHyperParameters();
 	
-	FILE_LOG(logDEBUG) << "Computing kernel parameters. Initial: " 
-			   << optimalTheta;
-
+	FILE_LOG(logDEBUG) << "Initial kernel parameters: " << optimalTheta;
 	kOptimizer->run(optimalTheta);
-	error = mKernel.setHyperParameters(optimalTheta);
-
-	if (error)
-	  {
-	    FILE_LOG(logERROR) << "Error updating kernel parameters.";
-	    exit(EXIT_FAILURE);
-	  }   
-
+	mKernel.setHyperParameters(optimalTheta);
 	FILE_LOG(logDEBUG) << "Final kernel parameters: " << optimalTheta;	
-	return error;
       }
   };
 
-  double EmpiricalBayesProcess::evaluateKernelParams(const vectord& query)
-  { 
-    if (mLearnType == L_FIXED)
-      {
-	FILE_LOG(logERROR) << "Fixed hyperparameters should not change.";
-	return -1;
-      }
-    else
-      {
-	int error = mKernel.setHyperParameters(query);
-	if (error) 
-	  {
-	    FILE_LOG(logERROR) << "Problem optimizing kernel parameters."; 
-	    exit(EXIT_FAILURE);	
-	  }
-	return evaluateKernelParams();
-      }
-  };
-  
   double EmpiricalBayesProcess::evaluateKernelParams()
   { 
     double result;
 	utils::erase_column(mMean.mFeatM,0);
 
 	// Compute the cross validation
-	precomputeSurrogate();
+	computeCholeskyCorrelation();
+	precomputePrediction(); 
 	ProbabilityDistribution* pd = prediction(x);
 	sum += log(pd->pdf(y));
 

src/fullbayesprocess.cpp

       }
   }
 
-  int FullBayesProcess::precomputePrediction()
+  void FullBayesProcess::precomputePrediction()
   {
     updateKernelParameters()
     for(size_t i=0;i<N_PROC;++i)
       { 
 	vectord th = column(kTheta,ii);
 	std::copy(th.begin(),th.end(),newParams.kernel.hp_mean);
-	mVProc.push_back(NonParametricProcess::create(dim_,newParams));
+	mVProc.push_back(BayesianRegressor::create(dim_,newParams));
       }
     
     return 0;
   }
 
-  int FullBayesProcess::updateKernelParameters()
+  void FullBayesProcess::updateKernelParameters()
   {
     double sum = 0.0;
     for (size_t ii = 0; ii < N_PROC; ++ii)

src/gaussian_process.cpp

   }
 
 
-  int GaussianProcess::precomputePrediction()
+  void GaussianProcess::precomputePrediction()
   {
     const size_t n = mData.getNSamples();
   
     mAlphaV.resize(n,false);
     mAlphaV = mData.mY-mMean.muTimesFeat();
     inplace_solve(mL,mAlphaV,ublas::lower_tag());
-
-    return 0; 
   }
 	
 } //namespace bayesopt

src/gaussian_process_ml.cpp

     return d_;
   }
 
-  int GaussianProcessML::precomputePrediction()
+  void GaussianProcessML::precomputePrediction()
   {
     size_t n = mData.getNSamples();
     size_t p = mMean.getMeanFunc()->nFeatures();
     mAlphaF = mData.mY - prod(mWML,mMean.mFeatM);
     inplace_solve(mL,mAlphaF,ublas::lower_tag());
     mSigma = inner_prod(mAlphaF,mAlphaF)/(n-p);
-  
-    return 0;
   }
 
 } //namespace bayesopt

src/gaussian_process_normal.cpp

 
 
 
-  int GaussianProcessNormal::precomputePrediction()
+  void GaussianProcessNormal::precomputePrediction()
   {
     size_t n = mData.getNSamples();
     size_t p = mMean.getMeanFunc()->nFeatures();
     mVf = mData.mY - prod(trans(mMean.mFeatM),mWMap);
     inplace_solve(mL,mVf,ublas::lower_tag());
 
-    if (boost::math::isnan(mWMap(0)))
+    if ((boost::math::isnan(mWMap(0))) || (boost::math::isnan(mSigma)))
       {
 	FILE_LOG(logERROR) << "Error in precomputed prediction. NaN found.";
-	return -1;
+	throw std::runtime_error("Error in precomputed prediction. NaN found.");
       }
-    return 0;
   }
 
 } //namespace bayesopt

src/kernel_functors.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 <stdexcept>
 #include "log.hpp"
 #include "parser.hpp"
 #include "ublas_extra.hpp"
     if (it == registry.end()) 
       {
 	FILE_LOG(logERROR) << "Error: Fatal error while parsing "
-			   << "kernel function: " << os 
-			   << " not found" << std::endl;
+			   << "kernel function: " << os << " not found";
+	throw std::invalid_argument("Kernel not found " + os);
 	return NULL;
       } 
     kFunc = it->second();
       {
 	kFunc->init(input_dim);
       } 
-    else 
+    else // Combined kernel
       {
 	kFunc->init(input_dim, create(os1,input_dim), create(os2,input_dim));
       }
     return kFunc;
-
   };
 
 
   //////////////////////////////////////////////////////////////////////
 
   KernelModel::KernelModel(size_t dim, bopt_params parameters)
-  {
-    int errorK = setKernel(parameters.kernel,dim);
-    if (errorK)
-      {
-	FILE_LOG(logERROR) << "Error initializing nonparametric process.";
-	exit(EXIT_FAILURE);
-      }
-  }
+  { setKernel(parameters.kernel,dim);  }
 
-  int KernelModel::setKernel (const vectord &thetav, 
+  void 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);
-    
-    if (mKernel == NULL || error)   return -1;
-
-    return mKernel->setHyperParameters(thetav);
+    setKernelPrior(thetav,stheta);
+    mKernel->setHyperParameters(thetav);
   }
 
-  int KernelModel::setKernel (kernel_parameters kernel, 
+  void KernelModel::setKernel (kernel_parameters kernel, 
 			      size_t dim)
   {
     size_t n = kernel.n_hp;
     vectord th = utils::array2vector(kernel.hp_mean,n);
     vectord sth = utils::array2vector(kernel.hp_std,n);
-    int error = setKernel(th, sth, kernel.name, dim);
-    return 0;
+    setKernel(th, sth, kernel.name, dim);
   };
 
-  int KernelModel::setKernelPrior (const vectord &theta, 
-				   const vectord &s_theta)
-  {
-    size_t n_theta = theta.size();
-    for (size_t i = 0; i<n_theta; ++i)
-      {
-	boost::math::normal n(theta(i),s_theta(i));
-	priorKernel.push_back(n);
-      }
-    return 0;
-  };
 
-  int KernelModel::computeCorrMatrix(const vecOfvec& XX, matrixd& corrMatrix, 
+  void KernelModel::computeCorrMatrix(const vecOfvec& XX, matrixd& corrMatrix, 
 				     double nugget)
   {
     assert(corrMatrix.size1() == XX.size());
 	  }
 	corrMatrix(ii,ii) = (*mKernel)(XX[ii],XX[ii]) + nugget;
       }
-    return 0;
   }
 
-  int KernelModel::computeDerivativeCorrMatrix(const vecOfvec& XX, 
+  void KernelModel::computeDerivativeCorrMatrix(const vecOfvec& XX, 
 					       matrixd& corrMatrix,
 					       int dth_index)
   {
 	  }
 	corrMatrix(ii,ii) = mKernel->gradient(XX[ii],XX[ii],dth_index);
       }
-    return 0;
   }
 
   

src/nonparametricprocess.cpp

-
 /*
 -------------------------------------------------------------------------
    This file is part of BayesOpt, an efficient C++ library for 
 
 #include <cstdio>
 #include <cstdlib>
+#include <stdexcept>
+#include <boost/lexical_cast.hpp>
+
 #include "nonparametricprocess.hpp"
+
 #include "log.hpp"
-#include "cholesky.hpp"
 #include "ublas_extra.hpp"
 
-#include "gaussian_process.hpp"
-#include "gaussian_process_ml.hpp"
-#include "gaussian_process_normal.hpp"
-#include "student_t_process_jef.hpp"
-#include "student_t_process_nig.hpp"
-
 
 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)
+    BayesianRegressor(dim,parameters), mRegularizer(parameters.noise),
+    mKernel(dim, parameters)
   { 
     setLearnType(parameters.l_type);
   }
 
-  NonParametricProcess::~NonParametricProcess()
-  {}
+  NonParametricProcess::~NonParametricProcess(){}
 
 
-  NonParametricProcess* NonParametricProcess::create(size_t dim, 
-						     bopt_params parameters)
-  {
-    NonParametricProcess* s_ptr;
 
-    std::string name = parameters.surr_name;
-
-    if (!name.compare("sGaussianProcess"))
-      s_ptr = new GaussianProcess(dim,parameters);
-    else  if(!name.compare("sGaussianProcessML"))
-      s_ptr = new GaussianProcessML(dim,parameters);
-    else  if(!name.compare("sGaussianProcessNormal"))
-      s_ptr = new GaussianProcessNormal(dim,parameters);
-    else if (!name.compare("sStudentTProcessJef"))
-      s_ptr = new StudentTProcessJeffreys(dim,parameters); 
-    else if (!name.compare("sStudentTProcessNIG"))
-      s_ptr = new StudentTProcessNIG(dim,parameters); 
-    else
-      {
-	FILE_LOG(logERROR) << "Error: surrogate function not supported.";
-	return NULL;
-      }
-    return s_ptr;
-  };
-
-
-  int NonParametricProcess::fitSurrogateModel()
-  {
-    int error = updateKernelParameters();
-    error += precomputeSurrogate();
-    return error;
-  };
-
-
-  int NonParametricProcess::precomputeSurrogate()
-  {
-    int error = -1;
-    error = computeCholeskyCorrelation();
-    if (error)
-      {
-	FILE_LOG(logERROR) << "Error computing the correlation matrix";
-	exit(EXIT_FAILURE);
-      }   
-
-    error = precomputePrediction(); 
-    if (error)
-      {
-	FILE_LOG(logERROR) << "Error pre-computing the prediction distribution";
-	exit(EXIT_FAILURE);
-      }   
-
-    return error; 
-  } // fitSurrogateModel
-
-
-  int NonParametricProcess::updateSurrogateModel( const vectord &Xnew,
+  void NonParametricProcess::updateSurrogateModel( const vectord &Xnew,
 						  double Ynew, bool retrain)
   {
     assert( mData.mX[1].size() == Xnew.size() );
 
     if (retrain)
       {
+	FILE_LOG(logDEBUG) << "Retraining model parameters";
 	addSample(Xnew,Ynew);
-	FILE_LOG(logDEBUG) << "Retraining model parameters";
-	return fitSurrogateModel();	
+	fitSurrogateModel();	
       }
     else
       {
 	double selfCorrelation = computeSelfCorrelation(Xnew) + mRegularizer;
 	addSample(Xnew,Ynew);
 	addNewPointToCholesky(newK,selfCorrelation);
-
-	int error = precomputePrediction(); 
-	if (error < 0)
-	  {
-	    FILE_LOG(logERROR) << "Error pre-computing the prediction distribution";
-	    exit(EXIT_FAILURE);
-	  }   
-	return error; 
+	precomputePrediction(); 
       }
-    return 0; //JIC
   } // updateSurrogateModel
 
 
-  //////////////////////////////////////////////////////////////////////////////
-  //// Getters and Setters
-  void NonParametricProcess::setSamples(const matrixd &x, const vectord &y)
-  {
-    mData.setSamples(x,y);
-    mMean.setPoints(mData.mX);  //Because it expects a vecOfvec instead of a matrixd
-  }
-
-  void NonParametricProcess::addSample(const vectord &x, double y)
-  {
-    mData.addSample(x,y);
-    mMean.addNewPoint(x);
-  };
-
-  
-  int NonParametricProcess::addNewPointToCholesky(const vectord& correlation,
-						  double selfcorrelation)
-  {
-    vectord newK(correlation);
-    utils::append(newK, selfcorrelation);
-    utils::cholesky_add_row(mL,newK);
-    return 0;
-  }
-
-
-  int NonParametricProcess::computeCholeskyCorrelation()
+  void NonParametricProcess::computeCholeskyCorrelation()
   {
     size_t nSamples = mData.getNSamples();
     mL.resize(nSamples,nSamples);
     //  const matrixd K = computeCorrMatrix();
     matrixd K(nSamples,nSamples);
     computeCorrMatrix(K);
-    return utils::cholesky_decompose(K,mL);
+    size_t line_error = utils::cholesky_decompose(K,mL);
+    if (line_error) 
+      throw std::runtime_error("Cholesky decomposition error at line " + 
+			       boost::lexical_cast<std::string>(line_error));
   }
 
   matrixd NonParametricProcess::computeDerivativeCorrMatrix(int dth_index)
   {
     const size_t nSamples = mData.getNSamples();
     matrixd corrMatrix(nSamples,nSamples);
-    int error = mKernel.computeDerivativeCorrMatrix(mData.mX,corrMatrix,dth_index);
+    mKernel.computeDerivativeCorrMatrix(mData.mX,corrMatrix,dth_index);
     return corrMatrix;
   }
 

src/student_t_process_jef.cpp

     return d_;
   }
 
-  int StudentTProcessJeffreys::precomputePrediction()
+  void StudentTProcessJeffreys::precomputePrediction()
   {
     size_t n = mData.getNSamples();
     size_t p = mMean.nFeatures();
     mSigma = inner_prod(mAlphaF,mAlphaF)/(n-p);
     
     d_->setDof(n-p);  
-    return 0;
   }
 
 } //namespace bayesopt

src/student_t_process_nig.cpp

 
 
 
-  int StudentTProcessNIG::precomputePrediction()
+  void StudentTProcessNIG::precomputePrediction()
   {
     size_t n = mData.getNSamples();
     size_t p = mMean.getMeanFunc()->nFeatures();
     if ((boost::math::isnan(mWMap(0))) || (boost::math::isnan(mSigma)))
       {
 	FILE_LOG(logERROR) << "Error in precomputed prediction. NaN found.";
-	return -1;
+	throw std::runtime_error("Error in precomputed prediction. NaN found.");
       }
 
 
     if (dof <= 0)  
       {
+	dof = n;
 	FILE_LOG(logERROR) << "ERROR: Incorrect alpha. Dof invalid."
 			   << "Forcing Dof <= num of points.";
-	dof = n;
       }
 
     d_->setDof(dof);  
-    return 0;
   }
 
 } //namespace bayesopt

tests/FastDelegate.h

+//						FastDelegate.h 
+//	Efficient delegates in C++ that generate only two lines of asm code!
+//  Documentation is found at http://www.codeproject.com/cpp/FastDelegate.asp
+//
+//						- Don Clugston, Mar 2004.
+//		Major contributions were made by Jody Hagins.
+// History:
+// 24-Apr-04 1.0  * Submitted to CodeProject. 
+// 28-Apr-04 1.1  * Prevent most unsafe uses of evil static function hack.
+//				  * Improved syntax for horrible_cast (thanks Paul Bludov).
+//				  * Tested on Metrowerks MWCC and Intel ICL (IA32)
+//				  * Compiled, but not run, on Comeau C++ and Intel Itanium ICL.
+//	27-Jun-04 1.2 * Now works on Borland C++ Builder 5.5
+//				  * Now works on /clr "managed C++" code on VC7, VC7.1
+//				  * Comeau C++ now compiles without warnings.
+//				  * Prevent the virtual inheritance case from being used on 
+//					  VC6 and earlier, which generate incorrect code.
+//				  * Improved warning and error messages. Non-standard hacks
+//					 now have compile-time checks to make them safer.
+//				  * implicit_cast used instead of static_cast in many cases.
+//				  * If calling a const member function, a const class pointer can be used.
+//				  * MakeDelegate() global helper function added to simplify pass-by-value.
+//				  * Added fastdelegate.clear()
+// 16-Jul-04 1.2.1* Workaround for gcc bug (const member function pointers in templates)
+// 30-Oct-04 1.3  * Support for (non-void) return values.
+//				  * No more workarounds in client code!
+//					 MSVC and Intel now use a clever hack invented by John Dlugosz:
+//				     - The FASTDELEGATEDECLARE workaround is no longer necessary.
+//					 - No more warning messages for VC6
+//				  * Less use of macros. Error messages should be more comprehensible.
+//				  * Added include guards
+//				  * Added FastDelegate::empty() to test if invocation is safe (Thanks Neville Franks).
+//				  * Now tested on VS 2005 Express Beta, PGI C++
+// 24-Dec-04 1.4  * Added DelegateMemento, to allow collections of disparate delegates.
+//                * <,>,<=,>= comparison operators to allow storage in ordered containers.
+//				  * Substantial reduction of code size, especially the 'Closure' class.
+//				  * Standardised all the compiler-specific workarounds.
+//                * MFP conversion now works for CodePlay (but not yet supported in the full code).
+//                * Now compiles without warnings on _any_ supported compiler, including BCC 5.5.1
+//				  * New syntax: FastDelegate< int (char *, double) >. 
+// 14-Feb-05 1.4.1* Now treats =0 as equivalent to .clear(), ==0 as equivalent to .empty(). (Thanks elfric).
+//				  * Now tested on Intel ICL for AMD64, VS2005 Beta for AMD64 and Itanium.
+// 30-Mar-05 1.5  * Safebool idiom: "if (dg)" is now equivalent to "if (!dg.empty())"
+//				  * Fully supported by CodePlay VectorC
+//                * Bugfix for Metrowerks: empty() was buggy because a valid MFP can be 0 on MWCC!
+//                * More optimal assignment,== and != operators for static function pointers.
+
+#ifndef FASTDELEGATE_H
+#define FASTDELEGATE_H
+#if _MSC_VER > 1000
+#pragma once
+#endif // _MSC_VER > 1000
+
+#include <memory.h> // to allow <,> comparisons
+
+////////////////////////////////////////////////////////////////////////////////
+//						Configuration options
+//
+////////////////////////////////////////////////////////////////////////////////
+
+// Uncomment the following #define for optimally-sized delegates.
+// In this case, the generated asm code is almost identical to the code you'd get
+// if the compiler had native support for delegates.
+// It will not work on systems where sizeof(dataptr) < sizeof(codeptr). 
+// Thus, it will not work for DOS compilers using the medium model.
+// It will also probably fail on some DSP systems.
+#define FASTDELEGATE_USESTATICFUNCTIONHACK
+
+// Uncomment the next line to allow function declarator syntax.
+// It is automatically enabled for those compilers where it is known to work.
+//#define FASTDELEGATE_ALLOW_FUNCTION_TYPE_SYNTAX
+
+////////////////////////////////////////////////////////////////////////////////
+//						Compiler identification for workarounds
+//
+////////////////////////////////////////////////////////////////////////////////
+
+// Compiler identification. It's not easy to identify Visual C++ because
+// many vendors fraudulently define Microsoft's identifiers.
+#if defined(_MSC_VER) && !defined(__MWERKS__) && !defined(__VECTOR_C) && !defined(__ICL) && !defined(__BORLANDC__)
+#define FASTDLGT_ISMSVC
+
+#if (_MSC_VER <1300) // Many workarounds are required for VC6.
+#define FASTDLGT_VC6
+#pragma warning(disable:4786) // disable this ridiculous warning
+#endif
+
+#endif
+
+// Does the compiler uses Microsoft's member function pointer structure?
+// If so, it needs special treatment.
+// Metrowerks CodeWarrior, Intel, and CodePlay fraudulently define Microsoft's 
+// identifier, _MSC_VER. We need to filter Metrowerks out.
+#if defined(_MSC_VER) && !defined(__MWERKS__)
+#define FASTDLGT_MICROSOFT_MFP
+
+#if !defined(__VECTOR_C)
+// CodePlay doesn't have the __single/multi/virtual_inheritance keywords
+#define FASTDLGT_HASINHERITANCE_KEYWORDS
+#endif
+#endif
+
+// Does it allow function declarator syntax? The following compilers are known to work:
+#if defined(FASTDLGT_ISMSVC) && (_MSC_VER >=1310) // VC 7.1
+#define FASTDELEGATE_ALLOW_FUNCTION_TYPE_SYNTAX
+#endif
+
+// Gcc(2.95+), and versions of Digital Mars, Intel and Comeau in common use.
+#if defined (__DMC__) || defined(__GNUC__) || defined(__ICL) || defined(__COMO__)
+#define FASTDELEGATE_ALLOW_FUNCTION_TYPE_SYNTAX
+#endif
+
+// It works on Metrowerks MWCC 3.2.2. From boost.Config it should work on earlier ones too.
+#if defined (__MWERKS__)
+#define FASTDELEGATE_ALLOW_FUNCTION_TYPE_SYNTAX
+#endif
+
+#ifdef __GNUC__ // Workaround GCC bug #8271 
+	// At present, GCC doesn't recognize constness of MFPs in templates
+#define FASTDELEGATE_GCC_BUG_8271
+#endif
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+//						General tricks used in this code
+//
+// (a) Error messages are generated by typdefing an array of negative size to
+//     generate compile-time errors.
+// (b) Warning messages on MSVC are generated by declaring unused variables, and
+//	    enabling the "variable XXX is never used" warning.
+// (c) Unions are used in a few compiler-specific cases to perform illegal casts.
+// (d) For Microsoft and Intel, when adjusting the 'this' pointer, it's cast to
+//     (char *) first to ensure that the correct number of *bytes* are added.
+//
+////////////////////////////////////////////////////////////////////////////////
+//						Helper templates
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
+namespace fastdelegate {
+namespace detail {	// we'll hide the implementation details in a nested namespace.
+
+//		implicit_cast< >
+// I believe this was originally going to be in the C++ standard but 
+// was left out by accident. It's even milder than static_cast.
+// I use it instead of static_cast<> to emphasize that I'm not doing
+// anything nasty. 
+// Usage is identical to static_cast<>
+template <class OutputClass, class InputClass>
+inline OutputClass implicit_cast(InputClass input){
+	return input;
+}
+
+//		horrible_cast< >
+// This is truly evil. It completely subverts C++'s type system, allowing you 
+// to cast from any class to any other class. Technically, using a union 
+// to perform the cast is undefined behaviour (even in C). But we can see if
+// it is OK by checking that the union is the same size as each of its members.
+// horrible_cast<> should only be used for compiler-specific workarounds. 
+// Usage is identical to reinterpret_cast<>.
+
+// This union is declared outside the horrible_cast because BCC 5.5.1
+// can't inline a function with a nested class, and gives a warning.
+template <class OutputClass, class InputClass>
+union horrible_union{
+	OutputClass out;
+	InputClass in;
+};
+
+template <class OutputClass, class InputClass>
+inline OutputClass horrible_cast(const InputClass input){
+	horrible_union<OutputClass, InputClass> u;
+	// Cause a compile-time error if in, out and u are not the same size.
+	// If the compile fails here, it means the compiler has peculiar
+	// unions which would prevent the cast from working.
+	typedef int ERROR_CantUseHorrible_cast[sizeof(InputClass)==sizeof(u) 
+		&& sizeof(InputClass)==sizeof(OutputClass) ? 1 : -1];
+	u.in = input;
+	return u.out;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+//						Workarounds
+//
+////////////////////////////////////////////////////////////////////////////////
+
+// Backwards compatibility: This macro used to be necessary in the virtual inheritance
+// case for Intel and Microsoft. Now it just forward-declares the class.
+#define FASTDELEGATEDECLARE(CLASSNAME)	class CLASSNAME;
+
+// Prevent use of the static function hack with the DOS medium model.
+#ifdef __MEDIUM__
+#undef FASTDELEGATE_USESTATICFUNCTIONHACK
+#endif
+
+//			DefaultVoid - a workaround for 'void' templates in VC6.
+//
+//  (1) VC6 and earlier do not allow 'void' as a default template argument.
+//  (2) They also doesn't allow you to return 'void' from a function.
+//
+// Workaround for (1): Declare a dummy type 'DefaultVoid' which we use
+//   when we'd like to use 'void'. We convert it into 'void' and back
+//   using the templates DefaultVoidToVoid<> and VoidToDefaultVoid<>.
+// Workaround for (2): On VC6, the code for calling a void function is
+//   identical to the code for calling a non-void function in which the
+//   return value is never used, provided the return value is returned
+//   in the EAX register, rather than on the stack. 
+//   This is true for most fundamental types such as int, enum, void *.
+//   Const void * is the safest option since it doesn't participate 
+//   in any automatic conversions. But on a 16-bit compiler it might
+//   cause extra code to be generated, so we disable it for all compilers
+//   except for VC6 (and VC5).
+#ifdef FASTDLGT_VC6
+// VC6 workaround
+typedef const void * DefaultVoid;
+#else
+// On any other compiler, just use a normal void.
+typedef void DefaultVoid;
+#endif
+
+// Translate from 'DefaultVoid' to 'void'.
+// Everything else is unchanged
+template <class T>
+struct DefaultVoidToVoid { typedef T type; };
+
+template <>
+struct DefaultVoidToVoid<DefaultVoid> {	typedef void type; };
+
+// Translate from 'void' into 'DefaultVoid'
+// Everything else is unchanged
+template <class T>
+struct VoidToDefaultVoid { typedef T type; };
+
+template <>
+struct VoidToDefaultVoid<void> { typedef DefaultVoid type; };
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+//						Fast Delegates, part 1:
+//
+//		Conversion of member function pointer to a standard form
+//
+////////////////////////////////////////////////////////////////////////////////
+
+// GenericClass is a fake class, ONLY used to provide a type.
+// It is vitally important that it is never defined, so that the compiler doesn't
+// think it can optimize the invocation. For example, Borland generates simpler
+// code if it knows the class only uses single inheritance.
+
+// Compilers using Microsoft's structure need to be treated as a special case.
+#ifdef  FASTDLGT_MICROSOFT_MFP
+
+#ifdef FASTDLGT_HASINHERITANCE_KEYWORDS
+	// For Microsoft and Intel, we want to ensure that it's the most efficient type of MFP 
+	// (4 bytes), even when the /vmg option is used. Declaring an empty class 
+	// would give 16 byte pointers in this case....
+	class __single_inheritance GenericClass;
+#endif
+	// ...but for Codeplay, an empty class *always* gives 4 byte pointers.
+	// If compiled with the /clr option ("managed C++"), the JIT compiler thinks
+	// it needs to load GenericClass before it can call any of its functions,
+	// (compiles OK but crashes at runtime!), so we need to declare an 
+	// empty class to make it happy.
+	// Codeplay and VC4 can't cope with the unknown_inheritance case either.
+	class GenericClass {};
+#else
+	class GenericClass;
+#endif
+
+// The size of a single inheritance member function pointer.
+const int SINGLE_MEMFUNCPTR_SIZE = sizeof(void (GenericClass::*)());
+
+//						SimplifyMemFunc< >::Convert()
+//
+//	A template function that converts an arbitrary member function pointer into the 
+//	simplest possible form of member function pointer, using a supplied 'this' pointer.
+//  According to the standard, this can be done legally with reinterpret_cast<>.
+//	For (non-standard) compilers which use member function pointers which vary in size 
+//  depending on the class, we need to use	knowledge of the internal structure of a 
+//  member function pointer, as used by the compiler. Template specialization is used
+//  to distinguish between the sizes. Because some compilers don't support partial 
+//	template specialisation, I use full specialisation of a wrapper struct.
+
+// general case -- don't know how to convert it. Force a compile failure
+template <int N>
+struct SimplifyMemFunc {
+	template <class X, class XFuncType, class GenericMemFuncType>
+	inline static GenericClass *Convert(X *pthis, XFuncType function_to_bind, 
+		GenericMemFuncType &bound_func) { 
+		// Unsupported member function type -- force a compile failure.
+	    // (it's illegal to have a array with negative size).
+		typedef char ERROR_Unsupported_member_function_pointer_on_this_compiler[N-100];
+		return 0; 
+	}
+};
+
+// For compilers where all member func ptrs are the same size, everything goes here.
+// For non-standard compilers, only single_inheritance classes go here.
+template <>
+struct SimplifyMemFunc<SINGLE_MEMFUNCPTR_SIZE>  {	
+	template <class X, class XFuncType, class GenericMemFuncType>
+	inline static GenericClass *Convert(X *pthis, XFuncType function_to_bind, 
+			GenericMemFuncType &bound_func) {
+#if defined __DMC__  
+		// Digital Mars doesn't allow you to cast between abitrary PMF's, 
+		// even though the standard says you can. The 32-bit compiler lets you
+		// static_cast through an int, but the DOS compiler doesn't.
+		bound_func = horrible_cast<GenericMemFuncType>(function_to_bind);
+#else 
+        bound_func = reinterpret_cast<GenericMemFuncType>(function_to_bind);
+#endif
+        return reinterpret_cast<GenericClass *>(pthis);
+	}
+};
+
+////////////////////////////////////////////////////////////////////////////////
+//						Fast Delegates, part 1b:
+//
+//					Workarounds for Microsoft and Intel
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
+// Compilers with member function pointers which violate the standard (MSVC, Intel, Codeplay),
+// need to be treated as a special case.
+#ifdef FASTDLGT_MICROSOFT_MFP
+
+// We use unions to perform horrible_casts. I would like to use #pragma pack(push, 1)
+// at the start of each function for extra safety, but VC6 seems to ICE
+// intermittently if you do this inside a template.
+
+// __multiple_inheritance classes go here
+// Nasty hack for Microsoft and Intel (IA32 and Itanium)
+template<>
+struct SimplifyMemFunc< SINGLE_MEMFUNCPTR_SIZE + sizeof(int) >  {
+	template <class X, class XFuncType, class GenericMemFuncType>
+	inline static GenericClass *Convert(X *pthis, XFuncType function_to_bind, 
+		GenericMemFuncType &bound_func) { 
+		// We need to use a horrible_cast to do this conversion.
+		// In MSVC, a multiple inheritance member pointer is internally defined as:
+        union {
+			XFuncType func;
+			struct {	 
+				GenericMemFuncType funcaddress; // points to the actual member function
+				int delta;	     // #BYTES to be added to the 'this' pointer
+			}s;
+        } u;
+		// Check that the horrible_cast will work
+		typedef int ERROR_CantUsehorrible_cast[sizeof(function_to_bind)==sizeof(u.s)? 1 : -1];
+        u.func = function_to_bind;
+		bound_func = u.s.funcaddress;
+		return reinterpret_cast<GenericClass *>(reinterpret_cast<char *>(pthis) + u.s.delta); 
+	}
+};
+
+// virtual inheritance is a real nuisance. It's inefficient and complicated.
+// On MSVC and Intel, there isn't enough information in the pointer itself to
+// enable conversion to a closure pointer. Earlier versions of this code didn't
+// work for all cases, and generated a compile-time error instead.
+// But a very clever hack invented by John M. Dlugosz solves this problem.
+// My code is somewhat different to his: I have no asm code, and I make no 
+// assumptions about the calling convention that is used.
+
+// In VC++ and ICL, a virtual_inheritance member pointer 
+// is internally defined as:
+struct MicrosoftVirtualMFP {
+	void (GenericClass::*codeptr)(); // points to the actual member function
+	int delta;		// #bytes to be added to the 'this' pointer
+	int vtable_index; // or 0 if no virtual inheritance
+};
+// The CRUCIAL feature of Microsoft/Intel MFPs which we exploit is that the
+// m_codeptr member is *always* called, regardless of the values of the other
+// members. (This is *not* true for other compilers, eg GCC, which obtain the
+// function address from the vtable if a virtual function is being called).
+// Dlugosz's trick is to make the codeptr point to a probe function which
+// returns the 'this' pointer that was used.
+
+// Define a generic class that uses virtual inheritance.
+// It has a trival member function that returns the value of the 'this' pointer.
+struct GenericVirtualClass : virtual public GenericClass
+{
+	typedef GenericVirtualClass * (GenericVirtualClass::*ProbePtrType)();
+	GenericVirtualClass * GetThis() { return this; }
+};
+
+// __virtual_inheritance classes go here
+template <>
+struct SimplifyMemFunc<SINGLE_MEMFUNCPTR_SIZE + 2*sizeof(int) >
+{
+
+	template <class X, class XFuncType, class GenericMemFuncType>
+	inline static GenericClass *Convert(X *pthis, XFuncType function_to_bind, 
+		GenericMemFuncType &bound_func) {
+		union {
+			XFuncType func;
+			GenericClass* (X::*ProbeFunc)();
+			MicrosoftVirtualMFP s;
+		} u;
+		u.func = function_to_bind;
+		bound_func = reinterpret_cast<GenericMemFuncType>(u.s.codeptr);
+		union {
+			GenericVirtualClass::ProbePtrType virtfunc;
+			MicrosoftVirtualMFP s;
+		} u2;
+		// Check that the horrible_cast<>s will work
+		typedef int ERROR_CantUsehorrible_cast[sizeof(function_to_bind)==sizeof(u.s)
+			&& sizeof(function_to_bind)==sizeof(u.ProbeFunc)
+			&& sizeof(u2.virtfunc)==sizeof(u2.s) ? 1 : -1];
+   // Unfortunately, taking the address of a MF prevents it from being inlined, so 
+   // this next line can't be completely optimised away by the compiler.
+		u2.virtfunc = &GenericVirtualClass::GetThis;
+		u.s.codeptr = u2.s.codeptr;
+		return (pthis->*u.ProbeFunc)();
+	}
+};
+
+#if (_MSC_VER <1300)
+
+// Nasty hack for Microsoft Visual C++ 6.0
+// unknown_inheritance classes go here
+// There is a compiler bug in MSVC6 which generates incorrect code in this case!!
+template <>
+struct SimplifyMemFunc<SINGLE_MEMFUNCPTR_SIZE + 3*sizeof(int) >
+{
+	template <class X, class XFuncType, class GenericMemFuncType>
+	inline static GenericClass *Convert(X *pthis, XFuncType function_to_bind, 
+		GenericMemFuncType &bound_func) {
+		// There is an apalling but obscure compiler bug in MSVC6 and earlier:
+		// vtable_index and 'vtordisp' are always set to 0 in the 
+		// unknown_inheritance case!
+		// This means that an incorrect function could be called!!!
+		// Compiling with the /vmg option leads to potentially incorrect code.
+		// This is probably the reason that the IDE has a user interface for specifying
+		// the /vmg option, but it is disabled -  you can only specify /vmg on 
+		// the command line. In VC1.5 and earlier, the compiler would ICE if it ever
+		// encountered this situation.
+		// It is OK to use the /vmg option if /vmm or /vms is specified.
+
+		// Fortunately, the wrong function is only called in very obscure cases.
+		// It only occurs when a derived class overrides a virtual function declared 
+		// in a virtual base class, and the member function 
+		// points to the *Derived* version of that function. The problem can be
+		// completely averted in 100% of cases by using the *Base class* for the 
+		// member fpointer. Ie, if you use the base class as an interface, you'll
+		// stay out of trouble.
+		// Occasionally, you might want to point directly to a derived class function
+		// that isn't an override of a base class. In this case, both vtable_index 
+		// and 'vtordisp' are zero, but a virtual_inheritance pointer will be generated.
+		// We can generate correct code in this case. To prevent an incorrect call from
+		// ever being made, on MSVC6 we generate a warning, and call a function to 
+		// make the program crash instantly. 
+		typedef char ERROR_VC6CompilerBug[-100];
+		return 0; 
+	}
+};
+
+
+#else 
+
+// Nasty hack for Microsoft and Intel (IA32 and Itanium)
+// unknown_inheritance classes go here 
+// This is probably the ugliest bit of code I've ever written. Look at the casts!
+// There is a compiler bug in MSVC6 which prevents it from using this code.
+template <>
+struct SimplifyMemFunc<SINGLE_MEMFUNCPTR_SIZE + 3*sizeof(int) >
+{
+	template <class X, class XFuncType, class GenericMemFuncType>
+	inline static GenericClass *Convert(X *pthis, XFuncType function_to_bind, 
+			GenericMemFuncType &bound_func) {
+		// The member function pointer is 16 bytes long. We can't use a normal cast, but
+		// we can use a union to do the conversion.
+		union {
+			XFuncType func;
+			// In VC++ and ICL, an unknown_inheritance member pointer 
+			// is internally defined as:
+			struct {
+				GenericMemFuncType m_funcaddress; // points to the actual member function
+				int delta;		// #bytes to be added to the 'this' pointer
+				int vtordisp;		// #bytes to add to 'this' to find the vtable
+				int vtable_index; // or 0 if no virtual inheritance
+			} s;
+		} u;
+		// Check that the horrible_cast will work
+		typedef int ERROR_CantUsehorrible_cast[sizeof(XFuncType)==sizeof(u.s)? 1 : -1];
+		u.func = function_to_bind;
+		bound_func = u.s.funcaddress;
+		int virtual_delta = 0;
+		if (u.s.vtable_index) { // Virtual inheritance is used
+			// First, get to the vtable. 
+			// It is 'vtordisp' bytes from the start of the class.
+			const int * vtable = *reinterpret_cast<const int *const*>(
+				reinterpret_cast<const char *>(pthis) + u.s.vtordisp );
+
+			// 'vtable_index' tells us where in the table we should be looking.
+			virtual_delta = u.s.vtordisp + *reinterpret_cast<const int *>( 
+				reinterpret_cast<const char *>(vtable) + u.s.vtable_index);
+		}
+		// The int at 'virtual_delta' gives us the amount to add to 'this'.
+        // Finally we can add the three components together. Phew!
+        return reinterpret_cast<GenericClass *>(
+			reinterpret_cast<char *>(pthis) + u.s.delta + virtual_delta);
+	};
+};
+#endif // MSVC 7 and greater
+
+#endif // MS/Intel hacks
+
+}  // namespace detail
+
+////////////////////////////////////////////////////////////////////////////////
+//						Fast Delegates, part 2:
+//
+//	Define the delegate storage, and cope with static functions
+//
+////////////////////////////////////////////////////////////////////////////////
+
+// DelegateMemento -- an opaque structure which can hold an arbitary delegate.
+// It knows nothing about the calling convention or number of arguments used by
+// the function pointed to.
+// It supplies comparison operators so that it can be stored in STL collections.
+// It cannot be set to anything other than null, nor invoked directly: 
+//   it must be converted to a specific delegate.
+
+// Implementation:
+// There are two possible implementations: the Safe method and the Evil method.
+//				DelegateMemento - Safe version
+//
+// This implementation is standard-compliant, but a bit tricky.
+// A static function pointer is stored inside the class. 
+// Here are the valid values:
+// +-- Static pointer --+--pThis --+-- pMemFunc-+-- Meaning------+
+// |   0				|  0       |   0        | Empty          |
+// |   !=0              |(dontcare)|  Invoker   | Static function|
+// |   0                |  !=0     |  !=0*      | Method call    |
+// +--------------------+----------+------------+----------------+
+//  * For Metrowerks, this can be 0. (first virtual function in a 
+//       single_inheritance class).
+// When stored stored inside a specific delegate, the 'dontcare' entries are replaced
+// with a reference to the delegate itself. This complicates the = and == operators
+// for the delegate class.
+
+//				DelegateMemento - Evil version
+//
+// For compilers where data pointers are at least as big as code pointers, it is 
+// possible to store the function pointer in the this pointer, using another 
+// horrible_cast. In this case the DelegateMemento implementation is simple:
+// +--pThis --+-- pMemFunc-+-- Meaning---------------------+
+// |    0     |  0         | Empty                         |
+// |  !=0     |  !=0*      | Static function or method call|
+// +----------+------------+-------------------------------+
+//  * For Metrowerks, this can be 0. (first virtual function in a 
+//       single_inheritance class).
+// Note that the Sun C++ and MSVC documentation explicitly state that they 
+// support static_cast between void * and function pointers.
+
+class DelegateMemento {
+protected: 
+	// the data is protected, not private, because many
+	// compilers have problems with template friends.
+	typedef void (detail::GenericClass::*GenericMemFuncType)(); // arbitrary MFP.
+	detail::GenericClass *m_pthis;
+	GenericMemFuncType m_pFunction;
+
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+	typedef void (*GenericFuncPtr)(); // arbitrary code pointer
+	GenericFuncPtr m_pStaticFunction;
+#endif
+
+public:
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+	DelegateMemento() : m_pthis(0), m_pFunction(0), m_pStaticFunction(0) {};
+	void clear() {
+		m_pthis=0; m_pFunction=0; m_pStaticFunction=0;
+	}
+#else
+	DelegateMemento() : m_pthis(0), m_pFunction(0) {};
+	void clear() {	m_pthis=0; m_pFunction=0;	}
+#endif
+public:
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+	inline bool IsEqual (const DelegateMemento &x) const{
+	    // We have to cope with the static function pointers as a special case
+		if (m_pFunction!=x.m_pFunction) return false;
+		// the static function ptrs must either both be equal, or both be 0.
+		if (m_pStaticFunction!=x.m_pStaticFunction) return false;
+		if (m_pStaticFunction!=0) return m_pthis==x.m_pthis;
+		else return true;
+	}
+#else // Evil Method
+	inline bool IsEqual (const DelegateMemento &x) const{
+		return m_pthis==x.m_pthis && m_pFunction==x.m_pFunction;
+	}
+#endif
+	// Provide a strict weak ordering for DelegateMementos.
+	inline bool IsLess(const DelegateMemento &right) const {
+		// deal with static function pointers first
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+		if (m_pStaticFunction !=0 || right.m_pStaticFunction!=0) 
+				return m_pStaticFunction < right.m_pStaticFunction;
+#endif
+		if (m_pthis !=right.m_pthis) return m_pthis < right.m_pthis;
+	// There are no ordering operators for member function pointers, 
+	// but we can fake one by comparing each byte. The resulting ordering is
+	// arbitrary (and compiler-dependent), but it permits storage in ordered STL containers.
+		return memcmp(&m_pFunction, &right.m_pFunction, sizeof(m_pFunction)) < 0;
+
+	}
+	// BUGFIX (Mar 2005):
+	// We can't just compare m_pFunction because on Metrowerks,
+	// m_pFunction can be zero even if the delegate is not empty!
+	inline bool operator ! () const		// Is it bound to anything?
+	{ return m_pthis==0 && m_pFunction==0; }
+	inline bool empty() const		// Is it bound to anything?
+	{ return m_pthis==0 && m_pFunction==0; }
+public:
+	DelegateMemento & operator = (const DelegateMemento &right)  {
+		SetMementoFrom(right); 
+		return *this;
+	}
+	inline bool operator <(const DelegateMemento &right) {
+		return IsLess(right);
+	}
+	inline bool operator >(const DelegateMemento &right) {
+		return right.IsLess(*this);
+	}
+	DelegateMemento (const DelegateMemento &right)  : 
+		m_pFunction(right.m_pFunction), m_pthis(right.m_pthis)
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+		, m_pStaticFunction (right.m_pStaticFunction)
+#endif
+		{}
+protected:
+	void SetMementoFrom(const DelegateMemento &right)  {
+		m_pFunction = right.m_pFunction;
+		m_pthis = right.m_pthis;
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+		m_pStaticFunction = right.m_pStaticFunction;
+#endif
+	}
+};
+
+
+//						ClosurePtr<>
+//
+// A private wrapper class that adds function signatures to DelegateMemento.
+// It's the class that does most of the actual work.
+// The signatures are specified by:
+// GenericMemFunc: must be a type of GenericClass member function pointer. 
+// StaticFuncPtr:  must be a type of function pointer with the same signature 
+//                 as GenericMemFunc.
+// UnvoidStaticFuncPtr: is the same as StaticFuncPtr, except on VC6
+//                 where it never returns void (returns DefaultVoid instead).
+
+// An outer class, FastDelegateN<>, handles the invoking and creates the
+// necessary typedefs.
+// This class does everything else.
+
+namespace detail {
+
+template < class GenericMemFunc, class StaticFuncPtr, class UnvoidStaticFuncPtr>
+class ClosurePtr : public DelegateMemento {
+public:
+	// These functions are for setting the delegate to a member function.
+
+	// Here's the clever bit: we convert an arbitrary member function into a 
+	// standard form. XMemFunc should be a member function of class X, but I can't 
+	// enforce that here. It needs to be enforced by the wrapper class.
+	template < class X, class XMemFunc >
+	inline void bindmemfunc(X *pthis, XMemFunc function_to_bind ) {
+		m_pthis = SimplifyMemFunc< sizeof(function_to_bind) >
+			::Convert(pthis, function_to_bind, m_pFunction);
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+		m_pStaticFunction = 0;
+#endif
+	}
+	// For const member functions, we only need a const class pointer.
+	// Since we know that the member function is const, it's safe to 
+	// remove the const qualifier from the 'this' pointer with a const_cast.
+	// VC6 has problems if we just overload 'bindmemfunc', so we give it a different name.
+	template < class X, class XMemFunc>
+	inline void bindconstmemfunc(const X *pthis, XMemFunc function_to_bind) {
+		m_pthis= SimplifyMemFunc< sizeof(function_to_bind) >
+			::Convert(const_cast<X*>(pthis), function_to_bind, m_pFunction);
+#if !defined(FASTDELEGATE_USESTATICFUNCTIONHACK)
+		m_pStaticFunction = 0;
+#endif
+	}