WSL/SLF GitLab Repository

Commit 9604a0f1 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A new test has been created to cover basic statistical functions on vectors....

A new test has been created to cover basic statistical functions on vectors. Two bugs have been found and fixed when computing the derivative on a vector of data (the right derivative was in fact a left derivative and the case of X containing nodata was not covered).

Documentation fixes in libinterpol2D and ResamplingAlgorithms. The cmake macro for finding MeteoIO now looks first in the home directory of the user on Mac (as is done on Linux).
parent aec7a348
......@@ -106,7 +106,7 @@ void ResamplingAlgorithms::getNearestValidPts(const size_t& pos, const size_t& p
* @param y1 y-coordinate of first point
* @param x2 x-coordinate of second point
* @param y2 y-coordinate of second point
* @param x3 x-coordinate of desired point
* @param x x-coordinate of desired point
* @return y-coordinate of desired point
*/
double ResamplingAlgorithms::linearInterpolation(const double& x1, const double& y1,
......
......@@ -80,6 +80,11 @@ std::vector<double> Interpol1D::quantiles(const std::vector<double>& X, const st
return vecResults;
}
//small helper function for checking if a point can be used when computing a vector derivative
bool Interpol1D::ptOK(const double& x, const double& y) {
return (x!=IOUtils::nodata && y!=IOUtils::nodata);
}
/**
* @brief This function returns the vector of local derivatives, given a vector of abscissae and ordinates.
* The vectors must be sorted by ascending x. The derivatives will be centered if possible, left or right otherwise or nodata
......@@ -90,28 +95,37 @@ std::vector<double> Interpol1D::quantiles(const std::vector<double>& X, const st
*/
std::vector<double> Interpol1D::derivative(const std::vector<double>& X, const std::vector<double>& Y)
{
const size_t Xsize = X.size();
if(Xsize!=Y.size()) {
const size_t n = X.size();
if(n!=Y.size()) {
stringstream ss;
ss << "X vector and Y vector don't match! " << Xsize << "!=" << Y.size() << "\n";
ss << "X vector and Y vector don't match! " << n << "!=" << Y.size() << "\n";
throw InvalidArgumentException(ss.str(), AT);
}
if(n<2) {
stringstream ss;
ss << "X and Y vector only contain " << n << "points, it is not possible to compute a derivative!\n";
throw InvalidArgumentException(ss.str(), AT);
}
std::vector<double> der(Xsize, IOUtils::nodata);
double right, centered, left, Dx_r, Dx_c, Dx_l;
size_t i=0;
std::vector<double> der(n, IOUtils::nodata);
//right hand derivative
Dx_r=X[i+1]-X[i];
if(Y[i]!=IOUtils::nodata && Y[i+1]!=IOUtils::nodata && Dx_r!=0.) der[i] = (Y[i+1]-Y[i]) / Dx_r;
{
const double x=X[0], x1=X[1];
const double y=Y[0], y1=Y[1];
const double Dx_r=x1-x;
if(ptOK(x,y) && ptOK(x1,y1) && Dx_r!=0.)
der[0] = (y1-y) / Dx_r;
}
//centered derivative if possible
i++;
for(; i<(Xsize-1); i++) {
Dx_r=X[i-1]-X[i]; Dx_c=X[i+1]-X[i-1]; Dx_l=X[i]-X[i-1];
for(size_t i=1; i<(n-1); i++) {
const double x0=X[i-1], x=X[i], x1=X[i+1];
const double y0=Y[i-1], y=Y[i], y1=Y[i+1];
const double Dx_r=x1-x, Dx_c=x1-x0, Dx_l=x-x0;
if(Y[i]!=IOUtils::nodata && Y[i+1]!=IOUtils::nodata && Dx_r!=0.) right=(Y[i+1]-Y[i])/Dx_r; else right=IOUtils::nodata;
if(Y[i]!=IOUtils::nodata && Y[i-1]!=IOUtils::nodata && Dx_l!=0.) left=(Y[i]-Y[i-1])/Dx_l; else left=IOUtils::nodata;
if(Y[i-1]!=IOUtils::nodata && Y[i+1]!=IOUtils::nodata && Dx_c!=0.) centered=(Y[i+1]-Y[i-1])/Dx_c; else centered=IOUtils::nodata;
const double right = (ptOK(x,y) && ptOK(x1,y1) && Dx_r!=0.)? (y1-y) / Dx_r : IOUtils::nodata;
const double left = (ptOK(x,y) && ptOK(x0,y0) && Dx_l!=0.)? (y-y0) / Dx_l : IOUtils::nodata;
const double centered = (ptOK(x0,y0) && ptOK(x1,y1) && Dx_c!=0.)? (y1-y0) / Dx_c : IOUtils::nodata;
if(centered!=IOUtils::nodata) der[i] = centered;
else if(right!=IOUtils::nodata) der[i] = right;
......@@ -119,8 +133,14 @@ std::vector<double> Interpol1D::derivative(const std::vector<double>& X, const s
}
//left hand derivative
Dx_l=X[i]-X[i-1];
if(Y[i]!=IOUtils::nodata && Y[i-1]!=IOUtils::nodata && Dx_l!=0.) der[i] = (Y[i]-Y[i-1]) / Dx_l;
{
const size_t last = n-1;
const double x0=X[last-1], x=X[last];
const double y0=Y[last-1], y=Y[last];
const double Dx_l=x-x0;
if(ptOK(x,y) && ptOK(x0,y0) && Dx_l!=0.)
der[last] = (y-y0) / Dx_l;
}
return der;
}
......
......@@ -56,6 +56,7 @@ class Interpol1D {
static void ExpRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r, std::stringstream& mesg);
private:
static bool ptOK(const double& x, const double& y);
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);
};
......
......@@ -331,7 +331,7 @@ void Interpol2D::IDW(const std::vector<double>& vecData_in, const std::vector<St
* @brief Grid filling function:
* This implementation fills a grid using a curvature and slope algorithm, as described in "A Meteorological
* Distribution System for High-Resolution Terrestrial Modeling (MicroMet)", Liston and Elder, 2006.
* @param dem array of elevations (dem). The slope must have been updated as it is required for the DEM analysis.
* @param i_dem array of elevations (dem). The slope must have been updated as it is required for the DEM analysis.
* @param VW 2D array of Wind Velocity to fill
* @param DW 2D array of Wind Direction to fill
*/
......
......@@ -34,3 +34,4 @@ ADD_SUBDIRECTORY(dem_reading)
ADD_SUBDIRECTORY(2D_interpolations)
ADD_SUBDIRECTORY(arrays)
ADD_SUBDIRECTORY(coords)
ADD_SUBDIRECTORY(stats)
## Test statistical functions
# generate executable
ADD_EXECUTABLE(stats stats.cc)
TARGET_LINK_LIBRARIES(stats ${LIBRARIES})
# add the tests
ADD_TEST(stats.smoke stats)
SET_TESTS_PROPERTIES(stats.smoke PROPERTIES LABELS smoke)
//#include <stdlib.h>
#include <cstdlib>
#include <time.h>
#include <meteoio/MeteoIO.h>
using namespace std;
using namespace mio;
const double rand_range = 1000.;
void cr_fixed_vectors(vector<double> &X, vector<double> &Y) {
X.clear(); X.resize(10);
Y.clear(); Y.resize(10);
X[0] = -499.882; Y[0] = 421.516;
X[1] = -114.998; Y[1] = 184.937;
X[2] = -103.033; Y[2] = 241.216;
X[3] = -86.3743; Y[3] = -140.725;
X[4] = 57.4737; Y[4] = -377.073;
X[5] = 127.179; Y[5] = -376.651;
X[6] = 250.781; Y[6] = 24.0065;
X[7] = -999.; Y[7] = 496.798;
X[8] = 393.442; Y[8] = -999.;
X[9] = 408.396; Y[9] = 105.45;
}
void cr_rand_vectors(vector<double> &X, vector<double> &Y) {
const size_t N = 10;
srand( time(NULL) );
X.clear(); X.resize(N);
Y.clear(); Y.resize(N);
for(size_t ii=0; ii<N; ++ii) {
X[ii] = rand()/(double)RAND_MAX*rand_range - rand_range/2.;
Y[ii] = rand()/(double)RAND_MAX*rand_range - rand_range/2.;
}
}
void print_vector(const vector<double>& X) {
for(size_t ii=0; ii<X.size(); ++ii)
std::cout << "X[" << ii << "] = " << X[ii] << "\n";
}
void print_vectors(const vector<double>& X, const vector<double>& Y) {
for(size_t ii=0; ii<X.size(); ++ii)
std::cout << "X[" << ii << "] = " << X[ii] << "; Y[" << ii << "] = " << Y[ii] << ";\n";
}
////////////////////////////////// start real tests
bool check_sort(const vector<double>& x, const vector<double>& y) {
vector<double> X(x), Y(y);
Interpol1D::sort(X, Y);
bool status = true;
for(size_t ii=1; ii<X.size(); ++ii) {
const double value = X[ii];
if(value < X[ii-1]) status=false;
const size_t index = find(x.begin(), x.end(), value) - x.begin();
if(y[index] != Y[ii]) status=false;
}
if(!status) std::cerr << "Sorting of vectors failed\n";
return status;
}
bool check_quantiles(const vector<double>& X) {
static const double arr[] = {.1, .2, .4, .5, .75, .95};
static const double results[] = {-191.9748, -107.819, -57.6047, 57.4737, 250.781, 402.4144};
const vector<double> quartiles(arr, arr + sizeof(arr) / sizeof(arr[0]));
const vector<double> quantiles = Interpol1D::quantiles(X, quartiles);
bool status = true;
for(size_t ii=0; ii<quantiles.size(); ++ii) {
if(!IOUtils::checkEpsilonEquality(quantiles[ii], results[ii], 1e-6)) {
std::cerr << setprecision(18) << "Quantile[" << ii << "] should be " << results[ii] << ", computed " << quantiles[ii] << " instead\n";
status=false;
}
}
return status;
}
bool check_basics(const vector<double>& X) {
//median
const double median = Interpol1D::getMedian(X);
const double median_ref = 57.4737;
const bool median_status = IOUtils::checkEpsilonEquality(median, median_ref, 1e-6);
if(!median_status)
std::cerr << setprecision(12) << "Median should be " << median_ref << ", computed " << median << " instead\n";
//MAD
const double mad = Interpol1D::getMedianAverageDeviation(X);
const double mad_ref = 172.4717;
const bool mad_status = IOUtils::checkEpsilonEquality(mad, mad_ref, 1e-6);
if(!mad_status)
std::cerr << setprecision(12) << "MAD should be " << mad_ref << ", computed " << mad << " instead\n";
//variance
const double variance = Interpol1D::variance(X);
const double variance_ref = 83038.1246275;
const bool variance_status = IOUtils::checkEpsilonEquality(variance, variance_ref, 1e-6);
if(!variance_status)
std::cerr << setprecision(12) << "Variance should be " << variance_ref << ", computed " << variance << " instead\n";
//stddev
const double stddev = Interpol1D::std_dev(X);
const double stddev_ref = 288.163364478;
const bool stddev_status = IOUtils::checkEpsilonEquality(stddev, stddev_ref, 1e-6);
if(!stddev_status)
std::cerr << setprecision(12) << "Stddev should be " << stddev_ref << ", computed " << stddev << " instead\n";
//weighted mean
const double d1 = 288.1643545;
const double d2 = 384.1562055;
const double w1 = 0.232326, w2 = 0.68125;
const double m1_ref = 310.465757275, m2_ref = 353.558802994;
const double m1 = Interpol1D::weightedMean(d1, d2, w1);
const double m2 = Interpol1D::weightedMean(d1, d2, w2);
double wmean_status = true;
if( !IOUtils::checkEpsilonEquality(m1, m1_ref, 1e-6) || !IOUtils::checkEpsilonEquality(m2, m2_ref, 1e-6)) {
std::cerr << setprecision(12) << "Weighted means should be " << m1_ref << " and " << m2_ref;
std::cerr << ", computed " << m1 << " and " << m2 << " instead\n";
wmean_status = false;
}
//vector weighted mean
const vector<double> weights(X.size(), 1./(X.size()));
const double vector_mean = Interpol1D::weightedMean(X, weights);
const double mean = Interpol1D::arithmeticMean(X);
const double vector_mean_status = IOUtils::checkEpsilonEquality(vector_mean, mean, 1e-6);
if(!vector_mean_status) {
std::cerr << setprecision(12) << "Vector mean should be " << mean << ", conputed " << vector_mean << " instead\n";
}
return (mad_status || variance_status || median_status || stddev_status || wmean_status || vector_mean_status);
}
bool check_covariance(const vector<double>& x, const vector<double>& y) {
const double cov = Interpol1D::covariance(x, y);
const double cov_ref = -35272.1266148;
if(!IOUtils::checkEpsilonEquality(cov, cov_ref, 1e-6)) {
std::cerr << setprecision(12) << "Covariance should be " << cov_ref << ", conputed " << cov << " instead\n";
return false;
}
return true;
}
bool check_derivative(const vector<double>& x, const vector<double>& y) {
static const double results[] = {-999., -0.6146761, -0.454329, -11.377355, -3.8521071, -1.104764, 2.0748285, 3.241513, 0.516724, -999.};
vector<double> X(x), Y(y);
Interpol1D::sort(X, Y);
const vector<double> der = Interpol1D::derivative(X, Y);
bool der_status = true;
for(size_t ii=0; ii<der.size(); ++ii) {
if(!IOUtils::checkEpsilonEquality(der[ii], results[ii], 1e-6)) {
std::cerr << setprecision(12) << "Derivative[" << ii << "] should be " << results[ii] << ", conputed " << der[ii] << " instead\n";
der_status = false;
}
}
return der_status;
}
int main() {
vector<double> x,y;
//cr_rand_vectors(x, y);
cr_fixed_vectors(x, y);
const bool basics_status = check_basics(x);
const bool sort_status = check_sort(x,y);
const bool quantiles_status = check_quantiles(x);
const bool covariance_status = check_covariance(x,y);
const bool der_status = check_derivative(x,y);
if(!basics_status || !sort_status || !quantiles_status || !covariance_status || !der_status)
throw IOException("Statistical functions error!", AT);
return 0;
}
\ No newline at end of file
......@@ -38,10 +38,10 @@ ELSE(WIN32)
FIND_LIBRARY(METEOIO_LIBRARY
NAMES meteoio${POPC_EXT}
PATHS
"/Applications/MeteoIO/lib"
ENV LD_LIBRARY_PATH
ENV DYLD_FALLBACK_LIBRARY_PATH
"~/usr/lib"
"/Applications/MeteoIO/lib"
"/usr/local/lib"
"/usr/lib"
"/opt/lib"
......
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