WSL/SLF GitLab Repository

Commit 597e6f48 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Moving forward with the trend/residuals approach for the spatial...

Moving forward with the trend/residuals approach for the spatial interpolations: here are the necessary changes in the supporting infrastructure: basically, better handling of forced lapse rates and fractional lapse rates, proper toString() methods for the fit1D object and components
parent 3b17026d
......@@ -51,7 +51,6 @@ void Fit1D::setModel(const std::string& i_regType, const std::vector<double>& in
if(i_regType=="ZERO") regType=ZERO;
else if(i_regType=="SIMPLE_LINEAR") regType=SIMPLE_LINEAR;
else if(i_regType=="NOISY_LINEAR") regType=NOISY_LINEAR;
else if(i_regType=="FRAC_LINEAR") regType=FRAC_LINEAR;
else if(i_regType=="LINVARIO") regType=LINVARIO;
else if(i_regType=="EXPVARIO") regType=EXPVARIO;
else if(i_regType=="SPHERICVARIO") regType=SPHERICVARIO;
......@@ -72,7 +71,6 @@ void Fit1D::setModel(const regression& regType, const std::vector<double>& in_X,
if(regType==ZERO) model=new Zero;
if(regType==SIMPLE_LINEAR) model=new SimpleLinear;
if(regType==NOISY_LINEAR) model=new NoisyLinear;
if(regType==FRAC_LINEAR) model=new FracLinear;
if(regType==LINVARIO) model=new LinVario;
if(regType==EXPVARIO) model=new ExpVario;
if(regType==SPHERICVARIO) model=new SphericVario;
......@@ -133,7 +131,7 @@ bool SimpleLinear::fit()
mesg << "Computed regression with " << regname << " model - r=" << r;
} else {
a = fixed_lapse_rate;
Interpol1D::LinRegressionFixedRate(X, Y, a, b, r, mesg);
Interpol1D::LinRegression(X, Y, a, b, r, mesg, true);
mesg << "Computed regression with " << regname << " model ";
mesg << "(fixed lapse rate=" << a << ") - r=" << r;
}
......@@ -144,11 +142,7 @@ bool SimpleLinear::fit()
return true;
}
double FracLinear::f(const double& x) {
return Lambda.at(0)*x*x + Lambda.at(1);
}
bool FracLinear::fit()
bool NoisyLinear::fit()
{
if(!checkInputs())
return false;
......@@ -157,15 +151,12 @@ bool FracLinear::fit()
double a,b,r;
std::stringstream mesg;
std::vector<double> X_sq(X);
std::transform(X_sq.begin(), X_sq.end(), X_sq.begin(), X_sq.begin(), multiplies<double>());
if(fixed_lapse_rate==IOUtils::nodata) {
Interpol1D::LinRegression(X_sq, Y, a, b, r, mesg);
Interpol1D::NoisyLinRegression(X, Y, a, b, r, mesg);
mesg << "Computed regression with " << regname << " model - r=" << r;
} else {
a = fixed_lapse_rate;
Interpol1D::LinRegressionFixedRate(X_sq, Y, a, b, r, mesg);
Interpol1D::NoisyLinRegression(X, Y, a, b, r, mesg, true);
mesg << "Computed regression with " << regname << " model ";
mesg << "(fixed lapse rate=" << a << ") - r=" << r;
}
......@@ -176,23 +167,6 @@ bool FracLinear::fit()
return true;
}
bool NoisyLinear::fit()
{
if(!checkInputs())
return false;
Lambda.clear();
double a,b,r;
std::stringstream mesg;
Interpol1D::NoisyLinRegression(X, Y, a, b, r, mesg);
Lambda.push_back(a);
Lambda.push_back(b);
mesg << "Computed regression with " << regname << " model - r=" << r;
infoString = mesg.str();
fit_ready = true;
return true;
}
//regression models using the standard least square algorithm
double SphericVario::f(const double& x) {
//c0>=0, cs>=0, as>=0
......
......@@ -45,13 +45,6 @@ class SimpleLinear : public FitModel {
double fixed_lapse_rate;
};
class FracLinear : public SimpleLinear {
public:
FracLinear() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "SimpleLinear";};
bool fit();
double f(const double& x);
};
class NoisyLinear : public SimpleLinear {
public:
NoisyLinear() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "NoisyLinear";};
......@@ -108,7 +101,6 @@ class Quadratic : public FitLeastSquare {
* - Specific fits:
* - SimpleLinear
* - NoisyLinear
* - fracLinear
* - Least Square fits:
* - SphericVario
* - LinVario
......@@ -132,7 +124,6 @@ class Fit1D {
ZERO, ///< always return zero (this is a way to disable detrending)
SIMPLE_LINEAR, ///< basic, cheap linear fit
NOISY_LINEAR, ///< same as SIMPLE_LINEAR but trying to remove outliers
FRAC_LINEAR, ///< "linear" but with a fractional lapse rate (ultimately, it is quadratic)
LINVARIO, ///< linear variogram
EXPVARIO, ///< exponential variogram
SPHERICVARIO, ///< spherical variogram
......@@ -250,6 +241,8 @@ class Fit1D {
*/
double operator ()(const double& x) const {return model->f(x);};
std::string toString() const {return model->toString();};
private:
FitModel *model;
};
......
......@@ -17,6 +17,7 @@
*/
#include <meteoio/meteostats/libfit1DCore.h>
#include <meteoio/meteostats/libinterpol1D.h>
#include <algorithm>
#include <cmath>
using namespace std;
......@@ -62,6 +63,27 @@ FitModel& FitModel::operator =(const FitModel& source) {
return *this;
}
std::string FitModel::toString() const {
std::ostringstream os;
os << "<FitModel>\n";
os << regname << " model with " << nParam << " (npts >= " << min_nb_pts << ")\n";
const double Xmin = *std::min_element(X.begin(),X.end());
const double Xmax = *std::max_element(X.begin(),X.end());
const double Ymin = *std::min_element(Y.begin(),Y.end());
const double Ymax = *std::max_element(Y.begin(),Y.end());
os << nPts << " calibration point(s) in\t";
os << "X[" << Xmin << "-" << Xmax << "]\tY[" << Ymin << "-" << Ymax << "]\n";
if(!Lambda.empty()) {
os << "Model parameters: \t";
for(size_t ii=0; ii<Lambda.size(); ii++) os << Lambda[ii] << " ";
os << "\n";
}
os << "</FitModel>\n";
return os.str();
}
////////////////////////////////////////////////////////////
//// Least square fit class
const double FitLeastSquare::lambda_init = 1.; //initial default guess
......
......@@ -26,7 +26,7 @@
namespace mio {
//purely virtual class to use as an interface
//class to use as an interface
class FitModel {
public:
FitModel() : Lambda(), X(), Y(), infoString(), regname(), nPts(0), nParam(0), min_nb_pts(0), fit_ready(false) {};
......@@ -40,6 +40,7 @@ class FitModel {
std::string getName() {return regname;};
std::string getInfo();
FitModel& operator =(const FitModel& source);
std::string toString() const;
protected:
virtual bool checkInputs() {return true;};
......
......@@ -351,9 +351,16 @@ double Interpol1D::covariance(const std::vector<double>& X, const std::vector<do
* @param b origin of the linear regression
* @param r absolute value of linear regression coefficient
* @param mesg information message if something fishy is detected
* @param fixed_rate force the lapse rate? (default=false)
*/
void Interpol1D::LinRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg)
{ //check arguments
void Interpol1D::LinRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg, const bool& fixed_rate)
{
if(fixed_rate) {
LinRegressionFixedRate(X, Y, a, b, r, mesg);
return;
}
//check arguments
const size_t n = X.size();
if(n==0)
throw NoAvailableDataException("Trying to calculate linear regression with no data points", AT);
......@@ -466,9 +473,10 @@ void Interpol1D::LinRegressionFixedRate(const std::vector<double>& X, const std:
* @param B origin of the linear regression
* @param R linear regression coefficient
* @param mesg information message if something fishy is detected
* @param fixed_rate force the lapse rate? (default=false)
* @return EXIT_SUCCESS or EXIT_FAILURE
*/
int Interpol1D::NoisyLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, double& A, double& B, double& R, std::stringstream& mesg)
int Interpol1D::NoisyLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, double& A, double& B, double& R, std::stringstream& mesg, const bool& fixed_rate)
{
//finds the linear regression for points (x,y,z,Value)
const double r_thres = 0.7;
......@@ -476,20 +484,9 @@ int Interpol1D::NoisyLinRegression(const std::vector<double>& in_X, const std::v
//we want at least 4 points AND 85% of the initial data set kept in the regression
const unsigned int min_dataset = (unsigned int)Optim::floor( 0.85*(double)nb_pts );
const unsigned int min_pts = (min_dataset>4)? min_dataset : 4;
double a,b,r;
if (nb_pts==2) {
mesg << "[W] only two points for linear regression!\n";
}
if (nb_pts<2) { //this should not be needed, we should have refrained from calling LinRegression in such a case
mesg << "[E] Not enough data points for linear regression!\n";
A=0.;
if (!in_X.empty()) B=in_X[0]; else B=0.;
R=1.;
return EXIT_FAILURE;
}
double a=A,b,r; //a needs to be initiallized to A in case of fixed_rate
LinRegression(in_X, in_Y, A, B, R, mesg);
LinRegression(in_X, in_Y, A, B, R, mesg, fixed_rate);
if(R>=r_thres)
return EXIT_SUCCESS;
......@@ -503,11 +500,11 @@ int Interpol1D::NoisyLinRegression(const std::vector<double>& in_X, const std::v
for (size_t i=0; i<nb_pts; i++) {
//invalidating alternatively each point
const double Y_tmp=Y[i]; Y[i]=IOUtils::nodata;
LinRegression(X, Y, a, b, r, mesg);
LinRegression(X, Y, a, b, r, mesg, fixed_rate);
Y[i]=Y_tmp;
if (fabs(r)>fabs(R)) {
A=a;
A=a; //if fixed_rate, a=A anyway...
B=b;
R=r;
index_bad=i;
......
......@@ -49,14 +49,14 @@ class Interpol1D {
static double std_dev(const std::vector<double>& X);
static double covariance(const std::vector<double>& z1, const std::vector<double>& z2);
static void LinRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg);
static void LinRegressionFixedRate(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg);
static int NoisyLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, double& A, double& B, double& R, std::stringstream& mesg);
static void LinRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg, const bool& fixed_rate=false);
static int NoisyLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, double& A, double& B, double& R, std::stringstream& mesg, const bool& fixed_rate=false);
static int twoLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, const double& bilin_inflection, std::vector<double>& coeffs);
static void LogRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg);
static void ExpRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg);
private:
static void LinRegressionFixedRate(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg);
static bool pair_comparator(const std::pair<double, double>& l, const std::pair<double, double>& r);
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment