add python bindings for GPs

src/limbo/model.hpp

@@ -58,5 +58,6 @@
 #include <limbo/model/gp/kernel_mean_lf_opt.hpp>
 #include <limbo/model/gp/mean_lf_opt.hpp>
 #include <limbo/model/gp/no_lf_opt.hpp>
+#include <limbo/model/multi_gp/parallel_lf_opt.hpp>
 
 #endif

src/python/limbo.cc

@@ -0,0 +1,113 @@
+#include <iostream>
+#include <limbo/model.hpp>
+#include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <vector>
+
+//  clang++ -bundle -I /usr/local/include/eigen3 -I../ -I /usr/local/include/ -L /usr/local/lib   -std=c++14 -L/usr/local/Cellar/python/3.7.0/Frameworks/Python.framework/Versions/3.7/lib/python3.7/config-3.7m-darwin/  -lpython3.7m -ldl -framework CoreFoundation -framework CoreFoundation  `python3.7-config --cflags --libs` -lpython3.7m   ./limbo.cc -o limbo.cpython-37m-darwin.so
+
+using namespace limbo;
+
+namespace py = pybind11;
+
+struct Params {
+    struct kernel : public defaults::kernel {
+        BO_DYN_PARAM(double, noise);
+        BO_DYN_PARAM(bool, optimize_noise);
+    };
+    struct kernel_squared_exp_ard : public defaults::kernel_squared_exp_ard {
+    };
+    struct opt_rprop : public defaults::opt_rprop {
+        BO_DYN_PARAM(int, iterations);
+        BO_DYN_PARAM(double, eps_stop);
+    };
+    struct mean_constant {
+        BO_DYN_PARAM(double, constant);
+    };
+};
+
+BO_DECLARE_DYN_PARAM(double, Params::kernel, noise);
+BO_DECLARE_DYN_PARAM(bool, Params::kernel, optimize_noise);
+BO_DECLARE_DYN_PARAM(int, Params::opt_rprop, iterations);
+BO_DECLARE_DYN_PARAM(double, Params::opt_rprop, eps_stop);
+BO_DECLARE_DYN_PARAM(double, Params::mean_constant, constant);
+
+using mean_t = mean::Constant<Params>;
+using kernel_t = kernel::SquaredExpARD<Params>;
+using GP_t = model::GP<Params, kernel_t, mean_t, model::gp::KernelLFOpt<Params>>;
+using MultiGP_t = model::MultiGP<Params, model::GP, kernel_t, mean_t,
+    model::multi_gp::ParallelLFOpt<Params, model::gp::KernelLFOpt<Params>>>;
+bool verbose = false;
+
+void _set_params(py::kwargs kwargs)
+{
+    Params::kernel::set_noise(0.01);
+    Params::kernel::set_optimize_noise(false);
+    Params::opt_rprop::set_iterations(300);
+    Params::opt_rprop::set_eps_stop(0.0);
+    Params::mean_constant::set_constant(0.0);
+
+    for (auto item : kwargs) {
+        auto key = py::cast<std::string>(item.first);
+        if (key == "noise") {
+            Params::kernel::set_noise(py::cast<double>(item.second));
+        }
+        else if (key == "optimize_noise") {
+            Params::kernel::set_optimize_noise(py::cast<bool>(item.second));
+        }
+        else if (key == "iterations") {
+            Params::opt_rprop::set_iterations(py::cast<int>(item.second));
+        }
+        else if (key == "eps_stop") {
+            Params::opt_rprop::set_eps_stop(py::cast<double>(item.second));
+        }
+        else if (key == "mean") {
+            Params::mean_constant::set_constant(py::cast<double>(item.second));
+        }
+        else if (key == "verbose") {
+            verbose = py::cast<bool>(item.second);
+        }
+        else {
+            std::cerr << "Unrecognized parameter:" << item.first << std::endl;
+        }
+    }
+    if (verbose) {
+        std::cout << "opt_rprop::iterations (iterations) => " << Params::opt_rprop::iterations()
+                  << std::endl
+                  << "opt_rprop::eps_stop (eps_stop)=> " << Params::opt_rprop::eps_stop()
+                  << std::endl
+                  << "kernel::noise (noise) => " << Params::kernel::noise() << std::endl
+                  << "WARNING: if not optimized, noise is set for all the GPs !" << std::endl
+                  << "kernel::optimize_noise (optimize_noise)=> "
+                  << Params::kernel::optimize_noise() << std::endl;
+    }
+}
+
+template <typename GP>
+GP make_gp(py::args args, py::kwargs kwargs)
+{
+    auto train_x = py::cast<std::vector<Eigen::VectorXd>>(args[0]);
+    auto train_y = py::cast<std::vector<Eigen::VectorXd>>(args[1]);
+    assert(train_x.size() == train_y.size());
+
+    _set_params(kwargs);
+
+    GP gp;
+    gp.compute(train_x, train_y, true);
+    gp.optimize_hyperparams();
+    if (verbose) {
+        std::cout << "Learning done. data points => " << train_x.size() << std::endl;
+    }
+    return gp;
+}
+
+PYBIND11_MODULE(limbo, m)
+{
+    m.doc() = "Simplified Limbo (GP only)";
+    m.def("make_gp", &make_gp<GP_t>, "Create a GP");
+    m.def("make_multi_gp", &make_gp<MultiGP_t>, "Create a Multi GP (multi-dimensional output)");
+
+    py::class_<GP_t>(m, "GP").def("query", &GP_t::query).def("get_log_lik", &GP_t::get_log_lik);
+    py::class_<MultiGP_t>(m, "MultiGP").def("query", &MultiGP_t::query);
+}

