WSL/SLF GitLab Repository

libfit1D.h 8.57 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/***********************************************************************************/
/*  Copyright 2011 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/>.
*/
#ifndef __LIBFIT1D_H__
#define __LIBFIT1D_H__

#include <meteoio/IOExceptions.h>
22
#include <meteoio/meteostats/libfit1DCore.h>
23
24
25
26
27

#include <vector>

namespace mio {

28
29
30
31
32
class Zero : public FitModel {
	public:
		Zero() {fit_ready = true; nParam = 0; min_nb_pts = 0; regname = "Zero";};
		void setData(const std::vector<double>& /*in_X*/, const std::vector<double>& /*in_Y*/) { };
		bool fit() { return true;};
33
		double f(const double& /*x*/) const {return 0.;};
34
35
};

36
37
class SimpleLinear : public FitModel {
	public:
38
		SimpleLinear() : fixed_lapse_rate(IOUtils::nodata) {fit_ready = false; nParam = 2; min_nb_pts = 2; regname = "SimpleLinear";};
39
		void setData(const std::vector<double>& in_X, const std::vector<double>& in_Y);
40
		bool fit();
41
		double f(const double& x) const;
42
43
44
45
46
47
		void setLapseRate(const double& in_lapse_rate) {fixed_lapse_rate = in_lapse_rate; fit_ready = false; min_nb_pts=1;};
	protected:
		bool checkInputs();
		double fixed_lapse_rate;
};

48
49
class NoisyLinear : public SimpleLinear {
	public:
50
		NoisyLinear() {fit_ready = false; nParam = 2; min_nb_pts = 2; regname = "NoisyLinear";};
51
		bool fit();
52
53
};

54
55
56
57
class SphericVario : public FitLeastSquare {
	public:
		SphericVario() {fit_ready = false; nParam = 3; min_nb_pts = 4; regname = "SphericVario";};
		void setDefaultGuess();
58
		double f(const double& x) const;
59
};
60

61
class LinVario : public FitLeastSquare {
62
	public:
63
64
		LinVario() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "LinVario";};
		void setDefaultGuess();
65
		double f(const double& x) const;
66
};
67

68
69
70
71
class ExpVario : public FitLeastSquare {
	public:
		ExpVario() {fit_ready = false; nParam = 3; min_nb_pts = 4; regname = "ExpVario";};
		void setDefaultGuess();
72
		double f(const double& x) const;
73
74
75
76
77
78
};

class RatQuadVario : public FitLeastSquare {
	public:
		RatQuadVario() {fit_ready = false; nParam = 3; min_nb_pts = 4; regname = "RatQuadVario";};
		void setDefaultGuess();
79
		double f(const double& x) const;
80
81
};

82
83
84
85
class LinearLS : public FitLeastSquare {
	public:
		LinearLS() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "LinearLS";};
		void setDefaultGuess();
86
		double f(const double& x) const;
87
};
88

89
90
91
92
class Quadratic : public FitLeastSquare {
	public:
		Quadratic() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "Quadratic";};
		void setDefaultGuess();
93
		double f(const double& x) const;
94
95
};

96
97
/**
 * @class Fit1D
98
 * @brief A class to perform 1D regressions
99
100
 * It works on a time serie and uses either ad-hoc methods or matrix arithmetic to perform an arbitrary fit.
 * Currently, the following models are supported:
101
 * - Specific fits:
102
 *    - SimpleLinear
103
104
105
106
107
108
109
110
 *    - NoisyLinear
 * - Least Square fits:
 *    - SphericVario
 *    - LinVario
 *    - ExpVario
 *    - RatQuadVario
 *    - LinearLS
 *    - Quadratic
111
 *
112
113
114
 * The various variogram models can be found in
 * <i>"Statistics for spatial data"</i>, Noel A. C. Cressie, John Wiley & Sons, revised edition, 1993, pp63.
 *
115
116
117
118
119
120
121
 *
 * @ingroup stats
 * @author Mathias Bavay
 * @date   2011-01-20
 */
