1-d and n-d Functions
Mathwrist C++ API Series

(October 5, 2023)

Abstract

This document details the C++ programming interface of 1-d and n-d mathematical functions, 1-d interpolation, 1-d integration, 1-d root finding and 1-d minimization features in Mathwrist’s Numerical Programming Library (NPL). Please refer to of our 1-d and n-d function white paper for a technical overview of these features.

1 1-d Functions

1.1 Function1D Base Class

Function1D is the pure virtual base class for all 1-dimensional functions. It is defined in header file “function_1d.h”.

1.1.1 Clone

virtual Function1D* clone() const = 0;

Description:

For concrete derived classes to create a cloned object.

1.1.2 Evaluation

virtual double eval(

double x,

double* _derv1 = 0,

double* _derv2 = 0) const = 0;

Description:

Return the function value $f(x)$ at a given point $x$ . Parameter _derv1 and _derv2 are optional. If non-zero pointers are provided, it is expected to compute the first and second derivatives of $f(x)$ respectively and save the result to the supplied pointers.

1.1.3 Integration

virtual double integral(double lo, double hi) const;

Description:

Return the finite integral $\int_{lo}^{hi}f(x)dx$ of the function. The Function1D base class by default throws an un-supported exception in this integral() function. It is upto the derived classes to override this default behavior if necessary.

1.1.4 Finite Difference Approximation

1 double derv1_fwd_diff(double x) const;

2 double derv1_ctr_diff(double x) const;

3 double derv2_fwd_diff_derv1(double x) const;

4 double derv2_diff_func(double x, double* _base_value = 0) const;

Description:

For convenience, the Function1D base class provides these derivative approximation functions using finite difference method. A smooth function derived from the Function1D base class may use these approximation methods as appropriate. The default bump size is $\Delta x=\sqrt{\epsilon}$ , where $\epsilon$ is the machine epsilon. A derived class may change this bump size by overriding a protected member function Function1D::pertubed_size().

1.

Return the approximation of the first derivative $f^{\prime}(x)$ by forward finite difference. $f^{\prime}(x)\approx\frac{f(x+\Delta x)-f(x)}{\Delta x}$ .
2.

Return the approximation of the first derivative $f^{\prime}(x)$ by center finite difference. $f^{\prime}(x)\approx\frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}$ .
3.

Return the approximation of the second derivative $f^{\prime\prime}(x)$ by forward finite difference on the first derivatives. $f^{\prime\prime}(x)\approx\frac{f^{\prime}(x+\Delta x)-f^{\prime}(x)}{\Delta x}$ . The implementation of function $f(x)$ is expected to compute first derivatives $f^{\prime}(x+\Delta x)$ and $f^{\prime}(x)$ .
4.

Return the approximation of the second derivative $f^{\prime\prime}(x)$ by center finite difference. $f^{\prime\prime}(x)\approx\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta x^{2}}$ . Optional parameter _base_value holds the function value $f(x)$ , if a non-zero pointer is supplied.

1.2 1-d Piecewise Function

Class FunctPiecewise is an intermediate pure virtual class to represent 1-d piecewise functions, defined in C++ header file “funtion_pw.h”. It is derived from Function1D base class. Concrete piecewise functions are further derived from this FunctPiecewise intermediate class.

class FunctPiecewise : public Function1D

{

};

1.2.1 Piece Lookup

static size_t offset(

const double* _x_begin, size_t n_x, double x);

Description:

Given a buffer of knot points $(x_{0},x_{1},\cdots,x_{n-1})$ in ascending order, return the offset set $i$ such that $x_{i}\leq x<x_{i+1}$ or return $n-2$ if $x$ is numerically equal to $x_{n-1}$ .

Parameter _x_begin is a pointer to the buffer address. Parameter n_x specifies the number of elements $n$ in the buffer. Parameter x is the lookup value.

REQUIRED: $x\in[x_{0},x_{n-1}]$ . An out of range exception will be thrown otherwise.

1.2.2 Knot Point Functions

1 size_t n_breaks() const;

2 const double* x() const;

Description:

1.

Return the number of knot points.
2.

Return the address of the knot point buffer.

1.2.3 Integration

1 double integral(double lo, double hi) const;

2 virtual void refresh() = 0;

Description:

1.

Implement the virtual integral() function declared in Function1D base class by summing up piecewise integrals.
2.

Refresh the internal cache of piecewise integrals, i.e after function value changes.

1.3 Piecewise Constant Function

Class FunctPiecewiseConst is a concrete type to represent 1-d piecewise constant functions. It is derived from the intermediate FunctPiecewise class.

class FunctPiecewiseConst : public FunctPiecewise

{

};

1.3.1 Constructor

1 FunctPiecewiseConst();

2 FunctPiecewiseConst(size_t n, const double* _x, const double* _y);

Description:

1.

Default contructor, create an empty object.
2.

Initialize the piecewise constant function with $n$ number of knot points $(x_{0},\cdots,x_{n-1})$ and $n-1$ function values at $y_{i}=f(x_{i}),i=0,\cdots,n-2$ .

Parameter n is the number of knot points. Parameter _x points to the buffer of knot points sorted in ascending order. Parameter _y points to the buffer of $n-1$ function values.

REQUIRED: The piecewise function internally just saves the pointer address _x and _y. Hence, the memory of these knot points and function values must be valid during the entire life of the piecewise function. The memory of _x and _y are not freed by the destructor.

1.3.2 Piecewise Values

const double* y() const;

Description:

This function returns the pointer of piecewise function values $(y_{0},\cdots,y_{n-2})$ .

1.3.3 Evaluation

double eval(

double x,

double* _derv1 = 0,

double* _derv2 = 0) const;

Description:

This function implements the public interface of the Function1D base class. It returns a value $f(x)=y_{i}$ , if $x_{i}\leq x<x_{i+1}$ . The first and second derivatives are 0. See also Function1D::eval().

1.3.4 Refresh Integral Cache

void refresh();

This function implements the public interface of the FunctPiecewise base class. Whenever the piecewise function values are changed, this refresh() function needs be called to recompute the piecewise integral cache. Please refer to FunctPiecewise::refresh() for API details.