src/python/test.py

@@ -0,0 +1,35 @@
+import numpy as np
+import limbo
+from timeit import default_timer as timer
+
+train_x = [np.random.randn(1) for i in range(0, 150)]
+train_y = [np.sin(train_x[i]) + np.random.randn(1) * 0.1 for i in range(0, 150)]
+
+print("init & opt GP...")
+gp = limbo.make_gp(train_x, train_y, optimize_noise=True, iterations=300, verbose=True)
+print("log likelihood:", gp.get_log_lik())
+
+print('starting queries (1000)')
+start = timer()
+for i in range(0, 1000):
+    x = np.random.randn(6)
+    m, sigma = gp.query(x)
+
+end = timer()
+print('time for eval:', end-start)
+print('', '')
+
+print("init & opt Multi GP...")
+train_x = [np.random.randn(2) for i in range(0, 150)]
+train_y = [np.sin(train_x[i]) + np.random.randn(2) * 0.1 for i in range(0, 150)]
+gp = limbo.make_multi_gp(train_x, train_y, optimize_noise=True, iterations=50, verbose=True)
+
+print('starting queries (1000)')
+start = timer()
+for i in range(0, 1000):
+    x = np.random.randn(6)
+    m, sigma = gp.query(x)
+
+end = timer()
+print('time for eval:', end-start)
+

src/python/wscript

@@ -0,0 +1,55 @@
+#!/usr/bin/env python
+# encoding: utf-8
+#| Copyright Inria May 2015
+#| This project has received funding from the European Research Council (ERC) under
+#| the European Union's Horizon 2020 research and innovation programme (grant
+#| agreement No 637972) - see http://www.resibots.eu
+#|
+#| Contributor(s):
+#|   - Jean-Baptiste Mouret ([email protected])
+#|   - Antoine Cully ([email protected])
+#|   - Konstantinos Chatzilygeroudis ([email protected])
+#|   - Federico Allocati ([email protected])
+#|   - Vaios Papaspyros ([email protected])
+#|   - Roberto Rama ([email protected])
+#|
+#| This software is a computer library whose purpose is to optimize continuous,
+#| black-box functions. It mainly implements Gaussian processes and Bayesian
+#| optimization.
+#| Main repository: http://github.com/resibots/limbo
+#| Documentation: http://www.resibots.eu/limbo
+#|
+#| This software is governed by the CeCILL-C license under French law and
+#| abiding by the rules of distribution of free software.  You can  use,
+#| modify and/ or redistribute the software under the terms of the CeCILL-C
+#| license as circulated by CEA, CNRS and INRIA at the following URL
+#| "http://www.cecill.info".
+#|
+#| As a counterpart to the access to the source code and  rights to copy,
+#| modify and redistribute granted by the license, users are provided only
+#| with a limited warranty  and the software's author,  the holder of the
+#| economic rights,  and the successive licensors  have only  limited
+#| liability.
+#|
+#| In this respect, the user's attention is drawn to the risks associated
+#| with loading,  using,  modifying and/or developing or reproducing the
+#| software by the user in light of its specific status of free software,
+#| that may mean  that it is complicated to manipulate,  and  that  also
+#| therefore means  that it is reserved for developers  and  experienced
+#| professionals having in-depth computer knowledge. Users are therefore
+#| encouraged to load and test the software's suitability as regards their
+#| requirements in conditions enabling the security of their systems and/or
+#| data to be ensured and,  more generally, to use and operate it in the
+#| same conditions as regards security.
+#|
+#| The fact that you are presently reading this means that you have had
+#| knowledge of the CeCILL-C license and that you accept its terms.
+#|
+
+def build(bld):
+    if bld.options.pybind:
+        bld(features  = 'c cshlib pyext',
+            includes = '. .. ../../',
+            uselib   = 'EIGEN TBB BOOST PYBIND11',
+            source   = 'limbo.cc',
+            target   = 'limbo')

