Kstars

tools/math_tools.cpp
Go to the documentation of this file.
1/*
2 SPDX-FileCopyrightText: 2014-2017 Max Planck Society.
3 All rights reserved.
4
5 SPDX-License-Identifier: BSD-3-Clause
6*/
7
8/**
9 * @file
10 * @date 2014-2017
11 * @copyright Max Planck Society
12 *
13 * @author Edgar D. Klenske <edgar.klenske@tuebingen.mpg.de>
14 * @author Stephan Wenninger <stephan.wenninger@tuebingen.mpg.de>
15 * @author Raffi Enficiaud <raffi.enficiaud@tuebingen.mpg.de>
16 *
17 * @brief Provides mathematical tools needed for the Gaussian process toolbox.
18 */
19
20#include "math_tools.h"
21#include <stdexcept>
22#include <cmath>
23#include <cstdint>
24
25namespace math_tools
26{
27
28Eigen::MatrixXd squareDistance(const Eigen::MatrixXd &a, const Eigen::MatrixXd &b)
29{
30 int aRows = a.rows(); // bRows must be identical to aRows
31 int aCols = a.cols();
32 int bCols = b.cols();
33
34 Eigen::MatrixXd am(aRows, aCols); // mean-corrected a
35 Eigen::MatrixXd bm(aRows, bCols); // mean-corrected b
36 // final result, aCols x bCols
37 Eigen::MatrixXd result(aCols, bCols);
38
39 Eigen::VectorXd mean(aRows);
40
41 /* If the two Matrices have the same address, it means the function was
42 called from the overloaded version, and thus the mean only has to be
43 computed once.
44 */
45 if (&a == &b) // Same address?
46 {
47 mean = a.rowwise().mean();
48 am = a.colwise() - mean;
49 bm = am;
50 }
51 else
52 {
53 if (aRows != b.rows())
54 {
55 throw std::runtime_error("Matrix dimension incorrect.");
56 }
57
58 mean = static_cast<double>(aCols) / (aCols + bCols) * a.rowwise().mean() +
59 static_cast<double>(bCols) / (bCols + aCols) * b.rowwise().mean();
60
61 // The mean of the two Matrices is subtracted beforehand, because the
62 // squared error is independent of the mean and this makes the squares
63 // smaller.
64 am = a.colwise() - mean;
65 bm = b.colwise() - mean;
66 }
67
68 // The square distance calculation (a - b)^2 is calculated as a^2 - 2*ab * b^2
69 // (using the binomial formula) because of numerical stability.
70
71 // fast version
72 return
73 ((am.array().square().colwise().sum().transpose().rowwise().replicate(bCols).matrix()
74 + bm.array().square().colwise().sum().colwise().replicate(aCols).matrix())
75 - 2 * (am.transpose()) * bm).array().max(0);
76
77 /* // verbose version
78 Eigen::MatrixXd a_square =
79 am.array().square().colwise().sum().transpose().rowwise().replicate(bCols);
80
81 Eigen::MatrixXd b_square =
82 bm.array().square().colwise().sum().colwise().replicate(aCols);
83
84 Eigen::MatrixXd twoab = 2 * (am.transpose()) * bm;
85
86 return ((a_square.matrix() + b_square.matrix()) - twoab).array().max(0);
87 */
88}
89
90Eigen::MatrixXd squareDistance(const Eigen::MatrixXd &a)
91{
92 return squareDistance(a, a);
93}
94
96 const size_t n,
97 const size_t m)
98{
99 Eigen::MatrixXd result = Eigen::MatrixXd(n, m);
100 result.setRandom(); // here we get uniform random values between -1 and 1
101 Eigen::MatrixXd temp = result.array() + 1; // shift to positive
102 result = temp / 2.0; // divide to obtain 0/1 interval
103 // prevent numerical problems by enforcing a particular interval
104 result = result.array().max(1e-10); // eliminate too small values
105 result = result.array().min(1.0); // eliminiate too large values
106 return result;
107}
108
109Eigen::MatrixXd box_muller(const Eigen::VectorXd &vRand)
110{
111 size_t n = vRand.rows();
112 size_t m = n / 2; // Box-Muller transforms pairs of numbers
113
114 Eigen::ArrayXd rand1 = vRand.head(m);
115 Eigen::ArrayXd rand2 = vRand.tail(m);
116
117 /* Implemented according to
118 * http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
119 */
120
121 // enforce interval to avoid numerical issues
122 rand1 = rand1.max(1e-10);
123 rand1 = rand1.min(1.0);
124
125 rand1 = -2 * rand1.log();
126 rand1 = rand1.sqrt(); // this is an amplitude
127
128 rand2 = rand2 * 2 * M_PI; // this is a random angle in the complex plane
129
130 Eigen::MatrixXd result(2 * m, 1);
131 Eigen::MatrixXd res1 = (rand1 * rand2.cos()).matrix(); // first elements
132 Eigen::MatrixXd res2 = (rand1 * rand2.sin()).matrix(); // second elements
133 result << res1, res2; // combine the pairs into one vector
134
135 return result;
136}
137
138Eigen::MatrixXd generate_normal_random_matrix(const size_t n, const size_t m)
139{
140 // if n*m is odd, we need one random number extra!
141 // therefore, we have to round up here.
142 size_t N = static_cast<size_t>(std::ceil(n * m / 2.0));
143
144 Eigen::MatrixXd result(2 * N, 1);
145 // push random samples through the Box-Muller transform
147 result.conservativeResize(n, m);
148 return result;
149}
150
151std::pair<Eigen::VectorXd, Eigen::VectorXd> compute_spectrum(Eigen::VectorXd &data, int N)
152{
153
154 int N_data = data.rows();
155
156 if (N < N_data)
157 {
158 N = N_data;
159 }
160 N = static_cast<int>(std::pow(2, std::ceil(std::log(N) / std::log(2)))); // map to nearest power of 2
161
162 Eigen::VectorXd padded_data = Eigen::VectorXd::Zero(N);
163 padded_data.head(N_data) = data;
164
165 Eigen::FFT<double> fft;
166
167 // initialize the double vector from Eigen vector. This works by initializing
168 // with two pointers: 1) the first element of the data, 2) the last element of the data
169 std::vector<double> vec_data(padded_data.data(), padded_data.data() + padded_data.rows() * padded_data.cols());
170 std::vector<std::complex<double> > vec_result;
171 fft.fwd(vec_result, vec_data); // this is the forward-FFT, from time domain to Fourier domain
172
173 // map back to Eigen vector with the address of the first element and the size
174 Eigen::VectorXcd result = Eigen::Map<Eigen::VectorXcd>(&vec_result[0], vec_result.size());
175
176 // the low_index is the lowest useful frequency, depending on the number of actual datapoints
177 int low_index = static_cast<int>(std::ceil(static_cast<double>(N) / static_cast<double>(N_data)));
178
179 // prepare amplitudes and frequencies, don't return frequencies introduced by padding
180 Eigen::VectorXd spectrum = result.segment(low_index, N / 2 - low_index + 1).array().abs().pow(2);
181 Eigen::VectorXd frequencies = Eigen::VectorXd::LinSpaced(N / 2 - low_index + 1, low_index, N / 2);
182 frequencies /= N;
183
184 return std::make_pair(spectrum, frequencies);
185}
186
187Eigen::VectorXd hamming_window(int N)
188{
189 double alpha = 0.54;
190 double beta = 0.46;
191
192 Eigen::VectorXd range = Eigen::VectorXd::LinSpaced(N, 0, 1);
193 Eigen::VectorXd window = alpha - beta * (2 * M_PI * range.array()).cos();
194
195 return window;
196}
197
198double stdandard_deviation(Eigen::VectorXd &input)
199{
200 Eigen::ArrayXd centered = input.array() - input.array().mean();
201 return std::sqrt(centered.pow(2).sum() / (centered.size() - 1));
202}
203
204
205} // namespace math_tools
206
207
208
209
210
QWidget * window(QObject *job)
Eigen::VectorXd hamming_window(int N)
Eigen::MatrixXd squareDistance(const Eigen::MatrixXd &a, const Eigen::MatrixXd &b)
Eigen::MatrixXd generate_normal_random_matrix(const size_t n, const size_t m)
Eigen::MatrixXd box_muller(const Eigen::VectorXd &vRand)
Eigen::MatrixXd generate_uniform_random_matrix_0_1(const size_t n, const size_t m)
std::pair< Eigen::VectorXd, Eigen::VectorXd > compute_spectrum(Eigen::VectorXd &data, int N)
double stdandard_deviation(Eigen::VectorXd &input)
Provides mathematical tools needed for the Gaussian process toolbox.
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.