1.4 Piecewise Linear Function

Class FunctPiecewiseLinear is a concrete type to represent 1-d continuous piecewise linear functions. It is derived from the intermediate FunctPiecewise class.

class FunctPiecewiseLinear : public FunctPiecewise

{

};

1.4.1 Constructor

1 FunctPiecewiseLinear();

2 FunctPiecewiseLinear(size_t n, const double* _x, const double* _y);

Description:

1.

Default contructor, create an empty object.
2.

Initialize the piecewise linear function with $n$ number of knot points $(x_{0},\cdots,x_{n-1})$ and $n$ function values at $y_{i}=f(x_{i}),i=0,\cdots,n-1$ .

Parameter n is the number of knot points. Parameter _x points to the buffer of knot points sorted in ascending order. Parameter _y points to the buffer of $n$ function values.

REQUIRED: The piecewise function internally just saves the pointer address _x and _y. Hence, the memory of these knot points and function values must be valid during the entire life of the piecewise function. The memory of _x and _y are not freed by the destructor.

1.4.2 Piecewise Values

const double* y() const;

Description:

This function returns the pointer of piecewise function values $(y_{0},\cdots,y_{n-1})$ .

1.4.3 Evaluation

double eval(

double x,

double* _derv1 = 0,

double* _derv2 = 0) const;

Description:

This function implements the public interface of the Function1D base class. It returns a linearly interpolated value between $y_{i}$ and $y_{i+1}$ , if $x_{i}\leq x<x_{i+1}$ . The first derivative is the linear slope. The second derivative is 0. See also Function1D::eval().

1.4.4 Refresh Integral Cache

void refresh();

This function implements the public interface of the FunctPiecewise base class. Whenever the piecewise function values are changed, this refresh() function needs be called to recompute the piecewise integral cache. Please refer to FunctPiecewise::refresh() for API details.

1.5 Piecewise Linear with Jump Function

Class FunctPiecewiseLinearJump is a concrete type to represent non-continuous 1-d piecewise linear functions. It is derived from the intermediate FunctPiecewise class.

class FunctPiecewiseLinearJump : public FunctPiecewise

{

};

1.5.1 Constructor

1 FunctPiecewiseLinearJump();

2 FunctPiecewiseLinearJump(size_t n,

3 const double* _x, const double* _y);

Description:

1.

Default contructor, create an empty object.
2.

Initialize the piecewise linear non-continuous function with $n$ number of knot points $(x_{0},\cdots,x_{n-1})$ and $2n-2$ function values.

Parameter n is the number of knot points. Parameter _x points to the buffer of knot points sorted in ascending order. Parameter _y points to the buffer of $2n-2$ function values. The function values are arranged by the start and end values of the first, second piece etc.

REQUIRED: The piecewise function internally just saves the pointer address _x and _y. Hence, the memory of these knot points and function values must be valid during the entire life of the piecewise function. The memory of _x and _y are not freed by the destructor.

1.5.2 Piecewise Values

const double* y() const;

Description:

This function returns the pointer of piecewise function values $(y_{0},\cdots,y_{2n-3})$ .

1.5.3 Evaluation

double eval(

double x,

double* _derv1 = 0,

double* _derv2 = 0) const;

Description:

This function implements the public interface of the Function1D base class. It returns a linearly interpolated value between $y_{2i}$ and $y_{2i+1}$ , if $x_{i}\leq x<x_{i+1}$ . The first derivative is the linear slope. The second derivative is 0. See also Function1D::eval().

1.5.4 Refresh Integral Cache

void refresh();

This function implements the public interface of the FunctPiecewise base class. Whenever the piecewise function values are changed, this refresh() function needs be called to recompute the piecewise integral cache. Please refer to FunctPiecewise::refresh() for API details.

1.6 B-spline Curve

Class BSplineCurve is a concrete type to represent continuous (smooth) functions using piecewise B-spline basis. It is derived from the intermediate FunctPiecewise class. BSplineCurve class is defined in the header file “b_spline_curve.h”.

class BSplineCurve : public FunctPiecewise

{

};

1.6.1 Constructor, Destructor and Copy

1 explicit BSplineCurve(size_t degree = 0);

2 BSplineCurve(const BSplineCurve& other);

3 BSplineCurve& operator= (const BSplineCurve& other);

4 ~BSplineCurve();

Description:

1.

Create a B-spline curve object of certain degree. Parameter degree is the degree of the B-spline basis. The degree is equal to B-spline’s order less 1.
2.

Copy constructor.
3.

Assignment operator to deep copy a B-spline object.
4.

Destructor, free all resources used by a B-spline object.

1.6.2 Knot Point Functions

1 void set_x_even(double x_lo, double x_hi, size_t n);

2 void set_x_customize(const double* x, size_t n);

Description:

1.

Set the knot points by evenly breaking the domain $(x_{lo},x_{hi})$ to $n$ pieces.
2.

Set the knot points with a given sequence of break points. Parameter x is the pointer address of the knot points. Parameter n is the total number of knot points.

1.6.3 Basis Functions

1 size_t n_basis() const;

2 void set_basis_coeff(const double* _beta = 0);

3 void set_basis_coeff(const Matrix& beta);

4 const ftl::Buffer<double>& basis_coeff() const;

Description:

1.

Return the total number of basis functions used by the B-spline curve.
2.

Set the coefficients of B-spline basis functions. Parameter _beta is the pointer address to the coefficients. If a zero pointer is passed, the curve sets all coefficients equal to 0.
3.

Set the coefficients of B-spline basis functions. Parameter beta is a vector of the coefficients.
4.

Return a constant reference of the B-spline basis coefficients stored in the B-spline curve object.

1.6.4 Evaluation

double eval(

double x,

double* _derv1 = 0,

double* _derv2 = 0) const;

Description:

This function implements the public interface of the Function1D base class. It returns the B-spline curve function value and computes derivatives. See also Function1D::eval().

1.7 Chebyshev Curve

