Merge pull request #167 from resibots/constrained_ei

Bo with constraints

src/examples/experimental/constrained_ei.cpp

+#include <limbo/limbo.hpp>
+
+#include <limbo/experimental/bayes_opt/cboptimizer.hpp>
+#include <limbo/experimental/acqui/cei.hpp>
+
+using namespace limbo;
+
+struct Params {
+    struct cbayes_opt_boptimizer : public defaults::cbayes_opt_boptimizer {
+        BO_PARAM(double, noise, 0.01);
+    };
+
+    struct init_randomsampling {
+        BO_PARAM(int, samples, 10);
+    };
+
+    struct kernel_exp : public defaults::kernel_exp {
+    };
+
+    struct bayes_opt_bobase {
+        BO_PARAM(bool, stats_enabled, true);
+    };
+
+    struct stop_maxiterations {
+        BO_PARAM(int, iterations, 100);
+    };
+
+    struct acqui_cei : public defaults::acqui_cei {
+    };
+
+    struct mean_constant {
+        BO_PARAM(double, constant, 1.0);
+    };
+
+#ifdef USE_NLOPT
+    struct opt_nloptnograd : public defaults::opt_nloptnograd {
+    };
+#elif defined(USE_LIBCMAES)
+    struct opt_cmaes : public defaults::opt_cmaes {
+    };
+#else
+    struct opt_gridsearch : public defaults::opt_gridsearch {
+    };
+#endif
+};
+
+struct cosine {
+    static constexpr size_t dim_in = 1;
+    static constexpr size_t dim_out = 1;
+    static constexpr size_t nb_constraints = 1;
+
+    Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
+    {
+        Eigen::VectorXd res(2);
+        // we _maximize in [0:1]
+        Eigen::VectorXd xx = -4.0 + 8.0 * x.array();
+        res(0) = std::cos(xx.array()(0));
+
+        // testing the constraints
+        std::string feas = "infeasible";
+        res(1) = 0;
+        if (res(0) < 0.5) {
+            res(1) = 1;
+            feas = "feasible";
+        }
+        std::cout << xx(0) << ": " << res(0) << " --> " << feas << std::endl;
+        return res;
+    }
+};
+
+int main()
+{
+    tools::par::init();
+    typedef cosine func_t;
+
+    using Stop_t = boost::fusion::vector<stop::MaxIterations<Params>>;
+    using Stat_t = boost::fusion::vector<stat::Samples<Params>,
+        stat::BlSamples<Params>,
+        stat::BestObservations<Params>,
+        stat::AggregatedObservations<Params>>;
+    using Mean_t = mean::Constant<Params>;
+    using Kernel_t = kernel::Exp<Params>;
+    using GP_t = model::GP<Params, Kernel_t, Mean_t>;
+    using Constrained_GP_t = model::GP<Params, Kernel_t, Mean_t>;
+
+    using Acqui_t = experimental::acqui::CEI<Params, GP_t, Constrained_GP_t>;
+    using Init_t = init::RandomSampling<Params>;
+
+    experimental::bayes_opt::CBOptimizer<Params,
+        modelfun<GP_t>,
+        acquifun<Acqui_t>,
+        statsfun<Stat_t>,
+        initfun<Init_t>,
+        stopcrit<Stop_t>,
+        experimental::constraint_modelfun<Constrained_GP_t>>
+        opt;
+
+    opt.optimize(func_t());
+
+    size_t n = 0;
+    double best = -100;
+    for (size_t i = 0; i < opt.samples().size(); i++) {
+        Eigen::VectorXd res = func_t()(opt.samples()[i]);
+        if (res(0) > best && res(0) < 0.5)
+            best = res(0);
+        if (res(0) >= 0.5)
+            n++;
+    }
+    std::cout << "Infeasible points tested: " << n << std::endl;
+    std::cout << "Best feasible observation: " << best << std::endl;
+
+    return 0;
+}

src/examples/experimental/wscript

