Math

Math provides various useful functions and constants. A number of useful function objects that operate with certain procedures in this part of the library is available from the Functors section.

Interface

Scalar Functions

  template <typename T, typename U>
  T round_cast(const U &u)  
A rounding type-cast, useful when converting floating-point types to integer types. When casting between floating-point types, rounding is not applied.
  template <typename T>
  T log2(const T n)  
Returns the largest result such that (2^result <= n). This version is intended for use with integer types, use std::tr1::log2() for floating point types.
  template <typename T>
  bool isPower2(const T x)  
Returns true if x is a power of two. Only for integral types.
  template <typename T>
  T clamp(const T &x, const T &low, const T &high)  
Returns the value of in the range [low, high] closest to x.
  template <typename T>
  T factorial(const T &x)  
Returns x!.
  template <typename T>
  T gcd(T u, T v)  
Greatest common divisor.
    template <typename T>
  T binomial(const T n, const T k)  
Binomial coefficient, choose k out of n.
    template <typename T, typename U>
  U binopmf(const T n, const T k, const U p)  
Probability mass function for the binomial distribution.
    template <typename T, typename U>
  U binocdf(const T n, const T k, const U p)  
Cumulative distribution function for the binomial distribution.

Vector Functions

  template <typename T>
  ValueType<T>::value_type modulus(const T &x)  
Computes the modulus of scalars and Vectors.
  template <typename Iterator, typename T>
  T average(Iterator begin, Iterator end, const T init)  
Returns the average of the values in the range [begin, end].
  template <typename Iterator, typename T>
  T median(Iterator begin, Iterator end, const T init)  
Returns the median of the values in the range [begin, end].
  template <typename Iterator, typename T>
  T deviation(Iterator begin, Iterator end, T init)  
Returns init plus the sample standard deviation using the samples in the given range, computed as the square root of the sample variance.
  template <typename Iterator, typename T>
  T deviation(Iterator begin, Iterator end, T avg, T init)  
As the previous, except that the argument avg must be the mean value of the given samples, which makes it faster to compute than the version without this argument.
  template <typename Iterator, typename T>
  T variance(Iterator begin, Iterator end, T init)  
Returns init plus the unbiased sample variance, i.e. normalized by N-1 rather than N, where N is the number of samples in the range.
  template <typename Iterator, typename T>
  T variance(Iterator begin, Iterator end, T avg, T init)  
As the previous, except that the argument avg must be the mean value of the given samples, which makes it faster to compute than the version without this argument.
  template <typename Iterator, typename T>
  T generalizedGaussianShape(Iterator begin, Iterator end,
			     T init)  
Estimates the ``shape''-parameter (alpha) of the Generalized Gaussian Distribution(GGD), based on the given samples, plus init.
  template <typename Iterator, typename T>
  T generalizedGaussianShape(Iterator begin, Iterator end,
			     T avg, T init)  
As the previous, except that the argument avg must be the mean value of the given samples, which makes it faster to compute than the version without this argument.

Constants

Constants of different types are defined through a template class:

  template <typename T> struct Constants;  

  value_type  
Alias for template parameter type T.
  value_type pi()  
Pi, the ratio of a circle's circumference to its diameter, i.e. 3.14.....
  value_type e()  
Base of the natural logarithm.
  std::complex<value_type> i()  
Complex component i is defined as sqrt(-1).

Compatibility Macros

Macro defined for backwards compatibility.

  PI  
Equal to Constants<double>::pi().

Procedures

Optimization

template <class Function, class Derivative>
typename Function::value_type optimize(const Function &f, const Derivative &d,
		typename Function::value_type low,
		typename Function::value_type high,
		const std::size_t N = 1024);
template <class Function>
typename Function::value_type optimize(const Function &f,
		typename Function::value_type low,
		typename Function::value_type high,
		const std::size_t N = 1024);

Newton-Raphson

The Newton-Raphson procedure is a so-called "root-finding" method, which, given a function f(x), tries to find an x such that f(x) = 0. The procedure is iterative, and computation stops upon convergence, i.e. when it appears that the solution is no longer improving. The implementation can handle multi-dimensional function by using Vectors.

As described here, the procedure is not entirely safe, and although the implementation addresses these weaknesses, it may fail. As a rule of the thumb, try to provide a starting point close to the solution and avoid having a local minimum or maximum in the function between the starting point and the solution.

There are two definitions of the classic method:

template <class Function, class Derivative>
bool doNewtonRaphson(const Function &f, const Derivative &d, Function::value_type &x,
		     unsigned N = 100000u);
template <class Function>
bool doNewtonRaphson(const Function &f, Function::value_type &x, unsigned N = 100000u);
And there are two that involve a fixed range into which x is clamped:
template <class Function, class Derivative>
bool doNewtonRaphson(const Function &f, const Derivative &d, Function::value_type &x,
		     const Function::value_type &low, const Function::value_type &high,
		     unsigned N = 100000u);
template <class Function>
bool doNewtonRaphson(const Function &f, Function::value_type &x,
		     const Function::value_type &low, const Function::value_type &high,
		     unsigned N = 100000u);
