Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/dev-PR137' into dev-tmp
Browse files Browse the repository at this point in the history
  • Loading branch information
@Chrismarsh
Chrismarsh committed Jan 16, 2024
2 parents 57e059b + 5b472b6 commit 38d8dd6
Show file tree
Hide file tree
Showing 6 changed files with 72 additions and 84 deletions.
25 changes: 25 additions & 0 deletions src/interpolation/TPSBasis.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
#pragma once
#include <boost/math/special_functions/expint.hpp>
#include <gsl/gsl_sf_expint.h>
#include <cmath>

/* Basis function used for building a thin plate spline. boost::math::expint is used here instead
* of gsl_sf_expint because it automatically works with Boost's automatic differentiation */
template <typename T>
T TPSBasis(T x){
T gamma = 0.5772156649015328606;
// gsl_sf_expint_E1(x) < epsilon_double and ln(x) > 1 for x > 32
// so -(log(x) + c + gsl_sf_expint_E1(x)) == -(log(x) + c) is true for x > 32
return (x <= 32) ? -(log(x) + gamma + boost::math::expint(1,x)) : -(log(x) + gamma);
}


/* TODO OLD_TPSBasis exists solely for timing FunC's LUTs against direct evaluation. It can be safely removed.
* The following two lines are needed if x>800 for OLD_TPSBasis because gsl_sf_expint_E1 overflows for x>800;
* however, old_error_handler is set to gsl_set_error_handler_off() somewhere else in CHM */
//#include <gsl/gsl_errno.h>
//static gsl_error_handler_t *old_error_handler=gsl_set_error_handler_off();
template <typename T>
T OLD_TPSBasis(T x){
return -(log(x) + 0.577215 + gsl_sf_expint_E1(x));
}
43 changes: 31 additions & 12 deletions src/interpolation/TPSpline.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,31 @@

#include "TPSpline.hpp"

//#include <iostream>
//#define FUNC_DEBUG
#include <func/func.hpp>
#include "TPSBasis.hpp"

// Build FunC lookup table for -(log(x)+gamma+gsl_sf_expint_E1(x))

// hardcoded error of 1e-8
//static func::FailureProofTable<func::UniformEqSpaceInterpTable<3,double>,double> TPSBasis_LUT({TPSBasis<double>}, {1e-5, 32, 0.11111});
static func::FailureProofTable<func::NonUniformEqSpaceInterpTable<3,double>,double> TPSBasis_LUT({FUNC_SET_F(TPSBasis,double)}, {1e-5, 32, 0.154207});
//static func::DirectEvaluation<double> TPSBasis_LUT({OLD_TPSBasis<double>}, 1e-8, 6500);


/* Each MPI process will print that text from the destructor after it finishes */
//static struct Counter{
// double total = 0.0;
// ~Counter(){ std::cerr << "Spent " << total << " seconds in spline::operator()" << std::endl; }
//} counter;


