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