Functions
AMA_MltvGrdMonoInterp.c File Reference
#include <AMA.h>

Functions

long int AMA_MltvGrdMonoInterp (AMA_OPTIONS *options, long int nind, long int *ng, double **x, double *z, long int *degree, AMA_SPLINE **spline)
 Monotonic Interpolation of Multivariate Gridded Data More...
 

Function Documentation

long int AMA_MltvGrdMonoInterp ( AMA_OPTIONS options,
long int  nind,
long int *  ng,
double **  x,
double *  z,
long int *  degree,
AMA_SPLINE **  spline 
)

Monotonic Interpolation of Multivariate Gridded Data

This function employs cnspla to compute a spline interpolant of independent variable data ${\bf X}\in{\rm R}^n$ and dependent variable data ${\bf Z}\in{\rm R}^1$. The independent variable data is given as $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$ and $1\le k\le n$, and it defines the rectilinear grid

\[ (x_{1,1},x_{1,2},\cdots,x_{1,N^g_1-1},x_{1,N^g_1})\times\cdots\times(x_{n,1},x_{n,2},\cdots,x_{n,N^g_n-1},x_{n,N^g_n}). \]

The dependent variable data lies on the grid and is given as $z_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$; that is, there are $N=\prod_{k=1}^nN^g_k$ dependent variable values.

For a given set of data this function employs cnspla to compute the spline

\[ s({\bf X}) = \sum_{j_n=1}^{m_n}\cdots\sum_{j_1=1}^{m_1}\alpha_{j_1,\ldots,j_n}B_{d_1,j_1}(x_1|{\bf\Lambda}_1)\cdots B_{d_n,j_n}(x_n|{\bf\Lambda_n}) \]

that minimizes

\[ F^{(2)}(s({\bf X})) = \sum_{k=1}^n\left[\int_R\left(\partial^2 s({\bf X})\over\partial^2 x_k\right)^2dR \right] \]

subject to the interpolation constraints

\[ s(x_{\ell_1},\cdots,x_{\ell_n}) = z_{\ell_1,\cdots,\ell_n} \]

for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$. Additionally, for each $1\le k_o \le n$, this function provides the option to impose the local monotonicity constraints

\[ {\partial s({\bf X})\over\partial x_{k_o}} \hspace{5pt} \cases { \ge 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_{k_o}\in[x_{k_o,\ell_{k_o}},x_{k_o,\ell_{k_o}+1}], & if $z_{\ell_1,\cdots,\ell_{k_o},\cdots,\ell_n} < z_{\ell_1,\cdots,\ell_{k_o}+1,\cdots,\ell_n}$; \cr \le 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_{k_o}\in[x_{k_o,\ell_{k_o}},x_{k_o,\ell_{k_o}+1}], & if $z_{\ell_1,\cdots,\ell_{k_o},\cdots,\ell_n} > z_{\ell_1,\cdots,\ell_{k_o}+1,\cdots,\ell_n}$; \cr = 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_{k_o}\in[x_{k_o,\ell_{k_o}},x_{k_o,\ell_{k_o}+1}], & otherwise. \cr } \]

These constraints are imposed at $x_{\ell_1,\cdots,\ell_{k_o},\cdots,\ell_n}$, for all $1\le \ell_k \le N^g_k$ and $1\le k\le n$ with $k\ne k_o$. By default the local monotonicity constraints are imposed on the spline in all of its independent variables but the constraints can be disabled in one or more of its independent variables with AMA_OptionsSetMonotonicity().

In the above definition of $s({\bf X})$ the $\alpha_{j_1,\ldots,j_n}$, for $j_k=1,\ldots,m_k$ and $1\le k\le n$, are the coefficients of the tensor product B-splines

\[ \Phi_{j_1,\ldots,j_n}({\bf X}) = B_{d_1,j_1}(x_1|{\bf\Lambda}_1)\cdots B_{d_n,j_n}(x_n|{\bf\Lambda}_n) \]

where $B_{d_k,j_k}(x_k|{\bf\Lambda}_k)$ are the $m_k$ univariate B-splines of degree $1 \le d_k\le 5$ defined by the knot vector ${\bf\Lambda}_k\in{\rm R}^{M_k}$. The knot vectors ${\bf\Lambda}_k$ are

