The MeteoData == operator now compares with an epsilon. The...

The MeteoData == operator now compares with an epsilon. The chekEpsilonEquality function has been rewritten, according to best practise and inlined. The powN functions have been moved to MathOptim and inlined. This makes the benchmarks 2% faster, that's such a victory...

meteoio/Coords.cc

@@ -18,6 +18,7 @@
 #include <cmath>
 
 #include <meteoio/Coords.h>
+#include <meteoio/MathOptim.h>
 #include <meteoio/meteolaws/Meteoconst.h> //for math constants
 
 #ifdef PROJ4
@@ -872,7 +873,7 @@ double Coords::lat_degree_lenght(const double& latitude) {
 	const double b = ellipsoids[E_WGS84].b;	//minor ellipsoid semi-axis
 	const double e2 = (a*a-b*b) / (a*a);	//ellispoid eccentricity, squared
 
-	const double degree_length = (Cst::PI*a*(1.-e2)) / ( 180.*pow(1.-e2*IOUtils::pow2(sin(latitude*Cst::to_rad)), 1.5) );
+	const double degree_length = (Cst::PI*a*(1.-e2)) / ( 180.*pow(1.-e2*Optim::pow2(sin(latitude*Cst::to_rad)), 1.5) );
 	return fabs( degree_length );
 }
 
@@ -888,7 +889,7 @@ double Coords::lon_degree_lenght(const double& latitude) {
 	const double b = ellipsoids[E_WGS84].b;	//minor ellipsoid semi-axis
 	const double e2 = (a*a-b*b) / (a*a);	//ellispoid eccentricity, squared
 
-	const double degree_length = (Cst::PI*a*cos(latitude*Cst::to_rad)) / ( 180.*sqrt(1.-e2*IOUtils::pow2(sin(latitude*Cst::to_rad))) );
+	const double degree_length = (Cst::PI*a*cos(latitude*Cst::to_rad)) / ( 180.*sqrt(1.-e2*Optim::pow2(sin(latitude*Cst::to_rad))) );
 	return fabs( degree_length );
 }
 
