diff options
Diffstat (limited to 'numpy/polynomial/hermite.py')
| -rw-r--r-- | numpy/polynomial/hermite.py | 235 |
1 files changed, 22 insertions, 213 deletions
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py index 93c9fc564..3870d1054 100644 --- a/numpy/polynomial/hermite.py +++ b/numpy/polynomial/hermite.py @@ -286,21 +286,7 @@ def hermfromroots(roots): array([0.+0.j, 0.+0.j]) """ - if len(roots) == 0: - return np.ones(1) - else: - [roots] = pu.as_series([roots], trim=False) - roots.sort() - p = [hermline(-r, 1) for r in roots] - n = len(p) - while n > 1: - m, r = divmod(n, 2) - tmp = [hermmul(p[i], p[i+m]) for i in range(m)] - if r: - tmp[0] = hermmul(tmp[0], p[-1]) - p = tmp - n = m - return p[0] + return pu._fromroots(hermline, hermmul, roots) def hermadd(c1, c2): @@ -340,15 +326,7 @@ def hermadd(c1, c2): array([2., 4., 6., 4.]) """ - # c1, c2 are trimmed copies - [c1, c2] = pu.as_series([c1, c2]) - if len(c1) > len(c2): - c1[:c2.size] += c2 - ret = c1 - else: - c2[:c1.size] += c1 - ret = c2 - return pu.trimseq(ret) + return pu._add(c1, c2) def hermsub(c1, c2): @@ -388,16 +366,7 @@ def hermsub(c1, c2): array([0., 0., 0., 4.]) """ - # c1, c2 are trimmed copies - [c1, c2] = pu.as_series([c1, c2]) - if len(c1) > len(c2): - c1[:c2.size] -= c2 - ret = c1 - else: - c2 = -c2 - c2[:c1.size] += c1 - ret = c2 - return pu.trimseq(ret) + return pu._sub(c1, c2) def hermmulx(c): @@ -564,26 +533,7 @@ def hermdiv(c1, c2): (array([1., 2., 3.]), array([1., 1.])) """ - # c1, c2 are trimmed copies - [c1, c2] = pu.as_series([c1, c2]) - if c2[-1] == 0: - raise ZeroDivisionError() - - lc1 = len(c1) - lc2 = len(c2) - if lc1 < lc2: - return c1[:1]*0, c1 - elif lc2 == 1: - return c1/c2[-1], c1[:1]*0 - else: - quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype) - rem = c1 - for i in range(lc1 - lc2, - 1, -1): - p = hermmul([0]*i + [1], c2) - q = rem[-1]/p[-1] - rem = rem[:-1] - q*p[:-1] - quo[i] = q - return quo, pu.trimseq(rem) + return pu._div(hermmul, c1, c2) def hermpow(c, pow, maxpower=16): @@ -620,24 +570,7 @@ def hermpow(c, pow, maxpower=16): array([81., 52., 82., 12., 9.]) """ - # c is a trimmed copy - [c] = pu.as_series([c]) - power = int(pow) - if power != pow or power < 0: - raise ValueError("Power must be a non-negative integer.") - elif maxpower is not None and power > maxpower: - raise ValueError("Power is too large") - elif power == 0: - return np.array([1], dtype=c.dtype) - elif power == 1: - return c - 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): - prd = hermmul(prd, c) - return prd + return pu._pow(hermmul, c, pow, maxpower) def hermder(c, m=1, scl=1, axis=0): @@ -698,14 +631,10 @@ def hermder(c, m=1, scl=1, axis=0): c = np.array(c, ndmin=1, copy=1) if c.dtype.char in '?bBhHiIlLqQpP': c = c.astype(np.double) - cnt, iaxis = [int(t) for t in [m, axis]] - - if cnt != m: - raise ValueError("The order of derivation must be integer") + cnt = pu._deprecate_as_int(m, "the order of derivation") + iaxis = pu._deprecate_as_int(axis, "the axis") if cnt < 0: raise ValueError("The order of derivation must be non-negative") - if iaxis != axis: - raise ValueError("The axis must be integer") iaxis = normalize_axis_index(iaxis, c.ndim) if cnt == 0: @@ -815,10 +744,8 @@ def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0): c = c.astype(np.double) if not np.iterable(k): k = [k] - cnt, iaxis = [int(t) for t in [m, axis]] - - if cnt != m: - raise ValueError("The order of integration must be integer") + cnt = pu._deprecate_as_int(m, "the order of integration") + iaxis = pu._deprecate_as_int(axis, "the axis") if cnt < 0: raise ValueError("The order of integration must be non-negative") if len(k) > cnt: @@ -827,8 +754,6 @@ def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0): raise ValueError("lbnd must be a scalar.") if np.ndim(scl) != 0: raise ValueError("scl must be a scalar.") - if iaxis != axis: - raise ValueError("The axis must be integer") iaxis = normalize_axis_index(iaxis, c.ndim) if cnt == 0: @@ -995,14 +920,7 @@ def hermval2d(x, y, c): .. versionadded:: 1.7.0 """ - try: - x, y = np.array((x, y), copy=0) - except Exception: - raise ValueError('x, y are incompatible') - - c = hermval(x, c) - c = hermval(y, c, tensor=False) - return c + return pu._valnd(hermval, c, x, y) def hermgrid2d(x, y, c): @@ -1055,9 +973,7 @@ def hermgrid2d(x, y, c): .. versionadded:: 1.7.0 """ - c = hermval(x, c) - c = hermval(y, c) - return c + return pu._gridnd(hermval, c, x, y) def hermval3d(x, y, z, c): @@ -1108,15 +1024,7 @@ def hermval3d(x, y, z, c): .. versionadded:: 1.7.0 """ - try: - x, y, z = np.array((x, y, z), copy=0) - except Exception: - raise ValueError('x, y, z are incompatible') - - c = hermval(x, c) - c = hermval(y, c, tensor=False) - c = hermval(z, c, tensor=False) - return c + return pu._valnd(hermval, c, x, y, z) def hermgrid3d(x, y, z, c): @@ -1172,10 +1080,7 @@ def hermgrid3d(x, y, z, c): .. versionadded:: 1.7.0 """ - c = hermval(x, c) - c = hermval(y, c) - c = hermval(z, c) - return c + return pu._gridnd(hermval, c, x, y, z) def hermvander(x, deg): @@ -1222,9 +1127,7 @@ def hermvander(x, deg): [ 1., 2., 2., -4.]]) """ - ideg = int(deg) - if ideg != deg: - raise ValueError("deg must be integer") + ideg = pu._deprecate_as_int(deg, "deg") if ideg < 0: raise ValueError("deg must be non-negative") @@ -1291,17 +1194,7 @@ def hermvander2d(x, y, deg): .. versionadded:: 1.7.0 """ - ideg = [int(d) for d in deg] - is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)] - if is_valid != [1, 1]: - raise ValueError("degrees must be non-negative integers") - degx, degy = ideg - x, y = np.array((x, y), copy=0) + 0.0 - - vx = hermvander(x, degx) - vy = hermvander(y, degy) - v = vx[..., None]*vy[..., None,:] - return v.reshape(v.shape[:-2] + (-1,)) + return pu._vander2d(hermvander, x, y, deg) def hermvander3d(x, y, z, deg): @@ -1355,18 +1248,7 @@ def hermvander3d(x, y, z, deg): .. versionadded:: 1.7.0 """ - ideg = [int(d) for d in deg] - is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)] - if is_valid != [1, 1, 1]: - raise ValueError("degrees must be non-negative integers") - degx, degy, degz = ideg - x, y, z = np.array((x, y, z), copy=0) + 0.0 - - vx = hermvander(x, degx) - vy = hermvander(y, degy) - vz = hermvander(z, degz) - v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:] - return v.reshape(v.shape[:-3] + (-1,)) + return pu._vander3d(hermvander, x, y, z, deg) def hermfit(x, y, deg, rcond=None, full=False, w=None): @@ -1493,81 +1375,7 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None): array([1.0218, 1.9986, 2.9999]) # may vary """ - x = np.asarray(x) + 0.0 - y = np.asarray(y) + 0.0 - deg = np.asarray(deg) - - # check arguments. - if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0: - raise TypeError("deg must be an int or non-empty 1-D array of int") - if deg.min() < 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: - raise TypeError("expected 1D or 2D array for y") - if len(x) != len(y): - raise TypeError("expected x and y to have same length") - - if deg.ndim == 0: - lmax = deg - order = lmax + 1 - van = hermvander(x, lmax) - else: - deg = np.sort(deg) - lmax = deg[-1] - order = len(deg) - van = hermvander(x, lmax)[:, deg] - - # set up the least squares matrices in transposed form - lhs = van.T - rhs = y.T - if w is not None: - w = np.asarray(w) + 0.0 - if w.ndim != 1: - raise TypeError("expected 1D vector for w") - if len(x) != len(w): - raise TypeError("expected x and w to have same length") - # apply weights. Don't use inplace operations as they - # can cause problems with NA. - lhs = lhs * w - rhs = rhs * w - - # set rcond - if rcond is None: - rcond = len(x)*np.finfo(x.dtype).eps - - # Determine the norms of the design matrix columns. - if issubclass(lhs.dtype.type, np.complexfloating): - scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1)) - else: - scl = np.sqrt(np.square(lhs).sum(1)) - scl[scl == 0] = 1 - - # Solve the least squares problem. - c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond) - c = (c.T/scl).T - - # Expand c to include non-fitted coefficients which are set to zero - if deg.ndim > 0: - if c.ndim == 2: - cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype) - else: - cc = np.zeros(lmax+1, dtype=c.dtype) - cc[deg] = c - c = cc - - # warn on rank reduction - if rank != order and not full: - msg = "The fit may be poorly conditioned" - warnings.warn(msg, pu.RankWarning, stacklevel=2) - - if full: - return c, [resids, rank, s, rcond] - else: - return c + return pu._fit(hermvander, x, y, deg, rcond, full, w) def hermcompanion(c): @@ -1668,7 +1476,8 @@ def hermroots(c): if len(c) == 2: return np.array([-.5*c[0]/c[1]]) - m = hermcompanion(c) + # rotated companion matrix reduces error + m = hermcompanion(c)[::-1,::-1] r = la.eigvals(m) r.sort() return r @@ -1753,9 +1562,9 @@ def hermgauss(deg): the right value when integrating 1. """ - ideg = int(deg) - if ideg != deg or ideg < 1: - raise ValueError("deg must be a non-negative integer") + ideg = pu._deprecate_as_int(deg, "deg") + if ideg <= 0: + raise ValueError("deg must be a positive integer") # first approximation of roots. We use the fact that the companion # matrix is symmetric in this case in order to obtain better zeros. |