Tutorial¶
This section shows how to create polynomial-like approximants to multi-dimensional real-valued functions. A few reminders to the theory follow hereafter; complete information about the methods is available in De Boor’s literature [2][3][4].
Define a function¶
PPPACK objects allow to represent a multivariate real-valued function belonging to a multi-dimensional Cartesian space \(F^{(d)}(V^{(d)})\) as \(f^{(d)}: V^{(d)} \rightarrow \mathbb{R}\), with \(V^{(d)} = V_1 \times \ldots \times V_d\), \(V_i=[a_i,b_i] \in \mathbb{R},\, \forall i\).
\(F^{(d)}\) is constructed by tensor products with \(d\) one-dimensional polynomial-like function spaces, \(F^{(d)} = F_1(V_1) \times \ldots \times F_d(V_d)\), where \(F_i=F_i^{(1)}\). The superscript (d) drops to access the class names of the objects, see the section Class definitions, where for instance \(F^{(d)}\) becomes simply Fd.
Two types of one-dimensional function spaces \(F(=F^{(1)})\) are available:
PP, the piece-wise polynomial space of (polynomial) order \(k\), breakpoints \(\vec{\xi}=[\xi_1,\ldots,\xi_{l+1}]\) and continuities of derivativies \(\vec{\nu}=[\nu_2,\ldots,\nu_{l}]\), \(\mathbb{P}_{k,\vec{\xi}(,\vec{\nu})}\), where \(l\) is the number of polynomial pieces;S, the space constructed by B-splines of (polynomial) order \(k\) and knot sequence \(\vec{t}\), \($_{k,\vec{t}}\).
The two spaces are equivalent if the requirements of the theorem of Curry and Schoenberg (1966) are fulfilled, see PP2S(). Simple polynomial spaces can be obtained as well with \(l=1\), being simply \(\vec{\xi} = [a,b]\).
Provided a base of functions on each dimension, i.e. \(F_j = \operatorname{span} ( \{\psi_{i_j},\,i_j=1,\ldots, I_j\} ),\; \forall j\), it is possible to write the multivariate function \(f\) as:
with \(\vec{i}=[i_1, \ldots, i_d]\) and \(\psi_{\vec{i}} = \prod_j{ \psi_{i_j} }\).
The coefficients \(C_{\vec{i}}\), instead of the heavier notation \(C_{i_1}C_{i_2} \ldots C_{i_d}\), are stored as d-dimensional FORTRAN array, see fd.
Possible one-dimensional (polynomial) functions \(f \in F\) are:
pp\(\;\in \mathbb{P}_{k,\vec{\xi}}\), a piece-wise polynomial function;s\(\;\in $_{k,\vec{t}}\), a function constructed by B-splines.
We use the convention that functions are in lower case letters, whereas function spaces follow in upper case.
The tensor product operator is implemented by overloading the special method mult to build multi-dimensional entities.
| Example: |
|---|
>>> from pppack import *
>>> k,t=2,[1., 1., 2., 3., 4., 4.]
>>> S1 = S(k,t)
>>> xi = [1., 4.]
>>> P1 = PP(k,xi)
>>> F1 = S1 * P1
>>> print(F1)
F^(2) = F0 x F1
with,
F0 = S_{ k=2, t }, with
t = [1.000,1.000,2.000,3.000,4.000,4.000]
F1 = PP_{ k=2, xi, nu }, with
xi = [1.000,4.000], nu = [-,,-]
>>> s0 = s.inS(S1) # = S(k,t) = s(k,t)
>>> print(s0)
s \in S(k=2,t=[1.0, 1.0, 2.0, 3.0, 4.0, 4.0])
coef (to be determined)
---
In addition to inheriting the init method from their space, functions have also another class method for clearer usage (see the example above): s.inS(), pp.inPP() and fd.inFd().
It is possible to assign a name to any object at init; default name is the class name plus an integer counter, increasing internally as far as new instances are created. The creation of multi-dimensional spaces with muliple occurrencies of the same elements is not allowed.
Compute the coefficients¶
The computation of the fitting coefficients requires the introduction of new linear functionals, \(\lambda_{\vec{j}}\), to apply to a target function \(g\) [2]. As well, the functionals apply also to the approximant \(f\) of \(g\), defining the kind of approximation:
by linearity. Interpolation requires as many conditions as unknown coefficients to determine. The functionals can be integral or derivative conditions, or divided differences. The last form covers osculatory interpolation too:
\((,*)\) means optional additional points for the directional divided differences. Simple evaluations of \(g\) at the data points \(x_{j_i}\), now called \(\tau_{j_i}\) as in the references, are
Other kinds of functionals are deferred to later versions.
The inversion of the Gramian matrices building \(\lambda_{\vec{j}}\psi_{\vec{i}}\) and their inner product with \(\lambda_{\vec{j}}g\) (index by index) yields the coefficients \(C_{\vec{i}}\). About the former, inversion is treated separately on the different dimensions: \(\lambda_{\vec{j}}\psi_{\vec{i}} = \lambda_{j_1}\psi_{i_1} \otimes \ldots \otimes \lambda_{j_d}\psi_{i_d}\), see s.cmpcoef() and pp.cmpcoef(). The following inner products are implemented according to De Boor’s intuition in [3], see fd.cmpcoef(), in order to achieve high computational performances.
Evaluate a function¶
Once the coefficients are known, function evaluation is accomplished by the call methods, for all classes of functions. Again, De Boor’s intuition allows to handle smartly the many inner products on successive indices, see fd.__call__(). One-dimensional functions call back the low-level optimized routines of the F90 extension library.
| Example: |
|---|
>>> import numpy as np
>>> g2= lambda x,y: np.sin(x)/(x+.1)+x*np.exp(y)+.25*y
>>> P1.name="linear"
>>> V2=Fd(S1,P1, name="2D")
>>> h2=fd.inFd(V2)
>>> tau=[]
>>> tau.append( S1.cmptau() )
>>> tau.append([1., 4.])
>>> Ni=[len(t) for t in tau]
>>> N=np.prod(Ni)
>>> #
>>> # Use FORTRAN order for continuous elements in memory
>>> # when creating arrays. Double precision arithmetics
>>> # is used as default in the fortran libraries. If simple
>>> # precision is requested, specify dtype=np.float32. However,
>>> # this will cause copies of data due to element type
>>> # mismatch.
>>> gtau=np.zeros((N,), order='F')
>>> for x,i in zip(itt.product(*tau), \
... itt.product(*[range(nx) for nx in Ni])):
... gtau[getidx(i,Ni)]=g2(*x)
>>> h2.cmpcoef(tau,gtau)
>>> # verify interpolation at data sites
>>> for x in itt.product(*tau):
... fx,gx=h2(x),g2(*x)
... print("%+13.6e %+13.6e %+7.2f"%(fx,gx,(1.-fx/gx)*100.))
+3.733255e+00 +3.733255e+00 +0.00
+5.636312e+01 +5.636312e+01 +0.00
+6.119562e+00 +6.119562e+00 +0.00
+1.106293e+02 +1.106293e+02 +0.00
+8.450368e+00 +8.450368e+00 +0.00
+1.648400e+02 +1.648400e+02 +0.00
+1.093854e+01 +1.093854e+01 +0.00
+2.192080e+02 +2.192080e+02 +0.00
>>> print(h2)
fd \in Fd(S(k=2,t=[1.0, 1.0, 2.0, 3.0, 4.0, 4.0]), PP(k=2,xi=[1.0, 4.0],nu=[]))
coef =
[ 3.73325545 6.11956243 8.45036807 10.93854134 17.5432894
34.8365788 52.1298682 69.42315761]
tau = [array([ 1., 2., 3., 4.]), [1.0, 4.0]]
gtau = [ 3.73325545 6.11956243 8.45036807 10.93854134 56.36312366
110.62929884 164.83997268 219.20801416]
---
>>> print("%13.6e"%h2([1.89,3.05]))
7.337242e+01
Array management¶
FORTRAN passes conventionally dummy arguments by reference, that is by passing a pointer to the data. This avoids internal copies and intermediate conversion steps before addressing data, so providing efficient memory management.
This feature is largely used in this package to exchange the multi-dimensional arrays: arrays may appear with indices of different sizes w.r.t. to the calling routine. For instance, a matrix mat(M,N) can be transfered as mat(N,M) to an auxiliary function or subroutine. The ordering of mat changes from mat(m + M(n-1)) to mat(n + N(m-1)) in the called module.
In order to keep the pass-by-reference behavior when C-wrapping, it is necessary to use f2py directives in the FORTRAN code. Here is an example of the chosen implementation strategy for the subroutine spli2d().
subroutine spli2d ( tau, gtau, t, n, k, m, work, q, bcoef, iflag )
! ...
!f2py integer, depend(gtau) :: n = shape(gtau,0)
!f2py integer, check(n*m == shape(gtau,0)*shape(gtau,1)), depend(gtau,n) :: m = shape(gtau,1)
integer ( kind = 4 ), intent(in) :: m, n
!f2py real ( kind = 8 ), intent(in), dimension(*,*) :: gtau
real ( kind = 8 ), intent(in) :: gtau(n,m)
! ...
end subroutine
With the directives above there is basic check of array dimensions, and an exception is raised in case of mismatch of gtau element types or gtau is not contiguous. Pearu Peterson, the main f2py developer, is greatly acknowledged for suggesting this implementation scheme.
Please, remind that copies are avoided only if the dummy argument is FORTRAN contiguous and the element types match exactly. In all other cases, data is copied. So use numpy.ndarray with parameter order=’F’ and dtype=np.float64 for floats and dtype=np.int32 for integers, according to the FORTRAN libraries of pppack (real in double precision with kind = 8, and integer of kind = 4). However, internal type cast is implemented wherever possible to help the users. When array copies are produced by f2py, a message is prompted as a warning.
Deprecated since version 1.0.0: A less efficient alternative is available in fd.cmpcoef() and fd.__call__() with explicit calls to array vectorization and matricization (by vectorize() and matricize() respectively) to implement correctly the calculation of the tensor coefficients and the evaluation of the multi-dimensional functions. Code changes are kept commented as reminder of previous implementation (in pure Python) and internal data organization.
Topical tests¶
Here are a few examples to understand the usage of the package.
Esteban, I’m all ears with the tests.