1 #include "ScatterDraw.h"
2
3 namespace Upp {
4 using namespace Eigen;
5
6 #define Membercall(fun) (this->*fun)
7
Min(Getdatafun getdata,int64 & id)8 double DataSource::Min(Getdatafun getdata, int64& id) {
9 double minVal = -DOUBLE_NULL;
10 for (int64 i = 0; i < GetCount(); ++i) {
11 double d = Membercall(getdata)(i);
12 if (!IsNull(d) && minVal > d) {
13 minVal = d;
14 id = i;
15 }
16 }
17 if (minVal == -DOUBLE_NULL)
18 return Null;
19 return minVal;
20 }
21
Max(Getdatafun getdata,int64 & id)22 double DataSource::Max(Getdatafun getdata, int64& id) {
23 double maxVal = DOUBLE_NULL;
24 for (int64 i = 0; i < GetCount(); ++i) {
25 double d = Membercall(getdata)(i);
26 if (!IsNull(d) && maxVal < d) {
27 maxVal = d;
28 id = i;
29 }
30 }
31 if (maxVal == DOUBLE_NULL)
32 return Null;
33 return maxVal;
34 }
35
MaxList(Getdatafun getdataY,Getdatafun getdataX,Vector<int64> & id,double width)36 void DataSource::MaxList(Getdatafun getdataY, Getdatafun getdataX, Vector<int64> &id, double width) {
37 id.Clear();
38 for (int64 i = 1; i < GetCount() - 1; ++i) {
39 double d = Membercall(getdataY)(i);
40 if (IsNull(d))
41 continue;
42 int64 ii;
43 for (ii = i-1; ii >= 0; --ii) {
44 if (!IsNull(Membercall(getdataY)(ii)))
45 break;
46 }
47 if (ii < 0)
48 continue;
49 double d_1 = Membercall(getdataY)(ii);
50 for (ii = i+1; ii < GetCount(); ++ii) {
51 if (!IsNull(Membercall(getdataY)(ii)))
52 break;
53 }
54 if (ii >= GetCount())
55 continue;
56 double d1 = Membercall(getdataY)(ii);
57 if (d >= d_1 && d > d1) {
58 if (id.GetCount() == 0 ||
59 (Membercall(getdataX)(i) - Membercall(getdataX)(id.GetCount() - 1) >= width))
60 id << i;
61 }
62 }
63 }
64
MaxSubDataImp(Getdatafun getdataY,Getdatafun getdataX,int64 maxId,int64 width)65 Pointf DataSource::MaxSubDataImp(Getdatafun getdataY, Getdatafun getdataX, int64 maxId, int64 width) {
66 Vector<Pointf> p;
67
68 int iw;
69 int64 idMin, idMax;
70 iw = 0;
71 for (idMin = maxId - 1; idMin >= 0 && iw < width; idMin--) {
72 if (IsNull(Membercall(getdataY)(idMin)) || IsNull(Membercall(getdataX)(idMin)))
73 continue;
74 iw++;
75 }
76 if (idMin < 0)
77 idMin = 0;
78 iw = 0;
79 for (idMax = maxId + 1; idMax < GetCount() && iw < width; idMax++) {
80 if (IsNull(Membercall(getdataY)(idMax)) || IsNull(Membercall(getdataX)(idMax)))
81 continue;
82 iw++;
83 }
84 if (idMax >= GetCount())
85 idMax = GetCount() - 1;
86 for (int64 i = idMin; i <= idMax; ++i) {
87 if (IsNull(Membercall(getdataY)(i)) || IsNull(Membercall(getdataX)(i)))
88 continue;
89 p << Pointf(Membercall(getdataX)(i), Membercall(getdataY)(i));
90 }
91 VectorPointf pf(p);
92 PolynomialEquation2 polyFit;
93 if (ExplicitEquation::NoError != polyFit.Fit(pf))
94 return Null;
95 const Vector<double> &coeff = polyFit.GetCoeff();
96 double b = coeff[1];
97 double a = coeff[2];
98 if (IsNull(a) || fabs(a) < 1E-10)
99 return Null;
100 return Pointf(-b/2/a, polyFit.f(-b/2/a));
101 }
102
Avg(Getdatafun getdata)103 double DataSource::Avg(Getdatafun getdata) {
104 double ret = 0;
105 int count = 0;
106 for (int64 i = 0; i < GetCount(); ++i) {
107 double d = Membercall(getdata)(i);
108 if (!IsNull(d)) {
109 ret += d;
110 count++;
111 }
112 }
113 if (count == 0)
114 return Null;
115 return ret/count;
116 }
117
Closest(Getdatafun getdata,double dat)118 int64 DataSource::Closest(Getdatafun getdata, double dat) {
119 double minD = DBL_MAX;
120 int64 minDat;
121 for (int64 i = 0; i < GetCount(); ++i) {
122 double d = Membercall(getdata)(i);
123 if (!IsNull(d)) {
124 if (minD > abs(d - dat)) {
125 minD = abs(d - dat);
126 minDat = i;
127 }
128 }
129 }
130 if (minD == DBL_MAX)
131 return Null;
132 return minDat;
133 }
134
Closest(Getdatafun getdataX,Getdatafun getdataY,double x,double y)135 int64 DataSource::Closest(Getdatafun getdataX, Getdatafun getdataY, double x, double y) {
136 double minD = DBL_MAX;
137 int64 minDat;
138 for (int64 i = 0; i < GetCount(); ++i) {
139 double dx = Membercall(getdataX)(i);
140 double dy = Membercall(getdataY)(i);
141 if (!IsNull(dx) && !IsNull(dy)) {
142 double d = sqr(dx - x) + sqr(dy - y);
143 if (minD > d) {
144 minD = d;
145 minDat = i;
146 }
147 }
148 }
149 if (minD == DBL_MAX)
150 return Null;
151 return minDat;
152 }
153
RMS(Getdatafun getdata)154 double DataSource::RMS(Getdatafun getdata) {
155 double ret = 0;
156 int count = 0;
157 for (int64 i = 0; i < GetCount(); ++i) {
158 double d = Membercall(getdata)(i);
159 if (!IsNull(d)) {
160 ret += d*d;
161 count++;
162 }
163 }
164 if (count == 0)
165 return Null;
166 return sqrt(ret/count);
167 }
168
IsSorted(Getdatafun getdata)169 double DataSource::IsSorted(Getdatafun getdata) {
170 int64 num = GetCount();
171 if (num == 0)
172 return false;
173 if (num == 1)
174 return 1;
175 for (int i = 1; i < num; ++i) {
176 if (Membercall(getdata)(i) < Membercall(getdata)(i - 1))
177 return false;
178 }
179 return true;
180 }
181
Variance(Getdatafun getdata,double avg)182 double DataSource::Variance(Getdatafun getdata, double avg) {
183 if (IsNull(avg))
184 avg = Avg(getdata);
185 if (IsNull(avg))
186 return Null;
187 double ret = 0;
188 int count = 0;
189 for (int64 i = 0; i < GetCount(); ++i) {
190 double d = Membercall(getdata)(i);
191 if (!IsNull(d)) {
192 d -= avg;
193 ret += d*d;
194 count++;
195 }
196 }
197 if (count <= 1)
198 return Null;
199 return ret/(count - 1);
200 }
201
Envelope(Getdatafun getdataY,Getdatafun getdataX,double width,bool (* fun)(double a,double b))202 Vector<int64> DataSource::Envelope(Getdatafun getdataY, Getdatafun getdataX, double width, bool (*fun)(double a, double b)) {
203 Vector<int64> ret;
204 double width_2 = width/2.;
205
206 for (int i = 0; i < GetCount(); ++i) {
207 double y = Membercall(getdataY)(i);
208 double x = Membercall(getdataX)(i);
209 if (IsNull(x) || IsNull(y))
210 continue;
211 int numComparisons = 0;
212 for (int j = i-1; j >= 0; --j) {
213 double ynext = Membercall(getdataY)(j);
214 double xnext = Membercall(getdataX)(j);
215 if (IsNull(xnext) || IsNull(ynext))
216 continue;
217 if ((x - xnext) > width_2)
218 break;
219 if (!fun(y, ynext)) {
220 numComparisons = Null;
221 break;
222 }
223 numComparisons++;
224 }
225 if (IsNull(numComparisons))
226 continue;
227 for (int j = i+1; j < GetCount(); ++j) {
228 double ynext = Membercall(getdataY)(j);
229 double xnext = Membercall(getdataX)(j);
230 if (IsNull(xnext) || IsNull(ynext))
231 continue;
232 if ((xnext - x) > width_2)
233 break;
234 if (!fun(y, ynext)) {
235 numComparisons = Null;
236 break;
237 }
238 numComparisons++;
239 }
240 if (IsNull(numComparisons))
241 continue;
242 if (numComparisons > 2) {
243 if (!ret.IsEmpty()) {
244 int64 prev_i = ret[ret.GetCount() - 1];
245 if (Membercall(getdataX)(prev_i) != x)
246 ret << i;
247 } else
248 ret << i;
249 }
250 }
251 return ret;
252 }
253
GreaterEqualThan(double a,double b)254 inline bool GreaterEqualThan(double a, double b) {return a >= b;}
LowerEqualThan(double a,double b)255 inline bool LowerEqualThan(double a, double b) {return a <= b;}
256
UpperEnvelope(Getdatafun getdataY,Getdatafun getdataX,double width)257 Vector<int64> DataSource::UpperEnvelope(Getdatafun getdataY, Getdatafun getdataX, double width) {return Envelope(getdataY, getdataX, width, GreaterEqualThan);}
LowerEnvelope(Getdatafun getdataY,Getdatafun getdataX,double width)258 Vector<int64> DataSource::LowerEnvelope(Getdatafun getdataY, Getdatafun getdataX, double width) {return Envelope(getdataY, getdataX, width, LowerEqualThan);}
259
Cumulative(Getdatafun getdataY,Getdatafun getdataX)260 Vector<Pointf> DataSource::Cumulative(Getdatafun getdataY, Getdatafun getdataX) {
261 Vector<Pointf> ret;
262
263 double acum = 0;
264 for (int i = 0; i < GetCount(); ++i) {
265 double y = Membercall(getdataY)(i);
266 double x = Membercall(getdataX)(i);
267
268 if (IsNull(x) || IsNull(y))
269 continue;
270 acum += y;
271 ret << Pointf(x, acum);
272 }
273 return ret;
274 }
275
CumulativeAverage(Getdatafun getdataY,Getdatafun getdataX)276 Vector<Pointf> DataSource::CumulativeAverage(Getdatafun getdataY, Getdatafun getdataX) {
277 Vector<Pointf> ret;
278
279 double acum = 0;
280 int num = 0;
281 for (int i = 0; i < GetCount(); ++i) {
282 double y = Membercall(getdataY)(i);
283 double x = Membercall(getdataX)(i);
284
285 if (IsNull(x) || IsNull(y))
286 continue;
287 acum += y;
288 num++;
289 ret << Pointf(x, acum/num);
290 }
291 return ret;
292 }
293
MovingAverage(Getdatafun getdataY,Getdatafun getdataX,double width)294 Vector<Pointf> DataSource::MovingAverage(Getdatafun getdataY, Getdatafun getdataX, double width) {
295 Vector<Pointf> ret;
296 double width_2 = width/2.;
297
298 for (int i = 0; i < GetCount(); ++i) {
299 double y = Membercall(getdataY)(i);
300 double x = Membercall(getdataX)(i);
301 if (IsNull(x) || IsNull(y))
302 continue;
303 int numAvg = 1;
304 double sum = y;
305 int j;
306 for (j = i-1; j >= 0; --j) {
307 double ynext = Membercall(getdataY)(j);
308 double xnext = Membercall(getdataX)(j);
309 if (IsNull(xnext) || IsNull(ynext))
310 continue;
311 if ((x - xnext) > width_2)
312 break;
313 sum += ynext;
314 numAvg++;
315 }
316 if (j < 0)
317 continue;
318 for (j = i+1; j < GetCount(); ++j) {
319 double ynext = Membercall(getdataY)(j);
320 double xnext = Membercall(getdataX)(j);
321 if (IsNull(xnext))
322 continue;
323 if ((xnext - x) > width_2)
324 break;
325 if (IsNull(ynext))
326 continue;
327 sum += ynext;
328 numAvg++;
329 }
330 if (j == GetCount())
331 continue;
332 ret << Pointf(x, sum/numAvg);
333 }
334 return ret;
335 }
336
SectorAverage(Getdatafun getdataY,Getdatafun getdataX,double width)337 Vector<Pointf> DataSource::SectorAverage(Getdatafun getdataY, Getdatafun getdataX, double width) {
338 Vector<Pointf> ret;
339
340 for (int i = 0; i < GetCount();) {
341 double y = Membercall(getdataY)(i);
342 double x = Membercall(getdataX)(i);
343 if (IsNull(x) || IsNull(y))
344 continue;
345
346 int numAvg = 1;
347 double sum = y;
348 double sumX = x;
349 int j;
350 for (j = i+1; j < GetCount(); ++j) {
351 double ynext = Membercall(getdataY)(j);
352 double xnext = Membercall(getdataX)(j);
353 if (IsNull(xnext))
354 continue;
355 if ((xnext - x) > width) {
356 --j;
357 break;
358 }
359 if (IsNull(ynext))
360 continue;
361 sumX += xnext;
362 sum += ynext;
363 numAvg++;
364 }
365 ret << Pointf(sumX/numAvg, sum/numAvg);
366 if (j == GetCount())
367 break;
368 i = j+1;
369 }
370 return ret;
371 }
372
ZeroCrossing(Getdatafun getdataY,Getdatafun getdataX,bool ascending,bool descending,Vector<double> & zeros,Vector<int64> & ids)373 void DataSource::ZeroCrossing(Getdatafun getdataY, Getdatafun getdataX, bool ascending, bool descending,
374 Vector<double> &zeros, Vector<int64> &ids) {
375 zeros.Clear();
376 ids.Clear();
377
378 double y_prev = 0, x_prev = 0;
379 int i0;
380 for (i0 = 0; i0 < GetCount(); ++i0) {
381 y_prev = Membercall(getdataY)(i0);
382 x_prev = Membercall(getdataX)(i0);
383 if (!IsNull(x_prev) && !IsNull(y_prev))
384 break;
385 }
386 for (int i = i0; i < GetCount(); ++i) {
387 double y = Membercall(getdataY)(i);
388 double x = Membercall(getdataX)(i);
389 if (IsNull(x) || IsNull(y))
390 continue;
391
392 if (((y >= 0 && y_prev < 0) && ascending) || ((y <= 0 && y_prev > 0) && descending)) {
393 ids << i;
394 zeros << (x_prev - (x - x_prev)*y_prev/(y - y_prev));
395 }
396 x_prev = x;
397 y_prev = y;
398 }
399 }
400
StdDev(Getdatafun getdata,double avg)401 double DataSource::StdDev(Getdatafun getdata, double avg) {
402 double var = Variance(getdata, avg);
403 return IsNull(var) ? Null : sqrt(var);
404 }
405
SinEstim_Amplitude(double avg)406 double DataSource::SinEstim_Amplitude(double avg) {
407 return sqrt(2.*VarianceY(avg));
408 }
409
SinEstim_FreqPhase(double & frequency,double & phase,double avg)410 bool DataSource::SinEstim_FreqPhase(double &frequency, double &phase, double avg) {
411 if (GetCount() < 4)
412 return false;
413 if (IsNull(avg))
414 avg = AvgY();
415 int64 firstId;
416 for (firstId = 0; firstId < GetCount(); ++firstId)
417 if (!IsNull(x(firstId)) && !IsNull(y(firstId)))
418 break;
419 bool firstIsToPositive = (y(firstId) - avg) < 0;
420 bool isPossitive = !firstIsToPositive;
421 double T = 0;
422 int numT = 0;
423 double lastZero = Null;
424 double firstZero = Null;
425 firstId++;
426 for (int64 id = firstId; id < GetCount(); ++id) {
427 if (IsNull(x(id)) || IsNull(y(id)))
428 continue;
429 if (((y(id) - avg) > 0) != isPossitive) {
430 isPossitive = !isPossitive;
431 double zero = x(id-1) - (y(id-1) - avg)*(x(id) - x(id-1))/(y(id) - y(id-1));
432 if (IsNull(lastZero))
433 firstZero = zero;
434 else {
435 T += zero - lastZero;
436 numT++;
437 }
438 lastZero = zero;
439 }
440 }
441 if (T == 0 || numT == 0)
442 return false;
443 T = 2*T/numT;
444 frequency = 2*M_PI/T;
445 phase = -frequency*firstZero;
446 if (!firstIsToPositive)
447 phase += M_PI;
448 phase = phase - 2*M_PI*int(phase/(2*M_PI));
449 if (phase > M_PI)
450 phase = phase - 2*M_PI;
451 if (phase < -M_PI)
452 phase = 2*M_PI + phase;
453 return true;
454 }
455
456
znFixed(int n,int64 id)457 double CArray::znFixed(int n, int64 id) {
458 if (n == 0)
459 return zData[id];
460 NEVER();
461 return Null;
462 }
463
FFTSimple(VectorXd & data,double tSample,bool frequency,int type,int window,int numOver)464 Vector<Pointf> FFTSimple(VectorXd &data, double tSample, bool frequency, int type,
465 int window, int numOver) {
466 int numData = int(data.size());
467
468 double numDataFact = 0;
469 switch (window) {
470 case DataSource::HAMMING:
471 for (int i = 0; i < numData; ++i) {
472 double windowFact = 0.54 - 0.46*cos(2*M_PI*i/numData);
473 numDataFact += windowFact;
474 data[i] *= windowFact;
475 }
476 break;
477 case DataSource::COS:
478 for (int i = 0; i < numOver; ++i) {
479 double windowFact = 0.5*(1 - cos(M_PI*i/numOver));
480 numDataFact += windowFact;
481 data[i] *= windowFact;
482 }
483 for (int i = numOver; i < numData - numOver; ++i) {
484 double windowFact = 1; // 1.004
485 numDataFact += windowFact;
486 //data[i] *= windowFact;
487 }
488 for (int i = numData - numOver; i < numData; ++i) {
489 double windowFact = 0.5*(1 + cos(M_PI*(numData - i - numOver)/numOver));
490 numDataFact += windowFact;
491 data[i] *= windowFact;
492 }
493 break;
494 default: numDataFact = numData;
495 }
496 Vector<Pointf> res;
497 VectorXcd freqbuf;
498 try {
499 FFT<double> fft;
500 fft.SetFlag(fft.HalfSpectrum);
501 fft.fwd(freqbuf, data);
502 } catch(...) {
503 return res;
504 }
505 double threshold = 0;
506 if (type == DataSource::T_PHASE) {
507 for (int i = 0; i < int(freqbuf.size()); ++i) {
508 if (threshold < std::abs(freqbuf[i]))
509 threshold = std::abs(freqbuf[i]);
510 }
511 }
512 threshold /= 10000.;
513 if (frequency) {
514 for (int i = 0; i < int(freqbuf.size()); ++i) {
515 double xdata = i/(tSample*numData);
516
517 switch (type) {
518 case DataSource::T_PHASE:
519 if (std::abs(freqbuf[i]) > threshold)
520 res << Pointf(xdata, std::arg(freqbuf[i]));
521 else
522 res << Pointf(xdata, 0);
523 break;
524 case DataSource::T_FFT: // Amplitude spectrum
525 res << Pointf(xdata, 2*std::abs(freqbuf[i])/numDataFact); break;
526 case DataSource::T_PSD: // Variance density spectrum
527 res << Pointf(xdata, 2*sqr(std::abs(freqbuf[i]))/(numDataFact/tSample)); // 1/2*FFT^2
528 }
529 }
530 } else {
531 for (int i = int(freqbuf.size()) - 1; i > 0; --i) {
532 double xdata = (tSample*numData)/i;
533
534 switch (type) {
535 case DataSource::T_PHASE:
536 if (std::abs(freqbuf[i]) > threshold)
537 res << Pointf(xdata, std::arg(freqbuf[i]));
538 else
539 res << Pointf(xdata, 0);
540 break;
541 case DataSource::T_FFT:
542 res << Pointf(xdata, 2*std::abs(freqbuf[i])/numDataFact); break;
543 case DataSource::T_PSD:
544 res << Pointf(xdata, 2*sqr(std::abs(freqbuf[i]))/(numDataFact/tSample));
545 }
546 }
547 }
548 return res;
549 }
550
FFT(Getdatafun getdata,double tSample,bool frequency,int type,int window,int numSub,double overlapping)551 Vector<Pointf> DataSource::FFT(Getdatafun getdata, double tSample, bool frequency, int type,
552 int window, int numSub, double overlapping) {
553 int numData = int(GetCount());
554
555 VectorXd data(numData);
556 int num = 0;
557 for (int i = 0; i < numData; ++i) {
558 double d = Membercall(getdata)(i);
559 if (!IsNull(d)) {
560 data[i] = d;
561 num++;
562 }
563 }
564 numData = num;
565
566 Vector<Pointf> res;
567 if (num < 3)
568 return res;
569 data.resize(numData);
570 double numOver;
571
572 if (numSub == 1) {
573 numOver = numData*overlapping;
574 return FFTSimple(data, tSample, frequency, type, window, int(numOver));
575 } else { // solve v t=2*(v-f*v/2) + (n-2)*(v-f*v) ==> v=t/(f + n -f*n)
576 double numDataPart = numData/(overlapping + numSub - overlapping*numSub);
577 int inumDataPart = int(numDataPart);
578 numOver = numDataPart*overlapping;
579 VectorXd dataPart(inumDataPart);
580 double izero = 0;
581 int izerod = 0;
582 Vector<Pointf> fft;
583 for (int iPart = 0; iPart < numSub; ++iPart) {
584 if (iPart > 0) {
585 izero += int(numDataPart - numOver);
586 izerod = int(izero);
587 }
588 for (int i = 0; i < inumDataPart; ++i)
589 dataPart[i] = data[izerod + i];
590 Vector<Pointf> fftPart;
591 fftPart = FFTSimple(dataPart, tSample, frequency, type, window, int(numOver));
592 if (iPart == 0)
593 fft = clone(fftPart); // pick()
594 else {
595 for (int i = 0; i < fftPart.GetCount(); ++i) {
596 fft[i].y += fftPart[i].y;
597 ASSERT(fft[i].x == fftPart[i].x);
598 }
599 }
600 }
601 for (int i = 0; i < fft.GetCount(); ++i)
602 fft[i].y /= numSub;
603 return fft;
604 }
605 }
606
Integral(double from,double to,double n)607 double DataSource::Integral(double from, double to, double n) {
608 double h = (to - from)/n;
609 double h_2 = h/2;
610
611 double sum1 = 0, sum2 = 0;
612 for(int i = 0; i < n; i++) {
613 sum1 += f(from + h*i + h_2);
614 sum2 += f(from + h*i);
615 }
616 return h/6*(f(from) + f(to) + 4*sum1 + 2*sum2);
617 }
618
Integral(Getdatafun getdataY,Getdatafun getdataX)619 double DataSource::Integral(Getdatafun getdataY, Getdatafun getdataX) {
620 double prevx = Membercall(getdataX)(0);
621 double prevy = Membercall(getdataY)(0);
622 double sum = 0;
623 for (int i = 1; i < GetCount(); ++i) {
624 double x = Membercall(getdataX)(i);
625 double y = Membercall(getdataY)(i);
626 sum += (x - prevx)*(y + prevy);
627 prevx = x;
628 prevy = y;
629 }
630 return sum/2;
631 }
632
GetSpectralMoments(double from,double to,double n,bool frequency,double & m_1,double & m0,double & m1,double & m2)633 void DataSource::GetSpectralMoments(double from, double to, double n, bool frequency,
634 double &m_1, double &m0, double &m1, double &m2) {
635 if (!frequency) {
636 from = 1/to;
637 to = 1/from;
638 }
639 double h = (to - from)/n;
640 double h2 = h/2;
641
642 double sum1_m_1 = 0, sum2_m_1 = 0;
643 double sum1_m0 = 0, sum2_m0 = 0;
644 double sum1_m1 = 0, sum2_m1 = 0;
645 double sum1_m2 = 0, sum2_m2 = 0;
646 double f1, f2, x1, x2;
647 for(int i = 0; i < n; i++) {
648 if (frequency) {
649 x1 = from + h*i + h2;
650 x2 = from + h*i;
651 f1 = f(x1);
652 f2 = f(x2);
653 } else {
654 x1 = 1/(from + h*i + h2);
655 x2 = 1/(from + h*i);
656 f1 = f(1/x1);
657 f2 = f(1/x2);
658 }
659 sum1_m_1 += f1/x1;
660 sum2_m_1 += f2/x2;
661 sum1_m0 += f1;
662 sum2_m0 += f2;
663 sum1_m1 += f1*x1;
664 sum2_m1 += f2*x2;
665 sum1_m2 += f1*x1*x1;
666 sum2_m2 += f2*x2*x2;
667 }
668 double f_from, f_to;
669 if (frequency) {
670 f_from = f(from);
671 f_to = f(to);
672 } else {
673 f_from = f(1/from);
674 f_to = f(1/to);
675 }
676 m_1 = h/6*(f_from/from + f_to/to + 4*sum1_m_1 + 2*sum2_m_1);
677 m0 = h/6*(f_from + f(to) + 4*sum1_m0 + 2*sum2_m0);
678 m1 = h/6*(f_from*from + f_to*to + 4*sum1_m1 + 2*sum2_m1);
679 m2 = h/6*(f_from*from*from + f_to*to*to + 4*sum1_m2 + 2*sum2_m2);
680 }
681
GetSpectralMoments(Getdatafun getdataY,Getdatafun getdataX,bool frequency,double & m_1,double & m0,double & m1,double & m2)682 void DataSource::GetSpectralMoments(Getdatafun getdataY, Getdatafun getdataX, bool frequency,
683 double &m_1, double &m0, double &m1, double &m2) {
684 double prevx = Membercall(getdataX)(0);
685 double Si_1 = Membercall(getdataY)(0);
686 m_1 = m0 = m1 = m2 = 0;
687 for (int i = 1; i < GetCount(); ++i) {
688 double x = Membercall(getdataX)(i);
689 double Si = Membercall(getdataY)(i);
690
691 double deltaX;
692 double fi, fi_1;
693 if (frequency) {
694 fi = x;
695 fi_1 = prevx;
696 deltaX = fi - fi_1;
697 } else {
698 fi = 1/x;
699 fi_1 = 1/prevx;
700 deltaX = fi_1 - fi;
701 }
702 if (fi != 0 && fi_1 != 0) {
703 m_1 += (Si/fi + Si_1/fi_1)*deltaX;
704 m0 += (Si + Si_1)*deltaX;
705 m1 += (Si*fi + Si_1*fi_1)*deltaX;
706 m2 += (Si*fi*fi + Si_1*fi_1*fi_1)*deltaX;
707 }
708 prevx = x;
709 Si_1 = Si;
710 }
711 m_1 /= 2;
712 m0 /= 2;
713 m1 /= 2;
714 m2 /= 2;
715 }
716
SameX(DataSource & data)717 bool DataSource::SameX(DataSource &data) {
718 int64 num = GetCount();
719 if (num == 0)
720 return false;
721 if (data.GetCount() != num)
722 return false;
723 for (int64 i = 0; i < num; ++i) {
724 if (data.x(i) != x(i))
725 return false;
726 }
727 return true;
728 }
729
SortData(Getdatafun getdata)730 Vector<double> DataSource::SortData(Getdatafun getdata) {
731 Vector<double> ret;
732 int count = 0;
733 for (int64 i = 0; i < GetCount(); ++i) {
734 double d = Membercall(getdata)(i);
735 if (!IsNull(d)) {
736 ret << d;
737 count++;
738 }
739 }
740 if (count == 0)
741 return ret;
742
743 Sort(ret);
744
745 return ret;
746 }
747
Percentile(Getdatafun getdata,double rate)748 Vector<double> DataSource::Percentile(Getdatafun getdata, double rate) {
749 ASSERT(rate >= 0 && rate <= 1);
750 Vector<double> data = SortData(getdata);
751 int num = int(data.GetCount()*rate) + 1;
752 if (num < data.GetCount())
753 data.Remove(num, data.GetCount()-num);
754 return data;
755 }
756
PercentileVal(Getdatafun getdata,double rate)757 double DataSource::PercentileVal(Getdatafun getdata, double rate) {
758 ASSERT(rate >= 0 && rate <= 1);
759
760 Vector<double> data = SortData(getdata);
761 int num = int(data.GetCount()*rate);
762 return LinearInterpolate<double>(data.GetCount()*rate, num, num+1, data[num-1], data[num]);
763 }
764
Derivative(Getdatafun getdataY,Getdatafun getdataX,int orderDer,int orderAcc)765 Vector<Pointf> DataSource::Derivative(Getdatafun getdataY, Getdatafun getdataX, int orderDer, int orderAcc) {
766 ASSERT(orderDer >= 1 && orderDer <= 2);
767 ASSERT(orderAcc == 2 || orderAcc == 4 || orderAcc == 6 || orderAcc == 8);
768
769 int numData = int(GetCount());
770
771 Vector<Pointf> data;
772 data.SetCount(numData);
773 int num = 0;
774 for (int i = 0; i < numData; ++i) {
775 double y = Membercall(getdataY)(i);
776 double x = Membercall(getdataX)(i);
777 if (!IsNull(y) && !IsNull(x)) {
778 data[i].y = y;
779 data[i].x = x;
780 num++;
781 }
782 }
783 numData = num;
784
785 Vector<Pointf> res;
786 if (numData < orderAcc+1)
787 return res;
788 data.SetCount(numData);
789
790 double minD = -DOUBLE_NULL_LIM, maxD = DOUBLE_NULL_LIM;
791 for (int i = 0; i < numData-1; ++i) {
792 double d = data[i+1].x - data[i].x;
793 minD = min(minD, d);
794 maxD = max(maxD, d);
795 }
796 if ((maxD - minD)/minD > 0.0001)
797 return res;
798
799 double h = (minD + maxD)/2;
800
801 VectorXd origE(numData);
802 for (int i = 0; i < numData; ++i)
803 origE(i) = data[i].y;
804
805 // From https://en.wikipedia.org/wiki/Finite_difference_coefficient
806 double kernels1[4][9] = {{-1/2., 0, 1/2.},
807 {1/12., -2/3., 0, 2/3., -1/12.},
808 {-1/60., 3/20., -3/4., 0, 3/4., -3/20., 1/60.},
809 {1/280., -4/105., 1/5., -4/5., 0., 4/5., -1/5., 4/105., -1/280.}};
810 double kernels2[4][9] = {{1, -2, 1},
811 {-1/12., 4/3., -5/2., 4/3., -1/12.},
812 {1/90., -3/20., 3/2., -49/18., 3/2., -3/20., 1/90.},
813 {-1/560., 8/315., -1/5., 8/5., -205/72., 8/5., -1/5., 8/315., -1/560.}};
814
815 int idkernel = orderAcc/2-1;
816
817 double factor;
818 VectorXd kernel;
819 if (orderDer == 1) {
820 factor = 1/h;
821 kernel = Map<MatrixXd>(kernels1[idkernel], orderAcc+1, 1);
822 } else if (orderDer == 2) {
823 factor = 1/h/h;
824 kernel = Map<MatrixXd>(kernels2[idkernel], orderAcc+1, 1);
825 } else
826 return res;
827
828 VectorXd resE = Convolution(origE, kernel, factor);
829
830 res.SetCount(numData-orderAcc);
831 int frame = orderAcc/2 + 1;
832 for (int i = 0; i < numData-orderAcc; ++i) {
833 res[i].y = resE(i);
834 res[i].x = data[i+frame].x;
835 }
836 return res;
837 }
838
SavitzkyGolay_CheckParams(int nleft,int nright,int deg,int der)839 bool SavitzkyGolay_CheckParams(int nleft, int nright, int deg, int der) {
840 return nleft >= 0 && nright >= 0 && der <= deg && nleft + nright >= deg;
841 }
842
SavitzkyGolay_Coeff(int nleft,int nright,int deg,int der)843 VectorXd SavitzkyGolay_Coeff(int nleft, int nright, int deg, int der) {
844 ASSERT(SavitzkyGolay_CheckParams(nleft, nright, deg, der));
845
846 int cols = deg + 1;
847 MatrixXd A(cols, cols);
848
849 for(int ipj = 0; ipj <= (deg << 1); ipj++) {
850 double sum = ipj ? 0 : 1;
851 for(int k = 1; k <= nright; k++)
852 sum += pow(k, ipj);
853 for(int k = 1; k <= nleft; k++)
854 sum += pow(-k, ipj);
855 int mm = min(ipj, 2 * deg - ipj);
856 for(int imj = -mm; imj <= mm; imj += 2)
857 A((ipj+imj)/2, (ipj-imj)/2) = sum;
858 }
859
860 MatrixXd At = A.transpose();
861 VectorXd B(VectorXd::Zero(cols));
862 B[der] = 1;
863 VectorXd y = (At * A).inverse() * (At * B);
864
865 VectorXd coeff(nleft + nright + 1);
866 double factor = der > 0 ? pow(2, der-1) : 1;
867 coeff.setZero();
868 int ic = 0;
869 for(int k = -nleft; k <= nright; k++) {
870 double sum = y[0];
871 double fac = 1;
872 for(int mm = 1; mm <= deg; mm++)
873 sum += y[mm]*(fac *= k);
874 coeff[ic++] = sum*factor;
875 }
876 return coeff;
877 }
878
SavitzkyGolay_Check(const VectorXd & coeff)879 bool SavitzkyGolay_Check(const VectorXd &coeff) {
880 double unity = coeff.sum();
881 return abs(1-unity) < 0.0000001 || abs(unity) < 0.0000001;
882 }
883
SavitzkyGolay(Getdatafun getdataY,Getdatafun getdataX,int deg,int size,int der)884 Vector<Pointf> DataSource::SavitzkyGolay(Getdatafun getdataY, Getdatafun getdataX, int deg, int size, int der) {
885 int numData = int(GetCount());
886
887 Vector<Pointf> data;
888 data.SetCount(numData);
889 int num = 0;
890 for (int i = 0; i < numData; ++i) {
891 double y = Membercall(getdataY)(i);
892 double x = Membercall(getdataX)(i);
893 if (!IsNull(y) && !IsNull(x)) {
894 data[i].y = y;
895 data[i].x = x;
896 num++;
897 }
898 }
899 numData = num;
900
901 Vector<Pointf> res;
902 if (numData < size)
903 return res;
904 data.SetCount(numData);
905
906 double minD = -DOUBLE_NULL_LIM, maxD = DOUBLE_NULL_LIM;
907 for (int i = 0; i < numData-1; ++i) {
908 double d = data[i+1].x - data[i].x;
909 minD = min(minD, d);
910 maxD = max(maxD, d);
911 }
912 if ((maxD - minD)/minD > 0.0001)
913 return res;
914
915 double h = (minD + maxD)/2;
916
917 VectorXd coeff = SavitzkyGolay_Coeff(size/2, size/2, deg, der);
918 if (!SavitzkyGolay_Check(coeff))
919 return res;
920
921 VectorXd origE(numData);
922 for (int i = 0; i < numData; ++i)
923 origE[i] = data[i].y;
924
925 VectorXd resE = Convolution(origE, coeff, 1/pow(h, der));
926
927 res.SetCount(numData-size);
928 int frame = size/2;
929 for (int i = 0; i < numData-size; ++i) {
930 res[i].y = resE(i);
931 res[i].x = data[i+frame].x;
932 }
933 return res;
934 }
935
936
Init(Function<double (double x,double y)> _funz,double _minX,double _maxX,double _minY,double _maxY)937 void ExplicitData::Init(Function<double (double x, double y)> _funz, double _minX, double _maxX, double _minY, double _maxY) {
938 ASSERT(maxX >= minX && maxY >= minY);
939 this->funz = _funz;
940 this->minX = _minX;
941 this->maxX = _maxX;
942 this->minY = _minY;
943 this->maxY = _maxY;
944
945 minZ = -DOUBLE_NULL_LIM;
946 maxZ = DOUBLE_NULL_LIM;
947 double deltax = (_maxX - _minX)/100.;
948 double deltay = (_maxY - _minY)/100.;
949 for (double x = _minX; x <= _maxX; x += deltax) {
950 for (double y = _minY; y <= _maxY; y += deltay) {
951 double z = _funz(x, y);
952 if (!IsNull(z)) {
953 minZ = min(minZ, z);
954 maxZ = max(maxZ, z);
955 }
956 }
957 }
958 }
959
FindClosest(Pointf & p,Vector<Pointf> & points,double deltaX,double deltaY,double & d)960 int FindClosest(Pointf &p, Vector<Pointf> &points, double deltaX, double deltaY, double &d) {
961 double dxmin = -DOUBLE_NULL, dymin = -DOUBLE_NULL;
962 int imin = -1;
963 for (int i = 0; i < points.GetCount(); ++i) {
964 double dx = abs(p.x - points[i].x);
965 double dy = abs(p.y - points[i].y);
966 if (dx*dx + dy*dy < dxmin*dxmin + dymin*dymin) {
967 dxmin = dx;
968 dymin = dy;
969 imin = i;
970 }
971 }
972 if ((dxmin > deltaX*1.00000000001) || (dymin > deltaY*1.00000000001))
973 return Null;
974 d = sqrt(dxmin*dxmin + dymin*dymin);
975 return imin;
976 }
977
GetIsolines(const Vector<double> & vals,const Rectf & area,double deltaX,double deltaY)978 Vector<Pointf> DataSourceSurf::GetIsolines(const Vector<double> &vals, const Rectf &area, double deltaX, double deltaY) {
979 Vector<Pointf> isolines;
980 for (int i = 0; i < vals.GetCount(); ++i) {
981 if (i > 0)
982 isolines << Null;
983 Vector<Pointf> isoaux = GetIsoline(vals[i], area, deltaX, deltaY);
984 isolines.Append(isoaux);
985 }
986 return isolines;
987 }
988
GetIsoline(double thres,const Rectf & area,double deltaX,double deltaY)989 Vector<Pointf> DataSourceSurf::GetIsoline(double thres, const Rectf &area, double deltaX, double deltaY) {
990 Vector<double> zp;
991
992 int width = (int)(area.GetWidth()/deltaX) + 1;
993 int height = -(int)(area.GetHeight()/deltaY) + 1;
994 zp.SetCount(width*height);
995 int iy = 0;
996 for (double y = area.bottom; iy < height; y += deltaY, iy++) {
997 int ix = 0;
998 for (double x = area.left; ix < width; x += deltaX, ix++)
999 zp[ix + iy*width] = z(x, y);
1000 }
1001
1002 Vector<Pointf> points;
1003 for (int iy = 0; iy < height; iy++) {
1004 for (int ix = 0; ix < width-1; ix++) {
1005 double z0 = zp[ix + iy*width];
1006 double z1 = zp[ix+1 + iy*width];
1007 if (IsNull(z0) || IsNull(z1))
1008 continue;
1009 if ((z1 > thres && z0 <= thres) || (z0 > thres && z1 <= thres)) {
1010 double delta = abs(thres - z0)/abs(z1 - z0);
1011 points << Pointf(area.left + (ix + delta)*deltaX, area.bottom + iy*deltaY);
1012 }
1013 }
1014 }
1015 for (int ix = 0; ix < width; ix++) {
1016 for (int iy = 0; iy < height-1; iy++) {
1017 double z0 = zp[ix + iy*width];
1018 double z1 = zp[ix + (iy+1)*width];
1019 if (IsNull(z0) || IsNull(z1))
1020 continue;
1021 if ((z1 > thres && z0 <= thres) || (z0 > thres && z1 <= thres)) {
1022 double delta = abs(thres - z0)/abs(z1 - z0);
1023 points << Pointf(area.left + ix*deltaX, area.bottom + (iy + delta)*deltaY);
1024 }
1025 }
1026 }
1027 if (points.IsEmpty())
1028 return points;
1029
1030 Vector<Pointf> isoline;
1031 isoline << points[0];
1032 points.Remove(0);
1033 while (!points.IsEmpty()) {
1034 int imin;
1035 double dt, d0;
1036 int iminT = FindClosest(isoline.Top(), points, deltaX, deltaY, dt);
1037 int imin0 = FindClosest(isoline[0], points, deltaX, deltaY, d0);
1038 if (IsNull(iminT) && IsNull(imin0)) {
1039 isoline << Null;
1040 imin = 0;
1041 } else if (IsNull(iminT)) {
1042 Reverse(isoline);
1043 imin = imin0;
1044 } else if (IsNull(imin0))
1045 imin = iminT;
1046 else {
1047 if (dt > d0) {
1048 Reverse(isoline);
1049 imin = imin0;
1050 } else
1051 imin = iminT;
1052 }
1053 isoline << points[imin];
1054 points.Remove(imin);
1055 }
1056 return isoline;
1057 }
1058
Intersection(Pointf & a,Pointf & b,Pointf & c,Pointf & d)1059 Pointf Intersection(Pointf &a, Pointf &b, Pointf &c, Pointf &d) {
1060 // Line AB represented as a1x + b1y = c1
1061 double a1 = b.y - a.y;
1062 double b1 = a.x - b.x;
1063 double c1 = a1*(a.x) + b1*(a.y);
1064
1065 // Line CD represented as a2x + b2y = c2
1066 double a2 = d.y - c.y;
1067 double b2 = c.x - d.x;
1068 double c2 = a2*(c.x)+ b2*(c.y);
1069
1070 double det = a1*b2 - a2*b1;
1071
1072 if (det == 0) // Parallel
1073 return Null;
1074 else {
1075 double x = (b2*c1 - b1*c2)/det;
1076 double y = (a1*c2 - a2*c1)/det;
1077 return Pointf(x, y);
1078 }
1079 }
1080
SegIntersection(Pointf & a,Pointf & b,Pointf & c,Pointf & d)1081 Pointf SegIntersection(Pointf &a, Pointf &b, Pointf &c, Pointf &d) {
1082 Pointf inter = Intersection(a, b, c, d);
1083 if (IsNull(inter))
1084 return Null;
1085 if (((a.x <= inter.x && b.x >= inter.x) || (b.x <= inter.x && a.x >= inter.x)) &&
1086 ((a.y <= inter.y && b.y >= inter.y) || (b.y <= inter.y && a.y >= inter.y)) &&
1087 ((c.x <= inter.x && d.x >= inter.x) || (d.x <= inter.x && c.x >= inter.x)) &&
1088 ((c.y <= inter.y && d.y >= inter.y) || (d.y <= inter.y && c.y >= inter.y)))
1089 return inter;
1090 else
1091 return Null;
1092 }
1093
Intersection(Vector<Pointf> & poly1,Vector<Pointf> & poly2)1094 Vector<Pointf> Intersection(Vector<Pointf> &poly1, Vector<Pointf> &poly2) {
1095 Vector<Pointf> listInter;
1096 for (int i1 = 0; i1 < poly1.GetCount() - 1; ++i1) {
1097 for (int i2 = 0; i2 < poly2.GetCount() - 1; ++i2) {
1098 Pointf inter = SegIntersection(poly1[i1], poly1[i1+1], poly2[i2], poly2[i2+1]);
1099 if (!IsNull(inter)) {
1100 bool found = false;
1101 for (int i = 0; i < listInter.GetCount(); ++i) {
1102 if (abs((inter.x - listInter[i].x)/inter.x) < 0.000001 &&
1103 abs((inter.y - listInter[i].y)/inter.y) < 0.000001) {
1104 found = true;
1105 break;
1106 }
1107 }
1108 if (!found)
1109 listInter << inter;
1110 }
1111 }
1112 }
1113 return listInter;
1114 }
1115
Simplify(Vector<Pointf> & poly,double dx,double dy)1116 void Simplify(Vector<Pointf> &poly, double dx, double dy) {
1117 for (int i = 1; i < poly.GetCount(); ++i) {
1118 if (abs(poly[i].x - poly[i-1].x) < dx && abs(poly[i].y - poly[i-1].y) < dy) {
1119 poly.Remove(i, 1);
1120 i--;
1121 }
1122 }
1123 }
1124
1125 }
1126