Miscellaneous mathematical functions that aren't in std::math More...

Classes
class	Approx_equal
	functor to compare whether values are approximately equal More...

struct	Ball
	minimal representation of an `n_dim`-dimensional ball More...

struct	Objective_finite_diff
	Estimates the derivatives of an optimization objective function by finite difference. More...

struct	Root_options
	provides a convenient way to pass options to root-finding algorithms More...

struct	Tolerance
	Specifies relative and absolute tolerances for comparing a value to truth data. More...

Functions
double	limit_abs (double value, double limit)

double	smooth_limit_abs (double value, double limit)

double	random_normal (double mean, double std_dev)

template<typename T >
void	set_random_normal (T &mat, double mean, double std_dev)
	Sets the entries of the provided matrix-like object to random values with a normal distribution. The mean and standard deviation of the distribution are controlled by `mean` and `std_dev`, respectively. `T` must support the methods `.rows()`, `.cols()`, and `operator()(int row, int col)` (for entry access).

template<typename number_t >
constexpr number_t	pow (number_t base, int exponent)
	Raises an arbitrary arithmetic type to an integer (not necessarily positive) power.

constexpr Int	log (Int base, Int arg)
	Integer logarithm.

template<typename T >
constexpr T	mod (T i, T j)
	Modulo operator.

template<typename T >
constexpr T	extreme (bool minmax, T arg)

template<typename T , typename... arg_ts>
constexpr T	extreme (bool minmax, T arg, arg_ts... args)
	if `minmax` is true, returns the maximum of remaining arguments, else the minimum

template<typename T , typename... arg_ts>
constexpr T	max (T arg, arg_ts... args)
	returns the maximum of all arguments

template<typename T , typename... arg_ts>
constexpr T	min (T arg, arg_ts... args)
	returns the minimum of all arguments

constexpr int	sign (bool condition)
	returns 1 if `condition` is true, otherwise -1

constexpr Int	stride (int n_dim, Int row_size, int i_dim)

constexpr Int	row_coordinate (int n_dim, Int row_size, int i_dim, Int index)
	Finds the `i_dim`th array index of a point with flat index `index` in an `n_dim` dimensional array of `row_size` on each side.

double	angle_diff (double angle0, double angle1)
	returns `angle0 - angle1`, where the difference is in \( [0, 2\pi) \)

double	broyden (std::function< double(double)> error, double init_guess, Root_options opts, double init_diff=1e-3)
	Finds a root of a scalar function with Broyden's method.

double	bisection (std::function< double(double)> error, std::array< double, 2 > bounds, Root_options opts)
	Finds a root of a scalar function with the bisection method.

Mat	newton (std::function< Mat< dyn, dyn >(Mat<>)> error_jacobian, Mat<> guess, Root_options options)
	Finds a root of a vector function with Newton's method.

Eigen::VectorXd	hypercube_matvec (const Eigen::MatrixXd &, const Eigen::VectorXd &)
	Multiply every dimension of a (flattened) N-dimensional array by a matrix.

Eigen::VectorXd	dimension_matvec (const Eigen::MatrixXd &, const Eigen::VectorXd &, int i_dim)
	Multiply a single dimension of a hypercubic ND array by a matrix.

Eigen::VectorXd	pow_outer (const Eigen::VectorXd &, int n_dim)
	Raises a vector to a power via ND outer products.

template<int n_dim>
Eigen::Matrix< double, n_dim, n_dim >	orthonormal (Eigen::Matrix< double, n_dim, n_dim > basis, int i_dim)
	Orthonormalize a vector basis (with dimension \(\le 3\)).

Eigen::MatrixXd	orthonormal (Eigen::MatrixXd basis, int i_dim)

int	stretched_ind (int n_dim, int ind, std::array< bool, 2 > stretch)
	for indexing faces/vertices in Refined_face s with possible stretching

template<int n_dim>
double	interp (Mat< pow(2, n_dim)> values, Mat< n_dim > coords)
	\(n\)-linear interpolation of `values`.

double	interp (Mat<> values, Mat<> coords)
	This is an overloaded member function, provided for convenience. It differs from the above function only in what argument(s) it accepts.

template<typename vec_t >
vec_t	proj_to_segment (std::array< vec_t, 2 > endpoints, vec_t target)
	Finds the nearest point to `target` on the line segment defined by `endpoints`.

