v2.0.0
Loading...
Searching...
No Matches
numerics.cpp
Go to the documentation of this file.
1//=============================================================================================================
26
27//=============================================================================================================
28// INCLUDES
29//=============================================================================================================
30
31// MSVC builds already define _USE_MATH_DEFINES globally (see src/CMakeLists.txt),
32// so define it here only for the toolchains that do not.
33#ifndef _USE_MATH_DEFINES
34#define _USE_MATH_DEFINES
35#endif
36#include <cmath>
37#include <iostream>
38#include "numerics.h"
39
40//=============================================================================================================
41// QT INCLUDES
42//=============================================================================================================
43
44#include <QDebug>
45#include <QStringList>
46
47//=============================================================================================================
48// USED NAMESPACES
49//=============================================================================================================
50
51using namespace UTILSLIB;
52using namespace Eigen;
53
54//=============================================================================================================
55// DEFINE MEMBER METHODS
56//=============================================================================================================
57
58int Numerics::gcd(int iA, int iB)
59{
60 if (iB == 0) {
61 return iA;
62 }
63
64 return gcd(iB, iA % iB);
65}
66
67//=============================================================================================================
68
69bool Numerics::issparse(VectorXd& v)
70{
71 qint32 c = 0;
72 qint32 n = v.rows();
73 qint32 t = n / 2;
74
75 for (qint32 i = 0; i < n; ++i) {
76 if (v(i) == 0)
77 ++c;
78 if (c > t)
79 return true;
80 }
81
82 return false;
83}
84
85//=============================================================================================================
86
87int Numerics::nchoose2(int n)
88{
89 int t_iNumOfCombination = static_cast<int>(n * (n - 1) * 0.5);
90
91 return t_iNumOfCombination;
92}
93
94//=============================================================================================================
95
96double Numerics::chi2Isf(double p, int dof)
97{
98 // P(a, x) summed from its power series (A&S 6.5.29), the quantile bisected; Q = 1 - P keeps
99 // ~1e-13 relative accuracy for the p = 1e-3 the inverse SNR estimates use.
100 const double a = 0.5 * dof;
101 auto upperTail = [a](double chi2) {
102 const double x = 0.5 * chi2;
103 double term = 1.0;
104 double sum = 1.0;
105 for (int n = 1; n < 100000 && term > 1e-17 * sum; ++n) {
106 term *= x / (a + n);
107 sum += term;
108 }
109 return 1.0 - sum * std::exp(a * std::log(x) - x - std::lgamma(a + 1.0));
110 };
111 double lo = 0.0;
112 double hi = dof + 10.0 * std::sqrt(2.0 * dof) + 50.0;
113 while (upperTail(hi) > p)
114 hi *= 2.0;
115 for (int i = 0; i < 200 && hi - lo > 1e-14 * hi; ++i) {
116 const double mid = 0.5 * (lo + hi);
117 (upperTail(mid) > p ? lo : hi) = mid;
118 }
119 return 0.5 * (lo + hi);
120}
121
122//=============================================================================================================
123
124MatrixXd Numerics::rescale(const MatrixXd& data,
125 const RowVectorXf& times,
126 const QPair<float, float>& baseline,
127 QString mode)
128{
129 MatrixXd data_out = data;
130 QStringList valid_modes;
131 valid_modes << "logratio" << "ratio" << "zscore" << "mean" << "percent";
132 if (!valid_modes.contains(mode)) {
133 qWarning().noquote() << "[Numerics::rescale] Mode" << mode << "is not supported. Supported modes are:" << valid_modes << "Returning input data.";
134 return data_out;
135 }
136
137 qInfo().noquote() << QString("[Numerics::rescale] Applying baseline correction ... (mode: %1)").arg(mode);
138
139 qint32 imin = 0;
140 qint32 imax = times.size();
141
142 if (baseline.second == baseline.first) {
143 imin = 0;
144 } else {
145 float bmin = baseline.first;
146 for (qint32 i = 0; i < times.size(); ++i) {
147 if (times[i] >= bmin) {
148 imin = i;
149 break;
150 }
151 }
152 }
153
154 float bmax = baseline.second;
155
156 if (baseline.second == baseline.first) {
157 bmax = 0;
158 }
159
160 for (qint32 i = times.size() - 1; i >= 0; --i) {
161 if (times[i] <= bmax) {
162 imax = i + 1;
163 break;
164 }
165 }
166
167 if (imax < imin) {
168 qWarning() << "[Numerics::rescale] imax < imin. Returning input data.";
169 return data_out;
170 }
171
172 VectorXd mean = data_out.block(0, imin, data_out.rows(), imax - imin).rowwise().mean();
173 if (mode.compare("mean") == 0) {
174 data_out -= mean.rowwise().replicate(data.cols());
175 } else if (mode.compare("logratio") == 0) {
176 for (qint32 i = 0; i < data_out.rows(); ++i)
177 for (qint32 j = 0; j < data_out.cols(); ++j)
178 data_out(i, j) = log10(data_out(i, j) / mean[i]);
179 } else if (mode.compare("ratio") == 0) {
180 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
181 } else if (mode.compare("zscore") == 0) {
182 MatrixXd std_mat = data.block(0, imin, data.rows(), imax - imin) - mean.rowwise().replicate(imax - imin);
183 std_mat = std_mat.cwiseProduct(std_mat);
184 VectorXd std_v = std_mat.rowwise().mean();
185 for (qint32 i = 0; i < std_v.size(); ++i)
186 std_v[i] = sqrt(std_v[i] / static_cast<float>(imax - imin));
187
188 data_out -= mean.rowwise().replicate(data_out.cols());
189 data_out = data_out.cwiseQuotient(std_v.rowwise().replicate(data_out.cols()));
190 } else if (mode.compare("percent") == 0) {
191 data_out -= mean.rowwise().replicate(data_out.cols());
192 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
193 }
194
195 return data_out;
196}
197
198//=============================================================================================================
199
200MatrixXd Numerics::rescale(const MatrixXd& data,
201 const RowVectorXf& times,
202 const std::pair<float, float>& baseline,
203 const std::string& mode)
204{
205 MatrixXd data_out = data;
206 std::vector<std::string> valid_modes{"logratio", "ratio", "zscore", "mean", "percent"};
207 if (std::find(valid_modes.begin(), valid_modes.end(), mode) == valid_modes.end()) {
208 qWarning().noquote() << "[Numerics::rescale] Mode" << mode.c_str() << "is not supported. Supported modes are:";
209 for (auto& m : valid_modes) {
210 std::cout << m << " ";
211 }
212 std::cout << "\n"
213 << "Returning input data.\n";
214 return data_out;
215 }
216
217 qInfo().noquote() << QString("[Numerics::rescale] Applying baseline correction ... (mode: %1)").arg(mode.c_str());
218
219 qint32 imin = 0;
220 qint32 imax = times.size();
221
222 if (baseline.second == baseline.first) {
223 imin = 0;
224 } else {
225 float bmin = baseline.first;
226 for (qint32 i = 0; i < times.size(); ++i) {
227 if (times[i] >= bmin) {
228 imin = i;
229 break;
230 }
231 }
232 }
233
234 float bmax = baseline.second;
235
236 if (baseline.second == baseline.first) {
237 bmax = 0;
238 }
239
240 for (qint32 i = times.size() - 1; i >= 0; --i) {
241 if (times[i] <= bmax) {
242 imax = i + 1;
243 break;
244 }
245 }
246
247 if (imax < imin) {
248 qWarning() << "[Numerics::rescale] imax < imin. Returning input data.";
249 return data_out;
250 }
251
252 VectorXd mean = data_out.block(0, imin, data_out.rows(), imax - imin).rowwise().mean();
253 if (mode.compare("mean") == 0) {
254 data_out -= mean.rowwise().replicate(data.cols());
255 } else if (mode.compare("logratio") == 0) {
256 for (qint32 i = 0; i < data_out.rows(); ++i)
257 for (qint32 j = 0; j < data_out.cols(); ++j)
258 data_out(i, j) = log10(data_out(i, j) / mean[i]);
259 } else if (mode.compare("ratio") == 0) {
260 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
261 } else if (mode.compare("zscore") == 0) {
262 MatrixXd std_mat = data.block(0, imin, data.rows(), imax - imin) - mean.rowwise().replicate(imax - imin);
263 std_mat = std_mat.cwiseProduct(std_mat);
264 VectorXd std_v = std_mat.rowwise().mean();
265 for (qint32 i = 0; i < std_v.size(); ++i)
266 std_v[i] = sqrt(std_v[i] / static_cast<float>(imax - imin));
267
268 data_out -= mean.rowwise().replicate(data_out.cols());
269 data_out = data_out.cwiseQuotient(std_v.rowwise().replicate(data_out.cols()));
270 } else if (mode.compare("percent") == 0) {
271 data_out -= mean.rowwise().replicate(data_out.cols());
272 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
273 }
274
275 return data_out;
276}
numerics.h
General numerical helpers: GCD, log2, histogram binning, baseline rescaling, sparsity tests.
UTILSLIB
Shared utilities (I/O helpers, spectral analysis, layout management, warp algorithms).
Definition rtfiffrawview.h:60
UTILSLIB::Numerics::gcd
static int gcd(int iA, int iB)
Definition numerics.cpp:58
UTILSLIB::Numerics::chi2Isf
static double chi2Isf(double p, int dof)
Definition numerics.cpp:96
UTILSLIB::Numerics::issparse
static bool issparse(Eigen::VectorXd &v)
Definition numerics.cpp:69
UTILSLIB::Numerics::rescale
static Eigen::MatrixXd rescale(const Eigen::MatrixXd &data, const Eigen::RowVectorXf &times, const QPair< float, float > &baseline, QString mode)
UTILSLIB::Numerics::nchoose2
static int nchoose2(int n)
Definition numerics.cpp:87
Generated on for MNE-CPP by doxygen 1.16.1