Class ChebyshevCurve is a concrete type to represent continuous (smooth) functions using Chebyshev basis polynomials. It is a directly derived class from the Function1D base class. Class ChebyshevCurve is defined in header file “chebyshev_curve.h”.

class ChebyshevCurve : public Function1D

{

};

1.7.1 Constructor

ChebyshevCurve(double a, double b, size_t n_degree);

Description:

Constructor, create a Chebyshev curve object with working domain $[a,b]$ and a given degree of the approximation polynomial. Parameter a and b defines the domain $(a,b)$ . Parameter n_degree is the highest degree of Chebyshev basis polynomial to be used in the curve.

1.7.2 Query Functions

1 double x_lo() const;

2 double x_hi() const;

3 size_t degree() const;

Description:

1.

Return the lower bound of the working domain.
2.

Return the upper bound of the working domain.
3.

Return the highest degree of Chebyshev basis polynomial used in the curve.

1.7.3 Basis Functions

1 void set_basis_coeff(const double* _beta = 0);

2 void set_basis_coeff(const Matrix& beta);

3 const ftl::Buffer<double>& basis_coeff() const;

Description:

1.

Set the coefficients of Chebyshev basis polynomials. Parameter _beta is the pointer address to the coefficients. If a zero pointer is passed, the curve sets all coefficients equal to 0.
2.

Set the coefficients of Chebyshev basis polynomials. Parameter beta is a vector of the coefficients.
3.

Return a constant reference of the Chebyshev basis coefficients stored in the Chebyshev curve object.

1.7.4 Evaluation

double eval(

double x,

double* _derv1 = 0,

double* _derv2 = 0) const;

Description:

This function implements the public interface of the Function1D base class. It returns the Chebyshev curve function value and computes derivatives. See also Function1D::eval().

1.7.5 Integration

double integral(double a, double b) const;

Description:

This function implements the public interface of the Function1D base class. It computes the finite integral over $(a,b)$ . See also Function1D::eval().

1.7.6 Approximation

void approx(const Function1D& f);

Description:

This function approximates a given function $f(x)$ and sets the coefficients of Chebyshev expansion basis by matching the function values $f(x_{i}),i=0,\cdots,n$ at $n+1$ number of node points, where

•

$n$ is the highest degree of the Chebyshev polynomial;
•

node point $x_{i}=a+\frac{b-a}{2}\left[\cos\left(\frac{\pi i}{n}\right)+1\right]$ .

2 n-d Functions

2.1 FunctionND Base Class

FunctionND is the pure virtual base class for all multi-dimensional functions. It is defined in header file “function_nd.h”.

2.1.1 Clone

virtual FunctionND* clone() const = 0;

Description:

For concrete derived classes to create a cloned object.

2.1.2 Evaluation

virtual double eval(

const Matrix& x,

Matrix* _g = 0,

Matrix* _H = 0) const = 0;

Description:

Return the function value $f(x)$ at a given point $x$ . Parameter _g and _H are optional. If non-zero pointers are provided, it is expected to compute the gradient and Hessian of $f(x)$ respectively and save the result to the supplied pointers.

2.1.3 Integration

virtual double integral(

const Matrix& X_LO, const Matrix& X_HI) const;

Description:

Return the multi-dimensional finite integral $\iint_{lo}^{hi}f(x)dx$ of the function. The FunctionND base class by default throws an un-supported exception in this integral() function. It is upto the derived classes to override this default behavior if necessary. Parameters X_LO and X_HI are the vectors of lower and upper limits of the integral respectively.

2.1.4 Dimension Size

virtual size_t dimension_size() const = 0;

Description:

This function returns the actual dimension size $n$ of the n-d function.

2.1.5 Finite Difference Approximation

void grad_fwd_diff(

const Matrix& x,

Matrix& g,

const double* _f_base = 0) const;

void grad_ctr_diff(

const Matrix& x,

Matrix& g) const;

void hess_fwd_diff_grad(

const Matrix& x,

Matrix& H,

const Matrix* _g_base = 0) const;

void hess_diff_func(

const Matrix& x,

Matrix& H,

const double* _f_base = 0) const;

Description:

For convenience, the FunctionND base class provides gradient and Hessian approximation functions using finite difference method. A smooth n-d function derived from FunctionND base class may use these approximation funcitons as appropriate.

The default bump size is $\Delta x=\sqrt{\epsilon}$ , where $\epsilon$ is the machine epsilon. A derived class may change this bump size by overriding a protected member function FunctionND::pertubed_size(). Let $\Delta\mathbf{x}_{i}$ be a vector of zero elements except its $i$ -th element is equal to $\Delta x$ .

In the same order of the above listed functions,

1.

Return the approximation of the gradient by forward finite difference. $\frac{\partial}{\partial x_{i}}f(\mathbf{x})\approx\frac{f(\mathbf{x}+\Delta% \mathbf{x}_{i})-f(\mathbf{x})}{\Delta x}$ .
2.

Return the gradient approximation by center finite difference. $\frac{\partial}{\partial x_{i}}f(\mathbf{x})\approx\frac{f(\mathbf{x}+\Delta% \mathbf{x}_{i})-f(\mathbf{x}-\Delta\mathbf{x}_{i})}{2\Delta x}$ .
3.

This approximation method requires the derived n-d function to compute gradient vectors. It then uses the gradient forward difference to approximate the Hessian. $\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}f(\mathbf{x})\approx\frac{g_{% i}(\mathbf{x}+\Delta x_{j})-g_{i}(\mathbf{x})}{2\Delta x}+\frac{g_{j}(\mathbf{% x}+\Delta x_{i})-g_{j}(\mathbf{x})}{2\Delta x}$ , where $g_{i}(\mathbf{x})$ is the partial derivative $\frac{\partial}{\partial x_{i}}f(\mathbf{x})$ .
4.

This approximation doesn’t require gradient calculation. It approximate the Hessian by center difference of function values only. $\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}f(\mathbf{x})\approx\frac{f(% \mathbf{x}+\Delta\mathbf{x}_{i}+\Delta\mathbf{x}_{j})-f(\mathbf{x}+\Delta% \mathbf{x}_{i})-f(\mathbf{x}+\Delta\mathbf{x}_{j})+f(\mathbf{x})}{\Delta x^{2}}$ .

