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