\[ {\bf\Lambda}_k = ( \lambda_{k,1},\cdots,\lambda_{k,d_k+1}, \lambda_{k,d_k+2}, \cdots, \lambda_{k,m_k}, \lambda_{k,m_k+1},\cdots,\lambda_{k,m_k+d_k+1} )^T \]

and they depend on the independent variable data and the degree $d_k$. If the local monotonicity constraints are requested in the $k$-th independent variable, then the knot vector ${\bf\Lambda}_k$ is defined by AMA_LamdaMonoInterp(); otherwise, it is defined by AMA_LamdaInterp(). In either case the knot vector ${\bf\Lambda}_k$ is based on the points $(x_{k,1},x_{k,2},\cdots,x_{k,N^g_k-1},x_{k,N^g_k})$.

This function computes a spline which satisfies either the full or reduced continuity condition at the interior grid points, as well as, the interpolation and local monotonicity constraints. The following table summarizes the continuity of $s({\bf X})$ with respect to $x_k$ at the interior grid points $x_{k,l_k}$, for $l_k=2,\ldots,n^g_k-1$, and the number of coefficients, $m_k$, for a spline that satisfies either the full or reduced continuity condition. In the table the term

\[ s^{(p)}(x_{k,l_k}) = {\partial^p s(x_1,\cdots,x_{k,l_k},\cdots,x_n)\over \partial x^p_k} \]

is the $p$-th order partial of $s({\bf X})$ with respect to $x_k$ evaluated at $x_{k,l_k}$ and defined over the region $[x_{1,1},x_{1,n^g_1}]\times\cdots\times[x_{k-1,1},x_{k-1,n^g_{k-1}}]\times[x_{k+1,1},x_{k+1,n^g_{k+1}}]\times\cdots\times[x_{n,1},x_{n,n^g_n}]$.

Degree Full Continuity $m_k$ Reduced Continuity $m_k$
Linear $ s^{(p)}(x_{k,l_k})$ is continuous for $p=0$ $n^g_k$ $s^{(p)}(x_{k,l_k})$ is continuous for $p=0$ $n^g_k$
Quadratic $s^{(p)}(x_{k,l_k})$ is continuous for $0\le p\le 1$ $2n^g_k$ $s^{(p)}(x_{k,l_k})$ is continuous for $p=0$ $2n^g_k-1$
Cubic $s^{(p)}(x_{k,l_k})$ is continuous for $0\le p\le 2$ $3n^g_k$ $s^{(p)}(x_{k,l_k})$ is continuous for $0\le p\le 1$ $2n^g_k$
Quartic $s^{(p)}(x_{k,l_k})$ is continuous for $0\le p\le 3$ $4n^g_k$ $s^{(p)}(x_{k,l_k})$ is continuous for $0\le p\le 2$ $3n^g_k$
Quintic $s^{(p)}(x_{k,l_k})$ is continuous for $0\le p\le 4$ $5n^g_k$ $s^{(p)}(x_{k,l_k})$ is continuous for $0\le p\le 3$ $4n^g_k$

This function employs an efficient algorithm for computing monotonic spline interpolants of multivariate gridded data. The algorithm consists of computing the splines

\[ \bar s_k({\bf X}) = \sum_{q_n=1}^{\bar m_n}\cdots\sum_{q_1=1}^{\bar m_1}\bar\alpha^k_{q_1,\ldots,q_n}B_{\bar d_1,q_1}(x_1|{\bf\bar\Lambda}_1)\cdots B_{\bar d_n,q_n}(x_n|{\bf\bar\Lambda_n}) \]

for $1\le k\le n$. The coefficients $\bar\alpha^k_{q_1,\ldots,q_n}$ of $\bar s_k({\bf X})$ are the coefficients of the tensor product B-splines

\[ \Phi_{q_1,\ldots,q_n}({\bf X}) = B_{\bar d_1,q_1}(x_1|{\bf\bar\Lambda}_1)\cdots B_{\bar d_n,q_n}(x_n|{\bf\bar\Lambda}_n) \]

where $B_{\bar d_k,q_k}(x_k|{\bf\bar\Lambda}_k)$ are the $\bar m_k$ univariate B-splines of degree $\bar d_k$ defined by the knot vector ${\bf\bar\Lambda}_k\in{\rm R}^{\bar m_k + \bar d_k + 1}$. The degree $\bar d_{k_o}$ is

