Commit a6ff5c28 authored by Jean-Baptiste Mouret's avatar Jean-Baptiste Mouret
Browse files

refactoring: new class bo_multi

parent bf39c52e
#include "limbo/parego.hpp"
#include "limbo/stat_pareto.hpp"
using namespace limbo;
......@@ -58,27 +59,6 @@ struct mop2 {
}
};
struct Pareto {
template<typename BO>
void operator()(const BO& opt) {
if (opt.iteration() % 10 != 0)
return;
auto p_model = opt.model_pareto_front(0, 1.0, 0.005);
auto p_data = opt.data_pareto_front();
std::string it = std::to_string(opt.iteration());
std::string model = it + "_pareto_model.dat";
std::string data = it + "_pareto_data.dat";
std::ofstream pareto_model(model.c_str()),
pareto_data(data.c_str());
for (auto x : p_model)
pareto_model << std::get<1>(x).transpose() << " "
<< (std::get<1>(x).array() + std::get<2>(x)).transpose() << " "
<< (std::get<1>(x).array() - std::get<2>(x)).transpose() << " "
<< std::endl;
for (auto x : p_data)
pareto_data << std::get<1>(x).transpose() << std::endl;
}
};
......@@ -89,18 +69,19 @@ int main() {
// typedef model::GP<Params, kernel_t, mean_t> gp_t;
// typedef acquisition_functions::UCB<Params, gp_t> ucb_t;
//Parego<Params, model_fun<gp_t>, acq_fun<ucb_t> > opt;
Parego<Params, stat_fun<Pareto> > opt;
Parego<Params> opt;
opt.optimize(mop2());
auto p_model = opt.model_pareto_front(0, 1.0, 0.001);
std::cout << "modeling pareto front..." << std::endl;
auto p_model = opt.model_pareto_front();
std::cout << "computing data pareto front" << std::endl;
auto p_data = opt.data_pareto_front();
std::ofstream pareto_model("mop2_pareto_model.dat"),
pareto_data("mop2_pareto_data.dat");
std::cout << "writing..." << std::endl;
for (auto x : p_model)
pareto_model << std::get<1>(x).transpose() << " "
<< (std::get<1>(x).array() + std::get<2>(x)).transpose() << " "
<< (std::get<1>(x).array() - std::get<2>(x)).transpose() << " "
<< std::endl;
for (auto x : p_data)
pareto_data << std::get<1>(x).transpose() << std::endl;
......
......@@ -13,7 +13,7 @@ def build(bld):
source = 'parego.cpp',
includes = '. .. ../../',
target = 'parego',
uselib = 'BOOST EIGEN TBB',
uselib = 'BOOST EIGEN TBB SFERES',
use = 'limbo')
obj = bld.program(features = 'cxx',
source = 'ns_ego.cpp',
......
......@@ -2,72 +2,11 @@
#define NSEGO_HPP_
#include <algorithm>
#define VERSION "xxx"
#include <sferes/phen/parameters.hpp>
#include <sferes/gen/evo_float.hpp>
#include <sferes/eval/parallel.hpp>
#include <sferes/modif/dummy.hpp>
#include <sferes/ea/nsga2.hpp>
#include "parego.hpp"
#include "bo_multi.hpp"
namespace limbo {
namespace ns_ego {
struct SferesParams {
struct evo_float {
typedef sferes::gen::evo_float::mutation_t mutation_t;
typedef sferes::gen::evo_float::cross_over_t cross_over_t;
SFERES_CONST float cross_rate = 0.5f;
SFERES_CONST float mutation_rate = 0.1f;
SFERES_CONST float eta_m = 15.0f;
SFERES_CONST float eta_c = 10.0f;
SFERES_CONST mutation_t mutation_type =
sferes::gen::evo_float::polynomial;
SFERES_CONST cross_over_t cross_over_type =
sferes::gen::evo_float::sbx;
};
struct pop {
SFERES_CONST unsigned size = 100;
SFERES_CONST unsigned nb_gen = 1000;
SFERES_CONST int dump_period = -1;
SFERES_CONST int initial_aleat = 1;
};
struct parameters {
SFERES_CONST float min = 0.0f;
SFERES_CONST float max = 1.0f;
};
};
template<typename M>
class SferesFit {
public:
SferesFit(const std::vector<M>& models) : _models(models) {
}
SferesFit() {}
const std::vector<float>& objs() const {
return _objs;
}
float obj(size_t i) const {
return _objs[i];
}
template<typename Indiv>
void eval(const Indiv& indiv) {
this->_objs.resize(_models.size());
Eigen::VectorXd v(indiv.size());
for (size_t j = 0; j < indiv.size(); ++j)
v[j] = indiv.data(j);
for (size_t i = 0; i < _models.size(); ++i)
this->_objs[i] = _models[i].mu(v);
}
protected:
std::vector<M> _models;
std::vector<float> _objs;
};
}
template <
class Params
......@@ -78,74 +17,23 @@ namespace limbo {
, class A6 = boost::parameter::void_
, class A7 = boost::parameter::void_
>
class NsEgo : public Parego<Params, A2, A3, A4, A5, A6, A7> {
class NsEgo : public BoMulti<Params, A2, A3, A4, A5, A6, A7> {
public:
typedef BoBase<Params, obs_type<Eigen::VectorXd>, A2, A3, A4, A5, A6, A7> base_t;
typedef typename base_t::obs_t obs_t;
typedef typename base_t::model_t model_t;
typedef typename base_t::inner_optimization_t inner_optimization_t;
typedef typename base_t::acquisition_function_t acquisition_function_t;
typedef std::tuple<Eigen::VectorXd, Eigen::VectorXd, double> pareto_point_t;
typedef std::vector<pareto_point_t> pareto_t;
typedef std::tuple<Eigen::VectorXd, Eigen::VectorXd, Eigen::VectorXd> pareto_point_t;
template<typename EvalFunction>
void optimize(const EvalFunction& feval, bool reset = true) {
this->_init(feval, reset);
size_t nb_objs = this->_observations[0].size();
size_t dim = this->_samples[0].size();
using namespace ns_ego;
typedef sferes::gen::EvoFloat<EvalFunction::dim, SferesParams> gen_t;
typedef sferes::phen::Parameters<gen_t, SferesFit<model_t>, SferesParams> phen_t;
typedef sferes::eval::Parallel<SferesParams> eval_t;
typedef boost::fusion::vector<> stat_t;
typedef sferes::modif::Dummy<> modifier_t;
typedef sferes::ea::Nsga2<phen_t, eval_t, stat_t, modifier_t, SferesParams> nsga2_t;
while (this->_samples.size() == 0 || this->_pursue()) {
// compute a model for each dimension
std::vector<model_t> models(nb_objs, model_t(dim));
std::vector<std::vector<double> > uni_obs(nb_objs);
for (size_t i = 0; i < this->_observations.size(); ++i)
for (size_t j = 0; j < this->_observations[i].size(); ++j)
uni_obs[j].push_back(this->_observations[i][j]);
for (size_t i = 0; i < uni_obs.size(); ++i) {
for (size_t j = 0; j < uni_obs[i].size(); ++j)
std::cout << uni_obs[i][j] << " ";
std::cout << std::endl;
models[i].compute(this->_samples, uni_obs[i], 0.0);
}
nsga2_t ea;
ea.set_fit_proto(SferesFit<model_t>(models));
ea.run();
auto pareto_front = ea.pareto_front();
std::sort(pareto_front.begin(), pareto_front.end(),
sferes::fit::compare_objs_lex());
std::ofstream ofs("pareto.dat");
Eigen::VectorXd best_v(EvalFunction::dim);
Eigen::VectorXd x(best_v.size());
double best_sigma = -1;
for (size_t i = 0; i < pareto_front.size(); ++i) {
double sigma = 0;
for (size_t j = 0; j < x.size(); ++j)
x[j] = pareto_front[i]->data(j);
for (size_t j = 0; j < models.size(); ++j)
sigma += models[j].sigma(x);
auto s = Eigen::VectorXd(2);
s << sqrt(2.0), sqrt(2.0);
Eigen::VectorXd p(2);
p << pareto_front[i]->fit().obj(0), pareto_front[i]->fit().obj(1);
ofs << p.transpose()
<< " " << (p + s * sigma).transpose()
<< " " << (p - s * sigma).transpose() << std::endl;
if (sigma > best_sigma) {
best_sigma = sigma;
best_v = x;
}
}
this->_update_models();
this->_update_pareto_front(feval);
auto pareto = this->model_pareto_front();
auto best = std::max_element(pareto.begin(), pareto.end(),
[](const pareto_point_t& x1, const pareto_point_t& x2) {
return std::get<2>(x1).sum() > std::get<2>(x2).sum();
});
Eigen::VectorXd best_v = std::get<0>(*best);
this->add_new_sample(best_v, feval(best_v));
this->_iteration++;
std::cout << this->_iteration << " | " << best_v.transpose() << std::endl;
......
......@@ -3,8 +3,7 @@
#include <algorithm>
#include "limbo.hpp"
#include "bo_base.hpp"
#include "pareto.hpp"
#include "bo_multi.hpp"
namespace limbo {
namespace defaults {
......@@ -15,17 +14,15 @@ namespace limbo {
template <
class Params
, 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_
>
class Parego : public BoBase<Params, obs_type<Eigen::VectorXd>,
A2, A3, A4, A5, A6, A7> {
class Parego : public BoMulti<Params, A3, A4, A5, A6, A7> {
public:
typedef BoBase<Params, obs_type<Eigen::VectorXd>, A2, A3, A4, A5, A6, A7> base_t;
typedef BoBase<Params, obs_type<Eigen::VectorXd>, A3, A4, A5, A6, A7> base_t;
typedef typename base_t::obs_t obs_t;
typedef typename base_t::model_t model_t;
typedef typename base_t::inner_optimization_t inner_optimization_t;
......@@ -57,97 +54,13 @@ namespace limbo {
this->_iteration++;
}
}
pareto_t data_pareto_front() const {
std::vector<double> v(this->_samples.size());
std::fill(v.begin(), v.end(), 0.0f);
return pareto::pareto_set(_pack_data(this->_samples, this->_observations, v));
}
pareto_t model_pareto_front(double min, double max, double inc) const {
assert(this->_observations.size());
size_t nb_objs = this->_observations[0].size();
size_t dim = this->_samples[0].size();
// compute a model for each dimension
std::vector<std::vector<double> > uni_obs(nb_objs);
for (size_t i = 0; i < this->_observations.size(); ++i)
for (size_t j = 0; j < this->_observations[i].size(); ++j)
uni_obs[j].push_back(this->_observations[i][j]);
std::vector<model_t> models(nb_objs, model_t(dim));
for (size_t i = 0; i < uni_obs.size(); ++i) {
for (size_t j = 0; j < uni_obs[i].size(); ++j)
std::cout << uni_obs[i][j] << " ";
std::cout << std::endl;
models[i].compute(this->_samples, uni_obs[i], 0.0);
}
std::cout << "enumerating..." << std::endl;
std::vector<Eigen::VectorXd> points;
_enumerate(dim - 1, min, max, inc, Eigen::VectorXd::Zero(dim), points);
std::cout << "evaluating models..." << std::endl;
std::vector<Eigen::VectorXd> objs;
std::vector<double> sigma;
tie(objs, sigma) = _eval_models(points, models);
std::ofstream ofs("model.dat");
for (size_t i = 0; i < objs.size(); ++i)
ofs << objs[i].transpose() << std::endl;
std::cout << "building Pareto set" << std::endl;
this->_update_pareto_front(feval);
return pareto::pareto_set(_pack_data(points, objs, sigma));
}
protected:
void _update_stats() {
boost::fusion::for_each(this->_stat, RefreshStat_f<Parego>(*this));
}
pareto_t _pack_data(const std::vector<Eigen::VectorXd>& points,
const std::vector<Eigen::VectorXd>& objs,
const std::vector<double>& sigma) const {
assert(points.size() == objs.size());
assert(sigma.size() == objs.size());
pareto_t p(points.size());
par::loop (0, p.size(), [&](size_t k) {
p[k] = std::make_tuple(points[k], objs[k], sigma[k]);
});
return p;
}
std::tuple<std::vector<Eigen::VectorXd>, std::vector<double> >
_eval_models(const std::vector<Eigen::VectorXd>& points, const std::vector<model_t>& models) const {
std::vector<Eigen::VectorXd> objs(points.size());
std::vector<double> sigma(points.size());
std::fill(sigma.begin(), sigma.end(), 0.0);
int nb_objs = models.size();
par::loop(size_t(0), points.size(), [&](size_t k) {
if (k % 10000 == 0) {
std::cout << k << " [" << points.size() << "] ";
std::cout.flush();
}
Eigen::VectorXd p(nb_objs);
for (size_t i = 0; i < p.size(); ++i) {
p[i] = models[i].mu(points[k]);
sigma[k] = std::max(sigma[k], models[i].sigma(points[k]));
}
objs[k] = p;
});
return std::make_tuple(objs, sigma);
}
void _enumerate(int dim, double min, double max, double inc,
Eigen::VectorXd p,
std::vector<Eigen::VectorXd>& res) const {
for (double x = min; x < max; x += inc) {
p[dim] = x;
if (dim == 0)
res.push_back(p);
else
_enumerate(dim - 1, min, max, inc, p, res);
}
}
std::vector<double> _scalarize_obs() {
assert(this->_observations.size() != 0);
......
......@@ -40,7 +40,7 @@ def configure(conf):
common_flags += " -DUSE_SFERES "
# release
opt_flags = common_flags + ' -O3 -msse2 -ggdb3'
opt_flags = common_flags + ' -O3 -msse2 -ggdb3 -g'
conf.env['CXXFLAGS'] = cxxflags + opt_flags.split(' ')
print conf.env['CXXFLAGS']
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment