WSL/SLF GitLab Repository

libinterpol1D.cc 21.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
18
#include <meteoio/meteostats/libinterpol1D.h>
19
#include <meteoio/MathOptim.h>
20
#include <algorithm>
21
#include <cmath>
22
23
24
25
26

using namespace std;

namespace mio {

27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
/**
 * @brief This function returns a vector of quantiles.
 * The vector does not have to be sorted. See https://secure.wikimedia.org/wikipedia/en/wiki/Quartile for more.
 * This code is heavily inspired by Ken Wilder, https://sites.google.com/site/jivsoft/Home/compute-ranks-of-elements-in-a-c---array-or-vector
 * (quantile method, replacing the nth-element call by direct access to a sorted vector).
 * @param X vector to classify
 * @param quartiles vector of quartiles, between 0 and 1
 * @return vector of ordinates of the quantiles
 */
std::vector<double> Interpol1D::quantiles(const std::vector<double>& X, const std::vector<double>& quartiles)
{
	const size_t Xsize = X.size();
	const size_t Qsize = quartiles.size();
	if (Xsize == 0)
		throw NoAvailableDataException("Trying to calculate quantiles with no data points", AT);
	if (Qsize == 0)
		throw NoAvailableDataException("No quantiles specified", AT);

	//in order to properly escape nodata points, we need to copy in a temporary vector
	vector<double> vecTemp;
	for(size_t i=0; i<Xsize; i++) {
		const double& value=X[i];
		if(value!=IOUtils::nodata)
			vecTemp.push_back(value);
	}

	//we will store results in a new vector
	std::vector<double> vecResults(Qsize, IOUtils::nodata);

	const size_t vecSize = vecTemp.size();
	if (vecSize == 0) {
		return vecResults; //ie: nodata values
	}
	if(vecSize == 1) {
		std::fill(vecResults.begin(), vecResults.end(), vecTemp[0]);
		return vecResults;
	}

	//compute quantiles
66
	std::sort( vecTemp.begin(), vecTemp.end()); //since we will process several values, we sort the vector
67
68
	for(size_t ii=0; ii<Qsize; ii++) {
		const double q = quartiles[ii];
69
70
		if(q<=0.) vecResults[ii] = vecTemp.front();
		else if(q>=1.) vecResults[ii] = vecTemp.back();
71
72
		else {
			const double pos = static_cast<double>(vecSize - 1) * q;
73
			const size_t ind = static_cast<size_t>(pos);
74
75
76
77
78
79
80
81
82
			const double delta = pos - static_cast<double>(ind);
			const double i1 = vecTemp[ind];
			const double i2 = vecTemp[ind+1];
			vecResults[ii] = i1 * (1.0 - delta) + i2 * delta;
		}
	}
	return vecResults;
}

83
//small helper function for checking if a point can be used when computing a vector derivative
84
inline bool Interpol1D::ptOK(const double& x, const double& y) {
85
86
87
	return (x!=IOUtils::nodata && y!=IOUtils::nodata);
}

88
89
90
91
92
93
94
95
96
97
/**
 * @brief This function returns the vector of local derivatives, given a vector of abscissae and ordinates.
 * The vectors must be sorted by ascending x. The derivatives will be centered if possible, left or right otherwise or nodata
 * if nothing else can be computed.
 * @param X vector of abscissae
 * @param Y vector of ordinates
 * @return vector of local derivatives
 */
std::vector<double> Interpol1D::derivative(const std::vector<double>& X, const std::vector<double>& Y)
{
98
99
	const size_t n = X.size();
	if(n!=Y.size()) {
100
		ostringstream ss;
101
102
103
104
		ss << "X vector and Y vector don't match! " << n << "!=" << Y.size() << "\n";
		throw InvalidArgumentException(ss.str(), AT);
	}
	if(n<2) {
105
		ostringstream ss;
106
		ss << "X and Y vector only contain " << n << "points, it is not possible to compute a derivative!\n";
107
108
109
		throw InvalidArgumentException(ss.str(), AT);
	}

110
	std::vector<double> der(n, IOUtils::nodata);
111
	//right hand derivative
112
113
114
115
116
117
118
	{
		const double x=X[0], x1=X[1];
		const double y=Y[0], y1=Y[1];
		const double Dx_r=x1-x;
		if(ptOK(x,y) && ptOK(x1,y1) && Dx_r!=0.)
			der[0] = (y1-y) / Dx_r;
	}
119
120

	//centered derivative if possible
121
122
123
124
	for(size_t i=1; i<(n-1); i++) {
		const double x0=X[i-1], x=X[i], x1=X[i+1];
		const double y0=Y[i-1], y=Y[i], y1=Y[i+1];
		const double Dx_r=x1-x, Dx_c=x1-x0, Dx_l=x-x0;
125

126
127
128
		const double right = (ptOK(x,y) && ptOK(x1,y1) && Dx_r!=0.)? (y1-y) / Dx_r : IOUtils::nodata;
		const double left = (ptOK(x,y) && ptOK(x0,y0) && Dx_l!=0.)? (y-y0) / Dx_l : IOUtils::nodata;
		const double centered = (ptOK(x0,y0) && ptOK(x1,y1) && Dx_c!=0.)? (y1-y0) / Dx_c : IOUtils::nodata;
129

130
131
132
		if(centered!=IOUtils::nodata) der[i] = centered;
		else if(right!=IOUtils::nodata) der[i] = right;
		else if(left!=IOUtils::nodata) der[i] = left;
133
134
135
	}

	//left hand derivative
136
137
138
139
140
141
142
143
	{
		const size_t last = n-1;
		const double x0=X[last-1], x=X[last];
		const double y0=Y[last-1], y=Y[last];
		const double Dx_l=x-x0;
		if(ptOK(x,y) && ptOK(x0,y0) && Dx_l!=0.)
			der[last] = (y-y0) / Dx_l;
	}
144
145
146
147
148
149

	return der;
}

void Interpol1D::sort(std::vector<double>& X, std::vector<double>& Y)
{
150
151
	const size_t Xsize = X.size();
	if(Xsize!=Y.size()) {
152
		ostringstream ss;
153
		ss << "X vector and Y vector don't match! " << Xsize << "!=" << Y.size() << "\n";
154
155
156
		throw InvalidArgumentException(ss.str(), AT);
	}

157
158
	std::vector< std::pair<double,double> > new_vec( Xsize );
	for(size_t i=0; i<Xsize; i++) {
159
		const std::pair<double,double> tmp(X[i],Y[i]);
160
		new_vec[i] = tmp;
161
162
163
164
	}

	std::sort( new_vec.begin(), new_vec.end(), pair_comparator );

165
	for(size_t i=0; i<Xsize; i++) {
166
167
168
169
170
		X[i] = new_vec[i].first;
		Y[i] = new_vec[i].second;
	}
}

171
inline bool Interpol1D::pair_comparator(const std::pair<double, double>& l, const std::pair<double, double>& r) {
172
173
174
	return l.first < r.first;
}

175
176
177
178
179
180
181
182
183
184
185
/**
 * @brief This function returns the weighted aritmetic mean of two numbers.
 * A weight of 0 returns d1, a weight of 1 returns d2, a weight of 0.5 returns a centered mean.
 * See https://secure.wikimedia.org/wikipedia/en/wiki/Weighted_mean for more...
 * @param d1 first value
 * @param d2 second value
 * @param weight weight to apply to the mean
 * @return weighted aritmetic mean
 */
double Interpol1D::weightedMean(const double& d1, const double& d2, const double& weight)
{
186
	const double tmp = abs(d1 - d2);
187
188
189
190
191
192
193
194
	if (d1 < d2) {
		return (d1 + tmp*weight);
	} else {
		return (d1 - tmp*weight);
	}
}

/**
195
 * @brief This function returns the weighted aritmetic mean of a vector.
196
197
198
199
200
201
202
203
204
205
206
 * See https://secure.wikimedia.org/wikipedia/en/wiki/Weighted_mean for more...
 * @param vecData vector of values
 * @param weight weights to apply to the mean
 * @return weighted aritmetic mean
 */
double Interpol1D::weightedMean(const std::vector<double>& vecData, const std::vector<double>& weight)
{
	const size_t nPts = vecData.size();
	if (nPts == 0)
		throw NoAvailableDataException("Trying to calculate an arithmetic mean with no data points", AT);
	if(nPts != weight.size()) {
207
		std::ostringstream ss;
208
209
210
211
212
213
214
		ss << "Computing weighted mean of a vector of size " << nPts;
		ss << " with vector of weights of size " << weight.size();
		throw InvalidArgumentException(ss.str(), AT);
	}

	double sum = 0., count = 0.;
	for (size_t ii=0; ii<nPts; ii++){
215
		const double value = vecData[ii];
216
217
		if(value!=IOUtils::nodata) {
			const double w = weight[ii];
218
			sum += value*w;
219
220
221
222
223
224
225
226
227
			count += w;
		}
	}

	if(count>0.)
		return (sum/count);
	else
		return IOUtils::nodata;
}
228
229

double Interpol1D::arithmeticMean(const std::vector<double>& vecData)
230
{
231
232
	const size_t nPts = vecData.size();
	if (nPts == 0)
233
234
		throw NoAvailableDataException("Trying to calculate an arithmetic mean with no data points", AT);

235
	unsigned int count=0;
236
	double sum = 0.0;
237
	for (size_t ii=0; ii<nPts; ii++){
238
		const double value = vecData[ii];
239
		if(value!=IOUtils::nodata) {
240
			sum += value;
241
242
			count++;
		}
243
244
	}

245
246
247
248
	if(count>0)
		return (sum/(double)count);
	else
		return IOUtils::nodata;
249
250
251
}

double Interpol1D::getMedian(const std::vector<double>& vecData)
252
{
253
//This uses a sorting algorithm for getting middle element
254
// as much more efficient than full sorting (O(n) compared to O(n log(n))
255
	if (vecData.empty())
256
		throw NoAvailableDataException("Trying to calculate a median with no data points", AT);
257

258
	vector<double> vecTemp;
259
	for(size_t i=0; i<vecData.size(); i++) {
260
261
262
263
		const double& value=vecData[i];
		if(value!=IOUtils::nodata)
			vecTemp.push_back(value);
	}
264

265
266
	const size_t vecSize = vecTemp.size();
	if (vecSize == 0)
267
		return IOUtils::nodata;
268

269
	if ((vecSize % 2) == 1){ //uneven
270
271
		const int middle = (int)(vecSize/2);
		nth_element(vecTemp.begin(), vecTemp.begin()+middle, vecTemp.end());
272
		return *(vecTemp.begin()+middle);
273
	} else { //use arithmetic mean of element n/2 and n/2-1
274
275
276
277
278
279
		const int middle = (int)(vecSize/2);
		nth_element(vecTemp.begin(), vecTemp.begin()+middle-1, vecTemp.end());
		const double m1 = *(vecTemp.begin()+middle-1);
		nth_element(vecTemp.begin(), vecTemp.begin()+middle, vecTemp.end());
		const double m2 = *(vecTemp.begin()+middle);
		return weightedMean( m1, m2, 0.5);
280
281
282
283
	}
}

double Interpol1D::getMedianAverageDeviation(const std::vector<double>& vecData)
284
{
285
	if (vecData.empty())
286
		throw NoAvailableDataException("Trying to calculate MAD with no data points", AT);
287

288
289
	vector<double> vecWindow(vecData);

290
	const double median = Interpol1D::getMedian(vecWindow);
291
292
	if(median==IOUtils::nodata)
		return IOUtils::nodata;
293
294

	//Calculate vector of deviations and write each value back into the vecWindow
295
	for(size_t ii=0; ii<vecWindow.size(); ii++){
296
		double& value = vecWindow[ii];
297
298
		if(value!=IOUtils::nodata)
			value = std::abs(value - median);
299
300
301
	}

	//Calculate the median of the deviations
302
	const double mad = Interpol1D::getMedian(vecWindow);
303
304
305
306

	return mad;
}

307
double Interpol1D::variance(const std::vector<double>& X)
308
309
310
{//The variance is computed using a compensated variance algorithm,
//(see https://secure.wikimedia.org/wikipedia/en/wiki/Algorithms_for_calculating_variance)
//in order to be more robust to small variations around the mean.
311
312
	const size_t n = X.size();
	size_t count=0;
313
	double sum=0.;
314

315
	for(size_t i=0; i<n; i++) {
316
		const double value = X[i];
317
318
319
320
321
322
323
		if(value!=IOUtils::nodata) {
			sum += value;
			count++;
		}
	}

	if(count<=1) return IOUtils::nodata;
324

325
	const double mean = sum/(double)count;
326
	double sum2=0., sum3=0.;
327
	for(size_t i=0; i<n; i++) {
328
		const double value = X[i];
329
		if(value!=IOUtils::nodata) {
330
331
332
			const double delta = value - mean;
			sum2 += delta*delta;
			sum3 += delta;
333
334
		}
	}
335
	const double variance = (sum2 - sum3*sum3/static_cast<double>(count)) / static_cast<double>(count - 1);
336
	return variance;
337
338
}

339
double Interpol1D::std_dev(const std::vector<double>& X) {
340
341
342
343
	return sqrt(variance(X));
}

double Interpol1D::covariance(const std::vector<double>& X, const std::vector<double>& Y)
344
{//this is a simple but still compensated covariance computation (see the notes on the variance)
345
346
	const size_t Xsize = X.size();
	if(Xsize!=Y.size())
347
		throw IOException("Vectors should have the same size for covariance!", AT);
348
	if(Xsize==0) return IOUtils::nodata;
349
350
351
352
353
354

	const double X_mean = Interpol1D::arithmeticMean(X);
	const double Y_mean = Interpol1D::arithmeticMean(Y);
	if(X_mean==IOUtils::nodata || Y_mean==IOUtils::nodata)
		return IOUtils::nodata;

355
	size_t count=0;
356
	double sum=0.;
357
	for(size_t i=0; i<Xsize; i++) {
358
359
360
361
362
363
364
365
366
		if(X[i]!=IOUtils::nodata && Y[i]!=IOUtils::nodata) {
			sum += (X[i] - X_mean) * (Y[i] - Y_mean);
			count++;
		}
	}
	if(count<=1) return IOUtils::nodata;
	return sum/((double)count-1.);
}

367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
/**
* @brief Computes the distance between a point (x,y) and a line y=ax+b
* @param x x coordinates of the point
* @param y y coordinates of the point
* @param a slope of the line
* @param c origin of the line
* @return distance of the point to the line
*/
double Interpol1D::pt_line_distance(const double& x, const double& y, const double& a, const double& c) {
	if(a==0.) return abs(y-c); //horizontal line

	//for ax+by+c=0; for us, b=-1
	const double b = -1.;
	const double d = abs(a*x +b*y + c) * Optim::invSqrt( a*a + b*b );
	return d;
}

384
385
386
387
388
389
390
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* the linear regression has the form Y = aX + b with a regression coefficient r (it is nodata safe)
* @param X vector of X coordinates
* @param Y vector of Y coordinates (same order as X)
* @param a slope of the linear regression
* @param b origin of the linear regression
391
* @param r absolute value of linear regression coefficient
392
* @param mesg information message if something fishy is detected
393
* @param fixed_rate force the lapse rate? (default=false)
394
*/
395
void Interpol1D::LinRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::ostringstream& mesg, const bool& fixed_rate)
396
397
398
399
400
401
402
{
	if(fixed_rate) {
		LinRegressionFixedRate(X, Y, a, b, r, mesg);
		return;
	}

	//check arguments
403
	const size_t n = X.size();
404
405
	if(n<2)
		throw NoAvailableDataException("Trying to calculate linear regression with too few data points", AT);
406
407
408
409
410
	if(n!=Y.size())
		throw IOException("Vectors should have the same size for linear regression!", AT);

	//computing x_avg and y_avg
	double x_avg=0., y_avg=0.;
411
	size_t count=0;
412
413
414
415
	for (size_t ii=0; ii<n; ii++) {
		if(X[ii]!=IOUtils::nodata && Y[ii]!=IOUtils::nodata) {
			x_avg += X[ii];
			y_avg += Y[ii];
416
417
418
			count++;
		}
	}
419
420
	if(count<2)
		throw NoAvailableDataException("Trying to calculate linear regression with too few valid data points", AT);
421
422
423
424
425
	x_avg /= (double)count;
	y_avg /= (double)count;

	//computing sx, sy, sxy
	double sx=0., sy=0., sxy=0.;
426
427
428
429
430
	for (size_t ii=0; ii<n; ii++) {
		if(X[ii]!=IOUtils::nodata && Y[ii]!=IOUtils::nodata) {
			sx += (X[ii]-x_avg) * (X[ii]-x_avg);
			sy += (Y[ii]-y_avg) * (Y[ii]-y_avg);
			sxy += (X[ii]-x_avg) * (Y[ii]-y_avg);
431
432
433
434
435
436
437
438
439
440
		}
	}

	//computing the regression line
	const double epsilon = 1e-6;
	if(sx <= abs(x_avg)*epsilon) { //sx and sy are always positive
		//all points have same X -> we return a constant value that is the average
		a = 0.;
		b = y_avg;
		r = 1.;
441
		mesg << "[W] Computing linear regression on data at identical X\n";
442
443
444
445
446
447
448
449
450
		return;
	}
	a = sxy / sx;
	b = y_avg - a*x_avg;
	if(sy==0) {
		//horizontal line: all y's are equals
		r = 1.;
	} else {
		//any other line
451
		r = abs( sxy / sqrt(sx*sy) );
452
453
454
	}
}

