WSL/SLF GitLab Repository

Commit 377d6e57 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The Krigging is now working. The current limitations are: 1)there is no...

The Krigging is now working. The current limitations are: 1)there is no maximum range, so each available station would take part to the variogram fit. 2)the covariance is currently NOT computed on past time series, limiting the relevance of the variogram. Practically, the variogram fit that always ends up being used is LINVARIO.

All matrix and fit methods that used to return void but could also throw exceptions now return a boolean to indicate if things went well (ie a matrix could be inverted, etc). They still throw exceptions for logic errors (incompatible dimensions, etc).
parent 84602c43
......@@ -53,6 +53,8 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& i_algo
return new SimpleWindInterpolationAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "ODKRIG"){// ordinary kriging
return new OrdinaryKrigingAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "ODKRIG_LAPSE"){// ordinary kriging with lapse rate
return new LapseOrdinaryKrigingAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "USER"){// read user provided grid
return new USERInterpolation(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "HNW_SNOW"){// precipitation interpolation according to (Magnusson, 2010)
......@@ -688,6 +690,58 @@ void SnowHNWInterpolation::calculate(Grid2DObject& grid)
if(new_mean!=0.) grid.grid2D *= orig_mean/new_mean;
}
void OrdinaryKrigingAlgorithm::getDataForVariogram(std::vector<double> &distData, std::vector<double> &variData)
{
distData.clear();
variData.clear();
for(unsigned int j=0; j<nrOfMeasurments; j++) {
const Coords& st1 = vecMeta[j].position;
const double x1 = st1.getEasting();
const double y1 = st1.getNorthing();
const double val1 = vecData[j];
for(unsigned int i=0; i<j; i++) {
//compute distance between stations
const Coords& st2 = vecMeta[i].position;
const double val2 = vecData[i];
const double DX = x1-st2.getEasting();
const double DY = y1-st2.getNorthing();
const double distance = Optim::fastSqrt_Q3(Optim::pow2(DX) + Optim::pow2(DY));
distData.push_back(distance);
variData.push_back( 0.5*Optim::pow2(val1-val2) );
}
}
}
bool OrdinaryKrigingAlgorithm::computeVariogram()
{//return variogram fit of covariance between stations i and j
std::vector<double> distData, variData;
getDataForVariogram(distData, variData);
std::vector<string> vario_types( vecArgs );
if(vario_types.empty()) vario_types.push_back("LINVARIO");
/*std::vector<string> vario_types;
vario_types.push_back("SPHERICVARIO");
vario_types.push_back("EXPVARIO");
vario_types.push_back("LINVARIO");*/
size_t args_index=0;
do {
const string vario_model = IOUtils::strToUpper( vario_types[args_index] );
const bool status = variogram.setModel(vario_model, distData, variData);
if(status) {
info << " - " << vario_model;
return true;
}
args_index++;
} while (args_index<vario_types.size());
return false;
}
void OrdinaryKrigingAlgorithm::initialize(const MeteoData::Parameters& in_param) {
param = in_param;
nrOfMeasurments = getData(param, vecData, vecMeta);
......@@ -695,11 +749,23 @@ void OrdinaryKrigingAlgorithm::initialize(const MeteoData::Parameters& in_param)
double OrdinaryKrigingAlgorithm::getQualityRating() const
{
if(nrOfMeasurments>=20) return 0.9;
if(nrOfMeasurments>=7) return 0.9;
return 0.1;
}
void OrdinaryKrigingAlgorithm::calculate(Grid2DObject& grid)
{
//optimization: getRange (from variogram fit -> exclude stations that are at distances > range (-> smaller matrix)
//or, get max range from io.ini, build variogram from this user defined max range
if (nrOfMeasurments==0)
throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
if(!computeVariogram()) //only refresh once a month, or once a week, etc
throw IOException("The variogram for parameter " + MeteoData::getParameterName(param) + " could not be computed!", AT);
Interpol2D::ODKriging(vecData, vecMeta, dem, variogram, grid);
}
void LapseOrdinaryKrigingAlgorithm::calculate(Grid2DObject& grid)
{
//optimization: getRange (from variogram fit -> exclude stations that are at distances > range (-> smaller matrix)
//or, get max range from io.ini, build variogram from this user defined max range
......@@ -708,61 +774,18 @@ void OrdinaryKrigingAlgorithm::calculate(Grid2DObject& grid)
if (vecAltitudes.empty() || nrOfMeasurments==0)
throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
Fit1D trend;
getTrend(vecAltitudes, vecData, trend);
Fit1D trend(Fit1D::NOISY_LINEAR, vecAltitudes, vecData, false);
if(!trend.fit())
throw IOException("Interpolation FAILED for parameter " + MeteoData::getParameterName(param) + ": " + trend.getInfo(), AT);
info << trend.getInfo();
detrend(trend, vecAltitudes, vecData);
computeVariogram(); //only refresh once a month, or once a week, etc
if(!computeVariogram()) //only refresh once a month, or once a week, etc
throw IOException("The variogram for parameter " + MeteoData::getParameterName(param) + " could not be computed!", AT);
Interpol2D::ODKriging(vecData, vecMeta, dem, variogram, grid);
retrend(trend, grid);
}
bool OrdinaryKrigingAlgorithm::computeVariogram()
{//return variogram fit of covariance between stations i and j
std::string vario_model;
if (vecArgs.empty()){
vario_model=std::string("LINVARIO");
} else if (vecArgs.size() == 1){
IOUtils::convertString(vario_model, vecArgs[0]);
} else { //incorrect arguments, throw an exception
throw InvalidArgumentException("Wrong number of arguments supplied for the ODKRIG algorithm", AT);
}
//give data to compute the variogram
std::vector<double> dist, vari;
for(unsigned int j=0; j<nrOfMeasurments; j++) {
const Coords& st1 = vecMeta[j].position;
const double x1 = st1.getEasting();
const double y1 = st1.getNorthing();
const double val1 = vecData[j];
for(unsigned int i=0; i<j; i++) {
//compute distance between stations
const Coords& st2 = vecMeta[i].position;
const double val2 = vecData[i];
const double DX = x1-st2.getEasting();
const double DY = y1-st2.getNorthing();
//const double distance = fastSqrt_Q3(DX*DX + DY*DY);
const double distance = sqrt(DX*DX + DY*DY);
dist.push_back(distance);
vari.push_back( 0.5*(val1-val2)*(val1-val2) );
}
}
try {
variogram.setModel(vario_model, dist, vari);
} catch(const IOException&) {
cerr << "[E] Variogram data points:\n";
for(unsigned int ii=0; ii<dist.size(); ii++) {
cerr << dist[ii] << " " << vari[ii] << endl;
}
throw IOException("The variogram model "+variogram.getName()+" could not be fitted to the available data!", AT);
}
return true;
}
} //namespace
......@@ -82,6 +82,7 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
* - WIND_CURV: the wind field (VW and DW) is interpolated using IDW_LAPSE and then altered depending on the local curvature and slope (taken from the DEM, see SimpleWindInterpolationAlgorithm)
* - HNW_SNOW: precipitation interpolation according to (Magnusson, 2011) (see SnowHNWInterpolation)
* - ODKRIG: ordinary kriging (see OrdinaryKrigingAlgorithm)
* - ODKRIG_LAPSE: ordinary kriging with lapse rate (see LapseOrdinaryKrigingAlgorithm)
* - USER: user provided grids to be read from disk (if available, see USERInterpolation)
* <!-- - LIDW_LAPSE: IDW_LAPSE restricted to a local scale (n neighbor stations, see LocalIDWLapseAlgorithm) -->
*
......@@ -447,22 +448,21 @@ class SnowHNWInterpolation : public InterpolationAlgorithm {
/**
* @class OrdinaryKrigingAlgorithm
* @brief Ordinary kriging.
* this implements ordinary krigging (see https://secure.wikimedia.org/wikipedia/en/wiki/Kriging)
* This implements ordinary krigging (see https://secure.wikimedia.org/wikipedia/en/wiki/Kriging)
* with user-selectable variogram model (see https://secure.wikimedia.org/wikipedia/en/wiki/Variogram).
* The variogram and krigging
* coefficients are re-computed fresh for each new grid (or time step). There is currently no fallback
* mechanism for variogram fit failures, which makes the whole thing not very robust.
* Therefore, USE IT AT YOUR OWN RISK!
* coefficients are re-computed fresh for each new grid (or time step).
*
* The variogram is currently computed with the current data (as 1/2*(X1-X2)^2), which makes it quite
* uninteresting... The next improvement will consist in calculating the covariances (used to build the
* variogram) from time series (thus reflecting the time-correlation between stations).
*
* The available variogram models are found in Fit1D::regression and written capitalized as optional arguments
* (by default, LINVARIO is used):
* The available variogram models are found in Fit1D::regression and given as optional arguments
* (by default, LINVARIO is used). Several models can be given, the first that can fit the data will be used
* for the current timestep:
* @code
* TA::algorithms = ODKRIG
* TA::odkrig = SPHERICVARIO
* TA::odkrig = SPHERICVARIO linvario
* @endcode
*
* @author Mathias Bavay
......@@ -478,11 +478,35 @@ class OrdinaryKrigingAlgorithm : public InterpolationAlgorithm {
virtual void initialize(const MeteoData::Parameters& in_param);
virtual double getQualityRating() const;
virtual void calculate(Grid2DObject& grid);
private:
protected:
void getDataForVariogram(std::vector<double> &distData, std::vector<double> &variData);
bool computeVariogram();
Fit1D variogram;
};
/**
* @class LapseOrdinaryKrigingAlgorithm
* @brief Ordinary kriging with detrending.
* This is very similar to OrdinaryKrigingAlgorithm but performs detrending on the data.
* @code
* TA::algorithms = ODKRIG_LAPSE
* TA::odkrig = SPHERICVARIO
* @endcode
*
* @author Mathias Bavay
*/
class LapseOrdinaryKrigingAlgorithm : public OrdinaryKrigingAlgorithm {
public:
LapseOrdinaryKrigingAlgorithm(Meteo2DInterpolator& i_mi,
const Date& i_date,
const DEMObject& i_dem,
const std::vector<std::string>& i_vecArgs,
const std::string& i_algo, IOManager& iom)
: OrdinaryKrigingAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algo, iom) {}
virtual void calculate(Grid2DObject& grid);
};
} //end namespace mio
#endif
......@@ -158,7 +158,7 @@ namespace Optim {
/**
* @brief Optimized version of c++ pow()
* This version works with positive and negative exponents and handles exponents bigger than 1.
* The relative error remains less than 6% for the benachmarks that we ran (argument between 0 and 500
* The relative error remains less than 6% for the benchmarks that we ran (argument between 0 and 500
* and exponent between -10 and +10). It is ~3.3 times faster than cmath's pow().
* Source: http://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp/
*
......
......@@ -25,12 +25,7 @@ namespace mio {
const double Matrix::epsilon = 1e-9; //for considering a determinant to be zero, etc
const double Matrix::epsilon_mtr = 1e-6; //for comparing two matrix
Matrix::Matrix() : vecData(), ncols(0), nrows(0)
{
}
Matrix::Matrix(const int& rows, const int& cols) : vecData(), ncols(0), nrows(0)
{
Matrix::Matrix(const int& rows, const int& cols) : vecData(), ncols(0), nrows(0) {
if(rows<0 || cols<0) {
std::stringstream tmp;
tmp << "Trying construct a matrix with negative dimensions: ";
......@@ -40,26 +35,10 @@ Matrix::Matrix(const int& rows, const int& cols) : vecData(), ncols(0), nrows(0)
resize((unsigned)rows,(unsigned)cols);
}
Matrix::Matrix(const unsigned int& rows, const unsigned int& cols) : vecData(rows*cols), ncols(cols), nrows(rows)
{
//resize(rows,cols);
}
Matrix::Matrix(const unsigned int& rows, const unsigned int& cols, const double& init) : vecData(rows*cols, init), ncols(cols), nrows(rows)
{
//resize(rows,cols,init);
}
Matrix::Matrix(const unsigned int& n, const double& init) : vecData(n*n, 0.), ncols(n), nrows(n)
{
//identity(n, init);
Matrix::Matrix(const unsigned int& n, const double& init) : vecData(n*n, 0.), ncols(n), nrows(n) {
for(unsigned int ii=1; ii<=n; ii++) operator()(ii,ii) = init;
}
Matrix::Matrix(const Matrix& init) : vecData(init.vecData), ncols(init.ncols), nrows(init.nrows)
{
}
void Matrix::identity(const unsigned int& n, const double& init) {
resize(n,n,0.);
for(unsigned int ii=1; ii<=n; ii++) operator()(ii,ii) = init;
......@@ -67,26 +46,16 @@ void Matrix::identity(const unsigned int& n, const double& init) {
void Matrix::resize(const unsigned int& rows, const unsigned int& cols) {
clear();
if ((rows > 0) && (cols > 0)) {
vecData.resize(rows*cols);
ncols = cols;
nrows = rows;
} else {
throw IndexOutOfBoundsException("Can not resize a matrix to negative sizes!", AT);
}
vecData.resize(rows*cols);
nrows = rows;
ncols = cols;
}
void Matrix::resize(const unsigned int& rows, const unsigned int& cols, const double& init) {
clear();
if ((rows > 0) && (cols > 0)) {
vecData.resize(rows*cols, init);
ncols = cols;
nrows = rows;
} else {
throw IndexOutOfBoundsException("Can not resize a matrix to negative sizes!", AT);
}
vecData.resize(rows*cols, init);
ncols = cols;
nrows = rows;
}
void Matrix::size(unsigned int& rows, unsigned int& cols) const{
......@@ -156,7 +125,7 @@ bool Matrix::operator==(const Matrix& in) const {
if(nrows!=in_nrows || ncols!=in_ncols)
return false;
for(unsigned int ii=0; ii<vecData.size(); ii++)
for(size_t ii=0; ii<vecData.size(); ii++)
if( !IOUtils::checkEpsilonEquality( vecData[ii] , in.vecData[ii], epsilon_mtr) ) return false;
return true;
......@@ -177,7 +146,7 @@ Matrix& Matrix::operator+=(const Matrix& rhs) {
}
//fill sum matrix
for(unsigned int ii=0; ii<vecData.size(); ii++)
for(size_t ii=0; ii<vecData.size(); ii++)
vecData[ii] += rhs.vecData[ii];
return *this;
......@@ -192,7 +161,7 @@ const Matrix Matrix::operator+(const Matrix& rhs) const {
Matrix& Matrix::operator+=(const double& rhs) {
//fill sum matrix
for(unsigned int ii=0; ii<vecData.size(); ii++)
for(size_t ii=0; ii<vecData.size(); ii++)
vecData[ii] += rhs;
return *this;
......@@ -216,7 +185,7 @@ Matrix& Matrix::operator-=(const Matrix& rhs) {
}
//fill sum matrix
for(unsigned int ii=0; ii<vecData.size(); ii++)
for(size_t ii=0; ii<vecData.size(); ii++)
vecData[ii] -= rhs.vecData[ii];
return *this;
......@@ -278,7 +247,7 @@ const Matrix Matrix::operator*(const Matrix& rhs) const {
}
Matrix& Matrix::operator*=(const double& rhs) {
for(unsigned int ii=0; ii<vecData.size(); ii++)
for(size_t ii=0; ii<vecData.size(); ii++)
vecData[ii] *= rhs;
return *this;
......@@ -473,7 +442,7 @@ Matrix Matrix::getInv() const {
return X;
}
void Matrix::inv() {
bool Matrix::inv() {
//same as getInv() const but we write the final result on top of the input matrix
if(nrows!=ncols) {
std::stringstream tmp;
......@@ -486,7 +455,7 @@ void Matrix::inv() {
Matrix U;
Matrix L;
if(LU(L, U)==false) {
throw IOException("LU decomposition of given matrix not possible", AT);
return false;
}
//we solve AX=I with X=A-1. Since A=LU, then LUX = I
......@@ -494,7 +463,7 @@ void Matrix::inv() {
Matrix Y(n, n);
for(unsigned int i=1; i<=n; i++) {
if(IOUtils::checkEpsilonEquality(L(i,i), 0., epsilon)) {
throw IOException("The given matrix can not be inversed", AT);
return false;
}
Y(i,i) = 1./L(i,i); //j==i
for(unsigned int j=1; j<i; j++) { //j<i
......@@ -512,8 +481,8 @@ void Matrix::inv() {
//now, we backward solve UX=Y
Matrix& X = *this; //we write the solution over the input matrix
for(unsigned int i=n; i>=1; i--) { //lines
if(IOUtils::checkEpsilonEquality(U(i,i), 0., epsilon)) { //HACK: actually, only U(n,n) needs checking
throw IOException("The given matrix is singular and can not be inversed", AT);
if(IOUtils::checkEpsilonEquality(U(i,i), 0., epsilon)) { //actually, only U(n,n) needs checking
return false;
}
for(unsigned int j=1; j<=n; j++) { //lines
double sum=0.;
......@@ -523,9 +492,11 @@ void Matrix::inv() {
X(i,j) = (Y(i,j) - sum) / U(i,i);
}
}
return true;
}
void Matrix::solve(const Matrix& A, const Matrix& B, Matrix& X) {
bool Matrix::solve(const Matrix& A, const Matrix& B, Matrix& X) {
//This uses an LU decomposition followed by backward and forward solving for A·X=B
unsigned int Anrows,Ancols, Bnrows, Bncols;
A.size(Anrows, Ancols);
......@@ -549,7 +520,7 @@ void Matrix::solve(const Matrix& A, const Matrix& B, Matrix& X) {
Matrix U;
Matrix L;
if(A.LU(L, U)==false) {
throw IOException("LU decomposition of A matrix not possible", AT);
return false;
}
//we solve AX=B. Since A=LU, then LUX = B
......@@ -557,7 +528,7 @@ void Matrix::solve(const Matrix& A, const Matrix& B, Matrix& X) {
Matrix Y(n, m);
for(unsigned int i=1; i<=n; i++) {
if(IOUtils::checkEpsilonEquality(L(i,i), 0., epsilon)) {
throw IOException("The given matrix can not be inversed", AT);
return false;
}
for(unsigned int j=1; j<=m; j++) {
double sum=0.;
......@@ -571,10 +542,9 @@ void Matrix::solve(const Matrix& A, const Matrix& B, Matrix& X) {
//now, we backward solve UX=Y
X.resize(n,m); //we need to ensure that X has the correct dimensions
for(unsigned int i=n; i>=1; i--) { //lines
if(IOUtils::checkEpsilonEquality(U(i,i), 0., epsilon)) { //HACK: actually, only U(n,n) needs checking
std::stringstream ss;
ss << "The given " << Ancols << "*" << Anrows << " matrix is singular and can not be inversed";
throw IOException(ss.str(), AT);
if(IOUtils::checkEpsilonEquality(U(i,i), 0., epsilon)) { //actually, only U(n,n) needs checking
//singular matrix
return false;
}
for(unsigned int j=1; j<=m; j++) {
double sum = 0.;
......@@ -584,16 +554,19 @@ void Matrix::solve(const Matrix& A, const Matrix& B, Matrix& X) {
X(i,j) = (Y(i,j) - sum) / U(i,i);
}
}
return true;
}
Matrix Matrix::solve(const Matrix& A, const Matrix& B) {
//This uses an LU decomposition followed by backward and forward solving for A·X=B
Matrix X;
solve(A, B, X);
if(!solve(A, B, X))
throw IOException("Matrix inversion failed!", AT);
return X;
}
void Matrix::TDMA_solve(const Matrix& A, const Matrix& B, Matrix& X)
bool Matrix::TDMA_solve(const Matrix& A, const Matrix& B, Matrix& X)
{ //Thomas algorithm for tridiagonal matrix solving of A·X=B
unsigned int Anrows,Ancols, Bnrows, Bncols;
A.size(Anrows, Ancols);
......@@ -624,7 +597,7 @@ void Matrix::TDMA_solve(const Matrix& A, const Matrix& B, Matrix& X)
b[1] = A(1,1); v[1] = B(1,1); //otherwise they would never be defined
for(unsigned int i=2; i<=n; i++) {
if(IOUtils::checkEpsilonEquality(b[i-1], 0., epsilon))
throw IOException("The given matrix can not be inversed", AT);
return false;
const double b_i = A(i,i);
const double v_i = B(i,1);
const double a_i = A(i,i-1);
......@@ -639,13 +612,17 @@ void Matrix::TDMA_solve(const Matrix& A, const Matrix& B, Matrix& X)
for(unsigned int i=n-1; i>=1; i--) {
X(i,1) = ( v[i] - c[i]*X(i+1,1) ) / b[i];
}
return true;
}
Matrix Matrix::TDMA_solve(const Matrix& A, const Matrix& B) {
//This uses the Thomas algorithm for tridiagonal matrix solving of A·X=B
Matrix X;
TDMA_solve(A, B, X);
return X;
if(TDMA_solve(A, B, X))
return X;
else
throw IOException("Matrix inversion failed!", AT);
}
bool Matrix::isIdentity() const {
......
......@@ -43,7 +43,7 @@ namespace mio {
*/
class Matrix {
public:
Matrix();
Matrix() : vecData(), ncols(0), nrows(0) {};
/**
* @brief A constructor that creates a matrix of a given size
......@@ -51,7 +51,7 @@ class Matrix {
* @param cols number of columns of the new matrix
*/
Matrix(const int& rows, const int& cols);
Matrix(const unsigned int& rows, const unsigned int& cols);
Matrix(const unsigned int& rows, const unsigned int& cols) : vecData(rows*cols), ncols(cols), nrows(rows) {};
/**
* @brief A constructor that creates a matrix filled with constant values
......@@ -59,7 +59,7 @@ class Matrix {
* @param cols number of columns of the new matrix
* @param init initial value to fill the matrix with
*/
Matrix(const unsigned int& rows, const unsigned int& cols, const double& init);
Matrix(const unsigned int& rows, const unsigned int& cols, const double& init) : vecData(rows*cols, init), ncols(cols), nrows(rows) {};
/**
* @brief A constructor that creates a diagonal matrix of size n
......@@ -72,7 +72,7 @@ class Matrix {
* @brief Copy constructor
* @param init matrix to copy
*/
Matrix(const Matrix& init);
Matrix(const Matrix& init) : vecData(init.vecData), ncols(init.ncols), nrows(init.nrows) {};
/**
* @brief Convert the current matrix to a identity matrix of size n
......@@ -134,7 +134,7 @@ class Matrix {
* @return inversed matrix
*/
Matrix getInv() const;
void inv();
bool inv();
/**
* @brief matrix solving for A·X=B.
......@@ -153,8 +153,9 @@ class Matrix {
* @param A A matrix
* @param B B matrix
* @param X solution matrix
* @return true is success
*/
static void solve(const Matrix& A, const Matrix& B, Matrix& X);
static bool solve(const Matrix& A, const Matrix& B, Matrix& X);
/**
* @brief Solving system of equations using Thomas Algorithm
......@@ -177,8 +178,9 @@ class Matrix {
* @param A A matrix
* @param B B matrix
* @param X solution matrix
* @return true is success
*/
static void TDMA_solve(const Matrix& A, const Matrix& B, Matrix& X);
static bool TDMA_solve(const Matrix& A, const Matrix& B, Matrix& X);
/**
* @brief matrix determinant
......
......@@ -45,7 +45,7 @@ Fit1D& Fit1D::operator=(const Fit1D& source) { //HACK: the pointer could not be
return *this;
}
void Fit1D::setModel(const std::string& i_regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit)
bool Fit1D::setModel(const std::string& i_regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit)
{
regression regType;
if(i_regType=="ZERO") regType=ZERO;
......@@ -61,10 +61,10 @@ void Fit1D::setModel(const std::string& i_regType, const std::vector<double>& in
throw IOException("The regression algorithm '"+i_regType+"' is not implemented" , AT);
}
setModel(regType, in_X, in_Y, updatefit);
return setModel(regType, in_X, in_Y, updatefit);
}
void Fit1D::setModel(const regression& regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit)
bool Fit1D::setModel(const regression& regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit)
{
if(model!=NULL) delete model;
......@@ -79,7 +79,10 @@ void Fit1D::setModel(const regression& regType, const std::vector<double>& in_X,
if(regType==QUADRATIC) model=new Quadratic;
model->setData(in_X, in_Y);
if(updatefit) fit();
if(updatefit)
return fit();
else
return true;
}
//////////////////////////////////////////////////////////////
......
......@@ -170,8 +170,9 @@ class Fit1D {
* @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())
* @return false if could not compute the parameters. if !updatefit, always return true
*/
void setModel(const regression& i_regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit=true);
bool setModel(const regression& i_regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit=true);
/**
* @brief Set or reset the regression model.
......@@ -179,8 +180,9 @@ class Fit1D {
* @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())
* @return false if could not compute the parameters. if !updatefit, always return true
*/
void setModel(const std::string& i_regType, const std::vector<double>& in_X, const std::vector<double>& in_Y, const bool& updatefit=true);