template<int n_dim, int n_point>
Ball< n_dim >	bounding_ball (Mat< n_dim, n_point > points)
	computes a bounding ball of the convex hull of a set of points

template<int n_dim>
bool	intersects (Ball< n_dim > b, Mat< n_dim > endpoint0, Mat< n_dim > endpoint1)
	returns true if the ball `b` intersects the line through `endpoint0` and `endpoint1`

double	chebyshev_step (int n_steps, int i_step, double safety=.9)
	the scaling factor for Chebyshev polynomial multistge time stepping

std::vector< double >	correct_values (std::vector< double > estimates, std::vector< double > exacts, double tol=huge)

std::mt19937	global_random_generator (global_random_device())

Variables
std::random_device	global_random_device

std::mt19937	global_random_generator

Detailed Description

Miscellaneous mathematical functions that aren't in std::math

Function Documentation

◆ bisection()

double hexed::math::bisection	(	std::function< double(double)>	error,
		std::array< double, 2 >	bounds,
		Root_options	opts )

Finds a root of a scalar function with the bisection method.

This is slower than broyden() but very robust.

Parameters

error	function to find the root of
bounds	Lower and upper bounds for a root.
opts	Use this to set the termination condition

◆ bounding_ball()

template<int n_dim, int n_point>

Ball< n_dim > hexed::math::bounding_ball ( Mat< n_dim, n_point > points )

computes a bounding ball of the convex hull of a set of points

Returns the smallest bounding ball centered at the center of mass of the points. This bounding ball is non-optimal but quick to compute. Note that a simplex is the convex hull of its vertices.

◆ broyden()

double hexed::math::broyden	(	std::function< double(double)>	error,
		double	init_guess,
		Root_options	opts,
		double	init_diff = 1e-3 )

Finds a root of a scalar function with Broyden's method.

Parameters

error	error function to find the root of
init_guess	Initial guess for the root.
opts	Use this to set the termination condition
init_diff	How far away the second point used to initialize the derivative estimate should be.

◆ chebyshev_step()

double hexed::math::chebyshev_step	(	int	n_steps,
		int	i_step,
		double	safety = .9 )

the scaling factor for Chebyshev polynomial multistge time stepping