455
456
457
458
459
460
461
462
463
464
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* the linear regression has the form Y = aX + b with a regression coefficient r (it is nodata safe) while forcing the value of a
* @param X vector of X coordinates
* @param Y vector of Y coordinates (same order as X)
* @param a slope of the linear regression (forced)
* @param b origin of the linear regression
* @param r absolute value of linear regression coefficient
* @param mesg information message if something fishy is detected
*/
465
void Interpol1D::LinRegressionFixedRate(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::ostringstream& mesg)
466
467
468
469
470
471
472
473
474
475
{	//check arguments
	const size_t n = X.size();
	if(n==0)
		throw NoAvailableDataException("Trying to calculate linear regression with no data points", AT);
	if(n!=Y.size())
		throw IOException("Vectors should have the same size for linear regression!", AT);

	//computing x_avg and y_avg
	int count=0;
	double x_avg=0, y_avg=0;
476
477
478
479
	for (size_t ii=0; ii<n; ii++) {
		if(X[ii]!=IOUtils::nodata && Y[ii]!=IOUtils::nodata) {
			x_avg += X[ii];
			y_avg += Y[ii];
480
481
482
483
484
485
486
487
488
489
490
491
			count++;
		}
	}
	if(count==0)
		throw NoAvailableDataException("Trying to calculate linear regression with no valid data points", AT);
	x_avg /= (double)count;
	y_avg /= (double)count;

	//computing the regression line
	b = y_avg - a*x_avg;

	double TSS=0, SSR=0; //Total Sum of Squares and Sum of Squared Residuals
492
493
494
495
	for (size_t ii=0; ii<n; ii++) {
		if(X[ii]!=IOUtils::nodata && Y[ii]!=IOUtils::nodata) {
			SSR += Optim::pow2( Y[ii] - (a*X[ii]+b) );
			TSS += Optim::pow2( Y[ii] );
496
497
498
499
500
501
502
503
504
505
		}
	}
	if(TSS!=0) {
		r = 1. - SSR/TSS;
	} else {
		r = 1.; //when all Y[i]=0 we automatically pick up the best possible fit. But r does not mean anything...
		mesg << "[W] Computing fixed lapse rate linear regression on data all at Y=0\n";
	}
}