src/wscript

@@ -50,6 +50,8 @@
 def build(bld):
     bld.recurse('examples')
     bld.recurse('tutorials')
+    bld.recurse('python')
+
     if bld.options.tests:
         bld.recurse('tests')
 

waf_tools/pybind11.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+# encoding: utf-8
+#| Copyright Inria May 2015
+#| This project has received funding from the European Research Council (ERC) under
+#| the European Union's Horizon 2020 research and innovation programme (grant
+#| agreement No 637972) - see http://www.resibots.eu
+#|
+#| Contributor(s):
+#|   - Jean-Baptiste Mouret ([email protected])
+#|   - Antoine Cully ([email protected])
+#|   - Konstantinos Chatzilygeroudis ([email protected])
+#|   - Federico Allocati ([email protected])
+#|   - Vaios Papaspyros ([email protected])
+#|   - Roberto Rama ([email protected])
+#|
+#| This software is a computer library whose purpose is to optimize continuous,
+#| black-box functions. It mainly implements Gaussian processes and Bayesian
+#| optimization.
+#| Main repository: http://github.com/resibots/limbo
+#| Documentation: http://www.resibots.eu/limbo
+#|
+#| This software is governed by the CeCILL-C license under French law and
+#| abiding by the rules of distribution of free software.  You can  use,
+#| modify and/ or redistribute the software under the terms of the CeCILL-C
+#| license as circulated by CEA, CNRS and INRIA at the following URL
+#| "http://www.cecill.info".
+#|
+#| As a counterpart to the access to the source code and  rights to copy,
+#| modify and redistribute granted by the license, users are provided only
+#| with a limited warranty  and the software's author,  the holder of the
+#| economic rights,  and the successive licensors  have only  limited
+#| liability.
+#|
+#| In this respect, the user's attention is drawn to the risks associated
+#| with loading,  using,  modifying and/or developing or reproducing the
+#| software by the user in light of its specific status of free software,
+#| that may mean  that it is complicated to manipulate,  and  that  also
+#| therefore means  that it is reserved for developers  and  experienced
+#| professionals having in-depth computer knowledge. Users are therefore
+#| encouraged to load and test the software's suitability as regards their
+#| requirements in conditions enabling the security of their systems and/or
+#| data to be ensured and,  more generally, to use and operate it in the
+#| same conditions as regards security.
+#|
+#| The fact that you are presently reading this means that you have had
+#| knowledge of the CeCILL-C license and that you accept its terms.
+#|
+#! /usr/bin/env python
+# JB Mouret - 2018
+
+"""
+Quick n dirty pybind11 detection
+"""
+
+import os, glob, types
+import subprocess
+from waflib.Configure import conf
+
+
+def options(opt):
+    opt.add_option('--pybind11', type='string', help='path to pybind11', dest='pybind11')
+
+@conf
+def check_pybind11(conf, *k, **kw):
+    required = kw.get('required', False)
+
+    conf.start_msg('Checking for pybind11')
+    includes_check = ['/usr/include/', '/usr/local/include/']
+
+    if conf.options.pybind11:
+        includes_check = [conf.options.pybind11]
+    try:
+        res = conf.find_file('pybind11/pybind11.h', includes_check)
+        incl = res[:-len('pybind11/pybind11.h')-1]
+        conf.env.INCLUDES_PYBIND11 = [incl]
+        conf.end_msg(incl)
+    except:
+        if required:
+            conf.fatal('Not found in %s' % str(includes_check))
+        conf.end_msg('Not found in %s' % str(includes_check), 'RED')
+    return 1