2.2 2-d Piecewise Function

Class FunctPiecewise2D is an intermediate pure virtual class to represent 2-d piecewise functions, defined in C++ header file “funtion_pw.h”. It is derived from FunctionND base class. Concrete 2-d piecewise functions are further derived from this FunctPiecewise2D intermediate class.

class FunctPiecewise2D : public FunctionND

{

};

2.2.1 Knot Point Function

1 size_t n_breaks_x() const;

2 size_t n_breaks_y() const;

3 const double* x() const;

4 const double* y() const;

Description:

1.

Return the number of knot points along $x$ axis.
2.

Return the number of knot points along $y$ axis.
3.

Return the buffer address of knot points along $x$ axis.
4.

Return the buffer address of knot points along $y$ axis.

2.2.2 Integration

double integral(const Matrix& X_LO, const Matrix& X_HI) const;

Description:

This function implements the virtual integral() function declared in FunctionND base class by summing up piecewise integrals.

2.3 2-d Piecewise Constant Function

Class FunctPiecewise2DConst is a concrete type to represent 2-d piecewise constant functions. It is derived from the intermediate FunctPiecewise2D class.

class FunctPiecewise2DConst : public FunctPiecewise2D

{

};

2.3.1 Constructors

FunctPiecewise2DConst();

FunctPiecewise2DConst(

size_t n_x, size_t n_y,

const double* _x, const double* _y,

const double* _z);

Description:

1.

Default contructor, create an empty object.
2.

Initialize the 2-d piecewise constant function with $n_{x}$ number of knot points $(x_{0},\cdots,x_{n_{x}-1})$ in $x$ axis, $n_{y}$ number of knot points $(y_{0},\cdots,y_{n_{y}-1})$ in $y$ axis and $(n_{x}-1)\times(n_{y}-1)$ number of piecewise function values at $z_{i,j}=f((x_{i},y_{j}))$ , $i=0,\cdots,n_{x}-2$ , $j=0,\cdots,n_{y}-2$ .

Parameter n_x and n_y are the number of knot points along $x$ and $y$ respectively. Parameter _x and _y point to the buffer of knot points sorted in ascending order along $x$ and $y$ respectively. Parameter _z points to the buffer of function values.

REQUIRED: The piecewise function internally just saves the pointer address _x, _y and _z. Hence, the memory of these knot points and function values must be valid during the entire life of the piecewise function. The memory of _x, _y and _z are not freed by the destructor.

2.3.2 Piecewise Values

const double* z() const;

Description:

This function returns the pointer of piecewise function values $(z_{0},\cdots)$ .

2.3.3 Evaluation

double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const;

Description:

This function implements the public interface of the FunctionND base class. It returns a value $f((x_{i},y_{j}))=z_{i,j}$ , if $x_{i}\leq x<x_{i+1}$ and $y_{j}\leq y<y_{j+1}$ . The gradient and Hessian are set to 0 values, if non-zero pointers are supplied. See also FunctionND::eval().

2.4 B-spline Surface

Class BSplineSurface is a concrete type to represent continuous (smooth) surface using the tensor product of two sets of B-spline basis along $x$ and $y$ . It is derived from the intermediate FunctPiecewise2D class. BSplineSurface class is defined in the header file “b_spline_surface.h”.

class BSplineSurface : public FunctPiecewise2D

{

};

2.4.1 Constructor, Destructor and Copy

1 BSplineSurface(size_t degree_x = 0, size_t degree_y = 0);

2 BSplineSurface(const BSplineSurface& other);

3 BSplineSurface& operator=(const BSplineSurface& other);

4 ~BSplineSurface();

Description:

1.

Create a B-spline surface of certain degrees. Parameter degree_x and degree_y are the degrees of the B-spline basis in $x$ and $y$ dimension. The degree is equal to B-spline’s order less 1.
2.

Copy constructor.
3.

Assignment operator to deep copy a B-spline surface object.
4.

Destructor, free all resources used by a B-spline surface object.

2.4.2 Knot Point Functions

1 void set_x_even(double x_lo, double x_hi, size_t n);

2 void set_y_even(double y_lo, double y_hi, size_t n);

3 void set_x_customize(const double* x, size_t n);

4 void set_y_customize(const double* y, size_t n);

Description:

1.

Set the knot points along $x$ dimension by evenly breaking the domain $(x_{lo},x_{hi})$ to $n$ pieces.
2.

Set the knot points along $y$ dimension by evenly breaking the domain $(y_{lo},y_{hi})$ to $n$ pieces.
3.

Set the knot points along $x$ dimension with a given sequence of break points. Parameter x is the pointer address of the knot points. Parameter n is the total number of knot points.
4.

Set the knot points along $y$ dimension with a given sequence of break points. Parameter y is the pointer address of the knot points. Parameter n is the total number of knot points.

2.4.3 Basis Functions

1 size_t n_basis_x() const;

2 size_t n_basis_y() const;

3 void set_basis_coeff(const double* _theta = 0);

4 void set_basis_coeff(const Matrix& theta);

5 const ftl::Buffer<double>& basis_coeff() const;

Description:

1.

Return the total number of basis functions along $x$ used by the surface.
2.

Return the total number of basis functions along $y$ used by the surface.
3.

Set the tensor product coefficients of B-spline basis functions. Parameter _theta is the pointer address to the $m\times n$ basis coefficients, where $m$ and $n$ are the number of basis in $x$ and $y$ coordinate respectively. If a zero pointer is passed, the surface sets all coefficients equal to 0.
4.

Set the tensor product coefficient matrix of B-spline basis functions. Parameter theta is the coefficient matrix.
5.

Return a constant reference of the B-spline basis tensor product coefficients. Coefficient are stored in the surface object in row-wise major.

2.4.4 Evaluation

double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const;

Description:

This function implements the public interface of the FunctionND base class. It returns the B-spline surface value and computes gradient and Hessian if required. See also FunctionND::eval().

2.5 Bi-linear Surface

