MeLOn
matern.h
Go to the documentation of this file.
1 /**********************************************************************************
2 * Copyright (c) 2020 Process Systems Engineering (AVT.SVT), RWTH Aachen University
3 *
4 * This program and the accompanying materials are made available under the
5 * terms of the Eclipse Public License 2.0 which is available at
6 * http://www.eclipse.org/legal/epl-2.0.
7 *
8 * SPDX-License-Identifier: EPL-2.0
9 *
10 * @file gp_data.h
11 *
12 * @brief File containing declaration of matern kernel classes
13 *
14 **********************************************************************************/
15 
16 #pragma once
17 
18 #include "kernel.h"
19 
20 #define USE_TAYLORMADE_RELAXATIONS
21 //#undef USE_TAYLORMADE_RELAXATIONS
22 
23 #ifdef USE_TAYLORMADE_RELAXATIONS
24 #include "ffunc.hpp"
25 #endif
26 
27 
28 namespace melon {
29  namespace kernel {
30 
31 
32 #ifdef USE_TAYLORMADE_RELAXATIONS
33 
34  using mc::covariance_function;
35 
36  template <typename T> inline T
37  covar_matern_1(const T& x) {
38  return covariance_function(x, 1);
39  }
40 
41  template <typename T> inline T
43  (const T& x) {
44  return covariance_function(x, 2);
45  }
46 
47  template <typename T> inline T
49  (const T& x) {
50  return covariance_function(x, 3);
51  }
52 
53  template <typename T> inline T
55  (const T& x) {
56  return covariance_function(x, 4);
57  }
58 
59 #endif
60 
64  struct KernelData {
65  double sf2;
66  std::vector<double> ell;
67  };
68 
73  template <typename T, typename V>
74  class Matern : public StationaryKernel<T, V> {
75  public:
76  Matern();
77  Matern(std::shared_ptr<const KernelData> data);
78 
79  protected:
80  using typename Kernel<T, V>::RET;
81 
82  std::shared_ptr<const KernelData> _data;
91  RET _quadratic_distance(std::vector<T> x1, std::vector<V> x2) override {
92  RET distance = 0;
93 
94  for (size_t i = 0; i < x1.size(); i++) { // i the demension of X and X_test
95  distance += pow(x1.at(i) - x2.at(i), 2) / pow(Matern<T, V>::_data->ell.at(i), 2);
96  }
97 
98  return distance;
99  };
100 
106  auto _euclidian_distance(RET quadraticDistance) {
107 
108  double const epsilon = 1e-16;
109  // Small number to avoid gradient calculation of sqrt(x) function at x = 0.
110  // Otherwise this leads to longer computation time of IPOPT iterations in the preprocessing. (SLSQP works very quick with epsilon = 0).
111  // In an example, epsilon = 1e-16 leads to an error of the prediction in the order of 1e-8.
112  // We accept this error.
113 
114  auto distance = quadraticDistance + epsilon;
115  distance = sqrt(distance);
116 
117  return distance;
118  };
119  };
120 
125  template <typename T, typename V>
126  class Matern12 : public Matern<T, V> {
127  public:
128  using Matern<T, V>::Matern;
129  using typename Kernel<T, V>::RET;
130 
131  RET evaluate_kernel(std::vector<T> x1, std::vector<V> x2);
132  RET evaluate_kernel(RET distance);
133  RET calculate_distance(std::vector<T> x1, std::vector<V> x2);
134  };
135 
140  template <typename T, typename V>
141  class Matern32 : public Matern<T, V> {
142  public:
143  using Matern<T, V>::Matern;
144  using typename Kernel<T, V>::RET;
145 
146  RET evaluate_kernel(std::vector<T> x1, std::vector<V> x2);
147  RET evaluate_kernel(RET distance);
148  RET calculate_distance(std::vector<T> x1, std::vector<V> x2);
149  };
150 
155  template <typename T, typename V>
156  class Matern52 : public Matern<T, V> {
157  public:
158  using Matern<T, V>::Matern;
159  using typename Kernel<T, V>::RET;
160 
161  RET evaluate_kernel(std::vector<T> x1, std::vector<V> x2);
162  RET evaluate_kernel(RET distance);
163  RET calculate_distance(std::vector<T> x1, std::vector<V> x2);
164  };
165 
170  template <typename T, typename V>
171  class MaternInf : public Matern<T, V> {
172  public:
173  using Matern<T, V>::Matern;
174  using typename Kernel<T, V>::RET;
175 
176  RET evaluate_kernel(std::vector<T> x1, std::vector<V> x2);
177  RET evaluate_kernel(RET distance);
178  RET calculate_distance(std::vector<T> x1, std::vector<V> x2);
179  };
180 
185  template<typename T, typename V>
187 
188  template<typename T, typename V>
189  Matern<T, V>::Matern(std::shared_ptr<const KernelData> data) : Matern() {
190  Matern<T, V>::_data = data;
191  }
192 
193 
194  // ---------------------------------------------------------------------------------
195  // Matern1/2
196  // ---------------------------------------------------------------------------------
197 
198  template<typename T, typename V>
200 
201 #ifdef USE_TAYLORMADE_RELAXATIONS
202  return Matern<T, V>::_data->sf2 * covar_matern_1(distance);
203 #else
204  RET euclidianDistance = sqrt(distance + 1e-16);
205  return Matern<T, V>::_data->sf2 * exp(-euclidianDistance);
206 #endif
207  }
208 
209  template<typename T, typename V>
210  typename Matern12<T, V>::RET Matern12<T, V>::evaluate_kernel(std::vector<T> x1, std::vector<V> x2) {
211  RET distances = calculate_distance(x1, x2);
212  return evaluate_kernel(distances);
213  }
214 
215  template<typename T, typename V>
216  typename Matern12<T, V>::RET Matern12<T, V>::calculate_distance(std::vector<T> x1, std::vector<V> x2) {
217  return this->_quadratic_distance(x1, x2);
218  }
219 
220 
221  // ---------------------------------------------------------------------------------
222  // Matern3/2
223  // ---------------------------------------------------------------------------------
224 
225  template<typename T, typename V>
227 
228 #ifdef USE_TAYLORMADE_RELAXATIONS
229  return Matern<T, V>::_data->sf2 * covar_matern_3(distance);
230 #else
231  RET euclidianDistance = sqrt(distance + 1e-16);
232  RET sqrtK = sqrt(3.0) * euclidianDistance;
233  // we multiply the covariance function out because we found that it has tighter relaxations
234  return Matern<T, V>::_data->sf2 * (exp(-sqrtK) + sqrtK * exp(-sqrtK));
235 #endif
236 
237  }
238 
239 
240  template<typename T, typename V>
241  typename Matern32<T, V>::RET Matern32<T, V>::evaluate_kernel(std::vector<T> x1, std::vector<V> x2) {
242  RET distances = calculate_distance(x1, x2);
243  return evaluate_kernel(distances);
244  }
245 
246  template<typename T, typename V>
247  typename Matern32<T, V>::RET Matern32<T, V>::calculate_distance(std::vector<T> x1, std::vector<V> x2) {
248  return this->_quadratic_distance(x1, x2);
249  }
250 
251 
252  // ---------------------------------------------------------------------------------
253  // Matern5/2
254  // ---------------------------------------------------------------------------------
255 
256  template<typename T, typename V>
258 
259 #ifdef USE_TAYLORMADE_RELAXATIONS
260  // use xexpax as envelope is implemented. return _data->sf2 * (1 + sqrtK) * exp(-sqrtK);
261  return Matern<T, V>::_data->sf2 * covar_matern_5(distance);
262 #else
263  RET euclidianDistance = sqrt(distance + 1e-16);
264  RET sqrtK = sqrt(5.0) * euclidianDistance;
265  RET sqrtK2_3 = 5. / 3. * distance;
266  // we multiply the covariance function out because we found that it has tighter relaxations
267  // return Matern<T, V>::_data->sf2 * (1 + sqrtK + sqrtK2_3) * exp(-sqrtK);
268  return Matern<T, V>::_data->sf2 * (exp(-sqrtK) + sqrtK * exp(-sqrtK) + sqrtK2_3 * exp(-sqrtK));
269 #endif
270 
271  }
272 
273  template<typename T, typename V>
274  typename Matern52<T, V>::RET Matern52<T, V>::evaluate_kernel(std::vector<T> x1, std::vector<V> x2) {
275  RET distances = calculate_distance(x1, x2);
276  return evaluate_kernel(distances);
277  }
278 
279  template<typename T, typename V>
280  typename Matern52<T, V>::RET Matern52<T, V>::calculate_distance(std::vector<T> x1, std::vector<V> x2) {
281  return this->_quadratic_distance(x1, x2);
282  }
283 
284 
285  // ---------------------------------------------------------------------------------
286  // MaternInf (squared exponential)
287  // ---------------------------------------------------------------------------------
288 
289  template<typename T, typename V>
291 
292 #ifdef USE_TAYLORMADE_RELAXATIONS
293  return Matern<T, V>::_data->sf2 * covar_sqrexp(distance);
294 #else
295  return Matern<T, V>::_data->sf2 * exp(-0.5 * distance);
296 #endif
297 
298  }
299 
300  template<typename T, typename V>
301  typename MaternInf<T, V>::RET MaternInf<T, V>::evaluate_kernel(std::vector<T> x1, std::vector<V> x2) {
302  RET distances = calculate_distance(x1, x2);
303  return evaluate_kernel(distances);
304  }
305 
306 
307  template<typename T, typename V>
308  typename MaternInf<T, V>::RET MaternInf<T, V>::calculate_distance(std::vector<T> x1, std::vector<V> x2) {
309  return this->_quadratic_distance(x1, x2);
310  }
311  }
312 }
kernel.h
melon::kernel::Matern::Matern
Matern()
Definition: matern.h:186
melon::kernel::covar_sqrexp
T covar_sqrexp(const T &x)
Definition: matern.h:55
melon::kernel::MaternInf::calculate_distance
RET calculate_distance(std::vector< T > x1, std::vector< V > x2)
Function for calculating the distance used in the kernel (type of distance used can vary among kernel...
Definition: matern.h:308
melon::kernel::Matern::_euclidian_distance
auto _euclidian_distance(RET quadraticDistance)
Calculates the euclidian distance from the quadratic distance.
Definition: matern.h:106
melon::kernel::Matern52
Class implementation of Matern52 kernel.
Definition: matern.h:156
melon::kernel::Matern12::evaluate_kernel
RET evaluate_kernel(std::vector< T > x1, std::vector< V > x2)
Function for evalualting the kernel for the points x1 and x2.
Definition: matern.h:210
melon::kernel::Matern52::calculate_distance
RET calculate_distance(std::vector< T > x1, std::vector< V > x2)
Function for calculating the distance used in the kernel (type of distance used can vary among kernel...
Definition: matern.h:280
melon::kernel::KernelData
struct containing kernel parameters
Definition: matern.h:64
melon::kernel::KernelData::sf2
double sf2
Definition: matern.h:65
melon::kernel::StationaryKernel
Definition: kernel.h:84
melon::kernel::MaternInf
Class implementation of MaternInf kernel.
Definition: matern.h:171
melon::kernel::Matern12::calculate_distance
RET calculate_distance(std::vector< T > x1, std::vector< V > x2)
Function for calculating the distance used in the kernel (type of distance used can vary among kernel...
Definition: matern.h:216
melon::kernel::covar_matern_5
T covar_matern_5(const T &x)
Definition: matern.h:49
melon::kernel::Matern52::evaluate_kernel
RET evaluate_kernel(std::vector< T > x1, std::vector< V > x2)
Function for evalualting the kernel for the points x1 and x2.
Definition: matern.h:274
melon::kernel::KernelData::ell
std::vector< double > ell
Definition: matern.h:66
melon::kernel::Matern12
Class implementation of Matern12 kernel.
Definition: matern.h:126
melon::kernel::Kernel::RET
decltype(std::declval< T >()+std::declval< V >()) RET
Definition: kernel.h:59
melon::kernel::Matern::_quadratic_distance
RET _quadratic_distance(std::vector< T > x1, std::vector< V > x2) override
Calculates the quadratic distance between two points x1 and x2.
Definition: matern.h:91
melon::kernel::Matern32
Class implementation of Matern32 kernel.
Definition: matern.h:141
melon::kernel::Matern
Abstract parent class for matern kernels.
Definition: matern.h:74
melon::kernel::covar_matern_1
T covar_matern_1(const T &x)
Definition: matern.h:37
melon::kernel::covar_matern_3
T covar_matern_3(const T &x)
Definition: matern.h:43
melon::kernel::Matern32::evaluate_kernel
RET evaluate_kernel(std::vector< T > x1, std::vector< V > x2)
Function for evalualting the kernel for the points x1 and x2.
Definition: matern.h:241
melon::kernel::MaternInf::evaluate_kernel
RET evaluate_kernel(std::vector< T > x1, std::vector< V > x2)
Function for evalualting the kernel for the points x1 and x2.
Definition: matern.h:301
melon
Definition: kernel.h:21
melon::kernel::Matern32::calculate_distance
RET calculate_distance(std::vector< T > x1, std::vector< V > x2)
Function for calculating the distance used in the kernel (type of distance used can vary among kernel...
Definition: matern.h:247
melon::kernel::Matern::_data
std::shared_ptr< const KernelData > _data
Definition: matern.h:82