\[ \bar d_{k_o} = \cases { d_{k_o}, & for $1\le k_o\le k$, \cr 1 , & for $k < k_o\le n $; \cr } \]

the number of coefficients $\bar m_{k_o}$ is

\[ \bar m_{k_o} = \cases { m_{k_o} , & for $1\le k_o\le k$, \cr n^g_{k_o}, & for $k < k_o\le n $; \cr } \]

and the knot vector ${\bf\bar\Lambda}_{k_o}$ is

\[ {\bf\bar\Lambda}_{k_o} = \cases { {\bf\Lambda}_{k_o} , & for $1\le k_o\le k$, \cr ( x_{k_o,1}, x_{k_o,1}, x_{k_o,2}, \cdots, x_{k_o,n^g_{k_o}-1}, x_{k_o,n^g_{k_o}}, x_{k_o,n^g_{k_o}})^T & for $k < k_o\le n $. \cr } \]

The spline $s({\bf X})$ is based on the spline $\bar s_n({\bf X})$ and is defined by setting its coefficents $\alpha_{j_1,\ldots,j_n}=\bar\alpha^n_{j_1,\ldots,j_n}$ for $j_k=1,\ldots,m_k$ and $1\le k\le n$. The spline $\bar s_n({\bf X})$ is defined through multiple invocations of AMA_UnvMonoInterp() in the following manner.

First, the spline $\bar s_1({\bf X})$ is defined by invoking AMA_UnvMonoInterp() $\prod_{k=2}^nN^g_{k}$ times to compute the univariate spline

\[ \bar s_{q_2,\cdots,q_n}(x_1) = \sum_{q_1=1}^{m_1}\bar\alpha^1_{q_1,\ldots,q_n}B_{d_1,q_1}(x_1|{\bf\Lambda}_1) \]

that minimizes

\[ F^{(2)}(\bar s_{q_2,\cdots,q_n}(x_1)) = \int_{x_{1,1}}^{x_{1,N^g_1}}\left( \bar s^{(2)}_{q_2,\cdots,q_n}(x_1) \right)^2 dx_1 \]

subject to the interpolation constraints

\[ \bar s_{q_2,\cdots,q_n}(x_{1,\ell_1}) = z_{\ell_1,q_2,\cdots,q_n} \]

and the local monotonicity constraints

\[ {\partial s_{q_2,\cdots,q_n}(x_1)\over\partial x_1} \hspace{5pt} \cases { \ge 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_1\in[x_{1,\ell_1},x_{1,\ell_1+1}], & if $z_{\ell_1,q_2,\cdots,q_n} < z_{\ell_1+1,q_2,\cdots,q_n}$; \cr \le 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_{1}\in[x_{1,\ell_1},x_{1,\ell_1+1}], & if $z_{\ell_1,q_2,\cdots,q_n} > z_{\ell_1+1,q_2,\cdots,q_n}$; \cr = 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_1\in[x_{1,\ell_1},x_{1,\ell_1+1}], & otherwise; \cr } \]

for $\ell_1=1,\ldots,N^g_1$. Computing the univariate splines $\bar s_{q_2,\cdots,q_n}(x_1)$, for $q_k=1,\ldots,N^g_k$ and $2\le k\le n$, defines the coefficients $\bar\alpha^1_{q_1,\ldots,q_n}$ of $\bar s_1({\bf X})$.

After computing $\bar s_1({\bf X})$ the splines $\bar s_k({\bf X})$, for $2\le k\le n$, are defined by invoking AMA_UnvMonoInterp() $\prod_{k_o=1}^{k-1}m_{k_o}\prod_{k_o=k+1}^nn^g_{k_o}$ times to compute the univariate spline

\[ \bar s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k) = \sum_{q_k=1}^{m_k}\bar\alpha^k_{q_1,\ldots,q_n}B_{d_k,q_k}(x_k|{\bf\Lambda}_k) \]

that minimizes

\[ F^{(2)}(\bar s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k)) = \int_{x_{k,1}}^{x_{k,n^g_k}}\left( s^{(2)}_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k) \right)^2 dx_k \]

subject to the interpolation constraints