Class BilinearSurface is a concrete type to represent bi-linear surfaces using the tensor product of two sets of linear basis along $x$ and $y$ . It is derived from the intermediate FunctPiecewise2D class. BilinearSurface class is defined in the header file “bi_linear_surface.h”.

class BilinearSurface : public FunctPiecewise2D

{

};

2.5.1 Knot Point Functions

1 void set_x_even(double x_lo, double x_hi, size_t n);

2 void set_y_even(double y_lo, double y_hi, size_t n);

3 void set_x_customize(const double* x, size_t n);

4 void set_y_customize(const double* y, size_t n);

Description:

1.

Set the knot points along $x$ dimension by evenly breaking the domain $(x_{lo},x_{hi})$ to $n$ pieces.
2.

Set the knot points along $y$ dimension by evenly breaking the domain $(y_{lo},y_{hi})$ to $n$ pieces.
3.

Set the knot points along $x$ dimension with a given sequence of break points. Parameter x is the pointer address of the knot points. Parameter n is the total number of knot points.
4.

Set the knot points along $y$ dimension with a given sequence of break points. Parameter y is the pointer address of the knot points. Parameter n is the total number of knot points.

2.5.2 Basis Functions

1 size_t n_basis_x() const;

2 size_t n_basis_y() const;

3 void set_basis_coeff(const double* _theta = 0);

4 void set_basis_coeff(const Matrix& theta);

5 const ftl::Buffer<double>& basis_coeff() const;

Description:

1.

Return the total number of basis functions along $x$ used by the surface.
2.

Return the total number of basis functions along $y$ used by the surface.
3.

Set the tensor product coefficients of bi-linear basis functions. Parameter _theta is the pointer address to the $m\times n$ basis coefficients, where $m$ and $n$ are the number of basis in $x$ and $y$ coordinate respectively. If a zero pointer is passed, the surface sets all coefficients equal to 0.
4.

Set the tensor product coefficient matrix of bi-linear basis functions. Parameter theta is the coefficient matrix.
5.

Return a constant reference of the bi-linear basis tensor product coefficients. Coefficient are stored in the surface object in row-wise major.

2.5.3 Evaluation

double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const;

Description:

This function implements the public interface of the FunctionND base class. It returns the bi-linear surface value and computes gradient and Hessian if required. See also FunctionND::eval().

2.6 Chebyshev Surface

Class ChebyshevSurface is a concrete type to represent continuous (smooth) surface using the tensor product of two sets of Chebyshev basis polynomials along $x$ and $y$ . It is a directly derived class from the FunctionND base class. Class ChebyshevSurface is defined in header file “chebyshev_surface.h”.

class ChebyshevSurface : public FunctionND

{

};

2.6.1 Constructor

ChebyshevSurface(

double a,

double b,

double c,

double d,

size_t degree_x,

size_t degree_y);

Description:

Constructor, create a Chebyshev surface over with working domain $[a,b]\times[c,d]$ and highest degrees of the polynomial. Parameter a and b defines the domain $[a,b]$ along $x$ . Parameter c and d defines the domain $[c,d]$ along $y$ . Parameter degree_x and degree_y are the highest degrees of Chebyshev basis polynomial along $x$ and $y$ respectively.

2.6.2 Query Functions

1 double x_lo() const;

2 double x_hi() const;

3 double y_lo() const;

4 double y_hi() const;

5 size_t degree_x() const;

6 size_t degree_y() const;

Description:

1.

Return the lower bound of the working domain along $x$ .
2.

Return the upper bound of the working domain along $x$ .
3.

Return the lower bound of the working domain along $y$ .
4.

Return the upper bound of the working domain along $y$ .
5.

Return the highest degree of Chebyshev basis polynomial along $x$ .
6.

Return the highest degree of Chebyshev basis polynomial along $y$ .

2.6.3 Basis Functions

1 void set_basis_coeff(const double* _theta = 0);

2 void set_basis_coeff(const Matrix& theta);

3 const Matrix& basis_coeff() const;

Description:

1.

Set the tensor product coefficients of Chebyshev basis polynomials. Parameter _theta is the pointer address to the $m\times n$ basis coefficients, where $m$ and $n$ are the number of basis in $x$ and $y$ coordinate respectively. In other words, $m-1$ and $n-1$ are the highest degree of Chebyshev basis in $x$ and $y$ respectively. If a zero pointer is passed, the surface sets all coefficients equal to 0.
2.

Set the tensor product coefficients of Chebyshev basis polynomials. Parameter theta is the coefficient matrix.
3.

Return a constant reference of the tensor product basis coefficient matrix stored in the Chebyshev surface object.

2.6.4 Evaluation

double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const;

Description:

This function implements the public interface of the FunctionND base class. It returns the Chebyshev surface value and computes gradient and Hessian if required. See also FunctionND::eval().

2.6.5 Integration

double integral(const Matrix& X_LO, const Matrix& X_HI) const;

Description:

This function implements the public interface of the FunctionND base class. It computes the finite integral with lower and upper limits given in parameter X_LO and X_HI respectively. See also FunctionND::eval().

3 n-d Vector Valued Functions

3.1 VtrValueFunctionND Base Class

VtrValueFunctionND is the pure virtual base class for vector valued multi-dimensional functions. It is defined in header file “function_nd.h”.

3.1.1 Evaluation

1 virtual void eval(

2 const Matrix& x,

3 Matrix* _v = 0,

4 Matrix* _J = 0) const = 0;

Description:

Given a $n$ -dimensional unknown variable $\mathbf{x}$ , this function evaluates $\mathbf{v}=f(\mathbf{x})$ where the result $\mathbf{v}$ is a $m$ -element vector. Parameter x holds the unknown variable $\mathbf{x}$ . Parameter _v is optional. If a non-zero pointer is supplied, it will hold the result of $\mathbf{v}$ on return. Parameter _J is optional. If a non-zero pointer is supplied, it will hold the Jacobian matrix on return.

When users need supply an implementation of the VtrValueFunctionND interface, i.e. a vector-valued residual function for model fitting, they could choose to just compute the value vector $\mathbf{v}$ and call the protected function jacobian_fwd_diff implemented in VtrValueFunctionND to approximate the Jacobian matrix by finite difference. An example of such implementation could be,

