Commits

Ruben Martinez-Cantin committed ca1e395

Separate compilation of examples. Solved bug while finding python libs.

Comments (0)

Files changed (19)

 SET(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/lib)
 SET(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/bin)
 
-option(BUILD_PYTHON_INTERFACE "Build Python Interface?" OFF)
+option(BUILD_EXAMPLES "Build examples and demos?" ON)
+option(BUILD_PYTHON_INTERFACE "Build Python interface?" OFF)
 option(MATLAB_COMPATIBLE "Build library compatible with Matlab?" ON)
 
 if(BUILD_PYTHON_INTERFACE)
 TARGET_LINK_LIBRARIES(bayesopt
   ${EXT_LIBS} ${PYTHON_LIB} )
 
-#Continuous test
-ADD_EXECUTABLE(bo_cont ./app/bo_cont.cpp)
-add_dependencies(bo_cont bayesopt)
-TARGET_LINK_LIBRARIES(bo_cont bayesopt)
 
-#Discrete test
-ADD_EXECUTABLE(bo_disc ./app/bo_disc.cpp)
-add_dependencies(bo_disc bayesopt)
-TARGET_LINK_LIBRARIES(bo_disc bayesopt)
-
-#1D test
-ADD_EXECUTABLE(bo_oned ./app/bo_oned.cpp)
-add_dependencies(bo_oned bayesopt)
-TARGET_LINK_LIBRARIES(bo_oned bayesopt)
-
-#Display test
-find_package(GLUT)
-find_package(OpenGL)
-if(OPENGL_FOUND AND GLUT_FOUND)
-  INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/matplotpp)
-
-  include_directories(${GLUT_INCLUDE_DIRS})
-  link_directories(${GLUT_LIBRARY_DIRS})
-  add_definitions(${GLUT_DEFINITIONS})
-
-  include_directories(${OpenGL_INCLUDE_DIRS})
-  link_directories(${OpenGL_LIBRARY_DIRS})
-  add_definitions(${OpenGL_DEFINITIONS})
-
-  enable_language(C)
-  ADD_EXECUTABLE(bo_display 
-    ./app/bo_display.cpp 
-    ./matplotpp/matplotpp.cc 
-    ./matplotpp/gl2ps.c)
-  add_dependencies(bo_display bayesopt)
-  TARGET_LINK_LIBRARIES(bo_display bayesopt ${GLUT_LIBRARY} ${OPENGL_LIBRARY})
-endif()
-
-#Branin
-ADD_EXECUTABLE(bo_branin ./app/bo_branin.cpp )
-add_dependencies(bo_branin bayesopt)
-TARGET_LINK_LIBRARIES(bo_branin bayesopt)
-
-#Test for random number generator
-#ADD_EXECUTABLE(randtest ./app/testrand.cpp)
-
-#Test for parsers
-#ADD_EXECUTABLE(parsetest ./utils/parser.cpp ./app/testparser.cpp)
-
+IF(BUILD_EXAMPLES)
+  ADD_SUBDIRECTORY(examples)
+endif(BUILD_EXAMPLES)
 
 INSTALL(FILES 
   ./include/bayesoptcont.hpp 

app/bo_branin.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/>.
-------------------------------------------------------------------------
-*/
-
-#define _USE_MATH_DEFINES
-#include <cmath>
-#include <valarray>
-#include "bayesoptcont.hpp"
-
-#ifndef M_PI
-#define M_PI           3.14159265358979323846
-#endif
-
-class ExampleBranin: public bayesopt::ContinuousModel
-{
-public:
-  ExampleBranin(size_t dim,bopt_params par):
-    ContinuousModel(dim,par) {}
-
-  double evaluateSample( const vectord& xin)
-  {
-     if (xin.size() != 2)
-      {
-	std::cout << "WARNING: This only works for 2D inputs." << std::endl
-		  << "WARNING: Using only first two components." << std::endl;
-      }
-
-    double x = xin(0) * 15 - 5;
-    double y = xin(1) * 15;
-
-    return sqr(y-(5.1/(4*sqr(M_PI)))*sqr(x)+5*x/M_PI-6)+10*(1-1/(8*M_PI))*cos(x)+10;
-  };
-
-  bool checkReachability(const vectord &query)
-  {return true;};
-
-  inline double sqr( double x ){ return x*x; };
-
-  void printOptimal()
-  {
-    vectord sv(2);  
-    sv(0) = 0.1239; sv(1) = 0.8183;
-    std::cout << "Solutions: " << sv << "->" 
-	      << evaluateSample(sv) << std::endl;
-    sv(0) = 0.5428; sv(1) = 0.1517;
-    std::cout << "Solutions: " << sv << "->" 
-	      << evaluateSample(sv) << std::endl;
-    sv(0) = 0.9617; sv(1) = 0.1650;
-    std::cout << "Solutions: " << sv << "->" 
-	      << evaluateSample(sv) << std::endl;
-  }
-
-};
-
-int main(int nargs, char *args[])
-{
-  bopt_params par = initialize_parameters_to_default();
-  par.n_iterations = 200;
-  par.n_init_samples = 50;
-  par.kernel.hp_mean[0] = 1.0;
-  par.kernel.n_hp = 1;
-  par.crit_name = "cHedge(cEI,cLCB,cPOI)";
-  
-  ExampleBranin branin(2,par);
-  vectord result(2);
-
-  branin.optimize(result);
-  std::cout << "Result:" << result << std::endl;
-  branin.printOptimal();
-
-  return 0;
-}

app/bo_cont.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 <ctime>
-#include "bayesoptwpr.h"                 // For the C API
-#include "bayesoptcont.hpp"              // For the C++ API
-
-
-/* Function to be used for C-API testing */
-double testFunction(unsigned int n, const double *x,
-		    double *gradient, /* NULL if not needed */
-		    void *func_data)
-{
-  double f = 10.;
-  for (unsigned int i = 0; i < n; ++i)
-    {
-      f += (x[i] - .53) * (x[i] - .53);
-    }
-  return f;
-}
-
-/* Class to be used for C++-API testing */
-class ExampleQuadratic: public bayesopt::ContinuousModel
-{
- public:
-
-  ExampleQuadratic(size_t dim,bopt_params param):
-    ContinuousModel(dim,param) {}
-
-  double evaluateSample( const vectord &Xi ) 
-  {
-    double x[100];
-    for (size_t i = 0; i < Xi.size(); ++i)
-	x[i] = Xi(i);
-    return testFunction(Xi.size(),x,NULL,NULL);
-  };
-
-
-  bool checkReachability( const vectord &query )
-  { return true; };
- 
-};
-
-
-int main(int nargs, char *args[])
-{    
-  int n = 10;                   // Number of dimensions
-
-  // Common configuration
-  // See parameters.h for the available options
-  // If we initialize the struct with the DEFAUL_PARAMS,
-  // the we can optionally change only few of them 
-  bopt_params par = initialize_parameters_to_default();
-
-  par.kernel.hp_mean[0] = KERNEL_THETA;
-  par.kernel.hp_mean[1] = KERNEL_THETA;
-  par.kernel.hp_std[0] = 1;
-  par.kernel.hp_std[1] = 1;
-  par.kernel.n_hp = 2;
-  par.mean.coef_mean[0] = 1.0;
-  par.mean.coef_std[0] = MEAN_SIGMA;
-  par.mean.n_coef = 1;
-  par.alpha = PRIOR_ALPHA;
-  par.beta = PRIOR_BETA;
-  par.noise = DEFAULT_NOISE;
-  par.surr_name = "sStudentTProcessJef";
-  par.kernel.name = "kSum(kSEISO,kConst)";
-  par.mean.name = "mConst";
-  par.l_type = L_ML;
-  par.n_iterations = 200;       // Number of iterations
-  par.n_init_samples = 50;
-  par.verbose_level = 2;
-  /*******************************************/
-
-  clock_t start, end;
-  double diff,diff2;
-
-  std::cout << "Running C++ interface" << std::endl;
-  // Configure C++ interface
-
-  ExampleQuadratic opt(n,par);
-  vectord result(n);
-
-  // Run C++ interface
-  start = clock();
-  opt.optimize(result);
-  end = clock();
-  diff = (double)(end-start) / (double)CLOCKS_PER_SEC;
-  /*******************************************/
-
-
-  std::cout << "Running C inferface" << std::endl;
-  
-  // Configure C interface
-  double u[128], x[128], l[128], fmin[128];
-
-  for (int i = 0; i < n; ++i) {
-    l[i] = 0.;    u[i] = 1.;
-  }
-
-  // Run C interface
-  start = clock();
-  bayes_optimization(n,&testFunction,NULL,l,u,x,fmin,par);
-  end = clock();
-  diff2 = (double)(end-start) / (double)CLOCKS_PER_SEC;
-  /*******************************************/
-
-
-  // Results
-  std::cout << "Final result C++: " << result << std::endl;
-  std::cout << "Elapsed time in C++: " << diff << " seconds" << std::endl;
-
-  std::cout << "Final result C: [" << n <<"](" << x[0];
-  for (int i = 1; i < n; ++i )
-    std::cout << "," << x[i];
-  
-  std::cout << ")" << std::endl;
-  std::cout << "Elapsed time in C: " << diff2 << " seconds" << std::endl;
-
-}
-

app/bo_disc.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 "bayesoptwpr.h"               // For the C API
-#include "bayesoptdisc.hpp"            // For the C++ API
-#include "lhs.hpp"
-
-
-/* Function to be used for C-API testing */
-double testFunction(unsigned int n, const double *x,
-		    double *gradient, /* NULL if not needed */
-		    void *func_data)
-{
-  double f = 10.;
-  for (unsigned int i = 0; i < n; ++i)
-    {
-      f += (x[i] - .53) * (x[i] - .53);
-    }
-  return f;
-}
-
-/* Class to be used for C++-API testing */
-class ExampleDisc: public bayesopt::DiscreteModel
-{
- public:
-
-  ExampleDisc(const vecOfvec & validSet, bopt_params param):
-    DiscreteModel(validSet,param) {}
-
-  double evaluateSample( const vectord &Xi ) 
-  {
-    double x[100];
-    for (size_t i = 0; i < Xi.size(); ++i)
-	x[i] = Xi(i);
-    return testFunction(Xi.size(),x,NULL,NULL);
-  };
-
-  bool checkReachability( const vectord &query )
-  { return true; }; 
-};
-
-
-int main(int nargs, char *args[])
-{    
-  int n = 6;                   // Number of dimensions
-
-  // Common configuration
-  // See parameters.h for the available options
-  // If we initialize the struct with the DEFAUL_PARAMS,
-  // the we can optionally change only few of them 
-  bopt_params par = initialize_parameters_to_default();
-
-  par.kernel.hp_mean[0] = KERNEL_THETA;
-  par.kernel.hp_std[0] = 1.0;
-  par.kernel.n_hp = 1;
-  par.mean.coef_mean[0] = 0.0;
-  par.mean.coef_std[0] = MEAN_SIGMA;
-  par.mean.n_coef = 1;
-  par.noise = DEFAULT_NOISE;
-  par.surr_name = "sStudentTProcessJef";
-  par.n_iterations = 20;       // Number of iterations
-  par.n_init_samples = 20;
-  /*******************************************/
-  
-  size_t nPoints = 1000;
-
-  randEngine mtRandom;
-  matrixd xPoints(nPoints,n);
-  vecOfvec xP;
-
-  //Thanks to the quirks of Visual Studio, the following expression is invalid,
-  //so we have to replace it by a literal.
-  //WARNING: Make sure to update the size of the array if the number of points
-  //or dimensions change.
-#ifdef _MSC_VER
-  double xPointsArray[6000];
-#else
-  const size_t nPinArr = n*nPoints;
-  double xPointsArray[nPinArr];
-#endif
-  
-  bayesopt::utils::lhs(xPoints,mtRandom);
-
-  for(size_t i = 0; i<nPoints; ++i)
-    {
-      vectord point = row(xPoints,i);  
-      xP.push_back(point);
-      for(size_t j=0; j<n; ++j)
-	{
-	  xPointsArray[i*n+j] = point(j);	  
-	}
-    }
-    
-
-  
-  // Run C++ interface
-  std::cout << "Running C++ interface" << std::endl;
-  ExampleDisc opt(xP,par);
-  vectord result(n);
-  opt.optimize(result);
-
-  
-  // Run C interface
-  std::cout << "Running C interface" << std::endl;
-  double x[128], fmin[128];
-  bayes_optimization_disc(n, &testFunction, NULL, xPointsArray, nPoints,
-			  x, fmin, par);
-
-  // Find the optimal value
-  size_t min = 0;
-  double minvalue = opt.evaluateSample(row(xPoints,min));
-  for(size_t i = 1; i<nPoints; ++i)
-    {
-      vectord point = row(xPoints,i);  
-      if (opt.evaluateSample(point) < minvalue)
-	{
-	  min = i;
-	  minvalue = opt.evaluateSample(row(xPoints,min));
-	  std::cout << i << "," << minvalue << std::endl;
-	}
-    }
-
-  std::cout << "Final result C++: " << result << std::endl;
-  std::cout << "Final result C: (";
-  for (int i = 0; i < n; i++ )
-    std::cout << x[i] << ", ";
-  std::cout << ")" << std::endl;
-  std::cout << "Optimal: " << row(xPoints,min) << std::endl;
-}

app/bo_display.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 <valarray>
-#include "bayesoptcont.hpp"
-
-class ExampleOneD: public bayesopt::ContinuousModel
-{
-public:
-  ExampleOneD(size_t dim, bopt_params par):
-    ContinuousModel(dim,par) {}
-
-  double evaluateSample(const vectord& xin)
-  {
-    if (xin.size() > 1)
-      {
-	std::cout << "WARNING: This only works for 1D inputs." << std::endl
-		  << "WARNING: Using only first component." << std::endl;
-      }
-
-    double x = xin(0);
-    return sqr(x-0.3) + sin(20*x)*0.2;
-  };
-
-  bool checkReachability(const vectord &query)
-  {return true;};
-
-  inline double sqr( double x ){ return x*x; }
-
-  void printOptimal()
-  {
-    std::cout << "Optimal:" << 0.5 << std::endl;
-  }
-
-};
-
-//#include "unistd.h"
-//using namespace std;
-#include "matplotpp.h"  
-
-using namespace bayesopt;
-
-int is_run=0;
-int is_step=0;
-size_t state_ii = 0;
-BayesOptBase* GLOBAL_MODEL;
-std::vector<double> lx,ly;
-
-class MP :public MatPlot{ 
-void DISPLAY(){
-  size_t nruns = GLOBAL_MODEL->getParameters()->n_iterations;
-  if ((is_run) && (state_ii < nruns))
-    {
-      ++state_ii;
-      GLOBAL_MODEL->stepOptimization(state_ii); 
-      vectord last(1);
-      double res = GLOBAL_MODEL->getSurrogateModel()->getLastSample(last);
-      ly.push_back(res);
-      lx.push_back(last(0));
-	  
-      if (is_step) { is_run = 0; is_step = 0; }
-    }
-    
-  int n=1000;
-  std::vector<double> x,y,z,su,sl,c;
-  x=linspace(0,1,n);
-  y = x; z = x; su = x; sl = x; c= x;
-  vectord q(1);
-  for(int i=0;i<n;++i)
-    {
-      q(0) = x[i];
-      ProbabilityDistribution* pd = GLOBAL_MODEL->getSurrogateModel()->prediction(q);
-      y[i] = pd->getMean();
-      su[i] = y[i] + 2*pd->getStd();
-      sl[i] = y[i] - 2*pd->getStd();
-      c[i] = -GLOBAL_MODEL->evaluateCriteria(q);
-      z[i] = GLOBAL_MODEL->evaluateSample(q);
-    }
- 
-  //plot
-  subplot(2,1,1);
-  title("press r to run and stop");
-  plot(x,y); set(3);
-  plot(lx,ly);set("k");set("*");
-  plot(x,su);set("g"); set(2);
-  plot(x,sl);set("g"); set(2);
-  plot(x,z);set("r"); set(3);
-  
-  subplot(2,1,2);
-  plot(x,c); set(3);
-}
-}mp;
-
-void display(){mp.display(); }
-void reshape(int w,int h){ mp.reshape(w,h); }
-void idle( void )
-{
-  glutPostRedisplay();
-}
-
-void mouse(int button, int state, int x, int y ){ mp.mouse(button,state,x,y); }
-void motion(int x, int y ){mp.motion(x,y); }
-void passive(int x, int y ){mp.passivemotion(x,y); }
-void keyboard(unsigned char key, int x, int y){
-    mp.keyboard(key, x, y); 
-    if(key=='r'){ if(is_run==0){is_run=1;}else{is_run=0;}}
-    if(key=='s'){ is_run=1; is_step=1; } 
-}
-
-int main(int nargs, char *args[])
-{
-  size_t dim = 1;
-  bopt_params parameters = initialize_parameters_to_default();
-  parameters.n_init_samples = 7;
-  parameters.n_iter_relearn = 0;
-  parameters.n_iterations = 300;
-  parameters.verbose_level = 2;
-
-  // Surrogate models
-  //  parameters.surr_name = "sStudentTProcessNIG";
-  parameters.surr_name = "sGaussianProcessNormal";
-
-  // Criterion model
-  // parameters.crit_name = "cAopt";
-  // parameters.n_crit_params = 0;
-
-  parameters.crit_name = "cEI";
-  parameters.crit_params[0] = 1;
-  parameters.n_crit_params = 1;
-
-  // parameters.crit_name = "cLCB";
-  // parameters.crit_params[0] = 5;
-  // parameters.n_crit_params = 1;
-
-  // Kernel models
-  // parameters.kernel.name = "kSum(kPoly3,kRQISO)";
-  // double mean[128] = {1, 1, 1, 1};
-  // double std[128] = {10, 10, 10, 10};
-  // size_t nhp = 4;
-  // memcpy(parameters.kernel.hp_mean, mean, nhp * sizeof(double));
-  // memcpy(parameters.kernel.hp_std,std, nhp * sizeof(double));
-  // parameters.kernel.n_hp = nhp;
-
-  parameters.kernel.name = "kMaternISO3";
-  parameters.kernel.hp_mean[0] = 1;
-  parameters.kernel.hp_std[0] = 5;
-  parameters.kernel.n_hp = 1;
-
-  state_ii = 0;
-
-  ExampleOneD* opt = new ExampleOneD(dim,parameters);
-  vectord result(dim);
-  GLOBAL_MODEL = opt;
-  opt->initializeOptimization();
-  vectord last(1);
-  size_t n_points = GLOBAL_MODEL->getSurrogateModel()->getNSamples();
-  for (size_t i = 0; i<n_points;++i)
-    {
-      double res = GLOBAL_MODEL->getSurrogateModel()->getSample(i,last);
-      ly.push_back(res);
-      lx.push_back(last(0));
-    }
-
-  glutInit(&nargs, args);
-  glutCreateWindow(50,50,800,650);
-  glutDisplayFunc( display );
-  glutReshapeFunc( reshape );
-  glutIdleFunc( idle );
-  glutMotionFunc( motion );
-  glutMouseFunc( mouse );
-  glutPassiveMotionFunc(passive);    
-  glutKeyboardFunc( keyboard );        
-  glutMainLoop();    
-
-  delete opt;
-  return 0;
-}

app/bo_oned.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 <valarray>
-#include "bayesoptcont.hpp"
-
-class ExampleOneD: public bayesopt::ContinuousModel
-{
-public:
-  ExampleOneD(size_t dim, bopt_params par):
-    ContinuousModel(dim,par) {}
-
-  double evaluateSample(const vectord& xin)
-  {
-    if (xin.size() > 1)
-      {
-	std::cout << "WARNING: This only works for 1D inputs." << std::endl
-		  << "WARNING: Using only first component." << std::endl;
-      }
-
-    double x = xin(0);
-    return sqr(x-0.5) + sin(20*x)*0.2;
-  };
-
-  bool checkReachability(const vectord &query)
-  {return true;};
-
-  inline double sqr( double x ){ return x*x; }
-
-  void printOptimal()
-  {
-    std::cout << "Optimal:" << 0.5 << std::endl;
-  }
-
-};
-
-int main(int nargs, char *args[])
-{
-  size_t dim = 1;
-  bopt_params parameters = initialize_parameters_to_default();
-  parameters.n_init_samples = 10;
-  parameters.n_iterations = 300;
-  parameters.surr_name = "sGaussianProcessML";
-  /*  parameters.kernel.hp_mean[0] = 1.0;
-  parameters.kernel.hp_std[0] = 100.0;
-  parameters.kernel.n_hp = 1;
-  parameters.crit_name = "cHedge(cEI,cLCB,cExpReturn,cOptimisticSampling)";
-  parameters.epsilon = 0.0;*/
-
-  ExampleOneD opt(dim,parameters);
-  vectord result(dim);
-  opt.optimize(result);
-  
-  std::cout << "Result:" << result << std::endl;
-  opt.printOptimal();
-
-  return 0;
-}

app/testchol.cpp

-/** -*- c++ -*- \file cholesky_test.cpp \brief test cholesky decomposition */
-/*
- -   begin                : 2005-08-24
- -   copyright            : (C) 2005 by Gunter Winkler, Konstantin Kutzkow
- -   email                : guwi17@gmx.de
-
-    This library is free software; you can redistribute it and/or
-    modify it under the terms of the GNU Lesser General Public
-    License as published by the Free Software Foundation; either
-    version 2.1 of the License, or (at your option) any later version.
-
-    This library 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
-    Lesser General Public License for more details.
-
-    You should have received a copy of the GNU Lesser General Public
-    License along with this library; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-
-*/
-
-#include <cassert>
-#include <limits>
-
-#include <boost/timer.hpp>
-
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/vector_proxy.hpp>
-
-#include <boost/numeric/ublas/triangular.hpp>
-#include <boost/numeric/ublas/banded.hpp>
-#include <boost/numeric/ublas/symmetric.hpp>
-
-#include <boost/numeric/ublas/matrix.hpp>
-#include <boost/numeric/ublas/matrix_proxy.hpp>
-
-#include <boost/numeric/ublas/io.hpp>
-
-#include "cholesky.hpp"
-
-namespace ublas = boost::numeric::ublas;
-
-/** \brief make a immutable triangular adaptor from a matrix
- *
- * \usage: 
-<code>
- A = triangular< lower >(B);
- A = triangular(B, lower());
-</code>
- */
-template < class TYPE, class MATRIX >
-ublas::triangular_adaptor<const MATRIX, TYPE>
-triangular(const MATRIX & A, const TYPE& uplo = TYPE())
-{
-  return ublas::triangular_adaptor<const MATRIX, TYPE>(A);
-}
-
-/** \brief make a immutable banded adaptor from a matrix
- *
- * \usage: 
-<code>
- A = banded(B, lower, upper);
-</code>
- */
-template < class MATRIX >
-ublas::banded_adaptor<const MATRIX>
-banded(const MATRIX & A, const size_t lower, const size_t upper)
-{
-  return ublas::banded_adaptor<const MATRIX>(A, lower, upper);
-}
-
-/** \brief make a immutable symmetric adaptor from a matrix
- *
- * \usage: 
-<code>
- A = symmetric< lower >(B);
- A = symmetric(B, lower());
-</code>
- */
-template < class TYPE, class MATRIX >
-ublas::symmetric_adaptor<const MATRIX, TYPE>
-symmetric(const MATRIX & A, const TYPE& uplo = TYPE())
-{
-  return ublas::symmetric_adaptor<const MATRIX, TYPE>(A);
-}
-
-
-/** \brief fill lower triangular matrix L 
- */
-template < class MATRIX >
-void fill_symm(MATRIX & L, const size_t bands = std::numeric_limits<size_t>::max() )
-{
-  typedef typename MATRIX::size_type size_type;
-  
-  assert(L.size1() == L.size2());
-
-  size_type size = L.size1();
-  for (size_type i=0; i<size; i++) {
-    for (size_type j = ((i>bands)?(i-bands):0); j<i; j++) {
-      L(i,j) = 1 + (1.0 + i)/(1.0 + j) + 1.0 / (0.5 + i - j);
-    }
-    L(i,i) = 1+i+size;
-  }
-
-  return;
-}
-
-int main(int argc, char * argv[] )
-{
-  size_t size = 10;
-  if (argc > 1)
-    size = ::atoi (argv [1]);
-  boost::timer  t1;
-  double pr, de, sv;
-
-  typedef double DBL;
-  typedef ublas::row_major  ORI;
-  {
-    // use dense matrix
-    ublas::matrix<DBL, ORI> A (size, size);
-    ublas::matrix<DBL, ORI> T (size, size);
-    ublas::matrix<DBL, ORI> L (size, size);
-
-    A = ublas::zero_matrix<DBL>(size, size);
-
-    ublas::vector<DBL> b (size);
-    ublas::vector<DBL> x (size);
-    ublas::vector<DBL> y (size);
-
-    std::fill(y.begin(), y.end(), 1.0);
-
-    fill_symm(T);
-    t1.restart();
-    A = ublas::prod(T, trans(T));
-    pr = t1.elapsed();
-
-    b = prod( A, y);
-
-    t1.restart();
-    size_t res = cholesky_decompose(A, L);
-    de = t1.elapsed();
-
-    t1.restart();
-    x = b;
-    cholesky_solve(L, x, ublas::lower());
-    sv = t1.elapsed();
-
-    std::cout << res << ": " 
-              << ublas::norm_inf(L-T) << " "
-              << ublas::norm_2(x-y) << " "
-              << " (deco: " << de << " sec)"
-              << " (prod: " << pr << " sec)"
-              << " (solve: " << sv << " sec)"
-              << " " << size
-              << std::endl;
-  }
-
-  {
-    // use dense triangular matrices
-    ublas::triangular_matrix<DBL, ublas::lower, ORI> A (size, size);
-    ublas::triangular_matrix<DBL, ublas::lower, ORI> T (size, size);
-    ublas::triangular_matrix<DBL, ublas::lower, ORI> L (size, size);
-    
-    A = ublas::zero_matrix<DBL> (size, size) ;
-    A = triangular<ublas::lower>( ublas::zero_matrix<DBL> (size, size) );
-    A = triangular( ublas::zero_matrix<DBL> (size, size), ublas::lower() );
-
-    ublas::vector<DBL> b (size);
-    ublas::vector<DBL> x (size);
-    ublas::vector<DBL> y (size);
-
-    std::fill(y.begin(), y.end(), 1.0);
-
-    fill_symm(T);
-    t1.restart();
-    A = triangular<ublas::lower>( ublas::prod(T, trans(T)) );
-    pr = t1.elapsed();
-
-    b = prod( symmetric<ublas::lower>(A), y);
-
-    t1.restart();
-    size_t res = cholesky_decompose(A, L);
-    de = t1.elapsed();
-
-    t1.restart();
-    x = b;
-    cholesky_solve(L, x, ublas::lower());
-    sv = t1.elapsed();
-
-    std::cout << res << ": " 
-              << ublas::norm_inf(L-T) << " "
-              << ublas::norm_2(x-y) << " "
-              << " (deco: " << de << " sec)"
-              << " (prod: " << pr << " sec)"
-              << " (solve: " << sv << " sec)"
-              << " " << size
-              << std::endl;
-//     std::cout << L << std::endl;
-//     std::cout << b << std::endl;
-//     std::cout << x << std::endl;
-//     std::cout << y << std::endl;
-//     std::cout << (b - prod(symmetric<ublas::lower>(A), x)) << std::endl;
-  }
-
-  {
-    // use banded matrices matrix
-    typedef ublas::banded_matrix<DBL, ORI> MAT;
-
-    size_t bands = std::min<size_t>( size, 50 );
-    MAT A (size, size, bands, 0);
-    MAT T (size, size, bands, 0);
-    MAT L (size, size, bands, 0);
-    
-    A = ublas::zero_matrix<DBL> (size, size) ;
-    A = banded( ublas::zero_matrix<DBL> (size, size), bands, 0 );
-
-    ublas::vector<DBL> b (size);
-    ublas::vector<DBL> x (size);
-    ublas::vector<DBL> y (size);
-
-    std::fill(y.begin(), y.end(), 1.0);
-
-    fill_symm(T, bands);
-    t1.restart();
-    A = banded( ublas::prod(T, trans(T)), bands, 0 );
-    pr = t1.elapsed();
-
-    b = prod( symmetric<ublas::lower>(A), y);
-
-    t1.restart();
-    size_t res = cholesky_decompose(A, L);
-    de = t1.elapsed();
-
-    t1.restart();
-    x = b;
-    cholesky_solve(L, x, ublas::lower());
-    sv = t1.elapsed();
-
-    std::cout << res << ": " 
-              << ublas::norm_inf(L-T) << " "
-              << ublas::norm_2(x-y) << " "
-              << " (deco: " << de << " sec)"
-              << " (prod: " << pr << " sec)"
-              << " (solve: " << sv << " sec)"
-              << " " << size
-              << std::endl;
-  }
-
-  return EXIT_SUCCESS;
-}
-
-
-/****************
-
-(note: dense + tria + 50 banded)
-
-$ g++-4.0 -I $HOME/include -o cholesky_test cholesky_test.cpp -Wall -g -O2 -DNDEBUG
-
-(column major, double)
-
-0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
-0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0.01 sec) (solve: 0 sec) 100
-0: 1.97176e-13 9.0801e-15  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
-0: 1.18172e-12 5.82623e-14  (deco: 0 sec) (prod: 0.06 sec) (solve: 0 sec) 200
-0: 1.18172e-12 5.82623e-14  (deco: 0.03 sec) (prod: 0.03 sec) (solve: 0 sec) 200
-0: 1.15463e-13 1.01481e-14  (deco: 0.01 sec) (prod: 0.01 sec) (solve: 0 sec) 200
-0: 4.80105e-12 2.13833e-13  (deco: 0.1 sec) (prod: 0.58 sec) (solve: 0 sec) 400
-0: 4.80105e-12 2.13833e-13  (deco: 0.22 sec) (prod: 0.22 sec) (solve: 0.01 sec) 400
-0: 6.39488e-14 1.14428e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0.01 sec) 400
-0: 1.80402e-11 3.72721e-13  (deco: 1.27 sec) (prod: 9.87 sec) (solve: 0 sec) 800
-0: 1.80402e-11 3.72721e-13  (deco: 2.16 sec) (prod: 2.13 sec) (solve: 0.04 sec) 800
-0: 4.04121e-14 1.52388e-14  (deco: 0.05 sec) (prod: 0.05 sec) (solve: 0.01 sec) 800
-0: 7.85829e-11 1.51858e-12  (deco: 33.5 sec) (prod: 279.15 sec) (solve: 0.05 sec) 1600
-0: 7.85829e-11 1.51858e-12  (deco: 21.8 sec) (prod: 22.45 sec) (solve: 0.18 sec) 1600
-0: 2.26485e-14 1.98783e-14  (deco: 0.19 sec) (prod: 0.09 sec) (solve: 0.01 sec) 1600
-
-(row major, double)
-
-0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
-0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
-0: 1.97176e-13 9.0801e-15  (deco: 0.01 sec) (prod: 0 sec) (solve: 0 sec) 100
-0: 1.18172e-12 5.82623e-14  (deco: 0.01 sec) (prod: 0.02 sec) (solve: 0 sec) 200
-0: 1.18172e-12 5.82623e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0 sec) 200
-0: 1.15463e-13 1.01481e-14  (deco: 0.01 sec) (prod: 0.01 sec) (solve: 0 sec) 200
-0: 4.80105e-12 2.13833e-13  (deco: 0.03 sec) (prod: 0.21 sec) (solve: 0.01 sec) 400
-0: 4.80105e-12 2.13833e-13  (deco: 0.08 sec) (prod: 0.09 sec) (solve: 0.01 sec) 400
-0: 6.39488e-14 1.14428e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0 sec) 400
-0: 1.80402e-11 3.72721e-13  (deco: 0.37 sec) (prod: 1.51 sec) (solve: 0 sec) 800
-0: 1.80402e-11 3.72721e-13  (deco: 0.69 sec) (prod: 0.57 sec) (solve: 0.03 sec) 800
-0: 4.04121e-14 1.52388e-14  (deco: 0.06 sec) (prod: 0.04 sec) (solve: 0 sec) 800
-0: 7.85829e-11 1.51858e-12  (deco: 2.57 sec) (prod: 11.26 sec) (solve: 0.05 sec) 1600
-0: 7.85829e-11 1.51858e-12  (deco: 4.98 sec) (prod: 4.43 sec) (solve: 0.15 sec) 1600
-0: 2.26485e-14 1.98783e-14  (deco: 0.21 sec) (prod: 0.06 sec) (solve: 0.01 sec) 1600
-
-$ g++-3.3.5 -I $HOME/include -o cholesky_test cholesky_test.cpp -Wall -g -O2 -DNDEBUG
-
-0: 6.72351e-13 7.10356e-14  (deco: 0.01 sec) (prod: 0.02 sec) (solve: 0 sec) 100
-0: 6.72351e-13 7.10356e-14  (deco: 0.01 sec) (prod: 0.01 sec) (solve: 0 sec) 100
-0: 1.12355e-13 1.96079e-14  (deco: 0.01 sec) (prod: 0 sec) (solve: 0 sec) 100
-0: 2.61746e-12 2.75266e-13  (deco: 0.02 sec) (prod: 0.15 sec) (solve: 0 sec) 200
-0: 2.61746e-12 2.75266e-13  (deco: 0.05 sec) (prod: 0.07 sec) (solve: 0 sec) 200
-0: 7.32747e-14 1.38831e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0 sec) 200
-0: 1.02713e-11 1.33167e-12  (deco: 0.18 sec) (prod: 1.04 sec) (solve: 0 sec) 400
-0: 1.02713e-11 1.33167e-12  (deco: 0.39 sec) (prod: 0.4 sec) (solve: 0.01 sec) 400
-0: 3.33067e-14 1.4327e-14  (deco: 0.03 sec) (prod: 0.04 sec) (solve: 0 sec) 400
-0: 4.21783e-11 3.23897e-12  (deco: 1.55 sec) (prod: 8.1 sec) (solve: 0.01 sec) 800
-0: 4.21783e-11 3.23897e-12  (deco: 3.17 sec) (prod: 3.07 sec) (solve: 0.02 sec) 800
-0: 1.42109e-14 1.87039e-14  (deco: 0.09 sec) (prod: 0.06 sec) (solve: 0 sec) 800
-0: 1.6946e-10 9.93422e-12  (deco: 12.09 sec) (prod: 64.96 sec) (solve: 0.07 sec) 1600
-0: 1.6946e-10 9.93422e-12  (deco: 24.72 sec) (prod: 24.28 sec) (solve: 0.1 sec) 1600
-0: 1.15463e-14 2.46404e-14  (deco: 0.28 sec) (prod: 0.1 sec) (solve: 0 sec) 1600
-
-
- ****************/

app/testparser.cpp

-#include <string>
-#include <iostream>
-#include "parser.hpp"
-
-int main()
-{
-  std::string test = "One(Two, Three, Four(Five, Six))";
-  std::string one;
-  std::vector<std::string> vs;
-  bayesopt::utils::parseExpresion(test,one,vs);
-
-  std::cout << one << std::endl;
-  for(std::vector<std::string>::iterator it = vs.begin();
-      it != vs.end(); ++it)
-    std::cout << *it << std::endl;
- 
-  return 0;
-}

app/testrand.cpp

-#include "randgen.hpp"
-
-int main()
-{
-  randEngine eng;
-  randFloat sample(eng, realUniformDist(0,1));
-
-  for (int i = 0; i<1000; i++)
-    std::cout << sample() << std::endl;
-
-  return 1;
-}

cmake/PythonMagic.cmake

 
 include(FindPythonInterp)
 ## check python
-find_package(PythonLibs) # using PYTHON_INCLUDE_PATH instead of PYTHON_INCLUDE_DIR
+find_package(PythonLibs 2 EXACT) # using PYTHON_INCLUDE_PATH instead of PYTHON_INCLUDE_DIR
 if( NOT PYTHON_EXECUTABLE )
   # look specifically for 2.7
   find_program(PYTHON_EXECUTABLE NAMES python2.7 python27 python PATHS [HKEY_LOCAL_MACHINE\\SOFTWARE\\Python\\PythonCore\\2.7\\InstallPath])

examples/bo_branin.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/>.
+------------------------------------------------------------------------
+*/
+
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <valarray>
+#include "bayesoptcont.hpp"
+
+#ifndef M_PI
+#define M_PI           3.14159265358979323846
+#endif
+
+class ExampleBranin: public bayesopt::ContinuousModel
+{
+public:
+  ExampleBranin(size_t dim,bopt_params par):
+    ContinuousModel(dim,par) {}
+
+  double evaluateSample( const vectord& xin)
+  {
+     if (xin.size() != 2)
+      {
+	std::cout << "WARNING: This only works for 2D inputs." << std::endl
+		  << "WARNING: Using only first two components." << std::endl;
+      }
+
+    double x = xin(0) * 15 - 5;
+    double y = xin(1) * 15;
+
+    return sqr(y-(5.1/(4*sqr(M_PI)))*sqr(x)+5*x/M_PI-6)+10*(1-1/(8*M_PI))*cos(x)+10;
+  };
+
+  bool checkReachability(const vectord &query)
+  {return true;};
+
+  inline double sqr( double x ){ return x*x; };
+
+  void printOptimal()
+  {
+    vectord sv(2);  
+    sv(0) = 0.1239; sv(1) = 0.8183;
+    std::cout << "Solutions: " << sv << "->" 
+	      << evaluateSample(sv) << std::endl;
+    sv(0) = 0.5428; sv(1) = 0.1517;
+    std::cout << "Solutions: " << sv << "->" 
+	      << evaluateSample(sv) << std::endl;
+    sv(0) = 0.9617; sv(1) = 0.1650;
+    std::cout << "Solutions: " << sv << "->" 
+	      << evaluateSample(sv) << std::endl;
+  }
+
+};
+
+int main(int nargs, char *args[])
+{
+  bopt_params par = initialize_parameters_to_default();
+  par.n_iterations = 200;
+  par.n_init_samples = 50;
+  par.kernel.hp_mean[0] = 1.0;
+  par.kernel.n_hp = 1;
+  par.crit_name = "cHedge(cEI,cLCB,cPOI)";
+  
+  ExampleBranin branin(2,par);
+  vectord result(2);
+
+  branin.optimize(result);
+  std::cout << "Result:" << result << std::endl;
+  branin.printOptimal();
+
+  return 0;
+}

examples/bo_cont.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 <ctime>
+#include "bayesoptwpr.h"                 // For the C API
+#include "bayesoptcont.hpp"              // For the C++ API
+
+
+/* Function to be used for C-API testing */
+double testFunction(unsigned int n, const double *x,
+		    double *gradient, /* NULL if not needed */
+		    void *func_data)
+{
+  double f = 10.;
+  for (unsigned int i = 0; i < n; ++i)
+    {
+      f += (x[i] - .53) * (x[i] - .53);
+    }
+  return f;
+}
+
+/* Class to be used for C++-API testing */
+class ExampleQuadratic: public bayesopt::ContinuousModel
+{
+ public:
+
+  ExampleQuadratic(size_t dim,bopt_params param):
+    ContinuousModel(dim,param) {}
+
+  double evaluateSample( const vectord &Xi ) 
+  {
+    double x[100];
+    for (size_t i = 0; i < Xi.size(); ++i)
+	x[i] = Xi(i);
+    return testFunction(Xi.size(),x,NULL,NULL);
+  };
+
+
+  bool checkReachability( const vectord &query )
+  { return true; };
+ 
+};
+
+
+int main(int nargs, char *args[])
+{    
+  int n = 10;                   // Number of dimensions
+
+  // Common configuration
+  // See parameters.h for the available options
+  // If we initialize the struct with the DEFAUL_PARAMS,
+  // the we can optionally change only few of them 
+  bopt_params par = initialize_parameters_to_default();
+
+  par.kernel.hp_mean[0] = KERNEL_THETA;
+  par.kernel.hp_mean[1] = KERNEL_THETA;
+  par.kernel.hp_std[0] = 1;
+  par.kernel.hp_std[1] = 1;
+  par.kernel.n_hp = 2;
+  par.mean.coef_mean[0] = 1.0;
+  par.mean.coef_std[0] = MEAN_SIGMA;
+  par.mean.n_coef = 1;
+  par.alpha = PRIOR_ALPHA;
+  par.beta = PRIOR_BETA;
+  par.noise = DEFAULT_NOISE;
+  par.surr_name = "sStudentTProcessJef";
+  par.kernel.name = "kSum(kSEISO,kConst)";
+  par.mean.name = "mConst";
+  par.l_type = L_ML;
+  par.n_iterations = 200;       // Number of iterations
+  par.n_init_samples = 50;
+  par.verbose_level = 2;
+  /*******************************************/
+
+  clock_t start, end;
+  double diff,diff2;
+
+  std::cout << "Running C++ interface" << std::endl;
+  // Configure C++ interface
+
+  ExampleQuadratic opt(n,par);
+  vectord result(n);
+
+  // Run C++ interface
+  start = clock();
+  opt.optimize(result);
+  end = clock();
+  diff = (double)(end-start) / (double)CLOCKS_PER_SEC;
+  /*******************************************/
+
+
+  std::cout << "Running C inferface" << std::endl;
+  
+  // Configure C interface
+  double u[128], x[128], l[128], fmin[128];
+
+  for (int i = 0; i < n; ++i) {
+    l[i] = 0.;    u[i] = 1.;
+  }
+
+  // Run C interface
+  start = clock();
+  bayes_optimization(n,&testFunction,NULL,l,u,x,fmin,par);
+  end = clock();
+  diff2 = (double)(end-start) / (double)CLOCKS_PER_SEC;
+  /*******************************************/
+
+
+  // Results
+  std::cout << "Final result C++: " << result << std::endl;
+  std::cout << "Elapsed time in C++: " << diff << " seconds" << std::endl;
+
+  std::cout << "Final result C: [" << n <<"](" << x[0];
+  for (int i = 1; i < n; ++i )
+    std::cout << "," << x[i];
+  
+  std::cout << ")" << std::endl;
+  std::cout << "Elapsed time in C: " << diff2 << " seconds" << std::endl;
+
+}
+

examples/bo_disc.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 "bayesoptwpr.h"               // For the C API
+#include "bayesoptdisc.hpp"            // For the C++ API
+#include "lhs.hpp"
+
+
+/* Function to be used for C-API testing */
+double testFunction(unsigned int n, const double *x,
+		    double *gradient, /* NULL if not needed */
+		    void *func_data)
+{
+  double f = 10.;
+  for (unsigned int i = 0; i < n; ++i)
+    {
+      f += (x[i] - .53) * (x[i] - .53);
+    }
+  return f;
+}
+
+/* Class to be used for C++-API testing */
+class ExampleDisc: public bayesopt::DiscreteModel
+{
+ public:
+
+  ExampleDisc(const vecOfvec & validSet, bopt_params param):
+    DiscreteModel(validSet,param) {}
+
+  double evaluateSample( const vectord &Xi ) 
+  {
+    double x[100];
+    for (size_t i = 0; i < Xi.size(); ++i)
+	x[i] = Xi(i);
+    return testFunction(Xi.size(),x,NULL,NULL);
+  };
+
+  bool checkReachability( const vectord &query )
+  { return true; }; 
+};
+
+
+int main(int nargs, char *args[])
+{    
+  int n = 6;                   // Number of dimensions
+
+  // Common configuration
+  // See parameters.h for the available options
+  // If we initialize the struct with the DEFAUL_PARAMS,
+  // the we can optionally change only few of them 
+  bopt_params par = initialize_parameters_to_default();
+
+  par.kernel.hp_mean[0] = KERNEL_THETA;
+  par.kernel.hp_std[0] = 1.0;
+  par.kernel.n_hp = 1;
+  par.mean.coef_mean[0] = 0.0;
+  par.mean.coef_std[0] = MEAN_SIGMA;
+  par.mean.n_coef = 1;
+  par.noise = DEFAULT_NOISE;
+  par.surr_name = "sStudentTProcessJef";
+  par.n_iterations = 20;       // Number of iterations
+  par.n_init_samples = 20;
+  /*******************************************/
+  
+  size_t nPoints = 1000;
+
+  randEngine mtRandom;
+  matrixd xPoints(nPoints,n);
+  vecOfvec xP;
+
+  //Thanks to the quirks of Visual Studio, the following expression is invalid,
+  //so we have to replace it by a literal.
+  //WARNING: Make sure to update the size of the array if the number of points
+  //or dimensions change.
+#ifdef _MSC_VER
+  double xPointsArray[6000];
+#else
+  const size_t nPinArr = n*nPoints;
+  double xPointsArray[nPinArr];
+#endif
+  
+  bayesopt::utils::lhs(xPoints,mtRandom);
+
+  for(size_t i = 0; i<nPoints; ++i)
+    {
+      vectord point = row(xPoints,i);  
+      xP.push_back(point);
+      for(size_t j=0; j<n; ++j)
+	{
+	  xPointsArray[i*n+j] = point(j);	  
+	}
+    }
+    
+
+  
+  // Run C++ interface
+  std::cout << "Running C++ interface" << std::endl;
+  ExampleDisc opt(xP,par);
+  vectord result(n);
+  opt.optimize(result);
+
+  
+  // Run C interface
+  std::cout << "Running C interface" << std::endl;
+  double x[128], fmin[128];
+  bayes_optimization_disc(n, &testFunction, NULL, xPointsArray, nPoints,
+			  x, fmin, par);
+
+  // Find the optimal value
+  size_t min = 0;
+  double minvalue = opt.evaluateSample(row(xPoints,min));
+  for(size_t i = 1; i<nPoints; ++i)
+    {
+      vectord point = row(xPoints,i);  
+      if (opt.evaluateSample(point) < minvalue)
+	{
+	  min = i;
+	  minvalue = opt.evaluateSample(row(xPoints,min));
+	  std::cout << i << "," << minvalue << std::endl;
+	}
+    }
+
+  std::cout << "Final result C++: " << result << std::endl;
+  std::cout << "Final result C: (";
+  for (int i = 0; i < n; i++ )
+    std::cout << x[i] << ", ";
+  std::cout << ")" << std::endl;
+  std::cout << "Optimal: " << row(xPoints,min) << std::endl;
+}

examples/bo_display.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 <valarray>
+#include "bayesoptcont.hpp"
+
+class ExampleOneD: public bayesopt::ContinuousModel
+{
+public:
+  ExampleOneD(size_t dim, bopt_params par):
+    ContinuousModel(dim,par) {}
+
+  double evaluateSample(const vectord& xin)
+  {
+    if (xin.size() > 1)
+      {
+	std::cout << "WARNING: This only works for 1D inputs." << std::endl
+		  << "WARNING: Using only first component." << std::endl;
+      }
+
+    double x = xin(0);
+    return sqr(x-0.3) + sin(20*x)*0.2;
+  };
+
+  bool checkReachability(const vectord &query)
+  {return true;};
+
+  inline double sqr( double x ){ return x*x; }
+
+  void printOptimal()
+  {
+    std::cout << "Optimal:" << 0.5 << std::endl;
+  }
+
+};
+
+//#include "unistd.h"
+//using namespace std;
+#include "matplotpp.h"  
+
+using namespace bayesopt;
+
+int is_run=0;
+int is_step=0;
+size_t state_ii = 0;
+BayesOptBase* GLOBAL_MODEL;
+std::vector<double> lx,ly;
+
+class MP :public MatPlot{ 
+void DISPLAY(){
+  size_t nruns = GLOBAL_MODEL->getParameters()->n_iterations;
+  if ((is_run) && (state_ii < nruns))
+    {
+      ++state_ii;
+      GLOBAL_MODEL->stepOptimization(state_ii); 
+      vectord last(1);
+      double res = GLOBAL_MODEL->getSurrogateModel()->getLastSample(last);
+      ly.push_back(res);
+      lx.push_back(last(0));
+	  
+      if (is_step) { is_run = 0; is_step = 0; }
+    }
+    
+  int n=1000;
+  std::vector<double> x,y,z,su,sl,c;
+  x=linspace(0,1,n);
+  y = x; z = x; su = x; sl = x; c= x;
+  vectord q(1);
+  for(int i=0;i<n;++i)
+    {
+      q(0) = x[i];
+      ProbabilityDistribution* pd = GLOBAL_MODEL->getSurrogateModel()->prediction(q);
+      y[i] = pd->getMean();
+      su[i] = y[i] + 2*pd->getStd();
+      sl[i] = y[i] - 2*pd->getStd();
+      c[i] = -GLOBAL_MODEL->evaluateCriteria(q);
+      z[i] = GLOBAL_MODEL->evaluateSample(q);
+    }
+ 
+  //plot
+  subplot(2,1,1);
+  title("press r to run and stop");
+  plot(x,y); set(3);
+  plot(lx,ly);set("k");set("*");
+  plot(x,su);set("g"); set(2);
+  plot(x,sl);set("g"); set(2);
+  plot(x,z);set("r"); set(3);
+  
+  subplot(2,1,2);
+  plot(x,c); set(3);
+}
+}mp;
+
+void display(){mp.display(); }
+void reshape(int w,int h){ mp.reshape(w,h); }
+void idle( void )
+{
+  glutPostRedisplay();
+}
+
+void mouse(int button, int state, int x, int y ){ mp.mouse(button,state,x,y); }
+void motion(int x, int y ){mp.motion(x,y); }
+void passive(int x, int y ){mp.passivemotion(x,y); }
+void keyboard(unsigned char key, int x, int y){
+    mp.keyboard(key, x, y); 
+    if(key=='r'){ if(is_run==0){is_run=1;}else{is_run=0;}}
+    if(key=='s'){ is_run=1; is_step=1; } 
+}
+
+int main(int nargs, char *args[])
+{
+  size_t dim = 1;
+  bopt_params parameters = initialize_parameters_to_default();
+  parameters.n_init_samples = 7;
+  parameters.n_iter_relearn = 0;
+  parameters.n_iterations = 300;
+  parameters.verbose_level = 2;
+
+  // Surrogate models
+  //  parameters.surr_name = "sStudentTProcessNIG";
+  parameters.surr_name = "sGaussianProcessNormal";
+
+  // Criterion model
+  // parameters.crit_name = "cAopt";
+  // parameters.n_crit_params = 0;
+
+  parameters.crit_name = "cEI";
+  parameters.crit_params[0] = 1;
+  parameters.n_crit_params = 1;
+
+  // parameters.crit_name = "cLCB";
+  // parameters.crit_params[0] = 5;
+  // parameters.n_crit_params = 1;
+
+  // Kernel models
+  // parameters.kernel.name = "kSum(kPoly3,kRQISO)";
+  // double mean[128] = {1, 1, 1, 1};
+  // double std[128] = {10, 10, 10, 10};
+  // size_t nhp = 4;
+  // memcpy(parameters.kernel.hp_mean, mean, nhp * sizeof(double));
+  // memcpy(parameters.kernel.hp_std,std, nhp * sizeof(double));
+  // parameters.kernel.n_hp = nhp;
+
+  parameters.kernel.name = "kMaternISO3";
+  parameters.kernel.hp_mean[0] = 1;
+  parameters.kernel.hp_std[0] = 5;
+  parameters.kernel.n_hp = 1;
+
+  state_ii = 0;
+
+  ExampleOneD* opt = new ExampleOneD(dim,parameters);
+  vectord result(dim);
+  GLOBAL_MODEL = opt;
+  opt->initializeOptimization();
+  vectord last(1);
+  size_t n_points = GLOBAL_MODEL->getSurrogateModel()->getNSamples();
+  for (size_t i = 0; i<n_points;++i)
+    {
+      double res = GLOBAL_MODEL->getSurrogateModel()->getSample(i,last);
+      ly.push_back(res);
+      lx.push_back(last(0));
+    }
+
+  glutInit(&nargs, args);
+  glutCreateWindow(50,50,800,650);
+  glutDisplayFunc( display );
+  glutReshapeFunc( reshape );
+  glutIdleFunc( idle );
+  glutMotionFunc( motion );
+  glutMouseFunc( mouse );
+  glutPassiveMotionFunc(passive);    
+  glutKeyboardFunc( keyboard );        
+  glutMainLoop();    
+
+  delete opt;
+  return 0;
+}

examples/bo_oned.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 <valarray>
+#include "bayesoptcont.hpp"
+
+class ExampleOneD: public bayesopt::ContinuousModel
+{
+public:
+  ExampleOneD(size_t dim, bopt_params par):
+    ContinuousModel(dim,par) {}
+
+  double evaluateSample(const vectord& xin)
+  {
+    if (xin.size() > 1)
+      {
+	std::cout << "WARNING: This only works for 1D inputs." << std::endl
+		  << "WARNING: Using only first component." << std::endl;
+      }
+
+    double x = xin(0);
+    return sqr(x-0.5) + sin(20*x)*0.2;
+  };
+
+  bool checkReachability(const vectord &query)
+  {return true;};
+
+  inline double sqr( double x ){ return x*x; }
+
+  void printOptimal()
+  {
+    std::cout << "Optimal:" << 0.5 << std::endl;
+  }
+
+};
+
+int main(int nargs, char *args[])
+{
+  size_t dim = 1;
+  bopt_params parameters = initialize_parameters_to_default();
+  parameters.n_init_samples = 10;
+  parameters.n_iterations = 300;
+  parameters.surr_name = "sGaussianProcessML";
+  /*  parameters.kernel.hp_mean[0] = 1.0;
+  parameters.kernel.hp_std[0] = 100.0;
+  parameters.kernel.n_hp = 1;
+  parameters.crit_name = "cHedge(cEI,cLCB,cExpReturn,cOptimisticSampling)";
+  parameters.epsilon = 0.0;*/
+
+  ExampleOneD opt(dim,parameters);
+  vectord result(dim);
+  opt.optimize(result);
+  
+  std::cout << "Result:" << result << std::endl;
+  opt.printOptimal();
+
+  return 0;
+}

examples/testchol.cpp

+/** -*- c++ -*- \file cholesky_test.cpp \brief test cholesky decomposition */
+/*
+ -   begin                : 2005-08-24
+ -   copyright            : (C) 2005 by Gunter Winkler, Konstantin Kutzkow
+ -   email                : guwi17@gmx.de
+
+    This library is free software; you can redistribute it and/or
+    modify it under the terms of the GNU Lesser General Public
+    License as published by the Free Software Foundation; either
+    version 2.1 of the License, or (at your option) any later version.
+
+    This library 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
+    Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public
+    License along with this library; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+*/
+
+#include <cassert>
+#include <limits>
+
+#include <boost/timer.hpp>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+
+#include <boost/numeric/ublas/triangular.hpp>
+#include <boost/numeric/ublas/banded.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
+
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+
+#include <boost/numeric/ublas/io.hpp>
+
+#include "cholesky.hpp"
+
+namespace ublas = boost::numeric::ublas;
+
+/** \brief make a immutable triangular adaptor from a matrix
+ *
+ * \usage: 
+<code>
+ A = triangular< lower >(B);
+ A = triangular(B, lower());
+</code>
+ */
+template < class TYPE, class MATRIX >
+ublas::triangular_adaptor<const MATRIX, TYPE>
+triangular(const MATRIX & A, const TYPE& uplo = TYPE())
+{
+  return ublas::triangular_adaptor<const MATRIX, TYPE>(A);
+}
+
+/** \brief make a immutable banded adaptor from a matrix
+ *
+ * \usage: 
+<code>
+ A = banded(B, lower, upper);
+</code>
+ */
+template < class MATRIX >
+ublas::banded_adaptor<const MATRIX>
+banded(const MATRIX & A, const size_t lower, const size_t upper)
+{
+  return ublas::banded_adaptor<const MATRIX>(A, lower, upper);
+}
+
+/** \brief make a immutable symmetric adaptor from a matrix
+ *
+ * \usage: 
+<code>
+ A = symmetric< lower >(B);
+ A = symmetric(B, lower());
+</code>
+ */
+template < class TYPE, class MATRIX >
+ublas::symmetric_adaptor<const MATRIX, TYPE>
+symmetric(const MATRIX & A, const TYPE& uplo = TYPE())
+{
+  return ublas::symmetric_adaptor<const MATRIX, TYPE>(A);
+}
+
+
+/** \brief fill lower triangular matrix L 
+ */
+template < class MATRIX >
+void fill_symm(MATRIX & L, const size_t bands = std::numeric_limits<size_t>::max() )
+{
+  typedef typename MATRIX::size_type size_type;
+  
+  assert(L.size1() == L.size2());
+
+  size_type size = L.size1();
+  for (size_type i=0; i<size; i++) {
+    for (size_type j = ((i>bands)?(i-bands):0); j<i; j++) {
+      L(i,j) = 1 + (1.0 + i)/(1.0 + j) + 1.0 / (0.5 + i - j);
+    }
+    L(i,i) = 1+i+size;
+  }
+
+  return;
+}
+
+int main(int argc, char * argv[] )
+{
+  size_t size = 10;
+  if (argc > 1)
+    size = ::atoi (argv [1]);
+  boost::timer  t1;
+  double pr, de, sv;
+
+  typedef double DBL;
+  typedef ublas::row_major  ORI;
+  {
+    // use dense matrix
+    ublas::matrix<DBL, ORI> A (size, size);
+    ublas::matrix<DBL, ORI> T (size, size);
+    ublas::matrix<DBL, ORI> L (size, size);
+
+    A = ublas::zero_matrix<DBL>(size, size);
+
+    ublas::vector<DBL> b (size);
+    ublas::vector<DBL> x (size);
+    ublas::vector<DBL> y (size);
+
+    std::fill(y.begin(), y.end(), 1.0);
+
+    fill_symm(T);
+    t1.restart();
+    A = ublas::prod(T, trans(T));
+    pr = t1.elapsed();
+
+    b = prod( A, y);
+
+    t1.restart();
+    size_t res = cholesky_decompose(A, L);
+    de = t1.elapsed();
+
+    t1.restart();
+    x = b;
+    cholesky_solve(L, x, ublas::lower());
+    sv = t1.elapsed();
+
+    std::cout << res << ": " 
+              << ublas::norm_inf(L-T) << " "
+              << ublas::norm_2(x-y) << " "
+              << " (deco: " << de << " sec)"
+              << " (prod: " << pr << " sec)"
+              << " (solve: " << sv << " sec)"
+              << " " << size
+              << std::endl;
+  }
+
+  {
+    // use dense triangular matrices
+    ublas::triangular_matrix<DBL, ublas::lower, ORI> A (size, size);
+    ublas::triangular_matrix<DBL, ublas::lower, ORI> T (size, size);
+    ublas::triangular_matrix<DBL, ublas::lower, ORI> L (size, size);
+    
+    A = ublas::zero_matrix<DBL> (size, size) ;
+    A = triangular<ublas::lower>( ublas::zero_matrix<DBL> (size, size) );
+    A = triangular( ublas::zero_matrix<DBL> (size, size), ublas::lower() );
+
+    ublas::vector<DBL> b (size);
+    ublas::vector<DBL> x (size);
+    ublas::vector<DBL> y (size);
+
+    std::fill(y.begin(), y.end(), 1.0);
+
+    fill_symm(T);
+    t1.restart();
+    A = triangular<ublas::lower>( ublas::prod(T, trans(T)) );
+    pr = t1.elapsed();
+
+    b = prod( symmetric<ublas::lower>(A), y);
+
+    t1.restart();
+    size_t res = cholesky_decompose(A, L);
+    de = t1.elapsed();
+
+    t1.restart();
+    x = b;
+    cholesky_solve(L, x, ublas::lower());
+    sv = t1.elapsed();
+
+    std::cout << res << ": " 
+              << ublas::norm_inf(L-T) << " "
+              << ublas::norm_2(x-y) << " "
+              << " (deco: " << de << " sec)"
+              << " (prod: " << pr << " sec)"
+              << " (solve: " << sv << " sec)"
+              << " " << size
+              << std::endl;
+//     std::cout << L << std::endl;
+//     std::cout << b << std::endl;
+//     std::cout << x << std::endl;
+//     std::cout << y << std::endl;
+//     std::cout << (b - prod(symmetric<ublas::lower>(A), x)) << std::endl;
+  }
+
+  {
+    // use banded matrices matrix
+    typedef ublas::banded_matrix<DBL, ORI> MAT;
+
+    size_t bands = std::min<size_t>( size, 50 );
+    MAT A (size, size, bands, 0);
+    MAT T (size, size, bands, 0);
+    MAT L (size, size, bands, 0);
+    
+    A = ublas::zero_matrix<DBL> (size, size) ;
+    A = banded( ublas::zero_matrix<DBL> (size, size), bands, 0 );
+
+    ublas::vector<DBL> b (size);
+    ublas::vector<DBL> x (size);
+    ublas::vector<DBL> y (size);
+
+    std::fill(y.begin(), y.end(), 1.0);
+
+    fill_symm(T, bands);
+    t1.restart();
+    A = banded( ublas::prod(T, trans(T)), bands, 0 );
+    pr = t1.elapsed();
+
+    b = prod( symmetric<ublas::lower>(A), y);
+
+    t1.restart();
+    size_t res = cholesky_decompose(A, L);
+    de = t1.elapsed();
+
+    t1.restart();
+    x = b;
+    cholesky_solve(L, x, ublas::lower());
+    sv = t1.elapsed();
+
+    std::cout << res << ": " 
+              << ublas::norm_inf(L-T) << " "
+              << ublas::norm_2(x-y) << " "
+              << " (deco: " << de << " sec)"
+              << " (prod: " << pr << " sec)"
+              << " (solve: " << sv << " sec)"
+              << " " << size
+              << std::endl;
+  }
+
+  return EXIT_SUCCESS;
+}
+
+
+/****************
+
+(note: dense + tria + 50 banded)
+
+$ g++-4.0 -I $HOME/include -o cholesky_test cholesky_test.cpp -Wall -g -O2 -DNDEBUG
+
+(column major, double)
+
+0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
+0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0.01 sec) (solve: 0 sec) 100
+0: 1.97176e-13 9.0801e-15  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
+0: 1.18172e-12 5.82623e-14  (deco: 0 sec) (prod: 0.06 sec) (solve: 0 sec) 200
+0: 1.18172e-12 5.82623e-14  (deco: 0.03 sec) (prod: 0.03 sec) (solve: 0 sec) 200
+0: 1.15463e-13 1.01481e-14  (deco: 0.01 sec) (prod: 0.01 sec) (solve: 0 sec) 200
+0: 4.80105e-12 2.13833e-13  (deco: 0.1 sec) (prod: 0.58 sec) (solve: 0 sec) 400
+0: 4.80105e-12 2.13833e-13  (deco: 0.22 sec) (prod: 0.22 sec) (solve: 0.01 sec) 400
+0: 6.39488e-14 1.14428e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0.01 sec) 400
+0: 1.80402e-11 3.72721e-13  (deco: 1.27 sec) (prod: 9.87 sec) (solve: 0 sec) 800
+0: 1.80402e-11 3.72721e-13  (deco: 2.16 sec) (prod: 2.13 sec) (solve: 0.04 sec) 800
+0: 4.04121e-14 1.52388e-14  (deco: 0.05 sec) (prod: 0.05 sec) (solve: 0.01 sec) 800
+0: 7.85829e-11 1.51858e-12  (deco: 33.5 sec) (prod: 279.15 sec) (solve: 0.05 sec) 1600
+0: 7.85829e-11 1.51858e-12  (deco: 21.8 sec) (prod: 22.45 sec) (solve: 0.18 sec) 1600
+0: 2.26485e-14 1.98783e-14  (deco: 0.19 sec) (prod: 0.09 sec) (solve: 0.01 sec) 1600
+
+(row major, double)
+
+0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
+0: 3.05533e-13 1.78207e-14  (deco: 0 sec) (prod: 0 sec) (solve: 0 sec) 100
+0: 1.97176e-13 9.0801e-15  (deco: 0.01 sec) (prod: 0 sec) (solve: 0 sec) 100
+0: 1.18172e-12 5.82623e-14  (deco: 0.01 sec) (prod: 0.02 sec) (solve: 0 sec) 200
+0: 1.18172e-12 5.82623e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0 sec) 200
+0: 1.15463e-13 1.01481e-14  (deco: 0.01 sec) (prod: 0.01 sec) (solve: 0 sec) 200
+0: 4.80105e-12 2.13833e-13  (deco: 0.03 sec) (prod: 0.21 sec) (solve: 0.01 sec) 400
+0: 4.80105e-12 2.13833e-13  (deco: 0.08 sec) (prod: 0.09 sec) (solve: 0.01 sec) 400
+0: 6.39488e-14 1.14428e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0 sec) 400
+0: 1.80402e-11 3.72721e-13  (deco: 0.37 sec) (prod: 1.51 sec) (solve: 0 sec) 800
+0: 1.80402e-11 3.72721e-13  (deco: 0.69 sec) (prod: 0.57 sec) (solve: 0.03 sec) 800
+0: 4.04121e-14 1.52388e-14  (deco: 0.06 sec) (prod: 0.04 sec) (solve: 0 sec) 800
+0: 7.85829e-11 1.51858e-12  (deco: 2.57 sec) (prod: 11.26 sec) (solve: 0.05 sec) 1600
+0: 7.85829e-11 1.51858e-12  (deco: 4.98 sec) (prod: 4.43 sec) (solve: 0.15 sec) 1600
+0: 2.26485e-14 1.98783e-14  (deco: 0.21 sec) (prod: 0.06 sec) (solve: 0.01 sec) 1600
+
+$ g++-3.3.5 -I $HOME/include -o cholesky_test cholesky_test.cpp -Wall -g -O2 -DNDEBUG
+
+0: 6.72351e-13 7.10356e-14  (deco: 0.01 sec) (prod: 0.02 sec) (solve: 0 sec) 100
+0: 6.72351e-13 7.10356e-14  (deco: 0.01 sec) (prod: 0.01 sec) (solve: 0 sec) 100
+0: 1.12355e-13 1.96079e-14  (deco: 0.01 sec) (prod: 0 sec) (solve: 0 sec) 100
+0: 2.61746e-12 2.75266e-13  (deco: 0.02 sec) (prod: 0.15 sec) (solve: 0 sec) 200
+0: 2.61746e-12 2.75266e-13  (deco: 0.05 sec) (prod: 0.07 sec) (solve: 0 sec) 200
+0: 7.32747e-14 1.38831e-14  (deco: 0.02 sec) (prod: 0.02 sec) (solve: 0 sec) 200
+0: 1.02713e-11 1.33167e-12  (deco: 0.18 sec) (prod: 1.04 sec) (solve: 0 sec) 400
+0: 1.02713e-11 1.33167e-12  (deco: 0.39 sec) (prod: 0.4 sec) (solve: 0.01 sec) 400
+0: 3.33067e-14 1.4327e-14  (deco: 0.03 sec) (prod: 0.04 sec) (solve: 0 sec) 400
+0: 4.21783e-11 3.23897e-12  (deco: 1.55 sec) (prod: 8.1 sec) (solve: 0.01 sec) 800
+0: 4.21783e-11 3.23897e-12  (deco: 3.17 sec) (prod: 3.07 sec) (solve: 0.02 sec) 800
+0: 1.42109e-14 1.87039e-14  (deco: 0.09 sec) (prod: 0.06 sec) (solve: 0 sec) 800
+0: 1.6946e-10 9.93422e-12  (deco: 12.09 sec) (prod: 64.96 sec) (solve: 0.07 sec) 1600
+0: 1.6946e-10 9.93422e-12  (deco: 24.72 sec) (prod: 24.28 sec) (solve: 0.1 sec) 1600
+0: 1.15463e-14 2.46404e-14  (deco: 0.28 sec) (prod: 0.1 sec) (solve: 0 sec) 1600
+
+
+ ****************/

examples/testparser.cpp

+#include <string>
+#include <iostream>
+#include "parser.hpp"
+
+int main()
+{
+  std::string test = "One(Two, Three, Four(Five, Six))";
+  std::string one;
+  std::vector<std::string> vs;
+  bayesopt::utils::parseExpresion(test,one,vs);
+
+  std::cout << one << std::endl;
+  for(std::vector<std::string>::iterator it = vs.begin();
+      it != vs.end(); ++it)
+    std::cout << *it << std::endl;
+ 
+  return 0;
+}

examples/testrand.cpp

+#include "randgen.hpp"
+
+int main()
+{
+  randEngine eng;
+  randFloat sample(eng, realUniformDist(0,1));
+
+  for (int i = 0; i<1000; i++)
+    std::cout << sample() << std::endl;
+
+  return 1;
+}

python/bayesopt.cpp

-/* Generated by Cython 0.16 on Wed May 15 16:39:41 2013 */
+/* Generated by Cython 0.16 on Tue Jul  2 14:48:59 2013 */
 
 #define PY_SSIZE_T_CLEAN
 #include "Python.h"
     }
 }
 
