WSL/SLF GitLab Repository

Commit 5f61e94a authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Preparing the way for the new spatial interpolation architecture: the 1D...

Preparing the way for the new spatial interpolation architecture: the 1D regressions should now offer the necessary features (but this is a first shot at it, it will need to be redone in a more generic way later).
parent f8b0d568
......@@ -25,9 +25,6 @@ using namespace std;
namespace mio {
Fit1D::Fit1D() : model(NULL) {}
//default constructor
Fit1D::Fit1D(const regression& regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit) : model(NULL) {
setModel(regType, in_X, in_Y, updatefit);
}
......@@ -53,7 +50,8 @@ void Fit1D::setModel(const std::string& i_regType, const std::vector<double>& in
regression regType;
if(i_regType=="ZERO") regType=ZERO;
else if(i_regType=="SIMPLE_LINEAR") regType=SIMPLE_LINEAR;
else if(i_regType=="NOISYLINEAR") regType=NOISYLINEAR;
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;
......@@ -73,7 +71,8 @@ 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==NOISYLINEAR) model=new NoisyLinear;
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;
......@@ -91,6 +90,11 @@ void SimpleLinear::setData(const std::vector<double>& in_X, const std::vector<do
X = in_X;
Y = in_Y;
fit_ready = false;
}
bool SimpleLinear::checkInputs()
{
//check input data consistency
nPts=X.size();
......@@ -104,30 +108,79 @@ void SimpleLinear::setData(const std::vector<double>& in_X, const std::vector<do
stringstream ss;
ss << "Only " << nPts << " data points for " << regname << " regression model.";
ss << " Expecting at least " << min_nb_pts << " for this model!\n";
throw InvalidArgumentException(ss.str(), AT);
infoString = ss.str();
return false;
}
fit_ready = false;
return true;
}
double SimpleLinear::f(const double& x) {
return Lambda.at(0)*x + Lambda.at(1);
}
bool SimpleLinear::fit() {
bool SimpleLinear::fit()
{
if(!checkInputs())
return false;
Lambda.clear();
double a,b,r;
std::stringstream mesg;
Interpol1D::LinRegression(X, Y, a, b, r, mesg);
if(fixed_lapse_rate==IOUtils::nodata) {
Interpol1D::LinRegression(X, Y, a, b, r, mesg);
mesg << "Computed regression with " << regname << " model - r=" << r;
} else {
a = fixed_lapse_rate;
Interpol1D::LinRegressionFixedRate(X, Y, a, b, r, mesg);
mesg << "Computed regression with " << regname << " model ";
mesg << "(fixed lapse rate=" << a << ") - r=" << r;
}
Lambda.push_back(a);
Lambda.push_back(b);
infoString = mesg.str();
fit_ready = true;
return true;
}
double FracLinear::f(const double& x) {
return Lambda.at(0)*x*x + Lambda.at(1);
}
bool FracLinear::fit()
{
if(!checkInputs())
return false;
Lambda.clear();
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);
mesg << "Computed regression with " << regname << " model - r=" << r;
} else {
a = fixed_lapse_rate;
Interpol1D::LinRegressionFixedRate(X_sq, Y, a, b, r, mesg);
mesg << "Computed regression with " << regname << " model ";
mesg << "(fixed lapse rate=" << a << ") - r=" << r;
}
Lambda.push_back(a);
Lambda.push_back(b);
mesg << "Computed regression with " << regname << " model - r=" << r;
infoString = mesg.str();
fit_ready = true;
return true;
}
bool NoisyLinear::fit() {
bool NoisyLinear::fit()
{
if(!checkInputs())
return false;
Lambda.clear();
double a,b,r;
std::stringstream mesg;
......
......@@ -35,10 +35,21 @@ class Zero : public FitModel {
class SimpleLinear : public FitModel {
public:
SimpleLinear() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "SimpleLinear";};
SimpleLinear() : fixed_lapse_rate(IOUtils::nodata) {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "SimpleLinear";};
void setData(const std::vector<double>& in_X, const std::vector<double>& in_Y);
bool fit();
double f(const double& x);
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;
};
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 {
......@@ -95,8 +106,9 @@ class Quadratic : public FitLeastSquare {
* 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:
* - Specific fits:
* - SimpleLinear:
* - SimpleLinear
* - NoisyLinear
* - fracLinear
* - Least Square fits:
* - SphericVario
* - LinVario
......@@ -119,7 +131,8 @@ class Fit1D {
typedef enum REGRESSION {
ZERO, ///< always return zero (this is a way to disable detrending)
SIMPLE_LINEAR, ///< basic, cheap linear fit
NOISYLINEAR, ///< same as SIMPLE_LINEAR but trying to remove outliers
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
......@@ -132,7 +145,7 @@ class Fit1D {
* @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.
*/
Fit1D();
Fit1D() : model(NULL) {};
/**
* @brief Constructor.
......@@ -186,6 +199,13 @@ class Fit1D {
*/
void setGuess(const std::vector<double>& lambda_in) {model->setGuess(lambda_in);};
/**
* @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);};
/**
* @brief Compute the regression parameters
* @return false if could not find the parameters
......@@ -222,6 +242,14 @@ class Fit1D {
Fit1D& operator =(const Fit1D& source);
/**
* @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);};
private:
FitModel *model;
};
......
......@@ -31,11 +31,7 @@ void FitModel::getParams(std::vector<double>& o_coefficients) {
}
std::string FitModel::getInfo() {
if(fit_ready==true) {
return infoString;
} else {
throw InvalidArgumentException("The regression has not yet being computed!", AT);
}
return infoString;
}
void FitModel::setGuess(const std::vector<double>& lambda_in) {
......@@ -74,14 +70,9 @@ const double FitLeastSquare::delta_init_rel = 0.2; //initial delta, relative
const double FitLeastSquare::eps_conv = 1e-6; //convergence criteria
const unsigned int FitLeastSquare::max_iter = 50; //maximum number of iterations
//default constructor
FitLeastSquare::FitLeastSquare() {
}
void FitLeastSquare::setData(const std::vector<double>& in_X, const std::vector<double>& in_Y) {
X = in_X;
Y = in_Y;
checkInputs();
Interpol1D::sort(X, Y); //sort the data by increasing X
setDefaultGuess();
fit_ready = false;
......@@ -94,13 +85,14 @@ void FitLeastSquare::setDefaultGuess() {
}
bool FitLeastSquare::fit() {
checkInputs();
return computeFit();
}
////////////////////////////////////////////////////////////
//// End of public methods
void FitLeastSquare::checkInputs()
bool FitLeastSquare::checkInputs()
{
nPts=X.size();
......@@ -114,8 +106,11 @@ void FitLeastSquare::checkInputs()
stringstream ss;
ss << "Only " << nPts << " data points for " << regname << " regression model.";
ss << " Expecting at least " << min_nb_pts << " for this model!\n";
throw InvalidArgumentException(ss.str(), AT);
infoString = ss.str();
return false;
}
return true;
}
bool FitLeastSquare::computeFit() {
......
......@@ -33,6 +33,7 @@ class FitModel {
virtual ~FitModel() {};
virtual void setData(const std::vector<double>& in_X, const std::vector<double>& in_Y) = 0;
void setGuess(const std::vector<double>& lambda_in);
virtual void setLapseRate(const double& /*lapse_rate*/) {throw InvalidArgumentException("Lapse rates can only be forced for linear regressions!", AT);};
virtual bool fit() = 0;
virtual double f(const double& x) = 0;
void getParams(std::vector<double>& o_coefficients);
......@@ -40,6 +41,8 @@ class FitModel {
std::string getInfo();
FitModel& operator =(const FitModel& source);
protected:
virtual bool checkInputs() {return true;};
std::vector<double> Lambda; //parameters of the fit
std::vector<double> X; //X of input data set to fit
std::vector<double> Y; //Y of input data set to fit
......@@ -48,7 +51,7 @@ class FitModel {
size_t nPts; //number of data points
size_t nParam; //number of parameters
size_t min_nb_pts; //minimum number of data points
bool fit_ready;
bool fit_ready; //set to false if fit() must be called before using the fit
};
/**
......@@ -62,7 +65,7 @@ class FitModel {
*/
class FitLeastSquare : public FitModel {
public:
FitLeastSquare();
FitLeastSquare() {};
void setData(const std::vector<double>& in_X, const std::vector<double>& in_Y);
bool fit();
virtual double f(const double& x) = 0;
......@@ -71,7 +74,7 @@ class FitLeastSquare : public FitModel {
virtual void setDefaultGuess(); //set defaults guess values. Called by setData
private:
void checkInputs();
bool checkInputs();
void initLambda();
void initDLambda(Matrix& dLambda) const;
double getDelta(const double& var) const;
......
......@@ -406,6 +406,57 @@ void Interpol1D::LinRegression(const std::vector<double>& X, const std::vector<d
}
}
/**
* @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
*/
void Interpol1D::LinRegressionFixedRate(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg)
{ //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;
for (size_t i=0; i<n; i++) {
if(X[i]!=IOUtils::nodata && Y[i]!=IOUtils::nodata) {
x_avg += X[i];
y_avg += Y[i];
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
for (size_t i=0; i<n; i++) {
if(X[i]!=IOUtils::nodata && Y[i]!=IOUtils::nodata) {
SSR += Optim::pow2( Y[i] - (a*X[i]+b) );
TSS += Optim::pow2( Y[i] );
}
}
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";
}
}
/**
* @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)
......
......@@ -50,6 +50,7 @@ class Interpol1D {
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 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);
......
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