wscript

@@ -74,7 +74,8 @@ def options(opt):
         opt.load('openmp')
         opt.load('nlopt')
         opt.load('libcmaes')
-        opt.load('xcode')
+        opt.load("python")
+        opt.load("pybind11")
 
         opt.add_option('--create', type='string', help='create a new exp', dest='create_exp')
         limbo.add_create_options(opt)
@@ -89,23 +90,23 @@ def options(opt):
         opt.add_option('--write_params', type='string', help='write all the default values of parameters in a file (used by the documentation system)', dest='write_params')
         opt.add_option('--regression_benchmarks', type='string', help='config file (json) to compile benchmark for regression', dest='regression_benchmarks')
         opt.add_option('--cpp14', action='store_true', default=False, help='force c++-14 compilation [--cpp14]', dest='cpp14')
+        opt.add_option('--pybind', action='store_true', default=False, help='build python bindings [--pybind]', dest='pybind')
         opt.add_option('--no-native', action='store_true', default=False, help='disable -march=native, which can cause some troubles [--no-native]', dest='no_native')
 
-
         try:
-                os.mkdir(blddir)# because this is not always created at that stage
+            os.mkdir(blddir) #because this is not always created at that stage
         except:
-                print("build dir not created (it probably already exists, this is fine)")
+            print("build dir not created (it probably already exists, this is fine)")
         opt.logger = Logs.make_logger(blddir + '/options.log', 'mylogger')
 
         for i in glob.glob('exp/*'):
-                if os.path.isdir(i):
-                    opt.start_msg('command-line options for [%s]' % i)
-                    try:
-                        opt.recurse(i)
-                        opt.end_msg(' -> OK')
-                    except WafError:
-                        opt.end_msg(' -> no options found')
+            if os.path.isdir(i):
+                opt.start_msg('command-line options for [%s]' % i)
+                try:
+                    opt.recurse(i)
+                    opt.end_msg(' -> OK')
+                except WafError:
+                    opt.end_msg(' -> no options found')
 
         opt.recurse('src/benchmarks')
 
@@ -117,9 +118,11 @@ def configure(conf):
         conf.load('sferes')
         conf.load('openmp')
         conf.load('mkl')
-        conf.load('xcode')
         conf.load('nlopt')
         conf.load('libcmaes')
+        if conf.options.pybind:
+            conf.load("python")
+            conf.load("pybind11")
 
         native_flags = "-march=native"
 
@@ -163,17 +166,31 @@ def configure(conf):
         conf.check_mkl()
         conf.check_nlopt()
         conf.check_libcmaes()
+        py_flags = ''
+        if conf.options.pybind:
+            conf.check_python_version((3, 0))
+            conf.check_python_headers(features='pyext')
+            conf.check_python_module('numpy')
+            conf.check_pybind11(required=True)
+            if conf.env.CXX_NAME in ["gcc", "g++"]:
+                # we need -fPIC in some Linux/gcc combinations
+                # also -undefined dynamic_lookup to work with Anaconda in macOS
+                py_flags = ' -fPIC'# -undefined dynamic_lookup'
+
 
         conf.env.INCLUDES_LIMBO = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) + "/src"
 
-        all_flags = common_flags + opt_flags
+        all_flags = common_flags + opt_flags + py_flags
         conf.env['CXXFLAGS'] = conf.env['CXXFLAGS'] + all_flags.split(' ')
         Logs.pprint('NORMAL', 'CXXFLAGS: %s' % conf.env['CXXFLAGS'])
+        if conf.options.pybind:
+            Logs.pprint('NORMAL', 'PYTHONDIR: %s' % conf.env['PYTHONDIR'])
+            Logs.pprint('NORMAL', 'PYTHON_PYEXT_LDFLAGS %s' % conf.env['PYEXT_LDFLAGS'])
 
         if conf.options.exp:
-                for i in conf.options.exp.split(','):
-                        Logs.pprint('NORMAL', 'configuring for exp: %s' % i)
-                        conf.recurse('exp/' + i)
+            for i in conf.options.exp.split(','):
+                Logs.pprint('NORMAL', 'configuring for exp: %s' % i)
+                conf.recurse('exp/' + i)
         conf.recurse('src/benchmarks')
         Logs.pprint('NORMAL', '')
         Logs.pprint('NORMAL', 'WHAT TO DO NOW?')