diff options
Diffstat (limited to 'numpy/polynomial/hermite_e.py')
| -rw-r--r-- | numpy/polynomial/hermite_e.py | 178 |
1 files changed, 110 insertions, 68 deletions
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index 874b42470..fce13a84e 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -66,18 +66,19 @@ import numpy.linalg as la from . import polyutils as pu from ._polybase import ABCPolyBase -__all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline', - 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermpow', - 'hermeval', - 'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots', - 'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE', - 'hermeval2d', 'hermeval3d', 'hermegrid2d', 'hermegrid3d', 'hermevander2d', - 'hermevander3d', 'hermecompanion', 'hermegauss', 'hermeweight'] +__all__ = [ + 'hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline', + 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', + 'hermepow', 'hermeval', 'hermeder', 'hermeint', 'herme2poly', + 'poly2herme', 'hermefromroots', 'hermevander', 'hermefit', 'hermetrim', + 'hermeroots', 'HermiteE', 'hermeval2d', 'hermeval3d', 'hermegrid2d', + 'hermegrid3d', 'hermevander2d', 'hermevander3d', 'hermecompanion', + 'hermegauss', 'hermeweight'] hermetrim = pu.trimcoef -def poly2herme(pol) : +def poly2herme(pol): """ poly2herme(pol) @@ -118,12 +119,12 @@ def poly2herme(pol) : [pol] = pu.as_series([pol]) deg = len(pol) - 1 res = 0 - for i in range(deg, -1, -1) : + for i in range(deg, -1, -1): res = hermeadd(hermemulx(res), pol[i]) return res -def herme2poly(c) : +def herme2poly(c): """ Convert a Hermite series to a polynomial. @@ -173,7 +174,7 @@ def herme2poly(c) : c0 = c[-2] c1 = c[-1] # i is the current degree of c1 - for i in range(n - 1, 1, -1) : + for i in range(n - 1, 1, -1): tmp = c0 c0 = polysub(c[i - 2], c1*(i - 1)) c1 = polyadd(tmp, polymulx(c1)) @@ -197,7 +198,7 @@ hermeone = np.array([1]) hermex = np.array([0, 1]) -def hermeline(off, scl) : +def hermeline(off, scl): """ Hermite series whose graph is a straight line. @@ -228,13 +229,13 @@ def hermeline(off, scl) : 5.0 """ - if scl != 0 : + if scl != 0: return np.array([off, scl]) - else : + else: return np.array([off]) -def hermefromroots(roots) : +def hermefromroots(roots): """ Generate a HermiteE series with given roots. @@ -284,9 +285,9 @@ def hermefromroots(roots) : array([ 0.+0.j, 0.+0.j]) """ - if len(roots) == 0 : + if len(roots) == 0: return np.ones(1) - else : + else: [roots] = pu.as_series([roots], trim=False) roots.sort() p = [hermeline(-r, 1) for r in roots] @@ -340,10 +341,10 @@ def hermeadd(c1, c2): """ # c1, c2 are trimmed copies [c1, c2] = pu.as_series([c1, c2]) - if len(c1) > len(c2) : + if len(c1) > len(c2): c1[:c2.size] += c2 ret = c1 - else : + else: c2[:c1.size] += c1 ret = c2 return pu.trimseq(ret) @@ -388,10 +389,10 @@ def hermesub(c1, c2): """ # c1, c2 are trimmed copies [c1, c2] = pu.as_series([c1, c2]) - if len(c1) > len(c2) : + if len(c1) > len(c2): c1[:c2.size] -= c2 ret = c1 - else : + else: c2 = -c2 c2[:c1.size] += c1 ret = c2 @@ -501,13 +502,13 @@ def hermemul(c1, c2): elif len(c) == 2: c0 = c[0]*xs c1 = c[1]*xs - else : + else: nd = len(c) c0 = c[-2]*xs c1 = c[-1]*xs - for i in range(3, len(c) + 1) : + for i in range(3, len(c) + 1): tmp = c0 - nd = nd - 1 + nd = nd - 1 c0 = hermesub(c[-i]*xs, c1*(nd - 1)) c1 = hermeadd(tmp, hermemulx(c1)) return hermeadd(c0, hermemulx(c1)) @@ -558,16 +559,16 @@ def hermediv(c1, c2): """ # c1, c2 are trimmed copies [c1, c2] = pu.as_series([c1, c2]) - if c2[-1] == 0 : + if c2[-1] == 0: raise ZeroDivisionError() lc1 = len(c1) lc2 = len(c2) - if lc1 < lc2 : + if lc1 < lc2: return c1[:1]*0, c1 - elif lc2 == 1 : + elif lc2 == 1: return c1/c2[-1], c1[:1]*0 - else : + else: quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype) rem = c1 for i in range(lc1 - lc2, - 1, -1): @@ -578,7 +579,7 @@ def hermediv(c1, c2): return quo, pu.trimseq(rem) -def hermepow(c, pow, maxpower=16) : +def hermepow(c, pow, maxpower=16): """Raise a Hermite series to a power. Returns the Hermite series `c` raised to the power `pow`. The @@ -615,24 +616,24 @@ def hermepow(c, pow, maxpower=16) : # c is a trimmed copy [c] = pu.as_series([c]) power = int(pow) - if power != pow or power < 0 : + if power != pow or power < 0: raise ValueError("Power must be a non-negative integer.") - elif maxpower is not None and power > maxpower : + elif maxpower is not None and power > maxpower: raise ValueError("Power is too large") - elif power == 0 : + elif power == 0: return np.array([1], dtype=c.dtype) - elif power == 1 : + elif power == 1: return c - else : + else: # This can be made more efficient by using powers of two # in the usual way. prd = c - for i in range(2, power + 1) : + for i in range(2, power + 1): prd = hermemul(prd, c) return prd -def hermeder(c, m=1, scl=1, axis=0) : +def hermeder(c, m=1, scl=1, axis=0): """ Differentiate a Hermite_e series. @@ -710,7 +711,7 @@ def hermeder(c, m=1, scl=1, axis=0) : n = len(c) if cnt >= n: return c[:1]*0 - else : + else: for i in range(cnt): n = n - 1 c *= scl @@ -814,9 +815,9 @@ def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0): if cnt != m: raise ValueError("The order of integration must be integer") - if cnt < 0 : + if cnt < 0: raise ValueError("The order of integration must be non-negative") - if len(k) > cnt : + if len(k) > cnt: raise ValueError("Too many integration constants") if iaxis != axis: raise ValueError("The axis must be integer") @@ -830,7 +831,7 @@ def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0): c = np.rollaxis(c, iaxis) k = list(k) + [0]*(cnt - len(k)) - for i in range(cnt) : + for i in range(cnt): n = len(c) c *= scl if n == 1 and np.all(c[0] == 0): @@ -922,21 +923,21 @@ def hermeval(x, c, tensor=True): if isinstance(x, (tuple, list)): x = np.asarray(x) if isinstance(x, np.ndarray) and tensor: - c = c.reshape(c.shape + (1,)*x.ndim) + c = c.reshape(c.shape + (1,)*x.ndim) - if len(c) == 1 : + if len(c) == 1: c0 = c[0] c1 = 0 - elif len(c) == 2 : + elif len(c) == 2: c0 = c[0] c1 = c[1] - else : + else: nd = len(c) c0 = c[-2] c1 = c[-1] - for i in range(3, len(c) + 1) : + for i in range(3, len(c) + 1): tmp = c0 - nd = nd - 1 + nd = nd - 1 c0 = c[-i] - c1*(nd - 1) c1 = tmp + c1*x return c0 + c1*x @@ -1171,7 +1172,7 @@ def hermegrid3d(x, y, z, c): return c -def hermevander(x, deg) : +def hermevander(x, deg): """Pseudo-Vandermonde matrix of given degree. Returns the pseudo-Vandermonde matrix of degree `deg` and sample points @@ -1226,14 +1227,14 @@ def hermevander(x, deg) : dtyp = x.dtype v = np.empty(dims, dtype=dtyp) v[0] = x*0 + 1 - if ideg > 0 : + if ideg > 0: v[1] = x - for i in range(2, ideg + 1) : + for i in range(2, ideg + 1): v[i] = (v[i-1]*x - v[i-2]*(i - 1)) return np.rollaxis(v, 0, v.ndim) -def hermevander2d(x, y, deg) : +def hermevander2d(x, y, deg): """Pseudo-Vandermonde matrix of given degrees. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample @@ -1296,7 +1297,7 @@ def hermevander2d(x, y, deg) : return v.reshape(v.shape[:-2] + (-1,)) -def hermevander3d(x, y, z, deg) : +def hermevander3d(x, y, z, deg): """Pseudo-Vandermonde matrix of given degrees. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample @@ -1487,13 +1488,13 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): y = np.asarray(y) + 0.0 # check arguments. - if deg < 0 : + if deg < 0: raise ValueError("expected deg >= 0") if x.ndim != 1: raise TypeError("expected 1D vector for x") if x.size == 0: raise TypeError("expected non-empty vector for x") - if y.ndim < 1 or y.ndim > 2 : + if y.ndim < 1 or y.ndim > 2: raise TypeError("expected 1D or 2D array for y") if len(x) != len(y): raise TypeError("expected x and y to have same length") @@ -1513,7 +1514,7 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): rhs = rhs * w # set rcond - if rcond is None : + if rcond is None: rcond = len(x)*np.finfo(x.dtype).eps # Determine the norms of the design matrix columns. @@ -1532,9 +1533,9 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): msg = "The fit may be poorly conditioned" warnings.warn(msg, pu.RankWarning) - if full : + if full: return c, [resids, rank, s, rcond] - else : + else: return c @@ -1565,7 +1566,6 @@ def hermecompanion(c): .. versionadded::1.7.0 """ - accprod = np.multiply.accumulate # c is a trimmed copy [c] = pu.as_series([c]) if len(c) < 2: @@ -1575,13 +1575,13 @@ def hermecompanion(c): n = len(c) - 1 mat = np.zeros((n, n), dtype=c.dtype) - scl = np.hstack((1., np.sqrt(np.arange(1, n)))) - scl = np.multiply.accumulate(scl) + scl = np.hstack((1., 1./np.sqrt(np.arange(n - 1, 0, -1)))) + scl = np.multiply.accumulate(scl)[::-1] top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = np.sqrt(np.arange(1, n)) bot[...] = top - mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1]) + mat[:, -1] -= scl*c[:-1]/c[-1] return mat @@ -1633,9 +1633,9 @@ def hermeroots(c): """ # c is a trimmed copy [c] = pu.as_series([c]) - if len(c) <= 1 : + if len(c) <= 1: return np.array([], dtype=c.dtype) - if len(c) == 2 : + if len(c) == 2: return np.array([-c[0]/c[1]]) m = hermecompanion(c) @@ -1644,6 +1644,49 @@ def hermeroots(c): return r +def _normed_hermite_e_n(x, n): + """ + Evaluate a normalized HermiteE polynomial. + + Compute the value of the normalized HermiteE polynomial of degree ``n`` + at the points ``x``. + + + Parameters + ---------- + x : ndarray of double. + Points at which to evaluate the function + n : int + Degree of the normalized HermiteE function to be evaluated. + + Returns + ------- + values : ndarray + The shape of the return value is described above. + + Notes + ----- + .. versionadded:: 1.10.0 + + This function is needed for finding the Gauss points and integration + weights for high degrees. The values of the standard HermiteE functions + overflow when n >= 207. + + """ + if n == 0: + return np.ones(x.shape)/np.sqrt(np.sqrt(2*np.pi)) + + c0 = 0. + c1 = 1./np.sqrt(np.sqrt(2*np.pi)) + nd = float(n) + for i in range(n - 1): + tmp = c0 + c0 = -c1*np.sqrt((nd - 1.)/nd) + c1 = tmp + c1*x*np.sqrt(1./nd) + nd = nd - 1.0 + return c0 + c1*x + + def hermegauss(deg): """ Gauss-HermiteE quadrature. @@ -1688,20 +1731,19 @@ def hermegauss(deg): # matrix is symmetric in this case in order to obtain better zeros. c = np.array([0]*deg + [1]) m = hermecompanion(c) - x = la.eigvals(m) + x = la.eigvalsh(m) x.sort() # improve roots by one application of Newton - dy = hermeval(x, c) - df = hermeval(x, hermeder(c)) + dy = _normed_hermite_e_n(x, ideg) + df = _normed_hermite_e_n(x, ideg - 1) * np.sqrt(ideg) x -= dy/df # compute the weights. We scale the factor to avoid possible numerical # overflow. - fm = hermeval(x, c[1:]) + fm = _normed_hermite_e_n(x, ideg - 1) fm /= np.abs(fm).max() - df /= np.abs(df).max() - w = 1/(fm * df) + w = 1/(fm * fm) # for Hermite_e we can also symmetrize w = (w + w[::-1])/2 |