MeLOn
Loading...
Searching...
No Matches
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
28namespace 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
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}
decltype(std::declval< T >()+std::declval< V >()) RET
Definition kernel.h:33
Class implementation of Matern12 kernel.
Definition matern.h:126
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
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
Class implementation of Matern32 kernel.
Definition matern.h:141
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
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
Class implementation of Matern52 kernel.
Definition matern.h:156
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
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
Abstract parent class for matern kernels.
Definition matern.h:74
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
auto _euclidian_distance(RET quadraticDistance)
Calculates the euclidian distance from the quadratic distance.
Definition matern.h:106
Matern()
Definition matern.h:186
std::shared_ptr< const KernelData > _data
Definition matern.h:82
Class implementation of MaternInf kernel.
Definition matern.h:171
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
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
Definition kernel.h:58
T covar_matern_5(const T &x)
Definition matern.h:49
T covar_sqrexp(const T &x)
Definition matern.h:55
T covar_matern_1(const T &x)
Definition matern.h:37
T covar_matern_3(const T &x)
Definition matern.h:43
Definition kernel.h:21
struct containing kernel parameters
Definition matern.h:64
double sf2
Definition matern.h:65
std::vector< double > ell
Definition matern.h:66