Kstars

src/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 <[email protected]>
14  * @author Stephan Wenninger <[email protected]>
15  * @author Raffi Enficiaud <[email protected]>
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 
25 namespace math_tools
26 {
27 
28 Eigen::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 
90 Eigen::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 
109 Eigen::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 
138 Eigen::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 
151 std::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 
187 Eigen::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 
198 double 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 
double stdandard_deviation(Eigen::VectorXd &input)
Eigen::MatrixXd squareDistance(const Eigen::MatrixXd &a, const Eigen::MatrixXd &b)
std::pair< Eigen::VectorXd, Eigen::VectorXd > compute_spectrum(Eigen::VectorXd &data, int N)
Eigen::MatrixXd box_muller(const Eigen::VectorXd &vRand)
Eigen::VectorXd hamming_window(int N)
Eigen::MatrixXd generate_normal_random_matrix(const size_t n, const size_t m)
Eigen::MatrixXd generate_uniform_random_matrix_0_1(const size_t n, const size_t m)
Provides mathematical tools needed for the Gaussian process toolbox.
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Sun Oct 1 2023 04:02:41 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.