The first form requires a user-defined function Derivative d, which must be the derivative of Function f. The second form does not require a user-supplied derivative but uses a numerical technique, see Derivative. The first form is recommended as it is expected to be both faster and more accurate. On input, the parameter x is the start-point of the procedure; on output it contains the result if the procedure has succesfully terminated. If the procedure has succesfully converged onto a solution, the return value is true, otherwise it is false and the contents of x are unspecified.

The "clamped" versions are a bit safer, and potentially faster, as they make sure x does not leave the given range; and employ a bi-section scheme to speed up convergence if possible. The method works best when the function is smaller than zero on one side of the range and grater than zero on the other side.

The template parameter classes Function and Derivative must model a well-behaved mathematical function, meaning that for any given parameter x, the corresponding f(x) is always the same. Both Function and Derivative must offer a typedef for value_type and an operator of the following form:

  value_type operator()(const value_type &x) const  

An example:

template <typename T>
struct Square
{
	typedef T value_type;
	T operator()(const T &x) const { return x * x; }
};

template <typename T>
struct SquareDerivative
{
	typedef T value_type;
	T operator()(const T &x) const { return T(2)*x; }
};

// Automatic Derivative
float r1 = -10.0f;
bool ok = doNewtonRaphson(Square<float>(), r1);
assert(ok);

// With User-Supplied Derivative
double r2 = -10.0;
ok = doNewtonRaphson(Square<double>(), SquareDerivative<double>(), r2);
assert(ok);

// Multi-Dimensional through Vectors
StaticVector<double, 2u> r3 = -10.0;
ok = doNewtonRaphson(Square<StaticVector<double, 2u> >(), r3);
assert(ok);

Runge-Kutta 4

The Runge-Kutta 4 algorithm is an iterative method for the approximation of solutions of ordinary differential equations of the form:

y' = f(t, y)
y(t0) = y0
In this case y' is the derivative of y, and its value in point t is given by f(t, y). The so-called "initial value" is y0, and it is the value of y(t0), where t0 is a specified point. The procedure will produce a vector of N + 1 values, corresponding to solution of the above initial-value problem in N + 1 points in the range [t0, tN].

The method has the following definition:

template <class Function>
void doRungeKutta(const Function &f, const Function::value_type t0,
		  const Function::value_type tN, const unsigned N,
		  const Function::value_type y0,
		  std::vector<Function::value_type> &y);

The template parameter class Function must model a well-behaved mathematical function, meaning that for any given parameters t, y, the corresponding f(t, y) is always the same. Function must offer a typedef for value_type and an operator of the following form:

value_type operator()(const value_type &t, const value_type &y) const

The parameters are defined as follows:

  const Function &f  
A Functor giving the derivative y' = f(t, y).
  const Function::value_type t0  
The position of the initial condition y(t0) = y0.
  const Function::value_type tN  
The final position.
  const unsigned N  
The number of iterations.
  const Function::value_type y0  
The initial condition, i.e. y(t0)=y0.
  std::vector<Function::value_type> &y  
The vector where the values of y will be stored.
Upon return, y will contain N + 1 values.

An example:

template <typename T>
struct RKSquare
{
	typedef T value_type;
	typedef T result_type;
	T operator()(const T t, const T y) const { return t*t-y; }
};

RKSquare<double> f;
std::vector<double> result;

// Compute 10 values in the range [0.0, 1.0], where y(0.0) = 3.0
doRungeKutta(f, 0.0, 1.0, 10u, 3.0, result);

Least Squares

Implements the well-known least-squares method. Here, the Matrix A is the design matrix, vector y contains observations, and the estimated parameters will be placed in vector x.
template <template <typename Tm, std::size_t Dm, typename Aux> class Array_t,
	  typename Ta, typename A, class XVector_t, class YVector_t>
bool leastSquaresFit(const Array_t<Ta, 2, A> &A, const YVector_t &y, XVector_t &x)

Hungarian Algorithm

The hungarian algorithm finds a matching in a bi-partite graph that minimizes overall cost or maximizes overall benefit. That is, given two sets of items A and B, the best possible set of pairs having one item from set A and one item from set B is found.

The overall cost or benefit is calculated as the sum of the costs or benefits from the individual pairs. Should your overall cost function actually be a product of the individual costs or benefits, then simply provide the log of your cost function.

Note: the costs matrix is destroyed during the calculation!

template <template <typename Tm, std::size_t D, typename Aux> class Array_t,
	typename T, typename A>
void find_matching(Array_t<T, 2, A> &costs,
		   std::vector<std::pair<std::size_t, std::size_t> > &matches,
		   const bool minimize_costs = true,
		   const bool heuristic = true);

  costs  
A 2-dimensional matrix containing at each element costs[a][b] the cost or benefit of a match between elements a and b from their respective sets.
  matches  
On exit, contains a set of pairs of indices <a,b> that each represent a match.
  minimize_costs  
If true (default) the costs will be minimized. If false, the costs (benefits) will be maximized.
  heuristic  
Indicate whether a heuristic should be used.

Matrix Operations

  void transpose(Array_t<T, 2u, A> &m)  
Transpose the Matrix.
  bool invert(Array_t<T, 2u, A> &m)  
Invert the matrix. Returns false if inversion fails.