# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 30.0 `fit` Chapter.
``fit`` - Curve and Surface Fitting
The main aim of this module is to assist you in finding a function which approximates a set of data points.
Typically the data contain random errors, as of experimental measurement, which need to be smoothed out.
To seek an approximation to the data, it is first necessary to specify for the approximating function a mathematical form (a polynomial, for example) which contains a number of unspecified coefficients: the appropriate fitting function then derives for the coefficients the values which provide the best fit of that particular form.
The module deals mainly with curve and surface fitting (i.e., fitting with functions of one and of two variables) when a polynomial or a cubic spline is used as the fitting function, since these cover the most common needs.
However, fitting with other functions and/or more variables can be undertaken by means of general linear or nonlinear functions (some of which are contained in other modules) depending on whether the coefficients in the function occur linearly or nonlinearly.
Cases where a graph rather than a set of data points is given can be treated simply by first reading a suitable set of points from the graph.
The module also contains functions for evaluating, differentiating and integrating polynomial and spline curves and surfaces, once the numerical values of their coefficients have been determined.
There is also a function for computing a Padé approximant of a mathematical function (see `Background to the Problems <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html#background6>`__).
See Also
--------
``naginterfaces.library.examples.fit`` :
This subpackage contains examples for the ``fit`` module.
See also the :ref:`library_fit_ex` subsection.
Functionality Index
-------------------
**Automatic fitting**
with cubic splines: :meth:`dim1_spline_auto`
**Automatic knot placement**
with bicubic splines
data on rectangular mesh: :meth:`dim2_spline_grid`
Data on lines: :meth:`dim2_cheb_lines`
Data on rectangular mesh: :meth:`dim2_spline_grid`
**Differentiation**
of bicubic splines: :meth:`dim2_spline_derivm`
of polynomials: :meth:`dim1_cheb_deriv`
**Evaluation**
at a point
of cubic splines: :meth:`dim1_spline_eval`
of cubic splines and derivatives: :meth:`dim1_spline_deriv`
at vector of points
of bicubic splines at vector of points: :meth:`dim2_spline_evalv`
of :math:`C^1` scattered fit: :meth:`dim2_spline_ts_evalv`
of cubic splines and optionally derivatives: :meth:`dim1_spline_deriv_vector`
of polynomials
in one variable: :meth:`dim1_cheb_eval2`
in one variable (simple interface): :meth:`dim1_cheb_eval`
in two variables: :meth:`dim2_cheb_eval`
of rational functions: :meth:`pade_eval`
on mesh
of bicubic splines: :meth:`dim2_spline_evalm`
of :math:`C^1` scattered fit: :meth:`dim2_spline_ts_evalm`
**Integration**
of cubic splines (definite integral): :meth:`dim1_spline_integ`
of polynomials: :meth:`dim1_cheb_integ`
:math:`l_1` **fit**
with constraints: :meth:`glinc_l1sol`
with general linear function: :meth:`glin_l1sol`
**Least squares curve fit**
with cubic splines: :meth:`dim1_spline_knots`
with polynomials
arbitrary data points: :meth:`dim1_cheb_arb`
selected data points: :meth:`dim1_cheb_glp`
with constraints: :meth:`dim1_cheb_con`
**Least squares surface fit**
with bicubic splines: :meth:`dim2_spline_panel`
with polynomials: :meth:`dim2_cheb_lines`
**Minimax space fit**
with general linear function: :meth:`glin_linf`
with polynomials in one variable: :meth:`dim1_minimax_polynomial`
Padé approximants: :meth:`pade_app`
**Scattered data fit**
bicubic spline: :meth:`dim2_spline_sctr`
:math:`C^1` spline: :meth:`dim2_spline_ts_sctr`
**Service functions**
general option getting function: :meth:`opt_get`
general option setting function: :meth:`opt_set`
Sorting: :meth:`dim2_spline_sort`
For full information please refer to the NAG Library document
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html
"""
# NAG Copyright 2017-2024.
[docs]def dim1_cheb_arb(kplus1, x, y, w):
r"""
``dim1_cheb_arb`` computes weighted least squares polynomial approximations to an arbitrary set of data points.
.. _e02ad-py2-py-doc:
For full information please refer to the NAG Library document for e02ad
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02adf.html
.. _e02ad-py2-py-parameters:
**Parameters**
**kplus1** : int
:math:`k+1`, where :math:`k` is the maximum degree required.
**x** : float, array-like, shape :math:`\left(m\right)`
The values :math:`x_{\textit{r}}` of the independent variable, for :math:`\textit{r} = 1,2,\ldots,m`.
**y** : float, array-like, shape :math:`\left(m\right)`
The values :math:`y_{\textit{r}}` of the dependent variable, for :math:`\textit{r} = 1,2,\ldots,m`.
**w** : float, array-like, shape :math:`\left(m\right)`
The set of weights, :math:`w_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`. For advice on the choice of weights, see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html#background12>`__.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{kplus1}, \mathrm{kplus1}\right)`
The coefficients of :math:`T_{\textit{j}}\left(\bar{x}\right)` in the approximating polynomial of degree :math:`\textit{i}`. :math:`\mathrm{a}[\textit{i},\textit{j}]` contains the coefficient :math:`a_{{\textit{i}+1,\textit{j}+1}}`, for :math:`\textit{j} = 0,1,\ldots,i`, for :math:`\textit{i} = 0,1,\ldots,k`.
**s** : float, ndarray, shape :math:`\left(\mathrm{kplus1}\right)`
:math:`\mathrm{s}[\textit{i}]` contains the root-mean-square residual :math:`s_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,k`, as described in :ref:`Notes <e02ad-py2-py-notes>`. For the interpretation of the values of the :math:`s_{\textit{i}}` and their use in selecting an appropriate degree, see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html#available1>`__.
.. _e02ad-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{w}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{w}[\textit{i}-1] > 0.0`.
(`errno` :math:`2`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[\textit{i}-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[\textit{i}-1]\geq \mathrm{x}[\textit{i}-2]`.
(`errno` :math:`3`)
On entry, all :math:`\mathrm{x}[\textit{I}-1]` have the same value: :math:`\mathrm{x}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{kplus1} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{mdist} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{kplus1}\leq \textit{mdist}`, where mdist is the number of distinct :math:`\mathrm{x}` values.
(`errno` :math:`4`)
On entry, :math:`\mathrm{kplus1} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{kplus1}\leq m`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{kplus1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{kplus1}\geq 1`.
.. _e02ad-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_cheb_arb`` determines least squares polynomial approximations of degrees :math:`0,1,\ldots,k` to the set of data points :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)` with weights :math:`w_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`.
The approximation of degree :math:`i` has the property that it minimizes :math:`\sigma_i` the sum of squares of the weighted residuals :math:`\epsilon_r`, where
.. math::
\epsilon_r = w_r\left(y_r-f_r\right)
and :math:`f_r` is the value of the polynomial of degree :math:`i` at the :math:`r`\ th data point.
Each polynomial is represented in Chebyshev series form with normalized argument :math:`\bar{x}`.
This argument lies in the range :math:`-1` to :math:`+1` and is related to the original variable :math:`x` by the linear transformation
.. math::
\bar{x} = \frac{\left(2x-x_{\mathrm{max}}-x_{\mathrm{min}}\right)}{\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)}\text{.}
Here :math:`x_{\mathrm{max}}` and :math:`x_{\mathrm{min}}` are respectively the largest and smallest values of :math:`x_r`.
The polynomial approximation of degree :math:`i` is represented as
.. math::
\frac{1}{2}a_{{i+1,1}}T_0\left(\bar{x}\right)+a_{{i+1,2}}T_1\left(\bar{x}\right)+a_{{i+1,3}}T_2\left(\bar{x}\right)+ \cdots +a_{{i+1,i+1}}T_i\left(\bar{x}\right)\text{,}
where :math:`T_{\textit{j}}\left(\bar{x}\right)`, for :math:`\textit{j} = 0,1,\ldots,i`, are the Chebyshev polynomials of the first kind of degree :math:`j` with argument :math:`\left(\bar{x}\right)`.
For :math:`i = 0,1,\ldots,k`, the function produces the values of :math:`a_{{i+1,\textit{j}+1}}`, for :math:`\textit{j} = 0,1,\ldots,i`, together with the value of the root-mean-square residual :math:`s_i = \sqrt{\sigma_i/\left(m-i-1\right)}`.
In the case :math:`m = i+1` the function sets the value of :math:`s_i` to zero.
The method employed is due to Forsythe (1957) and is based on the generation of a set of polynomials orthogonal with respect to summation over the normalized dataset.
The extensions due to Clenshaw (1960) to represent these polynomials as well as the approximating polynomials in their Chebyshev series forms are incorporated.
The modifications suggested by Reinsch and Gentleman (see Gentleman (1969)) to the method originally employed by Clenshaw for evaluating the orthogonal polynomials from their Chebyshev series representations are used to give greater numerical stability.
For further details of the algorithm and its use see Cox (1974) and Cox and Hayes (1973).
Subsequent evaluation of the Chebyshev series representations of the polynomial approximations should be carried out using :meth:`dim1_cheb_eval`.
.. _e02ad-py2-py-references:
**References**
Clenshaw, C W, 1960, `Curve fitting with a digital computer`, Comput. J. (2), 170--173
Cox, M G, 1974, `A data-fitting package for the non-specialist user`, Software for Numerical Mathematics, (ed D J Evans), Academic Press
Cox, M G and Hayes, J G, 1973, `Curve fitting: a guide and suite of algorithms for the non-specialist user`, NPL Report NAC26, National Physical Laboratory
Forsythe, G E, 1957, `Generation and use of orthogonal polynomials for data fitting with a digital computer`, J. Soc. Indust. Appl. Math. (5), 74--88
Gentleman, W M, 1969, `An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients`, Comput. J. (12), 160--165
Hayes, J G (ed.), 1970, `Numerical Approximation to Functions and Data`, Athlone Press, London
"""
raise NotImplementedError
[docs]def dim1_cheb_eval(a, xcap):
r"""
``dim1_cheb_eval`` evaluates a polynomial from its Chebyshev series representation.
.. _e02ae-py2-py-doc:
For full information please refer to the NAG Library document for e02ae
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02aef.html
.. _e02ae-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(\textit{nplus1}\right)`
:math:`\mathrm{a}[\textit{i}-1]` must be set to the value of the :math:`\textit{i}`\ th coefficient in the series, for :math:`\textit{i} = 1,2,\ldots,n+1`.
**xcap** : float
:math:`\bar{x}`, the argument at which the polynomial is to be evaluated. It should lie in the range :math:`-1` to :math:`+1`, but a value just outside this range is permitted (see :ref:`Exceptions <e02ae-py2-py-errors>`) to allow for possible rounding errors committed in the transformation from :math:`x` to :math:`\bar{x}` discussed in :ref:`Notes <e02ae-py2-py-notes>`. Provided the recommended form of the transformation is used, a successful exit is thus assured whenever the value of :math:`x` lies in the range :math:`x_{\mathrm{min}}` to :math:`x_{\mathrm{max}}`.
**Returns**
**p** : float
The value of the polynomial.
.. _e02ae-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{xcap} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{EPS} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\left\lvert \mathrm{xcap}\right\rvert \leq 1+4\times \mathrm{EPS}`, where :math:`\mathrm{EPS}` is machine precision.
(`errno` :math:`2`)
On entry, :math:`\textit{nplus1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nplus1}\geq 1`.
.. _e02ae-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_cheb_eval`` evaluates the polynomial
.. math::
\frac{1}{2}a_1T_0\left(\bar{x}\right)+a_2T_1\left(\bar{x}\right)+a_3T_2\left(\bar{x}\right)+ \cdots +a_{{n+1}}T_n\left(\bar{x}\right)
for any value of :math:`\bar{x}` satisfying :math:`-1\leq \bar{x}\leq 1`.
Here :math:`T_j\left(\bar{x}\right)` denotes the Chebyshev polynomial of the first kind of degree :math:`j` with argument :math:`\bar{x}`.
The value of :math:`n` is prescribed by you.
In practice, the variable :math:`\bar{x}` will usually have been obtained from an original variable :math:`x`, where :math:`x_{\mathrm{min}}\leq x\leq x_{\mathrm{max}}` and
.. math::
\bar{x} = \frac{{\left(\left(x-x_{\mathrm{min}}\right)-\left(x_{\mathrm{max}}-x\right)\right)}}{\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)}
Note that this form of the transformation should be used computationally rather than the mathematical equivalent
.. math::
\bar{x} = \frac{\left(2x-x_{\mathrm{min}}-x_{\mathrm{max}}\right)}{\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)}
since the former guarantees that the computed value of :math:`\bar{x}` differs from its true value by at most :math:`4\epsilon`, where :math:`\epsilon` is the machine precision, whereas the latter has no such guarantee.
The method employed is based on the three-term recurrence relation due to Clenshaw (1955), with modifications to give greater numerical stability due to Reinsch and Gentleman (see Gentleman (1969)).
For further details of the algorithm and its use see Cox (1974) and Cox and Hayes (1973).
.. _e02ae-py2-py-references:
**References**
Clenshaw, C W, 1955, `A note on the summation of Chebyshev series`, Math. Tables Aids Comput. (9), 118--120
Cox, M G, 1974, `A data-fitting package for the non-specialist user`, Software for Numerical Mathematics, (ed D J Evans), Academic Press
Cox, M G and Hayes, J G, 1973, `Curve fitting: a guide and suite of algorithms for the non-specialist user`, NPL Report NAC26, National Physical Laboratory
Gentleman, W M, 1969, `An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients`, Comput. J. (12), 160--165
"""
raise NotImplementedError
[docs]def dim1_cheb_glp(f):
r"""
``dim1_cheb_glp`` computes the coefficients of a polynomial, in its Chebyshev series form, which interpolates (passes exactly through) data at a special set of points.
Least squares polynomial approximations can also be obtained.
.. _e02af-py2-py-doc:
For full information please refer to the NAG Library document for e02af
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02aff.html
.. _e02af-py2-py-parameters:
**Parameters**
**f** : float, array-like, shape :math:`\left(\textit{nplus1}\right)`
For :math:`r = 1,2,\ldots,n+1`, :math:`\mathrm{f}[r-1]` must contain :math:`f_r` the value of the dependent variable (ordinate) corresponding to the value
.. math::
\bar{x}_r = \cos\left(\pi \left(r-1\right)/n\right)
of the independent variable (abscissa) :math:`\bar{x}`, or equivalently to the value
.. math::
x\left(r\right) = \frac{1}{2}\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)\cos\left(\pi \left(r-1\right)/n\right)+\frac{1}{2}\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)
of your original variable :math:`x`. Here :math:`x_{\mathrm{max}}` and :math:`x_{\mathrm{min}}` are respectively the upper and lower ends of the range over which you wish to interpolate.
**Returns**
**a** : float, ndarray, shape :math:`\left(\textit{nplus1}\right)`
:math:`\mathrm{a}[\textit{j}-1]` is the coefficient :math:`a_{\textit{j}}` in the interpolating polynomial, for :math:`\textit{j} = 1,2,\ldots,n+1`.
.. _e02af-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{nplus1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nplus1}\geq 2`.
.. _e02af-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_cheb_glp`` computes the coefficients :math:`a_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,n+1`, in the Chebyshev series
.. math::
\frac{1}{2}a_1T_0\left(\bar{x}\right)+a_2T_1\left(\bar{x}\right)+a_3T_2\left(\bar{x}\right)+ \cdots +a_{{n+1}}T_n\left(\bar{x}\right)\text{,}
which interpolates the data :math:`f_r` at the points
.. math::
\bar{x}_r = \cos\left(\left(r-1\right)\pi /n\right)\text{, }\quad r = 1,2,\ldots,n+1\text{.}
Here :math:`T_j\left(\bar{x}\right)` denotes the Chebyshev polynomial of the first kind of degree :math:`j` with argument :math:`\bar{x}`.
The use of these points minimizes the risk of unwanted fluctuations in the polynomial and is recommended when the data abscissae can be chosen by you, e.g., when the data is given as a graph.
For further advantages of this choice of points, see Clenshaw (1962).
In terms of your original variables, :math:`x` say, the values of :math:`x` at which the data :math:`f_r` are to be provided are
.. math::
x_r = \frac{1}{2}\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)\cos\left(\pi \left(r-1\right)/n\right)+\frac{1}{2}\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)\text{, }\quad r = 1,2,\ldots,n+1
where :math:`x_{\mathrm{max}}` and :math:`x_{\mathrm{min}}` are respectively the upper and lower ends of the range of :math:`x` over which you wish to interpolate.
Truncation of the resulting series after the term involving :math:`a_{{i+1}}`, say, yields a least squares approximation to the data.
This approximation, :math:`p\left(\bar{x}\right)`, say, is the polynomial of degree :math:`i` which minimizes
.. math::
\frac{1}{2}\epsilon_1^2+\epsilon_2^2+\epsilon_3^2 + \cdots +\epsilon_n^2+\frac{1}{2}\epsilon_{{n+1}}^2\text{,}
where the residual :math:`\epsilon_{\textit{r}} = p\left(\bar{x}_{\textit{r}}\right)-f_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,n+1`.
The method employed is based on the application of the three-term recurrence relation due to Clenshaw (1955) for the evaluation of the defining expression for the Chebyshev coefficients (see, for example, Clenshaw (1962)).
The modifications to this recurrence relation suggested by Reinsch and Gentleman (see Gentleman (1969)) are used to give greater numerical stability.
For further details of the algorithm and its use see Cox (1974) and Cox and Hayes (1973).
Subsequent evaluation of the computed polynomial, perhaps truncated after an appropriate number of terms, should be carried out using :meth:`dim1_cheb_eval`.
.. _e02af-py2-py-references:
**References**
Clenshaw, C W, 1955, `A note on the summation of Chebyshev series`, Math. Tables Aids Comput. (9), 118--120
Clenshaw, C W, 1962, `Chebyshev Series for Mathematical Functions`, Mathematical tables, HMSO
Cox, M G, 1974, `A data-fitting package for the non-specialist user`, Software for Numerical Mathematics, (ed D J Evans), Academic Press
Cox, M G and Hayes, J G, 1973, `Curve fitting: a guide and suite of algorithms for the non-specialist user`, NPL Report NAC26, National Physical Laboratory
Gentleman, W M, 1969, `An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients`, Comput. J. (12), 160--165
"""
raise NotImplementedError
[docs]def dim1_cheb_con(k, xmin, xmax, x, y, w, xf, yf, ip):
r"""
``dim1_cheb_con`` computes constrained weighted least squares polynomial approximations in Chebyshev series form to an arbitrary set of data points.
The values of the approximations and any number of their derivatives can be specified at selected points.
.. _e02ag-py2-py-doc:
For full information please refer to the NAG Library document for e02ag
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02agf.html
.. _e02ag-py2-py-parameters:
**Parameters**
**k** : int
:math:`k`, the maximum degree required.
**xmin** : float
The lower and upper end points, respectively, of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. Unless there are specific reasons to the contrary, it is recommended that :math:`\mathrm{xmin}` and :math:`\mathrm{xmax}` be set respectively to the lowest and highest value among the :math:`x_r` and :math:`{xf}_r`. This avoids the danger of extrapolation provided there is a constraint point or data point with nonzero weight at each end point.
**xmax** : float
The lower and upper end points, respectively, of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. Unless there are specific reasons to the contrary, it is recommended that :math:`\mathrm{xmin}` and :math:`\mathrm{xmax}` be set respectively to the lowest and highest value among the :math:`x_r` and :math:`{xf}_r`. This avoids the danger of extrapolation provided there is a constraint point or data point with nonzero weight at each end point.
**x** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{x}[\textit{r}-1]` must contain the value :math:`x_{\textit{r}}` of the independent variable at the :math:`\textit{r}`\ th data point, for :math:`\textit{r} = 1,2,\ldots,m`.
**y** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{y}[\textit{r}-1]` must contain :math:`y_{\textit{r}}`, the value of the dependent variable at the :math:`\textit{r}`\ th data point, for :math:`\textit{r} = 1,2,\ldots,m`.
**w** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{w}[\textit{r}-1]` must contain the weight :math:`w_{\textit{r}}` to be applied to the data point :math:`x_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`. For advice on the choice of weights see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html>`__. Negative weights are treated as positive. A zero weight causes the corresponding data point to be ignored. Zero weight should be given to any data point whose :math:`x` and :math:`y` values both coincide with those of a constraint (otherwise the denominators involved in the root mean square residuals :math:`S_i` will be slightly in error).
**xf** : float, array-like, shape :math:`\left(\textit{mf}\right)`
:math:`\mathrm{xf}[\textit{r}-1]` must contain :math:`{xf}_{\textit{r}}`, the value of the independent variable at which a constraint is specified, for :math:`\textit{r} = 1,2,\ldots,\textit{mf}`.
**yf** : float, array-like, shape :math:`\left(\textit{mf}+\mathrm{sum}\left(\mathrm{ip}\right)\right)`
The values which the approximating polynomials and their derivatives are required to take at the points specified in :math:`\mathrm{xf}`. For each value of :math:`\mathrm{xf}[\textit{r}-1]`, :math:`\mathrm{yf}` contains in successive elements the required value of the approximation, its first derivative, second derivative, :math:`\ldots,p_{\textit{r}}`\ th derivative, for :math:`\textit{r} = 1,2,\ldots,{mf}`. Thus the value, :math:`{yf}_s`, which the :math:`k`\ th derivative of each approximation (:math:`k = 0` referring to the approximation itself) is required to take at the point :math:`\mathrm{xf}[r-1]` must be contained in :math:`\mathrm{yf}[s-1]`, where
.. math::
s = r+k+p_1+p_2 + \cdots +p_{{r-1}}\text{,}
where :math:`k = 0,1,\ldots,p_r` and :math:`r = 1,2,\ldots,{mf}`. The derivatives are with respect to the independent variable :math:`x`.
**ip** : int, array-like, shape :math:`\left(\textit{mf}\right)`
:math:`\mathrm{ip}[\textit{r}-1]` must contain :math:`p_{\textit{r}}`, the order of the highest-order derivative specified at :math:`\mathrm{xf}[\textit{r}-1]`, for :math:`\textit{r} = 1,2,\ldots,{mf}`. :math:`p_r = 0` implies that the value of the approximation at :math:`\mathrm{xf}[r-1]` is specified, but not that of any derivative.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{k}+1, \mathrm{k}+1\right)`
:math:`\mathrm{a}[\textit{i},\textit{j}]` contains the coefficient :math:`a_{{\textit{i}\textit{j}}}` in the approximating polynomial of degree :math:`\textit{i}`, for :math:`\textit{j} = 0,1,\ldots,\textit{i}`, for :math:`\textit{i} = \textit{n},\ldots,k`.
**s** : float, ndarray, shape :math:`\left(\mathrm{k}+1\right)`
:math:`\mathrm{s}[\textit{i}]` contains :math:`S_{\textit{i}}`, for :math:`\textit{i} = \textit{n},\ldots,k`, the root mean square residual corresponding to the approximating polynomial of degree :math:`i`. In the case where the number of data points with nonzero weight is equal to :math:`k+1-\textit{n}`, :math:`S_i` is indeterminate: the function sets it to zero. For the interpretation of the values of :math:`S_i` and their use in selecting an appropriate degree, see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html#available1>`__.
**n** : int
Contains the total number of constraint conditions imposed: :math:`\mathrm{n} = \textit{mf}+p_1+p_2 + \cdots +p_{\textit{mf}}`.
**resid** : float, ndarray, shape :math:`\left(m\right)`
Contains weighted residuals of the highest degree of fit determined :math:`\left(k\right)`. The residual at :math:`x_{\textit{r}}` is in element :math:`\mathrm{resid}[\textit{r}-1]`, for :math:`\textit{r} = 1,2,\ldots,m`.
.. _e02ag-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{k}\geq \mathrm{n}`.
(`errno` :math:`1`)
On entry, :math:`\textit{mf} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{mf}\geq 1`.
(`errno` :math:`2`)
On entry, :math:`I = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ip}[I-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ip}[I-1]\geq 0`.
(`errno` :math:`3`)
On entry, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{xf}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{J} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xf}[\textit{J}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xf}[\textit{I}-1]\neq \mathrm{xf}[\textit{J}-1]`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{xf}[\textit{I}-1]` lies outside interval :math:`\left[\mathrm{xmin}, \mathrm{xmax}\right]`: :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{xf}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xmin} < \mathrm{xmax}`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{x}[\textit{I}-1]` lies outside interval :math:`\left[\mathrm{xmin}, \mathrm{xmax}\right]` for some :math:`\textit{I}`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{x}[\textit{I}-1]` lies outside interval :math:`\left[\mathrm{xmin}, \mathrm{xmax}\right]`: :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[\textit{i}-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[\textit{i}-1]\geq \mathrm{x}[\textit{i}-2]`.
(`errno` :math:`6`)
On entry, :math:`\mathrm{k}+1 > m^{{\prime \prime }}+\mathrm{n}`, where :math:`m^{{\prime \prime }}` is the number of data points with nonzero weight and distinct abscissae different from all :math:`\mathrm{xf}`, and :math:`\mathrm{n}` is the total number of constraints: :math:`\mathrm{k}+1 = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`m^{{\prime \prime }} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`7`)
The polynomials :math:`\mu \left(x\right)` and/or :math:`\nu \left(x\right)` cannot be found. The problem is too ill-conditioned.
.. _e02ag-py2-py-notes:
**Notes**
``dim1_cheb_con`` determines least squares polynomial approximations of degrees up to :math:`k` to the set of data points :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)` with weights :math:`w_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`.
The value of :math:`k`, the maximum degree required, is to be prescribed by you.
At each of the values :math:`{xf}_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,{mf}`, of the independent variable :math:`x`, the approximations and their derivatives up to order :math:`p_{\textit{r}}` are constrained to have one of the values :math:`{yf}_{\textit{s}}`, for :math:`\textit{s} = 1,2,\ldots,\textit{n}`, specified by you, where :math:`\textit{n} = {mf}+\sum_{{r = 0}}^{{mf}}p_r`.
The approximation of degree :math:`i` has the property that, subject to the imposed constraints, it minimizes :math:`\sigma_i`, the sum of the squares of the weighted residuals :math:`\epsilon_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`, where
.. math::
\epsilon_r = w_r\left(y_r-f_i\left(x_r\right)\right)
and :math:`f_i\left(x_r\right)` is the value of the polynomial approximation of degree :math:`i` at the :math:`r`\ th data point.
Each polynomial is represented in Chebyshev series form with normalized argument :math:`\bar{x}`.
This argument lies in the range :math:`-1` to :math:`+1` and is related to the original variable :math:`x` by the linear transformation
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)}
where :math:`x_{\mathrm{min}}` and :math:`x_{\mathrm{max}}`, specified by you, are respectively the lower and upper end points of the interval of :math:`x` over which the polynomials are to be defined.
The polynomial approximation of degree :math:`i` can be written as
.. math::
\frac{1}{2}a_{{i,0}}+a_{{i,1}}T_1\left(\bar{x}\right)+ \cdots +a_{{ij}}T_j\left(\bar{x}\right)+ \cdots +a_{{ii}}T_i\left(\bar{x}\right)
where :math:`T_j\left(\bar{x}\right)` is the Chebyshev polynomial of the first kind of degree :math:`j` with argument :math:`\bar{x}`.
For :math:`i = \textit{n},\textit{n}+1,\ldots,k`, the function produces the values of the coefficients :math:`a_{{i\textit{j}}}`, for :math:`\textit{j} = 0,1,\ldots,i`, together with the value of the root mean square residual,
.. math::
S_i = \sqrt{\frac{{\sigma_i}}{\left(m^{\prime }+\textit{n}-i-1\right)}}\text{,}
where :math:`m^{\prime }` is the number of data points with nonzero weight.
Values of the approximations may subsequently be computed using :meth:`dim1_cheb_eval` or :meth:`dim1_cheb_eval2`.
First ``dim1_cheb_con`` determines a polynomial :math:`\mu \left(\bar{x}\right)`, of degree :math:`\textit{n}-1`, which satisfies the given constraints, and a polynomial :math:`\nu \left(\bar{x}\right)`, of degree :math:`\textit{n}`, which has value (or derivative) zero wherever a constrained value (or derivative) is specified.
It then fits :math:`y_{\textit{r}}-\mu \left(x_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`, with polynomials of the required degree in :math:`\bar{x}` each with factor :math:`\nu \left(\bar{x}\right)`.
Finally the coefficients of :math:`\mu \left(\bar{x}\right)` are added to the coefficients of these fits to give the coefficients of the constrained polynomial approximations to the data points :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`.
The method employed is given in Hayes (1970): it is an extension of Forsythe's orthogonal polynomials method (see Forsythe (1957)) as modified by Clenshaw (see Clenshaw (1960)).
.. _e02ag-py2-py-references:
**References**
Clenshaw, C W, 1960, `Curve fitting with a digital computer`, Comput. J. (2), 170--173
Forsythe, G E, 1957, `Generation and use of orthogonal polynomials for data fitting with a digital computer`, J. Soc. Indust. Appl. Math. (5), 74--88
Hayes, J G (ed.), 1970, `Numerical Approximation to Functions and Data`, Athlone Press, London
"""
raise NotImplementedError
[docs]def dim1_cheb_deriv(n, xmin, xmax, a, ia1, iadif1):
r"""
``dim1_cheb_deriv`` determines the coefficients in the Chebyshev series representation of the derivative of a polynomial given in Chebyshev series form.
.. _e02ah-py2-py-doc:
For full information please refer to the NAG Library document for e02ah
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ahf.html
.. _e02ah-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the degree of the given polynomial :math:`p\left(x\right)`.
**xmin** : float
The lower and upper end points respectively of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. The Chebyshev series representation is in terms of the normalized variable :math:`\bar{x}`, where
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}\text{.}
**xmax** : float
The lower and upper end points respectively of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. The Chebyshev series representation is in terms of the normalized variable :math:`\bar{x}`, where
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}\text{.}
**a** : float, array-like, shape :math:`\left(1+\mathrm{n}\times \mathrm{ia1}\right)`
The Chebyshev coefficients of the polynomial :math:`p\left(x\right)`. Specifically, element :math:`\textit{i}\times \mathrm{ia1}` of :math:`\mathrm{a}` must contain the coefficient :math:`a_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,n`. Only these :math:`n+1` elements will be accessed.
Unchanged on exit, but see :math:`\mathrm{adif}`, below.
**ia1** : int
The index increment of :math:`\mathrm{a}`. Most frequently the Chebyshev coefficients are stored in adjacent elements of :math:`\mathrm{a}`, and :math:`\mathrm{ia1}` must be set to :math:`1`. However, if for example, they are stored in :math:`\mathrm{a}[0],\mathrm{a}[3],\mathrm{a}[6],\ldots \text{}`, the value of :math:`\mathrm{ia1}` must be :math:`3`. See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ahf.html#fcomments>`__.
**iadif1** : int
The index increment of :math:`\mathrm{adif}`. Most frequently the Chebyshev coefficients are required in adjacent elements of :math:`\mathrm{adif}`, and :math:`\mathrm{iadif1}` must be set to :math:`1`. However, if, for example, they are to be stored in :math:`\mathrm{adif}[0],\mathrm{adif}[3],\mathrm{adif}[6],\ldots \text{}`, the value of :math:`\mathrm{iadif1}` must be :math:`3`. See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ahf.html#fcomments>`__.
**Returns**
**patm1** : float
The value of :math:`p\left(x_{\mathrm{min}}\right)`. If this value is passed to the integration function :meth:`dim1_cheb_integ` with the coefficients of :math:`q\left(x\right)`, the original polynomial :math:`p\left(x\right)` is recovered, including its constant coefficient.
**adif** : float, ndarray, shape :math:`\left(1+\mathrm{n}\times \mathrm{iadif1}\right)`
The Chebyshev coefficients of the derived polynomial :math:`q\left(x\right)`. (The differentiation is with respect to the variable :math:`x`.) Specifically, element :math:`\textit{i}\times \mathrm{iadif1}` of :math:`\mathrm{adif}` contains the coefficient :math:`\bar{a}_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,n-1`. Additionally, element :math:`n\times \mathrm{iadif1}` is set to zero. A call of the function may have the array name :math:`\mathrm{adif}` the same as :math:`\mathrm{a}`, provided that note is taken of the order in which elements are overwritten, when choosing the starting elements and increments :math:`\mathrm{ia1}` and :math:`\mathrm{iadif1}`, i.e., the coefficients :math:`a_0,a_1,\ldots,a_{{i-1}}` must be intact after coefficient :math:`\bar{a}_i` is stored. In particular, it is possible to overwrite the :math:`a_i` completely by having :math:`\mathrm{ia1} = \mathrm{iadif1}`, and the actual arrays for :math:`\mathrm{a}` and :math:`\mathrm{adif}` identical.
.. _e02ah-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{ladif} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iadif1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ladif} > \mathrm{n}\times \mathrm{iadif1}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iadif1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{iadif1}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{la} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ia1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{la} > \mathrm{n}\times \mathrm{ia1}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ia1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ia1}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xmax} > \mathrm{xmin}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
.. _e02ah-py2-py-notes:
**Notes**
``dim1_cheb_deriv`` forms the polynomial which is the derivative of a given polynomial.
Both the original polynomial and its derivative are represented in Chebyshev series form.
Given the coefficients :math:`a_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,n`, of a polynomial :math:`p\left(x\right)` of degree :math:`n`, where
.. math::
p\left(x\right) = \frac{1}{2}a_0+a_1T_1\left(\bar{x}\right)+ \cdots +a_nT_n\left(\bar{x}\right)
the function returns the coefficients :math:`\bar{a}_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,n-1`, of the polynomial :math:`q\left(x\right)` of degree :math:`n-1`, where
.. math::
q\left(x\right) = \frac{{dp\left(x\right)}}{{dx}} = \frac{1}{2}\bar{a}_0+\bar{a}_1T_1\left(\bar{x}\right)+ \cdots +\bar{a}_{{n-1}}T_{{n-1}}\left(\bar{x}\right)\text{.}
Here :math:`T_j\left(\bar{x}\right)` denotes the Chebyshev polynomial of the first kind of degree :math:`j` with argument :math:`\bar{x}`.
It is assumed that the normalized variable :math:`\bar{x}` in the interval :math:`\left[-1, +1\right]` was obtained from your original variable :math:`x` in the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]` by the linear transformation
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}
and that you require the derivative to be with respect to the variable :math:`x`.
If the derivative with respect to :math:`\bar{x}` is required, set :math:`x_{\mathrm{max}} = 1` and :math:`x_{\mathrm{min}} = -1`.
Values of the derivative can subsequently be computed, from the coefficients obtained, by using :meth:`dim1_cheb_eval2`.
The method employed is that of Chebyshev series (see Module 8 of Modern Computing Methods (1961)), modified to obtain the derivative with respect to :math:`x`.
Initially setting :math:`\bar{a}_{{n+1}} = \bar{a}_n = 0`, the function forms successively
.. math::
\bar{a}_{{i-1}} = \bar{a}_{{i+1}}+\frac{2}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}2ia_i\text{, }\quad i = n,n-1,\ldots,1\text{.}
.. _e02ah-py2-py-references:
**References**
Modern Computing Methods, 1961, `Chebyshev-series`, NPL Notes on Applied Science (16), (2nd Edition), HMSO
"""
raise NotImplementedError
[docs]def dim1_cheb_integ(n, xmin, xmax, a, ia1, qatm1, iaint1):
r"""
``dim1_cheb_integ`` determines the coefficients in the Chebyshev series representation of the indefinite integral of a polynomial given in Chebyshev series form.
.. _e02aj-py2-py-doc:
For full information please refer to the NAG Library document for e02aj
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ajf.html
.. _e02aj-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the degree of the given polynomial :math:`p\left(x\right)`.
**xmin** : float
The lower and upper end points respectively of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. The Chebyshev series representation is in terms of the normalized variable :math:`\bar{x}`, where
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}\text{.}
**xmax** : float
The lower and upper end points respectively of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. The Chebyshev series representation is in terms of the normalized variable :math:`\bar{x}`, where
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}\text{.}
**a** : float, array-like, shape :math:`\left(1+\mathrm{n}\times \mathrm{ia1}\right)`
The Chebyshev coefficients of the polynomial :math:`p\left(x\right)`. Specifically, element :math:`\textit{i}\times \mathrm{ia1}` of :math:`\mathrm{a}` must contain the coefficient :math:`a_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,n`. Only these :math:`n+1` elements will be accessed.
Unchanged on exit, but see :math:`\mathrm{aintc}`, below.
**ia1** : int
The index increment of :math:`\mathrm{a}`. Most frequently the Chebyshev coefficients are stored in adjacent elements of :math:`\mathrm{a}`, and :math:`\mathrm{ia1}` must be set to :math:`1`. However, if for example, they are stored in :math:`\mathrm{a}[0],\mathrm{a}[3],\mathrm{a}[6],\ldots \text{}`, the value of :math:`\mathrm{ia1}` must be :math:`3`. See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ajf.html#fcomments>`__.
**qatm1** : float
The value that the integrated polynomial is required to have at the lower end point of its interval of definition, i.e., at :math:`\bar{x} = -1` which corresponds to :math:`x = x_{\mathrm{min}}`. Thus, :math:`\mathrm{qatm1}` is a constant of integration and will normally be set to zero by you.
**iaint1** : int
The index increment of :math:`\mathrm{aintc}`. Most frequently the Chebyshev coefficients are required in adjacent elements of :math:`\mathrm{aintc}`, and :math:`\mathrm{iaint1}` must be set to :math:`1`. However, if, for example, they are to be stored in :math:`\mathrm{aintc}[0],\mathrm{aintc}[3],\mathrm{aintc}[6],\ldots \text{}`, the value of :math:`\mathrm{iaint1}` must be :math:`3`. See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ajf.html#fcomments>`__.
**Returns**
**aintc** : float, ndarray, shape :math:`\left(1+\left(\mathrm{n}+1\right)\times \mathrm{iaint1}\right)`
The Chebyshev coefficients of the integral :math:`q\left(x\right)`. (The integration is with respect to the variable :math:`x`, and the constant coefficient is chosen so that :math:`q\left(x_{\mathrm{min}}\right)` equals :math:`\mathrm{qatm1}`). Specifically, element :math:`i\times \mathrm{iaint1}` of :math:`\mathrm{aintc}` contains the coefficient :math:`a_{\textit{i}}^{\prime }`, for :math:`\textit{i} = 0,1,\ldots,n+1`. A call of the function may have the array name :math:`\mathrm{aintc}` the same as :math:`\mathrm{a}`, provided that note is taken of the order in which elements are overwritten when choosing starting elements and increments :math:`\mathrm{ia1}` and :math:`\mathrm{iaint1}`: i.e., the coefficients, :math:`a_0,a_1,\ldots,a_{{i-2}}` must be intact after coefficient :math:`a_i^{\prime }` is stored. In particular it is possible to overwrite the :math:`a_i` entirely by having :math:`\mathrm{ia1} = \mathrm{iaint1}`, and the actual array for :math:`\mathrm{a}` and :math:`\mathrm{aintc}` identical.
.. _e02aj-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{laint} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iaint1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{laint} > \left(\mathrm{n}+1\right)\times \mathrm{iaint1}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iaint1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{iaint1}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{la} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ia1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{la} > \mathrm{n}\times \mathrm{ia1}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ia1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ia1}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xmax} > \mathrm{xmin}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 0`.
.. _e02aj-py2-py-notes:
**Notes**
``dim1_cheb_integ`` forms the polynomial which is the indefinite integral of a given polynomial.
Both the original polynomial and its integral are represented in Chebyshev series form.
If supplied with the coefficients :math:`a_i`, for :math:`\textit{i} = 0,1,\ldots,n`, of a polynomial :math:`p\left(x\right)` of degree :math:`n`, where
.. math::
p\left(x\right) = \frac{1}{2}a_0+a_1T_1\left(\bar{x}\right)+ \cdots +a_nT_n\left(\bar{x}\right)\text{,}
the function returns the coefficients :math:`a_i^{\prime }`, for :math:`\textit{i} = 0,1,\ldots,n+1`, of the polynomial :math:`q\left(x\right)` of degree :math:`n+1`, where
.. math::
q\left(x\right) = \frac{1}{2}a_0^{\prime }+a_1^{\prime }T_1\left(\bar{x}\right)+ \cdots +a_{{n+1}}^{\prime }T_{{n+1}}\left(\bar{x}\right)\text{,}
and
.. math::
q\left(x\right) = \int p\left(x\right){dx}\text{.}
Here :math:`T_j\left(\bar{x}\right)` denotes the Chebyshev polynomial of the first kind of degree :math:`j` with argument :math:`\bar{x}`.
It is assumed that the normalized variable :math:`\bar{x}` in the interval :math:`\left[-1, +1\right]` was obtained from your original variable :math:`x` in the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]` by the linear transformation
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}
and that you require the integral to be with respect to the variable :math:`x`.
If the integral with respect to :math:`\bar{x}` is required, set :math:`x_{\mathrm{max}} = 1` and :math:`x_{\mathrm{min}} = -1`.
Values of the integral can subsequently be computed, from the coefficients obtained, by using :meth:`dim1_cheb_eval2`.
The method employed is that of Chebyshev series (see Module 8 of Modern Computing Methods (1961)), modified for integrating with respect to :math:`x`.
Initially taking :math:`a_{{n+1}} = a_{{n+2}} = 0`, the function forms successively
.. math::
a_i^{\prime } = \frac{{a_{{i-1}}-a_{{i+1}}}}{{2i}}\times \frac{{x_{\mathrm{max}}-x_{\mathrm{min}}}}{2}\text{, }\quad i = n+1,n,\ldots,1\text{.}
The constant coefficient :math:`a_0^{\prime }` is chosen so that :math:`q\left(x\right)` is equal to a specified value, :math:`\mathrm{qatm1}`, at the lower end point of the interval on which it is defined, i.e., :math:`\bar{x} = -1`, which corresponds to :math:`x = x_{\mathrm{min}}`.
.. _e02aj-py2-py-references:
**References**
Modern Computing Methods, 1961, `Chebyshev-series`, NPL Notes on Applied Science (16), (2nd Edition), HMSO
"""
raise NotImplementedError
[docs]def dim1_cheb_eval2(n, xmin, xmax, a, ia1, x):
r"""
``dim1_cheb_eval2`` evaluates a polynomial from its Chebyshev series representation, allowing an arbitrary index increment for accessing the array of coefficients.
.. _e02ak-py2-py-doc:
For full information please refer to the NAG Library document for e02ak
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02akf.html
.. _e02ak-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the degree of the given polynomial :math:`p\left(\bar{x}\right)`.
**xmin** : float
The lower and upper end points respectively of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. The Chebyshev series representation is in terms of the normalized variable :math:`\bar{x}`, where
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}\text{.}
**xmax** : float
The lower and upper end points respectively of the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]`. The Chebyshev series representation is in terms of the normalized variable :math:`\bar{x}`, where
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}\text{.}
**a** : float, array-like, shape :math:`\left(\mathrm{n}\times \mathrm{ia1}+1\right)`
The Chebyshev coefficients of the polynomial :math:`p\left(\bar{x}\right)`. Specifically, element :math:`\textit{i}\times \mathrm{ia1}` must contain the coefficient :math:`a_{\textit{i}}`, for :math:`\textit{i} = 0,1,\ldots,n`. Only these :math:`n+1` elements will be accessed.
**ia1** : int
The index increment of :math:`\mathrm{a}`. Most frequently, the Chebyshev coefficients are stored in adjacent elements of :math:`\mathrm{a}`, and :math:`\mathrm{ia1}` must be set to :math:`1`. However, if, for example, they are stored in :math:`\mathrm{a}[0],\mathrm{a}[3],\mathrm{a}[6],\ldots \text{}`, the value of :math:`\mathrm{ia1}` must be :math:`3`.
**x** : float
The argument :math:`x` at which the polynomial is to be evaluated.
**Returns**
**result** : float
The value of the polynomial :math:`p\left(\bar{x}\right)`.
.. _e02ak-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xmax} > \mathrm{xmin}`.
(`errno` :math:`1`)
On entry, :math:`\textit{la} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ia1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{la} > \mathrm{n}\times \mathrm{ia1}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ia1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ia1}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{n}+1 = \langle\mathit{\boldsymbol{value}}\rangle`>.
Constraint: :math:`\mathrm{n}\geq 0`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{x}` does not lie in :math:`\left[\mathrm{xmin}, \mathrm{xmax}\right]`: :math:`\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _e02ak-py2-py-notes:
**Notes**
If supplied with the coefficients :math:`a_i`, for :math:`\textit{i} = 0,1,\ldots,n`, of a polynomial :math:`p\left(\bar{x}\right)` of degree :math:`n`, where
.. math::
p\left(\bar{x}\right) = \frac{1}{2}a_0+a_1T_1\left(\bar{x}\right)+ \cdots +a_nT_n\left(\bar{x}\right)\text{,}
``dim1_cheb_eval2`` returns the value of :math:`p\left(\bar{x}\right)` at a user-specified value of the variable :math:`x`.
Here :math:`T_j\left(\bar{x}\right)` denotes the Chebyshev polynomial of the first kind of degree :math:`j` with argument :math:`\bar{x}`.
It is assumed that the independent variable :math:`\bar{x}` in the interval :math:`\left[-1, +1\right]` was obtained from your original variable :math:`x` in the interval :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]` by the linear transformation
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{{x_{\mathrm{max}}-x_{\mathrm{min}}}}\text{.}
The coefficients :math:`a_i` may be supplied in the array :math:`\mathrm{a}`, with any increment between the indices of array elements which contain successive coefficients.
This enables the function to be used in surface fitting and other applications, in which the array might have two or more dimensions.
The method employed is based on the three-term recurrence relation due to Clenshaw (see Clenshaw (1955)), with modifications due to Reinsch and Gentleman (see Gentleman (1969)).
For further details of the algorithm and its use see Cox (1973) and Cox and Hayes (1973).
.. _e02ak-py2-py-references:
**References**
Clenshaw, C W, 1955, `A note on the summation of Chebyshev series`, Math. Tables Aids Comput. (9), 118--120
Cox, M G, 1973, `A data-fitting package for the non-specialist user`, NPL Report NAC 40, National Physical Laboratory
Cox, M G and Hayes, J G, 1973, `Curve fitting: a guide and suite of algorithms for the non-specialist user`, NPL Report NAC26, National Physical Laboratory
Gentleman, W M, 1969, `An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients`, Comput. J. (12), 160--165
"""
raise NotImplementedError
[docs]def dim1_minimax_polynomial(x, y, m):
r"""
``dim1_minimax_polynomial`` calculates a minimax polynomial fit to a set of data points.
.. _e02al-py2-py-doc:
For full information please refer to the NAG Library document for e02al
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02alf.html
.. _e02al-py2-py-parameters:
**Parameters**
**x** : float, array-like, shape :math:`\left(n\right)`
The values of the :math:`x` coordinates, :math:`x_i`, for :math:`\textit{i} = 1,2,\ldots,n`.
**y** : float, array-like, shape :math:`\left(n\right)`
The values of the :math:`y` coordinates, :math:`y_i`, for :math:`\textit{i} = 1,2,\ldots,n`.
**m** : int
:math:`m`, where :math:`m` is the degree of the polynomial to be found.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{m}+1\right)`
The coefficients :math:`a_i` of the minimax polynomial, for :math:`\textit{i} = 0,1,\ldots,m`.
**ref** : float
The final reference deviation, i.e., the maximum deviation of the computed polynomial evaluated at :math:`x_{\textit{i}}` from the reference values :math:`y_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n`. :math:`\mathrm{ref}` may return a negative value which indicates that the algorithm started to cycle due to round-off errors.
.. _e02al-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 1`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}\geq 0`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m} < 100`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m} < n-1`.
(`errno` :math:`3`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{i}] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[\textit{i}] > \mathrm{x}[\textit{i}-1]`.
.. _e02al-py2-py-notes:
**Notes**
Given a set of data points :math:`\left(x_i, y_i\right)`, for :math:`\textit{i} = 1,2,\ldots,n`, ``dim1_minimax_polynomial`` uses the exchange algorithm to compute an :math:`m`\ th-degree polynomial
.. math::
\mathrm{P}\left(x\right) = a_0+a_1x+a_2x^2 + \cdots +a_{{m}}x^m
such that :math:`\mathrm{max}_i\left\lvert \mathrm{P}\left(x_i\right)-y_i\right\rvert` is a minimum.
The function also returns a number whose absolute value is the final reference deviation (see :ref:`Parameters <e02al-py2-py-parameters>`).
The function is an adaptation of Boothroyd (1967).
.. _e02al-py2-py-references:
**References**
Boothroyd, J B, 1967, `Algorithm 318`, Comm. ACM (10), 801
Stieffel, E, 1959, `Numerical methods of Tchebycheff approximation`, On Numerical Approximation, (ed R E Langer), 217--232, University of Wisconsin Press
"""
raise NotImplementedError
[docs]def dim1_spline_knots(x, y, w, lamda):
r"""
``dim1_spline_knots`` computes a weighted least squares approximation to an arbitrary set of data points by a cubic spline with knots prescribed by you.
Cubic spline interpolation can also be carried out.
.. _e02ba-py2-py-doc:
For full information please refer to the NAG Library document for e02ba
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02baf.html
.. _e02ba-py2-py-parameters:
**Parameters**
**x** : float, array-like, shape :math:`\left(m\right)`
The values :math:`x_{\textit{r}}` of the independent variable (abscissa), for :math:`\textit{r} = 1,2,\ldots,m`. The values must satisfy the Schoenberg--Whitney conditions (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02baf.html#fcomments>`__).
**y** : float, array-like, shape :math:`\left(m\right)`
The values :math:`y_{\textit{r}}` of the dependent variable (ordinate), for :math:`\textit{r} = 1,2,\ldots,m`.
**w** : float, array-like, shape :math:`\left(m\right)`
The values :math:`w_{\textit{r}}` of the weights, for :math:`\textit{r} = 1,2,\ldots,m`. For advice on the choice of weights, see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html>`__.
**lamda** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
:math:`\mathrm{lamda}[\textit{i}-1]` must be set to the :math:`\left(\textit{i}-4\right)`\ th (interior) knot, :math:`\lambda_{\textit{i}}`, for :math:`\textit{i} = 5,6,\ldots,\bar{n}+3`.
**Returns**
**lamda** : float, ndarray, shape :math:`\left(\textit{ncap7}\right)`
The input values are unchanged, and :math:`\mathrm{lamda}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,4`, :math:`\textit{ncap7}-3`, :math:`\textit{ncap7}-2`, :math:`\textit{ncap7}-1`, :math:`\textit{ncap7}` contains the additional (exterior) knots introduced by the function. For advice on the choice of knots, see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html#available3>`__.
**c** : float, ndarray, shape :math:`\left(\textit{ncap7}\right)`
The coefficient :math:`c_{\textit{i}}` of the B-spline :math:`N_{\textit{i}}\left(x\right)`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3`. The remaining elements of the array are not used.
**ss** : float
The residual sum of squares, :math:`\theta`.
.. _e02ba-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{J} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{lamda}[\textit{J}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{J}] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[\textit{J}-1]\leq \mathrm{lamda}[\textit{J}]`.
(`errno` :math:`1`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{lamda}[\textit{ncap7}-5] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[m-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[\textit{ncap7}-5] < \mathrm{x}[m-1]`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{lamda}[4] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[4] > \mathrm{x}[0]`.
(`errno` :math:`2`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{w}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{w}[\textit{i}-1] > 0.0`.
(`errno` :math:`3`)
On entry, the :math:`\mathrm{x}` values are not in nondecreasing order. :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{J} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{xdist}[\textit{J}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[\textit{I}-1]\geq \textit{xdist}[\textit{J}-1]`, where xdist is the set of distinct :math:`x`-values.
(`errno` :math:`4`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{mdist} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncap7}\leq \textit{mdist}+4`, where mdist is the number of distinct x-values.
(`errno` :math:`4`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncap7}\leq m+4`.
(`errno` :math:`4`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncap7}\geq 8`.
(`errno` :math:`5`)
On entry, the Schoenberg--Whitney conditions fail to hold for at least one subset of the distinct data abscissae. :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{xdist}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{J} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{J}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{xdist}[\textit{I}-1] < \mathrm{lamda}[\textit{J}-1]`, where xdist is the set of distinct :math:`x`-values.
(`errno` :math:`5`)
On entry, the Schoenberg--Whitney conditions fail to hold for at least one subset of the distinct data abscissae. :math:`L = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{xdist}[L-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{xdist}[\textit{L}-1] > \mathrm{lamda}[\textit{I}-1]`, where xdist is the set of distinct :math:`x`-values.
(`errno` :math:`5`)
On entry, the Schoenberg--Whitney conditions fail to hold for at least one subset of the distinct data abscissae. :math:`\textit{J} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{xdist}[\textit{J}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{J}+4 = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{J}+3] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{xdist}[\textit{J}-1] < \mathrm{lamda}[\textit{J}+3]`, where xdist is the set of distinct :math:`x`-values.
.. _e02ba-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_spline_knots`` determines a least squares cubic spline approximation :math:`s\left(x\right)` to the set of data points :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)` with weights :math:`w_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`.
The value of :math:`\textit{ncap7} = \bar{n}+7`, where :math:`\bar{n}` is the number of intervals of the spline (one greater than the number of interior knots), and the values of the knots :math:`\lambda_5,\lambda_6,\ldots,\lambda_{{\bar{n}+3}}`, interior to the data interval, are prescribed by you.
:math:`s\left(x\right)` has the property that it minimizes :math:`\theta`, the sum of squares of the weighted residuals :math:`\epsilon_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`, where
.. math::
\epsilon_r = w_r\left(y_r-s\left(x_r\right)\right)\text{.}
The function produces this minimizing value of :math:`\theta` and the coefficients :math:`c_1,c_2,\ldots,c_q`, where :math:`q = \bar{n}+3`, in the B-spline representation
.. math::
s\left(x\right) = \sum_{{i = 1}}^qc_iN_i\left(x\right)\text{.}
Here :math:`N_i\left(x\right)` denotes the normalized B-spline of degree :math:`3` defined upon the knots :math:`\lambda_i,\lambda_{{i+1}},\ldots,\lambda_{{i+4}}`.
In order to define the full set of B-splines required, eight additional knots :math:`\lambda_1,\lambda_2,\lambda_3,\lambda_4` and :math:`\lambda_{{\bar{n}+4}},\lambda_{{\bar{n}+5}},\lambda_{{\bar{n}+6}},\lambda_{{\bar{n}+7}}` are inserted automatically by the function.
The first four of these are set equal to the smallest :math:`x_r` and the last four to the largest :math:`x_r`.
The representation of :math:`s\left(x\right)` in terms of B-splines is the most compact form possible in that only :math:`\bar{n}+3` coefficients, in addition to the :math:`\bar{n}+7` knots, fully define :math:`s\left(x\right)`.
The method employed involves forming and then computing the least squares solution of a set of :math:`m` linear equations in the coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,{\bar{n}+3}`.
The equations are formed using a recurrence relation for B-splines that is unconditionally stable (see Cox (1972) and de Boor (1972)), even for multiple (coincident) knots.
The least squares solution is also obtained in a stable manner by using orthogonal transformations, viz. a variant of Givens rotations (see Gentleman (1974) and Gentleman (1973)).
This requires only one equation to be stored at a time.
Full advantage is taken of the structure of the equations, there being at most four nonzero values of :math:`N_i\left(x\right)` for any value of :math:`x` and hence at most four coefficients in each equation.
For further details of the algorithm and its use see Cox (1974), Cox (1975) and Cox and Hayes (1973).
Subsequent evaluation of :math:`s\left(x\right)` from its B-spline representation may be carried out using :meth:`dim1_spline_eval`.
If derivatives of :math:`s\left(x\right)` are also required, :meth:`dim1_spline_deriv` may be used. :meth:`dim1_spline_integ` can be used to compute the definite integral of :math:`s\left(x\right)`.
.. _e02ba-py2-py-references:
**References**
Cox, M G, 1972, `The numerical evaluation of B-splines`, J. Inst. Math. Appl. (10), 134--149
Cox, M G, 1974, `A data-fitting package for the non-specialist user`, Software for Numerical Mathematics, (ed D J Evans), Academic Press
Cox, M G, 1975, `Numerical methods for the interpolation and approximation of data by spline functions`, PhD Thesis, City University, London
Cox, M G and Hayes, J G, 1973, `Curve fitting: a guide and suite of algorithms for the non-specialist user`, NPL Report NAC26, National Physical Laboratory
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
Gentleman, W M, 1973, `Least squares computations by Givens transformations without square roots`, J. Inst. Math. Applic. (12), 329--336
Gentleman, W M, 1974, `Algorithm AS 75. Basic procedures for large sparse or weighted linear least squares problems`, Appl. Statist. (23), 448--454
Schoenberg, I J and Whitney, A, 1953, `On Polya frequency functions III`, Trans. Amer. Math. Soc. (74), 246--259
"""
raise NotImplementedError
[docs]def dim1_spline_eval(lamda, c, x):
r"""
``dim1_spline_eval`` evaluates a cubic spline from its B-spline representation.
.. _e02bb-py2-py-doc:
For full information please refer to the NAG Library document for e02bb
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bbf.html
.. _e02bb-py2-py-parameters:
**Parameters**
**lamda** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
:math:`\mathrm{lamda}[\textit{j}-1]` must be set to the value of the :math:`\textit{j}`\ th member of the complete set of knots, :math:`\lambda_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,\bar{n}+7`.
**c** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
The coefficient :math:`c_{\textit{i}}` of the B-spline :math:`N_{\textit{i}}\left(x\right)`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3`. The remaining elements of the array are not referenced.
**x** : float
The argument :math:`x` at which the cubic spline is to be evaluated.
**Returns**
**s** : float
The value of the spline, :math:`s\left(x\right)`.
.. _e02bb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{ncap7}-4] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}\leq \mathrm{lamda}[\textit{ncap7}-4]`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[3] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}\geq \mathrm{lamda}[3]`.
(`errno` :math:`2`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncap7}\geq 8`.
.. _e02bb-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_spline_eval`` evaluates the cubic spline :math:`s\left(x\right)` at a prescribed argument :math:`x` from its augmented knot set :math:`\lambda_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n+7`, (see :meth:`dim1_spline_knots`) and from the coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,q` in its B-spline representation
.. math::
s\left(x\right) = \sum_{{i = 1}}^qc_iN_i\left(x\right)\text{.}
Here :math:`q = \bar{n}+3`, where :math:`\bar{n}` is the number of intervals of the spline, and :math:`N_i\left(x\right)` denotes the normalized B-spline of degree :math:`3` defined upon the knots :math:`\lambda_i,\lambda_{{i+1}},\ldots,\lambda_{{i+4}}`.
The prescribed argument :math:`x` must satisfy :math:`\lambda_4\leq x\leq \lambda_{{\bar{n}+4}}`.
It is assumed that :math:`\lambda_{\textit{j}}\geq \lambda_{{\textit{j}-1}}`, for :math:`\textit{j} = 2,3,\ldots,\bar{n}+7`, and :math:`\lambda_{{\bar{n}+4}} > \lambda_4`.
If :math:`x` is a point at which :math:`4` knots coincide, :math:`s\left(x\right)` is discontinuous at :math:`x`; in this case, :math:`\mathrm{s}` contains the value defined as :math:`x` is approached from the right.
The method employed is that of evaluation by taking convex combinations due to de Boor (1972).
For further details of the algorithm and its use see Cox (1972) and Cox and Hayes (1973).
It is expected that a common use of ``dim1_spline_eval`` will be the evaluation of the cubic spline approximations produced by :meth:`dim1_spline_knots`.
A generalization of ``dim1_spline_eval`` which also forms the derivative of :math:`s\left(x\right)` is :meth:`dim1_spline_deriv`. :meth:`dim1_spline_deriv` takes about :math:`50\%` longer than ``dim1_spline_eval``.
.. _e02bb-py2-py-references:
**References**
Cox, M G, 1972, `The numerical evaluation of B-splines`, J. Inst. Math. Appl. (10), 134--149
Cox, M G, 1978, `The numerical evaluation of a spline from its B-spline representation`, J. Inst. Math. Appl. (21), 135--143
Cox, M G and Hayes, J G, 1973, `Curve fitting: a guide and suite of algorithms for the non-specialist user`, NPL Report NAC26, National Physical Laboratory
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
"""
raise NotImplementedError
[docs]def dim1_spline_deriv(lamda, c, x, left):
r"""
``dim1_spline_deriv`` evaluates a cubic spline and its first three derivatives from its B-spline representation.
.. _e02bc-py2-py-doc:
For full information please refer to the NAG Library document for e02bc
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bcf.html
.. _e02bc-py2-py-parameters:
**Parameters**
**lamda** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
:math:`\mathrm{lamda}[\textit{j}-1]` must be set to the value of the :math:`\textit{j}`\ th member of the complete set of knots, :math:`\lambda_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,\bar{n}+7`.
**c** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
The coefficient :math:`c_{\textit{i}}` of the B-spline :math:`N_{\textit{i}}\left(x\right)`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3`. The remaining elements of the array are not referenced.
**x** : float
The argument :math:`x` at which the cubic spline and its derivatives are to be evaluated.
**left** : int
Specifies whether left- or right-hand values of the spline and its derivatives are to be computed (see :ref:`Notes <e02bc-py2-py-notes>`). Left - or right-hand values are formed according to whether :math:`\mathrm{left}` is equal or not equal to :math:`1`.
If :math:`x` does not coincide with a knot, the value of :math:`\mathrm{left}` is immaterial.
If :math:`x = \mathrm{lamda}[3]`, right-hand values are computed.
If :math:`x = \mathrm{lamda}[\textit{ncap7}-4]`, left-hand values are formed, regardless of the value of :math:`\mathrm{left}`.
**Returns**
**s** : float, ndarray, shape :math:`\left(4\right)`
:math:`\mathrm{s}[\textit{j}-1]` contains the value of the :math:`\left(\textit{j}-1\right)`\ th derivative of the spline at the argument :math:`x`, for :math:`\textit{j} = 1,2,\ldots,4`. Note that :math:`\mathrm{s}[0]` contains the value of the spline.
.. _e02bc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncap7}\geq 8`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{ncap7}-4] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}\leq \mathrm{lamda}[\textit{ncap7}-4]`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[3] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}\geq \mathrm{lamda}[3]`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{lamda}[3] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{ncap7}-4] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[3] < \mathrm{lamda}[\textit{ncap7}-4]`.
.. _e02bc-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_spline_deriv`` evaluates the cubic spline :math:`s\left(x\right)` and its first three derivatives at a prescribed argument :math:`x`.
It is assumed that :math:`s\left(x\right)` is represented in terms of its B-spline coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3` and (augmented) ordered knot set :math:`\lambda_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+7`, (see :meth:`dim1_spline_knots`), i.e.,
.. math::
s\left(x\right) = \sum_{{i = 1}}^qc_iN_i\left(x\right)\text{.}
Here :math:`q = \bar{n}+3`, :math:`\bar{n}` is the number of intervals of the spline and :math:`N_i\left(x\right)` denotes the normalized B-spline of degree :math:`3` (order :math:`4`) defined upon the knots :math:`\lambda_i,\lambda_{{i+1}},\ldots,\lambda_{{i+4}}`.
The prescribed argument :math:`x` must satisfy
.. math::
\lambda_4\leq x\leq \lambda_{{\bar{n}+4}}\text{.}
At a simple knot :math:`\lambda_i` (i.e., one satisfying :math:`\lambda_{{i-1}} < \lambda_i < \lambda_{{i+1}}`), the third derivative of the spline is in general discontinuous.
At a multiple knot (i.e., two or more knots with the same value), lower derivatives, and even the spline itself, may be discontinuous.
Specifically, at a point :math:`x = u` where (exactly) :math:`r` knots coincide (such a point is termed a knot of multiplicity :math:`r`), the values of the derivatives of order :math:`4-\textit{j}`, for :math:`\textit{j} = 1,2,\ldots,r`, are in general discontinuous. (Here :math:`1\leq r\leq 4`; :math:`r > 4` is not meaningful.) You must specify whether the value at such a point is required to be the left- or right-hand derivative.
The method employed is based upon:
(i) carrying out a binary search for the knot interval containing the argument :math:`x` (see Cox (1978)),
(#) evaluating the nonzero B-splines of orders :math:`1`, :math:`2`, :math:`3` and :math:`4` by recurrence (see Cox (1972) and Cox (1978)),
(#) computing all derivatives of the B-splines of order :math:`4` by applying a second recurrence to these computed B-spline values (see de Boor (1972)),
(#) multiplying the fourth-order B-spline values and their derivative by the appropriate B-spline coefficients, and summing, to yield the values of :math:`s\left(x\right)` and its derivatives.
``dim1_spline_deriv`` can be used to compute the values and derivatives of cubic spline fits and interpolants produced by :meth:`dim1_spline_knots`.
If only values and not derivatives are required, :meth:`dim1_spline_eval` may be used instead of ``dim1_spline_deriv``, which takes about :math:`50\%` longer than :meth:`dim1_spline_eval`.
.. _e02bc-py2-py-references:
**References**
Cox, M G, 1972, `The numerical evaluation of B-splines`, J. Inst. Math. Appl. (10), 134--149
Cox, M G, 1978, `The numerical evaluation of a spline from its B-spline representation`, J. Inst. Math. Appl. (21), 135--143
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
"""
raise NotImplementedError
[docs]def dim1_spline_integ(lamda, c):
r"""
``dim1_spline_integ`` computes the definite integral of a cubic spline from its B-spline representation.
.. _e02bd-py2-py-doc:
For full information please refer to the NAG Library document for e02bd
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bdf.html
.. _e02bd-py2-py-parameters:
**Parameters**
**lamda** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
:math:`\mathrm{lamda}[\textit{j}-1]` must be set to the value of the :math:`\textit{j}`\ th member of the complete set of knots, :math:`\lambda_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,\bar{n}+7`.
**c** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
The coefficient :math:`c_{\textit{i}}` of the B-spline :math:`N_{\textit{i}}\left(x\right)`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3`. The remaining elements of the array are not referenced.
**Returns**
**dint** : float
The value of the definite integral of :math:`s\left(x\right)` between the limits :math:`x = a` and :math:`x = b`, where :math:`a = \lambda_4` and :math:`b = \lambda_{{\bar{n}+4}}`.
.. _e02bd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncap7}\geq 8`.
(`errno` :math:`2`)
On entry, :math:`J = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{lamda}[J-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[J-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[J-1]\geq \mathrm{lamda}[J-2]`.
(`errno` :math:`2`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{lamda}[\textit{ncap7}-2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{ncap7}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[\textit{ncap7}-2] = \mathrm{lamda}[\textit{ncap7}-1]`.
(`errno` :math:`2`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{lamda}[\textit{ncap7}-3] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{ncap7}-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[\textit{ncap7}-3] = \mathrm{lamda}[\textit{ncap7}-2]`.
(`errno` :math:`2`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{lamda}[\textit{ncap7}-4] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{ncap7}-4] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[\textit{ncap7}-4] = \mathrm{lamda}[\textit{ncap7}-3]`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{lamda}[2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[3] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[2] = \mathrm{lamda}[3]`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{lamda}[1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[1] = \mathrm{lamda}[2]`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{lamda}[0] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[0] = \mathrm{lamda}[1]`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{lamda}[3] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{NCAP4} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{NCAP4}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[3] < \mathrm{lamda}[\textit{NCAP4}-1]`.
.. _e02bd-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_spline_integ`` computes the definite integral of the cubic spline :math:`s\left(x\right)` between the limits :math:`x = a` and :math:`x = b`, where :math:`a` and :math:`b` are respectively the lower and upper limits of the range over which :math:`s\left(x\right)` is defined.
It is assumed that :math:`s\left(x\right)` is represented in terms of its B-spline coefficients :math:`c_i`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3` and (augmented) ordered knot set :math:`\lambda_i`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+7`, with :math:`\lambda_i = a`, for :math:`\textit{i} = 1,2,\ldots,4` and :math:`\lambda_i = b`, for :math:`\textit{i} = \bar{n}+4,\ldots,\bar{n}+7`, (see :meth:`dim1_spline_knots`), i.e.,
.. math::
s\left(x\right) = \sum_{{i = 1}}^qc_iN_i\left(x\right)\text{.}
Here :math:`q = \bar{n}+3`, :math:`\bar{n}` is the number of intervals of the spline and :math:`N_i\left(x\right)` denotes the normalized B-spline of degree :math:`3` (order :math:`4`) defined upon the knots :math:`\lambda_i,\lambda_{{i+1}},\ldots,\lambda_{{i+4}}`.
The method employed uses the formula given in Section 3 of Cox (1975).
``dim1_spline_integ`` can be used to determine the definite integrals of cubic spline fits and interpolants produced by :meth:`dim1_spline_knots`.
.. _e02bd-py2-py-references:
**References**
Cox, M G, 1975, `An algorithm for spline interpolation`, J. Inst. Math. Appl. (15), 95--108
"""
raise NotImplementedError
[docs]def dim1_spline_auto(start, x, y, w, s, n, lamda, comm):
r"""
``dim1_spline_auto`` computes a cubic spline approximation to an arbitrary set of data points.
The knots of the spline are located automatically, but a single argument must be specified to control the trade-off between closeness of fit and smoothness of fit.
.. _e02be-py2-py-doc:
For full information please refer to the NAG Library document for e02be
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bef.html
.. _e02be-py2-py-parameters:
**Parameters**
**start** : str, length 1
Must be set to 'C' or 'W'.
:math:`\mathrm{start} = \texttt{'C'}`
The function will build up the knot set starting with no interior knots. No values need be assigned to the arguments :math:`\mathrm{n}`, :math:`\mathrm{lamda}`, :math:`\mathrm{comm}`\ ['wrk'] or :math:`\mathrm{comm}`\ ['iwrk'].
:math:`\mathrm{start} = \texttt{'W'}`
The function will restart the knot-placing strategy using the knots found in a previous call of the function. In this case, the arguments :math:`\mathrm{n}`, :math:`\mathrm{lamda}`, :math:`\mathrm{comm}`\ ['wrk'], and :math:`\mathrm{comm}`\ ['iwrk'] must be unchanged from that previous call. This warm start can save much time in searching for a satisfactory value of :math:`\mathrm{s}`.
**x** : float, array-like, shape :math:`\left(m\right)`
The values :math:`x_{\textit{r}}` of the independent variable (abscissa) :math:`x`, for :math:`\textit{r} = 1,2,\ldots,m`.
**y** : float, array-like, shape :math:`\left(m\right)`
The values :math:`y_{\textit{r}}` of the dependent variable (ordinate) :math:`y`, for :math:`\textit{r} = 1,2,\ldots,m`.
**w** : float, array-like, shape :math:`\left(m\right)`
The values :math:`w_{\textit{r}}` of the weights, for :math:`\textit{r} = 1,2,\ldots,m`. For advice on the choice of weights, see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html#background12>`__.
**s** : float
The smoothing factor, :math:`S`.
If :math:`S = 0.0`, the function returns an interpolating spline.
If :math:`S` is smaller than machine precision, it is assumed equal to zero.
For advice on the choice of :math:`S`, see :ref:`Notes <e02be-py2-py-notes>` and `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bef.html#fcomments>`__.
**n** : int
If the warm start option is used, the value of :math:`\mathrm{n}` must be left unchanged from the previous call.
**lamda** : float, array-like, shape :math:`\left(\textit{nest}\right)`
If the warm start option is used, the values :math:`\mathrm{lamda}[0],\mathrm{lamda}[1],\ldots,\mathrm{lamda}[ \mathrm{n} -1]` must be left unchanged from the previous call.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**Returns**
**n** : int
The total number, :math:`n`, of knots of the computed spline.
**lamda** : float, ndarray, shape :math:`\left(\textit{nest}\right)`
The knots of the spline, i.e., the positions of the interior knots :math:`\mathrm{lamda}[4], \mathrm{lamda}[5],\ldots,\mathrm{lamda}[ \mathrm{n}-5]` as well as the positions of the additional knots
.. math::
\mathrm{lamda}[0] = \mathrm{lamda}[1] = \mathrm{lamda}[2] = \mathrm{lamda}[3] = x_1
and
.. math::
\mathrm{lamda}[ \mathrm{n}-4] = \mathrm{lamda}[ \mathrm{n}-3] = \mathrm{lamda}[ \mathrm{n}-2] = \mathrm{lamda}[ \mathrm{n} -1] = x_m
needed for the B-spline representation.
**c** : float, ndarray, shape :math:`\left(\textit{nest}\right)`
The coefficient :math:`c_{\textit{i}}` of the B-spline :math:`N_{\textit{i}}\left(x\right)` in the spline approximation :math:`s\left(x\right)`, for :math:`\textit{i} = 1,2,\ldots,n-4`.
**fp** : float
The sum of the squared weighted residuals, :math:`\theta`, of the computed spline approximation. If :math:`\mathrm{fp} = 0.0`, this is an interpolating spline. :math:`\mathrm{fp}` should equal :math:`\mathrm{s}` within a relative tolerance of :math:`0.001` unless :math:`n = 8` when the spline has no interior knots and so is simply a cubic polynomial. For knots to be inserted, :math:`\mathrm{s}` must be set to a value below the value of :math:`\mathrm{fp}` produced in this case.
.. _e02be-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`r = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{w}[r-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{w}[r-1] > 0.0` for all :math:`r`.
(`errno` :math:`1`)
On entry, :math:`\textit{nest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nest}\geq 8`.
(`errno` :math:`1`)
On entry, :math:`\textit{nest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nest}\geq \langle\mathit{\boldsymbol{value}}\rangle` when :math:`\mathrm{s} = 0.0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{s}\geq 0.0`.
(`errno` :math:`1`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 4`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{start} \neq \texttt{'W'}` or :math:`\texttt{'C'}`: :math:`\mathrm{start} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`2`)
On entry, :math:`r = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{w}[r-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{w}[r-1] > 0.0` for all :math:`r`.
(`errno` :math:`3`)
On entry, :math:`r = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[r-2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[r-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[r-2] < \mathrm{x}[r-1]` for all :math:`r`.
(`errno` :math:`4`)
The number of knots needed is greater than :math:`\textit{nest}`: :math:`\textit{nest} = \langle\mathit{\boldsymbol{value}}\rangle`. Possibly :math:`\mathrm{s}` is too small: :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
The iterative process has failed to converge. Possibly :math:`\mathrm{s}` is too small: :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _e02be-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim1_spline_auto`` determines a smooth cubic spline approximation :math:`s\left(x\right)` to the set of data points :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)`, with weights :math:`w_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`.
The spline is given in the B-spline representation
.. math::
s\left(x\right) = \sum_{{i = 1}}^{{n-4}}c_iN_i\left(x\right)\text{,}
where :math:`N_i\left(x\right)` denotes the normalized cubic B-spline defined upon the knots :math:`\lambda_i,\lambda_{{i+1}},\ldots,\lambda_{{i+4}}`.
The total number :math:`n` of these knots and their values :math:`\lambda_1,\ldots,\lambda_n` are chosen automatically by the function.
The knots :math:`\lambda_5,\ldots,\lambda_{{n-4}}` are the interior knots; they divide the approximation interval :math:`\left[x_1, x_m\right]` into :math:`n-7` sub-intervals.
The coefficients :math:`c_1,c_2,\ldots,c_{{n-4}}` are then determined as the solution of the following constrained minimization problem:
minimize
.. math::
\eta = \sum_{{i = 5}}^{{n-4}}\delta_i^2
subject to the constraint
.. math::
\theta = \sum_{{r = 1}}^m\epsilon_r^2\leq S\text{,}
.. rst-class:: nag-rules-none nag-align-left
+-----+------------------+----------------------------------------------------------------------------------------------------------------------------------+
|where|:math:`\delta_i` |stands for the discontinuity jump in the third order derivative of :math:`s\left(x\right)` at the interior knot :math:`\lambda_i`,|
+-----+------------------+----------------------------------------------------------------------------------------------------------------------------------+
| |:math:`\epsilon_r`|denotes the weighted residual :math:`w_r\left(y_r-s\left(x_r\right)\right)`, |
+-----+------------------+----------------------------------------------------------------------------------------------------------------------------------+
|and |:math:`S` |is a non-negative number to be specified by you. |
+-----+------------------+----------------------------------------------------------------------------------------------------------------------------------+
The quantity :math:`\eta` can be seen as a measure of the (lack of) smoothness of :math:`s\left(x\right)`, while closeness of fit is measured through :math:`\theta`.
By means of the argument :math:`S`, 'the smoothing factor', you can then control the balance between these two (usually conflicting) properties.
If :math:`S` is too large, the spline will be too smooth and signal will be lost (underfit); if :math:`S` is too small, the spline will pick up too much noise (overfit).
In the extreme cases the function will return an interpolating spline :math:`\left(\theta = 0\right)` if :math:`S` is set to zero, and the weighted least squares cubic polynomial :math:`\left(\eta = 0\right)` if :math:`S` is set very large.
Experimenting with :math:`S` values between these two extremes should result in a good compromise. (See `Choice of S <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bef.html#fc-choice>`__ for advice on choice of :math:`S`.)
The method employed is outlined in `Outline of Method Used <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bef.html#fc-method>`__ and fully described in Dierckx (1975), Dierckx (1981) and Dierckx (1982).
It involves an adaptive strategy for locating the knots of the cubic spline (depending on the function underlying the data and on the value of :math:`S`), and an iterative method for solving the constrained minimization problem once the knots have been determined.
Values of the computed spline, or of its derivatives or definite integral, can subsequently be computed by calling :meth:`dim1_spline_eval`, :meth:`dim1_spline_deriv` or :meth:`dim1_spline_integ`, as described in `Evaluation of Computed Spline <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bef.html#fc-evaluation>`__.
.. _e02be-py2-py-references:
**References**
Dierckx, P, 1975, `An algorithm for smoothing, differentiating and integration of experimental data using spline functions`, J. Comput. Appl. Math. (1), 165--184
Dierckx, P, 1981, `An improved algorithm for curve fitting with spline functions`, Report TW54, Department of Computer Science, Katholieke Univerciteit Leuven
Dierckx, P, 1982, `A fast algorithm for smoothing data on a rectangular grid while using spline functions`, SIAM J. Numer. Anal. (19), 1286--1304
Reinsch, C H, 1967, `Smoothing by spline functions`, Numer. Math. (10), 177--183
"""
raise NotImplementedError
[docs]def dim1_spline_deriv_vector(start, lamda, c, deriv, xord, x, ixloc, iwrk=None):
r"""
``dim1_spline_deriv_vector`` evaluates a cubic spline and up to its first three derivatives from its B-spline representation at a vector of points. ``dim1_spline_deriv_vector`` can be used to compute the values and derivatives of cubic spline fits and interpolants produced by reference to :meth:`interp.dim1_spline <naginterfaces.library.interp.dim1_spline>`, :meth:`dim1_spline_knots` and :meth:`dim1_spline_auto`.
.. _e02bf-py2-py-doc:
For full information please refer to the NAG Library document for e02bf
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bff.html
.. _e02bf-py2-py-parameters:
**Parameters**
**start** : int
Indicates the completion state of the first phase of the algorithm.
:math:`\mathrm{start} = 0`
The enclosing interval numbers :math:`\textit{ixloc}_j` for the abscissae :math:`x_j` contained in :math:`\mathrm{x}` have not been determined, and you wish to use the sorted mode of vectorization.
:math:`\mathrm{start} = 1`
The enclosing interval numbers :math:`\textit{ixloc}_j` have been determined and are provided in :math:`\mathrm{ixloc}`, however the required permutation and interval related information has not been determined and you wish to use the sorted mode of vectorization.
:math:`\mathrm{start} = 2`
You wish to use the sorted mode of vectorization, and the entire first phase has been completed, with the enclosing interval numbers supplied in :math:`\mathrm{ixloc}`, and the required permutation and interval related information provided in :math:`\mathrm{iwrk}` (from a previous call to ``dim1_spline_deriv_vector``).
:math:`\mathrm{start} = 10`
The enclosing interval numbers :math:`\textit{ixloc}_j` for the abscissae :math:`x_j` contained in :math:`\mathrm{x}` have not been determined, and you wish to use the unsorted mode of vectorization.
:math:`\mathrm{start} = 11`
The enclosing interval numbers :math:`\textit{ixloc}_j` for the abscissae :math:`x_j` contained in :math:`\mathrm{x}` have been supplied in :math:`\mathrm{ixloc}`, and you wish to use the unsorted mode of vectorization.
`Additional`: :math:`\mathrm{start} = 0` or :math:`10` should be used unless you are sure that the knot set is unchanged between calls.
**lamda** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
:math:`\mathrm{lamda}[\textit{j}-1]` must be set to the value of the :math:`\textit{j}`\ th member of the complete set of knots, :math:`\lambda_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,\bar{n}+7`.
**c** : float, array-like, shape :math:`\left(\textit{ncap7}\right)`
The coefficient :math:`c_{\textit{i}}` of the B-spline :math:`N_{\textit{i}}\left(x\right)`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3`. The remaining elements of the array are not referenced.
**deriv** : int
Determines the maximum order of derivatives required, :math:`D_{\mathrm{ord}}`.
If :math:`\mathrm{deriv} < 0` left derivatives are calculated, otherwise right derivatives are calculated. For abscissae satisfying :math:`x_j = \lambda_4` or :math:`x_j = \lambda_{{\bar{n}+4}}` only right-handed or left-handed computation will be used respectively.
For abscissae which do not coincide exactly with a knot, the handedness of the computation is immaterial.
:math:`\mathrm{deriv} = 0`
No derivatives required.
:math:`\mathrm{deriv} = \pm 1`
Only :math:`s\left(x\right)` and its first derivative are required.
:math:`\mathrm{deriv} = \pm 2`
Only :math:`s\left(x\right)` and its first and second derivatives are required.
:math:`\mathrm{deriv} = \pm 3`
:math:`s\left(x\right)` and its first, second and third derivatives are required.
Note: if :math:`\left\lvert \mathrm{deriv}\right\rvert` is greater than :math:`3` only the derivatives up to and including :math:`3` will be returned.
**xord** : int
Indicates whether :math:`\mathrm{x}` is supplied in a sufficiently ordered manner. If :math:`\mathrm{x}` is sufficiently ordered ``dim1_spline_deriv_vector`` will complete faster.
:math:`\mathrm{xord} \neq 0`
The abscissae in :math:`\mathrm{x}` are ordered at least by ascending interval, in that any two abscissae contained in the same interval are only separated by abscissae in the same interval, and the intervals are arranged in ascending order. For example, :math:`x_j < x_{{\textit{j}+1}}`, for :math:`\textit{j} = 1,2,\ldots,\textit{nx}-1`.
:math:`\mathrm{xord} = 0`
The abscissae in :math:`\mathrm{x}` are not sufficiently ordered.
**x** : float, array-like, shape :math:`\left(\textit{nx}\right)`
The abscissae :math:`x_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,n_x`. If :math:`\mathrm{start} = 0` or :math:`10` then evaluations will only be performed for these :math:`x_j` satisfying :math:`\lambda_4\leq x_j\leq \lambda_{{\bar{n}+4}}`. Otherwise evaluation will be performed unless the corresponding element of :math:`\mathrm{ixloc}` contains an invalid interval number. Please note that if the :math:`\mathrm{ixloc}[j]` is a valid interval number then no check is made that :math:`\mathrm{x}[j]` actually lies in that interval.
**ixloc** : int, array-like, shape :math:`\left(\textit{nx}\right)`
If :math:`\mathrm{start} = 1`, :math:`2` or :math:`11`, if you wish :math:`x_j` to be evaluated, :math:`\mathrm{ixloc}[j-1]` must be the enclosing interval number :math:`\textit{ixloc}_j` of the abscissae :math:`x_j` (see `[equation] <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02bff.html#ixloc_def>`__). If you do not wish :math:`x_j` to be evaluated, you may set the interval number to be either less than :math:`4` or greater than :math:`\bar{n}+4`.
Otherwise, :math:`\mathrm{ixloc}` need not be set.
**iwrk** : None or int, array-like, shape :math:`\left(\textit{liwrk}\right)`, optional
If :math:`\mathrm{start} = 2`, :math:`\mathrm{iwrk}` must be unchanged from a previous call to ``dim1_spline_deriv_vector`` with :math:`\mathrm{start} = 0` or :math:`1`.
Otherwise, :math:`\mathrm{iwrk}` need not be set.
**Returns**
**ixloc** : int, ndarray, shape :math:`\left(\textit{nx}\right)`
If :math:`\mathrm{start} = 1`, :math:`2` or :math:`11`, :math:`\mathrm{ixloc}` is unchanged on exit.
Otherwise, :math:`\mathrm{ixloc}[\textit{j}-1]`, contains the enclosing interval number :math:`\textit{ixloc}_{\textit{j}}`, for the abscissa supplied in :math:`\mathrm{x}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,n_x`.
Evaluations will only be performed for abscissae :math:`x_j` satisfying :math:`\lambda_4\leq x_j\leq \lambda_{{\bar{n}+4}}`.
If evaluation is not performed :math:`\mathrm{ixloc}[j-1]` is set to :math:`0` if :math:`x_j < \lambda_4` or :math:`\bar{n}+7` if :math:`x_j > \lambda_{{\bar{n}+4}}`.
**s** : float, ndarray, shape :math:`\left(\textit{nx}, \min\left(\left\lvert \mathrm{deriv}\right\rvert,3\right)+1\right)`
If :math:`x_j` is valid, :math:`\mathrm{s}[\textit{j}-1,\textit{d}-1]` will contain the (:math:`\textit{d}-1`)th derivative of :math:`s\left(x\right)`, for :math:`\textit{j} = 1,2,\ldots,n_x`, for :math:`\textit{d} = 1,2,\ldots,D_{\mathrm{ord}}+1`. In particular, :math:`\mathrm{s}[j-1,0]` will contain the approximation of :math:`s\left(x_j\right)` for all legal values in :math:`\mathrm{x}`.
**iwrk** : None or int, ndarray, shape :math:`\left(\textit{liwrk}\right)`
If :math:`\mathrm{start} = 10` or :math:`11`, :math:`\mathrm{iwrk}` is unchanged on exit.
Otherwise, :math:`\mathrm{iwrk}` contains the required permutation of elements of :math:`\mathrm{x}`, if any, and information related to the division of the abscissae :math:`x_j` between the intervals derived from :math:`\mathrm{lamda}`.
.. _e02bf-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`2`)
On entry, all elements of :math:`\mathrm{x}` had enclosing interval numbers in :math:`\mathrm{ixloc}` outside the domain allowed by the provided spline.
:math:`\langle\mathit{\boldsymbol{value}}\rangle` entries of :math:`\mathrm{x}` were indexed below the lower bound :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
:math:`\langle\mathit{\boldsymbol{value}}\rangle` entries of :math:`\mathrm{x}` were indexed above the upper bound :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`11`)
On entry, :math:`\mathrm{start} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{start} = 0`, :math:`1`, :math:`2`, :math:`10` or :math:`11`.
(`errno` :math:`12`)
On entry, :math:`\mathrm{start} = 2` and :math:`\textit{nx}` is not consistent with the previous call to ``dim1_spline_deriv_vector``.
On entry, :math:`\textit{nx} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nx} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`21`)
On entry, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ncap7}\geq 8`.
(`errno` :math:`31`)
On entry, :math:`\mathrm{lamda}[3] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{ncap7} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{ncap7}-4] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[3] < \mathrm{lamda}[\textit{ncap7}-4]`.
(`errno` :math:`91`)
On entry, :math:`\textit{nx} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nx}\geq 1`.
(`errno` :math:`131`)
On entry, :math:`\textit{liwrk} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{liwrk}\geq 3\times \textit{nx}+3 = \langle\mathit{\boldsymbol{value}}\rangle`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
On entry, at least one element of :math:`\mathrm{x}` has an enclosing interval number in :math:`\mathrm{ixloc}` outside the set allowed by the provided spline. The spline has been evaluated for all :math:`\mathrm{x}` with enclosing interval numbers inside the allowable set.
:math:`\langle\mathit{\boldsymbol{value}}\rangle` entries of :math:`\mathrm{x}` were indexed below the lower bound :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
:math:`\langle\mathit{\boldsymbol{value}}\rangle` entries of :math:`\mathrm{x}` were indexed above the upper bound :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _e02bf-py2-py-notes:
**Notes**
``dim1_spline_deriv_vector`` evaluates the cubic spline :math:`s\left(x\right)` and optionally derivatives up to order :math:`3` for a vector of points :math:`x_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,n_x`.
It is assumed that :math:`s\left(x\right)` is represented in terms of its B-spline coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+3`, and (augmented) ordered knot set :math:`\lambda_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\bar{n}+7`, (see :meth:`dim1_spline_knots` and :meth:`dim1_spline_auto`), i.e.,
.. math::
s\left(x\right) = \sum_{{i = 1}}^qc_iN_i\left(x\right)\text{.}
Here :math:`q = \bar{n}+3`, :math:`\bar{n}` is the number of intervals of the spline and :math:`N_i\left(x\right)` denotes the normalized B-spline of degree :math:`3` (order :math:`4`) defined upon the knots :math:`\lambda_i,\lambda_{{i+1}},\ldots,\lambda_{{i+4}}`.
The knots :math:`\lambda_5,\lambda_6,\ldots,\lambda_{{\bar{n}+3}}` are the interior knots.
The remaining knots, :math:`\lambda_1`, :math:`\lambda_2`, :math:`\lambda_3`, :math:`\lambda_4` and :math:`\lambda_{{\bar{n}+4}}`, :math:`\lambda_{{\bar{n}+5}}`, :math:`\lambda_{{\bar{n}+6}}`, :math:`\lambda_{{\bar{n+7}}}` are the exterior knots.
The knots :math:`\lambda_4` and :math:`\lambda_{{\bar{n}+4}}` are the boundaries of the spline.
Only abscissae satisfying,
.. math::
\lambda_4\leq x_j\leq \lambda_{{\bar{n}+4}}\text{,}
will be evaluated.
At a simple knot :math:`\lambda_i` (i.e., one satisfying :math:`\lambda_{{i-1}} < \lambda_i < \lambda_{{i+1}}`), the third derivative of the spline is, in general, discontinuous.
At a multiple knot (i.e., two or more knots with the same value), lower derivatives, and even the spline itself, may be discontinuous.
Specifically, at a point :math:`x = u` where (exactly) :math:`r` knots coincide (such a point is termed a knot of multiplicity :math:`r`), the values of the derivatives of order :math:`4-\textit{j}`, for :math:`\textit{j} = 1,2,\ldots,r`, are, in general, discontinuous. (Here :math:`1\leq r\leq 4`; :math:`r > 4` is not meaningful.) The maximum order of the derivatives to be evaluated :math:`D_{\mathrm{ord}}`, and the left- or right-handedness of the computation when an abscissa corresponds exactly to an interior knot, are determined by the value of :math:`\mathrm{deriv}`.
Each abscissa (point at which the spline is to be evaluated) :math:`x_j` contained in :math:`\mathrm{x}` has an associated enclosing interval number, :math:`\textit{ixloc}_j` either supplied or returned in :math:`\mathrm{ixloc}` (see argument :math:`\mathrm{start}`).
A simple call to ``dim1_spline_deriv_vector`` would set :math:`\mathrm{start} = 0` and the contents of :math:`\mathrm{ixloc}` need never be set nor referenced, and the following description on modes of operation can be ignored.
However, where efficiency is an important consideration, the following description will help to choose the appropriate mode of operation.
The interval numbers are used to determine which B-splines must be evaluated for a given abscissa, and are defined as
.. math::
\textit{ixloc}_j = \begin{pmatrix}\leq 0& x_j < \lambda_1 \\4& \lambda_4 = x_j \\k& \lambda_k < x_j < \lambda_{{k+1}} &\\k& \lambda_4 < \lambda_k = x_j &\text{left derivatives}\\k& x_j = \lambda_{{k+1}} < \lambda_{{\bar{n}+4}} &\text{right derivatives or no derivatives}\\\bar{n}+4& \lambda_{{\bar{n}+4}} = x_j \\ > \bar{n}+7& x_j > \lambda_{{\bar{n}+7}} \end{pmatrix}
The algorithm has two modes of vectorization, termed here sorted and unsorted, which are selectable by the argument :math:`\mathrm{start}`.
Furthermore, if the supplied abscissae are sufficiently ordered, as indicated by the argument :math:`\mathrm{xord}`, the algorithm will take advantage of significantly faster methods for the determination of both the interval numbers and the subsequent spline evaluations.
The sorted mode has two phases, a sorting phase and an evaluation phase.
This mode is recommended if there are many abscissae to evaluate relative to the number of intervals of the spline, or the abscissae are distributed relatively densely over a subsection of the spline.
In the first phase, :math:`\textit{ixloc}_j` is determined for each :math:`x_j` and a permutation is calculated to sort the :math:`x_j` by interval number.
The first phase may be either partially or completely by-passed using the argument :math:`\mathrm{start}` if the enclosing segments and/or the subsequent ordering are already known `a priori`, for example if multiple spline coefficients :math:`\mathrm{c}` are to be evaluated over the same set of knots :math:`\mathrm{lamda}`.
In the second phase of the sorted mode, spline approximations are evaluated by segment, so that non-abscissa dependent calculations over a segment may be reused in the evaluation for all abscissae belonging to a specific segment.
For example, all third derivatives of all abscissae in the same segment will be identical.
In the unsorted mode of vectorization, no `a priori` segment sorting is performed, and if the abscissae are not sufficiently ordered, the evaluation at an abscissa will be independent of evaluations at other abscissae; also non-abscissa dependent calculations over a segment will be repeated for each abscissa in a segment.
This may be quicker if the number of abscissa is small in comparison to the number of knots in the spline, and they are distributed sparsely throughout the domain of the spline.
This is effectively a direct vectorization of :meth:`dim1_spline_eval` and :meth:`dim1_spline_deriv`, although if the enclosing interval numbers :math:`\textit{ixloc}_j` are known, these may again be provided.
If the abscissae are sufficiently ordered, then once the first abscissa in a segment is known, an efficient algorithm will be used to determine the location of the final abscissa in this segment.
The spline will subsequently be evaluated in a vectorized manner for all the abscissae indexed between the first and last of the current segment.
If no derivatives are required, the spline evaluation is calculated by taking convex combinations due to de Boor (1972).
Otherwise, the calculation of :math:`s\left(x\right)` and its derivatives is based upon,
(i) evaluating the nonzero B-splines of orders :math:`1`, :math:`2`, :math:`3` and :math:`4` by recurrence (see Cox (1972) and Cox (1978)),
(#) computing all derivatives of the B-splines of order :math:`4` by applying a second recurrence to these computed B-spline values (see de Boor (1972)),
(#) multiplying the fourth-order B-spline values and their derivative by the appropriate B-spline coefficients, and summing, to yield the values of :math:`s\left(x\right)` and its derivatives.
The method of convex combinations is significantly faster than the recurrence based method.
If higher derivatives of order :math:`2` or :math:`3` are not required, as much computation as possible is avoided.
.. _e02bf-py2-py-references:
**References**
Cox, M G, 1972, `The numerical evaluation of B-splines`, J. Inst. Math. Appl. (10), 134--149
Cox, M G, 1978, `The numerical evaluation of a spline from its B-spline representation`, J. Inst. Math. Appl. (21), 135--143
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
"""
raise NotImplementedError
[docs]def dim2_cheb_lines(m, k, l, x, y, f, w, na, xmin, xmax, nux, nuy):
r"""
``dim2_cheb_lines`` forms an approximation to the weighted, least squares Chebyshev series surface fit to data arbitrarily distributed on lines parallel to one independent coordinate axis.
.. _e02ca-py2-py-doc:
For full information please refer to the NAG Library document for e02ca
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02caf.html
.. _e02ca-py2-py-parameters:
**Parameters**
**m** : int, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{m}[\textit{s}-1]` must be set to :math:`m_{\textit{s}}`, the number of data :math:`x` values on the line :math:`y = y_{\textit{s}}`, for :math:`\textit{s} = 1,2,\ldots,n`.
**k** : int
:math:`k`, the required degree of :math:`x` in the fit.
**l** : int
:math:`l`, the required degree of :math:`y` in the fit.
**x** : float, array-like, shape :math:`\left(\mathrm{sum}\left(\mathrm{m}\right)\right)`
The :math:`x` values of the data points. The sequence must be
all points on :math:`y = y_1`, followed by
all points on :math:`y = y_2`, followed by
:math:`\vdots`
all points on :math:`y = y_n`.
**y** : float, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{y}[\textit{s}-1]` must contain the :math:`y` value of line :math:`y = y_{\textit{s}}`, for :math:`\textit{s} = 1,2,\ldots,n`, on which data is given.
**f** : float, array-like, shape :math:`\left(\mathrm{sum}\left(\mathrm{m}\right)\right)`
:math:`f`, the data values of the dependent variable in the same sequence as the :math:`x` values.
**w** : float, array-like, shape :math:`\left(\mathrm{sum}\left(\mathrm{m}\right)\right)`
The weights to be assigned to the data points, in the same sequence as the :math:`x` values. These weights should be calculated from estimates of the absolute accuracies of the :math:`f_r`, expressed as standard deviations, probable errors or some other measure which is of the same dimensions as :math:`f_r`. Specifically, each :math:`w_r` should be inversely proportional to the accuracy estimate of :math:`f_r`. Often weights all equal to unity will be satisfactory. If a particular weight is zero, the corresponding data point is omitted from the fit.
**na** : int
The dimension of the array :math:`\mathrm{a}`.
**xmin** : float, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{xmin}[\textit{s}-1]` must contain :math:`x_{\mathrm{min}}^{\left(\textit{s}\right)}`, the lower end of the range of :math:`x` on the line :math:`y = y_{\textit{s}}`, for :math:`\textit{s} = 1,2,\ldots,n`. It must not be greater than the lowest data value of :math:`x` on the line. Each :math:`x_{\mathrm{min}}^{\left(s\right)}` is scaled to :math:`{-1.0}` in the fit. (See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02caf.html#fcomments>`__.)
**xmax** : float, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{xmax}[\textit{s}-1]` must contain :math:`x_{\mathrm{max}}^{\left(\textit{s}\right)}`, the upper end of the range of :math:`x` on the line :math:`y = y_{\textit{s}}`, for :math:`\textit{s} = 1,2,\ldots,n`. It must not be less than the highest data value of :math:`x` on the line. Each :math:`x_{\mathrm{max}}^{\left(s\right)}` is scaled to :math:`{+1.0}` in the fit. (See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02caf.html#fcomments>`__.)
**nux** : float, array-like, shape :math:`\left(\textit{inuxp1}\right)`
:math:`\mathrm{nux}[\textit{i}-1]` must contain the coefficient of the Chebyshev polynomial of degree :math:`\left(\textit{i}-1\right)` in :math:`\bar{x}`, in the Chebyshev series representation of the polynomial factor in :math:`\bar{x}` which you require the fit to contain, for :math:`\textit{i} = 1,2,\ldots,\textit{inuxp1}`. These coefficients are defined according to the standard convention of :ref:`Notes <e02ca-py2-py-notes>`.
**nuy** : float, array-like, shape :math:`\left(\textit{inuyp1}\right)`
:math:`\mathrm{nuy}[\textit{i}-1]` must contain the coefficient of the Chebyshev polynomial of degree :math:`\left(\textit{i}-1\right)` in :math:`\bar{y}`, in the Chebyshev series representation of the polynomial factor which you require the fit to contain, for :math:`\textit{i} = 1,2,\ldots,\textit{inuyp1}`. These coefficients are defined according to the standard convention of :ref:`Notes <e02ca-py2-py-notes>`.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{na}\right)`
Contains the Chebyshev coefficients of the fit. :math:`\mathrm{a}[i\times \left(\mathrm{l}+1\right)+j-1]` is the coefficient :math:`a_{{ij}}` of :ref:`Notes <e02ca-py2-py-notes>` defined according to the standard convention. These coefficients are used by :meth:`dim2_cheb_eval` to calculate values of the fitted function.
.. _e02ca-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{na}` is too small. :math:`\mathrm{na} = \langle\mathit{\boldsymbol{value}}\rangle`. Minimum possible dimension: :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\textit{mtot}` is too small. :math:`\textit{mtot} = \langle\mathit{\boldsymbol{value}}\rangle`. Minimum possible dimension: :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{m}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{inuxp1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}[\textit{I}-1]\geq \mathrm{k}-\textit{inuxp1}+2`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{l} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{inuyp1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq \mathrm{l}-\textit{inuyp1}+2`.
(`errno` :math:`1`)
On entry, :math:`\textit{inuyp1} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{l} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{inuyp1}\leq \mathrm{l}+1`.
(`errno` :math:`1`)
On entry, :math:`\textit{inuxp1} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{inuxp1}\leq \mathrm{k}+1`.
(`errno` :math:`1`)
On entry, :math:`\textit{inuyp1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{inuyp1}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{inuxp1} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{inuxp1}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{l} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{l}\geq 0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{k}\geq 0`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{xmin}[\textit{I}-1]` and :math:`\mathrm{xmax}[\textit{I}-1]` do not span the data :math:`\mathrm{x}` values on :math:`\mathrm{y} = \mathrm{y}[\textit{I}-1]`: :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{xmin}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{xmax}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{y}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`3`)
On entry, the data :math:`\mathrm{x}` values are not nondecreasing for :math:`\mathrm{y} = \mathrm{y}[\textit{I}-1]`: :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{y}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`3`)
On entry, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{y}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{y}[\textit{I}-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{y}[\textit{I}-1] > \mathrm{y}[\textit{I}-2]`.
(`errno` :math:`4`)
On entry, the number of distinct :math:`\mathrm{x}` values with nonzero weight on :math:`\mathrm{y} = \mathrm{y}[\textit{I}-1]` is less than :math:`\mathrm{k}-\textit{inuxp1}+2`: :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{y}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{inuxp1} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
On entry, :math:`\textit{inuxp1} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{nux}[\textit{inuxp1}-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{inuyp1} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nuy}[\textit{inuyp1}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{nux}[\textit{inuxp1}-1] = 0.0`, :math:`\textit{inuxp1} = 1`; if :math:`\mathrm{nuy}[\textit{inuyp1}-1] = 0.0`, :math:`\textit{inuyp1} = 1`.
.. _e02ca-py2-py-notes:
**Notes**
``dim2_cheb_lines`` determines a bivariate polynomial approximation of degree :math:`k` in :math:`x` and :math:`l` in :math:`y` to the set of data points :math:`\left(x_{{\textit{r},\textit{s}}}, y_{\textit{s}}, f_{{\textit{r},\textit{s}}}\right)`, with weights :math:`w_{{\textit{r},\textit{s}}}`, for :math:`\textit{r} = 1,2,\ldots,m_{\textit{s}}`, for :math:`\textit{s} = 1,2,\ldots,n`.
That is, the data points are on lines :math:`y = y_s`, but the :math:`x` values may be different on each line.
The values of :math:`k` and :math:`l` are prescribed by you (for guidance on their choice, see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02caf.html#fcomments>`__).
The function is based on the method described in Sections 5 and 6 of Clenshaw and Hayes (1965).
The polynomial is represented in double Chebyshev series form with arguments :math:`\bar{x}` and :math:`\bar{y}`.
The arguments lie in the range :math:`-1` to :math:`+1` and are related to the original variables :math:`x` and :math:`y` by the transformations
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)}\quad \text{ and }\quad \bar{y} = \frac{{2y-\left(y_{\mathrm{max}}+y_{\mathrm{min}}\right)}}{\left(y_{\mathrm{max}}-y_{\mathrm{min}}\right)}\text{.}
Here :math:`y_{\mathrm{max}}` and :math:`y_{\mathrm{min}}` are set by the function to, respectively, the largest and smallest value of :math:`y_s`, but :math:`x_{\mathrm{max}}` and :math:`x_{\mathrm{min}}` are functions of :math:`y` prescribed by you (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02caf.html#fcomments>`__).
For this function, only their values :math:`x_{\mathrm{max}}^{\left(s\right)}` and :math:`x_{\mathrm{min}}^{\left(s\right)}` at each :math:`y = y_s` are required.
For each :math:`s = 1,2,\ldots,n`, :math:`x_{\mathrm{max}}^{\left(s\right)}` must not be less than the largest :math:`x_{{r,s}}` on the line :math:`y = y_s`, and, similarly, :math:`x_{\mathrm{min}}^{\left(s\right)}` must not be greater than the smallest :math:`x_{{r,s}}`.
The double Chebyshev series can be written as
.. math::
\sum_{{i = 0}}^k\sum_{{j = 0}}^la_{{ij}}T_i\left(\bar{x}\right)T_j\left(\bar{y}\right)
where :math:`T_i\left(\bar{x}\right)` is the Chebyshev polynomial of the first kind of degree :math:`i` with argument :math:`\bar{x}`, and :math:`T_j\left(y\right)` is similarly defined.
However, the standard convention, followed in this function, is that coefficients in the above expression which have either :math:`i` or :math:`j` zero are written as :math:`\frac{1}{2}a_{{ij}}`, instead of simply :math:`a_{{ij}}`, and the coefficient with both :math:`i` and :math:`j` equal to zero is written as :math:`\frac{1}{4}a_{{0,0}}`.
The series with coefficients output by the function should be summed using this convention. :meth:`dim2_cheb_eval` is available to compute values of the fitted function from these coefficients.
The function first obtains Chebyshev series coefficients :math:`c_{{s,\textit{i}}}`, for :math:`\textit{i} = 0,1,\ldots,k`, of the weighted least squares polynomial curve fit of degree :math:`k` in :math:`\bar{x}` to the data on each line :math:`y = y_{\textit{s}}`, for :math:`\textit{s} = 1,2,\ldots,n`, in turn, using an auxiliary function.
The same function is then called :math:`k+1` times to fit :math:`c_{{\textit{s},i}}`, for :math:`\textit{s} = 1,2,\ldots,n`, by a polynomial of degree :math:`l` in :math:`\bar{y}`, for each :math:`i = 0,1,\ldots,k`.
The resulting coefficients are the required :math:`a_{{ij}}`.
You can force the fit to contain a given polynomial factor.
This allows for the surface fit to be constrained to have specified values and derivatives along the boundaries :math:`x = x_{\mathrm{min}}`, :math:`x = x_{\mathrm{max}}`, :math:`y = y_{\mathrm{min}}` and :math:`y = y_{\mathrm{max}}` or indeed along any lines :math:`\bar{x} = \text{}` constant or :math:`\bar{y} = \text{}` constant (see Section 8 of Clenshaw and Hayes (1965)).
.. _e02ca-py2-py-references:
**References**
Clenshaw, C W and Hayes, J G, 1965, `Curve and surface fitting`, J. Inst. Math. Appl. (1), 164--183
Hayes, J G (ed.), 1970, `Numerical Approximation to Functions and Data`, Athlone Press, London
"""
raise NotImplementedError
[docs]def dim2_cheb_eval(mfirst, k, l, x, xmin, xmax, y, ymin, ymax, a):
r"""
``dim2_cheb_eval`` evaluates a bivariate polynomial from the rectangular array of coefficients in its double Chebyshev series representation.
.. _e02cb-py2-py-doc:
For full information please refer to the NAG Library document for e02cb
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02cbf.html
.. _e02cb-py2-py-parameters:
**Parameters**
**mfirst** : int
The index of the first and last :math:`x` value in the array :math:`x` at which the evaluation is required respectively (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02cbf.html#fcomments>`__).
**k** : int
The degree :math:`k` of :math:`x` and :math:`l` of :math:`y`, respectively, in the polynomial.
**l** : int
The degree :math:`k` of :math:`x` and :math:`l` of :math:`y`, respectively, in the polynomial.
**x** : float, array-like, shape :math:`\left(\textit{mlast}\right)`
:math:`\mathrm{x}[\textit{i}-1]`, for :math:`\textit{i} = \mathrm{mfirst},\ldots,\textit{mlast}`, must contain the :math:`x` values at which the evaluation is required.
**xmin** : float
The lower and upper ends, :math:`x_{\mathrm{min}}` and :math:`x_{\mathrm{max}}`, of the range of the variable :math:`x` (see :ref:`Notes <e02cb-py2-py-notes>`).
The values of :math:`\mathrm{xmin}` and :math:`\mathrm{xmax}` may depend on the value of :math:`y` (e.g., when the polynomial has been derived using :meth:`dim2_cheb_lines`).
**xmax** : float
The lower and upper ends, :math:`x_{\mathrm{min}}` and :math:`x_{\mathrm{max}}`, of the range of the variable :math:`x` (see :ref:`Notes <e02cb-py2-py-notes>`).
The values of :math:`\mathrm{xmin}` and :math:`\mathrm{xmax}` may depend on the value of :math:`y` (e.g., when the polynomial has been derived using :meth:`dim2_cheb_lines`).
**y** : float
The value of the :math:`y` coordinate of all the points at which the evaluation is required.
**ymin** : float
The lower and upper ends, :math:`y_{\mathrm{min}}` and :math:`y_{\mathrm{max}}`, of the range of the variable :math:`y` (see :ref:`Notes <e02cb-py2-py-notes>`).
**ymax** : float
The lower and upper ends, :math:`y_{\mathrm{min}}` and :math:`y_{\mathrm{max}}`, of the range of the variable :math:`y` (see :ref:`Notes <e02cb-py2-py-notes>`).
**a** : float, array-like, shape :math:`\left(\left(\mathrm{k}+1\right)\times \left(\mathrm{l}+1\right)\right)`
The Chebyshev coefficients of the polynomial. The coefficient :math:`a_{{ij}}` defined according to the standard convention (see :ref:`Notes <e02cb-py2-py-notes>`) must be in :math:`\mathrm{a}[i\times \left(l+1\right)+j]`.
**Returns**
**ff** : float, ndarray, shape :math:`\left(\textit{mlast}\right)`
:math:`\mathrm{ff}[\textit{i}-1]` gives the value of the polynomial at the point :math:`\left(x_{\textit{i}}, y\right)`, for :math:`\textit{i} = \mathrm{mfirst},\ldots,\textit{mlast}`.
.. _e02cb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{na} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{l} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{na}\geq \left(\mathrm{k}+1\right)\times \left(\mathrm{l}+1\right)`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{l} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{l}\geq 0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{k}\geq 0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{mfirst} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{mlast} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mfirst}\leq \textit{mlast}`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{y} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ymax} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{y}\leq \mathrm{ymax}`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{y} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ymin} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{y}\geq \mathrm{ymin}`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{ymin} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ymax} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ymin} < \mathrm{ymax}`.
(`errno` :math:`3`)
On entry, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[\textit{I}-1]\leq \mathrm{xmax}`.
(`errno` :math:`3`)
On entry, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[\textit{I}-1]\geq \mathrm{xmin}`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{xmin} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xmax} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xmin} < \mathrm{xmax}`.
(`errno` :math:`4`)
Unexpected internal failure when evaluating the polynomial.
.. _e02cb-py2-py-notes:
**Notes**
This function evaluates a bivariate polynomial (represented in double Chebyshev form) of degree :math:`k` in one variable, :math:`\bar{x}`, and degree :math:`l` in the other, :math:`\bar{y}`.
The range of both variables is :math:`-1` to :math:`+1`.
However, these normalized variables will usually have been derived (as when the polynomial has been computed by :meth:`dim2_cheb_lines`, for example) from your original variables :math:`x` and :math:`y` by the transformations
.. math::
\bar{x} = \frac{{2x-\left(x_{\mathrm{max}}+x_{\mathrm{min}}\right)}}{\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)}\quad \text{ and }\quad \bar{y} = \frac{{2y-\left(y_{\mathrm{max}}+y_{\mathrm{min}}\right)}}{\left(y_{\mathrm{max}}-y_{\mathrm{min}}\right)}\text{.}
(Here :math:`x_{\mathrm{min}}` and :math:`x_{\mathrm{max}}` are the ends of the range of :math:`x` which has been transformed to the range :math:`-1` to :math:`+1` of :math:`\bar{x}`. :math:`y_{\mathrm{min}}` and :math:`y_{\mathrm{max}}` are correspondingly for :math:`y`.
See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02cbf.html#fcomments>`__).
For this reason, the function has been designed to accept values of :math:`x` and :math:`y` rather than :math:`\bar{x}` and :math:`\bar{y}`, and so requires values of :math:`x_{\mathrm{min}}`, etc. to be supplied by you.
In fact, for the sake of efficiency in appropriate cases, the function evaluates the polynomial for a sequence of values of :math:`x`, all associated with the same value of :math:`y`.
The double Chebyshev series can be written as
.. math::
\sum_{{i = 0}}^k\sum_{{j = 0}}^la_{{ij}}T_i\left(\bar{x}\right)T_j\left(\bar{y}\right)\text{,}
where :math:`T_i\left(\bar{x}\right)` is the Chebyshev polynomial of the first kind of degree :math:`i` and argument :math:`\bar{x}`, and :math:`T_j\left(\bar{y}\right)` is similarly defined.
However the standard convention, followed in this function, is that coefficients in the above expression which have either :math:`i` or :math:`j` zero are written :math:`\frac{1}{2}a_{{ij}}`, instead of simply :math:`a_{{ij}}`, and the coefficient with both :math:`i` and :math:`j` zero is written :math:`\frac{1}{4}a_{{0,0}}`.
The function first forms :math:`c_i = \sum_{{j = 0}}^la_{{ij}}T_j\left(\bar{y}\right)`, with :math:`a_{{i,0}}` replaced by :math:`\frac{1}{2}a_{{i,0}}`, for each of :math:`i = 0,1,\ldots,k`.
The value of the double series is then obtained for each value of :math:`x`, by summing :math:`c_i\times T_i\left(\bar{x}\right)`, with :math:`c_0` replaced by :math:`\frac{1}{2}c_0`, over :math:`i = 0,1,\ldots,k`.
The Clenshaw three term recurrence (see Clenshaw (1955)) with modifications due to Reinsch and Gentleman (1969) is used to form the sums.
.. _e02cb-py2-py-references:
**References**
Clenshaw, C W, 1955, `A note on the summation of Chebyshev series`, Math. Tables Aids Comput. (9), 118--120
Gentleman, W M, 1969, `An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients`, Comput. J. (12), 160--165
"""
raise NotImplementedError
[docs]def dim2_spline_panel(x, y, f, w, lamda, mu, point, eps):
r"""
``dim2_spline_panel`` forms a minimal, weighted least squares bicubic spline surface fit with prescribed knots to a given set of data points.
.. _e02da-py2-py-doc:
For full information please refer to the NAG Library document for e02da
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02daf.html
.. _e02da-py2-py-parameters:
**Parameters**
**x** : float, array-like, shape :math:`\left(m\right)`
The coordinates of the data point :math:`\left(x_{\textit{r}}, y_{\textit{r}}, f_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`. The order of the data points is immaterial, but see the array :math:`\mathrm{point}`.
**y** : float, array-like, shape :math:`\left(m\right)`
The coordinates of the data point :math:`\left(x_{\textit{r}}, y_{\textit{r}}, f_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`. The order of the data points is immaterial, but see the array :math:`\mathrm{point}`.
**f** : float, array-like, shape :math:`\left(m\right)`
The coordinates of the data point :math:`\left(x_{\textit{r}}, y_{\textit{r}}, f_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`. The order of the data points is immaterial, but see the array :math:`\mathrm{point}`.
**w** : float, array-like, shape :math:`\left(m\right)`
The weight :math:`w_r` of the :math:`r`\ th data point. It is important to note the definition of weight implied by the equation `(1) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02daf.html#eqn1>`__ in :ref:`Notes <e02da-py2-py-notes>`, since it is also common usage to define weight as the square of this weight. In this function, each :math:`w_r` should be chosen inversely proportional to the (absolute) accuracy of the corresponding :math:`f_r`, as expressed, for example, by the standard deviation or probable error of the :math:`f_r`. When the :math:`f_r` are all of the same accuracy, all the :math:`w_r` may be set equal to :math:`1.0`.
**lamda** : float, array-like, shape :math:`\left(\textit{px}\right)`
:math:`\mathrm{lamda}[\textit{i}+3]` must contain the :math:`\textit{i}`\ th interior knot :math:`\lambda_{{\textit{i}+4}}` associated with the variable :math:`x`, for :math:`\textit{i} = 1,2,\ldots,\textit{px}-8`. The knots must be in nondecreasing order and lie strictly within the range covered by the data values of :math:`x`. A knot is a value of :math:`x` at which the spline is allowed to be discontinuous in the third derivative with respect to :math:`x`, though continuous up to the second derivative. This degree of continuity can be reduced, if you require, by the use of coincident knots, provided that no more than four knots are chosen to coincide at any point. Two, or three, coincident knots allow loss of continuity in, respectively, the second and first derivative with respect to :math:`x` at the value of :math:`x` at which they coincide. Four coincident knots split the spline surface into two independent parts. For choice of knots see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02daf.html#fcomments>`__.
**mu** : float, array-like, shape :math:`\left(\textit{py}\right)`
:math:`\mathrm{mu}[\textit{i}+3]` must contain the :math:`\textit{i}`\ th interior knot :math:`\mu_{{\textit{i}+4}}` associated with the variable :math:`y`, for :math:`\textit{i} = 1,2,\ldots,\textit{py}-8`.
**point** : int, array-like, shape :math:`\left(\textit{npoint}\right)`
Indexing information usually provided by :meth:`dim2_spline_sort` which enables the data points to be accessed in the order which produces the advantageous matrix structure mentioned in :ref:`Notes <e02da-py2-py-notes>`. This order is such that, if the :math:`\left(x, y\right)` plane is thought of as being divided into rectangular panels by the two sets of knots, all data in a panel occur before data in succeeding panels, where the panels are numbered from bottom to top and then left to right with the usual arrangement of axes, as indicated in Figure [label omitted].
[figure omitted]
A data point lying exactly on one or more panel sides is considered to be in the highest numbered panel adjacent to the point. :meth:`dim2_spline_sort` should be called to obtain the array :math:`\mathrm{point}`, unless it is provided by other means.
**eps** : float
A threshold :math:`\epsilon` for determining the effective rank of the system of linear equations. The rank is determined as the number of elements of the array :math:`\mathrm{dl}` which are nonzero. An element of :math:`\mathrm{dl}` is regarded as zero if it is less than :math:`\epsilon`. Machine precision is a suitable value for :math:`\epsilon` in most practical applications which have only :math:`2` or :math:`3` decimals accurate in data. If some coefficients of the fit prove to be very large compared with the data ordinates, this suggests that :math:`\epsilon` should be increased so as to decrease the rank. The array :math:`\mathrm{dl}` will give a guide to appropriate values of :math:`\epsilon` to achieve this, as well as to the choice of :math:`\epsilon` in other cases where some experimentation may be needed to determine a value which leads to a satisfactory fit.
**Returns**
**lamda** : float, ndarray, shape :math:`\left(\textit{px}\right)`
The interior knots :math:`\mathrm{lamda}[4]` to :math:`\mathrm{lamda}[\textit{px}-5]` are unchanged, and the segments :math:`\mathrm{lamda}[0:4]` and :math:`\mathrm{lamda}[\textit{px}-4:\textit{px}]` contain additional (exterior) knots introduced by the function in order to define the full set of B-splines required. The four knots in the first segment are all set equal to the lowest data value of :math:`x` and the other four additional knots are all set equal to the highest value: there is experimental evidence that coincident end-knots are best for numerical accuracy. The complete array must be left undisturbed if :meth:`dim2_spline_evalv` or :meth:`dim2_spline_evalm` is to be used subsequently.
**mu** : float, ndarray, shape :math:`\left(\textit{py}\right)`
The same remarks apply to :math:`\mathrm{mu}` as to :math:`\mathrm{lamda}` above, with :math:`\mathrm{y}` replacing :math:`\mathrm{x}`, and :math:`y` replacing :math:`x`.
**dl** : float, ndarray, shape :math:`\left(\left(\textit{px}-4\right)\times \left(\textit{py}-4\right)\right)`
Gives the squares of the diagonal elements of the reduced triangular matrix, divided by the mean squared weight. It includes those elements, less than :math:`\epsilon`, which are treated as zero (see :ref:`Notes <e02da-py2-py-notes>`).
**c** : float, ndarray, shape :math:`\left(\left(\textit{px}-4\right)\times \left(\textit{py}-4\right)\right)`
Gives the coefficients of the fit. :math:`\mathrm{c}[ \left(\textit{py}-4\right)\times \left(\textit{i}-1\right) +\textit{j} -1]` is the coefficient :math:`c_{{\textit{i}\textit{j}}}` of :ref:`Notes <e02da-py2-py-notes>` and `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02daf.html#fcomments>`__, for :math:`\textit{j} = 1,2,\ldots,\textit{py}-4`, for :math:`\textit{i} = 1,2,\ldots,\textit{px}-4`. These coefficients are used by :meth:`dim2_spline_evalv` or :meth:`dim2_spline_evalm` to calculate values of the fitted function.
**sigma** : float
:math:`\Sigma`, the weighted sum of squares of residuals. This is not computed from the individual residuals but from the right-hand sides of the orthogonally-transformed linear equations. For further details see page 97 of Hayes and Halliday (1974). The two methods of computation are theoretically equivalent, but the results may differ because of rounding error.
**rank** : int
The rank of the system as determined by the value of the threshold :math:`\epsilon`.
:math:`\mathrm{rank} = \left(\textit{px}-4\right)\times \left(\textit{py}-4\right)`
The least squares solution is unique.
:math:`\mathrm{rank}\neq \left(\textit{px}-4\right)\times \left(\textit{py}-4\right)`
The minimal least squares solution is computed.
.. _e02da-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
At least one set of knots is not in nondecreasing order.
(`errno` :math:`2`)
More than four knots coincide at a single point.
(`errno` :math:`3`)
Array :math:`\mathrm{point}` does not indicate the data points in panel order.
(`errno` :math:`4`)
On entry, :math:`\textit{npoint} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{px} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{py} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{npoint}\geq m+\left(\textit{px}-7\right)\times \left(\textit{py}-7\right)`.
(`errno` :math:`4`)
On entry, :math:`\textit{nc} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{px} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{py} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nc} = \left(\textit{px}-4\right)\times \left(\textit{py}-4\right)`.
(`errno` :math:`4`)
On entry, :math:`\textit{py} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{py}\geq 8`.
(`errno` :math:`4`)
On entry, :math:`\textit{px} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{px}\geq 8`.
(`errno` :math:`4`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m > 1`.
(`errno` :math:`5`)
All the weights are zero, or rank determined as zero.
.. _e02da-py2-py-notes:
**Notes**
``dim2_spline_panel`` determines a bicubic spline fit :math:`s\left(x, y\right)` to the set of data points :math:`\left(x_r, y_r, f_r\right)` with weights :math:`w_r`, for :math:`\textit{r} = 1,2,\ldots,m`.
The two sets of internal knots of the spline, :math:`\left\{\lambda \right\}` and :math:`\left\{\mu \right\}`, associated with the variables :math:`x` and :math:`y` respectively, are prescribed by you.
These knots can be thought of as dividing the data region of the :math:`\left(x, y\right)` plane into panels (see Figure [label omitted] in :ref:`Parameters <e02da-py2-py-parameters>`).
A bicubic spline consists of a separate bicubic polynomial in each panel, the polynomials joining together with continuity up to the second derivative across the panel boundaries.
:math:`s\left(x, y\right)` has the property that :math:`\Sigma`, the sum of squares of its weighted residuals :math:`\rho_r`, for :math:`\textit{r} = 1,2,\ldots,m`, where
.. math::
\rho_r = w_r\left(s\left(x_r, y_r\right)-f_r\right)
is as small as possible for a bicubic spline with the given knot sets.
The function produces this minimized value of :math:`\Sigma` and the coefficients :math:`c_{{ij}}` in the B-spline representation of :math:`s\left(x, y\right)` -- see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02daf.html#fcomments>`__. :meth:`dim2_spline_evalv`, :meth:`dim2_spline_evalm` and :meth:`dim2_spline_derivm` are available to compute values and derivatives of the fitted spline from the coefficients :math:`c_{{ij}}`.
The least squares criterion is not always sufficient to determine the bicubic spline uniquely: there may be a whole family of splines which have the same minimum sum of squares.
In these cases, the function selects from this family the spline for which the sum of squares of the coefficients :math:`c_{{ij}}` is smallest: in other words, the minimal least squares solution.
This choice, although arbitrary, reduces the risk of unwanted fluctuations in the spline fit.
The method employed involves forming a system of :math:`m` linear equations in the coefficients :math:`c_{{ij}}` and then computing its least squares solution, which will be the minimal least squares solution when appropriate.
The basis of the method is described in Hayes and Halliday (1974).
The matrix of the equation is formed using a recurrence relation for B-splines which is numerically stable (see Cox (1972) and de Boor (1972) -- the former contains the more elementary derivation but, unlike de Boor (1972), does not cover the case of coincident knots).
The least squares solution is also obtained in a stable manner by using orthogonal transformations, viz. a variant of Givens rotation (see Gentleman (1973)).
This requires only one row of the matrix to be stored at a time.
Advantage is taken of the stepped-band structure which the matrix possesses when the data points are suitably ordered, there being at most sixteen nonzero elements in any row because of the definition of B-splines.
First the matrix is reduced to upper triangular form and then the diagonal elements of this triangle are examined in turn.
When an element is encountered whose square, divided by the mean squared weight, is less than a threshold :math:`\epsilon`, it is replaced by zero and the rest of the elements in its row are reduced to zero by rotations with the remaining rows.
The rank of the system is taken to be the number of nonzero diagonal elements in the final triangle, and the nonzero rows of this triangle are used to compute the minimal least squares solution.
If all the diagonal elements are nonzero, the rank is equal to the number of coefficients :math:`c_{{ij}}` and the solution obtained is the ordinary least squares solution, which is unique in this case.
.. _e02da-py2-py-references:
**References**
Cox, M G, 1972, `The numerical evaluation of B-splines`, J. Inst. Math. Appl. (10), 134--149
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
Gentleman, W M, 1973, `Least squares computations by Givens transformations without square roots`, J. Inst. Math. Applic. (12), 329--336
Hayes, J G and Halliday, J, 1974, `The least squares fitting of cubic spline surfaces to general data sets`, J. Inst. Math. Appl. (14), 89--103
"""
raise NotImplementedError
[docs]def dim2_spline_grid(start, x, y, f, s, nx, lamda, ny, mu, comm):
r"""
``dim2_spline_grid`` computes a bicubic spline approximation to a set of data values, given on a rectangular grid in the :math:`x`-:math:`y` plane.
The knots of the spline are located automatically, but a single argument must be specified to control the trade-off between closeness of fit and smoothness of fit.
.. _e02dc-py2-py-doc:
For full information please refer to the NAG Library document for e02dc
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dcf.html
.. _e02dc-py2-py-parameters:
**Parameters**
**start** : str, length 1
Determines whether calculations are to be performed afresh (Cold Start) or whether knots found in previous calls are to be used as an initial estimate of knot placement (Warm Start).
:math:`\mathrm{start} = \texttt{'C'}`
The function will build up the knot set starting with no interior knots. No values need be assigned to the arguments :math:`\mathrm{nx}`, :math:`\mathrm{ny}`, :math:`\mathrm{lamda}`, :math:`\mathrm{mu}`, :math:`\mathrm{comm}`\ ['wrk'] or :math:`\mathrm{comm}`\ ['iwrk'].
:math:`\mathrm{start} = \texttt{'W'}`
The function will restart the knot-placing strategy using the knots found in a previous call of the function. In this case, the arguments :math:`\mathrm{nx}`, :math:`\mathrm{ny}`, :math:`\mathrm{lamda}`, :math:`\mathrm{mu}`, :math:`\mathrm{comm}`\ ['wrk'] and :math:`\mathrm{comm}`\ ['iwrk'] must be unchanged from that previous call. This warm start can save much time in determining a satisfactory set of knots for the given value of :math:`\mathrm{s}`. This is particularly useful when different smoothing factors are used for the same dataset.
**x** : float, array-like, shape :math:`\left(\textit{mx}\right)`
:math:`\mathrm{x}[\textit{q}-1]` must be set to :math:`x_{\textit{q}}`, the :math:`x` coordinate of the :math:`\textit{q}`\ th grid point along the :math:`x` axis, for :math:`\textit{q} = 1,2,\ldots,m_x`.
**y** : float, array-like, shape :math:`\left(\textit{my}\right)`
:math:`\mathrm{y}[\textit{r}-1]` must be set to :math:`y_{\textit{r}}`, the :math:`y` coordinate of the :math:`\textit{r}`\ th grid point along the :math:`y` axis, for :math:`\textit{r} = 1,2,\ldots,m_y`.
**f** : float, array-like, shape :math:`\left(\textit{mx}\times \textit{my}\right)`
:math:`\mathrm{f}[m_y\times \left(\textit{q}-1\right)+\textit{r}-1]` must contain the data value :math:`f_{{\textit{q},\textit{r}}}`, for :math:`\textit{r} = 1,2,\ldots,m_y`, for :math:`\textit{q} = 1,2,\ldots,m_x`.
**s** : float
The smoothing factor, :math:`S`.
If :math:`\mathrm{s} = 0.0`, the function returns an interpolating spline.
If :math:`\mathrm{s}` is smaller than machine precision, it is assumed equal to zero.
For advice on the choice of :math:`\mathrm{s}`, see :ref:`Notes <e02dc-py2-py-notes>` and `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dcf.html#fcomments3>`__.
**nx** : int
If the warm start option is used, the value of :math:`\mathrm{nx}` must be left unchanged from the previous call.
**lamda** : float, array-like, shape :math:`\left(\textit{nxest}\right)`
If the warm start option is used, the values :math:`\mathrm{lamda}[0],\mathrm{lamda}[1],\ldots,\mathrm{lamda}[\mathrm{nx}-1]` must be left unchanged from the previous call.
**ny** : int
If the warm start option is used, the value of :math:`\mathrm{ny}` must be left unchanged from the previous call.
**mu** : float, array-like, shape :math:`\left(\textit{nyest}\right)`
If the warm start option is used, the values :math:`\mathrm{mu}[0],\mathrm{mu}[1],\ldots,\mathrm{mu}[\mathrm{ny}-1]` must be left unchanged from the previous call.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**Returns**
**nx** : int
The total number of knots, :math:`n_x`, of the computed spline with respect to the :math:`x` variable.
**lamda** : float, ndarray, shape :math:`\left(\textit{nxest}\right)`
Contains the complete set of knots :math:`\lambda_i` associated with the :math:`x` variable, i.e., the interior knots :math:`\mathrm{lamda}[4],\mathrm{lamda}[5],\ldots,\mathrm{lamda}[\mathrm{nx}-5]` as well as the additional knots
.. math::
\mathrm{lamda}[0] = \mathrm{lamda}[1] = \mathrm{lamda}[2] = \mathrm{lamda}[3] = \mathrm{x}[0]
and
.. math::
\mathrm{lamda}[\mathrm{nx}-4] = \mathrm{lamda}[\mathrm{nx}-3] = \mathrm{lamda}[\mathrm{nx}-2] = \mathrm{lamda}[\mathrm{nx}-1] = \mathrm{x}[\textit{mx}-1]
needed for the B-spline representation.
**ny** : int
The total number of knots, :math:`n_y`, of the computed spline with respect to the :math:`y` variable.
**mu** : float, ndarray, shape :math:`\left(\textit{nyest}\right)`
Contains the complete set of knots :math:`\mu_i` associated with the :math:`y` variable, i.e., the interior knots :math:`\mathrm{mu}[4],\mathrm{mu}[5],\ldots,\mathrm{mu}[\mathrm{ny}-5]` as well as the additional knots
.. math::
\mathrm{mu}[0] = \mathrm{mu}[1] = \mathrm{mu}[2] = \mathrm{mu}[3] = \mathrm{y}[0]
and
.. math::
\mathrm{mu}[\mathrm{ny}-4] = \mathrm{mu}[\mathrm{ny}-3] = \mathrm{mu}[\mathrm{ny}-2] = \mathrm{mu}[\mathrm{ny}-1] = \mathrm{y}[\textit{my}-1]
needed for the B-spline representation.
**c** : float, ndarray, shape :math:`\left(\left(\textit{nxest}-4\right)\times \left(\textit{nyest}-4\right)\right)`
The coefficients of the spline approximation. :math:`\mathrm{c}[\left(n_y-4\right) \times \left(i-1\right) +j -1]` is the coefficient :math:`c_{{ij}}` defined in :ref:`Notes <e02dc-py2-py-notes>`.
**fp** : float
The sum of squared residuals, :math:`\theta`, of the computed spline approximation. If :math:`\mathrm{fp} = 0.0`, this is an interpolating spline. :math:`\mathrm{fp}` should equal :math:`\mathrm{s}` within a relative tolerance of :math:`0.001` unless :math:`\mathrm{nx} = \mathrm{ny} = 8`, when the spline has no interior knots and so is simply a bicubic polynomial. For knots to be inserted, :math:`\mathrm{s}` must be set to a value below the value of :math:`\mathrm{fp}` produced in this case.
.. _e02dc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{nxest} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{mx} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{nyest} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{my} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{s} = 0.0`, :math:`\textit{nxest}\geq \textit{mx}+4` and :math:`\textit{nyest}\geq \textit{my}+4`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{s}\geq 0.0`.
(`errno` :math:`1`)
On entry, :math:`\textit{nyest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nyest}\geq \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\textit{my} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{my}\geq \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\textit{nxest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nxest}\geq \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\textit{mx} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{mx}\geq \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{start} \neq \texttt{'W'}` or :math:`\texttt{'C'}`: :math:`\mathrm{start} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`2`)
On entry, :math:`q = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[q-2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[q-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[q-2] < \mathrm{x}[q-1]` for all :math:`q`.
(`errno` :math:`3`)
On entry, :math:`r = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{y}[r-2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{y}[r-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{y}[r-2] < \mathrm{y}[r-1]` for all :math:`r`.
(`errno` :math:`4`)
The number of knots needed in one direction is greater than :math:`\textit{nxest}` or :math:`\textit{nyest}`: :math:`\textit{nxest}`, :math:`\textit{nyest} = \langle\mathit{\boldsymbol{value}}\rangle,\langle\mathit{\boldsymbol{value}}\rangle`. Possibly :math:`\mathrm{s}` is too small: :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
The iterative process has failed to converge. Possibly :math:`\mathrm{s}` is too small: :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _e02dc-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim2_spline_grid`` determines a smooth bicubic spline approximation :math:`s\left(x, y\right)` to the set of data points :math:`\left(x_{\textit{q}}, y_{\textit{r}}, f_{{\textit{q},\textit{r}}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m_y`, for :math:`\textit{q} = 1,2,\ldots,m_x`.
The spline is given in the B-spline representation
.. math::
s\left(x, y\right) = \sum_{{i = 1}}^{{n_x-4}}\sum_{{j = 1}}^{{n_y-4}}c_{{ij}}M_i\left(x\right)N_j\left(y\right)\text{,}
where :math:`M_i\left(x\right)` and :math:`N_j\left(y\right)` denote normalized cubic B-splines, the former defined on the knots :math:`\lambda_i` to :math:`\lambda_{{i+4}}` and the latter on the knots :math:`\mu_j` to :math:`\mu_{{j+4}}`.
For further details, see Hayes and Halliday (1974) for bicubic splines and de Boor (1972) for normalized B-splines.
The total numbers :math:`n_x` and :math:`n_y` of these knots and their values :math:`\lambda_1,\ldots,\lambda_{n_x}` and :math:`\mu_1,\ldots,\mu_{n_y}` are chosen automatically by the function.
The knots :math:`\lambda_5,\ldots,\lambda_{{n_x-4}}` and :math:`\mu_5,\ldots,\mu_{{n_y-4}}` are the interior knots; they divide the approximation domain :math:`\left[x_1, {x_{m_x}}\right]\times \left[y_1, {y_{m_y}}\right]` into :math:`\left(n_x-7\right)\times \left(n_y-7\right)` subpanels :math:`\left[{\lambda_{\textit{i}}}, {\lambda_{{\textit{i}+1}}}\right]\times \left[\mu_{\textit{j}}, \mu_{{\textit{j}+1}}\right]`, for :math:`\textit{j} = 4,5,\ldots,n_y-4`, for :math:`\textit{i} = 4,5,\ldots,n_x-4`.
Then, much as in the curve case (see :meth:`dim1_spline_auto`), the coefficients :math:`c_{{\textit{i}\textit{j}}}` are determined as the solution of the following constrained minimization problem:
.. math::
\mathrm{minimize}\left(\eta \right)\text{,}
subject to the constraint
.. math::
\theta = \sum_{{q = 1}}^{m_x}\sum_{{r = 1}}^{m_y}\epsilon_{{q,r}}^2\leq S\text{,}
.. rst-class:: nag-rules-none nag-align-left
+-----+------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|where|:math:`\eta` |is a measure of the (lack of) smoothness of :math:`s\left(x, y\right)`. Its value depends on the discontinuity jumps in :math:`s\left(x, y\right)` across the boundaries of the subpanels. It is zero only when there are no discontinuities and is positive otherwise, increasing with the size of the jumps (see Dierckx (1982) for details).|
+-----+------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| |:math:`\epsilon_{{q,r}}`|denotes the residual :math:`f_{{q,r}}-s\left(x_q, y_r\right)`, |
+-----+------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|and |:math:`S` |is a non-negative number specified by you. |
+-----+------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
By means of the argument :math:`S`, 'the smoothing factor', you will then control the balance between smoothness and closeness of fit, as measured by the sum of squares of residuals in `(3) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dcf.html#eqn3>`__.
If :math:`S` is too large, the spline will be too smooth and signal will be lost (underfit); if :math:`S` is too small, the spline will pick up too much noise (overfit).
In the extreme cases the function will return an interpolating spline :math:`\left(\theta = 0\right)` if :math:`S` is set to zero, and the least squares bicubic polynomial :math:`\left(\eta = 0\right)` if :math:`S` is set very large.
Experimenting with :math:`S`-values between these two extremes should result in a good compromise. (See `Choice of s <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dcf.html#fcomments3>`__ for advice on choice of :math:`S`.)
The method employed is outlined in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dcf.html#fcomments5>`__ and fully described in Dierckx (1981) and Dierckx (1982).
It involves an adaptive strategy for locating the knots of the bicubic spline (depending on the function underlying the data and on the value of :math:`S`), and an iterative method for solving the constrained minimization problem once the knots have been determined.
Values and derivatives of the computed spline can subsequently be computed by calling :meth:`dim2_spline_evalv`, :meth:`dim2_spline_evalm` or :meth:`dim2_spline_derivm` as described in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dcf.html#fcomments6>`__.
.. _e02dc-py2-py-references:
**References**
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
Dierckx, P, 1981, `An improved algorithm for curve fitting with spline functions`, Report TW54, Department of Computer Science, Katholieke Univerciteit Leuven
Dierckx, P, 1982, `A fast algorithm for smoothing data on a rectangular grid while using spline functions`, SIAM J. Numer. Anal. (19), 1286--1304
Hayes, J G and Halliday, J, 1974, `The least squares fitting of cubic spline surfaces to general data sets`, J. Inst. Math. Appl. (14), 89--103
Reinsch, C H, 1967, `Smoothing by spline functions`, Numer. Math. (10), 177--183
"""
raise NotImplementedError
[docs]def dim2_spline_sctr(start, x, y, f, w, s, nx, lamda, ny, mu, comm):
r"""
``dim2_spline_sctr`` computes a bicubic spline approximation to a set of scattered data.
The knots of the spline are located automatically, but a single argument must be specified to control the trade-off between closeness of fit and smoothness of fit.
.. _e02dd-py2-py-doc:
For full information please refer to the NAG Library document for e02dd
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ddf.html
.. _e02dd-py2-py-parameters:
**Parameters**
**start** : str, length 1
Determines whether calculations are to be performed afresh (Cold Start) or whether knots found in previous calls are to be used as an initial estimate of knot placement (Warm Start).
:math:`\mathrm{start} = \texttt{'C'}`
The function will build up the knot set starting with no interior knots. No values need be assigned to the arguments :math:`\mathrm{nx}`, :math:`\mathrm{ny}`, :math:`\mathrm{lamda}`, :math:`\mathrm{mu}` or :math:`\mathrm{comm}`\ ['wrk'].
:math:`\mathrm{start} = \texttt{'W'}`
The function will restart the knot-placing strategy using the knots found in a previous call of the function. In this case, the arguments :math:`\mathrm{nx}`, :math:`\mathrm{ny}`, :math:`\mathrm{lamda}`, :math:`\mathrm{mu}` and :math:`\mathrm{comm}`\ ['wrk'] must be unchanged from that previous call. This warm start can save much time in determining a satisfactory set of knots for the given value of :math:`\mathrm{s}`. This is particularly useful when different smoothing factors are used for the same dataset.
**x** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{x}[\textit{r}-1]`, :math:`\mathrm{y}[\textit{r}-1]`, :math:`\mathrm{f}[\textit{r}-1]` must be set to the coordinates of :math:`\left(x_{\textit{r}}, y_{\textit{r}}, f_{\textit{r}}\right)`, the :math:`\textit{r}`\ th data point, for :math:`\textit{r} = 1,2,\ldots,m`. The order of the data points is immaterial.
**y** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{x}[\textit{r}-1]`, :math:`\mathrm{y}[\textit{r}-1]`, :math:`\mathrm{f}[\textit{r}-1]` must be set to the coordinates of :math:`\left(x_{\textit{r}}, y_{\textit{r}}, f_{\textit{r}}\right)`, the :math:`\textit{r}`\ th data point, for :math:`\textit{r} = 1,2,\ldots,m`. The order of the data points is immaterial.
**f** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{x}[\textit{r}-1]`, :math:`\mathrm{y}[\textit{r}-1]`, :math:`\mathrm{f}[\textit{r}-1]` must be set to the coordinates of :math:`\left(x_{\textit{r}}, y_{\textit{r}}, f_{\textit{r}}\right)`, the :math:`\textit{r}`\ th data point, for :math:`\textit{r} = 1,2,\ldots,m`. The order of the data points is immaterial.
**w** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{w}[\textit{r}-1]` must be set to :math:`\textit{w}_{\textit{r}}`, the :math:`\textit{r}`\ th value in the set of weights, for :math:`\textit{r} = 1,2,\ldots,m`. Zero weights are permitted and the corresponding points are ignored, except when determining :math:`x_{\mathrm{min}}`, :math:`x_{\mathrm{max}}`, :math:`y_{\mathrm{min}}` and :math:`y_{\mathrm{max}}` (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ddf.html#fcomments4>`__). For advice on the choice of weights, see `the E02 Introduction <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02intro.html#background12>`__.
**s** : float
The smoothing factor, :math:`\mathrm{s}`.
For advice on the choice of :math:`\mathrm{s}`, see :ref:`Notes <e02dd-py2-py-notes>` and `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ddf.html#fcomments2>`__.
**nx** : int
If the warm start option is used, the value of :math:`\mathrm{nx}` must be left unchanged from the previous call.
**lamda** : float, array-like, shape :math:`\left(\textit{nxest}\right)`
If the warm start option is used, the values :math:`\mathrm{lamda}[0],\mathrm{lamda}[1],\ldots,\mathrm{lamda}[\mathrm{nx}-1]` must be left unchanged from the previous call.
**ny** : int
If the warm start option is used, the value of :math:`\mathrm{ny}` must be left unchanged from the previous call.
**mu** : float, array-like, shape :math:`\left(\textit{nyest}\right)`
If the warm start option is used, the values :math:`\mathrm{mu}[0],\mathrm{mu}[1],\ldots,\mathrm{mu}[\mathrm{ny}-1]` must be left unchanged from the previous call.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**Returns**
**nx** : int
The total number of knots, :math:`n_x`, of the computed spline with respect to the :math:`x` variable.
**lamda** : float, ndarray, shape :math:`\left(\textit{nxest}\right)`
Contains the complete set of knots :math:`\lambda_i` associated with the :math:`x` variable, i.e., the interior knots :math:`\mathrm{lamda}[4],\mathrm{lamda}[5],\ldots,\mathrm{lamda}[\mathrm{nx}-5]` as well as the additional knots
.. math::
\mathrm{lamda}[0] = \mathrm{lamda}[1] = \mathrm{lamda}[2] = \mathrm{lamda}[3] = x_{\mathrm{min}}
and
.. math::
\mathrm{lamda}[\mathrm{nx}-4] = \mathrm{lamda}[\mathrm{nx}-3] = \mathrm{lamda}[\mathrm{nx}-2] = \mathrm{lamda}[\mathrm{nx}-1] = x_{\mathrm{max}}
needed for the B-spline representation (where :math:`x_{\mathrm{min}}` and :math:`x_{\mathrm{max}}` are as described in :ref:`Notes <e02dd-py2-py-notes>`).
**ny** : int
The total number of knots, :math:`n_y`, of the computed spline with respect to the :math:`y` variable.
**mu** : float, ndarray, shape :math:`\left(\textit{nyest}\right)`
Contains the complete set of knots :math:`\mu_i` associated with the :math:`y` variable, i.e., the interior knots :math:`\mathrm{mu}[4],\mathrm{mu}[5],\ldots,\mathrm{mu}[\mathrm{ny}-5]` as well as the additional knots
.. math::
\mathrm{mu}[0] = \mathrm{mu}[1] = \mathrm{mu}[2] = \mathrm{mu}[3] = y_{\mathrm{min}}
and
.. math::
\mathrm{mu}[\mathrm{ny}-4] = \mathrm{mu}[\mathrm{ny}-3] = \mathrm{mu}[\mathrm{ny}-2] = \mathrm{mu}[\mathrm{ny}-1] = y_{\mathrm{max}}
needed for the B-spline representation (where :math:`y_{\mathrm{min}}` and :math:`y_{\mathrm{max}}` are as described in :ref:`Notes <e02dd-py2-py-notes>`).
**c** : float, ndarray, shape :math:`\left(\left(\textit{nxest}-4\right)\times \left(\textit{nyest}-4\right)\right)`
The coefficients of the spline approximation. :math:`\mathrm{c}[\left(n_y-4\right) \times \left(i-1\right) +j -1]` is the coefficient :math:`c_{{ij}}` defined in :ref:`Notes <e02dd-py2-py-notes>`.
**fp** : float
The weighted sum of squared residuals, :math:`\theta`, of the computed spline approximation. :math:`\mathrm{fp}` should equal :math:`\mathrm{s}` within a relative tolerance of :math:`0.001` unless :math:`\mathrm{nx} = \mathrm{ny} = 8`, when the spline has no interior knots and so is simply a bicubic polynomial. For knots to be inserted, :math:`\mathrm{s}` must be set to a value below the value of :math:`\mathrm{fp}` produced in this case.
**rank** : int
Gives the rank of the system of equations used to compute the final spline (as determined by a suitable machine-dependent threshold). When :math:`\mathrm{rank} = \left(\mathrm{nx}-4\right)\times \left(\mathrm{ny}-4\right)`, the solution is unique; otherwise the system is rank-deficient and the minimum-norm solution is computed. The latter case may be caused by too small a value of :math:`\mathrm{s}`.
.. _e02dd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{nyest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nyest}\geq \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\textit{nxest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nxest}\geq \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{start} \neq \texttt{'W'}` or :math:`\texttt{'C'}`: :math:`\mathrm{start} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{s} > 0.0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{MM} < 16`, where :math:`\mathrm{MM}` is the number of points with nonzero weight: :math:`\mathrm{MM} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`2`)
On entry, all the values in array :math:`\mathrm{y}` are equal.
(`errno` :math:`2`)
On entry, all the values in array :math:`\mathrm{x}` are equal.
(`errno` :math:`3`)
The number of knots needed in one direction is greater than :math:`\textit{nxest}` or :math:`\textit{nyest}`: :math:`\textit{nxest}`, :math:`\textit{nyest} = \langle\mathit{\boldsymbol{value}}\rangle,\langle\mathit{\boldsymbol{value}}\rangle`. Possibly :math:`\mathrm{s}` is too small: :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`4`)
No more knots added; the number of :math:`B`-spline coefficients already exceeds :math:`\textit{m}`. Either :math:`\textit{m}` or :math:`\mathrm{s}` is probably too small: :math:`m = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
No more knots added; the additional knot would coincide with an old one. Possibly an inaccurate data point has too large a weight, or :math:`\mathrm{s}` is too small: :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`6`)
The iterative process has failed to converge. Possibly :math:`\mathrm{s}` is too small: :math:`\mathrm{s} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _e02dd-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim2_spline_sctr`` determines a smooth bicubic spline approximation :math:`s\left(x, y\right)` to the set of data points :math:`\left(x_{\textit{r}}, y_{\textit{r}}, f_{\textit{r}}\right)` with weights :math:`\textit{w}_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`.
The approximation domain is considered to be the rectangle :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]\times \left[y_{\mathrm{min}}, y_{\mathrm{max}}\right]`, where :math:`x_{\mathrm{min}}\left(y_{\mathrm{min}}\right)` and :math:`x_{\mathrm{max}}\left(y_{\mathrm{max}}\right)` denote the lowest and highest data values of :math:`x\left(y\right)`.
The spline is given in the B-spline representation
.. math::
s\left(x, y\right) = \sum_{{i = 1}}^{{n_x-4}}\sum_{{j = 1}}^{{n_y-4}}c_{{ij}}M_i\left(x\right)N_j\left(y\right)\text{,}
where :math:`M_i\left(x\right)` and :math:`N_j\left(y\right)` denote normalized cubic B-splines, the former defined on the knots :math:`\lambda_i` to :math:`\lambda_{{i+4}}` and the latter on the knots :math:`\mu_j` to :math:`\mu_{{j+4}}`.
For further details, see Hayes and Halliday (1974) for bicubic splines and de Boor (1972) for normalized B-splines.
The total numbers :math:`n_x` and :math:`n_y` of these knots and their values :math:`\lambda_1,\ldots,\lambda_{n_x}` and :math:`\mu_1,\ldots,\mu_{n_y}` are chosen automatically by the function.
The knots :math:`\lambda_5,\ldots,\lambda_{{n_x-4}}` and :math:`\mu_5,\ldots,\mu_{{n_y-4}}` are the interior knots; they divide the approximation domain :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]\times \left[y_{\mathrm{min}}, y_{\mathrm{max}}\right]` into :math:`\left(n_x-7\right)\times \left(n_y-7\right)` subpanels :math:`\left[{\lambda_{\textit{i}}}, {\lambda_{{\textit{i}+1}}}\right]\times \left[\mu_{\textit{j}}, \mu_{{\textit{j}+1}}\right]`, for :math:`\textit{j} = 4,5,\ldots,n_y-4`, for :math:`\textit{i} = 4,5,\ldots,n_x-4`.
Then, much as in the curve case (see :meth:`dim1_spline_auto`), the coefficients :math:`c_{{ij}}` are determined as the solution of the following constrained minimization problem:
minimize
.. math::
\eta \text{,}
subject to the constraint
.. math::
\theta = \sum_{{r = 1}}^m\epsilon_r^2\leq S
.. rst-class:: nag-rules-none
+-----+------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|where|:math:`\eta` |is a measure of the (lack of) smoothness of :math:`s\left(x, y\right)`. Its value depends on the discontinuity jumps in :math:`s\left(x, y\right)` across the boundaries of the subpanels. It is zero only when there are no discontinuities and is positive otherwise, increasing with the size of the jumps (see Dierckx (1981b) for details).|
+-----+------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| |:math:`\epsilon_r`|denotes the weighted residual :math:`\textit{w}_r\left(f_r-s\left(x_r, y_r\right)\right)`, |
+-----+------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|and |:math:`\mathrm{s}`|is a non-negative number to be specified by you. |
+-----+------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
By means of the argument :math:`\mathrm{s}`, 'the smoothing factor', you will then control the balance between smoothness and closeness of fit, as measured by the sum of squares of residuals in `(3) <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ddf.html#eqn3>`__.
If :math:`\mathrm{s}` is too large, the spline will be too smooth and signal will be lost (underfit); if :math:`\mathrm{s}` is too small, the spline will pick up too much noise (overfit).
In the extreme cases the method would return an interpolating spline :math:`\left(\theta = 0\right)` if :math:`\mathrm{s}` were set to zero, and returns the least squares bicubic polynomial :math:`\left(\eta = 0\right)` if :math:`\mathrm{s}` is set very large.
Experimenting with :math:`\mathrm{s}`-values between these two extremes should result in a good compromise. (See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ddf.html#fcomments2>`__ for advice on choice of :math:`\mathrm{s}`.) Note however, that this function, unlike :meth:`dim1_spline_auto` and :meth:`dim2_spline_grid`, does not allow :math:`\mathrm{s}` to be set exactly to zero: to compute an interpolant to scattered data, :meth:`interp.dim2_scat <naginterfaces.library.interp.dim2_scat>` or :meth:`interp.dim2_scat_shep <naginterfaces.library.interp.dim2_scat_shep>` should be used.
The method employed is outlined in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ddf.html#fcomments5>`__ and fully described in Dierckx (1981a) and Dierckx (1981b).
It involves an adaptive strategy for locating the knots of the bicubic spline (depending on the function underlying the data and on the value of :math:`\mathrm{s}`), and an iterative method for solving the constrained minimization problem once the knots have been determined.
Values and derivatives of the computed spline can subsequently be computed by calling :meth:`dim2_spline_evalv`, :meth:`dim2_spline_evalm` or :meth:`dim2_spline_derivm` as described in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02ddf.html#fcomments6>`__.
.. _e02dd-py2-py-references:
**References**
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
Dierckx, P, 1981, `An algorithm for surface fitting with spline functions`, IMA J. Numer. Anal. (1), 267--283
Dierckx, P, 1981, `An improved algorithm for curve fitting with spline functions`, Report TW54, Department of Computer Science, Katholieke Univerciteit Leuven
Hayes, J G and Halliday, J, 1974, `The least squares fitting of cubic spline surfaces to general data sets`, J. Inst. Math. Appl. (14), 89--103
Peters, G and Wilkinson, J H, 1970, `The least squares problem and pseudo-inverses`, Comput. J. (13), 309--316
Reinsch, C H, 1967, `Smoothing by spline functions`, Numer. Math. (10), 177--183
"""
raise NotImplementedError
[docs]def dim2_spline_evalv(x, y, lamda, mu, c):
r"""
``dim2_spline_evalv`` calculates values of a bicubic spline from its B-spline representation.
.. _e02de-py2-py-doc:
For full information please refer to the NAG Library document for e02de
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02def.html
.. _e02de-py2-py-parameters:
**Parameters**
**x** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{x}` and :math:`\mathrm{y}` must contain :math:`x_{\textit{r}}` and :math:`y_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`, respectively. These are the coordinates of the points at which values of the spline are required. The order of the points is immaterial.
**y** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{x}` and :math:`\mathrm{y}` must contain :math:`x_{\textit{r}}` and :math:`y_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m`, respectively. These are the coordinates of the points at which values of the spline are required. The order of the points is immaterial.
**lamda** : float, array-like, shape :math:`\left(\textit{px}\right)`
:math:`\mathrm{lamda}` and :math:`\mathrm{mu}` must contain the complete sets of knots :math:`\left\{\lambda \right\}` and :math:`\left\{\mu \right\}` associated with the :math:`x` and :math:`y` variables respectively.
**mu** : float, array-like, shape :math:`\left(\textit{py}\right)`
:math:`\mathrm{lamda}` and :math:`\mathrm{mu}` must contain the complete sets of knots :math:`\left\{\lambda \right\}` and :math:`\left\{\mu \right\}` associated with the :math:`x` and :math:`y` variables respectively.
**c** : float, array-like, shape :math:`\left(\left(\textit{px}-4\right)\times \left(\textit{py}-4\right)\right)`
:math:`\mathrm{c}[ \left(\textit{py}-4\right) \times \left(\textit{i}-1\right) +\textit{j} -1]` must contain the coefficient :math:`c_{{\textit{i}\textit{j}}}` described in :ref:`Notes <e02de-py2-py-notes>`, for :math:`\textit{j} = 1,2,\ldots,\textit{py}-4`, for :math:`\textit{i} = 1,2,\ldots,\textit{px}-4`.
**Returns**
**ff** : float, ndarray, shape :math:`\left(m\right)`
:math:`\mathrm{ff}[\textit{r}-1]` contains the value of the spline at the point :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`.
.. _e02de-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{py} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{py}\geq 8`.
(`errno` :math:`1`)
On entry, :math:`\textit{px} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{px}\geq 8`.
(`errno` :math:`1`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 1`.
(`errno` :math:`2`)
On entry, the knots in :math:`\mathrm{mu}` are not in nondecreasing order.
(`errno` :math:`2`)
On entry, the knots in :math:`\mathrm{lamda}` are not in nondecreasing order.
(`errno` :math:`3`)
On entry, point :math:`\left({\mathrm{x}[\textit{K}-1]}, {\mathrm{y}[\textit{K}-1]}\right)` lies outside the rectangle bounded by :math:`\mathrm{lamda}[3]`, :math:`\mathrm{lamda}[\textit{px}-4]`, :math:`\mathrm{mu}[3]`, :math:`\mathrm{mu}[\textit{py}-4]`: :math:`\textit{K} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[\textit{K}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{y}[\textit{K}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _e02de-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim2_spline_evalv`` calculates values of the bicubic spline :math:`s\left(x, y\right)` at prescribed points :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`, from its augmented knot sets :math:`\left\{\lambda \right\}` and :math:`\left\{\mu \right\}` and from the coefficients :math:`c_{{ij}}`, for :math:`\textit{j} = 1,2,\ldots,\textit{py}-4`, for :math:`\textit{i} = 1,2,\ldots,\textit{px}-4`, in its B-spline representation
.. math::
s\left(x, y\right) = \sum_{{ij}}c_{{ij}}M_i\left(x\right)N_j\left(y\right)\text{.}
Here :math:`M_i\left(x\right)` and :math:`N_j\left(y\right)` denote normalized cubic B-splines, the former defined on the knots :math:`\lambda_i` to :math:`\lambda_{{i+4}}` and the latter on the knots :math:`\mu_j` to :math:`\mu_{{j+4}}`.
This function may be used to calculate values of a bicubic spline given in the form produced by :meth:`interp.dim2_spline_grid <naginterfaces.library.interp.dim2_spline_grid>`, :meth:`dim2_spline_panel`, :meth:`dim2_spline_grid` and :meth:`dim2_spline_sctr`.
It is derived from the function B2VRE in Anthony `et al.` (1982).
.. _e02de-py2-py-references:
**References**
Anthony, G T, Cox, M G and Hayes, J G, 1982, `DASL -- Data Approximation Subroutine Library`, National Physical Laboratory
Cox, M G, 1978, `The numerical evaluation of a spline from its B-spline representation`, J. Inst. Math. Appl. (21), 135--143
"""
raise NotImplementedError
[docs]def dim2_spline_evalm(x, y, lamda, mu, c):
r"""
``dim2_spline_evalm`` calculates values of a bicubic spline from its B-spline representation.
The spline is evaluated at all points on a rectangular grid.
.. _e02df-py2-py-doc:
For full information please refer to the NAG Library document for e02df
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dff.html
.. _e02df-py2-py-parameters:
**Parameters**
**x** : float, array-like, shape :math:`\left(\textit{mx}\right)`
:math:`\mathrm{x}` and :math:`\mathrm{y}` must contain :math:`x_{\textit{q}}`, for :math:`\textit{q} = 1,2,\ldots,m_x`, and :math:`y_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m_y`, respectively. These are the :math:`x` and :math:`y` coordinates that define the rectangular grid of points at which values of the spline are required.
**y** : float, array-like, shape :math:`\left(\textit{my}\right)`
:math:`\mathrm{x}` and :math:`\mathrm{y}` must contain :math:`x_{\textit{q}}`, for :math:`\textit{q} = 1,2,\ldots,m_x`, and :math:`y_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m_y`, respectively. These are the :math:`x` and :math:`y` coordinates that define the rectangular grid of points at which values of the spline are required.
**lamda** : float, array-like, shape :math:`\left(\textit{px}\right)`
:math:`\mathrm{lamda}` and :math:`\mathrm{mu}` must contain the complete sets of knots :math:`\left\{\lambda \right\}` and :math:`\left\{\mu \right\}` associated with the :math:`x` and :math:`y` variables respectively.
**mu** : float, array-like, shape :math:`\left(\textit{py}\right)`
:math:`\mathrm{lamda}` and :math:`\mathrm{mu}` must contain the complete sets of knots :math:`\left\{\lambda \right\}` and :math:`\left\{\mu \right\}` associated with the :math:`x` and :math:`y` variables respectively.
**c** : float, array-like, shape :math:`\left(\left(\textit{px}-4\right)\times \left(\textit{py}-4\right)\right)`
:math:`\mathrm{c}[ \left(\textit{py}-4\right) \times \left(\textit{i}-1\right) +\textit{j} -1]` must contain the coefficient :math:`c_{{\textit{i}\textit{j}}}` described in :ref:`Notes <e02df-py2-py-notes>`, for :math:`\textit{j} = 1,2,\ldots,\textit{py}-4`, for :math:`\textit{i} = 1,2,\ldots,\textit{px}-4`.
**Returns**
**ff** : float, ndarray, shape :math:`\left(\textit{mx}\times \textit{my}\right)`
:math:`\mathrm{ff}[\textit{my}\times \left(\textit{q}-1\right)+\textit{r}-1]` contains the value of the spline at the point :math:`\left(x_{\textit{q}}, y_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m_y`, for :math:`\textit{q} = 1,2,\ldots,m_x`.
.. _e02df-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{my} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{my}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{mx} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{mx}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{py} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{py}\geq 8`.
(`errno` :math:`1`)
On entry, :math:`\textit{px} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{px}\geq 8`.
(`errno` :math:`3`)
On entry, the knots in :math:`\mathrm{mu}` are not in nondecreasing order.
(`errno` :math:`3`)
On entry, the knots in :math:`\mathrm{lamda}` are not in nondecreasing order.
(`errno` :math:`4`)
On entry, :math:`\mathrm{mu}[3]\leq \mathrm{y}[0] < \cdots < \mathrm{y}[\textit{my}-1]\leq \mathrm{mu}[\textit{py}-4]` is violated.
(`errno` :math:`4`)
On entry, :math:`\mathrm{lamda}[3]\leq \mathrm{x}[0] < \cdots < \mathrm{x}[\textit{mx}-1]\leq \mathrm{lamda}[\textit{px}-4]` is violated.
.. _e02df-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``dim2_spline_evalm`` calculates values of the bicubic spline :math:`s\left(x, y\right)` on a rectangular grid of points in the :math:`x`-:math:`y` plane, from its augmented knot sets :math:`\left\{\lambda \right\}` and :math:`\left\{\mu \right\}` and from the coefficients :math:`c_{{ij}}`, for :math:`\textit{j} = 1,2,\ldots,\textit{py}-4`, for :math:`\textit{i} = 1,2,\ldots,\textit{px}-4`, in its B-spline representation
.. math::
s\left(x, y\right) = \sum_{{ij}}c_{{ij}}M_i\left(x\right)N_j\left(y\right)\text{.}
Here :math:`M_i\left(x\right)` and :math:`N_j\left(y\right)` denote normalized cubic B-splines, the former defined on the knots :math:`\lambda_i` to :math:`\lambda_{{i+4}}` and the latter on the knots :math:`\mu_j` to :math:`\mu_{{j+4}}`.
The points in the grid are defined by coordinates :math:`x_q`, for :math:`\textit{q} = 1,2,\ldots,m_x`, along the :math:`x` axis, and coordinates :math:`y_r`, for :math:`\textit{r} = 1,2,\ldots,m_y`, along the :math:`y` axis.
This function may be used to calculate values of a bicubic spline given in the form produced by :meth:`interp.dim2_spline_grid <naginterfaces.library.interp.dim2_spline_grid>`, :meth:`dim2_spline_panel`, :meth:`dim2_spline_grid` and :meth:`dim2_spline_sctr`.
It is derived from the function B2VRE in Anthony `et al.` (1982).
.. _e02df-py2-py-references:
**References**
Anthony, G T, Cox, M G and Hayes, J G, 1982, `DASL -- Data Approximation Subroutine Library`, National Physical Laboratory
Cox, M G, 1978, `The numerical evaluation of a spline from its B-spline representation`, J. Inst. Math. Appl. (21), 135--143
"""
raise NotImplementedError
[docs]def dim2_spline_derivm(x, y, lamda, mu, c, nux, nuy):
r"""
``dim2_spline_derivm`` computes the partial derivative (of order :math:`\nu_x`, :math:`\nu_y`), of a bicubic spline approximation to a set of data values, from its B-spline representation, at points on a rectangular grid in the :math:`x`-:math:`y` plane.
This function may be used to calculate derivatives of a bicubic spline given in the form produced by :meth:`interp.dim2_spline_grid <naginterfaces.library.interp.dim2_spline_grid>`, :meth:`dim2_spline_panel`, :meth:`dim2_spline_grid` and :meth:`dim2_spline_sctr`.
.. _e02dh-py2-py-doc:
For full information please refer to the NAG Library document for e02dh
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02dhf.html
.. _e02dh-py2-py-parameters:
**Parameters**
**x** : float, array-like, shape :math:`\left(\textit{mx}\right)`
:math:`\mathrm{x}[q-1]` must be set to :math:`x_{\textit{q}}`, the :math:`x` coordinate of the :math:`\textit{q}`\ th grid point along the :math:`x` axis, for :math:`\textit{q} = 1,2,\ldots,m_x`, on which values of the partial derivative are sought.
**y** : float, array-like, shape :math:`\left(\textit{my}\right)`
:math:`\mathrm{y}[\textit{r}-1]` must be set to :math:`y_{\textit{r}}`, the :math:`y` coordinate of the :math:`\textit{r}`\ th grid point along the :math:`y` axis, for :math:`\textit{r} = 1,2,\ldots,m_y` on which values of the partial derivative are sought.
**lamda** : float, array-like, shape :math:`\left(\textit{px}\right)`
Contains the position of the knots in the :math:`x`-direction of the bicubic spline approximation to be differentiated, e.g., :math:`\textit{lamda}` as returned by :meth:`dim2_spline_grid`.
**mu** : float, array-like, shape :math:`\left(\textit{py}\right)`
Contains the position of the knots in the :math:`y`-direction of the bicubic spline approximation to be differentiated, e.g., :math:`\textit{mu}` as returned by :meth:`dim2_spline_grid`.
**c** : float, array-like, shape :math:`\left(\left(\textit{px}-4\right)\times \left(\textit{py}-4\right)\right)`
The coefficients of the bicubic spline approximation to be differentiated, e.g., :math:`\textit{c}` as returned by :meth:`dim2_spline_grid`.
**nux** : int
Specifies the order, :math:`\nu_x` of the partial derivative in the :math:`x`-direction.
**nuy** : int
Specifies the order, :math:`\nu_y` of the partial derivative in the :math:`y`-direction.
**Returns**
**z** : float, ndarray, shape :math:`\left(\textit{mx}\times \textit{my}\right)`
:math:`\mathrm{z}[m_y\times \left(\textit{q}-1\right)+\textit{r}-1]` contains the derivative :math:`\frac{\partial^{{\nu_x+\nu_y}}}{{{\partial x}^{\nu_x}{\partial y}^{\nu_y}}}s\left(x_q, y_r\right)`, for :math:`\textit{r} = 1,2,\ldots,m_y`, for :math:`\textit{q} = 1,2,\ldots,m_x`.
.. _e02dh-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{nux} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{nux}\leq 2`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{nuy} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{nuy}\leq 2`.
(`errno` :math:`3`)
On entry, :math:`\textit{mx} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{mx}\geq 1`.
(`errno` :math:`4`)
On entry, :math:`\textit{my} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{my}\geq 1`.
(`errno` :math:`5`)
On entry, for :math:`i = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{x}[i-2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{x}[i-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{x}[\textit{i}-2]\leq \mathrm{x}[\textit{i}-1]`, for :math:`\textit{i} = 2,3,\ldots,\textit{mx}`.
(`errno` :math:`6`)
On entry, for :math:`i = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{y}[i-2] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{y}[i-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{y}[\textit{i}-2]\leq \mathrm{y}[\textit{i}-1]`, for :math:`\textit{i} = 2,3,\ldots,\textit{my}`.
.. _e02dh-py2-py-notes:
**Notes**
``dim2_spline_derivm`` determines the partial derivative :math:`\frac{\partial^{{\nu_x+\nu_y}}}{{\partial x^{\nu_x}\partial y^{\nu_y}}}` of a smooth bicubic spline approximation :math:`s\left(x, y\right)` at the set of data points :math:`\left(x_q, y_r\right)`.
The spline is given in the B-spline representation
.. math::
s\left(x, y\right) = \sum_{{i = 1}}^{{n_x-4}}\sum_{{j = 1}}^{{n_y-4}}c_{{ij}}M_i\left(x\right)N_j\left(y\right)\text{,}
where :math:`M_i\left(x\right)` and :math:`N_j\left(y\right)` denote normalized cubic B-splines, the former defined on the knots :math:`\lambda_i` to :math:`\lambda_{{i+4}}` and the latter on the knots :math:`\mu_j` to :math:`\mu_{{j+4}}`, with :math:`n_x` and :math:`n_y` the total numbers of knots of the computed spline with respect to the :math:`x` and :math:`y` variables respectively.
For further details, see Hayes and Halliday (1974) for bicubic splines and de Boor (1972) for normalized B-splines.
This function is suitable for B-spline representations returned by :meth:`interp.dim2_spline_grid <naginterfaces.library.interp.dim2_spline_grid>`, :meth:`dim2_spline_panel`, :meth:`dim2_spline_grid` and :meth:`dim2_spline_sctr`.
The partial derivatives can be up to order :math:`2` in each direction; thus the highest mixed derivative available is :math:`\frac{\partial^4}{{\partial x^2\partial y^2}}`.
The points in the grid are defined by coordinates :math:`x_{\textit{q}}`, for :math:`\textit{q} = 1,2,\ldots,m_x`, along the :math:`x` axis, and coordinates :math:`y_{\textit{r}}`, for :math:`\textit{r} = 1,2,\ldots,m_y`, along the :math:`y` axis.
.. _e02dh-py2-py-references:
**References**
de Boor, C, 1972, `On calculating with B-splines`, J. Approx. Theory (6), 50--62
Dierckx, P, 1981, `An improved algorithm for curve fitting with spline functions`, Report TW54, Department of Computer Science, Katholieke Univerciteit Leuven
Dierckx, P, 1982, `A fast algorithm for smoothing data on a rectangular grid while using spline functions`, SIAM J. Numer. Anal. (19), 1286--1304
Hayes, J G and Halliday, J, 1974, `The least squares fitting of cubic spline surfaces to general data sets`, J. Inst. Math. Appl. (14), 89--103
Reinsch, C H, 1967, `Smoothing by spline functions`, Numer. Math. (10), 177--183
"""
raise NotImplementedError
[docs]def glin_l1sol(a, b, toler=0.0):
r"""
``glin_l1sol`` calculates an :math:`l_1` solution to an overdetermined system of linear equations.
.. _e02ga-py2-py-doc:
For full information please refer to the NAG Library document for e02ga
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02gaf.html
.. _e02ga-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(m+2, \textit{nplus2}\right)`
:math:`\mathrm{a}[\textit{i}-1,\textit{j}-1]` must contain :math:`a_{{\textit{i}\textit{j}}}`, the element in the :math:`\textit{i}`\ th row and :math:`\textit{j}`\ th column of the matrix :math:`A`, for :math:`\textit{j} = 1,2,\ldots,n`, for :math:`\textit{i} = 1,2,\ldots,m`. The remaining elements need not be set.
**b** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{b}[\textit{i}-1]` must contain :math:`b_{\textit{i}}`, the :math:`\textit{i}`\ th element of the vector :math:`b`, for :math:`\textit{i} = 1,2,\ldots,m`.
**toler** : float, optional
A non-negative value. In general :math:`\mathrm{toler}` specifies a threshold below which numbers are regarded as zero. The recommended threshold value is :math:`\epsilon^{{2/3}}` where :math:`\epsilon` is the machine precision. The recommended value can be computed within the function by setting :math:`\mathrm{toler}` to zero. If premature termination occurs a larger value for :math:`\mathrm{toler}` may result in a valid solution.
**Returns**
**a** : float, ndarray, shape :math:`\left(m+2, \textit{nplus2}\right)`
Contains the last simplex tableau generated by the simplex method.
**b** : float, ndarray, shape :math:`\left(m\right)`
The :math:`\textit{i}`\ th residual :math:`r_{\textit{i}}` corresponding to the solution vector :math:`x`, for :math:`\textit{i} = 1,2,\ldots,m`.
**x** : float, ndarray, shape :math:`\left(\textit{nplus2}\right)`
:math:`\mathrm{x}[\textit{j}-1]` contains the :math:`\textit{j}`\ th element of the solution vector :math:`x`, for :math:`\textit{j} = 1,2,\ldots,n`. The elements :math:`\mathrm{x}[n]` and :math:`\mathrm{x}[n+1]` are unused.
**resid** : float
The sum of the absolute values of the residuals for the solution vector :math:`x`.
**irank** : int
The computed rank of the matrix :math:`A`.
**itera** : int
The number of iterations taken by the simplex method.
.. _e02ga-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`2`)
Premature termination due to rounding errors. Try using larger value of :math:`\mathrm{toler}`: :math:`\mathrm{toler} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`3`)
On entry, :math:`\textit{nplus2} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`3\leq \textit{nplus2}\leq m+2`.
(`errno` :math:`3`)
On entry, :math:`\textit{nplus2} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nplus2}\geq 3`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
An optimal solution has been obtained, but may not be unique.
(`errno` :math:`4`)
More than :math:`1000*n` iterations were performed. ``glin_l1sol`` has terminated without calculating a solution. The output data from the function is as computed on the last good iteration. Consider increasing the value of :math:`\mathrm{toler}`. Alternatively, :math:`A` may be ill conditioned---try scaling its columns.
.. _e02ga-py2-py-notes:
**Notes**
Given a matrix :math:`A` with :math:`m` rows and :math:`n` columns :math:`\left(m\geq n\right)` and a vector :math:`b` with :math:`m` elements, the function calculates an :math:`l_1` solution to the overdetermined system of equations
.. math::
Ax = b\text{.}
That is to say, it calculates a vector :math:`x`, with :math:`n` elements, which minimizes the :math:`l_1` norm (the sum of the absolute values) of the residuals
.. math::
r\left(x\right) = \sum_{{i = 1}}^m\left\lvert r_i\right\rvert \text{,}
where the residuals :math:`r_i` are given by
.. math::
r_i = b_i-\sum_{{j = 1}}^na_{{ij}}x_j\text{, }\quad i = 1,2,\ldots,m\text{.}
Here :math:`a_{{ij}}` is the element in row :math:`i` and column :math:`j` of :math:`A`, :math:`b_i` is the :math:`i`\ th element of :math:`b` and :math:`x_j` the :math:`j`\ th element of :math:`x`.
The matrix :math:`A` need not be of full rank.
Typically in applications to data fitting, data consisting of :math:`m` points with coordinates :math:`\left(t_i, y_i\right)` are to be approximated in the :math:`l_1` norm by a linear combination of known functions :math:`\phi_j\left(t\right)`,
.. math::
\alpha_1\phi_1\left(t\right)+\alpha_2\phi_2\left(t\right)+ \cdots +\alpha_n\phi_n\left(t\right)\text{.}
This is equivalent to fitting an :math:`l_1` solution to the overdetermined system of equations
.. math::
\sum_{{j = 1}}^n\phi_j\left(t_i\right)\alpha_j = y_i\text{, }\quad i = 1,2,\ldots,m\text{.}
Thus if, for each value of :math:`i` and :math:`j`, the element :math:`a_{{ij}}` of the matrix :math:`A` in the previous paragraph is set equal to the value of :math:`\phi_j\left(t_i\right)` and :math:`b_i` is set equal to :math:`y_i`, the solution vector :math:`x` will contain the required values of the :math:`\alpha_j`.
Note that the independent variable :math:`t` above can, instead, be a vector of several independent variables (this includes the case where each :math:`\phi_i` is a function of a different variable, or set of variables).
The algorithm is a modification of the simplex method of linear programming applied to the primal formulation of the :math:`l_1` problem (see Barrodale and Roberts (1973) and Barrodale and Roberts (1974)).
The modification allows several neighbouring simplex vertices to be passed through in a single iteration, providing a substantial improvement in efficiency.
.. _e02ga-py2-py-references:
**References**
Barrodale, I and Roberts, F D K, 1973, `An improved algorithm for discrete` :math:`l_1` `linear approximation`, SIAM J. Numer. Anal. (10), 839--848
Barrodale, I and Roberts, F D K, 1974, `Solution of an overdetermined system of equations in the` :math:`l_1` `-norm`, Comm. ACM (17(6)), 319--320
"""
raise NotImplementedError
[docs]def glinc_l1sol(m, e, f, x, mxs, monit, iprint, data=None):
r"""
``glinc_l1sol`` calculates an :math:`l_1` solution to an overdetermined system of linear equations, possibly subject to linear inequality constraints.
.. _e02gb-py2-py-doc:
For full information please refer to the NAG Library document for e02gb
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02gbf.html
.. _e02gb-py2-py-parameters:
**Parameters**
**m** : int
The number of equations in the overdetermined system, :math:`m` (i.e., the number of rows of the matrix :math:`A`).
**e** : float, array-like, shape :math:`\left(n, \textit{mpl}\right)`
The equation and constraint matrices stored in the following manner.
The first :math:`m` columns contain the :math:`m` rows of the matrix :math:`A`; element :math:`\mathrm{e}[\textit{i}-1,\textit{j}-1]` specifying the element :math:`a_{{\textit{j}\textit{i}}}` in the :math:`\textit{j}`\ th row and :math:`\textit{i}`\ th column of :math:`A` (the coefficient of the :math:`\textit{i}`\ th unknown in the :math:`\textit{j}`\ th equation), for :math:`\textit{j} = 1,2,\ldots,m`, for :math:`\textit{i} = 1,2,\ldots,n`.
The next :math:`l` columns contain the :math:`l` rows of the constraint matrix :math:`C`; element :math:`\mathrm{e}[\textit{i}-1,\textit{j}+m-1]` containing the element :math:`c_{{\textit{j}\textit{i}}}` in the :math:`\textit{j}`\ th row and :math:`\textit{i}`\ th column of :math:`C` (the coefficient of the :math:`\textit{i}`\ th unknown in the :math:`\textit{j}`\ th constraint), for :math:`\textit{j} = 1,2,\ldots,l`, for :math:`\textit{i} = 1,2,\ldots,n`.
**f** : float, array-like, shape :math:`\left(\textit{mpl}\right)`
:math:`\mathrm{f}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,m`, must contain :math:`b_{\textit{i}}` (the :math:`\textit{i}`\ th element of the right-hand side vector of the overdetermined system of equations) and :math:`\mathrm{f}[m+\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,l`, must contain :math:`d_i` (the :math:`i`\ th element of the right-hand side vector of the constraints), where :math:`l` is the number of constraints.
**x** : float, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{x}[\textit{i}-1]` must contain an estimate of the :math:`\textit{i}`\ th unknown, for :math:`\textit{i} = 1,2,\ldots,n`. If no better initial estimate for :math:`\mathrm{x}[i-1]` is available, set :math:`\mathrm{x}[i-1] = 0.0`.
**mxs** : int
The maximum number of steps to be allowed for the solution of the unconstrained problem. Typically this may be a modest multiple of :math:`n`. If, on entry, :math:`\mathrm{mxs}` is zero or negative, the value returned by :meth:`machine.integer_max <naginterfaces.library.machine.integer_max>` is used.
**monit** : callable monit(x, niter, k, el1n, data=None)
:math:`\mathrm{monit}` can be used to print out the current values of any selection of its arguments.
The frequency with which :math:`\mathrm{monit}` is called in ``glinc_l1sol`` is controlled by :math:`\mathrm{iprint}`.
**Parameters**
**x** : float, ndarray, shape :math:`\left(n\right)`
The latest estimate of the unknowns.
**niter** : int
The number of iterations so far carried out.
**k** : int
The total number of equations and constraints which are currently active (i.e., the number of equations with zero residuals plus the number of constraints which are satisfied as equations).
**el1n** : float
The :math:`l_1`-norm of the current residuals of the overdetermined system of equations.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**iprint** : int
The frequency of iteration print out.
:math:`\mathrm{iprint} > 0`
:math:`\mathrm{monit}` is called every :math:`\mathrm{iprint}` iterations and at the solution.
:math:`\mathrm{iprint} = 0`
Information is printed out at the solution only. Otherwise :math:`\mathrm{monit}` is not called.
**data** : arbitrary, optional
User-communication data for callback functions.
**Returns**
**e** : float, ndarray, shape :math:`\left(n, \textit{mpl}\right)`
Unchanged, except possibly to the extent of a small multiple of the machine precision. (See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02gbf.html#fcomments>`__.)
**x** : float, ndarray, shape :math:`\left(n\right)`
The latest estimate of the :math:`\textit{i}`\ th unknown, for :math:`\textit{i} = 1,2,\ldots,n`. If no exception or warning is raised on exit, these are the solution values.
**k** : int
The total number of equations and constraints which are then active (i.e., the number of equations with zero residuals plus the number of constraints which are satisfied as equalities).
**el1n** : float
The :math:`l_1`-norm (sum of absolute values) of the equation residuals.
**indx** : int, ndarray, shape :math:`\left(\textit{mpl}\right)`
Specifies which columns of :math:`\mathrm{e}` relate to the inactive equations and constraints. :math:`\mathrm{indx}[0]` up to :math:`\mathrm{indx}[\mathrm{k}-1]` number the active columns and :math:`\mathrm{indx}[\mathrm{k}]` up to :math:`\mathrm{indx}[\textit{mpl}-1]` number the inactive columns.
.. _e02gb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
The constraints cannot all be satisfied simultaneously: they are not compatible with one another. Hence no solution is possible.
(`errno` :math:`2`)
The limit imposed by :math:`\mathrm{mxs}` has been reached without finding a solution. Consider restarting from the current point by simply calling ``glinc_l1sol`` again without changing the arguments.
(`errno` :math:`3`)
The function has failed because of numerical difficulties; the problem is too ill-conditioned. Consider rescaling the unknowns.
(`errno` :math:`4`)
Elements :math:`1` to :math:`\mathrm{m}` of one of the first :math:`\textit{mpl}` columns of the array :math:`\mathrm{e}` are all zero -- this corresponds to a zero row in either of the matrices :math:`A` or :math:`C`.
(`errno` :math:`4`)
On entry, :math:`\textit{mpl} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{mpl}\geq \mathrm{m}`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}\geq n`.
(`errno` :math:`4`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n\geq 2`.
.. _e02gb-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
Given a matrix :math:`A` with :math:`m` rows and :math:`n` columns :math:`\left(m\geq n\right)` and a vector :math:`b` with :math:`m` elements, the function calculates an :math:`l_1` solution to the overdetermined system of equations
.. math::
Ax = b\text{.}
That is to say, it calculates a vector :math:`x`, with :math:`n` elements, which minimizes the :math:`l_1`-norm (the sum of the absolute values) of the residuals
.. math::
r\left(x\right) = \sum_{{i = 1}}^m\left\lvert r_i\right\rvert \text{,}
where the residuals :math:`r_i` are given by
.. math::
r_i = b_i-\sum_{{j = 1}}^na_{{ij}}x_j\text{, }\quad i = 1,2,\ldots,m\text{.}
Here :math:`a_{{ij}}` is the element in row :math:`i` and column :math:`j` of :math:`A`, :math:`b_i` is the :math:`i`\ th element of :math:`b` and :math:`x_j` the :math:`j`\ th element of :math:`x`.
If, in addition, a matrix :math:`C` with :math:`l` rows and :math:`n` columns and a vector :math:`d` with :math:`l` elements, are given, the vector :math:`x` computed by the function is such as to minimize the :math:`l_1`-norm :math:`r\left(x\right)` subject to the set of inequality constraints :math:`Cx\geq d`.
The matrices :math:`A` and :math:`C` need not be of full rank.
Typically in applications to data fitting, data consisting of :math:`m` points with coordinates :math:`\left(t_i, y_i\right)` is to be approximated by a linear combination of known functions :math:`\phi_i\left(t\right)`,
.. math::
\alpha_1\phi_1\left(t\right)+\alpha_2\phi_2\left(t\right)+ \cdots +\alpha_n\phi_n\left(t\right)\text{,}
in the :math:`l_1`-norm, possibly subject to linear inequality constraints on the coefficients :math:`\alpha_j` of the form :math:`C\alpha \geq d` where :math:`\alpha` is the vector of the :math:`\alpha_j` and :math:`C` and :math:`d` are as in the previous paragraph.
This is equivalent to finding an :math:`l_1` solution to the overdetermined system of equations
.. math::
\sum_{{j = 1}}^n\phi_j\left(t_i\right)\alpha_j = y_i\text{, }\quad i = 1,2,\ldots,m\text{,}
subject to :math:`C\alpha \geq d`.
Thus if, for each value of :math:`i` and :math:`j`, the element :math:`a_{{ij}}` of the matrix :math:`A` above is set equal to the value of :math:`\phi_j\left(t_i\right)` and :math:`b_i` is equal to :math:`y_i` and :math:`C` and :math:`d` are also supplied to the function, the solution vector :math:`x` will contain the required values of the :math:`\alpha_j`.
Note that the independent variable :math:`t` above can, instead, be a vector of several independent variables (this includes the case where each of :math:`\phi_i` is a function of a different variable, or set of variables).
The algorithm follows the Conn--Pietrzykowski approach (see Bartels `et al.` (1978) and Conn and Pietrzykowski (1977)), which is via an exact penalty function
.. math::
g\left(x\right) = \gamma r\left(x\right)-\sum_{{i = 1}}^l\mathrm{min}\left(0, {c_i^\mathrm{T}x-d_i}\right)\text{,}
where :math:`\gamma` is a penalty argument, :math:`c_i^\mathrm{T}` is the :math:`i`\ th row of the matrix :math:`C`, and :math:`d_i` is the :math:`i`\ th element of the vector :math:`d`.
It proceeds in a step-by-step manner much like the simplex method for linear programming but does not move from vertex to vertex and does not require the problem to be cast in a form containing only non-negative unknowns.
It uses stable procedures to update an orthogonal factorization of the current set of active equations and constraints.
.. _e02gb-py2-py-references:
**References**
Bartels, R H, Conn, A R and Charalambous, C, 1976, `Minimisation techniques for piecewise Differentiable functions -- the` :math:`l_{\infty }` `solution to an overdetermined linear system`, Technical Report No. 247, CORR 76/30, Mathematical Sciences Department, The John Hopkins University
Bartels, R H, Conn, A R and Sinclair, J W, 1976, `A Fortran program for solving overdetermined systems of linear equations in the` :math:`l_1` `Sense`, Technical Report No. 236, CORR 76/7, Mathematical Sciences Department, The John Hopkins University
Bartels, R H, Conn, A R and Sinclair, J W, 1978, `Minimisation techniques for piecewise differentiable functions -- the` :math:`l_1` `solution to an overdetermined linear system`, SIAM J. Numer. Anal. (15), 224--241
Conn, A R and Pietrzykowski, T, 1977, `A penalty-function method converging directly to a constrained optimum`, SIAM J. Numer. Anal. (14), 348--375
"""
raise NotImplementedError
[docs]def glin_linf(n, a, b, relerr, tol=0.0):
r"""
``glin_linf`` calculates an :math:`l_{\infty }` solution to an overdetermined system of linear equations.
.. _e02gc-py2-py-doc:
For full information please refer to the NAG Library document for e02gc
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02gcf.html
.. _e02gc-py2-py-parameters:
**Parameters**
**n** : int
The number of unknowns, :math:`n` (the number of columns of the matrix :math:`A`).
**a** : float, array-like, shape :math:`\left(\mathrm{n}+3, m+1\right)`
:math:`\mathrm{a}[\textit{j}-1,\textit{i}-1]` must contain :math:`a_{{\textit{i}\textit{j}}}`, the element in the :math:`\textit{i}`\ th row and :math:`\textit{j}`\ th column of the matrix :math:`A`, for :math:`\textit{j} = 1,2,\ldots,n`, for :math:`\textit{i} = 1,2,\ldots,m`, (that is, the **transpose** of the matrix). The remaining elements need not be set. Preferably, the columns of the matrix :math:`A` (rows of the argument :math:`\mathrm{a}`) should be scaled before entry: see `Accuracy <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02gcf.html#accuracy>`__.
**b** : float, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{b}[\textit{i}-1]` must contain :math:`b_{\textit{i}}`, the :math:`\textit{i}`\ th element of the vector :math:`b`, for :math:`\textit{i} = 1,2,\ldots,m`.
**relerr** : float
Must be set to a bound on the relative error acceptable in the maximum residual at the solution.
If :math:`\mathrm{relerr}\leq 0.0`, the :math:`l_{\infty }` solution is computed, and :math:`\mathrm{relerr}` is set to :math:`0.0` on exit.
If :math:`\mathrm{relerr} > 0.0`, the function obtains instead an approximate solution for which the largest residual is less than :math:`1.0+\mathrm{relerr}` times that of the :math:`l_{\infty }` solution; on exit, :math:`\mathrm{relerr}` contains a smaller value such that the above bound still applies. (The usual result of this option, say with :math:`\mathrm{relerr} = 0.1`, is a saving in the number of simplex iterations).
**tol** : float, optional
A threshold below which numbers are regarded as zero. The recommended threshold value is :math:`10.0\times \epsilon`, where :math:`\epsilon` is the machine precision. If :math:`\mathrm{tol}\leq 0.0` on entry, the recommended value is used within the function. If premature termination occurs, a larger value for :math:`\mathrm{tol}` may result in a valid solution.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{n}+3, m+1\right)`
Contains the last simplex tableau.
**b** : float, ndarray, shape :math:`\left(m\right)`
The :math:`\textit{i}`\ th residual :math:`\textit{r}_{\textit{i}}` corresponding to the solution vector :math:`x`, for :math:`\textit{i} = 1,2,\ldots,m`. Note however that these residuals may contain few significant figures, especially when :math:`\mathrm{resmax}` is within one or two orders of magnitude of :math:`\mathrm{tol}`. Indeed if :math:`\mathrm{resmax}\leq \mathrm{tol}`, the elements :math:`\mathrm{b}[i-1]` may all be set to zero. It is, therefore, often advisable to compute the residuals directly.
**relerr** : float
Is altered as described above.
**x** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
If the function exits successfully or :math:`\mathrm{errno}` = 1, :math:`\mathrm{x}[\textit{j}-1]` contains the :math:`\textit{j}`\ th element of the solution vector :math:`x`, for :math:`\textit{j} = 1,2,\ldots,n`. Whether this is an :math:`l_{\infty }` solution or an approximation to one, depends on the value of :math:`\mathrm{relerr}` on entry.
**resmax** : float
If the function exits successfully or :math:`\mathrm{errno}` = 1, :math:`\mathrm{resmax}` contains the absolute value of the largest residual(s) for the solution vector :math:`x`. (See :math:`\mathrm{b}`.)
**irank** : int
If the function exits successfully or :math:`\mathrm{errno}` = 1, :math:`\mathrm{irank}` contains the computed rank of the matrix :math:`A`.
**itera** : int
If the function exits successfully or :math:`\mathrm{errno}` = 1, :math:`\mathrm{itera}` contains the number of iterations taken by the simplex method.
.. _e02gc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`2`)
Premature termination due to rounding errors. Try using larger value of :math:`\mathrm{tol}`: :math:`\mathrm{tol} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 1`.
(`errno` :math:`3`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq \mathrm{n}`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
An optimal solution has been obtained, but may not be unique.
.. _e02gc-py2-py-notes:
**Notes**
Given a matrix :math:`A` with :math:`m` rows and :math:`n` columns :math:`\left(m\geq n\right)` and a vector :math:`b` with :math:`m` elements, the function calculates an :math:`l_{\infty }` solution to the overdetermined system of equations
.. math::
Ax = b\text{.}
That is to say, it calculates a vector :math:`x`, with :math:`n` elements, which minimizes the :math:`l_{\infty }` norm of the residuals (the absolutely largest residual)
.. math::
r\left(x\right) = \mathrm{max}_{{1\leq i\leq m}}\left\lvert r_i\right\rvert
where the residuals :math:`r_i` are given by
.. math::
r_i = b_i-\sum_{{j = 1}}^na_{{ij}}x_{{j}}\text{, }\quad i = 1,2,\ldots,m\text{.}
Here :math:`a_{{ij}}` is the element in row :math:`i` and column :math:`j` of :math:`A`, :math:`b_i` is the :math:`i`\ th element of :math:`b` and :math:`x_j` the :math:`j`\ th element of :math:`x`.
The matrix :math:`A` need not be of full rank.
The solution is not unique in this case, and may not be unique even if :math:`A` is of full rank.
Alternatively, in applications where a complete minimization of the :math:`l_{\infty }` norm is not necessary, you may obtain an approximate solution, usually in shorter time, by giving an appropriate value to the argument :math:`\mathrm{relerr}`.
Typically in applications to data fitting, data consisting of :math:`m` points with coordinates :math:`\left(t_i, y_i\right)` is to be approximated in the :math:`l_{\infty }` norm by a linear combination of known functions :math:`\phi_j\left(t\right)`,
.. math::
\alpha_1\phi_1\left(t\right)+\alpha_2\phi_2\left(t\right)+ \cdots +\alpha_n\phi_n\left(t\right)\text{.}
This is equivalent to finding an :math:`l_{\infty }` solution to the overdetermined system of equations
.. math::
\sum_{{j = 1}}^n\phi_j\left(t_i\right)\alpha_j = y_{{i}}\text{, }\quad i = 1,2,\ldots,m\text{.}
Thus if, for each value of :math:`i` and :math:`j` the element :math:`a_{{ij}}` of the matrix :math:`A` above is set equal to the value of :math:`\phi_j\left(t_i\right)` and :math:`b_i` is set equal to :math:`y_i`, the solution vector :math:`x` will contain the required values of the :math:`\alpha_j`.
Note that the independent variable :math:`t` above can, instead, be a vector of several independent variables (this includes the case where each :math:`\phi_i` is a function of a different variable, or set of variables).
The algorithm is a modification of the simplex method of linear programming applied to the dual formation of the :math:`l_{\infty }` problem (see Barrodale and Phillips (1974) and Barrodale and Phillips (1975)).
The modifications are designed to improve the efficiency and stability of the simplex method for this particular application.
.. _e02gc-py2-py-references:
**References**
Barrodale, I and Phillips, C, 1974, `An improved algorithm for discrete Chebyshev linear approximation`, Proc. 4th Manitoba Conf. Numerical Mathematics, 177--190, University of Manitoba, Canada
Barrodale, I and Phillips, C, 1975, `Solution of an overdetermined system of linear equations in the Chebyshev norm [F4] (Algorithm 495)`, ACM Trans. Math. Software (1(3)), 264--270
"""
raise NotImplementedError
[docs]def dim2_spline_ts_sctr(x, y, f, lsminp, lsmaxp, nxcels, nycels, comm):
r"""
``dim2_spline_ts_sctr`` computes a spline approximation to a set of scattered data using a two-stage approximation method.
The computational complexity of the method grows linearly with the number of data points; hence large datasets are easily accommodated.
Note: this function uses optional algorithmic parameters, see also: :meth:`opt_set`, :meth:`opt_get`.
.. _e02jd-py2-py-doc:
For full information please refer to the NAG Library document for e02jd
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02jdf.html
.. _e02jd-py2-py-parameters:
**Parameters**
**x** : float, array-like, shape :math:`\left(n\right)`
The :math:`\left(x_i, y_i, f_i\right)` data values to be fitted.
**y** : float, array-like, shape :math:`\left(n\right)`
The :math:`\left(x_i, y_i, f_i\right)` data values to be fitted.
**f** : float, array-like, shape :math:`\left(n\right)`
The :math:`\left(x_i, y_i, f_i\right)` data values to be fitted.
**lsminp** : int
:math:`\mathrm{lsminp}` and :math:`\mathrm{lsmaxp}` are control parameters for the local approximations.
Each local approximation is computed on a local domain containing one of the triangles in the discretization of the bounding box.
The size of each local domain will be adaptively chosen such that if it contains fewer than :math:`\mathrm{lsminp}` sample points it is expanded, else if it contains greater than :math:`\mathrm{lsmaxp}` sample points a thinning method is applied. :math:`\mathrm{lsmaxp}` mainly controls computational cost (in that working with a thinned set of points is cheaper and may be appropriate if the input data is densely distributed), while :math:`\mathrm{lsminp}` allows handling of different types of scattered data.
Setting :math:`\mathrm{lsmaxp} < \mathrm{lsminp}`, and, therefore, forcing either expansion or thinning, may be useful for computing initial coarse approximations.
In general smaller values for these arguments reduces cost.
A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate values for :math:`\mathrm{lsminp}` and :math:`\mathrm{lsmaxp}`.
**lsmaxp** : int
:math:`\mathrm{lsminp}` and :math:`\mathrm{lsmaxp}` are control parameters for the local approximations.
Each local approximation is computed on a local domain containing one of the triangles in the discretization of the bounding box.
The size of each local domain will be adaptively chosen such that if it contains fewer than :math:`\mathrm{lsminp}` sample points it is expanded, else if it contains greater than :math:`\mathrm{lsmaxp}` sample points a thinning method is applied. :math:`\mathrm{lsmaxp}` mainly controls computational cost (in that working with a thinned set of points is cheaper and may be appropriate if the input data is densely distributed), while :math:`\mathrm{lsminp}` allows handling of different types of scattered data.
Setting :math:`\mathrm{lsmaxp} < \mathrm{lsminp}`, and, therefore, forcing either expansion or thinning, may be useful for computing initial coarse approximations.
In general smaller values for these arguments reduces cost.
A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate values for :math:`\mathrm{lsminp}` and :math:`\mathrm{lsmaxp}`.
**nxcels** : int
:math:`\mathrm{nxcels}` (respectively :math:`\mathrm{nycels}`) is the number of cells in the :math:`x` (respectively :math:`y`) direction that will be used to create the triangulation of the bounding box of the domain of the function to be fitted.
Greater efficiency generally comes when :math:`\mathrm{nxcels}` and :math:`\mathrm{nycels}` are chosen to be of the same order of magnitude and are such that :math:`\textit{n}` is :math:`\mathrm{O}\left(\mathrm{nxcels}\times \mathrm{nycels}\right)`.
Thus for a 'square' triangulation --- when :math:`\mathrm{nxcels} = \mathrm{nycels}` --- the quantities :math:`\sqrt{n}` and :math:`\mathrm{nxcels}` should be of the same order of magnitude.
See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02jdf.html#fcomments>`__.
**nycels** : int
:math:`\mathrm{nxcels}` (respectively :math:`\mathrm{nycels}`) is the number of cells in the :math:`x` (respectively :math:`y`) direction that will be used to create the triangulation of the bounding box of the domain of the function to be fitted.
Greater efficiency generally comes when :math:`\mathrm{nxcels}` and :math:`\mathrm{nycels}` are chosen to be of the same order of magnitude and are such that :math:`\textit{n}` is :math:`\mathrm{O}\left(\mathrm{nxcels}\times \mathrm{nycels}\right)`.
Thus for a 'square' triangulation --- when :math:`\mathrm{nxcels} = \mathrm{nycels}` --- the quantities :math:`\sqrt{n}` and :math:`\mathrm{nxcels}` should be of the same order of magnitude.
See also `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02jdf.html#fcomments>`__.
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt_set`.
.. _e02jd-py2-py-other_params:
**Other Parameters**
**'Averaged Spline'** : str
Default :math:`\text{} = \texttt{'NO'}`
When the bounding box is triangulated there are :math:`8` equivalent configurations of the mesh.
Setting :math:`\text{‘Averaged Spline'} = \texttt{'YES'}` will use the averaged value of the :math:`8` possible local polynomial approximations over each triangle in the mesh.
This usually gives better results but at (about 8 times) higher computational cost.
Constraint: :math:`\text{‘Averaged Spline'} = \texttt{'YES'}` or :math:`\texttt{'NO'}`.
**'Global Smoothing Level'** : int
Default :math:`\text{} = 1`
The smoothness level for the global spline approximation.
:math:`\text{‘Global Smoothing Level'} = 1`
Will use :math:`C^1` piecewise cubics.
:math:`\text{‘Global Smoothing Level'} = 2`
Will use :math:`C^2` piecewise sextics.
**'Interpolation Only RBF'** : str
Default :math:`\text{} = \texttt{'YES'}`
If :math:`\text{‘Interpolation Only RBF'} = \texttt{'YES'}`, each local RBF approximation is computed by interpolation.
If :math:`\text{‘Interpolation Only RBF'} = \texttt{'NO'}`, each local RBF approximation is computed by a discrete least squares approach.
This is likely to be more accurate and more expensive than interpolation.
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` or :math:`\texttt{'POLYNOMIAL'}`, this option setting is ignored.
**'Local Method'** : str
Default :math:`\text{} = \texttt{'POLYNOMIAL'}`
The local approximation scheme to use.
:math:`\text{‘Local Method'} = \texttt{'POLYNOMIAL'}`
Uses least squares polynomial approximations.
:math:`\text{‘Local Method'} = \texttt{'HYBRID'}`
Uses hybrid polynomial and RBF approximations.
:math:`\text{‘Local Method'} = \texttt{'RBF'}`
Uses pure RBF approximations.
In general :math:`\texttt{'POLYNOMIAL'}` is less computationally expensive than :math:`\texttt{'HYBRID'}` is less computationally expensive than :math:`\texttt{'RBF'}` with the reverse ordering holding for accuracy of results.
**'Minimum Singular Value LHA'** : float
Default :math:`\text{} = 1.0`
A tolerance measure for accepting or rejecting a local hybrid approximation (LHA) as reliable.
The solution of a local least squares problem solved on each triangle subdomain is accepted as reliable if the minimum singular value :math:`\sigma` of the collocation matrix (of polynomial and radial basis function terms) associated with the least squares problem satisfies :math:`\text{‘Minimum Singular Value LHA'}\leq \sigma`.
In general the approximation power will be reduced as 'Minimum Singular Value LHA' is reduced. (A small :math:`\sigma` indicates that the local data has hidden redundancies which prevent it from carrying enough information for a good approximation to be made.) Setting 'Minimum Singular Value LHA' very large may have the detrimental effect that only approximations of low degree are deemed reliable.
A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate value for this parameter.
If :math:`\text{‘Local Method'} = \texttt{'POLYNOMIAL'}` or :math:`\texttt{'RBF'}`, this option setting is ignored.
**'Minimum Singular Value LPA'** : float
Default :math:`\text{} = 1.0`
A tolerance measure for accepting or rejecting a local polynomial approximation (LPA) as reliable.
Clearly this setting is relevant when :math:`\text{‘Local Method'} = \texttt{'POLYNOMIAL'}`, but it also may be used when :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02jdf.html#fcomments>`__.)
The solution of a local least squares problem solved on each triangle subdomain is accepted as reliable if the minimum singular value :math:`\sigma` of the matrix (of Bernstein polynomial values) associated with the least squares problem satisfies :math:`\text{‘Minimum Singular Value LPA'}\leq \sigma`.
In general the approximation power will be reduced as 'Minimum Singular Value LPA' is reduced. (A small :math:`\sigma` indicates that the local data has hidden redundancies which prevent it from carrying enough information for a good approximation to be made.) Setting 'Minimum Singular Value LPA' very large may have the detrimental effect that only approximations of low degree are deemed reliable.
'Minimum Singular Value LPA' will have no effect if :math:`\text{‘Polynomial Starting Degree'} = 0`, and it will have little effect if the input data is 'smooth' (e.g., from a known function).
A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate value for this parameter.
If :math:`\text{‘Local Method'} = \texttt{'RBF'}`, this option setting is ignored.
**'Polynomial Starting Degree'** : int
Default :math:`\text{} = 5` if :math:`\text{‘Local Method'} = \texttt{'HYBRID'}`,
Default :math:`\text{} = 1` otherwise
The degree to be used for the polynomial part in the initial step of each local approximation.
At this initial step the method will attempt to fit with a local approximation having polynomial part of degree 'Polynomial Starting Degree'.
If :math:`\text{‘Local Method'} = \texttt{'POLYNOMIAL'}` and the approximation is deemed unreliable (according to 'Minimum Singular Value LPA'), the degree will be decremented by one and a new local approximation computed, ending with a constant approximation if no other is reliable.
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` and the approximation is deemed unreliable (according to 'Minimum Singular Value LHA'), a pure polynomial approximation of this degree will be tried instead.
The method then proceeds as in the 'POLYNOMIAL' case.
'Polynomial Starting Degree' is bounded from above by the maximum possible spline degree, which is :math:`6` (when performing :math:`C^2` global super-smoothing).
Note that the best-case approximation error (see `Accuracy <https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02jdf.html#accuracy>`__) for :math:`C^2` smoothing with :math:`\text{‘Supersmooth C2'} = \texttt{'NO'}` is achieved for local polynomials of degree :math:`5`; that is, for this level of global smoothing no further benefit is gained by setting :math:`\text{‘Polynomial Starting Degree'} = 6`.
The default value gives a good compromise between efficiency and accuracy.
In general the best approximation can be obtained by setting:
If :math:`\text{‘Local Method'} = \texttt{'POLYNOMIAL'}`
if :math:`\text{‘Global Smoothing Level'} = 1`, :math:`\text{‘Polynomial Starting Degree'} = 3`;
if :math:`\text{‘Global Smoothing Level'} = 2`;
if :math:`\text{‘Supersmooth C2'} = \texttt{'NO'}`, :math:`\text{‘Polynomial Starting Degree'} = 5`;
otherwise :math:`\text{‘Polynomial Starting Degree'} = 6`.
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}`, 'Polynomial Starting Degree' as small as possible.
If :math:`\text{‘Local Method'} = \texttt{'RBF'}`, this option setting is ignored.
**'Radial Basis Function'** : str
Default :math:`\text{} = \texttt{'MQ'}`
'Radial Basis Function' selects the RBF to use in each local RBF approximation, while 'Scaling Coefficient RBF' selects the scale factor to use in its evaluation, as described below.
A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate scale factor and RBF.
If :math:`\text{‘Local Method'} = \texttt{'POLYNOMIAL'}`, these option settings are ignored.
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` or :math:`\texttt{'RBF'}`, the following (conditionally) positive definite functions may be chosen.
Define :math:`R = \sqrt{x^2+y^2}` and :math:`\rho = R/r`.
.. rst-class:: nag-rules-none nag-align-left
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'GAUSS'}` |Gaussian :math:`\mathrm{exp}\left(-\rho^2\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ'}` |inverse multiquadric :math:`1/\sqrt{r^2+R^2}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ2'}` |inverse multiquadric :math:`1/\left(r^2+R^2\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ3'}` |inverse multiquadric :math:`1/\left(r^2+R^2\right)^{\left(3/2\right)}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ0_5'}` |inverse multiquadric :math:`1/\left(r^2+R^2\right)^{\left(1/4\right)}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'WENDLAND31'}` |H. Wendland's :math:`C^2` function :math:`\mathrm{max}\left(0, {1-\rho }\right)^4\left(4\rho +1\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'WENDLAND32'}` |H. Wendland's :math:`C^4` function :math:`\mathrm{max}\left(0, {1-\rho }\right)^6\left(35\rho^2+18\rho +3\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'WENDLAND33'}` |H. Wendland's :math:`C^6` function :math:`\mathrm{max}\left(0, {1-\rho }\right)^8\left(32\rho^3+25\rho^2+8\rho +1\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'BUHMANNC3'}` |M. Buhmann's :math:`C^3` function :math:`112/45\rho^{\left(9/2\right)}+16/3\rho^{\left(7/2\right)}-7\rho^4-14/15\rho^2+1/9` if :math:`\rho \leq 1`, :math:`0` otherwise|
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'MQ'}` |multiquadric :math:`\sqrt{r^2+R^2}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'MQ1_5'}` |multiquadric :math:`\left(r^2+R^2\right)^{\left(1.5/2\right)}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC1_5'}` |polyharmonic spline :math:`\rho^{1.5}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC1_75'}`|polyharmonic spline :math:`\rho^{1.75}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` the following conditionally positive definite functions may also be chosen.
.. rst-class:: nag-rules-none nag-align-left
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'MQ2'}` |multiquadric :math:`\left(r^2+R^2\right)\log\left(r^2+R^2\right)`|
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'MQ3'}` |multiquadric :math:`\left(r^2+R^2\right)^{\left(3/2\right)}` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'TPS'}` |thin plate spline :math:`\rho^2\log\left(\rho^2\right)` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC3'}`|polyharmonic spline :math:`\rho^3` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'TPS4'}` |thin plate spline :math:`\rho^4\log\left(\rho^2\right)` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC5'}`|polyharmonic spline :math:`\rho^5` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'TPS6'}` |thin plate spline :math:`\rho^6\log\left(\rho^2\right)` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC7'}`|polyharmonic spline :math:`\rho^7` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC9'}`|polyharmonic spline :math:`\rho^9` |
+--------------------------------+-----------------------------------------------------------------+
**'Scaling Coefficient RBF'** : float
Default :math:`\text{} = 1.0`
'Radial Basis Function' selects the RBF to use in each local RBF approximation, while 'Scaling Coefficient RBF' selects the scale factor to use in its evaluation, as described below.
A calibration procedure (experimenting with a small subset of the data to be fitted and validating the results) may be needed to choose the most appropriate scale factor and RBF.
If :math:`\text{‘Local Method'} = \texttt{'POLYNOMIAL'}`, these option settings are ignored.
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` or :math:`\texttt{'RBF'}`, the following (conditionally) positive definite functions may be chosen.
Define :math:`R = \sqrt{x^2+y^2}` and :math:`\rho = R/r`.
.. rst-class:: nag-rules-none nag-align-left
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'GAUSS'}` |Gaussian :math:`\mathrm{exp}\left(-\rho^2\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ'}` |inverse multiquadric :math:`1/\sqrt{r^2+R^2}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ2'}` |inverse multiquadric :math:`1/\left(r^2+R^2\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ3'}` |inverse multiquadric :math:`1/\left(r^2+R^2\right)^{\left(3/2\right)}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'IMQ0_5'}` |inverse multiquadric :math:`1/\left(r^2+R^2\right)^{\left(1/4\right)}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'WENDLAND31'}` |H. Wendland's :math:`C^2` function :math:`\mathrm{max}\left(0, {1-\rho }\right)^4\left(4\rho +1\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'WENDLAND32'}` |H. Wendland's :math:`C^4` function :math:`\mathrm{max}\left(0, {1-\rho }\right)^6\left(35\rho^2+18\rho +3\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'WENDLAND33'}` |H. Wendland's :math:`C^6` function :math:`\mathrm{max}\left(0, {1-\rho }\right)^8\left(32\rho^3+25\rho^2+8\rho +1\right)` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'BUHMANNC3'}` |M. Buhmann's :math:`C^3` function :math:`112/45\rho^{\left(9/2\right)}+16/3\rho^{\left(7/2\right)}-7\rho^4-14/15\rho^2+1/9` if :math:`\rho \leq 1`, :math:`0` otherwise|
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'MQ'}` |multiquadric :math:`\sqrt{r^2+R^2}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'MQ1_5'}` |multiquadric :math:`\left(r^2+R^2\right)^{\left(1.5/2\right)}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC1_5'}` |polyharmonic spline :math:`\rho^{1.5}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC1_75'}`|polyharmonic spline :math:`\rho^{1.75}` |
+-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` the following conditionally positive definite functions may also be chosen.
.. rst-class:: nag-rules-none nag-align-left
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'MQ2'}` |multiquadric :math:`\left(r^2+R^2\right)\log\left(r^2+R^2\right)`|
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'MQ3'}` |multiquadric :math:`\left(r^2+R^2\right)^{\left(3/2\right)}` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'TPS'}` |thin plate spline :math:`\rho^2\log\left(\rho^2\right)` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC3'}`|polyharmonic spline :math:`\rho^3` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'TPS4'}` |thin plate spline :math:`\rho^4\log\left(\rho^2\right)` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC5'}`|polyharmonic spline :math:`\rho^5` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'TPS6'}` |thin plate spline :math:`\rho^6\log\left(\rho^2\right)` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC7'}`|polyharmonic spline :math:`\rho^7` |
+--------------------------------+-----------------------------------------------------------------+
|:math:`\texttt{'POLYHARMONIC9'}`|polyharmonic spline :math:`\rho^9` |
+--------------------------------+-----------------------------------------------------------------+
**'Separation LRBFA'** : float
Default :math:`\text{} = 16.0/\text{‘Scaling Coefficient RBF'}`
A knot-separation parameter used to control the condition number of the matrix used in each local RBF approximation (LRBFA).
A smaller value may mean greater numerical stability but fewer knots.
If :math:`\text{‘Local Method'} = \texttt{'HYBRID'}` or :math:`\texttt{'POLYNOMIAL'}`, this option setting is ignored.
**'Supersmooth C2'** : str
Default :math:`\text{} = \texttt{'NO'}`
If :math:`\text{‘Supersmooth C2'} = \texttt{'YES'}`, the :math:`C^2` spline is generated using additional smoothness constraints.
This usually gives better results but at higher computational cost.
If :math:`\text{‘Global Smoothing Level'} = 1` this option setting is ignored.
.. _e02jd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`2`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n > 1`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{lsminp} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{lsminp}\leq n`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{lsmaxp} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lsmaxp}\geq 1`.
(`errno` :math:`6`)
On entry, :math:`\mathrm{nxcels} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nxcels}\geq 1`.
(`errno` :math:`7`)
On entry, :math:`\mathrm{nycels} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nycels}\geq 1`.
(`errno` :math:`9`)
Option arrays are not initialized or are corrupted.
(`errno` :math:`11`)
An unexpected algorithmic failure was encountered. Please contact `NAG <https://www.nag.com>`__.
(`errno` :math:`12`)
On entry, all elements of :math:`\mathrm{x}` or of :math:`\mathrm{y}` are equal.
(`errno` :math:`20`)
The selected radial basis function cannot be used with the RBF local method.
(`errno` :math:`21`)
The value of option 'Polynomial Starting Degree' was invalid.
.. _e02jd-py2-py-notes:
**Notes**
``dim2_spline_ts_sctr`` determines a smooth bivariate spline approximation to a set of data points :math:`\left(x_{\textit{i}}, y_{\textit{i}}, f_{\textit{i}}\right)`, for :math:`\textit{i} = 1,2,\ldots,n`.
Here, 'smooth' means :math:`C^1` or :math:`C^2`. (You may select the degree of smoothing using the option 'Global Smoothing Level'.)
The approximation domain is the bounding box :math:`\left[x_{\mathrm{min}}, x_{\mathrm{max}}\right]\times \left[y_{\mathrm{min}}, y_{\mathrm{max}}\right]`, where :math:`x_{\mathrm{min}}` (respectively :math:`y_{\mathrm{min}}`) and :math:`x_{\mathrm{max}}` (respectively :math:`y_{\mathrm{max}}`) denote the lowest and highest data values of the :math:`\left(x_i\right)` (respectively :math:`\left(y_i\right)`).
The spline is computed by local approximations on a uniform triangulation of the bounding box.
These approximations are extended to a smooth spline representation of the surface over the domain.
The local approximation scheme is controlled by the option 'Local Method'.
The schemes provided are: by least squares polynomial approximation (Davydov and Zeilfelder (2004)); by hybrid polynomial and radial basis function (RBF) approximation (Davydov `et al.` (2006)); or by pure RBF approximation (Davydov `et al.` (2005)).
The two-stage approximation method employed by ``dim2_spline_ts_sctr`` is derived from the TSFIT package of O. Davydov and F. Zeilfelder.
Values of the computed spline can subsequently be computed by calling :meth:`dim2_spline_ts_evalv` or :meth:`dim2_spline_ts_evalm`.
.. _e02jd-py2-py-references:
**References**
Davydov, O, Morandi, R and Sestini, A, 2006, `Local hybrid approximation for scattered data fitting with bivariate splines`, Comput. Aided Geom. Design (23), 703--721
Davydov, O, Sestini, A and Morandi, R, 2005, `Local RBF approximation for scattered data fitting with bivariate splines`, `Trends and Applications in Constructive Approximation`, M. G. de Bruin, D. H. Mache, and J. Szabados, Eds (ISNM Vol. 151), Birkhauser, 91--102
Davydov, O and Zeilfelder, F, 2004, `Scattered data fitting by direct extension of local polynomials to bivariate splines`, Advances in Comp. Math. (21), 223--271
See Also
--------
:meth:`naginterfaces.library.examples.fit.dim2_spline_ts_sctr_ex.main`
"""
raise NotImplementedError
[docs]def dim2_spline_ts_evalv(xevalv, yevalv, comm):
r"""
``dim2_spline_ts_evalv`` calculates a vector of values of a spline computed by :meth:`dim2_spline_ts_sctr`.
.. _e02je-py2-py-doc:
For full information please refer to the NAG Library document for e02je
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02jef.html
.. _e02je-py2-py-parameters:
**Parameters**
**xevalv** : float, array-like, shape :math:`\left(\textit{nevalv}\right)`
The :math:`\left(x_i\right)` values at which the spline is to be evaluated.
**yevalv** : float, array-like, shape :math:`\left(\textit{nevalv}\right)`
The :math:`\left(y_i\right)` values at which the spline is to be evaluated.
**comm** : dict, communication object
Communication structure.
This argument must have been initialized by prior calls to :meth:`dim2_spline_ts_sctr` and :meth:`opt_set`.
**Returns**
**fevalv** : float, ndarray, shape :math:`\left(\textit{nevalv}\right)`
If no exception or warning is raised on exit :math:`\mathrm{fevalv}[i-1]` contains the computed spline value at :math:`\left(x_i, y_i\right)`.
.. _e02je-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`2`)
On entry, :math:`\textit{nevalv} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nevalv}\geq 1`.
(`errno` :math:`9`)
Option arrays are not initialized or are corrupted.
(`errno` :math:`10`)
The fitting routine has not been called, or the array of coefficients has been corrupted.
(`errno` :math:`13`)
On entry, :math:`\mathrm{xevalv}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle` was outside the bounding box.
Constraint: :math:`\langle\mathit{\boldsymbol{value}}\rangle\leq \mathrm{xevalv}[i-1]\leq \langle\mathit{\boldsymbol{value}}\rangle` for all :math:`i`.
(`errno` :math:`14`)
On entry, :math:`\mathrm{yevalv}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle` was outside the bounding box.
Constraint: :math:`\langle\mathit{\boldsymbol{value}}\rangle\leq \mathrm{yevalv}[i-1]\leq \langle\mathit{\boldsymbol{value}}\rangle` for all :math:`i`.
.. _e02je-py2-py-notes:
**Notes**
``dim2_spline_ts_evalv`` calculates values at prescribed points :math:`\left(x_i, y_i\right)`, for :math:`\textit{i} = 1,2,\ldots,n`, of a bivariate spline computed by :meth:`dim2_spline_ts_sctr`.
It is derived from the TSFIT package of O. Davydov and F. Zeilfelder.
.. _e02je-py2-py-references:
**References**
Davydov, O, Morandi, R and Sestini, A, 2006, `Local hybrid approximation for scattered data fitting with bivariate splines`, Comput. Aided Geom. Design (23), 703--721
Davydov, O, Sestini, A and Morandi, R, 2005, `Local RBF approximation for scattered data fitting with bivariate splines`, `Trends and Applications in Constructive Approximation`, M. G. de Bruin, D. H. Mache, and J. Szabados, Eds (ISNM Vol. 151), Birkhauser, 91--102
Davydov, O and Zeilfelder, F, 2004, `Scattered data fitting by direct extension of local polynomials to bivariate splines`, Advances in Comp. Math. (21), 223--271
Farin, G and Hansford, D, 2000, `The Essentials of CAGD`, Natic, MA: A K Peters, Ltd.
See Also
--------
:meth:`naginterfaces.library.examples.fit.dim2_spline_ts_sctr_ex.main`
"""
raise NotImplementedError
[docs]def dim2_spline_ts_evalm(xevalm, yevalm, comm):
r"""
``dim2_spline_ts_evalm`` calculates a mesh of values of a spline computed by :meth:`dim2_spline_ts_sctr`.
.. _e02jf-py2-py-doc:
For full information please refer to the NAG Library document for e02jf
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02jff.html
.. _e02jf-py2-py-parameters:
**Parameters**
**xevalm** : float, array-like, shape :math:`\left(\textit{nxeval}\right)`
The :math:`\left(x_i\right)` values forming the mesh on which the spline is to be evaluated.
**yevalm** : float, array-like, shape :math:`\left(\textit{nyeval}\right)`
The :math:`\left(y_j\right)` values forming the mesh on which the spline is to be evaluated.
**comm** : dict, communication object
Communication structure.
This argument must have been initialized by prior calls to :meth:`dim2_spline_ts_sctr` and :meth:`opt_set`.
**Returns**
**fevalm** : float, ndarray, shape :math:`\left(\textit{nxeval}, \textit{nyeval}\right)`
If no exception or warning is raised on exit :math:`\mathrm{fevalm}[i-1,j-1]` contains the computed spline value at :math:`\left(x_i, y_j\right)`.
.. _e02jf-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`6`)
On entry, :math:`\textit{nxeval} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nxeval}\geq 1`.
(`errno` :math:`7`)
On entry, :math:`\textit{nyeval} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nyeval}\geq 1`.
(`errno` :math:`9`)
Option arrays are not initialized or are corrupted.
(`errno` :math:`10`)
The fitting routine has not been called, or the array of coefficients has been corrupted.
(`errno` :math:`13`)
On entry, :math:`\mathrm{xevalm}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle` was outside the bounding box.
Constraint: :math:`\langle\mathit{\boldsymbol{value}}\rangle\leq \mathrm{xevalm}[i-1]\leq \langle\mathit{\boldsymbol{value}}\rangle` for all :math:`i`.
(`errno` :math:`14`)
On entry, :math:`\mathrm{yevalm}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle` was outside the bounding box.
Constraint: :math:`\langle\mathit{\boldsymbol{value}}\rangle\leq \mathrm{yevalm}[j-1]\leq \langle\mathit{\boldsymbol{value}}\rangle` for all :math:`j`.
.. _e02jf-py2-py-notes:
**Notes**
``dim2_spline_ts_evalm`` calculates values on a rectangular mesh of a bivariate spline computed by :meth:`dim2_spline_ts_sctr`.
The points in the mesh are defined by :math:`x` coordinates (:math:`x_{\textit{i}}`), for :math:`\textit{i} = 1,2,\ldots,n_x`, and :math:`y` coordinates (:math:`y_{\textit{j}}`), for :math:`\textit{j} = 1,2,\ldots,n_y`.
This function is derived from the TSFIT package of O. Davydov and F. Zeilfelder.
.. _e02jf-py2-py-references:
**References**
Davydov, O, Morandi, R and Sestini, A, 2006, `Local hybrid approximation for scattered data fitting with bivariate splines`, Comput. Aided Geom. Design (23), 703--721
Davydov, O, Sestini, A and Morandi, R, 2005, `Local RBF approximation for scattered data fitting with bivariate splines`, `Trends and Applications in Constructive Approximation`, M. G. de Bruin, D. H. Mache, and J. Szabados, Eds (ISNM Vol. 151), Birkhauser, 91--102
Davydov, O and Zeilfelder, F, 2004, `Scattered data fitting by direct extension of local polynomials to bivariate splines`, Advances in Comp. Math. (21), 223--271
Farin, G and Hansford, D, 2000, `The Essentials of CAGD`, Natic, MA: A K Peters, Ltd.
"""
raise NotImplementedError
[docs]def pade_app(ia, ib, c):
r"""
``pade_app`` calculates the coefficients in a Padé approximant to a function from its user-supplied Maclaurin expansion.
.. _e02ra-py2-py-doc:
For full information please refer to the NAG Library document for e02ra
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02raf.html
.. _e02ra-py2-py-parameters:
**Parameters**
**ia** : int
:math:`\mathrm{ia}` must specify :math:`l+1` and :math:`\mathrm{ib}` must specify :math:`m+1`, where :math:`l` and :math:`m` are the degrees of the numerator and denominator of the approximant, respectively.
**ib** : int
:math:`\mathrm{ia}` must specify :math:`l+1` and :math:`\mathrm{ib}` must specify :math:`m+1`, where :math:`l` and :math:`m` are the degrees of the numerator and denominator of the approximant, respectively.
**c** : float, array-like, shape :math:`\left(\mathrm{ia}+\mathrm{ib}-1\right)`
:math:`\mathrm{c}[\textit{i}-1]` must specify, for :math:`\textit{i} = 1,2,\ldots,l+m+1`, the coefficient of :math:`x^{{\textit{i}-1}}` in the given power series.
**Returns**
**a** : float, ndarray, shape :math:`\left(\mathrm{ia}\right)`
:math:`\mathrm{a}[\textit{j}]`, for :math:`\textit{j} = 1,2,\ldots,l+1`, contains the coefficient :math:`a_{\textit{j}}` in the numerator of the approximant.
**b** : float, ndarray, shape :math:`\left(\mathrm{ib}\right)`
:math:`\mathrm{b}[\textit{k}]`, for :math:`\textit{k} = 1,2,\ldots,m+1`, contains the coefficient :math:`b_{\textit{k}}` in the denominator of the approximant.
.. _e02ra-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{ic} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{ia} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ib} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ic}\geq \mathrm{ia}+\mathrm{ib}-1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ib} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ib}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ia} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ia}\geq 1`.
(`errno` :math:`2`)
The Pade approximant is degenerate.
.. _e02ra-py2-py-notes:
**Notes**
Given a power series
.. math::
c_0+c_1x+c_2x^2 + \cdots +c_{{l+m}}x^{{l+m}} + \cdots
``pade_app`` uses the coefficients :math:`c_i`, for :math:`\textit{i} = 0,1,\ldots,l+m`, to form the :math:`\left[l/m\right]` Padé approximant of the form
.. math::
\frac{{a_0+a_1x+a_2x^2 + \cdots +a_lx^l}}{{b_0+b_1x+b_2x^2 + \cdots +b_mx^m}}
with :math:`b_0` defined to be unity.
The two sets of coefficients :math:`a_j`, for :math:`\textit{j} = 0,1,\ldots,l`, and :math:`b_k`, for :math:`\textit{k} = 0,1,\ldots,m`, in the numerator and denominator are calculated by direct solution of the Padé equations (see Graves--Morris (1979)); these values are returned through the argument list unless the approximant is degenerate.
Padé approximation is a useful technique when values of a function are to be obtained from its Maclaurin expansion but convergence of the series is unacceptably slow or even nonexistent.
It is based on the hypothesis of the existence of a sequence of convergent rational approximations, as described in Baker and Graves--Morris (1981) and Graves--Morris (1979).
Unless there are reasons to the contrary (as discussed in Module 4, Section 2, Modules 5 and 6 of Baker and Graves--Morris (1981)), one normally uses the diagonal sequence of Padé approximants, namely
.. math::
\left\{\left[m/m\right],m = 0,1,2,\ldots \right\}\text{.}
Subsequent evaluation of the approximant at a given value of :math:`x` may be carried out using :meth:`pade_eval`.
.. _e02ra-py2-py-references:
**References**
Baker, G A Jr and Graves--Morris, P R, 1981, `Padé approximants, Part 1: Basic theory`, encyclopaedia of Mathematics and its Applications, Addison--Wesley
Graves--Morris, P R, 1979, `The numerical calculation of Padé approximants`, Padé Approximation and its Applications. Lecture Notes in Mathematics, (ed L Wuytack) (765), 231--245, Adison--Wesley
"""
raise NotImplementedError
[docs]def pade_eval(a, b, x):
r"""
``pade_eval`` evaluates a rational function at a user-supplied point, given the numerator and denominator coefficients.
.. _e02rb-py2-py-doc:
For full information please refer to the NAG Library document for e02rb
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02rbf.html
.. _e02rb-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(\textit{ia}\right)`
:math:`\mathrm{a}[\textit{j}]`, for :math:`\textit{j} = 1,2,\ldots,l+1`, must contain the value of the coefficient :math:`a_{\textit{j}}` in the numerator of the rational function.
**b** : float, array-like, shape :math:`\left(\textit{ib}\right)`
:math:`\mathrm{b}[\textit{k}]`, for :math:`\textit{k} = 1,2,\ldots,m+1`, must contain the value of the coefficient :math:`b_k` in the denominator of the rational function.
**x** : float
The point :math:`x` at which the rational function is to be evaluated.
**Returns**
**ans** : float
The result of evaluating the rational function at the given point :math:`x`.
.. _e02rb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
Evaluation at or near a pole.
(`errno` :math:`2`)
The first :math:`\textit{ib}` entries in :math:`\mathrm{b}` are zero: :math:`\textit{ib} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`2`)
On entry, :math:`\textit{ib} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ib}\geq 1`.
(`errno` :math:`2`)
On entry, :math:`\textit{ia} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ia}\geq 1`.
.. _e02rb-py2-py-notes:
**Notes**
Given a real value :math:`x` and the coefficients :math:`a_j`, for :math:`\textit{j} = 0,1,\ldots,l` and :math:`b_k`, for :math:`\textit{k} = 0,1,\ldots,m`, ``pade_eval`` evaluates the rational function
.. math::
\frac{{\sum_{{j = 0}}^la_jx^j}}{{\sum_{{k = 0}}^mb_kx^k}}\text{.}
using nested multiplication (see Conte and de Boor (1965)).
A particular use of ``pade_eval`` is to compute values of the Padé approximants determined by :meth:`pade_app`.
.. _e02rb-py2-py-references:
**References**
Conte, S D and de Boor, C, 1965, `Elementary Numerical Analysis`, McGraw--Hill
Peters, G and Wilkinson, J H, 1971, `Practical problems arising in the solution of polynomial equations`, J. Inst. Maths. Applics. (8), 16--35
"""
raise NotImplementedError
[docs]def dim2_spline_sort(lamda, mu, x, y):
r"""
``dim2_spline_sort`` sorts two-dimensional data into rectangular panels.
.. _e02za-py2-py-doc:
For full information please refer to the NAG Library document for e02za
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02zaf.html
.. _e02za-py2-py-parameters:
**Parameters**
**lamda** : float, array-like, shape :math:`\left(\textit{px}\right)`
:math:`\mathrm{lamda}[4]` to :math:`\mathrm{lamda}[\textit{px}-5]` must contain, in nondecreasing order, the intercepts on the :math:`x` axis of the sides of the panels parallel to the :math:`y` axis.
**mu** : float, array-like, shape :math:`\left(\textit{py}\right)`
:math:`\mathrm{mu}[4]` to :math:`\mathrm{mu}[\textit{py}-5]` must contain, in nondecreasing order, the intercepts on the :math:`y` axis of the sides of the panels parallel to the :math:`x` axis.
**x** : float, array-like, shape :math:`\left(m\right)`
The coordinates of the :math:`\textit{r}`\ th data point :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`.
**y** : float, array-like, shape :math:`\left(m\right)`
The coordinates of the :math:`\textit{r}`\ th data point :math:`\left(x_{\textit{r}}, y_{\textit{r}}\right)`, for :math:`\textit{r} = 1,2,\ldots,m`.
**Returns**
**point** : int, ndarray, shape :math:`\left(m+\left(\textit{px}-7\right)\times \left(\textit{py}-7\right)\right)`
For :math:`i = 1,2,\ldots,\textit{npoint}`, :math:`\mathrm{point}[m+i-1] = \mathrm{I1}` is the index of the first point in panel :math:`i`, :math:`\mathrm{point}[\mathrm{I1}-1] = \mathrm{I2}` is the index of the second point in panel :math:`i` and so on.
:math:`\mathrm{point}[\mathrm{In}-1] = 0` indicates that :math:`\mathrm{x}[\mathrm{In}-1],\mathrm{y}[\mathrm{In}-1]` was the last point in the panel.
The coordinates of points in panel :math:`i` can be accessed in turn by means of the following instructions:
::
i_n = point[m+i-1]
while i_n != 0:
xi = x[i_n-1]
yi = y[i_n-1]
.
.
.
i_n = point[i_n]
...
.. _e02za-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{mu}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{mu}[\textit{I}-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mu}[\textit{I}-1]\geq \mathrm{mu}[\textit{I}-2]`.
(`errno` :math:`1`)
On entry, :math:`\textit{I} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{lamda}[\textit{I}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{lamda}[\textit{I}-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lamda}[\textit{I}-1]\geq \mathrm{lamda}[\textit{I}-2]`.
(`errno` :math:`2`)
On entry, :math:`\textit{npoint} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{px} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{py} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{npoint}\geq m+\left(\textit{px}-7\right)\times \left(\textit{py}-7\right)`.
(`errno` :math:`2`)
On entry, :math:`\textit{py} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{py}\geq 8`.
(`errno` :math:`2`)
On entry, :math:`\textit{px} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{px}\geq 8`.
(`errno` :math:`2`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m > 0`.
.. _e02za-py2-py-notes:
**Notes**
A set of :math:`m` data points with rectangular Cartesian coordinates :math:`x_r,y_r` are sorted into panels defined by lines parallel to the :math:`y` and :math:`x` axes.
The intercepts of these lines on the :math:`x` and :math:`y` axes are given in :math:`\mathrm{lamda}[\textit{i}-1]`, for :math:`\textit{i} = 5,6,\ldots,\textit{px}-4` and :math:`\mathrm{mu}[\textit{j}-1]`, for :math:`\textit{j} = 5,6,\ldots,\textit{py}-4`, respectively.
The function orders the data so that all points in a panel occur before data in succeeding panels, where the panels are numbered from bottom to top and then left to right, with the usual arrangement of axes, as shown in the diagram.
Within a panel the points maintain their original order.
[figure omitted]
A data point lying exactly on one or more panel sides is taken to be in the highest-numbered panel adjacent to the point.
The function does not physically rearrange the data, but provides the array :math:`\mathrm{point}` which contains a linked list for each panel, pointing to the data in that panel.
The total number of panels is :math:`\left(\textit{px}-7\right)\times \left(\textit{py}-7\right)`.
"""
raise NotImplementedError
[docs]def opt_set(optstr, comm):
r"""
``opt_set`` either initializes or resets the option arrays or sets a single option for supported problem solving functions in submodule ``fit``.
Currently, only :meth:`dim2_spline_ts_sctr` is supported.
.. _e02zk-py2-py-doc:
For full information please refer to the NAG Library document for e02zk
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02zkf.html
.. _e02zk-py2-py-parameters:
**Parameters**
**optstr** : str
A string identifying the option to be set.
Initialize = :math:`\textit{function name}`
Initialize the option arrays :math:`\mathrm{comm}`\ ['iopts'] and :math:`\mathrm{comm}`\ ['opts'] for use with function :math:`\textit{function name}`, where :math:`\textit{function name}` is the name of the problem solving function you wish to use.
Defaults
Resets all options to their default values.
:math:`\textit{option} = \textit{optval}`
See :ref:`Other Parameters for dim2_spline_ts_sctr <e02jd-py2-py-other_params>` for details of valid values for :math:`\textit{option}` and :math:`\textit{optval}`. The equals sign (:math:`=`) delimiter must be used to separate the :math:`\textit{option}` from its :math:`\textit{optval}`.
The processing of :math:`\mathrm{optstr}` does not depend on its case.
Each token in the :math:`\textit{option}` and :math:`\textit{optval}` component must be separated by at least one space.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
.. _e02zk-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`11`)
On entry, the option in :math:`\mathrm{optstr}` was not recognized: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`12`)
On entry, the expected delimiter ':math:`=`' was not found in :math:`\mathrm{optstr}`: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`13`)
On entry, could not convert the specified :math:`\textit{optval}` to an integer: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`13`)
On entry, could not convert the specified :math:`\textit{optval}` to a real: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`14`)
On entry, attempting to initialize the option arrays but specified function name was not valid: :math:`\text{name} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`15`)
On entry, the :math:`\textit{optval}` supplied for the integer option is not valid.
:math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`16`)
On entry, the :math:`\textit{optval}` supplied for the real option is not valid.
:math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`17`)
On entry, the :math:`\textit{optval}` supplied for the character option is not valid.
:math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`21`)
On entry, either the option arrays have not been initialized or they have been corrupted.
.. _e02zk-py2-py-notes:
**Notes**
``opt_set`` has three purposes: to initialize option arrays, to reset all options to their default values or to set a single option to a user-supplied value.
Options and their values are, in general, presented as a character string, :math:`\mathrm{optstr}`, of the form ':math:`\textit{option} = \textit{optval}`'; alphabetic characters can be supplied in either upper or lower case.
Both :math:`\textit{option}` and :math:`\textit{optval}` may consist of one or more tokens separated by white space.
The tokens that comprise :math:`\textit{optval}` will normally be either an integer, real or character value as defined in the description of the specific optional argument.
In addition all options can take an :math:`\textit{optval}` 'DEFAULT' which resets the option to its default value.
It is imperative that option arrays are initialized before any options are set, before the relevant problem solving function is called and before any options are queried using :meth:`opt_get`.
Information relating to available option names and their corresponding valid values is given in :ref:`Other Parameters for dim2_spline_ts_sctr <e02jd-py2-py-other_params>`.
See Also
--------
:meth:`naginterfaces.library.examples.fit.dim2_spline_ts_sctr_ex.main`
"""
raise NotImplementedError
[docs]def opt_get(optstr, comm):
r"""
``opt_get`` is used to query the value of options available to supported problem solving functions in submodule ``fit``.
Currently, only :meth:`dim2_spline_ts_sctr` is supported.
.. _e02zl-py2-py-doc:
For full information please refer to the NAG Library document for e02zl
https://support.nag.com/numeric/nl/nagdoc_30/flhtml/e02/e02zlf.html
.. _e02zl-py2-py-parameters:
**Parameters**
**optstr** : str
A string identifying the option whose current value is required. See :ref:`Other Parameters for dim2_spline_ts_sctr <e02jd-py2-py-other_params>` for information on valid options. In addition, the following is a valid option:
Identify
``opt_get`` returns the function name supplied to :meth:`opt_set` when the option arrays :math:`\mathrm{comm}`\ ['iopts'] and :math:`\mathrm{comm}`\ ['opts'] were initialized.
**comm** : dict, communication object
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt_set`.
**Returns**
**optvalue** : dict
The option-value ``dict``, with the following keys:
``'value'`` : float, int or str
The value of the requested option.
``'annotation'`` : None or str
Possible additional information about the option value.
.. _e02zl-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`11`)
On entry, the option in :math:`\mathrm{optstr}` was not recognized: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`61`)
On entry, either the option arrays have not been initialized or they have been corrupted.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`41`)
On entry, :math:`\mathrm{optstr}` indicates a character option, but :math:`\textit{cvalue}` is too short to hold the stored value. The returned value will be truncated.
.. _e02zl-py2-py-notes:
**Notes**
``opt_get`` is used to query the current values of options.
It is necessary to initalize option arrays using :meth:`opt_set` before any options are queried.
Information on option names and whether these options are real, integer or character can be found in :ref:`Other Parameters for dim2_spline_ts_sctr <e02jd-py2-py-other_params>`.
"""
raise NotImplementedError
