v2.0.0
Loading...
Searching...
No Matches
numerics.cpp
Go to the documentation of this file.
1//=============================================================================================================
26
27//=============================================================================================================
28// INCLUDES
29//=============================================================================================================
30
31#define _USE_MATH_DEFINES
32#include <cmath>
33#include <iostream>
34#include "numerics.h"
35
36//=============================================================================================================
37// QT INCLUDES
38//=============================================================================================================
39
40#include <QDebug>
41#include <QStringList>
42
43//=============================================================================================================
44// USED NAMESPACES
45//=============================================================================================================
46
47using namespace UTILSLIB;
48using namespace Eigen;
49
50//=============================================================================================================
51// DEFINE MEMBER METHODS
52//=============================================================================================================
53
54int Numerics::gcd(int iA, int iB)
55{
56 if (iB == 0) {
57 return iA;
58 }
59
60 return gcd(iB, iA % iB);
61}
62
63//=============================================================================================================
64
65bool Numerics::issparse(VectorXd &v)
66{
67 qint32 c = 0;
68 qint32 n = v.rows();
69 qint32 t = n/2;
70
71 for(qint32 i = 0; i < n; ++i)
72 {
73 if(v(i) == 0)
74 ++c;
75 if(c > t)
76 return true;
77 }
78
79 return false;
80}
81
82//=============================================================================================================
83
84int Numerics::nchoose2(int n)
85{
86 int t_iNumOfCombination = static_cast<int>(n*(n-1)*0.5);
87
88 return t_iNumOfCombination;
89}
90
91//=============================================================================================================
92
93MatrixXd Numerics::rescale(const MatrixXd &data,
94 const RowVectorXf &times,
95 const QPair<float,float>& baseline,
96 QString mode)
97{
98 MatrixXd data_out = data;
99 QStringList valid_modes;
100 valid_modes << "logratio" << "ratio" << "zscore" << "mean" << "percent";
101 if(!valid_modes.contains(mode))
102 {
103 qWarning().noquote() << "[Numerics::rescale] Mode" << mode << "is not supported. Supported modes are:" << valid_modes << "Returning input data.";
104 return data_out;
105 }
106
107 qInfo().noquote() << QString("[Numerics::rescale] Applying baseline correction ... (mode: %1)").arg(mode);
108
109 qint32 imin = 0;
110 qint32 imax = times.size();
111
112 if (baseline.second == baseline.first) {
113 imin = 0;
114 } else {
115 float bmin = baseline.first;
116 for(qint32 i = 0; i < times.size(); ++i) {
117 if(times[i] >= bmin) {
118 imin = i;
119 break;
120 }
121 }
122 }
123
124 float bmax = baseline.second;
125
126 if (baseline.second == baseline.first) {
127 bmax = 0;
128 }
129
130 for(qint32 i = times.size()-1; i >= 0; --i) {
131 if(times[i] <= bmax) {
132 imax = i+1;
133 break;
134 }
135 }
136
137 if(imax < imin) {
138 qWarning() << "[Numerics::rescale] imax < imin. Returning input data.";
139 return data_out;
140 }
141
142 VectorXd mean = data_out.block(0, imin, data_out.rows(), imax-imin).rowwise().mean();
143 if(mode.compare("mean") == 0) {
144 data_out -= mean.rowwise().replicate(data.cols());
145 } else if(mode.compare("logratio") == 0) {
146 for(qint32 i = 0; i < data_out.rows(); ++i)
147 for(qint32 j = 0; j < data_out.cols(); ++j)
148 data_out(i,j) = log10(data_out(i,j)/mean[i]);
149 } else if(mode.compare("ratio") == 0) {
150 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
151 } else if(mode.compare("zscore") == 0) {
152 MatrixXd std_mat = data.block(0, imin, data.rows(), imax-imin) - mean.rowwise().replicate(imax-imin);
153 std_mat = std_mat.cwiseProduct(std_mat);
154 VectorXd std_v = std_mat.rowwise().mean();
155 for(qint32 i = 0; i < std_v.size(); ++i)
156 std_v[i] = sqrt(std_v[i] / static_cast<float>(imax-imin));
157
158 data_out -= mean.rowwise().replicate(data_out.cols());
159 data_out = data_out.cwiseQuotient(std_v.rowwise().replicate(data_out.cols()));
160 } else if(mode.compare("percent") == 0) {
161 data_out -= mean.rowwise().replicate(data_out.cols());
162 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
163 }
164
165 return data_out;
166}
167
168//=============================================================================================================
169
170MatrixXd Numerics::rescale(const MatrixXd& data,
171 const RowVectorXf& times,
172 const std::pair<float,float>& baseline,
173 const std::string& mode)
174{
175 MatrixXd data_out = data;
176 std::vector<std::string> valid_modes{"logratio", "ratio", "zscore", "mean", "percent"};
177 if(std::find(valid_modes.begin(), valid_modes.end(), mode) == valid_modes.end())
178 {
179 qWarning().noquote() << "[Numerics::rescale] Mode" << mode.c_str() << "is not supported. Supported modes are:";
180 for(auto& m : valid_modes){
181 std::cout << m << " ";
182 }
183 std::cout << "\n" << "Returning input data.\n";
184 return data_out;
185 }
186
187 qInfo().noquote() << QString("[Numerics::rescale] Applying baseline correction ... (mode: %1)").arg(mode.c_str());
188
189 qint32 imin = 0;
190 qint32 imax = times.size();
191
192 if (baseline.second == baseline.first) {
193 imin = 0;
194 } else {
195 float bmin = baseline.first;
196 for(qint32 i = 0; i < times.size(); ++i) {
197 if(times[i] >= bmin) {
198 imin = i;
199 break;
200 }
201 }
202 }
203
204 float bmax = baseline.second;
205
206 if (baseline.second == baseline.first) {
207 bmax = 0;
208 }
209
210 for(qint32 i = times.size()-1; i >= 0; --i) {
211 if(times[i] <= bmax) {
212 imax = i+1;
213 break;
214 }
215 }
216
217 if(imax < imin) {
218 qWarning() << "[Numerics::rescale] imax < imin. Returning input data.";
219 return data_out;
220 }
221
222 VectorXd mean = data_out.block(0, imin, data_out.rows(), imax-imin).rowwise().mean();
223 if(mode.compare("mean") == 0) {
224 data_out -= mean.rowwise().replicate(data.cols());
225 } else if(mode.compare("logratio") == 0) {
226 for(qint32 i = 0; i < data_out.rows(); ++i)
227 for(qint32 j = 0; j < data_out.cols(); ++j)
228 data_out(i,j) = log10(data_out(i,j)/mean[i]);
229 } else if(mode.compare("ratio") == 0) {
230 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
231 } else if(mode.compare("zscore") == 0) {
232 MatrixXd std_mat = data.block(0, imin, data.rows(), imax-imin) - mean.rowwise().replicate(imax-imin);
233 std_mat = std_mat.cwiseProduct(std_mat);
234 VectorXd std_v = std_mat.rowwise().mean();
235 for(qint32 i = 0; i < std_v.size(); ++i)
236 std_v[i] = sqrt(std_v[i] / static_cast<float>(imax-imin));
237
238 data_out -= mean.rowwise().replicate(data_out.cols());
239 data_out = data_out.cwiseQuotient(std_v.rowwise().replicate(data_out.cols()));
240 } else if(mode.compare("percent") == 0) {
241 data_out -= mean.rowwise().replicate(data_out.cols());
242 data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
243 }
244
245 return data_out;
246}
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:58
UTILSLIB::Numerics::gcd
static int gcd(int iA, int iB)
Definition numerics.cpp:54
UTILSLIB::Numerics::issparse
static bool issparse(Eigen::VectorXd &v)
Definition numerics.cpp:65
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:84
Generated on for MNE-CPP by doxygen 1.16.1