geocat.comp.ndpolyval

geocat.comp.ndpolyval(p: Iterable, x: Iterable, axis: int = 0, **kwargs)

Extended version of numpy.polyval to support multi-dimensional outputs provided by geocat.comp.ndpolyfit.

As the name suggest, this version supports a multi-dimensional p array. Let’s say p is of dimension (s0,s1,s2) and axis=1, then the output would be of dimension (s0, M, s2) where M depends on x. The same way, x could be a multi dimensional array or a single array. In another word, x is either of dimension (M, ), (M, 1), (1, M) or, in this example, of dimension (s0, M, s2). When x is not the vector, it must have the same dimension as of p except for the dimension that is defined by axis.

Parameters

  • p (Iterable) – the polynomial coeficients

  • x (Iterable) – the coordinates where polynomial must be evaluated

  • axis (int) –

  • axis where the polynomials are there. (The) –

  • **kwargs – Currently not used.

Returns (:class: xr.DataArray:

polynomial evaluated with the provided coordinates.

Examples

  • Evaluating a polynomial:

>>> p = [2, 3] # representing y = 2*x + 3
>>> x = np.arange(-5, 6, dtype="float")
>>> x
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.])
>>> from geocat.comp.polynomial import ndpolyval
>>> y = ndpolyval(p, x)
>>> y.shape
(11,)
>>> print(y)
<xarray.DataArray (dim_0: 11)>
array([-7., -5., -3., -1.,  1.,  3.,  5.,  7.,  9., 11., 13.])
Dimensions without coordinates: dim_0
>>> np.testing.assert_almost_equal(y, 2*x+3)

  • evaluating a multi-dimensional fitted polynomial:

>>> p = np.tile(np.asarray(p).reshape(1, 2, 1, 1), [3, 1, 4, 5])
>>> p.shape
(3, 2, 4, 5)
>>> p[1, :, 1, 1]
array([2, 3])
>>> p[1, :, 1, 2]
array([2, 3])
>>> y = ndpolyval(p, x, axis=1)
>>> y.shape
(3, 11, 4, 5)

  • Fitting a first degree polynomial and calculating the residual manually:

>>> from geocat.comp.polynomial import ndpolyfit, ndpolyval
>>> import numpy as np
>>> n = 10
>>> x = np.linspace(0, 1, n)
>>> print(x.shape)
(10,)
>>> y = np.random.random((3, 4, n, 2, 5))
>>> print(y.shape)
(3, 4, 10, 2, 5)
>>> # First fitting a first degree polynomial to our 5-Dimensional array
... p = ndpolyfit(x, y, deg=1, axis=2)
>>> print("p.shape: ", p.shape)
p.shape:  (3, 4, 2, 2, 5)
>>> # Now re-evaluating the polynomial at the same points:
... y_fitted = ndpolyval(p, x, axis=2)
>>> We are ready to manually calculate the residual now:
... signed_residual = y_fitted - y
>>> y_fitted.shape
(3, 4, 10, 2, 5)

  • when evaluating x can be another multi-dimensional array

>>> # Let's work on the p that was calculated in previous example
... # using the command:
... #           p = ndpolyfit(x, y, deg=1, axis=2)
... print("p.shape: ", p.shape)
p.shape:  (3, 4, 2, 2, 5)

>>> # Let's pass x which has 20 members this time:
... x = np.random.random(20)
>>> print("x.shape: ", x.shape)
x.shape:  (20,)
>>> new_y = ndpolyval(p, x, axis=2)
>>> print("new_y.shape: ", new_y.shape) # Note new_y.shape[2] = 20
new_y.shape:  (3, 4, 20, 2, 5)

>>> # Now let's make a multi-dimensional x
... # NOTE: all dimension, except axis=2 is the same as the one in p
... #       because during fitting step axis was set to 2
... x = np.random.random((3, 4, 42, 2, 5))
>>> new_y = ndpolyval(p, x, axis=2)
>>> print("new_y.shape: ", new_y.shape) # Note new_y.shape[2] = 42
new_y.shape:  (3, 4, 42, 2, 5)

>>> # now if a wrongly sized array is passed, you would get an error:
... x = np.random.random((5, 2, 42, 4, 3))
>>> new_y = ndpolyval(p, x, axis=2)
Traceback (most recent call last):
  ...
ValueError: x has invalid shape.