double thin_plate_spline::operator()(std::vector< boost::tuple<double,double,double> >& sample_points, boost::tuple<double,double,double>& query_point)
{
//func::Timer t;
//see if we can reuse our
if(sample_points.size() +1 != size)
if(sample_points.size() + 1 != size)
{
size = sample_points.size();
size++; // need to make room for the physics
Expand Down Expand Up @@ -78,8 +99,8 @@ double thin_plate_spline::operator()(std::vector< boost::tuple<double,double,dou
//And Hengl and Evans in geomorphometry p.52 do not, but have some undefined omega_0/omega_1 weights
//it is all rather confusing. But this follows Mitášová exactly, and produces essentially the same answer
//as the worked example in box 16.2 in Chang
Rd = -(log(dij) + c + gsl_sf_expint_E1(dij));

// set Rd = -(log(dij) + c + gsl_sf_expint_E1(dij))
Rd = TPSBasis_LUT(dij);
}

A(i, j + 1) = Rd;
Expand Down Expand Up @@ -132,18 +153,21 @@ double thin_plate_spline::operator()(std::vector< boost::tuple<double,double,dou
double ydiff = (sy - ey);
double dij = sqrt(xdiff*xdiff + ydiff*ydiff);
dij = (dij * weight/2.0) * (dij * weight/2.0);
double Rd = -(log(dij) + c + gsl_sf_expint_E1(dij));
// set Rd equal to -(log(dij) + c + gsl_sf_expint_E1(dij))

double Rd = TPSBasis_LUT(dij);

z0 = z0 + x(i)*Rd;
}

//t.stop();
//counter.total += t.duration();
return z0;
}

thin_plate_spline::thin_plate_spline(size_t sz, std::map<std::string,std::string> config )
: thin_plate_spline()
{

size = sz;
size++; // need to make room for the physics

Expand All @@ -157,10 +181,9 @@ thin_plate_spline::thin_plate_spline(size_t sz, std::map<std::string,std::string
if(config["reuse_LU"] == "true")
reuse_LU = true;
}

reuse_LU = false;

}

thin_plate_spline::thin_plate_spline()
{
pi = 3.14159;
Expand All @@ -170,10 +193,6 @@ thin_plate_spline::thin_plate_spline()

reuse_LU = false;
uninit_lu_decomp = true;


}
thin_plate_spline::~thin_plate_spline()
{

}
thin_plate_spline::~thin_plate_spline(){}
1 change: 0 additions & 1 deletion src/interpolation/TPSpline.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,5 +75,4 @@ class thin_plate_spline : public interp_base
double weight;
bool uninit_lu_decomp;
size_t size;

};
24 changes: 4 additions & 20 deletions third_party/snobal/chm_sati_function.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,24 +1,8 @@
#pragma once
#include <func/EvaluationFunctor.hpp>
#include <cmath>
#define FREEZE 2.7316e2
#define FUNC_I(x) (pow(1.e1, -9.09718 * ((FREEZE / x) - 1.) - 3.56654 * log(FREEZE / x) / log(1.e1) + 8.76793e-1 * (1. - (x / FREEZE)) + log(6.1071) / log(1.e1)) * 100.0)
#define FUNCNAME_I "pow(1.e1, -9.09718 * ((FREEZE / x) - 1.) - 3.56654 * log(FREEZE / x) / log(1.e1) + 8.76793e-1 * (1. - (x / FREEZE)) + log(6.1071) / log(1.e1)) * 100.0"
#define SATI_MACRO(x) (pow(1.e1, -9.09718 * ((FREEZE / x) - 1.) - 3.56654 * log(FREEZE / x) / log(1.e1) + 8.76793e-1 * (1. - (x / FREEZE)) + log(6.1071) / log(1.e1)) * 100.0)
#define FUNCNAME_I "sati"

class MyFunction_sati final : public EvaluationFunctor<double,double>
{
public:
double operator()(double x) override { return FUNC_I(x); }
double deriv(double x) override
{
return 100.0 * log(1.e1) * pow(1.e1, -9.09718 * ((FREEZE / x) - 1.) - 3.56654 * log(FREEZE / x) / log(1.e1) + 8.76793e-1 * (1. - (x / FREEZE)) + log(6.1071) / log(1.e1)) * (9.09718 * FREEZE / x / x + 3.56654 / x / log(1.e1) - 8.76793e-1 / FREEZE);
}
double deriv2(double x) override
{
return 100.0 * log(1.e1) * pow(1.e1, -9.09718 * ((FREEZE / x) - 1.) - 3.56654 * log(FREEZE / x) / log(1.e1) + 8.76793e-1 * (1. - (x / FREEZE)) + log(6.1071) / log(1.e1)) * ((9.09718 * FREEZE / x / x + 3.56654 / x / log(1.e1) - 8.76793e-1 / FREEZE) * (9.09718 * FREEZE / x / x + 3.56654 / x / log(1.e1) - 8.76793e-1 / FREEZE) * log(1.e1) + (-18.19436 * FREEZE / x / x / x - 3.56654 / x / x / log(1.e1)));
}
double deriv3(double x) override
{
return 100.0 * log(1.e1) * pow(1.e1, -9.09718 * ((FREEZE / x) - 1.) - 3.56654 * log(FREEZE / x) / log(1.e1) + 8.76793e-1 * (1. - (x / FREEZE)) + log(6.1071) / log(1.e1)) * ((9.09718 * FREEZE / x / x + 3.56654 / x / log(1.e1) - 8.76793e-1 / FREEZE) * (9.09718 * FREEZE / x / x + 3.56654 / x / log(1.e1) - 8.76793e-1 / FREEZE) * (9.09718 * FREEZE / x / x + 3.56654 / x / log(1.e1) - 8.76793e-1 / FREEZE) * log(1.e1) * log(1.e1) + 3.0 * log(1.e1) * (9.09718 * FREEZE / x / x + 3.56654 / x / log(1.e1) - 8.76793e-1 / FREEZE) * (-18.19436 * FREEZE / x / x / x - 3.56654 / x / x / log(1.e1)) + (54.58308 * FREEZE / x / x / x / x + 7.13308 / x / x / x / log(1.e1)));
}
};
template <typename T>
T MyFunction_sati(T x){ return SATI_MACRO(x); }
24 changes: 4 additions & 20 deletions third_party/snobal/chm_satw_function.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,25 +1,9 @@
#pragma once
#include <func/EvaluationFunctor.hpp>
#include <cmath>
#define SEA_LEVEL 1.013246e5
#define BOIL 3.7315e2
#define FUNC_W(x) (pow(1.e1, -7.90298 * (BOIL / x - 1.) + 5.02808 * log(BOIL / x) / log(1.e1) - 1.3816e-7 * (pow(1.e1, 1.1344e1 * (1. - x / BOIL)) - 1.) + 8.1328e-3 * (pow(1.e1, -3.49149 * (BOIL / x - 1.)) - 1.) + log(SEA_LEVEL) / log(1.e1)))
#define FUNCNAME_W "pow(1.e1, -7.90298 * (BOIL / x - 1.) + 5.02808 * log(BOIL / x) / log(1.e1) - 1.3816e-7 * (pow(1.e1, 1.1344e1 * (1. - x / BOIL)) - 1.) + 8.1328e-3 * (pow(1.e1, -3.49149 * (BOIL / x - 1.)) - 1.) + log(SEA_LEVEL) / log(1.e1))"
#define SATW_MACRO(x) (pow(1.e1, -7.90298 * (BOIL / x - 1.) + 5.02808 * log(BOIL / x) / log(1.e1) - 1.3816e-7 * (pow(1.e1, 1.1344e1 * (1. - x / BOIL)) - 1.) + 8.1328e-3 * (pow(1.e1, -3.49149 * (BOIL / x - 1.)) - 1.) + log(SEA_LEVEL) / log(1.e1)))
#define FUNCNAME_W "satw"

