WSL/SLF GitLab Repository

Commit 8a8dce1f authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The nodata points are now always properly filtered out when computing a least...

The nodata points are now always properly filtered out when computing a least square fit. The code has been a little reshaped in order to be closer to the reference document that has also been cited. Some comments have been added. The quadratic fit was wrong and has been fixed.

The statistical tests now cover the least square fits and ad-hoc linear regressions.

The ARC plugin properly builds the file names for reading and writing grids with the "A3D_VIEW" option. 
parent 5dd995f2
......@@ -78,7 +78,16 @@ bool Fit1D::setModel(const regression& regType, const std::vector<double>& in_X,
if(regType==LINEARLS) model=new LinearLS;
if(regType==QUADRATIC) model=new Quadratic;
model->setData(in_X, in_Y);
//remove nodata points
std::vector<double> X, Y;
for (size_t ii=0; ii<in_X.size(); ii++) {
if(in_X[ii]!=IOUtils::nodata && in_Y[ii]!=IOUtils::nodata) {
X.push_back( in_X[ii] );
Y.push_back( in_Y[ii] );
}
}
model->setData(X, Y);
if(updatefit)
return fit();
else
......@@ -294,7 +303,7 @@ double Quadratic::f(const double& x) const {
}
void Quadratic::setDefaultGuess() {
std::vector<double> der = Interpol1D::derivative(X, Y);
std::vector<double> der( Interpol1D::derivative(X, Y) );
const double acc = 0.5 * Interpol1D::arithmeticMean( Interpol1D::derivative(X, der) );
double xzero=der[0];
size_t xzero_idx=0;
......
......@@ -88,7 +88,7 @@ class LinearLS : public FitLeastSquare {
class Quadratic : public FitLeastSquare {
public:
Quadratic() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "Quadratic";};
Quadratic() {fit_ready = false; nParam = 3; min_nb_pts = 4; regname = "Quadratic";};
void setDefaultGuess();
double f(const double& x) const;
};
......
......@@ -135,6 +135,7 @@ bool FitLeastSquare::checkInputs()
return true;
}
//see http://mathworld.wolfram.com/NonlinearLeastSquaresFitting.html
bool FitLeastSquare::computeFit() {
double max_delta;
initLambda();
......@@ -147,6 +148,11 @@ bool FitLeastSquare::computeFit() {
unsigned int iter = 0;
do {
iter++;
//set dBeta matrix
for(size_t m=1; m<=nPts; m++) {
dBeta(m,1) = Y[m-1] - f(X[m-1]); //X and Y are vectors
}
//set A matrix
for(size_t m=1; m<=nPts; m++) {
for(size_t n=1; n<=nParam; n++) {
......@@ -155,11 +161,6 @@ bool FitLeastSquare::computeFit() {
}
}
//set dBeta matrix
for(size_t m=1; m<=nPts; m++) {
dBeta(m,1) = Y[m-1] - f(X[m-1]); //X and Y are vectors
}
//calculate parameters deltas
const Matrix a = A.getT() * A;
const Matrix b = A.getT() * dBeta;
......@@ -167,9 +168,9 @@ bool FitLeastSquare::computeFit() {
//apply the deltas to the parameters, record maximum delta
max_delta = 0.;
for(size_t m=1; m<=nParam; m++) {
Lambda[m-1] += dLambda(m,1); //Lambda is a vector
if( fabs(dLambda(m,1))>max_delta ) max_delta=fabs(dLambda(m,1));
for(size_t n=1; n<=nParam; n++) {
Lambda[n-1] += dLambda(n,1); //Lambda is a vector
if( fabs(dLambda(n,1))>max_delta ) max_delta=fabs(dLambda(n,1));
}
} while (max_delta>eps_conv && iter<max_iter);
......@@ -229,9 +230,10 @@ double FitLeastSquare::DDer(const double& x, const size_t& index) {
const double y1 = f(x);
Lambda[index-1] = v2;
const double y2 = f(x);
Lambda[index-1] = var;
return (y2-y1)/(v2-v1);
Lambda[index-1] = var; //restore initial value
return (y2-y1)/(2.*delta);
}
} //namespace
......@@ -58,7 +58,8 @@ class FitModel {
/**
* @class FitLeastSquare
* @brief A class to perform non-linear least square fitting.
* It works on a time serie and uses matrix arithmetic to perform an arbitrary fit.
* It works on a time serie and uses matrix arithmetic to perform an arbitrary fit
* (see http://mathworld.wolfram.com/NonlinearLeastSquaresFitting.html).
*
* @ingroup stats
* @author Mathias Bavay
......
......@@ -269,6 +269,8 @@ void ARCIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& par
ext="swr";
else if (parameter==MeteoGrids::ILWR)
ext="lwr";
else if(parameter==MeteoGrids::DEM)
ext="asc";
else {
ext = MeteoGrids::getParameterName(parameter);
IOUtils::toLower(ext);
......@@ -382,12 +384,18 @@ void ARCIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameter
{
//the path will be added by write2DGrid
if(a3d_view_out) {
// the A3D grid viewer looks for the following extensions:
//sdp, tss, swr, lwr, swe, alb, wet
string ext;
if(parameter==MeteoGrids::HS) {
if (parameter==MeteoGrids::HS)
ext="sdp";
} else if(parameter==MeteoGrids::DEM) {
else if (parameter==MeteoGrids::ISWR)
ext="swr";
else if (parameter==MeteoGrids::ILWR)
ext="lwr";
else if(parameter==MeteoGrids::DEM)
ext="asc";
} else {
else {
ext = MeteoGrids::getParameterName(parameter);
IOUtils::toLower(ext);
}
......
......@@ -237,6 +237,93 @@ bool check_derivative(const vector<double>& x, const vector<double>& y) {
return status;
}
bool check_regressions(const vector<double>& x, const vector<double>& y) {
bool status = true;
Fit1D fit(Fit1D::SIMPLE_LINEAR, x, y);
std::vector<double> coeff;
fit.getParams(coeff);
if (!IOUtils::checkEpsilonEquality(coeff[0], -0.500941, 1e-6) ||
!IOUtils::checkEpsilonEquality(coeff[1], 12.810614, 1e-6)) {
std::cout << "wrong results for simple_linear regression: returned ";
std::cout << "a=" << setprecision(12) << coeff[0] << " and b=" << coeff[1] << "\n";
status = false;
}
fit.setModel(Fit1D::LINEARLS, x, y); //this should return the same values
fit.getParams(coeff);
if (!IOUtils::checkEpsilonEquality(coeff[0], -0.500941, 1e-6) ||
!IOUtils::checkEpsilonEquality(coeff[1], 12.810614, 1e-6)) {
std::cout << "wrong results for linear_ls regression: returned ";
std::cout << "a=" << setprecision(12) << coeff[0] << " and b=" << coeff[1] << "\n";
status = false;
}
fit.setModel(Fit1D::NOISY_LINEAR, x, y);
const double y1 = fit.f(-700.), yres1 = 372.9869;
const double y2 = fit.f(-31.5), yres2 = 141.9284;
const double y3 = fit.f(45.), yres3 = 115.4871;
const double y4 = fit.f(250.781), yres4 = 44.3615;
const double y5 = fit.f(500.), yres5 = -41.7778;
if (!IOUtils::checkEpsilonEquality(y1, yres1, 1e-3) ||
!IOUtils::checkEpsilonEquality(y2, yres2, 1e-3) ||
!IOUtils::checkEpsilonEquality(y3, yres3, 1e-3) ||
!IOUtils::checkEpsilonEquality(y4, yres4, 1e-3) ||
!IOUtils::checkEpsilonEquality(y5, yres5, 1e-3) ) {
std::cout << "wrong results for noisy_linear regression:\n\texpected\treturned\n";
std::cout << setprecision(12);
std::cout << "\t" << yres1 << "\t\t" << y1 << "\n";
std::cout << "\t" << yres2 << "\t\t" << y2 << "\n";
std::cout << "\t" << yres3 << "\t\t" << y3 << "\n";
std::cout << "\t" << yres4 << "\t\t" << y4 << "\n";
std::cout << "\t" << yres5 << "\t\t" << y5 << "\n";
status = false;
}
fit.setModel(Fit1D::QUADRATIC, x, y);
fit.getParams(coeff);
if (!IOUtils::checkEpsilonEquality(coeff[0], 0.0016668, 1e-6) ||
!IOUtils::checkEpsilonEquality(coeff[1], -0.3605318, 1e-6) ||
!IOUtils::checkEpsilonEquality(coeff[2], -98.3815970, 1e-6)) {
std::cout << "wrong results for quadratic regression: returned ";
std::cout << "a=" << setprecision(12) << coeff[0] << " b=" << coeff[1] << " c=" << coeff[2] << "\n";
status = false;
}
//set the parameters and get the values
vector<double> lambda(3);
lambda[0] = 0.05; lambda[1] = 0.70; lambda[2] = 40.;
vector<double> X(11), Y(11);
X[0] = 2.; Y[0] = 0.102456;
X[1] = 5.; Y[1] = 0.180566;
X[2] = 7.; Y[2] = 0.231874;
X[3] = 10.; Y[3] = 0.307031;
X[4] = 14.; Y[4] = 0.402494;
X[5] = 19.; Y[5] = 0.51124;
X[6] = 22.; Y[6] = 0.569269;
X[7] = 39.; Y[7] = 0.749349;
X[8] = 60.; Y[8] = 0.75;
X[9] = 90.; Y[9] = 0.75;
X[10] = 100.; Y[10] = 0.75;
fit.setModel(Fit1D::SPHERICVARIO, X, Y, false);
fit.setGuess(lambda);
for (size_t ii=0; ii<X.size(); ii++) {
if ( !IOUtils::checkEpsilonEquality(fit.f(X[ii]), Y[ii], 1e-5) ) {
std::cout << "Wrong result when setting parameters: expected " << Y[ii];
std::cout << " received " << fit.f(X[ii]) << " instead\n";
status = false;
}
}
if (status)
std::cout << "Regressions: success\n";
else
std::cout << "Regressions: failed\n";
return status;
}
int main() {
vector<double> x,y;
//cr_rand_vectors(x, y);
......@@ -248,8 +335,9 @@ int main() {
const bool covariance_status = check_covariance(x,y);
const bool der_status = check_derivative(x,y);
const bool quantiles_status = check_quantiles(x);
const bool regressions_status = check_regressions(x, y);
if(!basics_status || !sort_status || !bin_status || !quantiles_status || !covariance_status || !der_status)
if(!basics_status || !sort_status || !bin_status || !quantiles_status || !covariance_status || !der_status || !regressions_status)
throw IOException("Statistical functions error!", AT);
......
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