@@ -4,32 +4,32 @@
 #| 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 (jean-baptiste.mouret@inria.fr)
 #|   - Antoine Cully (antoinecully@gmail.com)
 #|   - Kontantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
 #|   - Federico Allocati (fede.allocati@gmail.com)
 #|   - Vaios Papaspyros (b.papaspyros@gmail.com)
-#| 
+#|
 #| 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,
@@ -40,10 +40,10 @@
 #| 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):
@@ -59,10 +59,15 @@ def build(bld):
                       target = 'parego',
                       uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT',
                       use = 'limbo')
-
     obj = bld.program(features = 'cxx',
                       source = 'multi.cpp',
                       includes = '.. ../.. ../../../',
                       target = 'multi',
                       uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT',
                       use = 'limbo')
+    obj = bld.program(features = 'cxx',
+                      source = 'constrained_ei.cpp',
+                      includes = '.. ../.. ../../../',
+                      target = 'constrained_ei',
+                      uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT',
+                      use = 'limbo')

src/limbo/experimental/acqui/cei.hpp

+//| 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 (jean-baptiste.mouret@inria.fr)
+//|   - Antoine Cully (antoinecully@gmail.com)
+//|   - Kontantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
+//|   - Federico Allocati (fede.allocati@gmail.com)
+//|   - Vaios Papaspyros (b.papaspyros@gmail.com)
+//|
+//| 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.
+//|
+#ifndef LIMBO_ACQUI_CEI_HPP
+#define LIMBO_ACQUI_CEI_HPP
+
+#include <cmath>
+#include <vector>
+#include <Eigen/Core>
+
+#include <limbo/tools/macros.hpp>
+
+namespace limbo {
+    namespace defaults {
+        struct acqui_cei {
+            /// @ingroup acqui_defaults
+            BO_PARAM(double, jitter, 0.0);
+        };
+    }
+
+    namespace experimental {
+        namespace acqui {
+            template <typename Params, typename Model, typename ConstraintModel>
+            class CEI {
+            public:
+                CEI(const Model& model, const ConstraintModel& constraint_model, int iteration = 0)
+                    : _model(model), _constraint_model(constraint_model) {}
+
+                size_t dim_in() const { return _model.dim_in(); }
+
+                size_t dim_out() const { return _model.dim_out(); }
+
+                template <typename AggregatorFunction>
+                double operator()(const Eigen::VectorXd& v, const AggregatorFunction& afun) const
+                {
+                    Eigen::VectorXd mu;
+                    double sigma_sq;
+                    std::tie(mu, sigma_sq) = _model.query(v);
+                    double sigma = std::sqrt(sigma_sq);
+
+                    // If \sigma(x) = 0 or we do not have any observation yet we return 0
+                    if (sigma < 1e-10 || _model.samples().size() < 1)
+                        return 0.0;
+
+                    // Compute constrained EI(x)
+                    // First find the best (predicted) observation so far
+                    // (We are zeroing infeasible samples subject to the constraint value)
+                    std::vector<double> rewards;
+                    for (auto s : _model.samples())
+                        rewards.push_back(afun(_model.mu(s)));
+
+                    double f_max = *std::max_element(rewards.begin(), rewards.end());
+                    // Calculate Z and \Phi(Z) and \phi(Z)
+                    double X = afun(mu) - f_max - Params::acqui_cei::jitter();
+                    double Z = X / sigma;
+                    double phi = std::exp(-0.5 * std::pow(Z, 2.0)) / std::sqrt(2.0 * M_PI);
+                    double Phi = 0.5 * std::erfc(-Z / std::sqrt(2));
+
+                    return _pf(v, afun) * (X * Phi + sigma * phi);
+                }
+
+            protected:
+                const Model& _model;
+                const ConstraintModel& _constraint_model;
+
+                template <typename AggregatorFunction>
+                double _pf(const Eigen::VectorXd& v, const AggregatorFunction& afun) const
+                {
+                    Eigen::VectorXd mu;
+                    double sigma_sq;
+                    std::tie(mu, sigma_sq) = _constraint_model.query(v);
+                    double sigma = std::sqrt(sigma_sq);
+
+                    if (sigma < 1e-10 || _constraint_model.samples().size() < 1)
+                        return 1.0;
+
+                    double Z = (afun(mu) - 1.0) / sigma;
+                    double Phi = 0.5 * std::erfc(-Z / std::sqrt(2));
+
+                    return Phi;
+                }
+            };
+        }
+    }
+}
+
+#endif

src/limbo/experimental/bayes_opt/cboptimizer.hpp

+//| 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 (jean-baptiste.mouret@inria.fr)
+//|   - Antoine Cully (antoinecully@gmail.com)
+//|   - Kontantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
+//|   - Federico Allocati (fede.allocati@gmail.com)
+//|   - Vaios Papaspyros (b.papaspyros@gmail.com)
+//|
+//| 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.
+//|
+#ifndef LIMBO_BAYES_OPT_CBOPTIMIZER_HPP
+#define LIMBO_BAYES_OPT_CBOPTIMIZER_HPP
+
+#include <iostream>
+#include <algorithm>
+#include <iterator>
+
+#include <boost/parameter/aux_/void.hpp>
+
+#include <Eigen/Core>
+
+#include <limbo/tools/macros.hpp>
+#include <limbo/tools/random_generator.hpp>
+#include <limbo/bayes_opt/bo_base.hpp>
+#ifdef USE_NLOPT
+#include <limbo/opt/nlopt_no_grad.hpp>
+#elif defined USE_LIBCMAES
+#include <limbo/opt/cmaes.hpp>
+#else
+#include <limbo/opt/grid_search.hpp>
+#endif
+
+namespace limbo {
+    namespace defaults {
+        struct cbayes_opt_boptimizer {
+            BO_PARAM(double, noise, 1e-6);
+            BO_PARAM(int, hp_period, -1);
+        };
+    }
+
+    namespace experimental {
+        BOOST_PARAMETER_TEMPLATE_KEYWORD(constraint_modelfun)
+
+        namespace bayes_opt {
+
+            typedef boost::parameter::parameters<boost::parameter::optional<limbo::tag::acquiopt>,
+                boost::parameter::optional<limbo::tag::statsfun>,
+                boost::parameter::optional<limbo::tag::initfun>,
+                boost::parameter::optional<limbo::tag::acquifun>,
+                boost::parameter::optional<limbo::tag::stopcrit>,
+                boost::parameter::optional<limbo::tag::modelfun>,
+                boost::parameter::optional<limbo::experimental::tag::constraint_modelfun>>
+                cboptimizer_signature;
+
+            // clang-format off
+        /**
+        The classic Bayesian optimization algorithm.
+
+        \\rst
+        References: :cite:`brochu2010tutorial,Mockus2013`
+        \\endrst
+
+        This class takes the same template parameters as BoBase. It adds:
+        \\rst
+        +---------------------+------------+----------+---------------+
+        |type                 |typedef     | argument | default       |
+        +=====================+============+==========+===============+
+        |acqui. optimizer     |acqui_opt_t | acquiopt | see below     |
+        +---------------------+------------+----------+---------------+
+        \\endrst
+
+        The default value of acqui_opt_t is:
+        - ``opt::Cmaes<Params>`` if libcmaes was found in `waf configure`
+        - ``opt::NLOptNoGrad<Params, nlopt::GN_DIRECT_L_RAND>`` if NLOpt was found but libcmaes was not found
+        - ``opt::GridSearch<Params>`` otherwise (please do not use this: the algorithm will not work at all!)
+        */
+        template <class Params,
+          class A1 = boost::parameter::void_,
+          class A2 = boost::parameter::void_,
+          class A3 = boost::parameter::void_,
+          class A4 = boost::parameter::void_,
+          class A5 = boost::parameter::void_,
+          class A6 = boost::parameter::void_,
+          class A7 = boost::parameter::void_>
+            // clang-format on
+            class CBOptimizer : public limbo::bayes_opt::BoBase<Params, A1, A2, A3, A4, A5, A6> {
+            public:
+                // defaults
+                struct defaults {
+#ifdef USE_NLOPT
+                    typedef limbo::opt::NLOptNoGrad<Params, nlopt::GN_DIRECT_L_RAND> acquiopt_t;
+#elif defined(USE_LIBCMAES)
+                    typedef limbo::opt::Cmaes<Params> acquiopt_t;
+#else
+#warning NO NLOpt, and NO Libcmaes: the acquisition function will be optimized by a grid search algorithm (which is usually bad). Please install at least NLOpt or libcmaes to use limbo!.
+                    typedef limbo::opt::GridSearch<Params> acquiopt_t;
+#endif
+                    typedef limbo::kernel::Exp<Params> kf_t;
+                    typedef limbo::mean::Constant<Params> mean_t;
+                    typedef limbo::model::GP<Params, kf_t, mean_t> constraint_model_t;
+                };
+                /// link to the corresponding BoBase (useful for typedefs)
+                typedef limbo::bayes_opt::BoBase<Params, A1, A2, A3, A4, A5, A6> base_t;
+                typedef typename base_t::model_t model_t;
+
+                typedef typename base_t::acquisition_function_t acquisition_function_t;
+                // extract the types
+                typedef typename cboptimizer_signature::bind<A1, A2, A3, A4, A5, A6, A7>::type args;
+                typedef typename boost::parameter::binding<args, limbo::tag::acquiopt, typename defaults::acquiopt_t>::type acqui_optimizer_t;
+
+                typedef typename boost::parameter::binding<args, limbo::experimental::tag::constraint_modelfun, typename defaults::constraint_model_t>::type constraint_model_t;
+
+                /// The main function (run the Bayesian optimization algorithm)
+                template <typename StateFunction, typename AggregatorFunction = FirstElem>
+                void optimize(const StateFunction& sfun, const AggregatorFunction& afun = AggregatorFunction(), bool reset = true)
+                {
+                    _nb_constraints = StateFunction::nb_constraints;
+                    _dim_out = StateFunction::dim_out;
+
+                    this->_init(sfun, afun, reset);
+
+                    if (!this->_observations.empty()) {
+                        _split_observations();
+                        _model.compute(this->_samples, _obs[0], Eigen::VectorXd::Constant(_obs[0].size(), Params::cbayes_opt_boptimizer::noise()), this->_bl_samples, Eigen::VectorXd::Constant(this->_bl_samples.size(), Params::cbayes_opt_boptimizer::noise()));
+                        if (_nb_constraints > 0)
+                            _constraint_model.compute(this->_samples, _obs[1], Eigen::VectorXd::Constant(_obs[1].size(), Params::cbayes_opt_boptimizer::noise()), this->_bl_samples, Eigen::VectorXd::Constant(this->_bl_samples.size(), Params::cbayes_opt_boptimizer::noise()));
+                    }
+                    else {
+                        _model = model_t(StateFunction::dim_in, StateFunction::dim_out);
+                        if (_nb_constraints > 0)
+                            _constraint_model = constraint_model_t(StateFunction::dim_in, _nb_constraints);
+                    }
+
+                    acqui_optimizer_t acqui_optimizer;
+
+                    while (!this->_stop(*this, afun)) {
+                        acquisition_function_t acqui(_model, _constraint_model, this->_current_iteration);
+
+                        // we do not have gradient in our current acquisition function
+                        auto acqui_optimization =
+                            [&](const Eigen::VectorXd& x, bool g) { return opt::no_grad(acqui(x, afun)); };
+                        Eigen::VectorXd starting_point = tools::random_vector(StateFunction::dim_in);
+                        Eigen::VectorXd new_sample = acqui_optimizer(acqui_optimization, starting_point, true);
+                        bool blacklisted = !this->eval_and_add(sfun, new_sample);
+
+                        this->_update_stats(*this, afun, blacklisted);
+
+                        if (blacklisted) {
+                            _model.add_bl_sample(this->_bl_samples.back(), Params::cbayes_opt_boptimizer::noise());
+                            if (_nb_constraints > 0)
+                                _constraint_model.add_bl_sample(this->_bl_samples.back(), Params::cbayes_opt_boptimizer::noise());
+                        }
+                        else {
+                            _model.add_sample(this->_samples.back(), _obs[0].back(), Params::cbayes_opt_boptimizer::noise());
+                            if (_nb_constraints > 0)
+                                _constraint_model.add_sample(this->_samples.back(), _obs[1].back(), Params::cbayes_opt_boptimizer::noise());
+                        }
+
+                        if (Params::cbayes_opt_boptimizer::hp_period() > 0
+                            && (this->_current_iteration + 1) % Params::cbayes_opt_boptimizer::hp_period() == 0) {
+                            _model.optimize_hyperparams();
+                            if (_nb_constraints > 0)
+                                _constraint_model.optimize_hyperparams();
+                        }
+
+                        this->_current_iteration++;
+                        this->_total_iterations++;
+                    }
+                }
+
+                /// return the best observation so far (i.e. max(f(x)))
+                template <typename AggregatorFunction = FirstElem>
+                const Eigen::VectorXd& best_observation(const AggregatorFunction& afun = AggregatorFunction()) const
+                {
+                    _dim_out = _model.dim_out();
+                    _split_observations();
+
+                    std::vector<Eigen::VectorXd> obs = _feasible_observations();
+                    if (obs.size() > 0)
+                        _obs[0] = obs;
+
+                    auto rewards = std::vector<double>(_obs[0].size());
+                    std::transform(_obs[0].begin(), _obs[0].end(), rewards.begin(), afun);
+                    auto max_e = std::max_element(rewards.begin(), rewards.end());
+                    return _obs[0][std::distance(rewards.begin(), max_e)];
+                }
+
+                /// return the best sample so far (i.e. the argmax(f(x)))
+                template <typename AggregatorFunction = FirstElem>
+                const Eigen::VectorXd& best_sample(const AggregatorFunction& afun = AggregatorFunction()) const
+                {
+                    _dim_out = _model.dim_out();
+                    _split_observations();
+
+                    std::vector<Eigen::VectorXd> obs = _feasible_observations();
+                    if (obs.size() > 0)
+                        _obs[0] = obs;
+
+                    auto rewards = std::vector<double>(_obs[0].size());
+                    std::transform(_obs[0].begin(), _obs[0].end(), rewards.begin(), afun);
+                    auto max_e = std::max_element(rewards.begin(), rewards.end());
+                    return this->_samples[std::distance(rewards.begin(), max_e)];
+                }
+
+                const model_t& model() const { return _model; } // returns model for the objective function
+            protected:
+                model_t _model;
+                constraint_model_t _constraint_model;
+
+                size_t _nb_constraints;
+                mutable size_t _dim_out;
+                mutable std::vector<std::vector<Eigen::VectorXd>> _obs;
+
+                std::vector<Eigen::VectorXd> _feasible_observations() const
+                {
+                    std::vector<Eigen::VectorXd> feasible_obs;
+                    for (size_t i = 0; i < _obs[0].size(); ++i) {
+                        if (_obs[1][i].prod() > 0)
+                            feasible_obs.push_back(_obs[0][i]);
+                    }
+
+                    return feasible_obs;
+                }
+
+                void _split_observations() const
+                {
+                    _obs.clear();
+                    _obs.resize(2);
+
+                    for (size_t i = 0; i < this->_observations.size(); ++i) {
+                        assert(this->_observations[i].size() == _dim_out + _nb_constraints);
+
+                        Eigen::VectorXd vec_obj(_dim_out);
+                        for (size_t j = 0; j < _dim_out; ++j)
+                            vec_obj(j) = this->_observations[i](j);
+                        _obs[0].push_back(vec_obj);
+
+                        if (_nb_constraints > 0) {
+                            Eigen::VectorXd vec_con(_nb_constraints);
+                            for (size_t j = _dim_out, ind = 0; j < this->_observations[i].size(); ++j, ++ind)
+                                vec_con(ind) = this->_observations[i](j);
+                            _obs[1].push_back(vec_con);
+                        }
+                    }
+                }
+            };
+        }
+    }
+}
+
+#endif