Commit cafa4dc8 authored by Aneoshun's avatar Aneoshun
Browse files

GPs now deal with multi dimensional observations. Likelihood also works with it

parent 696e1c40
......@@ -243,7 +243,7 @@ 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!=NULL && 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");
}*/
......
......@@ -33,27 +33,38 @@ namespace limbo {
this->_init(feval, reset);
model_t model(EvalFunction::dim);
model.compute(this->_samples, this->_observations, Params::boptimizer::noise());
if(this->_observations.size())
model.compute(this->_samples, this->_observations, Params::boptimizer::noise());
inner_optimization_t inner_optimization;
while (this->_samples.size() == 0 || this->_pursue(*this)) {
acquisition_function_t acqui(model, this->_iteration);
Eigen::VectorXd new_sample = inner_optimization(acqui, acqui.dim());
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());
this->_update_stats(*this);
std::cout << this->_iteration << " new point: "
<< this->_samples[this->_samples.size() - 1].transpose()
<< " value: " << this->_observations[this->_observations.size() - 1].transpose()
// << " mu: "<< model.mu(this->_samples[this->_samples.size() - 1]).transpose()
// << " sigma: "<< model.sigma(this->_samples[this->_samples.size() - 1])
// << " acqui: "<< acqui(this->_samples[this->_samples.size() - 1])
<< " best:" << this->best_observation().transpose()
<< std::endl;
model.compute(this->_samples, this->_observations, Params::boptimizer::noise());
this->_update_stats(*this);
this->_iteration++;
}
}
......
......@@ -6,6 +6,7 @@
#include <Eigen/Core>
#include <Eigen/LU>
#include <Eigen/Cholesky>
#include <limits>
#include "kernel_functions.hpp"
#include "mean_functions.hpp"
......@@ -25,12 +26,14 @@ namespace limbo {
const std::vector<ObsType>& observations,
double noise) {
if (_dim == -1) {
assert(samples.size() != 0);
assert(observations.size() != 0);
assert(samples.size() == observations.size());
_dim = samples[0].size();
}
_samples = samples;
_observations.resize(observations.size(),observations[0].size());
_noise = noise;
......@@ -45,7 +48,7 @@ namespace limbo {
_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
_mean_vector.row(i) = _mean_function(_samples[i], *this);
_obs_mean = _observations - _mean_vector;
......@@ -91,6 +94,8 @@ namespace limbo {
ObsType mean_observation() const {
return _mean_observation;
}
const Eigen::MatrixXd& mean_vector() const{return _mean_vector;}
protected:
int _dim;
KernelFunction _kernel_function;
......@@ -111,11 +116,13 @@ namespace limbo {
Eigen::LLT<Eigen::MatrixXd> _llt;
void _compute_kernel() {
// O(n^2) [should be negligible]
_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;
_kernel(i, j) = _kernel_function(_samples[i], _samples[j]) + ((i==j)?_noise:0);
// O(n^3)
// _inverted_kernel = _kernel.inverse();
......@@ -129,16 +136,20 @@ namespace limbo {
}
ObsType _mu(const Eigen::VectorXd& v, const Eigen::VectorXd& k) const {
return (k.transpose() * _alpha) + _mean_function(v, *this);
return (k.transpose() * _alpha) + _mean_function(v, *this).transpose();
// return _mean_function(v)
// + (k.transpose() * _inverted_kernel * (_obs_mean))[0];
}
double _sigma(const Eigen::VectorXd& v, const Eigen::VectorXd& k) const {
Eigen::VectorXd z = _llt.matrixL().solve(k);
return _kernel_function(v, v) - z.dot(z);
double res= _kernel_function(v, v) - z.dot(z);
return (res<=std::numeric_limits<double>::epsilon())?0:res;
// return _kernel_function(v, v) - (k.transpose() * _inverted_kernel * k)[0];
}
Eigen::VectorXd _compute_k(const Eigen::VectorXd& v) const {
Eigen::VectorXd k(_samples.size());
for (int i = 0; i < k.size(); i++)
k[i] = _kernel_function(_samples[i], v);
......
......@@ -33,10 +33,16 @@ namespace limbo {
double noise) {
GP<Params, KernelFunction, MeanFunction>::compute(samples, observations, noise);
_optimize_likelihood();
_optimize_likelihood();
this->_compute_kernel();
}
Eigen::VectorXd check_inverse(){
return this->_kernel*this->_alpha.col(0) - this->_obs_mean;
}
// see Rasmussen and Williams, 2006 (p. 113)
double log_likelihood(const Eigen::VectorXd& h_params,
bool update_kernel = true) {
......@@ -51,8 +57,12 @@ namespace limbo {
long double det = 2 * l.diagonal().array().log().sum();
// 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);
// double a = this->_obs_mean.col(0).dot(this->_alpha.col(0));
double a = (this->_obs_mean.transpose() * this->_alpha).trace(); // generalization for multi dimensional observation
// std::cout<<" a: "<<a <<" det: "<< det<<std::endl;
double lik= -0.5 * a - 0.5 * det - 0.5 * n * log(2 * M_PI);
return lik;
}
// see Rasmussen and Williams, 2006 (p. 114)
......@@ -61,14 +71,14 @@ namespace limbo {
this->_kernel_function.set_h_params(h_params);
if (update_kernel)
this->_compute_kernel();
size_t n = this->_observations.size();
size_t n = this->_observations.rows();
/// what we should write, but it is less numerically stable than using the Cholesky decomposition
// Eigen::MatrixXd alpha = this->_inverted_kernel * this->_obs_mean;
// Eigen::MatrixXd w = alpha * alpha.transpose() - this->_inverted_kernel;
// alpha = K^{-1} * this->_obs_mean;
Eigen::VectorXd alpha = this->_llt.solve(this->_obs_mean);
Eigen::MatrixXd alpha = this->_llt.matrixL().solve(this->_obs_mean);
this->_llt.matrixL().adjoint().solveInPlace(alpha);
// K^{-1} using Cholesky decomposition
......@@ -107,7 +117,8 @@ namespace limbo {
return gp.log_likelihood_grad(v, false);
},
this->kernel_function().h_params_size(), Params::gp_auto::n_rprop());
double lik = this->log_likelihood(v);
double lik = gp.log_likelihood(v);//this->log_likelihood(v);
return std::make_pair(v, lik);
};
auto comp = [](const pair_t& v1, const pair_t& v2) {
......
......@@ -28,7 +28,7 @@ namespace limbo {
void operator()(const F& feval, Opt& opt) const {
for (int i = 0; i < Params::init::nb_samples(); i++) {
Eigen::VectorXd new_sample(F::dim);
for (size_t i = 0; i < F::dim; i++)
for (int i = 0; i < F::dim; i++)
new_sample[i] = misc::rand<double>(0, 1);
std::cout << "random sample:" << new_sample.transpose() << std::endl;
opt.add_new_sample(new_sample, feval(new_sample));
......
......@@ -19,7 +19,7 @@ namespace limbo {
namespace defaults {
struct cmaes {
BO_PARAM(int, nrestarts, 10);
BO_PARAM(int, nrestarts, 1);
BO_PARAM(double, max_fun_evals, -1);
};
}
......@@ -34,6 +34,8 @@ namespace limbo {
}
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;
......@@ -72,10 +74,10 @@ namespace limbo {
while (!(stop = cmaes_TestForTermination(&evo))) {
pop = cmaes_SamplePopulation(&evo);
par::loop(0, pop_size, [&](int i) {
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];
pop_eigen[i](j) = all_x_in_bounds[i][j];
fitvals[i] = -acqui(pop_eigen[i]);
});
cmaes_UpdateDistribution(&evo, fitvals);
......
......@@ -56,7 +56,7 @@ namespace limbo {
template<typename Params>
struct SquaredExpARD {
SquaredExpARD(int dim) : _sf2(0), _ell(dim), _input_dim(dim){
this->set_h_params(Eigen::VectorXd::Ones(_ell.size()+1));
this->set_h_params(Eigen::VectorXd::Ones(_ell.size()+1)*-1);
}
size_t h_params_size() const {
return _ell.size() + 1;
......@@ -68,7 +68,7 @@ namespace limbo {
_h_params = p;
for (size_t i = 0; i < _input_dim; ++i)
_ell(i) = exp(p(i));
_sf2 = exp(2 * p(_input_dim));
_sf2 = 1;//exp(2 * p(_input_dim));
}
Eigen::VectorXd grad(const Eigen::VectorXd &x1, const Eigen::VectorXd &x2) const {
......@@ -76,7 +76,7 @@ namespace limbo {
Eigen::VectorXd z = (x1 - x2).cwiseQuotient(_ell).array().square();
double k = _sf2 * exp(-0.5 * z.sum());
grad.head(_input_dim) = z * k;
grad(_input_dim) = 2.0 * k;
grad(_input_dim) = 0;//2.0 * k;
return grad;
}
......
......@@ -13,21 +13,21 @@ namespace limbo {
};
}
namespace mean_functions {
template<typename Params>
template<typename Params,typename ObsType=Eigen::VectorXd>
struct NullFunction {
NullFunction() {}
template<typename GP>
double operator()(const Eigen::VectorXd& v, const GP&)const {
return 0;
ObsType operator()(const Eigen::VectorXd& v, const GP&)const {
return ObsType::Zeros(Params::NullFunction::dim());
}
};
template<typename Params>
template<typename Params,typename ObsType=Eigen::VectorXd>
struct MeanConstant {
MeanConstant() {
}
template<typename GP>
double operator()(const Eigen::VectorXd& v, const GP&)const {
ObsType operator()(const Eigen::VectorXd& v, const GP&)const {
return Params::meanconstant::constant();
}
};
......@@ -44,7 +44,7 @@ namespace limbo {
template<typename Params>
template<typename Params,typename ObsType=Eigen::VectorXd>
struct MeanArchive {
MeanArchive() {
// create and open an archive for input
......@@ -57,7 +57,7 @@ namespace limbo {
}
template<typename GP>
double operator()(const Eigen::VectorXd& v, const GP&)const {
ObsType operator()(const Eigen::VectorXd& v, const GP&)const {
std::vector<double> key(v.size(), 0);
for (int i = 0; i < v.size(); i++)
key[i] = v[i];
......
......@@ -25,7 +25,7 @@ namespace rprop {
Eigen::VectorXd delta = Eigen::VectorXd::Ones(param_dim) * delta0;
Eigen::VectorXd grad_old = Eigen::VectorXd::Zero(param_dim);
Eigen::VectorXd params = (Eigen::VectorXd::Random(param_dim).array() + 1);
Eigen::VectorXd params = (Eigen::VectorXd::Random(param_dim).array() -1);
Eigen::VectorXd best_params = params;
double best = log(0);
......
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