Commit 696e1c40 authored by Aneoshun's avatar Aneoshun
Browse files

generalization over multi-dimensional observations

parent cc1672bb
......@@ -60,7 +60,7 @@ def qsub(conf_file):
# maximum execution time
#PBS -l walltime=@wall_time
# mail parameters
#PBS -m abe
# number of nodes
#PBS -l nodes=@nb_cores:ppn=@ppn
#PBS -l pmem=5200mb -l mem=5200mb
......
......@@ -243,10 +243,11 @@ cmaes_init(cmaes_t *t, /* "this" */
double dtest, trace;
static const char * version = "3.11.02.beta";
if (t->version && strcmp(version, t->version) == 0) {
/* if (t->version && strcmp(version, t->version) == 0) {
ERRORMESSAGE("cmaes_init called twice, which will lead to a memory leak, use cmaes_exit first",0,0,0);
printf("Warning: cmaes_init called twice, which will lead to a memory leak, use cmaes_exit first\n");
}
}*/
t->version = version;
/* assign_string(&t->signalsFilename, "cmaes_signals.par"); */
......
......@@ -10,8 +10,18 @@ using namespace limbo;
struct Params {
struct gp_ucb : public defaults::gp_ucb {};
struct cmaes : public defaults::cmaes {};
struct gp_auto : public defaults::gp_auto {};
struct meanconstant : public defaults::meanconstant {};
struct ucb
{
BO_PARAM(float, alpha,0.1);
};
struct kf_maternfivehalfs
{
BO_PARAM( float,sigma,1);
BO_PARAM( float, l, 0.2);
};
struct boptimizer {
BO_PARAM(double, noise, 0.001);
......@@ -28,10 +38,11 @@ struct Params {
struct fit_eval {
static constexpr size_t dim = 2;
double operator()(const Eigen::VectorXd& x) const {
double res = 0;
Eigen::VectorXd operator()(const Eigen::VectorXd& x) const {
Eigen::VectorXd res(1);
res(0) = 0;
for (int i = 0; i < x.size(); i++)
res += 1 - (x[i] - 0.3) * (x[i] - 0.3) + sin(10 * x[i]) * 0.2;
res(0) += 1 - (x[i] - 0.3) * (x[i] - 0.3) + sin(10 * x[i]) * 0.2;
return res;
}
......@@ -39,7 +50,12 @@ struct fit_eval {
int main() {
BOptimizer<Params> opt;
typedef kernel_functions::MaternFiveHalfs<Params> Kernel_t;
typedef mean_functions::MeanData<Params> Mean_t;
typedef model::GP<Params, Kernel_t, Mean_t> GP_t;
typedef acquisition_functions::UCB<Params,GP_t> Acqui_t;
BOptimizer<Params,model_fun<GP_t>, acq_fun<Acqui_t>> opt;
opt.optimize(fit_eval());
std::cout << opt.best_observation()
<< " res " << opt.best_sample().transpose()
......
//#define SHOW_TIMER
#include <boost/fusion/container/vector.hpp>
#include <boost/fusion/include/vector.hpp>
#include "limbo/limbo.hpp"
#include "limbo/inner_cmaes.hpp"
using namespace limbo;
struct Params {
struct gp_ucb : public defaults::gp_ucb {};
struct cmaes : public defaults::cmaes {};
struct meanconstant : public defaults::meanconstant {};
struct ucb
{
BO_PARAM(float, alpha,0.1);
};
struct kf_maternfivehalfs
{
BO_PARAM( float,sigma,1);
BO_PARAM( float, l, 0.2);
};
struct boptimizer {
BO_PARAM(double, noise, 0.001);
BO_PARAM(int, dump_period, 1);
};
struct init {
BO_PARAM(int, nb_samples, 5);
};
struct maxiterations {
BO_PARAM(int, n_iterations, 20);
};
};
template<typename Params, typename Model>
class UCB_multi {
public:
UCB_multi(const Model& model, int iteration = 0) : _model(model) {
}
size_t dim() const {
return _model.dim();
}
double operator()(const Eigen::VectorXd& v) const {
//double mu, sigma;
//std::tie(mu, sigma) = _model.query(v);
//return (mu + Params::ucb::alpha() * sqrt(sigma));
return (sqrt(_model.sigma(v)));
}
protected:
const Model& _model;
};
struct fit_eval {
static constexpr size_t dim = 2;
Eigen::VectorXd operator()(const Eigen::VectorXd& x) const {
Eigen::VectorXd res(2);
res(0) = 0;
for (int i = 0; i < x.size(); i++)
res(0) += 1 - (x[i] - 0.3) * (x[i] - 0.3) + sin(10 * x[i]) * 0.2;
res(1)=res(0);
return res;
}
};
int main() {
typedef kernel_functions::MaternFiveHalfs<Params> Kernel_t;
typedef mean_functions::MeanData<Params> Mean_t;
typedef model::GP<Params, Kernel_t, Mean_t> GP_t;
typedef UCB_multi<Params,GP_t> Acqui_t;
BOptimizer<Params,model_fun<GP_t>, acq_fun<Acqui_t>> opt;
opt.optimize(fit_eval());
std::cout << opt.best_observation()
<< " res " << opt.best_sample().transpose()
<< std::endl;
return 0;
}
......@@ -9,6 +9,12 @@ def build(bld):
target = 'mono_dim',
uselib = 'BOOST EIGEN TBB',
use = 'limbo')
obj = bld.program(features = 'cxx',
source = 'obs_multi.cpp',
includes = '. .. ../../',
target = 'obs_multi',
uselib = 'BOOST EIGEN TBB',
use = 'limbo')
obj = bld.program(features = 'cxx',
source = 'parego.cpp',
includes = '. .. ../../',
......
......@@ -31,9 +31,11 @@ namespace limbo {
return _model.dim();
}
double operator()(const Eigen::VectorXd& v) const {
double mu, sigma;
std::tie(mu, sigma) = _model.query(v);
return (mu + Params::ucb::alpha() * sqrt(sigma));
//double mu, sigma;
//std::tie(mu, sigma) = _model.query(v);
//return (mu + Params::ucb::alpha() * sqrt(sigma));
return (_model.mu(v)(0) + Params::ucb::alpha() * sqrt(_model.sigma(v)));
}
protected:
const Model& _model;
......@@ -55,9 +57,10 @@ namespace limbo {
return _model.dim();
}
double operator()(const Eigen::VectorXd& v) const {
double mu, sigma;
std::tie(mu, sigma) = _model.query(v);
return (mu + _beta * sqrt(sigma));
//double mu, sigma;
//std::tie(mu, sigma) = _model.query(v);
//return (mu + _beta * sqrt(sigma));
return (_model.mu(v)(0) + _beta * sqrt(_model.sigma(v)));
}
protected:
const Model& _model;
......
......@@ -109,7 +109,7 @@ namespace limbo {
typedef acquisition_functions::GP_UCB<Params, model_t> acqui_t; //4
typedef stat::Acquisitions<Params> stat_t; //5
typedef boost::fusion::vector<stopping_criterion::MaxIterations<Params> > stop_t; //6
typedef double obs_t; //7
typedef Eigen::VectorXd obs_t; //7
};
// extract the types
......
......@@ -7,6 +7,8 @@
namespace limbo {
bool compareVectorXd(Eigen::VectorXd i, Eigen::VectorXd j) { return i(0)<j(0); }
template <
class Params
, class A1 = boost::parameter::void_
......@@ -27,9 +29,11 @@ namespace limbo {
template<typename EvalFunction>
void optimize(const EvalFunction& feval, bool reset = true) {
static_assert(std::is_floating_point<obs_t>::value, "BOptimizer wants double/double for obs");
// static_assert(std::is_floating_point<obs_t>::value, "BOptimizer wants double/double for obs");
this->_init(feval, reset);
model_t model(EvalFunction::dim);
model.compute(this->_samples, this->_observations, Params::boptimizer::noise());
inner_optimization_t inner_optimization;
......@@ -37,6 +41,7 @@ namespace limbo {
acquisition_function_t acqui(model, this->_iteration);
Eigen::VectorXd new_sample = inner_optimization(acqui, acqui.dim());
this->add_new_sample(new_sample, feval(new_sample));
model.compute(this->_samples, this->_observations, Params::boptimizer::noise());
......@@ -44,8 +49,8 @@ namespace limbo {
std::cout << this->_iteration << " new point: "
<< this->_samples[this->_samples.size() - 1].transpose()
<< " value: " << this->_observations[this->_observations.size() - 1]
<< " best:" << this->best_observation()
<< " value: " << this->_observations[this->_observations.size() - 1].transpose()
<< " best:" << this->best_observation().transpose()
<< std::endl;
this->_iteration++;
......@@ -53,11 +58,11 @@ namespace limbo {
}
const obs_t& best_observation() const {
return *std::max_element(this->_observations.begin(), this->_observations.end());
return *std::max_element(this->_observations.begin(), this->_observations.end(),compareVectorXd);
}
const Eigen::VectorXd& best_sample() const {
auto max_e = std::max_element(this->_observations.begin(), this->_observations.end());
auto max_e = std::max_element(this->_observations.begin(), this->_observations.end(),compareVectorXd);
return this->_samples[std::distance(this->_observations.begin(), max_e)];
}
......
#ifndef GP_HPP_
#define GP_HPP_
......@@ -6,22 +7,24 @@
#include <Eigen/LU>
#include <Eigen/Cholesky>
#include "kernel_functions.hpp"
#include "mean_functions.hpp"
namespace limbo {
namespace model {
template<typename Params, typename KernelFunction, typename MeanFunction>
template<typename Params, typename KernelFunction, typename MeanFunction, typename ObsType= Eigen::VectorXd>
class GP {
public:
GP() : _dim(-1) {}
// useful because the model might created before having samples
// useful because the model might be created before having samples
GP(int d) : _dim(d), _kernel_function(d) {}
void compute(const std::vector<Eigen::VectorXd>& samples,
const std::vector<double>& observations,
const std::vector<ObsType>& observations,
double noise) {
if (_dim == -1) {
assert(samples.size() != 0);
assert(observations.size() != 0);
......@@ -29,23 +32,29 @@ namespace limbo {
_dim = samples[0].size();
}
_samples = samples;
_observations.resize(observations.size());
_observations.resize(observations.size(),observations[0].size());
_noise = noise;
for (int i = 0; i < _observations.size(); ++i)
_observations(i) = observations[i];
_mean_observation = _observations.sum() / _observations.size();
int obs_dim=observations[0].size();
_mean_vector.resize(_samples.size());
for (int i = 0; i < _mean_vector.size(); i++)
_mean_vector(i) = _mean_function(_samples[i], *this);
for (int i = 0; i < _observations.rows(); ++i)
_observations.row(i) = observations[i];
_mean_observation.resize(obs_dim);
for(int i=0; i< _observations.cols(); i++)
_mean_observation(i) = _observations.col(i).sum() / _observations.rows();
_mean_vector.resize(_samples.size(),obs_dim);
for (int i = 0; i < _mean_vector.rows(); i++)
_mean_vector.row(i) = ObsType::Zero(obs_dim)+_mean_function(_samples[i], *this); // small trick to accept either Double or Vector
_obs_mean = _observations - _mean_vector;
_compute_kernel();
}
// return mu, sigma (unormaliz)
std::tuple<double, double> query(const Eigen::VectorXd& v) const {
std::tuple<ObsType, double> query(const Eigen::VectorXd& v) const {
if (_samples.size() == 0)
return std::make_tuple(_mean_function(v, *this),
sqrt(_kernel_function(v, v)));
......@@ -54,7 +63,7 @@ namespace limbo {
return std::make_tuple(_mu(v, k), _sigma(v, k));
}
double mu(const Eigen::VectorXd& v) const {
ObsType mu(const Eigen::VectorXd& v) const {
if (_samples.size() == 0)
return _mean_function(v, *this);
return _mu(v, _compute_k(v));
......@@ -74,10 +83,12 @@ namespace limbo {
const MeanFunction& mean_function() const {
return _mean_function;
}
double max_observation() const {
ObsType max_observation() const {
if(_observations.cols()>1)
std::cout<<"WARNING max_observation with multi dimensional observations doesn't make sense"<<std::endl;
return _observations.maxCoeff();
}
double mean_observation() const {
ObsType mean_observation() const {
return _mean_observation;
}
protected:
......@@ -86,13 +97,13 @@ namespace limbo {
MeanFunction _mean_function;
std::vector<Eigen::VectorXd> _samples;
Eigen::VectorXd _observations;
Eigen::VectorXd _mean_vector;
Eigen::VectorXd _obs_mean;
Eigen::MatrixXd _observations;
Eigen::MatrixXd _mean_vector;
Eigen::MatrixXd _obs_mean;
double _noise;
Eigen::VectorXd _alpha;
double _mean_observation;
Eigen::MatrixXd _alpha;
ObsType _mean_observation;
Eigen::MatrixXd _kernel;
Eigen::MatrixXd _inverted_kernel;
......@@ -101,9 +112,9 @@ namespace limbo {
void _compute_kernel() {
// O(n^2) [should be negligible]
_kernel.resize(_observations.size(), _observations.size());
for (int i = 0; i < _observations.size(); i++)
for (int j = 0; j < _observations.size(); ++j)
_kernel.resize(_samples.size(), _samples.size());
for (size_t i = 0; i < _samples.size(); i++)
for (size_t j = 0; j < _samples.size(); ++j)
_kernel(i, j) = _kernel_function(_samples[i], _samples[j]) + _noise;
// O(n^3)
......@@ -114,10 +125,11 @@ namespace limbo {
// alpha = K^{-1} * this->_obs_mean;
_alpha = _llt.matrixL().solve(this->_obs_mean);
_llt.matrixL().adjoint().solveInPlace(_alpha);
}
double _mu(const Eigen::VectorXd& v, const Eigen::VectorXd& k) const {
return _mean_function(v, *this) + k.transpose() * _alpha;
ObsType _mu(const Eigen::VectorXd& v, const Eigen::VectorXd& k) const {
return (k.transpose() * _alpha) + _mean_function(v, *this);
// return _mean_function(v)
// + (k.transpose() * _inverted_kernel * (_obs_mean))[0];
}
......
......@@ -22,14 +22,14 @@ namespace limbo {
};
}
namespace model {
template<typename Params, typename KernelFunction, typename MeanFunction>
template<typename Params, typename KernelFunction, typename MeanFunction, typename ObsType=Eigen::VectorXd>
class GPAuto : public GP<Params, KernelFunction, MeanFunction> {
public:
// TODO : init KernelFunction with dim in GP
GPAuto(int d) : GP<Params, KernelFunction, MeanFunction>(d) {}
void compute(const std::vector<Eigen::VectorXd>& samples,
const std::vector<double>& observations,
const std::vector<ObsType>& observations,
double noise) {
GP<Params, KernelFunction, MeanFunction>::compute(samples, observations, noise);
......@@ -43,7 +43,7 @@ namespace limbo {
this->_kernel_function.set_h_params(h_params);
if (update_kernel)
this->_compute_kernel();
size_t n = this->_obs_mean.size();
size_t n = this->_obs_mean.rows();
// --- cholesky ---
// see: http://xcorr.net/2008/06/11/log-determinant-of-positive-definite-matrices-in-matlab/
......@@ -52,7 +52,6 @@ namespace limbo {
// alpha = K^{-1} * this->_obs_mean;
double a = this->_obs_mean.dot(this->_alpha);
return -0.5 * a - 0.5 * det - 0.5 * n * log(2 * M_PI);
}
......
......@@ -30,13 +30,13 @@ namespace limbo {
Cmaes() {}
template <typename AcquisitionFunction>
Eigen::VectorXd operator()(const AcquisitionFunction& acqui, int dim) const {
return this->operator()(acqui, dim, Eigen::VectorXd::Constant(dim, 0.5));
return this->operator()(acqui, dim, Eigen::VectorXd::Constant(dim, 0.5));
}
template <typename AcquisitionFunction>
Eigen::VectorXd operator()(const AcquisitionFunction& acqui, int dim, const Eigen::VectorXd& init) const {
int nrestarts = Params::cmaes::nrestarts();
double incpopsize = 2;
cmaes_t evo;
cmaes_t evo;
double *const*pop;
double *fitvals;
double fbestever = 0, *xbestever = NULL;
......@@ -53,8 +53,11 @@ namespace limbo {
double init_point[dim];
for (int i = 0; i < dim; ++i)
init_point[i] = init(i);
for (irun = 0; irun < nrestarts + 1; ++irun) {
fitvals = cmaes_init(&evo, acqui.dim(), init_point, NULL, 0, lambda, NULL);
evo.countevals = countevals;
evo.sp.stopMaxFunEvals =
Params::cmaes::max_fun_evals() < 0 ?
......@@ -70,13 +73,14 @@ namespace limbo {
while (!(stop = cmaes_TestForTermination(&evo))) {
pop = cmaes_SamplePopulation(&evo);
par::loop(0, pop_size, [&](int i) {
boundary_transformation(&boundaries, pop[i], all_x_in_bounds[i], dim);
for (int j = 0; j < dim; ++j)
pop_eigen[i](j) = x_in_bounds[j];
fitvals[i] = -acqui(pop_eigen[i]);
});
boundary_transformation(&boundaries, pop[i], all_x_in_bounds[i], dim);
for (int j = 0; j < dim; ++j)
pop_eigen[i](j) = x_in_bounds[j];
fitvals[i] = -acqui(pop_eigen[i]);
});
cmaes_UpdateDistribution(&evo, fitvals);
}
for (int i = 0; i < pop_size; ++i)
free(all_x_in_bounds[i]);
......@@ -97,7 +101,7 @@ namespace limbo {
xbestever = cmaes_GetInto(&evo, "xmean", xbestever);
}
cmaes_exit(&evo);
cmaes_exit(&evo);
if (stop) {
if (strncmp(stop, "Fitness", 7) == 0 || strncmp(stop, "MaxFunEvals", 11) == 0) {
......
......@@ -32,13 +32,13 @@ namespace limbo {
}
};
template<typename Params>
template<typename Params, typename ObsType=Eigen::VectorXd>
struct MeanData {
MeanData() {
}
template<typename GP>
double operator()(const Eigen::VectorXd& v, const GP& gp)const {
return gp.mean_observation();
ObsType operator()(const Eigen::VectorXd& v, const GP& gp)const {
return gp.mean_observation().array();
}
};
......
......@@ -18,6 +18,7 @@ def options(opt):
opt.load('tbb')
opt.load('sferes')
opt.load('limbo')
opt.load('ode')
opt.add_option('--exp', type='string', help='exp(s) to build, separate by comma', dest='exp')
for i in glob.glob('exp/*'):
opt.recurse(i)
......@@ -42,14 +43,18 @@ def configure(conf):
conf.check_eigen()
conf.check_tbb()
conf.check_sferes()
conf.check_ode()
if conf.is_defined('USE_TBB'):
common_flags += " -DUSE_TBB"
if conf.is_defined('USE_SFERES'):
common_flags += " -DUSE_SFERES -DSFERES_FAST_DOMSORT"
if conf.env['CXXFLAGS_ODE']:
common_flags += ' '+ conf.env['CXXFLAGS_ODE']
opt_flags = common_flags + ' -O3 -msse2 -g'
conf.env['CXXFLAGS'] = cxxflags + opt_flags.split(' ')
conf.env['CXXFLAGS']= cxxflags + opt_flags.split(' ')
print conf.env['CXXFLAGS']
if conf.options.exp:
......
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