506
507
508
509
510
511
512
513
514
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* the linear regression has the form Y = aX + b with a regression coefficient r. If the regression coefficient is not good enough, tries to remove bad points (up to 15% of the initial data set can be removed, keeping at least 4 points)
* @param in_X vector of X coordinates
* @param in_Y vector of Y coordinates (same order as X)
* @param A slope of the linear regression
* @param B origin of the linear regression
* @param R linear regression coefficient
* @param mesg information message if something fishy is detected
515
* @param fixed_rate force the lapse rate? (default=false)
516
*/
517
void Interpol1D::NoisyLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, double& A, double& B, double& R, std::ostringstream& mesg, const bool& fixed_rate)
518
519
{
	//finds the linear regression for points (x,y,z,Value)
520
	const double r_thres = 0.7;
521
	const size_t nb_pts = in_X.size();
522
523
	//we want at least 4 points AND 75% of the initial data set kept in the regression
	const size_t min_dataset = (size_t)Optim::floor( 0.75*(double)nb_pts );
524
	const size_t min_pts = (min_dataset>4)? min_dataset : 4;
525

526
	LinRegression(in_X, in_Y, A, B, R, mesg, fixed_rate);
527
	if(R>=r_thres) return;
528
529

	std::vector<double> X(in_X), Y(in_Y);
530
	size_t nb_valid_pts = nb_pts;
531

532
	while(R<r_thres && nb_valid_pts>min_pts) {
533
		//we try to remove the one point in the data set that is the worst
534
		size_t index_bad=0;
535
536
537
538
539
540
541
		double max_dist = -1.;
		for (size_t ii=0; ii<nb_pts; ii++) {
			if(Y[ii]==IOUtils::nodata) continue;
			const double dist = pt_line_distance(X[ii], Y[ii], A, B);
			if(dist>max_dist) {
				max_dist = dist;
				index_bad = ii;
542
543
			}
		}
544
		//the worst point has been found, we remove it
545
		Y[index_bad] = IOUtils::nodata;
546
		nb_valid_pts--;
547
		LinRegression(X, Y, A, B, R, mesg, fixed_rate);
548
549
550
	}

	//check if r is reasonnable
551
	if (R<r_thres) {
552
553
554
		mesg << "[W] Poor regression coefficient: " << std::setprecision(4) << R << "\n";
	}
}
555

556
557
558
559
560
561
562
563
/**
* @brief Computes the bi-linear regression coefficients fitting the points given as X and Y in two vectors
* We consider that the regression can be made with 2 linear segments with a fixed inflection point. It relies on Interpol1D::NoisyLinRegression.
* @param in_X vector of X coordinates
* @param in_Y vector of Y coordinates (same order as X)
* @param bilin_inflection inflection point absissa
* @param coeffs a,b,r coefficients in a vector
*/
564
void Interpol1D::twoLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, const double& bilin_inflection, std::vector<double>& coeffs)
565
566
567
{
	//build segments
	std::vector<double> X1, Y1, X2, Y2;
568
	for(size_t ii=0; ii<in_X.size(); ii++) {
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
		if(in_X[ii]<bilin_inflection) { //first segment
			X1.push_back( in_X[ii] );
			Y1.push_back( in_Y.at(ii) );
		} else if(in_X[ii]>bilin_inflection) { //second segment
			X2.push_back( in_X[ii] );
			Y2.push_back( in_Y.at(ii) );
		} else { //point belongs to both segments
			X1.push_back( in_X[ii] );
			Y1.push_back( in_Y.at(ii) );
			X2.push_back( in_X[ii] );
			Y2.push_back( in_Y.at(ii) );
		}
	}

	//first segment
584
	double a1, b1, r1;
585
	std::ostringstream mesg1;
586
	NoisyLinRegression(X1, Y1, a1, b1, r1, mesg1);
587
588

	//second segment
589
	double a2, b2, r2;
590
	std::ostringstream mesg2;
591
	NoisyLinRegression(X2, Y2, a2, b2, r2, mesg2);
592
593
594
595
596

	coeffs.push_back(a1); coeffs.push_back(b1);
	coeffs.push_back(a2); coeffs.push_back(b2);
}

597
598
599
600
601
602
603
604
/**
* @brief Computes the Log regression coefficients fitting the points given as X and Y in two vectors
* the log regression has the form Y = a*ln(X) + b with a regression coefficient r (it is nodata safe)
* @param X vector of X coordinates
* @param Y vector of Y coordinates (same order as X)
* @param a slope of the regression
* @param b origin of the regression
* @param r regression coefficient
605
* @param mesg information message if something fishy is detected
606
*/
607
void Interpol1D::LogRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::ostringstream& mesg)
608
{
609
	std::vector<double> x( X.size() );
610

611
	for(size_t i=0; i<X.size(); i++) {
612
		const double val = X[i];
613
		x[i] = (val!=IOUtils::nodata)? log(val) : IOUtils::nodata;
614
615
	}

616
	LinRegression(x, Y, a, b, r, mesg); //HACK: how should we transform r?
617
618
619
620
621
622
623
624
625
626
}

/**
* @brief Computes the power regression coefficients fitting the points given as X and Y in two vectors
* the power regression has the form Y = b*X^a with a regression coefficient r (it is nodata safe)
* @param X vector of X coordinates
* @param Y vector of Y coordinates (same order as X)
* @param a slope of the regression
* @param b origin of the regression
* @param r regression coefficient
627
* @param mesg information message if something fishy is detected
628
*/
629
void Interpol1D::ExpRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::ostringstream& mesg)
630
{
631
	std::vector<double> y( Y.size() );
632

633
	for(size_t i=0; i<Y.size(); i++) {
634
		const double val = Y[i];
635
		y[i] = (val!=IOUtils::nodata)? log(val) : IOUtils::nodata;
636
637
	}

638
	LinRegression(X, y, a, b, r, mesg); //HACK: how should we transform r?
639
640
}

641
} //namespace