diff options
Diffstat (limited to 'numpy/polynomial/polyutils.py')
| -rw-r--r-- | numpy/polynomial/polyutils.py | 126 |
1 files changed, 83 insertions, 43 deletions
diff --git a/numpy/polynomial/polyutils.py b/numpy/polynomial/polyutils.py index 35b24d1ab..ec7ba6f1d 100644 --- a/numpy/polynomial/polyutils.py +++ b/numpy/polynomial/polyutils.py @@ -43,9 +43,8 @@ Functions mapparms parameters of the linear map between domains. """ -from __future__ import division, absolute_import, print_function - import operator +import functools import warnings import numpy as np @@ -79,7 +78,7 @@ class PolyDomainError(PolyError): # Base class for all polynomial types # -class PolyBase(object): +class PolyBase: """ Base class for all polynomial types. @@ -415,45 +414,89 @@ def mapdomain(x, old, new): return off + scl*x -def _vander2d(vander_f, x, y, deg): - """ - Helper function used to implement the ``<type>vander2d`` functions. +def _nth_slice(i, ndim): + sl = [np.newaxis] * ndim + sl[i] = slice(None) + return tuple(sl) + + +def _vander_nd(vander_fs, points, degrees): + r""" + A generalization of the Vandermonde matrix for N dimensions + + The result is built by combining the results of 1d Vandermonde matrices, + + .. math:: + W[i_0, \ldots, i_M, j_0, \ldots, j_N] = \prod_{k=0}^N{V_k(x_k)[i_0, \ldots, i_M, j_k]} + + where + + .. math:: + N &= \texttt{len(points)} = \texttt{len(degrees)} = \texttt{len(vander\_fs)} \\ + M &= \texttt{points[k].ndim} \\ + V_k &= \texttt{vander\_fs[k]} \\ + x_k &= \texttt{points[k]} \\ + 0 \le j_k &\le \texttt{degrees[k]} + + Expanding the one-dimensional :math:`V_k` functions gives: + + .. math:: + W[i_0, \ldots, i_M, j_0, \ldots, j_N] = \prod_{k=0}^N{B_{k, j_k}(x_k[i_0, \ldots, i_M])} + + where :math:`B_{k,m}` is the m'th basis of the polynomial construction used along + dimension :math:`k`. For a regular polynomial, :math:`B_{k, m}(x) = P_m(x) = x^m`. Parameters ---------- - vander_f : function(array_like, int) -> ndarray - The 1d vander function, such as ``polyvander`` - x, y, deg : - See the ``<type>vander2d`` functions for more detail + vander_fs : Sequence[function(array_like, int) -> ndarray] + The 1d vander function to use for each axis, such as ``polyvander`` + points : Sequence[array_like] + Arrays of point coordinates, all of the same shape. The dtypes + will be converted to either float64 or complex128 depending on + whether any of the elements are complex. Scalars are converted to + 1-D arrays. + This must be the same length as `vander_fs`. + degrees : Sequence[int] + The maximum degree (inclusive) to use for each axis. + This must be the same length as `vander_fs`. + + Returns + ------- + vander_nd : ndarray + An array of shape ``points[0].shape + tuple(d + 1 for d in degrees)``. """ - degx, degy = deg - x, y = np.array((x, y), copy=False) + 0.0 + n_dims = len(vander_fs) + if n_dims != len(points): + raise ValueError( + f"Expected {n_dims} dimensions of sample points, got {len(points)}") + if n_dims != len(degrees): + raise ValueError( + f"Expected {n_dims} dimensions of degrees, got {len(degrees)}") + if n_dims == 0: + raise ValueError("Unable to guess a dtype or shape when no points are given") - vx = vander_f(x, degx) - vy = vander_f(y, degy) - v = vx[..., None]*vy[..., None,:] - return v.reshape(v.shape[:-2] + (-1,)) + # convert to the same shape and type + points = tuple(np.array(tuple(points), copy=False) + 0.0) + # produce the vandermonde matrix for each dimension, placing the last + # axis of each in an independent trailing axis of the output + vander_arrays = ( + vander_fs[i](points[i], degrees[i])[(...,) + _nth_slice(i, n_dims)] + for i in range(n_dims) + ) -def _vander3d(vander_f, x, y, z, deg): - """ - Helper function used to implement the ``<type>vander3d`` functions. + # we checked this wasn't empty already, so no `initial` needed + return functools.reduce(operator.mul, vander_arrays) - Parameters - ---------- - vander_f : function(array_like, int) -> ndarray - The 1d vander function, such as ``polyvander`` - x, y, z, deg : - See the ``<type>vander3d`` functions for more detail + +def _vander_nd_flat(vander_fs, points, degrees): """ - degx, degy, degz = deg - x, y, z = np.array((x, y, z), copy=False) + 0.0 + Like `_vander_nd`, but flattens the last ``len(degrees)`` axes into a single axis - vx = vander_f(x, degx) - vy = vander_f(y, degy) - vz = vander_f(z, degz) - v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:] - return v.reshape(v.shape[:-3] + (-1,)) + Used to implement the public ``<type>vander<n>d`` functions. + """ + v = _vander_nd(vander_fs, points, degrees) + return v.reshape(v.shape[:-len(degrees)] + (-1,)) def _fromroots(line_f, mul_f, roots): @@ -497,17 +540,15 @@ def _valnd(val_f, c, *args): c, args : See the ``<type>val<n>d`` functions for more detail """ - try: - args = tuple(np.array(args, copy=False)) - except Exception: - # preserve the old error message - if len(args) == 2: + args = [np.asanyarray(a) for a in args] + shape0 = args[0].shape + if not all((a.shape == shape0 for a in args[1:])): + if len(args) == 3: raise ValueError('x, y, z are incompatible') - elif len(args) == 3: + elif len(args) == 2: raise ValueError('x, y are incompatible') else: raise ValueError('ordinates are incompatible') - it = iter(args) x0 = next(it) @@ -745,12 +786,11 @@ def _deprecate_as_int(x, desc): else: if ix == x: warnings.warn( - "In future, this will raise TypeError, as {} will need to " - "be an integer not just an integral float." - .format(desc), + f"In future, this will raise TypeError, as {desc} will " + "need to be an integer not just an integral float.", DeprecationWarning, stacklevel=3 ) return ix - raise TypeError("{} must be an integer".format(desc)) + raise TypeError(f"{desc} must be an integer") |