1 void eval(const Matrix& x, Matrix* _v=0, Matrix* _J=0) const

2 {

3 if (_v != 0)

4 {

5 // compute value vector v = f(x)

6 }

8 if (_J != 0)

9 {

10 // Finite difference approximate of Jacobian

11 this->jacobian_fwd_diff(x, *_J, _v);

12 }

13 }

3.1.2 Query Function

1 virtual size_t vector_size() const = 0;

2 virtual size_t dimension_size() const = 0;

Description:

1.

Return the number of value elements to be computed in result vector.
2.

Return the dimension size of unknown variables.

3.2 VtrValueFunctionND_M

Class VtrValueFunctionND_M is a concrete type derived from VtrValueFunctionND base class. It holds $m$ number of FunctionND object handles in a buffer and evaluates them simultaneously to produce a $m$ -element result vector. A FunctionND::Handle is typedefed from shared pointer of FunctionND.

struct VtrValueFunctionND_M :

public VtrValueFunctionND,

public ftl::Buffer<FunctionND::Handle>

{

};

4 1-d Interpolation

Class InterpCubicSplineCurve in interp_cubic_spline_curve.h provides a generalized cubic spline interpolation method that can match one, two or all three of $[f(x),f^{\prime}(x),f^{\prime\prime}(x)]$ at a list of data points. This allows users to more effectively control the level, slope and curvature of the interpolated curve.

InterpCubicSplineCurve is a derived class from Function1D. Once fitted to the given interpolation points, it can be used just like a regular 1-d function.

class InterpCubicSplineCurve : public Function1D

{

};

4.1 Definition of Data Points

1 struct Point

2 {

3 enum Type

4 {

5 VALUE,

6 SLOPE,

7 CURVATURE,

8 };

10 double x;

11 double v;

12 Type type;

14 Point() : x(0), v(0), type(VALUE) {}

15 Point(double xx, double vv, Type t) : x(xx),v(vv),type(t) {}

16 };

Interpolation data points are represented by a nested struct, InterpCubicSplineCurve::Point, inside the InterpCubicSplineCurve class. The details of this Point are given below by the line numbers in the struct definition.

•

Line number 3-8, this is the type enumeration indicating whether the point holds a function value $f(x)$ , slope $f^{\prime}(x)$ or curvature $f^{\prime\prime}(x)$ .
•

Line number 10, member variable x holds the $x$ coordinate value of the data point.
•

Line number 11, member variable v holds the value of $f(x)$ , $f^{\prime}(x)$ or $f^{\prime\prime}(x)$ .
•

Line number 12, member variable type indicates the value type of member variable v.
•

Line number 14-15, contructors to initialize a data point object.

4.2 Interpolation Constructor

InterpCubicSplineCurve();

Description:

Create a cublic spline interpolation curve. This curve can be fit to a list of data points at a later time.

4.3 Clone

Function1D* clone() const;

Description:

Implement the interface function Function1D::clone(). It returns a new instance of cubic interpolation curve identical to the current curve.

4.4 Evaluation

double eval(double x, double* _derv1=0, double* _derv2=0) const;

Description:

Implement the interface function Function1D::eval(). It returns the cubic spline function value $f(x)$ and computes $f^{\prime}(x),f^{\prime\prime}(x)$ if non-zero parameters are supplied to _derv1 and _derv2.

4.5 Integration

double integral(double lo, double hi) const;

Description:

Implement the interface function Function1D::integral(). It returns the finite integral of the cubic spline function from lower limit lo to upper limit hi.

4.6 Interpolation

void interp(const ftl::Buffer<Point>& points);

Description:

Interpolate the cubic spline curve through a list of data points $(x_{0},\cdots,x_{n-1})$ and match $f(x_{i})$ , $f^{\prime}(x_{i})$ or $f^{\prime\prime}(x_{i})$ depending on the value type specified in the data point object. $(x_{0},\cdots,x_{n-1})$ doesn’t need be sorted. The buffer can also hold data points with same $x$ but different value types. This function can be repeated called, in which case the curve object is fitted to the new list of input data points.

5 1-d Integration

Class Integrate1D in integrate_1d.h implements an adaptive Gaussian quadrature method to numerically approximate a finite integral $\int_{a}^{b}f(x)dx$ .

class Integrate1D

{

};

5.1 Constructor

1 Integrate(const Function1D& f, double eps)

Description:

The first parameter f is a constant reference of the 1-d integrand function $f(x)$ . The quadrature object internally saves $f$ as a data member.

The second parameter eps is the error tolerance that controls whether the adaptive quadrature needs further refine a sub-interval.

REQUIRED: The 1-d integrand function object needs be alive in the lifecycle of the quadrature.

5.2 Integrand Function

1 const Function1D& integrand() const;

Description:

Return a constant reference to the integrand function. The reference is held as a data member of class Integrate1D.

5.3 Precision Control

1 double precision() const;

2 void set_precision(double eps);

Description:

1.

Return the error tolerance level used in adaptive quadrature method.
2.

Set the error tolerance used in adaptive quadrature method.

5.4 Integrate Operator

1 double operator() (double a, double b) const;

Description:

Return the quadrature approximation of finite integral $\int_{a}^{b}f(x)dx$ . Parameters a and b are the lower and upper limite of the integral.

6 Iterative Method

Most numerical problems are solved iteratively with some convergence requirement. The IterativeMethod class defined in iterative.h is the common base class of iterative solvers implemented in the library.

6.1 Constructor and Destructor

1 protected:

2 explicit IterativeMethod(size_t max_iter = 100);

3 public:

4 virtual ~IterativeMethod() = 0;

Description:

This class is a pure virtual class. It has a protected default constructor and public pure virtual destructor.

6.2 Iteration Control

1 void set_max_iter(size_t max_iter);

2 size_t max_iter() const;

3 size_t actual_iter() const;

4 void enable_exception_on_max_iter(bool flag = true);

Description:

1.

Set the max number of iterations allowed to run.
2.

Return the max number of iterations allowed.
3.