class MyFunction_satw final : public EvaluationFunctor<double,double>
{
public:
double operator()(double x) override { return FUNC_W(x); }
double deriv(double x) override
{
return log(1.e1) * pow(1.e1, -7.90298 * (BOIL / x - 1.) + 5.02808 * log(BOIL / x) / log(1.e1) - 1.3816e-7 * (pow(1.e1, 1.1344e1 * (1. - x / BOIL)) - 1.) + 8.1328e-3 * (pow(1.e1, -3.49149 * (BOIL / x - 1.)) - 1.) + log(SEA_LEVEL) / log(1.e1)) * (7.90298 * BOIL / x / x - 5.02808 / log(1.e1) / x - 1.3816e-7 * log(1.e1)* pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x));
}
double deriv2(double x) override
{
return log(1.e1) * pow(1.e1, -7.90298 * (BOIL / x - 1.) + 5.02808 * log(BOIL / x) / log(1.e1) - 1.3816e-7 * (pow(1.e1, 1.1344e1 * (1. - x / BOIL)) - 1.) + 8.1328e-3 * (pow(1.e1, -3.49149 * (BOIL / x - 1.)) - 1.) + log(SEA_LEVEL) / log(1.e1)) * (log(1.e1) * (7.90298 * BOIL / x / x - 5.02808 / log(1.e1) / x - 1.3816e-7 * log(1.e1)* pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x)) * (7.90298 * BOIL / x / x - 5.02808 / log(1.e1) / x - 1.3816e-7 * log(1.e1)* pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x)) + (-15.80596 * BOIL / x / x / x + 5.02808 / log(1.e1) / x / x - 1.3816e-7 * log(1.e1) * log(1.e1) * pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x) * (3.49149 * BOIL / x / x) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (-6.98298 * BOIL / x / x / x)));
}
double deriv3(double x) override
{
return log(1.e1) * pow(1.e1, -7.90298 * (BOIL / x - 1.) + 5.02808 * log(BOIL / x) / log(1.e1) - 1.3816e-7 * (pow(1.e1, 1.1344e1 * (1. - x / BOIL)) - 1.) + 8.1328e-3 * (pow(1.e1, -3.49149 * (BOIL / x - 1.)) - 1.) + log(SEA_LEVEL) / log(1.e1)) * (log(1.e1) * log(1.e1) * (7.90298 * BOIL / x / x - 5.02808 / log(1.e1) / x - 1.3816e-7 * log(1.e1)* pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x)) * (7.90298 * BOIL / x / x - 5.02808 / log(1.e1) / x - 1.3816e-7 * log(1.e1)* pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x)) * (7.90298 * BOIL / x / x - 5.02808 / log(1.e1) / x - 1.3816e-7 * log(1.e1)* pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x)) + 3 * log(1.e1) * (7.90298 * BOIL / x / x - 5.02808 / log(1.e1) / x - 1.3816e-7 * log(1.e1)* pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x)) * (-15.80596 * BOIL / x / x / x + 5.02808 / log(1.e1) / x / x - 1.3816e-7 * log(1.e1) * log(1.e1) * pow(1.e1, 1.1344e1 * (1. - x / BOIL)) * (-1.1344e1 / BOIL) * (-1.1344e1 / BOIL) + 8.1328e-3 * log(1.e1) * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (3.49149 * BOIL / x / x) * (3.49149 * BOIL / x / x) + 8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (-6.98298 * BOIL / x / x / x)) + (8.1328e-3 * log(1.e1) * pow(1.e1, -3.49149 * (BOIL / x - 1.)) * (log(1.e1) * log(1.e1) * (3.49149 * BOIL / x / x) * (3.49149 * BOIL / x / x) * (3.49149 * BOIL / x / x) + 3 * log(1.e1) * (3.49149 * BOIL / x / x) * (-6.98298 * BOIL / x / x / x) + (20.94894 * BOIL / x / x / x / x))));
}
};
template <typename T>
T MyFunction_satw(T x){ return SATW_MACRO(x); }
39 changes: 8 additions & 31 deletions third_party/snobal/sno.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,20 +26,14 @@
//