@@ -1122,7 +1123,7 @@ void Coords::WGS84_to_UTM(double lat_in, double long_in, double& east_out, doubl
 	const double long0 = (double)((zoneNumber - 1)*6 - 180 + 3) * Cst::to_rad; //+3 puts origin in middle of zone
 
 	//Geometrical parameters
-	const double nu = a / sqrt(1.-e2*IOUtils::pow2(sin(Lat))); //radius of curvature of the earth perpendicular to the meridian plane
+	const double nu = a / sqrt(1.-e2*Optim::pow2(sin(Lat))); //radius of curvature of the earth perpendicular to the meridian plane
 	const double p = (Long-long0);
 
 	//calculating first the coefficients of the series, then the Meridional Arc M itself
@@ -1138,9 +1139,9 @@ void Coords::WGS84_to_UTM(double lat_in, double long_in, double& east_out, doubl
 	//calculating the coefficients for the series
 	const double K1 = M*k0;
 	const double K2 = 1./4.*k0*nu*sin(2.*Lat);
-	const double K3 = (k0*nu*sin(Lat)*IOUtils::pow3(cos(Lat))*1./24.) * (5. - IOUtils::pow2(tan(Lat)) + 9.*eP2*IOUtils::pow2(cos(Lat)) + 4.*eP2*eP2*IOUtils::pow4(cos(Lat)));
+	const double K3 = (k0*nu*sin(Lat)*Optim::pow3(cos(Lat))*1./24.) * (5. - Optim::pow2(tan(Lat)) + 9.*eP2*Optim::pow2(cos(Lat)) + 4.*eP2*eP2*Optim::pow4(cos(Lat)));
 	const double K4 = k0*nu*cos(Lat);
-	const double K5 = (k0*nu*IOUtils::pow3(cos(Lat))*1./6.) * (1. - IOUtils::pow2(tan(Lat)) + eP2*IOUtils::pow2(cos(Lat)));
+	const double K5 = (k0*nu*Optim::pow3(cos(Lat))*1./6.) * (1. - Optim::pow2(tan(Lat)) + eP2*Optim::pow2(cos(Lat)));
 
 	north_out = K1 + K2*p*p + K3*p*p*p*p;
 	east_out = K4*p + K5*p*p*p + 500000.0;
@@ -1193,10 +1194,10 @@ void Coords::UTM_to_WGS84(double east_in, double north_in, double& lat_out, doub
 	const double fp = mu + J1*sin(2.*mu) + J2*sin(4.*mu) + J3*sin(6.*mu) + J4*sin(8.*mu);
 
 	//calculating the parameters
-	const double C1 = eP2 * IOUtils::pow2(cos(fp));
-	const double T1 = IOUtils::pow2( tan(fp) );
-	const double R1 = a*(1.-e2) / pow((1.-e2*IOUtils::pow2(sin(fp))), 1.5);
-	const double N1 = a / sqrt(1.-e2*IOUtils::pow2(sin(fp)));
+	const double C1 = eP2 * Optim::pow2(cos(fp));
+	const double T1 = Optim::pow2( tan(fp) );
+	const double R1 = a*(1.-e2) / pow((1.-e2*Optim::pow2(sin(fp))), 1.5);
+	const double N1 = a / sqrt(1.-e2*Optim::pow2(sin(fp)));
 	const double D = east_in / (N1*k0);
 
 	//calculating the coefficients of the series for latitude and longitude
@@ -1320,7 +1321,7 @@ void Coords::distance(const Coords& destination, double& distance, double& beari
 //HACK: this is the 2D distance, it does not work in 3D!!
 	if(isSameProj(destination)) {
 		//we can use simple cartesian grid arithmetic
-		distance = sqrt( IOUtils::pow2(easting - destination.getEasting()) + IOUtils::pow2(northing - destination.getNorthing()) );
+		distance = sqrt( Optim::pow2(easting - destination.getEasting()) + Optim::pow2(northing - destination.getNorthing()) );
 		bearing = atan2( northing - destination.getNorthing() , easting - destination.getEasting() );
 		bearing = fmod( bearing*Cst::to_deg+360. , 360. );
 	} else {
@@ -1380,7 +1381,7 @@ void Coords::WGS84_to_local(double lat_in, double long_in, double& east_out, dou
 */
 void Coords::local_to_WGS84(double east_in, double north_in, double& lat_out, double& long_out) const
 {
-	const double distance = sqrt( IOUtils::pow2(east_in) + IOUtils::pow2(north_in) );
+	const double distance = sqrt( Optim::pow2(east_in) + Optim::pow2(north_in) );
 	const double bearing = fmod( atan2(east_in, north_in)*Cst::to_deg+360. , 360.);
 
 	if((ref_latitude==IOUtils::nodata) || (ref_longitude==IOUtils::nodata)) {
@@ -1502,12 +1503,12 @@ double Coords::VincentyDistance(const double& lat1, const double& lon1, const do
 	double sin_sigma, cos_sigma, sigma, cos_alpha2, cos_2sigma_m;
 	int n=0;
 	do {
-		sin_sigma = sqrt( IOUtils::pow2(cos(U2)*sin(lambda)) + IOUtils::pow2(cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda)) );
+		sin_sigma = sqrt( Optim::pow2(cos(U2)*sin(lambda)) + Optim::pow2(cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda)) );
 		if(sin_sigma==0.) return 0.; //co-incident points
 		cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(lambda);
 		sigma = atan2(sin_sigma,cos_sigma);
 		const double sin_alpha = cos(U1)*cos(U2)*sin(lambda) / sin_sigma;
-		cos_alpha2 = 1. - IOUtils::pow2(sin_alpha);
+		cos_alpha2 = 1. - Optim::pow2(sin_alpha);
 		if(lat1==0. && lat2==0.) {
 			cos_2sigma_m = 0.;
 		} else {
@@ -1516,7 +1517,7 @@ double Coords::VincentyDistance(const double& lat1, const double& lon1, const do
 		const double C = f/16. * cos_alpha2*(4.+f*(4.-3.*cos_alpha2));
 		lambda_p = lambda;
 		lambda = L + (1.-C)*f*sin_alpha*(
-			sigma + C*sin_sigma*( cos_2sigma_m + C * cos_sigma * (-1.+2.*IOUtils::pow2(cos_2sigma_m)) )
+			sigma + C*sin_sigma*( cos_2sigma_m + C * cos_sigma * (-1.+2.*Optim::pow2(cos_2sigma_m)) )
 			);
 		n++;
 	} while ( (n<n_max) && (fabs(lambda - lambda_p) > thresh) );
@@ -1528,7 +1529,7 @@ double Coords::VincentyDistance(const double& lat1, const double& lon1, const do
 	const double u2 = cos_alpha2 * (a*a - b*b) / (b*b);
 	const double A = 1. + u2/16384. * ( 4096.+u2*(-768.+u2*(320.-175.*u2)) );
 	const double B = u2/1024. * ( 256.+u2*(-128.+u2*(74.-47.*u2)) );
-	const double delta_sigma = B*sin_sigma*( cos_2sigma_m+B/4.*( cos_sigma*(-1.+2.*IOUtils::pow2(cos_2sigma_m)) - B/6.*(cos_2sigma_m*(-3.+4.*IOUtils::pow2(sin_sigma))*(-3.+4.*IOUtils::pow2(cos_2sigma_m))) ) );
+	const double delta_sigma = B*sin_sigma*( cos_2sigma_m+B/4.*( cos_sigma*(-1.+2.*Optim::pow2(cos_2sigma_m)) - B/6.*(cos_2sigma_m*(-3.+4.*Optim::pow2(sin_sigma))*(-3.+4.*Optim::pow2(cos_2sigma_m))) ) );
 
 	const double s = b*A*(sigma - delta_sigma);	//distance between the two points
 
@@ -1584,14 +1585,14 @@ void Coords::VincentyInverse(const double& lat_ref, const double& lon_ref, const
 		cos2sigma_m = cos( 2.*sigma1 + sigma );
 		double delta_sigma = B*sin(sigma) * ( cos2sigma_m + B/4. * (
 			cos(sigma)*(-1.+2.*cos2sigma_m*cos2sigma_m)
-			-B/6. * cos2sigma_m * (-3.+4.*IOUtils::pow2(sin(sigma))) * (-3.+4.*cos2sigma_m*cos2sigma_m)
+			-B/6. * cos2sigma_m * (-3.+4.*Optim::pow2(sin(sigma))) * (-3.+4.*cos2sigma_m*cos2sigma_m)
 			) );
 		sigma_p = sigma;
 		sigma = distance / (b*A) + delta_sigma;
 	}
 
 	lat = atan2( sinU1*cos(sigma) + cosU1*sin(sigma)*cos(alpha1),
-		     (1.-f) * sqrt( sinAlpha*sinAlpha + IOUtils::pow2(sinU1*sin(sigma) - cosU1*cos(sigma)*cos(alpha1)) )
+		     (1.-f) * sqrt( sinAlpha*sinAlpha + Optim::pow2(sinU1*sin(sigma) - cosU1*cos(sigma)*cos(alpha1)) )
 		   );
 	const double lambda = atan2( sin(sigma)*sin(alpha1), cosU1*cos(sigma) - sinU1*sin(sigma)*cos(alpha1) );
 	const double C = f/16. * cos2alpha * (4.+f*(4.-3.*cos2alpha));

meteoio/DEMObject.cc

@@ -19,6 +19,7 @@
 #include <limits.h>
 
 #include <meteoio/DEMObject.h>
+#include <meteoio/MathOptim.h>
 #include <meteoio/meteolaws/Meteoconst.h> //for math constants
 
 /**
@@ -370,7 +371,7 @@ void DEMObject::sanitize() {
 */
 double DEMObject::horizontalDistance(const double& xcoord1, const double& ycoord1, const double& xcoord2, const double& ycoord2)
 {
-	return sqrt( IOUtils::pow2(xcoord2-xcoord1) + IOUtils::pow2(ycoord2-ycoord1) );
+	return sqrt( Optim::pow2(xcoord2-xcoord1) + Optim::pow2(ycoord2-ycoord1) );
 }
 
 /**
@@ -419,9 +420,9 @@ double DEMObject::terrainDistance(Coords point1, const Coords& point2) {
 				//distance += sqrt( pow2((ix2-ix1)*cellsize) + pow2((iy2-iy1)*cellsize) + pow2(grid2D(ix2,iy2)-grid2D(ix1,iy1)) );
 				const double z1=grid2D(ix1,iy1);
 				const double z2=grid2D(ix2,iy2);
-				const double tmpx=IOUtils::pow2((double)(ix2-ix1)*cellsize);
-				const double tmpy=IOUtils::pow2((double)(iy2-iy1)*cellsize);
-				const double tmpz=IOUtils::pow2(z2-z1);
+				const double tmpx=Optim::pow2((double)(ix2-ix1)*cellsize);
+				const double tmpy=Optim::pow2((double)(iy2-iy1)*cellsize);
+				const double tmpz=Optim::pow2(z2-z1);
 				distance += sqrt( tmpx + tmpy + tmpz );
 			}
 			last_point = ii;

meteoio/Date.cc

@@ -16,6 +16,7 @@
     along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
 */
 #include <meteoio/Date.h>
+#include <meteoio/IOUtils.h>
 #include <cmath>
 
 using namespace std;
@@ -633,7 +634,7 @@ bool Date::operator==(const Date& indate) const {
 		return( undef==true && indate.isUndef() );
 	}
 
-	return( fabs(indate.gmt_julian - gmt_julian) < epsilon );
+	return IOUtils::checkEpsilonEquality(gmt_julian, indate.gmt_julian, epsilon);
 }
 
 bool Date::operator!=(const Date& indate) const {

meteoio/IOUtils.cc

@@ -59,30 +59,6 @@ std::string getLibVersion() {
 	return ss.str();
 }
 
-bool IOUtils::checkEpsilonEquality(const double& val1, const double& val2, const double& epsilon)
-{
-	if (((val1-epsilon) < val2) && ((val1+epsilon) > val2)) {
-		return true;
-	}
-
-	return false;
-}
-
-double IOUtils::pow2(const double& val)
-{
-	return (val*val);
-}
-
-double IOUtils::pow3(const double& val)
-{
-	return (val*val*val);
-}
-
-double IOUtils::pow4(const double& val)
-{
-	return (val*val*val*val);
-}
-
 double IOUtils::bearing_to_angle(const double& bearing) {
 	return (fmod(360.-bearing+90., 360.)*Cst::to_rad);
 }

meteoio/IOUtils.h

@@ -34,6 +34,7 @@
 #include <fstream>
 #include <cctype>
 #include <limits>
+#include <cmath>
 
 #ifndef MAX
 #define MAX(x,y)    (((x) < (y)) ? (y) : (x))
@@ -86,11 +87,7 @@ namespace IOUtils {
 	* @param epsilon is a radius around val1
 	* @return true if val2 is within the radius around val1, false otherwise.
 	*/
-	bool checkEpsilonEquality(const double& val1, const double& val2, const double& epsilon);
-
-	double pow2(const double& val);
-	double pow3(const double& val);
-	double pow4(const double& val);
+	inline bool checkEpsilonEquality(const double& val1, const double& val2, const double& epsilon) {return (fabs(val1-val2) < epsilon);}
 
 	size_t seek(const Date& soughtdate, const std::vector<MeteoData>& vecM, const bool& exactmatch=true);
 

meteoio/InterpolationAlgorithms.cc

@@ -19,6 +19,7 @@
 #include <meteoio/meteolaws/Atmosphere.h>
 #include <meteoio/Meteo2DInterpolator.h>
 #include <meteoio/IOManager.h>
+#include <meteoio/MathOptim.h>
 
 using namespace std;
 
@@ -286,7 +287,7 @@ void ConstLapseRateAlgorithm::calculate(Grid2DObject& grid)
 	}
 
 	//run algorithm
-	info << "r^2=" << IOUtils::pow2( vecCoefficients[3] );
+	info << "r^2=" << Optim::pow2( vecCoefficients[3] );
 	Interpol2D::constantLapseGrid2DFill(avgData, avgAltitudes, dem, vecCoefficients, funcptr, grid);
 }
 
@@ -373,7 +374,7 @@ void IDWLapseAlgorithm::calculate(Grid2DObject& grid)
 	}
 
 	//run algorithm
-	info << "r^2=" << IOUtils::pow2( vecCoefficients[3] );
+	info << "r^2=" << Optim::pow2( vecCoefficients[3] );
 	Interpol2D::LapseIDW(vecData, vecMeta, dem, vecCoefficients, funcptr, grid);
 }
 
@@ -412,7 +413,7 @@ void LocalIDWLapseAlgorithm::calculate(Grid2DObject& grid)
 	//run algorithm
 	double r2=0.;
 	Interpol2D::LocalLapseIDW(vecData, vecMeta, dem, nrOfNeighbors, grid, r2);
-	info << "r^2=" << IOUtils::pow2( r2 );
+	info << "r^2=" << Optim::pow2( r2 );
 }
 
 
@@ -555,7 +556,7 @@ void ILWRAlgorithm::calculate(Grid2DObject& grid)
 	}
 
 	//run algorithm
-	info << "r^2=" << IOUtils::pow2( vecCoefficients[3] );
+	info << "r^2=" << Optim::pow2( vecCoefficients[3] );
 
 	//run algorithm
 	Interpol2D::LapseIDW(vecDataEA, vecMeta, dem, vecCoefficients, funcptr, grid);

meteoio/MathOptim.h

@@ -113,6 +113,9 @@ namespace Optim {
 		return x * invSqrt(x);
 	}
 
+	inline double pow2(const double& val) {return (val*val);}
+	inline double pow3(const double& val) {return (val*val*val);}
+	inline double pow4(const double& val) {return (val*val*val*val);}
 }
 
 } //end namespace

meteoio/MeteoData.cc

@@ -18,6 +18,9 @@
 #include <meteoio/MeteoData.h>
 #include <meteoio/StationData.h>
 
+#include <cmath>
+#include <limits>
+
 using namespace std;
 namespace mio {
 
@@ -213,8 +216,9 @@ bool MeteoData::operator==(const MeteoData& in) const
 	if (nrOfAllParameters != in.nrOfAllParameters) //the number of meteo parameters has to be consistent
 		return false;
 
-	for (size_t ii=0; ii<nrOfAllParameters; ii++){
-		if (data[ii] != in.data[ii])
+	for (size_t ii=0; ii<nrOfAllParameters; ii++) {
+		const double epsilon = (fabs(data[ii]) < fabs(in.data[ii]) ? fabs(in.data[ii]) : fabs(data[ii])) * std::numeric_limits<double>::epsilon();
+		if( !IOUtils::checkEpsilonEquality(data[ii], in.data[ii], epsilon) )
 			return false;
 	}
 

meteoio/plugins/ARPSIO.cc

@@ -16,6 +16,8 @@
     along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
 */
 #include <meteoio/meteolaws/Meteoconst.h> //for PI
+#include <meteoio/MathOptim.h>
+
 #include <string.h>
 #include <algorithm>
 
@@ -147,7 +149,7 @@ void ARPSIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& pa
 		readGridLayer("v", 2, V);
 		for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
 			for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
-				grid_out(ii,jj) = sqrt( IOUtils::pow2(grid_out(ii,jj)) + IOUtils::pow2(V(ii,jj)) );
+				grid_out(ii,jj) = sqrt( Optim::pow2(grid_out(ii,jj)) + Optim::pow2(V(ii,jj)) );
 			}
 		}
 	}

meteoio/plugins/GRIBIO.cc

@@ -20,6 +20,7 @@
 #include <meteoio/meteolaws/Atmosphere.h>
 #include <meteoio/meteolaws/Meteoconst.h> //for PI
 #include <meteoio/DEMObject.h>
+#include <meteoio/MathOptim.h>
 
 #include <cmath>
 #include <errno.h>
@@ -532,7 +533,7 @@ void GRIBIO::readWind(const std::string& filename, const Date& date)
 		DW.set(U.ncols, U.nrows, U.cellsize, U.llcorner);
 		for(unsigned int jj=0; jj<VW.nrows; jj++) {
 			for(unsigned int ii=0; ii<VW.ncols; ii++) {
-				VW(ii,jj) = sqrt( IOUtils::pow2(U(ii,jj)) + IOUtils::pow2(V(ii,jj)) );
+				VW(ii,jj) = sqrt( Optim::pow2(U(ii,jj)) + Optim::pow2(V(ii,jj)) );
 				DW(ii,jj) = fmod( atan2( U(ii,jj), V(ii,jj) ) * Cst::to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
 			}
 		}
@@ -1034,7 +1035,7 @@ void GRIBIO::readMeteoStep(std::vector<StationData> &stations, double *lats, dou
 		if(readMeteoValues(34.2, 105, 10, i_date, npoints, lats, lons, values) //V_10M
 		   && readMeteoValues(33.2, 105, 10, i_date, npoints, lats, lons, values2)) { //U_10M
 			for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
-				Meteo[ii](MeteoData::VW) =  sqrt( IOUtils::pow2(values[ii]) + IOUtils::pow2(values2[ii]) );
+				Meteo[ii](MeteoData::VW) =  sqrt( Optim::pow2(values[ii]) + Optim::pow2(values2[ii]) );
 			}
 		}
 	}