Returns the ratio of the i_stepth step of n_steps to the nominal time step. A value of safety = 1 gives you the maximum stable time step and safety = 0 gives you a time step of 0 (although it isn't a linear function).

◆ correct_values()

std::vector< double > hexed::math::correct_values	(	std::vector< double >	estimates,
		std::vector< double >	exacts,
		double	tol = huge )

Suppose that you compute some values, estimates, with a method that is very robust but has low accuracy. Suppose you also recompute these values as exacts with a method that is very accurate but not robust—It may give you extraneous values and/or fail to produce some of the correct values. The purpose of this function is to filter out the implausible values from exacts and give you only the ones that appear to be correct. The return value will have the same number of values as estimates, but some of them will be replaced with values from exacts in a way that obtains the best possible agreement. Replacements will be considered only if the difference between the "exact" value and the "estimate" is less than tol. The order of exacts does not matter. Values will be returned in the same order as estimates.

◆ dimension_matvec()

Eigen::VectorXd hexed::math::dimension_matvec	(	const Eigen::MatrixXd &	mat,
		const Eigen::VectorXd &	vec,
		int	i_dim )

Multiply a single dimension of a hypercubic ND array by a matrix.

C.f. hypercube_matvec. If matrix is square, the shape of the output array will match the input. If matrix is a row vector, the output will still be hypercubic, but with one less dimension than the input. Otherwise, the resulting array will no longer be hypercubic.

◆ extreme()

template<typename T , typename... arg_ts>

T hexed::math::extreme	(	bool	minmax,
		T	arg,
		arg_ts...	args )

constexpr

if minmax is true, returns the maximum of remaining arguments, else the minimum

Warning: always returns the type of the first argument, regardless of type conversion rules

◆ hypercube_matvec()

Eigen::VectorXd hexed::math::hypercube_matvec	(	const Eigen::MatrixXd &	mat,
		const Eigen::VectorXd &	vec )

Multiply every dimension of a (flattened) N-dimensional array by a matrix.

Size of array along each dimension must be equal (i.e. the array is hypercube-shaped, or in my terminology, "hypercubic"). Matrix does not have to be square. Dimensionality of the array is inferred automatically by comparing the number of matrix columns to the array size. Array size must be a power of the number of matrix columns.

◆ interp()

template<int n_dim>

double hexed::math::interp	(	Mat< pow(2, n_dim)>	values,
		Mat< n_dim >	coords )

\(n\)-linear interpolation of values.

ND generalization of bilinear interpolation.

Parameters

values	Values to interpolate. Assumed to be at corners of the unit hypercube.
coords	Coordinates to interpolate to.

◆ log()

Int hexed::math::log	(	Int	base,
		Int	arg )

constexpr

Integer logarithm.

If base < 2, returns -1 to indicate failure. Otherwise, if arg < 1, returns 0. In the usual case where neither of the above are true, returns \(\lceil\log_{\mathtt{base}}(\mathtt{arg})\rceil\).

◆ max()

template<typename T , typename... arg_ts>

T hexed::math::max	(	T	arg,
		arg_ts...	args )

constexpr

returns the maximum of all arguments

Warning: always returns the type of the first argument, regardless of type conversion rules

◆ min()

template<typename T , typename... arg_ts>

T hexed::math::min	(	T	arg,
		arg_ts...	args )

constexpr

returns the minimum of all arguments

Warning: always returns the type of the first argument, regardless of type conversion rules

◆ mod()

template<typename T >

T hexed::math::mod	(	T	i,
		T	j )

constexpr

Modulo operator.

Similar to the remainder operator (ij) but returns a nonnegative result even when i is negative.

◆ newton()

Mat hexed::math::newton	(	std::function< Mat< dyn, dyn >(Mat<>)>	error_jacobian,
		Mat<>	guess,
		Root_options	options )

Finds a root of a vector function with Newton's method.

Parameters

error_jacobian	Function to find the root of. If the input ( \( x \)) has \( n \) entries, this should be an \( n \times n + 1 \) matrix. The first column should be the error vector, and the remaing columns should be the Jacobian matrix, such that entry i, j of the Jacobian (i, j + 1 of the return matrix) should be the derivative of the ith component of the error w.r.t. the jth component of the input.
guess	initial guess for the root
options	use this to set the termination condition

◆ orthonormal()

template<int n_dim>

Eigen::Matrix< double, n_dim, n_dim > hexed::math::orthonormal	(	Eigen::Matrix< double, n_dim, n_dim >	basis,
		int	i_dim )

Orthonormalize a vector basis (with dimension \(\le 3\)).

Assumes basis is invertible. Returns a matrix with the following properties:

Unitary.
Span of columns excluding the i_dimth is the same as for basis.
Minimizes RMS difference between columns excluding i_dimth of return matrix and basis (sensitive to order).
Inner product of i_dimth columns of return matrix and basis is positive.

◆ pow()

template<typename number_t >

number_t hexed::math::pow	(	number_t	base,
		int	exponent )

constexpr

Raises an arbitrary arithmetic type to an integer (not necessarily positive) power.

Can return constexpr, which std::pow is not allowed to do according to the standard (although the GCC implementation can anyway).

◆ pow_outer()

Eigen::VectorXd hexed::math::pow_outer	(	const Eigen::VectorXd &	vec,
		int	n_dim )

Raises a vector to a power via ND outer products.

That is, takes an outer product with the vector {1} n_dim times along different dimensions to produce an n_dim-dimensional hypercubic array.

◆ proj_to_segment()

template<typename vec_t >

vec_t hexed::math::proj_to_segment	(	std::array< vec_t, 2 >	endpoints,
		vec_t	target )

Finds the nearest point to target on the line segment defined by endpoints.

Works for 2D or 3D.

Classes

Functions

Variables

Detailed Description

Function Documentation

◆ bisection()

◆ bounding_ball()

◆ broyden()

◆ chebyshev_step()

◆ correct_values()

◆ dimension_matvec()

◆ extreme()

◆ hypercube_matvec()

◆ interp()

◆ log()

◆ max()

◆ min()

◆ mod()

◆ newton()

◆ orthonormal()

◆ pow()

◆ pow_outer()

◆ proj_to_segment()