#include "sno.h"
//#include <func/func.hpp>
//#include "chm_satw_function.hpp"
//#include "chm_sati_function.hpp"
#include <func/func.hpp>
#include "chm_satw_function.hpp"
#include "chm_sati_function.hpp"
using namespace snobalMacros;

//MyFunction_satw func_w;
//MyFunction_sati func_i;
//
//UniformLookupTableGenerator gen_w(&func_w, FREEZE - 50.0, FREEZE + 50.0);
//UniformLookupTableGenerator gen_i(&func_i, 90.0, FREEZE);
//std::string implName = "UniformLinearInterpolationTable";
//double tableTol = 1e-8;
//static std::unique_ptr<EvaluationImplementation> impl_w = gen_w.generate_by_tol(implName, tableTol);
//static std::unique_ptr<EvaluationImplementation> impl_i = gen_i.generate_by_tol(implName, tableTol);
// These hard coded stepsizes satisfy a tolerance of 1e-8. LookupTableParameters takes {min, max, step}
static func::UniformLinearRawInterpTable<double> satw_lut({MyFunction_satw<double>}, {FREEZE - 50.0, FREEZE + 50.0, 0.00258598});
static func::UniformLinearRawInterpTable<double> sati_lut({MyFunction_sati<double>}, {90.0, FREEZE, 0.000388715});

double sno::ssxfr(
double k1, /* layer 1's thermal conductivity (J / (m K sec)) */
Expand All @@ -66,17 +60,7 @@ double sno::satw(double tk) /* air temperature (K) */
BOOST_THROW_EXCEPTION(module_error() << errstr_info ("tk < 0"));
}

l10 = log(1.e1);

x = -7.90298 * (BOIL / tk - 1.) + 5.02808 * log(BOIL / tk) / l10 -
1.3816e-7 * (pow(1.e1, 1.1344e1 * (1. - tk / BOIL)) - 1.) +
8.1328e-3 * (pow(1.e1, -3.49149 * (BOIL / tk - 1.)) - 1.) +
log(SEA_LEVEL) / l10;

x = pow(1.e1, x);

return (x);
// return ((*impl_w)(tk));
return (satw_lut)(tk);
}

double sno::sati(double tk) /* air temperature (K) */
Expand All @@ -95,19 +79,12 @@ double sno::sati(double tk) /* air temperature (K) */
return (x);
}

l10 = log(1.e1);

x = pow(1.e1, -9.09718 * ((FREEZE / tk) - 1.) - 3.56654 * log(FREEZE / tk) / l10 +
8.76793e-1 * (1. - (tk / FREEZE)) + log(6.1071) / l10);


return (x * 1.e2);
if (tk < 90.0)
{
return 0.0;
}

// return ((*impl_i)(tk));
return (sati_lut)(tk);
}

double sno::new_tsno(
Expand Down

0 comments on commit 38d8dd6

Please sign in to comment.