Consider the following inverse triangular formula
$$\left( \begin{array}{ccccc} & & & & N_{i-p,p}\left(u_0\right) \\ & & N_{i-2,2}\left(u_0\right) & & \\ & N_{i-1,1}\left(u_0\right) & & & \\ \color{red}{N_{i,0}\left(u_0\right)=1} & & N_{i-1,2}\left(u_0\right) & \cdots & \vdots \\ & N_{i,1}\left(u_0\right) & & & \\ & & N_{i,2}\left(u_0\right) & & \\ & & & & N_{i,p}\left(u_0\right) \end{array} \right)$$
where, $N_{i,0}=1$, and
In addition, $u_0 \in [u_i,u_{i+1})$ knots = $\{u_0,u_1, \cdots, u_m\}, 0\leq u_i \leq u_j$
Here is a procedural implementaion calculate $\color{blue}{N_{i-p,p}(u_0),B_{i-p+1,p}(u_0), \cdots, N_{i,p}(u_0)}$
Search the index $i$ by the auxiliary function
searchSpan
In the code, I use the following local array to store the values of $N_{m,n}$
$$ \left( \begin{array}{ccccc} & & & & N_{i-p,p}\left(u_0\right) \\ & & & & \\ & & N_{i-2,2}\left(u_0\right) & & \vdots \\ & N_{i-1,1}\left(u_0\right) & N_{i-1,2}\left(u_0\right) & \cdots & \\ N_{i,0}\left(u_0\right) & N_{i,1}\left(u_0\right) & N_{i,2}\left(u_0\right) & \cdots & N_{i,p}\left(u_0\right) \end{array} \right)_{(p+1)\times (p+1)} $$
where $N_{m,n}$ was stored in the position $(p+1-i+m,n+1)$ of local array
Search the index of span $[u_i,u_{i+1})$
searchSpan[{deg_, knots_}, u0_] :=
Module[{biSearch},
biSearch =
Function[{low, high},
With[{mid = Floor[(low + high)/2]},
If[u0 < knots[[mid]], {low, mid}, {mid, high}]]
];(*Do bisection search*)
First@
NestWhile[
biSearch[Sequence @@ #, u0] &,
{deg + 1, Length@knots - deg}, Subtract @@ # != -1 &] - 1
]
NonzeroBasis[{deg_, knots_}, u0_] :=
Module[{coeff, basis, i},
coeff =
(u0 - knots[[#1 + 1]])/(knots[[#1 + #2 + 1]] - knots[[#1 + 1]]) &;
basis = ConstantArray[1, {deg + 1, deg + 1}];
i = searchSpan[{deg, knots}, u0];
Do[
basis[[deg + 1 - k, k + 1]] =
(1 - coeff[i - k + 1, k]) basis[[deg + 2 - k, k]];
With[{m = deg + 1 - i},
Do[
basis[[m + j, k + 1]] =
{coeff[j, k], 1 - coeff[j + 1, k]}.{basis[[m + j, k]], basis[[m + j + 1, k]]},
{j, i - k + 1, i - 1}]
];
basis[[deg + 1, k + 1]] =
coeff[i, k] basis[[deg + 1, k]],
{k, deg}];
basis
]
###Test
knots = {0, 0, 0, 0, 0, 1, 2, 3, 4, 4, 5, 6, 7, 8, 9, 10, 10, 10, 10, 10};
deg = 4;
NonzeroBasis[{deg, knots}, 5/2] // MatrixForm
$\left( \begin{array}{ccccc} 1 & 1 & 1 & 1 & \frac{1}{288} \\ 1 & 1 & 1 & \frac{1}{48} & \frac{227}{1152} \\ 1 & 1 & \frac{1}{8} & \frac{23}{48} & \frac{205}{384} \\ 1 & \frac{1}{2} & \frac{3}{4} & \frac{15}{32} & \frac{25}{96} \\ 1 & \frac{1}{2} & \frac{1}{8} & \frac{1}{32} & \frac{1}{192} \end{array} \right)$
###Performance test
knots0 =
Join[ConstantArray[0, 3001], Range[1, 5000], ConstantArray[5001, 3001]];
deg0 = 3000;
NonzeroBasis[{deg0, knots0}, 2.5]; // AbsoluteTiming
###Question
- How to implement this
triangular formulain a non-procedural(like functional or rule-based) method? - How to improve the efficience of
NonzeroBasis?
###Update
Another example I found today was the calculation of Bernstein function
$$B_{n,i}(u)=\binom n i u^i(1-u)^{n-i}$$, where $0 \leq u \leq 1$
In addition, Bernstein function owns the following recursive relationship:
$$B_{n,i}(u)=(1-u) B_{n-1,i}(u)+uB_{n-1,i-1}(u)$$
where, $B_{n,i}(u)=0$ when $i<0$ or $i>n$
So we can use the following triangular schematic digram to calculate $\color{blue}{B_{n,0}(u),B_{n,0}(u), \cdots, B_{n,n}(u)}$
$$\left( \begin{array}{ccccc} \text{} & \text{} & \text{} & \text{} & B_{n,0} (u) \\ \text{} & \text{} & \text{} & .\cdot{}^{\cdot} & \text{} \\ \text{} & \text{} & B_{2,0}(u)& \text{} & \text{} \\ \text{} & B_{1,0} (u) & \text{} & \text{} & \text{} \\ \color{red}{B_{0,0} (u)=1} & \text{} & B_{2,1}(u)& \vdots & \vdots \\ \text{} & B_{1,0}(u) & \text{} & \text{} & \text{} \\ \text{} & \text{} & B_{2,2}(u)& \text{ } & \text{} \\ \text{} & \text{} & \text{} & \ddots & \text{} \\ \text{} & \text{} & \text{} & \text{} & B_{n,n} (u) \\ \end{array} \right)$$
Related question
