benchmark_template.cpp 8.61 KB
 1 2 3 4 5 ``````#include #include #include #include #include `````` Konstantinos Chatzilygeroudis committed Jun 12, 2017 6 7 ``````#include #include `````` 8 9 ``````#include "test_functions.hpp" `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 10 11 12 13 14 15 ``````#ifdef USE_LIBGP #include #include #include #endif `````` 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 ``````using namespace limbo; template inline T gaussian_rand(T m = 0.0, T v = 1.0) { static std::random_device rd; static std::mt19937 gen(rd()); std::normal_distribution gaussian(m, v); return gaussian(gen); } @PARAMS inline Eigen::MatrixXd sample_covariance(const std::vector& points) { assert(points.size()); // Get the sample means Eigen::VectorXd means = Eigen::VectorXd::Zero(points[0].size()); for (size_t i = 0; i < points.size(); i++) { means.array() += points[i].array(); } means = means.array() / double(points.size()); // Calculate the sample covariance matrix Eigen::MatrixXd cov = Eigen::MatrixXd::Zero(points[0].size(), points[0].size()); for (size_t i = 0; i < points.size(); i++) { cov = cov + points[i] * points[i].transpose(); } cov = (cov.array() - (double(points.size()) * means * means.transpose()).array()) / (double(points.size()) - 1.0); return cov; } template void benchmark(const std::string& name, std::vector dimensions, std::vector points) { Function func; std::vector dims; if (func.dims() > 0) dims.push_back(func.dims()); else { dims = dimensions; } // Whether to add noise or not bool add_noise = @NOISE; int N_test = 10000; std::ofstream ofs_res(name + ".dat"); `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 70 71 72 ``````#ifdef USE_LIBGP std::ofstream ofs_libgp(name + "_libgp.dat"); #endif `````` 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 `````` for (size_t dim = 0; dim < dims.size(); dim++) { std::vector bounds = func.bounds(); bool one_bound = (bounds.size() == 1); int D = dims[dim]; for (size_t n = 0; n < points.size(); n++) { int N = points[n]; if (N <= 0) { std::cerr << "Number of points less or equal to zero!" << std::endl; continue; } // Output number of models ofs_res << D << " " << N << " @NMODELS" << std::endl; `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 88 89 90 ``````#ifdef USE_LIBGP ofs_libgp << D << " " << N << " 1" << std::endl; #endif `````` 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 `````` std::cout << name << " in dim: " << D << " and # of points: " << N << std::endl; std::string file_name = name + "_" + std::to_string(D) + "_" + std::to_string(N); std::vector points, obs; for (int i = 0; i < N; i++) { Eigen::VectorXd p = limbo::tools::random_vector(D); //.array() * 10.24 - 5.12; if (one_bound) p = p.array() * (bounds[0](1) - bounds[0](0)) + bounds[0](0); else { for (int j = 0; j < D; j++) { p(j) = p(j) * (bounds[j](1) - bounds[j](0)) + bounds[j](0); } } points.push_back(p); Eigen::VectorXd ob(1); ob << func(p); obs.push_back(ob); } if (add_noise) { Eigen::MatrixXd cov = sample_covariance(obs); double sigma = std::sqrt(cov(0, 0)) / 20.0; std::cout << "Adding noise of: " << sigma << std::endl; for (int i = 0; i < N; i++) obs[i] = obs[i].array() + gaussian_rand(0.0, sigma); } `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 ``````#ifdef USE_LIBGP std::cout << "Training libGP GP..." << std::endl; // Learn libGP GP here libgp::GaussianProcess libgp_gp(D, "CovSum ( CovSEard, CovNoise)"); // set initial hyper-parameters Eigen::VectorXd libgp_params = Eigen::VectorXd::Zero(libgp_gp.covf().get_param_dim()); libgp_params(D - 1) = std::log(0.01); libgp_gp.covf().set_loghyper(libgp_params); // add observations to libGP for (int i = 0; i < N; i++) { libgp_gp.add_pattern(points[i].data(), obs[i](0)); } // optimize hyper-parameters of libGP libgp::RProp rprop_optimizer; double delta0 = 0.1; double deltamin = 1e-6; double deltamax = 50; double etaminus = 0.5; double etaplus = 1.2; `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 141 `````` rprop_optimizer.init(1e-2, delta0, deltamin, deltamax, etaminus, etaplus); `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 142 `````` auto start_libgp = std::chrono::high_resolution_clock::now(); `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 143 `````` rprop_optimizer.maximize(&libgp_gp, 50, false); `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 144 145 146 147 148 `````` auto time_libgp = std::chrono::duration_cast(std::chrono::high_resolution_clock::now() - start_libgp).count(); std::cout << std::setprecision(std::numeric_limits::digits10 + 1); std::cout << "Time of libGP in secs: " << time_libgp / double(1000000.0) << std::endl; #endif `````` 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 `````` // Learn the GPs code here @GPSLEARN // Generation of test data std::vector test_points, test_obs; for (int i = 0; i < N_test; i++) { Eigen::VectorXd p = limbo::tools::random_vector(D); if (one_bound) p = p.array() * (bounds[0](1) - bounds[0](0)) + bounds[0](0); else { for (int j = 0; j < D; j++) { p(j) = p(j) * (bounds[j](1) - bounds[j](0)) + bounds[j](0); } } test_points.push_back(p); Eigen::VectorXd ob(1); ob << func(p); test_obs.push_back(ob); } `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 ``````#ifdef USE_LIBGP // Predicition of the libGP GP std::vector predict_libgp(N_test); std::vector dummy_libgp(N_test); start_libgp = std::chrono::high_resolution_clock::now(); for (int i = 0; i < N_test; i++) { predict_libgp[i] = libgp_gp.f(test_points[i].data()); dummy_libgp[i] = libgp_gp.var(test_points[i].data()); // just here for the time comparisons } auto time2_libgp = std::chrono::duration_cast(std::chrono::high_resolution_clock::now() - start_libgp).count(); std::cout << "Time of libGP (query) in ms: " << time2_libgp * 1e-3 / double(N_test) << std::endl; double err_libgp = 0.0; for (int i = 0; i < N_test; i++) { Eigen::VectorXd p(1); p << predict_libgp[i]; err_libgp += (p - test_obs[i]).squaredNorm(); } err_libgp /= double(N_test); std::cout << "MSE(libGP): " << err_libgp << std::endl; // save results of libGP ofs_libgp << std::setprecision(std::numeric_limits::digits10 + 1); `````` Konstantinos Chatzilygeroudis committed Mar 30, 2018 194 `````` ofs_libgp << time_libgp / double(1000000.0) << " " << time2_libgp * 1e-3 / double(N_test) << " " << err_libgp << " SE-Full" << std::endl; `````` Konstantinos Chatzilygeroudis committed Mar 27, 2018 195 196 ``````#endif `````` 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 `````` // Prediction of the GPs code here @GPSQUERY std::cout << "Saving data..." << std::endl; // std::ofstream ofs(file_name + "_gp.dat"); // std::ofstream ofs_real(file_name + "_real.dat"); // int pp = 4000; // for (int i = 0; i < pp; ++i) { // Eigen::VectorXd v = limbo::tools::random_vector(D); // if (one_bound) // v = v.array() * (bounds[0](1) - bounds[0](0)) + bounds[0](0); // else { // for (int j = 0; j < D; j++) { // v(j) = v(j) * (bounds[j](1) - bounds[j](0)) + bounds[j](0); // } // } // Eigen::VectorXd mu; // double sigma; // std::tie(mu, sigma) = gp.query(v); // // ofs << v.transpose() << " " << mu[0] << " " << std::sqrt(sigma) << std::endl; // // double val = func(v); // ofs_real << v.transpose() << " " << val << " 0" << std::endl; // } std::ofstream ofs_data(file_name + "_data.dat"); for (size_t i = 0; i < points.size(); ++i) ofs_data << points[i].transpose() << " " << obs[i].transpose() << std::endl; std::ofstream ofs_test(file_name + "_test.dat"); for (size_t i = 0; i < test_points.size(); ++i) ofs_test << test_points[i].transpose() << " " << test_obs[i].transpose() << std::endl; std::cout << "Data saved...!" << std::endl; } } } int main(int argc, char* argv[]) { @FUNCS return 0; }``````