Return the actual number of iterations being used by an solver.
4.

Enable or disable throwing an exception when the max number of iterations is reached.

6.3 Convergence Control

1 void set_converge_tolerance(double t);

2 double converge_tolerance() const;

3 bool is_converged_to_zero(const Matrix& v) const;

4 bool is_converged(const Matrix& v, const Matrix& t) const;

Description:

1.

Set the convergence tolerance. It is upto the actual algorithm to perform convergency test.
2.

Return the convergence tolerance used by a solver.
3.

Return if a vector’s max norm $|\mathbf{v}|_{\infty}$ is smaller than the convergence tolerance.
4.

Return if two vector’s difference $|\mathbf{v}-\mathbf{t}|_{\infty}$ is smaller than the convergence tolerance. REQUIRED: v and t are vectors of same size.

7 Numerical Bounds

NPL represents the numerical boundary of variable $x\in[x_{\text{lo}},x_{\text{hi}}]$ by class Bound in header file bound.h. Across the library, two doubles $a$ and $b$ are considered as numerically equal if $\frac{|a-b|}{1+|b|}<\epsilon$ for some tolerance $\epsilon$ , default to $\epsilon=1.e^{-12}$ . Note that this default tolerance is not the machine precision. The Bound class uses this default tolerance to test if a value is equal to the boundary value.

7.1 Object Construction

1 Bound();

2 Bound(double _lo, double _hi);

3 static Bound lower(double _lo);

4 static Bound upper(double _hi);

Description:

1.

Default constructor, initialized the object as unbounded.
2.

Initialize the object with a given lower and upper bound. REQUIRED: upper bound is strictly greater than or numerically equal to lower bound.
3.

Create a lower bound object.
4.

Create an upper bound object.

7.2 Comparison Operator

1 bool operator == (const Bound& other) const;

2 bool operator != (const Bound& other) const;

Description:

1.

Return if two bound objects have numerically equal lower and upper bounds.
2.

Return if two bound objects do NOT have numerically equal lower and upper bounds.

7.3 Arithmetic Operator

1 Bound& operator += (double v);

2 Bound& operator -= (double v);

3 Bound& operator *= (double v);

4 Bound operator + (double v) const;

5 Bound operator - (double v) const;

6 Bound operator * (double v) const;

Description:

1.

Shift (addition) both lower and upper bounds of the current object by amount v.
2.

Shift (subtraction) both lower and upper bounds of the current object by amount v.
3.

Scale both lower and upper bounds of the current object by factor v.
4.

Return a new object that shifts (addition) both lower and upper bounds of the current object by amount v.
5.

Return a new object that shifts (subtraction) both lower and upper bounds of the current object by amount v.
6.

Return a new object that scales both lower and upper bounds of the current object by amount v.

7.4 Elastic Bounds

1 void set_elastic(bool flag = true);

2 void set_elastic_lo(bool flag = true);

3 void set_elastic_hi(bool flag = true);

4 bool elastic_lo() const;

5 bool elastic_hi() const;

Description:

This set of functions is only used for nonlinear programming at the moment.

1.

Allow both the lower and upper bound to be elastic.
2.

Explicitly set the lower bound to be elastic.
3.

Explicitly set the upper bound to be elastic.
4.

Return the current object has elastic lower bound.
5.

Return the current object has elastic upper bound.

7.5 Query Functions

1 double lo() const;

2 double hi() const;

3 bool has_lower_bound() const;

4 bool has_upper_bound() const;

5 bool has_both_bounds() const;

6 bool has_bounds() const;

7 bool is_equality_bound() const;

8 double violation(double x, double tolerance=EPSILON()) const;

9 Level status(double v, double eps=EPSILON()) const;

Description:

1.

Return the lower bound value.
2.

Return the upper bound value.
3.

Return if the object has a lower bound.
4.

Return if the object has a upper bound.
5.

Return if the object has both lower and upper bounds.
6.

Return if the object has either lower or upper bound.
7.

Return if the lower bound is same as upper bound, effective an equality bound.
8.

Return 0 if a value is in the range between lower and upper bounds; or a positive distance above upper bound; or a negative distance below lower bound.
9.

Return one of the level enumerations to indicate the status of a value relative to the bounds.

Level Enumerations:

1 enum Level

2 {

3 BELOW_LO, // a given value is below its lower bound

4 AT_LO, // a given value is bounded at its lower bound

5 IN_RANGE, // a given value is within bounded range

6 AT_HI, // a given value is bounded at its upper bound

7 ABOVE_HI // a given value is above its upper bound

8 };

8 1-d Root Finding

Class RootSolver1D in root_1d.h derives from IterativeMethod base class. It solves 1-dimensional root finding problem $f(x)=0$ using a safeguarded algorithm that combines bitsection, rational interpolation and Newton method. Please see our white paper for details of the algorithm.

class RootSolver1D : public IterativeMethod

{

};

8.1 Constructor

1 RootSolver1D(const Function1D& f, const Bound& domain=Bound())

Description:

The first parameter f is a constant reference of the 1-d objective function $f(x)$ . The root-finding solver internally saves $f$ as a data member.

REQUIRED: The 1-d function object needs be alive in the lifecycle of the root-finding solver.

The second parameter domain is optional. By default, $f(x)$ has unbounded working domain $x\in(-\infty,\infty)$ . When a bounded domain is provided, the solver will throw an exception when a bracketing operation exceeds the domain.

8.2 Query Functions

1 const Function1D& fx() const;

2 const Bound& domain() const;

Description:

1.

Return the constant reference to the underlying objective function $f(x)$ . This reference is saved as a class data member.
2.

Return the working domain of $f(x)$ .

8.3 Enable or Disable Derivatives $f^{\prime}(x)$

1 void enable_derivative(bool use_df = true);

Description:

Enable or disable the solver to use the derivatives $f^{\prime}(x)$ . Newton method is used to compute the search sequence if $f^{\prime}(x)$ is enabled. Otherwise, rational interpolation is used.

8.4 Solver Operator

1 double operator()(double a, double& x, double b);

Description:

This function computes the root $x^{*}$ such that $|f(x^{*})|<\delta$ , where $\delta$ is set by the set_converge_tolerance() function from base class IterativeMethod.

Parameter a and b is a guess of the bracketing interval $[a,b]$ . If it actually doesn’t bracket, the solver will enter bracketing mode to locate a bracketing interval.

Parameter x holds the initial guess $x_{0}$ on call. REQUIRED: $a<x_{0}<b$ . On return, parameter x holds the solution $x^{*}$ . Meanwhile, $f(x^{*})$ is returned.

9 1-d Minimization

Class MiniSolver1D in mini_1d.h derives from IterativeMethod base class. It solves 1-dimensional minimization problem

\min_{x\in[x_{l},x_{u}]}f(x)

using a safeguarded algorithm that combines polynomial interpolation with Golden section method. Please see our white paper for details of the algorithm.

class MiniSolver1D : public IterativeMethod

{

};

9.1 Constructor

1 MiniSolver1D(const Function1D& f, const Bound& domain=Bound())

Description:

The first parameter f is a constant reference of the 1-d objective function $f(x)$ . The minimization solver internally saves $f$ as a data member.

REQUIRED: The 1-d function object needs be alive in the lifecycle of the minimization solver.

The second parameter domain is optional. By default, $f(x)$ has unbounded working domain $x\in(-\infty,\infty)$ . When a bounded domain is provided, the solver will restrict the search strictly within the specified domain.

9.2 Query Functions

1 const Function1D& objective() const;

2 const Bound& domain() const;

Description:

1.

Return the constant reference to the underlying objective function $f(x)$ . This reference is saved as a class data member.
2.

Return the working domain of $f(x)$ .

9.3 Enable or Disable Derivatives $f^{\prime}(x)$

1 void enable_derivative(bool use_df = true);

Description:

Enable or disable the solver to use the derivatives $f^{\prime}(x)$ . The solver internally fits a quadratic model function $\hat{f}(x)$ to approximate $f(x)$ . If $f^{\prime}(x)$ is enabled, $\hat{f}(x)$ will match $f^{\prime}(u)$ , $f(u)$ and $f(v)$ at two points $u$ and $v$ . When only function values $f(x)$ is used, $\hat{f}(x)$ is fit to match function values at three points.

9.4 Step Length Control

1 void set_max_step(double step);

2 double max_step() const;

Description:

1.

Set the max step length that $x$ is allowed to change at each iteration.
2.

Return the max step that $x$ is allowed to change.

9.5 Solver Operator

1 double operator()(double a, double& x, double b);

Description:

This function computes a local minimum $x^{*}$ for function $f(x)$ . Parameter a and b is a guess of the bracketing interval $[a,b]$ . Parameter x is the initial guess $x_{0}$ on call. On termination, parameter x holds the solution $x^{*}$ and $f(x^{*})$ is returned. REQUIRED: $a<x_{0}<b$ .

The solver will first test if the initial interval $[a,b]$ is indeed a bracketing interval. If not, it will enter bracketing mode to locate a bracketing interval. Once the interval $[a,b]$ is found, the algorithm will enter safeguarded mode to compute $x^{*}\in[a,b]$ . It terminates when the max number of iterations $N$ is reached, or when $|a-b|<\delta$ . Both $N$ and $\delta$ can be set by calling relevant member functions of the base class IterativeMethod.

1-d and n-d Functions Mathwrist C++ API Series

Abstract

1 1-d Functions

1.1 Function1D Base Class

1.1.1 Clone

1.1.2 Evaluation

1.1.3 Integration

1.1.4 Finite Difference Approximation

1.2 1-d Piecewise Function

1.2.1 Piece Lookup

1.2.2 Knot Point Functions

1.2.3 Integration

1.3 Piecewise Constant Function

1.3.1 Constructor

1.3.2 Piecewise Values

1.3.3 Evaluation

1.3.4 Refresh Integral Cache

1.4 Piecewise Linear Function

1.4.1 Constructor

1.4.2 Piecewise Values

1.4.3 Evaluation

1.4.4 Refresh Integral Cache

1.5 Piecewise Linear with Jump Function

1.5.1 Constructor

1.5.2 Piecewise Values

1.5.3 Evaluation

1.5.4 Refresh Integral Cache

1.6 B-spline Curve

1.6.1 Constructor, Destructor and Copy

1.6.2 Knot Point Functions

1.6.3 Basis Functions

1.6.4 Evaluation

1.7 Chebyshev Curve

1.7.1 Constructor

1.7.2 Query Functions

1.7.3 Basis Functions

1.7.4 Evaluation

1.7.5 Integration

1.7.6 Approximation

2 n-d Functions

2.1 FunctionND Base Class

2.1.1 Clone

2.1.2 Evaluation

2.1.3 Integration

2.1.4 Dimension Size

2.1.5 Finite Difference Approximation

2.2 2-d Piecewise Function

2.2.1 Knot Point Function

2.2.2 Integration

2.3 2-d Piecewise Constant Function

2.3.1 Constructors

2.3.2 Piecewise Values

2.3.3 Evaluation

2.4 B-spline Surface

2.4.1 Constructor, Destructor and Copy

2.4.2 Knot Point Functions

2.4.3 Basis Functions

2.4.4 Evaluation

2.5 Bi-linear Surface

2.5.1 Knot Point Functions

2.5.2 Basis Functions

2.5.3 Evaluation

2.6 Chebyshev Surface

2.6.1 Constructor

2.6.2 Query Functions

2.6.3 Basis Functions

2.6.4 Evaluation

2.6.5 Integration

3 n-d Vector Valued Functions

3.1 VtrValueFunctionND Base Class

3.1.1 Evaluation

3.1.2 Query Function

3.2 VtrValueFunctionND_M

4 1-d Interpolation

4.1 Definition of Data Points

4.2 Interpolation Constructor

4.3 Clone

4.4 Evaluation

4.5 Integration

4.6 Interpolation

1-d and n-d Functions
Mathwrist C++ API Series

8.3 Enable or Disable Derivatives $f^{\prime}(x)$

9.3 Enable or Disable Derivatives $f^{\prime}(x)$