+#ifdef __FreeBSD__
+#include <floatingpoint.h>
+#endif
+#if PY_MAJOR_VERSION < 3
+int main(int argc, char** argv) {
+#elif defined(WIN32) || defined(MS_WINDOWS)
+int wmain(int argc, wchar_t **argv) {
+#else
+static int __Pyx_main(int argc, wchar_t **argv) {
+#endif
+    /* 754 requires that FP exceptions run in "no stop" mode by default,
+     * and until C vendors implement C99's ways to control FP exceptions,
+     * Python requires non-stop mode.  Alas, some platforms enable FP
+     * exceptions by default.  Here we disable them.
+     */
+#ifdef __FreeBSD__
+    fp_except_t m;
+    m = fpgetmask();
+    fpsetmask(m & ~FP_X_OFL);
+#endif
+    if (argc && argv)
+        Py_SetProgramName(argv[0]);
+    Py_Initialize();
+    if (argc && argv)
+        PySys_SetArgv(argc, argv);
+    { /* init module 'bayesopt' as '__main__' */
+      PyObject* m = NULL;
+      __pyx_module_is_main_bayesopt = 1;
+      #if PY_MAJOR_VERSION < 3
+          initbayesopt();
+      #else
+          m = PyInit_bayesopt();
+      #endif
+      if (PyErr_Occurred()) {
+          PyErr_Print(); /* This exits with the right code if SystemExit. */
+          #if PY_MAJOR_VERSION < 3
+          if (Py_FlushLine()) PyErr_Clear();