\[ \bar s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_{k,\ell_k}) = \bar\alpha^{k-1}_{q_1,\cdots,q_{k-1},\ell_k,q_{k+1},\cdots,q_n} \]

and the local monotonicity constraints

\[ {\partial s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k)\over\partial x_k} \hspace{5pt} \cases { \ge 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_k\in[x_{k,\ell_k},x_{k,\ell_k+1}], & if $\bar\alpha^{k-1}_{q_1,\cdots,q_{k-1},\ell_k,q_{k+1},\cdots,q_n} < \bar\alpha^{k-1}_{q_1,\cdots,q_{k-1},\ell_k+1,q_{k+1},\cdots,q_n}$; \cr \le 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_{k}\in[x_{k,\ell_k},x_{k,\ell_k+1}], & if $\bar\alpha^{k-1}_{q_1,\cdots,q_{k-1},\ell_k,q_{k+1},\cdots,q_n} > \bar\alpha^{k-1}_{q_1,\cdots,q_{k-1},\ell_k+1,q_{k+1},\cdots,q_n}$; \cr = 0.0 {\rm \hspace{5pt}for\hspace{5pt}all\hspace{5pt}} x_k\in[x_{k,\ell_k},x_{k,\ell_{k+1}}], & otherwise; \cr } \]

for $\ell_k=1,\ldots,N^g_k$. Computing the univariate splines $\bar s_{q_1,\cdots,q_{k-1},q_{k+1},\cdots,q_n}(x_k)$, for $q_{k_o}=1,\ldots,m_{k_o}$ when $1\le k_o\le k-1$ and for $q_{k_o}=1,\ldots,N^g_{k_o}$ when $k+1 \le k_o \le n$, defines the coefficients $\bar\alpha^k_{q_1,\ldots,q_n}$ of $\bar s_k({\bf X})$.

In total this algorithm consists of $\sum_{k=1}^n\left(\prod_{k_o=1}^{k-1}m_{k_o}\prod_{k_o=k+1}^nN^g_{k_o}\right)$ invocations of AMA_UnvMonoInterp() to compute a multivariate spline which satisfies the interpolation and local monotonicity constraints.

This function does the following:

  • Checks the input parameters for validity, see AMA_MltvGrdInputCheck().
  • Defines the spline $s({\bf X})$ based on the knot vectors ${\bf\Lambda}_k$, for $1\le k\le n$.
  • Defines the splines $\bar s_k({\bf X})$, for $1\le k\le n$, through multiple invocations of AMA_UnvMonoInterp().
  • Stores the monotonic interpolant $s({\bf X})=\bar s_n({\bf X})$ into spline.
Note
By default the spline coefficients are initialized to zero but a different initial value for the spline coefficients can be set with AMA_OptionsSetCoefficients().
Because this function employs the univariate spline interpolation function AMA_UnvMonoInterp(), it does not minimizes the cross partial terms of the penalty term utilized by AMA_MltvInterp().

Parameter Note: In the parameters description given below the limits on $k$ are $1\le k\le n$, k = $k-1$ and N $=\prod_{k=1}^nN^g_k$.

Parameters
options[in] Pointer to AMA_OPTIONS. Must be initialized with AMA_Options() prior to calling AMA_MltvGrdMonoInterp().
nind[in] The number of independent variables $n$. Must satisfy $2\le$ nind $\le$ AMA_MXNIND.
ng[in] Array of size nind containing the number of points $N^g_k$ where ng[k] $= N^g_k$. Must satisfy ng[k] $\ge 2$.
x[in] Array of size nind containing arrays of size ng[k] where x[k] contains the independent variable data $x_{k,\ell_k}$, for $\ell_k=1,\ldots,N^g_k$. The values of x[k] must be in strictly ascending order.
z[in] Array of size N containing the dependent variable data $z_{\ell_1,\cdots,\ell_n}$, for $\ell_1=1,\ldots,N^g_1,\cdots,\ell_n=1,\ldots,N^g_n$.
degree[in] Array of size nind containing the degree $d_k$ where degree[k] $= d_k$. Must satisfy $0\le$ degree[k] $\le 5$.
spline[out] Pointer to AMA_SPLINE pointer containing the monotonic spline interpolant. Must satisfy spline $\ne$ NULL.
Returns
Success/Error Code.

User Callable Function - Documented 030215