class Fit1D {
 	public:
122
123
		///Keywords for regression model
		typedef enum REGRESSION {
124
			ZERO, ///< always return zero (this is a way to disable detrending)
125
			SIMPLE_LINEAR, ///< basic, cheap linear fit
126
			NOISY_LINEAR, ///< same as SIMPLE_LINEAR but trying to remove outliers
127
128
129
130
131
132
133
134
			LINVARIO, ///< linear variogram
			EXPVARIO, ///< exponential variogram
			SPHERICVARIO, ///< spherical variogram
			RATQUADVARIO, ///< rational quadratic variogram
			LINEARLS, ///< linear, using least squares
			QUADRATIC ///< quadratic
		} regression;

135
136
137
138
		/**
		* @brief Empty Constructor. The model must be set afterwards.
		* If the model has not been set before calling other methods, a NULL pointer exception will be thrown.
		*/
139
		Fit1D() : model(NULL) {};
140

141
142
		/**
		* @brief Constructor.
143
		* @param regType regression model to use
144
145
		* @param in_X vector of data points abscissae
		* @param in_Y vector of data points ordinates
146
		* @param updatefit should the fit be redone? (default=true, otherwise you must manually call fit())
147
		*/
148
		Fit1D(const regression& regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit=true);
149

150
151
152
153
154
155
156
157
158
		/**
		* @brief Constructor for user provided model type.
		* @param regType regression model to use
		* @param in_X vector of data points abscissae
		* @param in_Y vector of data points ordinates
		* @param updatefit should the fit be redone? (default=true, otherwise you must manually call fit())
		*/
		Fit1D(const std::string& regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit=true);

159
160
161
162
163
164
165
		/**
		* @brief Copy constructor.
		* @param i_fit Object to copy
		*/
		Fit1D(const Fit1D& i_fit);

		~Fit1D() {delete model;};
166

167
168
169
170
171
172
		/**
		* @brief Set or reset the regression model.
		* @param i_regType regression model to use
		* @param in_X vector of data points abscissae
		* @param in_Y vector of data points ordinates
		* @param updatefit should the fit be redone? (default=true, otherwise you must manually call fit())
173
		* @return false if could not compute the parameters. if !updatefit, always return true
174
		*/
175
		bool setModel(const regression& i_regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit=true);
176

177
178
179
180
181
182
		/**
		* @brief Set or reset the regression model.
		* @param i_regType regression model to use
		* @param in_X vector of data points abscissae
		* @param in_Y vector of data points ordinates
		* @param updatefit should the fit be redone? (default=true, otherwise you must manually call fit())
183
		* @return false if could not compute the parameters. if !updatefit, always return true
184
		*/
185
		bool setModel(const std::string& i_regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit=true);
186

187
188
		/**
		* @brief Provide a set of initial values for the model parameters.
189
190
		* The model can be used right after providing the guesses, and it would use those guesses as parameters,
		* thus allowing the user to force his model parameters.
191
192
		* @param lambda_in one initial value per model parameter
		*/
193
		void setGuess(const std::vector<double>& lambda_in) {model->setGuess(lambda_in);};
194

195
196
197
198
199
200
201
		/**
		* @brief Set a forced lapse rate for linear regressions
		* This will throw an exception for all other regression models!
		* @param lapse_rate lapse rate to set
		*/
		void setLapseRate(const double& lapse_rate) {model->setLapseRate(lapse_rate);};

202
203
		/**
		* @brief Compute the regression parameters
204
		* @return false if could not compute the parameters
205
		*/
206
		bool fit() {return model->fit();};
207
208
209
210
211
212
213

		/**
		* @brief Calculate a value using the computed least square fit.
		* The fit has to be computed before.
		* @param x abscissa
		* @return f(x) using the computed least square fit
		*/
214
		double f(const double& x) const {return model->f(x);};
215
216
217
218
219
220

		/**
		* @brief Calculate the parameters of the fit.
		* The fit has to be computed before.
		* @param coefficients vector containing the coefficients
		*/
221
		void getParams(std::vector<double>& coefficients) const {model->getParams(coefficients);};
222
223
224
225
226

		/**
		* @brief Return the name of the fit model.
		* @return model name
		*/
227
		std::string getName() const {return model->getName();};
228
229
230
231
232
233

		/**
		* @brief Return a string of information about the fit.
		* The fit has to be computed before.
		* @return info string
		*/
234
		std::string getInfo() const {return model->getInfo();};
235

236
237
		Fit1D& operator =(const Fit1D& source);

238
239
240
241
242
243
244
245
		/**
		* @brief Calculate a value using the computed least square fit.
		* The fit has to be computed before.
		* @param x abscissa
		* @return f(x) using the computed least square fit
		*/
		double operator ()(const double& x) const {return model->f(x);};

246
247
		std::string toString() const {return model->toString();};

248
	private:
249
		FitModel *model;
250
251
252
253
254
};

} //end namespace

#endif