You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
389 lines
13 KiB
389 lines
13 KiB
#ifndef ALGOIM_INTERVAL_HPP
|
|
#define ALGOIM_INTERVAL_HPP
|
|
|
|
// algoim::Interval<N> methods which compute first-order Taylor series with bounded remainder.
|
|
// These methods are based on those described in the paper
|
|
// R. I. Saye, High-Order Quadrature Methods for Implicitly Defined Surfaces and Volumes
|
|
// in Hyperrectangles,SIAM Journal on Scientific Computing, 37(2), A993-A1019 (2015),
|
|
// http://dx.doi.org/10.1137/140966290
|
|
|
|
#include "real.hpp"
|
|
#include "uvector.hpp"
|
|
#include "utility.hpp"
|
|
|
|
namespace algoim
|
|
{
|
|
/* Interval arithmetic using a first-order Taylor series with remainder. A function's range of attainable values is
|
|
evaluated as
|
|
f(x_c + y) = alpha + beta.y + [-eps,eps]
|
|
where alpha is the value of the function at the centre x_c of an interval [x_c - delta, x_c + delta],
|
|
beta is the n-dimensional gradient of the function evaluated at the centre, y is a placeholder for first-order
|
|
variations, and eps bounds the remainder term.
|
|
|
|
delta is a static (thread-local) variable and should be appropriately initialised before a sequence of calculations
|
|
involving any Interval<N> objects. Thus, all intervals share the same delta bounds. This is thread-safe, however
|
|
note this technique is not appropriate when using a complicated composition of algorithms (e.g., if a method to
|
|
calculate a function's value internally depends on another method that applies separate interval arithmetic
|
|
calculations).
|
|
|
|
For a general C^2 function f, a first order Taylor series with remainder can be computed via
|
|
f(alpha + beta.y + h(y)) = f(alpha) + f'(alpha)*(beta.y + h(y)) + 0.5 f''(xi) (beta.y + h(y))^2
|
|
where xi is somewhere in the interval [alpha, alpha + beta.y + h(y)]. Thus, if C bounds |f''| on the same interval,
|
|
then one can deduce
|
|
f(alpha + beta.y + h(y)) = f(alpha) + f'(alpha)*beta.y + hh(y)
|
|
where
|
|
hh(y) = f'(alpha) h(y) + 0.5 f''(xi) (beta.y + h(y))^2
|
|
and
|
|
|hh(y)| <= |f'(alpha)| eps + 0.5 C (|beta|.delta + eps)^2.
|
|
*/
|
|
template<int N>
|
|
struct Interval
|
|
{
|
|
real alpha;
|
|
uvector<real,N> beta;
|
|
real eps;
|
|
|
|
// Default constructor corresponds to the zero function
|
|
Interval() : alpha(0.0), beta(0.0), eps(0.0) {}
|
|
|
|
// Constructor for a constant
|
|
explicit Interval(real alpha) : alpha(alpha), beta(0.0), eps(0.0) {}
|
|
|
|
// Constructor for the linear function f(y) = alpha + beta.y
|
|
explicit Interval(real alpha, const uvector<real,N>& beta) : alpha(alpha), beta(beta), eps(0.0) {}
|
|
|
|
// Constructor with all parameters given
|
|
explicit Interval(real alpha, const uvector<real,N>& beta, real eps) : alpha(alpha), beta(beta), eps(eps) {}
|
|
|
|
// Assignment from a constant
|
|
Interval& operator=(real alpha)
|
|
{
|
|
this->alpha = alpha;
|
|
beta = 0.0;
|
|
eps = 0.0;
|
|
return *this;
|
|
}
|
|
|
|
// Assigment from another Interval, i.e., copy constructor
|
|
Interval& operator=(const Interval& other)
|
|
{
|
|
// Shouldn't need this, but if not done manually (i.e., cannot use "= default;"), Intel compiler
|
|
// icpc 15 has uninitialised memory problem
|
|
alpha = other.alpha;
|
|
beta = other.beta;
|
|
eps = other.eps;
|
|
return *this;
|
|
}
|
|
|
|
// Maximum deviation of the interval's values from its value at the centre of the interval, i.e.,
|
|
// place bounds on the linear plus remainder term.
|
|
real maxDeviation() const
|
|
{
|
|
using std::abs;
|
|
real b = eps;
|
|
for (int dim = 0; dim < N; ++dim)
|
|
b += abs(beta(dim)) * delta(dim);
|
|
return b;
|
|
}
|
|
|
|
// Unary negation operator
|
|
Interval operator-() const
|
|
{
|
|
return Interval(-alpha, -beta, eps);
|
|
}
|
|
|
|
// Addition by a constant
|
|
Interval& operator+=(real c)
|
|
{
|
|
alpha += c;
|
|
return *this;
|
|
}
|
|
|
|
// Addition by another Interval
|
|
Interval& operator+=(const Interval& rhs)
|
|
{
|
|
alpha += rhs.alpha;
|
|
beta += rhs.beta;
|
|
eps += rhs.eps;
|
|
return *this;
|
|
}
|
|
|
|
// Negation by a constant
|
|
Interval& operator-=(real c)
|
|
{
|
|
alpha -= c;
|
|
return *this;
|
|
}
|
|
|
|
// Negation by another Interval
|
|
Interval& operator-=(const Interval& rhs)
|
|
{
|
|
// Note: if rhs == *this, the result could be considered the zero function, but this code may double eps
|
|
alpha -= rhs.alpha;
|
|
beta -= rhs.beta;
|
|
eps += rhs.eps;
|
|
return *this;
|
|
}
|
|
|
|
// Multiplication by a constant
|
|
Interval& operator*=(real c)
|
|
{
|
|
using std::abs;
|
|
alpha *= c;
|
|
beta *= c;
|
|
eps *= abs(c);
|
|
return *this;
|
|
}
|
|
|
|
// Multiplication by another Interval
|
|
Interval& operator*=(const Interval& rhs)
|
|
{
|
|
using std::abs;
|
|
// This interval's function is: (*this)(y) = alpha + beta . y + [-eps, eps]
|
|
// The rhs's interval's function is: rhs(y) = rhs.alpha + rhs.beta . y + [-rhs.eps, rhs.eps]
|
|
// The product of the two functions takes the form
|
|
// (*this * rhs)(y) = alpha*rhs.alpha + (alpha*rhs.beta + rhs.alpha*beta).y + (beta . y)(rhs.beta . y) + [-eps, eps]*(rhs.alpha + rhs.beta . y) + [-rhs.eps, rhs.eps]*(alpha + beta . y) + [-eps,eps] * [-rhs.eps,rhs.eps]
|
|
real ell_this = maxDeviation();
|
|
real ell_rhs = rhs.maxDeviation();
|
|
// "Reverse ordering" so old values are not prematurely overwritten; this works even if &rhs == this.
|
|
eps = ell_this * ell_rhs + abs(alpha) * rhs.eps + abs(rhs.alpha) * eps;
|
|
for (int dim = 0; dim < N; ++dim)
|
|
beta(dim) = alpha * rhs.beta(dim) + rhs.alpha * beta(dim);
|
|
alpha *= rhs.alpha;
|
|
return *this;
|
|
}
|
|
|
|
// Division by a constant, presumed != 0
|
|
Interval& operator/=(real c)
|
|
{
|
|
using std::abs;
|
|
alpha /= c;
|
|
beta /= c;
|
|
eps /= abs(c);
|
|
return *this;
|
|
}
|
|
|
|
// Division by another interval
|
|
Interval& operator/=(const Interval& rhs)
|
|
{
|
|
using std::abs;
|
|
real rhsxinv = 1.0 / rhs.alpha;
|
|
real tau = rhs.maxDeviation()*abs(rhsxinv);
|
|
if (tau >= 1.0)
|
|
throw std::domain_error("Unable to perform operator/= with supplied argument");
|
|
real r = util::sqr(tau)/util::cube(1.0 - tau);
|
|
real bnum = maxDeviation();
|
|
eps = (abs(alpha*rhsxinv)*rhs.eps + (abs(alpha) + bnum)*r + bnum*tau + eps) * abs(rhsxinv);
|
|
for (int dim = 0; dim < N; ++dim)
|
|
beta(dim) = beta(dim)*rhsxinv - alpha*rhs.beta(dim)*util::sqr(rhsxinv);
|
|
alpha *= rhsxinv;
|
|
return *this;
|
|
}
|
|
|
|
// If the interval's values are uniformly positive or negative, returns +1 or -1, respectively. If
|
|
// the values cannot be guaranteed to have uniform sign, returns 0. Internally, this function uses
|
|
// a small tolerance based on machine precision, so that, e.g., x -> x has uniform sign on the
|
|
// interval [0,1]. The function x -> 0 is declared to have sign 0.
|
|
int sign() const
|
|
{
|
|
real x = (real(1.0) - real(10.0)*real(std::numeric_limits<real>::epsilon()))*maxDeviation();
|
|
if (alpha > x)
|
|
return 1;
|
|
else if (alpha < -x)
|
|
return -1;
|
|
else
|
|
return 0;
|
|
}
|
|
|
|
// Returns whether the interval's value are guaranteed to have uniform sign
|
|
bool uniformSign() const
|
|
{
|
|
return sign() != 0;
|
|
}
|
|
|
|
// Half-width of interval calculations, across all dimensions
|
|
static uvector<real,N>& delta()
|
|
{
|
|
// thread_local to allow different threads to perform separate chains of Interval arithmetic.
|
|
// Originally, delta was a static member variable of Interval<N> but GCC has compiler bugs
|
|
// concerning thread_local variables, possibly related to
|
|
// https://gcc.gnu.org/bugzilla/show_bug.cgi?id=81880. Changed to a static member function
|
|
// instead.
|
|
static thread_local uvector<real,N> d{};
|
|
return d;
|
|
}
|
|
|
|
// Half-width of interval calculations, for a specific dimension
|
|
static real& delta(int dim)
|
|
{
|
|
return delta()(dim);
|
|
}
|
|
};
|
|
|
|
// Interval<N> binary operators
|
|
|
|
template<int N>
|
|
Interval<N> operator+(Interval<N> alpha, const Interval<N>& y)
|
|
{
|
|
return alpha += y;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator+(Interval<N> alpha, real y)
|
|
{
|
|
return alpha += y;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator+(real alpha, Interval<N> y)
|
|
{
|
|
return y += alpha;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator-(Interval<N> alpha, const Interval<N>& y)
|
|
{
|
|
return alpha -= y;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator-(Interval<N> alpha, real y)
|
|
{
|
|
return alpha -= y;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator-(real alpha, Interval<N> y)
|
|
{
|
|
y *= -1.0;
|
|
return y += alpha;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator*(Interval<N> alpha, const Interval<N>& y)
|
|
{
|
|
return alpha *= y;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator*(Interval<N> alpha, real y)
|
|
{
|
|
return alpha *= y;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator*(real alpha, Interval<N> y)
|
|
{
|
|
return y *= alpha;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator/(Interval<N> alpha, const Interval<N>& y)
|
|
{
|
|
return alpha /= y;
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> operator/(Interval<N> alpha, real y)
|
|
{
|
|
return alpha /= y;
|
|
}
|
|
|
|
// Streaming operator for inspection
|
|
template<int N>
|
|
std::ostream& operator<< (std::ostream& o, const Interval<N>& alpha)
|
|
{
|
|
return o << "Interval<" << N << ">(" << alpha.alpha << ", " << alpha.beta << ", " << alpha.eps << ')';
|
|
}
|
|
|
|
// Various unary function evaluations of Interval objects
|
|
|
|
template<int N>
|
|
Interval<N> sin(const Interval<N>& i)
|
|
{
|
|
using std::sin;
|
|
using std::cos;
|
|
using std::abs;
|
|
real c = cos(i.alpha);
|
|
return Interval<N>(sin(i.alpha), c * i.beta, abs(c)*i.eps + 0.5*util::sqr(i.maxDeviation()));
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> cos(const Interval<N>& i)
|
|
{
|
|
using std::sin;
|
|
using std::cos;
|
|
using std::abs;
|
|
real s = -sin(i.alpha);
|
|
return Interval<N>(cos(i.alpha), s * i.beta, abs(s)*i.eps + 0.5*util::sqr(i.maxDeviation()));
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> sqrt(const Interval<N>& i)
|
|
{
|
|
using std::sqrt;
|
|
real b = i.maxDeviation();
|
|
if (b > i.alpha)
|
|
throw std::domain_error("Unable to compute sqrt() with supplied argument");
|
|
real sqrtx = sqrt(i.alpha);
|
|
real sqrtxinv = 1.0 / sqrtx;
|
|
real C = 0.25/((i.alpha - b)*sqrt(i.alpha - b));
|
|
return Interval<N>(sqrtx, (0.5*sqrtxinv)*i.beta, 0.5 * C * util::sqr(b));
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> exp(const Interval<N>& i)
|
|
{
|
|
using std::exp;
|
|
real ex = exp(i.alpha);
|
|
real b = i.maxDeviation();
|
|
// Bounding (exp(alpha))'' using montonocity and thus argument [i.alpha + b]
|
|
return Interval<N>(ex, ex * i.beta, ex*i.eps + 0.5 * exp(i.alpha + b) * util::sqr(b));
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> cosh(const Interval<N>& z)
|
|
{
|
|
using std::sinh;
|
|
using std::cosh;
|
|
using std::abs;
|
|
real coshz = cosh(z.alpha);
|
|
real sinhz = sinh(z.alpha);
|
|
real b = z.maxDeviation();
|
|
// Second derivative of cosh(alpha) is cosh(alpha), and is even and monotone increasing for alpha >= 0. Therefore
|
|
// for z.alpha - z.maxDeviation <= alpha <= z.alpha + z.maxDeviation, the maximum value for cosh(alpha) is cosh(abs(z.alpha) + z.maxDeviation)
|
|
return Interval<N>(coshz, sinhz * z.beta, abs(sinhz)*z.eps + 0.5 * cosh(abs(z.alpha) + b) * util::sqr(b));
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> tanh(const Interval<N>& z)
|
|
{
|
|
using std::tanh;
|
|
using std::cosh;
|
|
using std::abs;
|
|
real tanhz = tanh(z.alpha);
|
|
real sechzsqr = util::sqr(1.0/cosh(z.alpha));
|
|
// The second derivative of tanh(alpha) is -2*tanh(alpha)*sech(alpha)^2, and has maximum absolute value 0.76980035892...
|
|
return Interval<N>(tanhz, sechzsqr * z.beta, abs(sechzsqr)*z.eps + 0.5 * 0.77 * util::sqr(z.maxDeviation()));
|
|
}
|
|
|
|
inline real sech(real alpha)
|
|
{
|
|
using std::cosh;
|
|
return real(1)/cosh(alpha);
|
|
}
|
|
|
|
template<int N>
|
|
Interval<N> sech(const Interval<N>& z)
|
|
{
|
|
using std::tanh;
|
|
using std::abs;
|
|
real sechz = sech(z.alpha);
|
|
real fprime = -sechz*tanh(z.alpha);
|
|
uvector<real,N> beta = fprime * z.beta;
|
|
// The second derivative of sech(alpha) is {sech(alpha)*tanh(alpha)^2 - sech(alpha)^3} and has maximum absolute value 1.0
|
|
return Interval<N>(sechz, beta, abs(fprime)*z.eps + 0.5 * 1.0 * util::sqr(z.maxDeviation()));
|
|
}
|
|
} // namespace algoim
|
|
|
|
#endif
|
|
|