comparison DEPENDENCIES/mingw32/Python27/Lib/site-packages/numpy/polynomial/hermite_e.py @ 87:2a2c65a20a8b

Add Python libs and headers
author Chris Cannam
date Wed, 25 Feb 2015 14:05:22 +0000
parents
children
comparison
equal deleted inserted replaced
86:413a9d26189e 87:2a2c65a20a8b
1 """
2 Objects for dealing with Hermite_e series.
3
4 This module provides a number of objects (mostly functions) useful for
5 dealing with Hermite_e series, including a `HermiteE` class that
6 encapsulates the usual arithmetic operations. (General information
7 on how this module represents and works with such polynomials is in the
8 docstring for its "parent" sub-package, `numpy.polynomial`).
9
10 Constants
11 ---------
12 - `hermedomain` -- Hermite_e series default domain, [-1,1].
13 - `hermezero` -- Hermite_e series that evaluates identically to 0.
14 - `hermeone` -- Hermite_e series that evaluates identically to 1.
15 - `hermex` -- Hermite_e series for the identity map, ``f(x) = x``.
16
17 Arithmetic
18 ----------
19 - `hermemulx` -- multiply a Hermite_e series in ``P_i(x)`` by ``x``.
20 - `hermeadd` -- add two Hermite_e series.
21 - `hermesub` -- subtract one Hermite_e series from another.
22 - `hermemul` -- multiply two Hermite_e series.
23 - `hermediv` -- divide one Hermite_e series by another.
24 - `hermeval` -- evaluate a Hermite_e series at given points.
25 - `hermeval2d` -- evaluate a 2D Hermite_e series at given points.
26 - `hermeval3d` -- evaluate a 3D Hermite_e series at given points.
27 - `hermegrid2d` -- evaluate a 2D Hermite_e series on a Cartesian product.
28 - `hermegrid3d` -- evaluate a 3D Hermite_e series on a Cartesian product.
29
30 Calculus
31 --------
32 - `hermeder` -- differentiate a Hermite_e series.
33 - `hermeint` -- integrate a Hermite_e series.
34
35 Misc Functions
36 --------------
37 - `hermefromroots` -- create a Hermite_e series with specified roots.
38 - `hermeroots` -- find the roots of a Hermite_e series.
39 - `hermevander` -- Vandermonde-like matrix for Hermite_e polynomials.
40 - `hermevander2d` -- Vandermonde-like matrix for 2D power series.
41 - `hermevander3d` -- Vandermonde-like matrix for 3D power series.
42 - `hermegauss` -- Gauss-Hermite_e quadrature, points and weights.
43 - `hermeweight` -- Hermite_e weight function.
44 - `hermecompanion` -- symmetrized companion matrix in Hermite_e form.
45 - `hermefit` -- least-squares fit returning a Hermite_e series.
46 - `hermetrim` -- trim leading coefficients from a Hermite_e series.
47 - `hermeline` -- Hermite_e series of given straight line.
48 - `herme2poly` -- convert a Hermite_e series to a polynomial.
49 - `poly2herme` -- convert a polynomial to a Hermite_e series.
50
51 Classes
52 -------
53 - `HermiteE` -- A Hermite_e series class.
54
55 See also
56 --------
57 `numpy.polynomial`
58
59 """
60 from __future__ import division, absolute_import, print_function
61
62 import warnings
63 import numpy as np
64 import numpy.linalg as la
65
66 from . import polyutils as pu
67 from ._polybase import ABCPolyBase
68
69 __all__ = [
70 'hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
71 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv',
72 'hermepow', 'hermeval', 'hermeder', 'hermeint', 'herme2poly',
73 'poly2herme', 'hermefromroots', 'hermevander', 'hermefit', 'hermetrim',
74 'hermeroots', 'HermiteE', 'hermeval2d', 'hermeval3d', 'hermegrid2d',
75 'hermegrid3d', 'hermevander2d', 'hermevander3d', 'hermecompanion',
76 'hermegauss', 'hermeweight']
77
78 hermetrim = pu.trimcoef
79
80
81 def poly2herme(pol):
82 """
83 poly2herme(pol)
84
85 Convert a polynomial to a Hermite series.
86
87 Convert an array representing the coefficients of a polynomial (relative
88 to the "standard" basis) ordered from lowest degree to highest, to an
89 array of the coefficients of the equivalent Hermite series, ordered
90 from lowest to highest degree.
91
92 Parameters
93 ----------
94 pol : array_like
95 1-D array containing the polynomial coefficients
96
97 Returns
98 -------
99 c : ndarray
100 1-D array containing the coefficients of the equivalent Hermite
101 series.
102
103 See Also
104 --------
105 herme2poly
106
107 Notes
108 -----
109 The easy way to do conversions between polynomial basis sets
110 is to use the convert method of a class instance.
111
112 Examples
113 --------
114 >>> from numpy.polynomial.hermite_e import poly2herme
115 >>> poly2herme(np.arange(4))
116 array([ 2., 10., 2., 3.])
117
118 """
119 [pol] = pu.as_series([pol])
120 deg = len(pol) - 1
121 res = 0
122 for i in range(deg, -1, -1):
123 res = hermeadd(hermemulx(res), pol[i])
124 return res
125
126
127 def herme2poly(c):
128 """
129 Convert a Hermite series to a polynomial.
130
131 Convert an array representing the coefficients of a Hermite series,
132 ordered from lowest degree to highest, to an array of the coefficients
133 of the equivalent polynomial (relative to the "standard" basis) ordered
134 from lowest to highest degree.
135
136 Parameters
137 ----------
138 c : array_like
139 1-D array containing the Hermite series coefficients, ordered
140 from lowest order term to highest.
141
142 Returns
143 -------
144 pol : ndarray
145 1-D array containing the coefficients of the equivalent polynomial
146 (relative to the "standard" basis) ordered from lowest order term
147 to highest.
148
149 See Also
150 --------
151 poly2herme
152
153 Notes
154 -----
155 The easy way to do conversions between polynomial basis sets
156 is to use the convert method of a class instance.
157
158 Examples
159 --------
160 >>> from numpy.polynomial.hermite_e import herme2poly
161 >>> herme2poly([ 2., 10., 2., 3.])
162 array([ 0., 1., 2., 3.])
163
164 """
165 from .polynomial import polyadd, polysub, polymulx
166
167 [c] = pu.as_series([c])
168 n = len(c)
169 if n == 1:
170 return c
171 if n == 2:
172 return c
173 else:
174 c0 = c[-2]
175 c1 = c[-1]
176 # i is the current degree of c1
177 for i in range(n - 1, 1, -1):
178 tmp = c0
179 c0 = polysub(c[i - 2], c1*(i - 1))
180 c1 = polyadd(tmp, polymulx(c1))
181 return polyadd(c0, polymulx(c1))
182
183 #
184 # These are constant arrays are of integer type so as to be compatible
185 # with the widest range of other types, such as Decimal.
186 #
187
188 # Hermite
189 hermedomain = np.array([-1, 1])
190
191 # Hermite coefficients representing zero.
192 hermezero = np.array([0])
193
194 # Hermite coefficients representing one.
195 hermeone = np.array([1])
196
197 # Hermite coefficients representing the identity x.
198 hermex = np.array([0, 1])
199
200
201 def hermeline(off, scl):
202 """
203 Hermite series whose graph is a straight line.
204
205
206
207 Parameters
208 ----------
209 off, scl : scalars
210 The specified line is given by ``off + scl*x``.
211
212 Returns
213 -------
214 y : ndarray
215 This module's representation of the Hermite series for
216 ``off + scl*x``.
217
218 See Also
219 --------
220 polyline, chebline
221
222 Examples
223 --------
224 >>> from numpy.polynomial.hermite_e import hermeline
225 >>> from numpy.polynomial.hermite_e import hermeline, hermeval
226 >>> hermeval(0,hermeline(3, 2))
227 3.0
228 >>> hermeval(1,hermeline(3, 2))
229 5.0
230
231 """
232 if scl != 0:
233 return np.array([off, scl])
234 else:
235 return np.array([off])
236
237
238 def hermefromroots(roots):
239 """
240 Generate a HermiteE series with given roots.
241
242 The function returns the coefficients of the polynomial
243
244 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
245
246 in HermiteE form, where the `r_n` are the roots specified in `roots`.
247 If a zero has multiplicity n, then it must appear in `roots` n times.
248 For instance, if 2 is a root of multiplicity three and 3 is a root of
249 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
250 roots can appear in any order.
251
252 If the returned coefficients are `c`, then
253
254 .. math:: p(x) = c_0 + c_1 * He_1(x) + ... + c_n * He_n(x)
255
256 The coefficient of the last term is not generally 1 for monic
257 polynomials in HermiteE form.
258
259 Parameters
260 ----------
261 roots : array_like
262 Sequence containing the roots.
263
264 Returns
265 -------
266 out : ndarray
267 1-D array of coefficients. If all roots are real then `out` is a
268 real array, if some of the roots are complex, then `out` is complex
269 even if all the coefficients in the result are real (see Examples
270 below).
271
272 See Also
273 --------
274 polyfromroots, legfromroots, lagfromroots, hermfromroots,
275 chebfromroots.
276
277 Examples
278 --------
279 >>> from numpy.polynomial.hermite_e import hermefromroots, hermeval
280 >>> coef = hermefromroots((-1, 0, 1))
281 >>> hermeval((-1, 0, 1), coef)
282 array([ 0., 0., 0.])
283 >>> coef = hermefromroots((-1j, 1j))
284 >>> hermeval((-1j, 1j), coef)
285 array([ 0.+0.j, 0.+0.j])
286
287 """
288 if len(roots) == 0:
289 return np.ones(1)
290 else:
291 [roots] = pu.as_series([roots], trim=False)
292 roots.sort()
293 p = [hermeline(-r, 1) for r in roots]
294 n = len(p)
295 while n > 1:
296 m, r = divmod(n, 2)
297 tmp = [hermemul(p[i], p[i+m]) for i in range(m)]
298 if r:
299 tmp[0] = hermemul(tmp[0], p[-1])
300 p = tmp
301 n = m
302 return p[0]
303
304
305 def hermeadd(c1, c2):
306 """
307 Add one Hermite series to another.
308
309 Returns the sum of two Hermite series `c1` + `c2`. The arguments
310 are sequences of coefficients ordered from lowest order term to
311 highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
312
313 Parameters
314 ----------
315 c1, c2 : array_like
316 1-D arrays of Hermite series coefficients ordered from low to
317 high.
318
319 Returns
320 -------
321 out : ndarray
322 Array representing the Hermite series of their sum.
323
324 See Also
325 --------
326 hermesub, hermemul, hermediv, hermepow
327
328 Notes
329 -----
330 Unlike multiplication, division, etc., the sum of two Hermite series
331 is a Hermite series (without having to "reproject" the result onto
332 the basis set) so addition, just like that of "standard" polynomials,
333 is simply "component-wise."
334
335 Examples
336 --------
337 >>> from numpy.polynomial.hermite_e import hermeadd
338 >>> hermeadd([1, 2, 3], [1, 2, 3, 4])
339 array([ 2., 4., 6., 4.])
340
341 """
342 # c1, c2 are trimmed copies
343 [c1, c2] = pu.as_series([c1, c2])
344 if len(c1) > len(c2):
345 c1[:c2.size] += c2
346 ret = c1
347 else:
348 c2[:c1.size] += c1
349 ret = c2
350 return pu.trimseq(ret)
351
352
353 def hermesub(c1, c2):
354 """
355 Subtract one Hermite series from another.
356
357 Returns the difference of two Hermite series `c1` - `c2`. The
358 sequences of coefficients are from lowest order term to highest, i.e.,
359 [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
360
361 Parameters
362 ----------
363 c1, c2 : array_like
364 1-D arrays of Hermite series coefficients ordered from low to
365 high.
366
367 Returns
368 -------
369 out : ndarray
370 Of Hermite series coefficients representing their difference.
371
372 See Also
373 --------
374 hermeadd, hermemul, hermediv, hermepow
375
376 Notes
377 -----
378 Unlike multiplication, division, etc., the difference of two Hermite
379 series is a Hermite series (without having to "reproject" the result
380 onto the basis set) so subtraction, just like that of "standard"
381 polynomials, is simply "component-wise."
382
383 Examples
384 --------
385 >>> from numpy.polynomial.hermite_e import hermesub
386 >>> hermesub([1, 2, 3, 4], [1, 2, 3])
387 array([ 0., 0., 0., 4.])
388
389 """
390 # c1, c2 are trimmed copies
391 [c1, c2] = pu.as_series([c1, c2])
392 if len(c1) > len(c2):
393 c1[:c2.size] -= c2
394 ret = c1
395 else:
396 c2 = -c2
397 c2[:c1.size] += c1
398 ret = c2
399 return pu.trimseq(ret)
400
401
402 def hermemulx(c):
403 """Multiply a Hermite series by x.
404
405 Multiply the Hermite series `c` by x, where x is the independent
406 variable.
407
408
409 Parameters
410 ----------
411 c : array_like
412 1-D array of Hermite series coefficients ordered from low to
413 high.
414
415 Returns
416 -------
417 out : ndarray
418 Array representing the result of the multiplication.
419
420 Notes
421 -----
422 The multiplication uses the recursion relationship for Hermite
423 polynomials in the form
424
425 .. math::
426
427 xP_i(x) = (P_{i + 1}(x) + iP_{i - 1}(x)))
428
429 Examples
430 --------
431 >>> from numpy.polynomial.hermite_e import hermemulx
432 >>> hermemulx([1, 2, 3])
433 array([ 2., 7., 2., 3.])
434
435 """
436 # c is a trimmed copy
437 [c] = pu.as_series([c])
438 # The zero series needs special treatment
439 if len(c) == 1 and c[0] == 0:
440 return c
441
442 prd = np.empty(len(c) + 1, dtype=c.dtype)
443 prd[0] = c[0]*0
444 prd[1] = c[0]
445 for i in range(1, len(c)):
446 prd[i + 1] = c[i]
447 prd[i - 1] += c[i]*i
448 return prd
449
450
451 def hermemul(c1, c2):
452 """
453 Multiply one Hermite series by another.
454
455 Returns the product of two Hermite series `c1` * `c2`. The arguments
456 are sequences of coefficients, from lowest order "term" to highest,
457 e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
458
459 Parameters
460 ----------
461 c1, c2 : array_like
462 1-D arrays of Hermite series coefficients ordered from low to
463 high.
464
465 Returns
466 -------
467 out : ndarray
468 Of Hermite series coefficients representing their product.
469
470 See Also
471 --------
472 hermeadd, hermesub, hermediv, hermepow
473
474 Notes
475 -----
476 In general, the (polynomial) product of two C-series results in terms
477 that are not in the Hermite polynomial basis set. Thus, to express
478 the product as a Hermite series, it is necessary to "reproject" the
479 product onto said basis set, which may produce "unintuitive" (but
480 correct) results; see Examples section below.
481
482 Examples
483 --------
484 >>> from numpy.polynomial.hermite_e import hermemul
485 >>> hermemul([1, 2, 3], [0, 1, 2])
486 array([ 14., 15., 28., 7., 6.])
487
488 """
489 # s1, s2 are trimmed copies
490 [c1, c2] = pu.as_series([c1, c2])
491
492 if len(c1) > len(c2):
493 c = c2
494 xs = c1
495 else:
496 c = c1
497 xs = c2
498
499 if len(c) == 1:
500 c0 = c[0]*xs
501 c1 = 0
502 elif len(c) == 2:
503 c0 = c[0]*xs
504 c1 = c[1]*xs
505 else:
506 nd = len(c)
507 c0 = c[-2]*xs
508 c1 = c[-1]*xs
509 for i in range(3, len(c) + 1):
510 tmp = c0
511 nd = nd - 1
512 c0 = hermesub(c[-i]*xs, c1*(nd - 1))
513 c1 = hermeadd(tmp, hermemulx(c1))
514 return hermeadd(c0, hermemulx(c1))
515
516
517 def hermediv(c1, c2):
518 """
519 Divide one Hermite series by another.
520
521 Returns the quotient-with-remainder of two Hermite series
522 `c1` / `c2`. The arguments are sequences of coefficients from lowest
523 order "term" to highest, e.g., [1,2,3] represents the series
524 ``P_0 + 2*P_1 + 3*P_2``.
525
526 Parameters
527 ----------
528 c1, c2 : array_like
529 1-D arrays of Hermite series coefficients ordered from low to
530 high.
531
532 Returns
533 -------
534 [quo, rem] : ndarrays
535 Of Hermite series coefficients representing the quotient and
536 remainder.
537
538 See Also
539 --------
540 hermeadd, hermesub, hermemul, hermepow
541
542 Notes
543 -----
544 In general, the (polynomial) division of one Hermite series by another
545 results in quotient and remainder terms that are not in the Hermite
546 polynomial basis set. Thus, to express these results as a Hermite
547 series, it is necessary to "reproject" the results onto the Hermite
548 basis set, which may produce "unintuitive" (but correct) results; see
549 Examples section below.
550
551 Examples
552 --------
553 >>> from numpy.polynomial.hermite_e import hermediv
554 >>> hermediv([ 14., 15., 28., 7., 6.], [0, 1, 2])
555 (array([ 1., 2., 3.]), array([ 0.]))
556 >>> hermediv([ 15., 17., 28., 7., 6.], [0, 1, 2])
557 (array([ 1., 2., 3.]), array([ 1., 2.]))
558
559 """
560 # c1, c2 are trimmed copies
561 [c1, c2] = pu.as_series([c1, c2])
562 if c2[-1] == 0:
563 raise ZeroDivisionError()
564
565 lc1 = len(c1)
566 lc2 = len(c2)
567 if lc1 < lc2:
568 return c1[:1]*0, c1
569 elif lc2 == 1:
570 return c1/c2[-1], c1[:1]*0
571 else:
572 quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
573 rem = c1
574 for i in range(lc1 - lc2, - 1, -1):
575 p = hermemul([0]*i + [1], c2)
576 q = rem[-1]/p[-1]
577 rem = rem[:-1] - q*p[:-1]
578 quo[i] = q
579 return quo, pu.trimseq(rem)
580
581
582 def hermepow(c, pow, maxpower=16):
583 """Raise a Hermite series to a power.
584
585 Returns the Hermite series `c` raised to the power `pow`. The
586 argument `c` is a sequence of coefficients ordered from low to high.
587 i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.``
588
589 Parameters
590 ----------
591 c : array_like
592 1-D array of Hermite series coefficients ordered from low to
593 high.
594 pow : integer
595 Power to which the series will be raised
596 maxpower : integer, optional
597 Maximum power allowed. This is mainly to limit growth of the series
598 to unmanageable size. Default is 16
599
600 Returns
601 -------
602 coef : ndarray
603 Hermite series of power.
604
605 See Also
606 --------
607 hermeadd, hermesub, hermemul, hermediv
608
609 Examples
610 --------
611 >>> from numpy.polynomial.hermite_e import hermepow
612 >>> hermepow([1, 2, 3], 2)
613 array([ 23., 28., 46., 12., 9.])
614
615 """
616 # c is a trimmed copy
617 [c] = pu.as_series([c])
618 power = int(pow)
619 if power != pow or power < 0:
620 raise ValueError("Power must be a non-negative integer.")
621 elif maxpower is not None and power > maxpower:
622 raise ValueError("Power is too large")
623 elif power == 0:
624 return np.array([1], dtype=c.dtype)
625 elif power == 1:
626 return c
627 else:
628 # This can be made more efficient by using powers of two
629 # in the usual way.
630 prd = c
631 for i in range(2, power + 1):
632 prd = hermemul(prd, c)
633 return prd
634
635
636 def hermeder(c, m=1, scl=1, axis=0):
637 """
638 Differentiate a Hermite_e series.
639
640 Returns the series coefficients `c` differentiated `m` times along
641 `axis`. At each iteration the result is multiplied by `scl` (the
642 scaling factor is for use in a linear change of variable). The argument
643 `c` is an array of coefficients from low to high degree along each
644 axis, e.g., [1,2,3] represents the series ``1*He_0 + 2*He_1 + 3*He_2``
645 while [[1,2],[1,2]] represents ``1*He_0(x)*He_0(y) + 1*He_1(x)*He_0(y)
646 + 2*He_0(x)*He_1(y) + 2*He_1(x)*He_1(y)`` if axis=0 is ``x`` and axis=1
647 is ``y``.
648
649 Parameters
650 ----------
651 c : array_like
652 Array of Hermite_e series coefficients. If `c` is multidimensional
653 the different axis correspond to different variables with the
654 degree in each axis given by the corresponding index.
655 m : int, optional
656 Number of derivatives taken, must be non-negative. (Default: 1)
657 scl : scalar, optional
658 Each differentiation is multiplied by `scl`. The end result is
659 multiplication by ``scl**m``. This is for use in a linear change of
660 variable. (Default: 1)
661 axis : int, optional
662 Axis over which the derivative is taken. (Default: 0).
663
664 .. versionadded:: 1.7.0
665
666 Returns
667 -------
668 der : ndarray
669 Hermite series of the derivative.
670
671 See Also
672 --------
673 hermeint
674
675 Notes
676 -----
677 In general, the result of differentiating a Hermite series does not
678 resemble the same operation on a power series. Thus the result of this
679 function may be "unintuitive," albeit correct; see Examples section
680 below.
681
682 Examples
683 --------
684 >>> from numpy.polynomial.hermite_e import hermeder
685 >>> hermeder([ 1., 1., 1., 1.])
686 array([ 1., 2., 3.])
687 >>> hermeder([-0.25, 1., 1./2., 1./3., 1./4 ], m=2)
688 array([ 1., 2., 3.])
689
690 """
691 c = np.array(c, ndmin=1, copy=1)
692 if c.dtype.char in '?bBhHiIlLqQpP':
693 c = c.astype(np.double)
694 cnt, iaxis = [int(t) for t in [m, axis]]
695
696 if cnt != m:
697 raise ValueError("The order of derivation must be integer")
698 if cnt < 0:
699 raise ValueError("The order of derivation must be non-negative")
700 if iaxis != axis:
701 raise ValueError("The axis must be integer")
702 if not -c.ndim <= iaxis < c.ndim:
703 raise ValueError("The axis is out of range")
704 if iaxis < 0:
705 iaxis += c.ndim
706
707 if cnt == 0:
708 return c
709
710 c = np.rollaxis(c, iaxis)
711 n = len(c)
712 if cnt >= n:
713 return c[:1]*0
714 else:
715 for i in range(cnt):
716 n = n - 1
717 c *= scl
718 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
719 for j in range(n, 0, -1):
720 der[j - 1] = j*c[j]
721 c = der
722 c = np.rollaxis(c, 0, iaxis + 1)
723 return c
724
725
726 def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
727 """
728 Integrate a Hermite_e series.
729
730 Returns the Hermite_e series coefficients `c` integrated `m` times from
731 `lbnd` along `axis`. At each iteration the resulting series is
732 **multiplied** by `scl` and an integration constant, `k`, is added.
733 The scaling factor is for use in a linear change of variable. ("Buyer
734 beware": note that, depending on what one is doing, one may want `scl`
735 to be the reciprocal of what one might expect; for more information,
736 see the Notes section below.) The argument `c` is an array of
737 coefficients from low to high degree along each axis, e.g., [1,2,3]
738 represents the series ``H_0 + 2*H_1 + 3*H_2`` while [[1,2],[1,2]]
739 represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) + 2*H_0(x)*H_1(y) +
740 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
741
742 Parameters
743 ----------
744 c : array_like
745 Array of Hermite_e series coefficients. If c is multidimensional
746 the different axis correspond to different variables with the
747 degree in each axis given by the corresponding index.
748 m : int, optional
749 Order of integration, must be positive. (Default: 1)
750 k : {[], list, scalar}, optional
751 Integration constant(s). The value of the first integral at
752 ``lbnd`` is the first value in the list, the value of the second
753 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
754 default), all constants are set to zero. If ``m == 1``, a single
755 scalar can be given instead of a list.
756 lbnd : scalar, optional
757 The lower bound of the integral. (Default: 0)
758 scl : scalar, optional
759 Following each integration the result is *multiplied* by `scl`
760 before the integration constant is added. (Default: 1)
761 axis : int, optional
762 Axis over which the integral is taken. (Default: 0).
763
764 .. versionadded:: 1.7.0
765
766 Returns
767 -------
768 S : ndarray
769 Hermite_e series coefficients of the integral.
770
771 Raises
772 ------
773 ValueError
774 If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or
775 ``np.isscalar(scl) == False``.
776
777 See Also
778 --------
779 hermeder
780
781 Notes
782 -----
783 Note that the result of each integration is *multiplied* by `scl`.
784 Why is this important to note? Say one is making a linear change of
785 variable :math:`u = ax + b` in an integral relative to `x`. Then
786 .. math::`dx = du/a`, so one will need to set `scl` equal to
787 :math:`1/a` - perhaps not what one would have first thought.
788
789 Also note that, in general, the result of integrating a C-series needs
790 to be "reprojected" onto the C-series basis set. Thus, typically,
791 the result of this function is "unintuitive," albeit correct; see
792 Examples section below.
793
794 Examples
795 --------
796 >>> from numpy.polynomial.hermite_e import hermeint
797 >>> hermeint([1, 2, 3]) # integrate once, value 0 at 0.
798 array([ 1., 1., 1., 1.])
799 >>> hermeint([1, 2, 3], m=2) # integrate twice, value & deriv 0 at 0
800 array([-0.25 , 1. , 0.5 , 0.33333333, 0.25 ])
801 >>> hermeint([1, 2, 3], k=1) # integrate once, value 1 at 0.
802 array([ 2., 1., 1., 1.])
803 >>> hermeint([1, 2, 3], lbnd=-1) # integrate once, value 0 at -1
804 array([-1., 1., 1., 1.])
805 >>> hermeint([1, 2, 3], m=2, k=[1, 2], lbnd=-1)
806 array([ 1.83333333, 0. , 0.5 , 0.33333333, 0.25 ])
807
808 """
809 c = np.array(c, ndmin=1, copy=1)
810 if c.dtype.char in '?bBhHiIlLqQpP':
811 c = c.astype(np.double)
812 if not np.iterable(k):
813 k = [k]
814 cnt, iaxis = [int(t) for t in [m, axis]]
815
816 if cnt != m:
817 raise ValueError("The order of integration must be integer")
818 if cnt < 0:
819 raise ValueError("The order of integration must be non-negative")
820 if len(k) > cnt:
821 raise ValueError("Too many integration constants")
822 if iaxis != axis:
823 raise ValueError("The axis must be integer")
824 if not -c.ndim <= iaxis < c.ndim:
825 raise ValueError("The axis is out of range")
826 if iaxis < 0:
827 iaxis += c.ndim
828
829 if cnt == 0:
830 return c
831
832 c = np.rollaxis(c, iaxis)
833 k = list(k) + [0]*(cnt - len(k))
834 for i in range(cnt):
835 n = len(c)
836 c *= scl
837 if n == 1 and np.all(c[0] == 0):
838 c[0] += k[i]
839 else:
840 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
841 tmp[0] = c[0]*0
842 tmp[1] = c[0]
843 for j in range(1, n):
844 tmp[j + 1] = c[j]/(j + 1)
845 tmp[0] += k[i] - hermeval(lbnd, tmp)
846 c = tmp
847 c = np.rollaxis(c, 0, iaxis + 1)
848 return c
849
850
851 def hermeval(x, c, tensor=True):
852 """
853 Evaluate an HermiteE series at points x.
854
855 If `c` is of length `n + 1`, this function returns the value:
856
857 .. math:: p(x) = c_0 * He_0(x) + c_1 * He_1(x) + ... + c_n * He_n(x)
858
859 The parameter `x` is converted to an array only if it is a tuple or a
860 list, otherwise it is treated as a scalar. In either case, either `x`
861 or its elements must support multiplication and addition both with
862 themselves and with the elements of `c`.
863
864 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
865 `c` is multidimensional, then the shape of the result depends on the
866 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
867 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
868 scalars have shape (,).
869
870 Trailing zeros in the coefficients will be used in the evaluation, so
871 they should be avoided if efficiency is a concern.
872
873 Parameters
874 ----------
875 x : array_like, compatible object
876 If `x` is a list or tuple, it is converted to an ndarray, otherwise
877 it is left unchanged and treated as a scalar. In either case, `x`
878 or its elements must support addition and multiplication with
879 with themselves and with the elements of `c`.
880 c : array_like
881 Array of coefficients ordered so that the coefficients for terms of
882 degree n are contained in c[n]. If `c` is multidimensional the
883 remaining indices enumerate multiple polynomials. In the two
884 dimensional case the coefficients may be thought of as stored in
885 the columns of `c`.
886 tensor : boolean, optional
887 If True, the shape of the coefficient array is extended with ones
888 on the right, one for each dimension of `x`. Scalars have dimension 0
889 for this action. The result is that every column of coefficients in
890 `c` is evaluated for every element of `x`. If False, `x` is broadcast
891 over the columns of `c` for the evaluation. This keyword is useful
892 when `c` is multidimensional. The default value is True.
893
894 .. versionadded:: 1.7.0
895
896 Returns
897 -------
898 values : ndarray, algebra_like
899 The shape of the return value is described above.
900
901 See Also
902 --------
903 hermeval2d, hermegrid2d, hermeval3d, hermegrid3d
904
905 Notes
906 -----
907 The evaluation uses Clenshaw recursion, aka synthetic division.
908
909 Examples
910 --------
911 >>> from numpy.polynomial.hermite_e import hermeval
912 >>> coef = [1,2,3]
913 >>> hermeval(1, coef)
914 3.0
915 >>> hermeval([[1,2],[3,4]], coef)
916 array([[ 3., 14.],
917 [ 31., 54.]])
918
919 """
920 c = np.array(c, ndmin=1, copy=0)
921 if c.dtype.char in '?bBhHiIlLqQpP':
922 c = c.astype(np.double)
923 if isinstance(x, (tuple, list)):
924 x = np.asarray(x)
925 if isinstance(x, np.ndarray) and tensor:
926 c = c.reshape(c.shape + (1,)*x.ndim)
927
928 if len(c) == 1:
929 c0 = c[0]
930 c1 = 0
931 elif len(c) == 2:
932 c0 = c[0]
933 c1 = c[1]
934 else:
935 nd = len(c)
936 c0 = c[-2]
937 c1 = c[-1]
938 for i in range(3, len(c) + 1):
939 tmp = c0
940 nd = nd - 1
941 c0 = c[-i] - c1*(nd - 1)
942 c1 = tmp + c1*x
943 return c0 + c1*x
944
945
946 def hermeval2d(x, y, c):
947 """
948 Evaluate a 2-D HermiteE series at points (x, y).
949
950 This function returns the values:
951
952 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * He_i(x) * He_j(y)
953
954 The parameters `x` and `y` are converted to arrays only if they are
955 tuples or a lists, otherwise they are treated as a scalars and they
956 must have the same shape after conversion. In either case, either `x`
957 and `y` or their elements must support multiplication and addition both
958 with themselves and with the elements of `c`.
959
960 If `c` is a 1-D array a one is implicitly appended to its shape to make
961 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
962
963 Parameters
964 ----------
965 x, y : array_like, compatible objects
966 The two dimensional series is evaluated at the points `(x, y)`,
967 where `x` and `y` must have the same shape. If `x` or `y` is a list
968 or tuple, it is first converted to an ndarray, otherwise it is left
969 unchanged and if it isn't an ndarray it is treated as a scalar.
970 c : array_like
971 Array of coefficients ordered so that the coefficient of the term
972 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
973 dimension greater than two the remaining indices enumerate multiple
974 sets of coefficients.
975
976 Returns
977 -------
978 values : ndarray, compatible object
979 The values of the two dimensional polynomial at points formed with
980 pairs of corresponding values from `x` and `y`.
981
982 See Also
983 --------
984 hermeval, hermegrid2d, hermeval3d, hermegrid3d
985
986 Notes
987 -----
988
989 .. versionadded::1.7.0
990
991 """
992 try:
993 x, y = np.array((x, y), copy=0)
994 except:
995 raise ValueError('x, y are incompatible')
996
997 c = hermeval(x, c)
998 c = hermeval(y, c, tensor=False)
999 return c
1000
1001
1002 def hermegrid2d(x, y, c):
1003 """
1004 Evaluate a 2-D HermiteE series on the Cartesian product of x and y.
1005
1006 This function returns the values:
1007
1008 .. math:: p(a,b) = \sum_{i,j} c_{i,j} * H_i(a) * H_j(b)
1009
1010 where the points `(a, b)` consist of all pairs formed by taking
1011 `a` from `x` and `b` from `y`. The resulting points form a grid with
1012 `x` in the first dimension and `y` in the second.
1013
1014 The parameters `x` and `y` are converted to arrays only if they are
1015 tuples or a lists, otherwise they are treated as a scalars. In either
1016 case, either `x` and `y` or their elements must support multiplication
1017 and addition both with themselves and with the elements of `c`.
1018
1019 If `c` has fewer than two dimensions, ones are implicitly appended to
1020 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
1021 x.shape.
1022
1023 Parameters
1024 ----------
1025 x, y : array_like, compatible objects
1026 The two dimensional series is evaluated at the points in the
1027 Cartesian product of `x` and `y`. If `x` or `y` is a list or
1028 tuple, it is first converted to an ndarray, otherwise it is left
1029 unchanged and, if it isn't an ndarray, it is treated as a scalar.
1030 c : array_like
1031 Array of coefficients ordered so that the coefficients for terms of
1032 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1033 greater than two the remaining indices enumerate multiple sets of
1034 coefficients.
1035
1036 Returns
1037 -------
1038 values : ndarray, compatible object
1039 The values of the two dimensional polynomial at points in the Cartesian
1040 product of `x` and `y`.
1041
1042 See Also
1043 --------
1044 hermeval, hermeval2d, hermeval3d, hermegrid3d
1045
1046 Notes
1047 -----
1048
1049 .. versionadded::1.7.0
1050
1051 """
1052 c = hermeval(x, c)
1053 c = hermeval(y, c)
1054 return c
1055
1056
1057 def hermeval3d(x, y, z, c):
1058 """
1059 Evaluate a 3-D Hermite_e series at points (x, y, z).
1060
1061 This function returns the values:
1062
1063 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * He_i(x) * He_j(y) * He_k(z)
1064
1065 The parameters `x`, `y`, and `z` are converted to arrays only if
1066 they are tuples or a lists, otherwise they are treated as a scalars and
1067 they must have the same shape after conversion. In either case, either
1068 `x`, `y`, and `z` or their elements must support multiplication and
1069 addition both with themselves and with the elements of `c`.
1070
1071 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1072 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1073 x.shape.
1074
1075 Parameters
1076 ----------
1077 x, y, z : array_like, compatible object
1078 The three dimensional series is evaluated at the points
1079 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
1080 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1081 to an ndarray, otherwise it is left unchanged and if it isn't an
1082 ndarray it is treated as a scalar.
1083 c : array_like
1084 Array of coefficients ordered so that the coefficient of the term of
1085 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1086 greater than 3 the remaining indices enumerate multiple sets of
1087 coefficients.
1088
1089 Returns
1090 -------
1091 values : ndarray, compatible object
1092 The values of the multidimensional polynomial on points formed with
1093 triples of corresponding values from `x`, `y`, and `z`.
1094
1095 See Also
1096 --------
1097 hermeval, hermeval2d, hermegrid2d, hermegrid3d
1098
1099 Notes
1100 -----
1101
1102 .. versionadded::1.7.0
1103
1104 """
1105 try:
1106 x, y, z = np.array((x, y, z), copy=0)
1107 except:
1108 raise ValueError('x, y, z are incompatible')
1109
1110 c = hermeval(x, c)
1111 c = hermeval(y, c, tensor=False)
1112 c = hermeval(z, c, tensor=False)
1113 return c
1114
1115
1116 def hermegrid3d(x, y, z, c):
1117 """
1118 Evaluate a 3-D HermiteE series on the Cartesian product of x, y, and z.
1119
1120 This function returns the values:
1121
1122 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * He_i(a) * He_j(b) * He_k(c)
1123
1124 where the points `(a, b, c)` consist of all triples formed by taking
1125 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1126 a grid with `x` in the first dimension, `y` in the second, and `z` in
1127 the third.
1128
1129 The parameters `x`, `y`, and `z` are converted to arrays only if they
1130 are tuples or a lists, otherwise they are treated as a scalars. In
1131 either case, either `x`, `y`, and `z` or their elements must support
1132 multiplication and addition both with themselves and with the elements
1133 of `c`.
1134
1135 If `c` has fewer than three dimensions, ones are implicitly appended to
1136 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1137 x.shape + y.shape + z.shape.
1138
1139 Parameters
1140 ----------
1141 x, y, z : array_like, compatible objects
1142 The three dimensional series is evaluated at the points in the
1143 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
1144 list or tuple, it is first converted to an ndarray, otherwise it is
1145 left unchanged and, if it isn't an ndarray, it is treated as a
1146 scalar.
1147 c : array_like
1148 Array of coefficients ordered so that the coefficients for terms of
1149 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1150 greater than two the remaining indices enumerate multiple sets of
1151 coefficients.
1152
1153 Returns
1154 -------
1155 values : ndarray, compatible object
1156 The values of the two dimensional polynomial at points in the Cartesian
1157 product of `x` and `y`.
1158
1159 See Also
1160 --------
1161 hermeval, hermeval2d, hermegrid2d, hermeval3d
1162
1163 Notes
1164 -----
1165
1166 .. versionadded::1.7.0
1167
1168 """
1169 c = hermeval(x, c)
1170 c = hermeval(y, c)
1171 c = hermeval(z, c)
1172 return c
1173
1174
1175 def hermevander(x, deg):
1176 """Pseudo-Vandermonde matrix of given degree.
1177
1178 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1179 `x`. The pseudo-Vandermonde matrix is defined by
1180
1181 .. math:: V[..., i] = He_i(x),
1182
1183 where `0 <= i <= deg`. The leading indices of `V` index the elements of
1184 `x` and the last index is the degree of the HermiteE polynomial.
1185
1186 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
1187 array ``V = hermevander(x, n)``, then ``np.dot(V, c)`` and
1188 ``hermeval(x, c)`` are the same up to roundoff. This equivalence is
1189 useful both for least squares fitting and for the evaluation of a large
1190 number of HermiteE series of the same degree and sample points.
1191
1192 Parameters
1193 ----------
1194 x : array_like
1195 Array of points. The dtype is converted to float64 or complex128
1196 depending on whether any of the elements are complex. If `x` is
1197 scalar it is converted to a 1-D array.
1198 deg : int
1199 Degree of the resulting matrix.
1200
1201 Returns
1202 -------
1203 vander : ndarray
1204 The pseudo-Vandermonde matrix. The shape of the returned matrix is
1205 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1206 corresponding HermiteE polynomial. The dtype will be the same as
1207 the converted `x`.
1208
1209 Examples
1210 --------
1211 >>> from numpy.polynomial.hermite_e import hermevander
1212 >>> x = np.array([-1, 0, 1])
1213 >>> hermevander(x, 3)
1214 array([[ 1., -1., 0., 2.],
1215 [ 1., 0., -1., -0.],
1216 [ 1., 1., 0., -2.]])
1217
1218 """
1219 ideg = int(deg)
1220 if ideg != deg:
1221 raise ValueError("deg must be integer")
1222 if ideg < 0:
1223 raise ValueError("deg must be non-negative")
1224
1225 x = np.array(x, copy=0, ndmin=1) + 0.0
1226 dims = (ideg + 1,) + x.shape
1227 dtyp = x.dtype
1228 v = np.empty(dims, dtype=dtyp)
1229 v[0] = x*0 + 1
1230 if ideg > 0:
1231 v[1] = x
1232 for i in range(2, ideg + 1):
1233 v[i] = (v[i-1]*x - v[i-2]*(i - 1))
1234 return np.rollaxis(v, 0, v.ndim)
1235
1236
1237 def hermevander2d(x, y, deg):
1238 """Pseudo-Vandermonde matrix of given degrees.
1239
1240 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1241 points `(x, y)`. The pseudo-Vandermonde matrix is defined by
1242
1243 .. math:: V[..., deg[1]*i + j] = He_i(x) * He_j(y),
1244
1245 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
1246 `V` index the points `(x, y)` and the last index encodes the degrees of
1247 the HermiteE polynomials.
1248
1249 If ``V = hermevander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1250 correspond to the elements of a 2-D coefficient array `c` of shape
1251 (xdeg + 1, ydeg + 1) in the order
1252
1253 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
1254
1255 and ``np.dot(V, c.flat)`` and ``hermeval2d(x, y, c)`` will be the same
1256 up to roundoff. This equivalence is useful both for least squares
1257 fitting and for the evaluation of a large number of 2-D HermiteE
1258 series of the same degrees and sample points.
1259
1260 Parameters
1261 ----------
1262 x, y : array_like
1263 Arrays of point coordinates, all of the same shape. The dtypes
1264 will be converted to either float64 or complex128 depending on
1265 whether any of the elements are complex. Scalars are converted to
1266 1-D arrays.
1267 deg : list of ints
1268 List of maximum degrees of the form [x_deg, y_deg].
1269
1270 Returns
1271 -------
1272 vander2d : ndarray
1273 The shape of the returned matrix is ``x.shape + (order,)``, where
1274 :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same
1275 as the converted `x` and `y`.
1276
1277 See Also
1278 --------
1279 hermevander, hermevander3d. hermeval2d, hermeval3d
1280
1281 Notes
1282 -----
1283
1284 .. versionadded::1.7.0
1285
1286 """
1287 ideg = [int(d) for d in deg]
1288 is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
1289 if is_valid != [1, 1]:
1290 raise ValueError("degrees must be non-negative integers")
1291 degx, degy = ideg
1292 x, y = np.array((x, y), copy=0) + 0.0
1293
1294 vx = hermevander(x, degx)
1295 vy = hermevander(y, degy)
1296 v = vx[..., None]*vy[..., None,:]
1297 return v.reshape(v.shape[:-2] + (-1,))
1298
1299
1300 def hermevander3d(x, y, z, deg):
1301 """Pseudo-Vandermonde matrix of given degrees.
1302
1303 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1304 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
1305 then Hehe pseudo-Vandermonde matrix is defined by
1306
1307 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = He_i(x)*He_j(y)*He_k(z),
1308
1309 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
1310 indices of `V` index the points `(x, y, z)` and the last index encodes
1311 the degrees of the HermiteE polynomials.
1312
1313 If ``V = hermevander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
1314 of `V` correspond to the elements of a 3-D coefficient array `c` of
1315 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
1316
1317 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
1318
1319 and ``np.dot(V, c.flat)`` and ``hermeval3d(x, y, z, c)`` will be the
1320 same up to roundoff. This equivalence is useful both for least squares
1321 fitting and for the evaluation of a large number of 3-D HermiteE
1322 series of the same degrees and sample points.
1323
1324 Parameters
1325 ----------
1326 x, y, z : array_like
1327 Arrays of point coordinates, all of the same shape. The dtypes will
1328 be converted to either float64 or complex128 depending on whether
1329 any of the elements are complex. Scalars are converted to 1-D
1330 arrays.
1331 deg : list of ints
1332 List of maximum degrees of the form [x_deg, y_deg, z_deg].
1333
1334 Returns
1335 -------
1336 vander3d : ndarray
1337 The shape of the returned matrix is ``x.shape + (order,)``, where
1338 :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will
1339 be the same as the converted `x`, `y`, and `z`.
1340
1341 See Also
1342 --------
1343 hermevander, hermevander3d. hermeval2d, hermeval3d
1344
1345 Notes
1346 -----
1347
1348 .. versionadded::1.7.0
1349
1350 """
1351 ideg = [int(d) for d in deg]
1352 is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
1353 if is_valid != [1, 1, 1]:
1354 raise ValueError("degrees must be non-negative integers")
1355 degx, degy, degz = ideg
1356 x, y, z = np.array((x, y, z), copy=0) + 0.0
1357
1358 vx = hermevander(x, degx)
1359 vy = hermevander(y, degy)
1360 vz = hermevander(z, degz)
1361 v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
1362 return v.reshape(v.shape[:-3] + (-1,))
1363
1364
1365 def hermefit(x, y, deg, rcond=None, full=False, w=None):
1366 """
1367 Least squares fit of Hermite series to data.
1368
1369 Return the coefficients of a HermiteE series of degree `deg` that is
1370 the least squares fit to the data values `y` given at points `x`. If
1371 `y` is 1-D the returned coefficients will also be 1-D. If `y` is 2-D
1372 multiple fits are done, one for each column of `y`, and the resulting
1373 coefficients are stored in the corresponding columns of a 2-D return.
1374 The fitted polynomial(s) are in the form
1375
1376 .. math:: p(x) = c_0 + c_1 * He_1(x) + ... + c_n * He_n(x),
1377
1378 where `n` is `deg`.
1379
1380 Parameters
1381 ----------
1382 x : array_like, shape (M,)
1383 x-coordinates of the M sample points ``(x[i], y[i])``.
1384 y : array_like, shape (M,) or (M, K)
1385 y-coordinates of the sample points. Several data sets of sample
1386 points sharing the same x-coordinates can be fitted at once by
1387 passing in a 2D-array that contains one dataset per column.
1388 deg : int
1389 Degree of the fitting polynomial
1390 rcond : float, optional
1391 Relative condition number of the fit. Singular values smaller than
1392 this relative to the largest singular value will be ignored. The
1393 default value is len(x)*eps, where eps is the relative precision of
1394 the float type, about 2e-16 in most cases.
1395 full : bool, optional
1396 Switch determining nature of return value. When it is False (the
1397 default) just the coefficients are returned, when True diagnostic
1398 information from the singular value decomposition is also returned.
1399 w : array_like, shape (`M`,), optional
1400 Weights. If not None, the contribution of each point
1401 ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
1402 weights are chosen so that the errors of the products ``w[i]*y[i]``
1403 all have the same variance. The default value is None.
1404
1405 Returns
1406 -------
1407 coef : ndarray, shape (M,) or (M, K)
1408 Hermite coefficients ordered from low to high. If `y` was 2-D,
1409 the coefficients for the data in column k of `y` are in column
1410 `k`.
1411
1412 [residuals, rank, singular_values, rcond] : list
1413 These values are only returned if `full` = True
1414
1415 resid -- sum of squared residuals of the least squares fit
1416 rank -- the numerical rank of the scaled Vandermonde matrix
1417 sv -- singular values of the scaled Vandermonde matrix
1418 rcond -- value of `rcond`.
1419
1420 For more details, see `linalg.lstsq`.
1421
1422 Warns
1423 -----
1424 RankWarning
1425 The rank of the coefficient matrix in the least-squares fit is
1426 deficient. The warning is only raised if `full` = False. The
1427 warnings can be turned off by
1428
1429 >>> import warnings
1430 >>> warnings.simplefilter('ignore', RankWarning)
1431
1432 See Also
1433 --------
1434 chebfit, legfit, polyfit, hermfit, polyfit
1435 hermeval : Evaluates a Hermite series.
1436 hermevander : pseudo Vandermonde matrix of Hermite series.
1437 hermeweight : HermiteE weight function.
1438 linalg.lstsq : Computes a least-squares fit from the matrix.
1439 scipy.interpolate.UnivariateSpline : Computes spline fits.
1440
1441 Notes
1442 -----
1443 The solution is the coefficients of the HermiteE series `p` that
1444 minimizes the sum of the weighted squared errors
1445
1446 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
1447
1448 where the :math:`w_j` are the weights. This problem is solved by
1449 setting up the (typically) overdetermined matrix equation
1450
1451 .. math:: V(x) * c = w * y,
1452
1453 where `V` is the pseudo Vandermonde matrix of `x`, the elements of `c`
1454 are the coefficients to be solved for, and the elements of `y` are the
1455 observed values. This equation is then solved using the singular value
1456 decomposition of `V`.
1457
1458 If some of the singular values of `V` are so small that they are
1459 neglected, then a `RankWarning` will be issued. This means that the
1460 coefficient values may be poorly determined. Using a lower order fit
1461 will usually get rid of the warning. The `rcond` parameter can also be
1462 set to a value smaller than its default, but the resulting fit may be
1463 spurious and have large contributions from roundoff error.
1464
1465 Fits using HermiteE series are probably most useful when the data can
1466 be approximated by ``sqrt(w(x)) * p(x)``, where `w(x)` is the HermiteE
1467 weight. In that case the weight ``sqrt(w(x[i])`` should be used
1468 together with data values ``y[i]/sqrt(w(x[i])``. The weight function is
1469 available as `hermeweight`.
1470
1471 References
1472 ----------
1473 .. [1] Wikipedia, "Curve fitting",
1474 http://en.wikipedia.org/wiki/Curve_fitting
1475
1476 Examples
1477 --------
1478 >>> from numpy.polynomial.hermite_e import hermefik, hermeval
1479 >>> x = np.linspace(-10, 10)
1480 >>> err = np.random.randn(len(x))/10
1481 >>> y = hermeval(x, [1, 2, 3]) + err
1482 >>> hermefit(x, y, 2)
1483 array([ 1.01690445, 1.99951418, 2.99948696])
1484
1485 """
1486 order = int(deg) + 1
1487 x = np.asarray(x) + 0.0
1488 y = np.asarray(y) + 0.0
1489
1490 # check arguments.
1491 if deg < 0:
1492 raise ValueError("expected deg >= 0")
1493 if x.ndim != 1:
1494 raise TypeError("expected 1D vector for x")
1495 if x.size == 0:
1496 raise TypeError("expected non-empty vector for x")
1497 if y.ndim < 1 or y.ndim > 2:
1498 raise TypeError("expected 1D or 2D array for y")
1499 if len(x) != len(y):
1500 raise TypeError("expected x and y to have same length")
1501
1502 # set up the least squares matrices in transposed form
1503 lhs = hermevander(x, deg).T
1504 rhs = y.T
1505 if w is not None:
1506 w = np.asarray(w) + 0.0
1507 if w.ndim != 1:
1508 raise TypeError("expected 1D vector for w")
1509 if len(x) != len(w):
1510 raise TypeError("expected x and w to have same length")
1511 # apply weights. Don't use inplace operations as they
1512 # can cause problems with NA.
1513 lhs = lhs * w
1514 rhs = rhs * w
1515
1516 # set rcond
1517 if rcond is None:
1518 rcond = len(x)*np.finfo(x.dtype).eps
1519
1520 # Determine the norms of the design matrix columns.
1521 if issubclass(lhs.dtype.type, np.complexfloating):
1522 scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
1523 else:
1524 scl = np.sqrt(np.square(lhs).sum(1))
1525 scl[scl == 0] = 1
1526
1527 # Solve the least squares problem.
1528 c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
1529 c = (c.T/scl).T
1530
1531 # warn on rank reduction
1532 if rank != order and not full:
1533 msg = "The fit may be poorly conditioned"
1534 warnings.warn(msg, pu.RankWarning)
1535
1536 if full:
1537 return c, [resids, rank, s, rcond]
1538 else:
1539 return c
1540
1541
1542 def hermecompanion(c):
1543 """
1544 Return the scaled companion matrix of c.
1545
1546 The basis polynomials are scaled so that the companion matrix is
1547 symmetric when `c` is an HermiteE basis polynomial. This provides
1548 better eigenvalue estimates than the unscaled case and for basis
1549 polynomials the eigenvalues are guaranteed to be real if
1550 `numpy.linalg.eigvalsh` is used to obtain them.
1551
1552 Parameters
1553 ----------
1554 c : array_like
1555 1-D array of HermiteE series coefficients ordered from low to high
1556 degree.
1557
1558 Returns
1559 -------
1560 mat : ndarray
1561 Scaled companion matrix of dimensions (deg, deg).
1562
1563 Notes
1564 -----
1565
1566 .. versionadded::1.7.0
1567
1568 """
1569 # c is a trimmed copy
1570 [c] = pu.as_series([c])
1571 if len(c) < 2:
1572 raise ValueError('Series must have maximum degree of at least 1.')
1573 if len(c) == 2:
1574 return np.array([[-c[0]/c[1]]])
1575
1576 n = len(c) - 1
1577 mat = np.zeros((n, n), dtype=c.dtype)
1578 scl = np.hstack((1., np.sqrt(np.arange(1, n))))
1579 scl = np.multiply.accumulate(scl)
1580 top = mat.reshape(-1)[1::n+1]
1581 bot = mat.reshape(-1)[n::n+1]
1582 top[...] = np.sqrt(np.arange(1, n))
1583 bot[...] = top
1584 mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])
1585 return mat
1586
1587
1588 def hermeroots(c):
1589 """
1590 Compute the roots of a HermiteE series.
1591
1592 Return the roots (a.k.a. "zeros") of the polynomial
1593
1594 .. math:: p(x) = \\sum_i c[i] * He_i(x).
1595
1596 Parameters
1597 ----------
1598 c : 1-D array_like
1599 1-D array of coefficients.
1600
1601 Returns
1602 -------
1603 out : ndarray
1604 Array of the roots of the series. If all the roots are real,
1605 then `out` is also real, otherwise it is complex.
1606
1607 See Also
1608 --------
1609 polyroots, legroots, lagroots, hermroots, chebroots
1610
1611 Notes
1612 -----
1613 The root estimates are obtained as the eigenvalues of the companion
1614 matrix, Roots far from the origin of the complex plane may have large
1615 errors due to the numerical instability of the series for such
1616 values. Roots with multiplicity greater than 1 will also show larger
1617 errors as the value of the series near such points is relatively
1618 insensitive to errors in the roots. Isolated roots near the origin can
1619 be improved by a few iterations of Newton's method.
1620
1621 The HermiteE series basis polynomials aren't powers of `x` so the
1622 results of this function may seem unintuitive.
1623
1624 Examples
1625 --------
1626 >>> from numpy.polynomial.hermite_e import hermeroots, hermefromroots
1627 >>> coef = hermefromroots([-1, 0, 1])
1628 >>> coef
1629 array([ 0., 2., 0., 1.])
1630 >>> hermeroots(coef)
1631 array([-1., 0., 1.])
1632
1633 """
1634 # c is a trimmed copy
1635 [c] = pu.as_series([c])
1636 if len(c) <= 1:
1637 return np.array([], dtype=c.dtype)
1638 if len(c) == 2:
1639 return np.array([-c[0]/c[1]])
1640
1641 m = hermecompanion(c)
1642 r = la.eigvals(m)
1643 r.sort()
1644 return r
1645
1646
1647 def hermegauss(deg):
1648 """
1649 Gauss-HermiteE quadrature.
1650
1651 Computes the sample points and weights for Gauss-HermiteE quadrature.
1652 These sample points and weights will correctly integrate polynomials of
1653 degree :math:`2*deg - 1` or less over the interval :math:`[-\inf, \inf]`
1654 with the weight function :math:`f(x) = \exp(-x^2/2)`.
1655
1656 Parameters
1657 ----------
1658 deg : int
1659 Number of sample points and weights. It must be >= 1.
1660
1661 Returns
1662 -------
1663 x : ndarray
1664 1-D ndarray containing the sample points.
1665 y : ndarray
1666 1-D ndarray containing the weights.
1667
1668 Notes
1669 -----
1670
1671 .. versionadded::1.7.0
1672
1673 The results have only been tested up to degree 100, higher degrees may
1674 be problematic. The weights are determined by using the fact that
1675
1676 .. math:: w_k = c / (He'_n(x_k) * He_{n-1}(x_k))
1677
1678 where :math:`c` is a constant independent of :math:`k` and :math:`x_k`
1679 is the k'th root of :math:`He_n`, and then scaling the results to get
1680 the right value when integrating 1.
1681
1682 """
1683 ideg = int(deg)
1684 if ideg != deg or ideg < 1:
1685 raise ValueError("deg must be a non-negative integer")
1686
1687 # first approximation of roots. We use the fact that the companion
1688 # matrix is symmetric in this case in order to obtain better zeros.
1689 c = np.array([0]*deg + [1])
1690 m = hermecompanion(c)
1691 x = la.eigvals(m)
1692 x.sort()
1693
1694 # improve roots by one application of Newton
1695 dy = hermeval(x, c)
1696 df = hermeval(x, hermeder(c))
1697 x -= dy/df
1698
1699 # compute the weights. We scale the factor to avoid possible numerical
1700 # overflow.
1701 fm = hermeval(x, c[1:])
1702 fm /= np.abs(fm).max()
1703 df /= np.abs(df).max()
1704 w = 1/(fm * df)
1705
1706 # for Hermite_e we can also symmetrize
1707 w = (w + w[::-1])/2
1708 x = (x - x[::-1])/2
1709
1710 # scale w to get the right value
1711 w *= np.sqrt(2*np.pi) / w.sum()
1712
1713 return x, w
1714
1715
1716 def hermeweight(x):
1717 """Weight function of the Hermite_e polynomials.
1718
1719 The weight function is :math:`\exp(-x^2/2)` and the interval of
1720 integration is :math:`[-\inf, \inf]`. the HermiteE polynomials are
1721 orthogonal, but not normalized, with respect to this weight function.
1722
1723 Parameters
1724 ----------
1725 x : array_like
1726 Values at which the weight function will be computed.
1727
1728 Returns
1729 -------
1730 w : ndarray
1731 The weight function at `x`.
1732
1733 Notes
1734 -----
1735
1736 .. versionadded::1.7.0
1737
1738 """
1739 w = np.exp(-.5*x**2)
1740 return w
1741
1742
1743 #
1744 # HermiteE series class
1745 #
1746
1747 class HermiteE(ABCPolyBase):
1748 """An HermiteE series class.
1749
1750 The HermiteE class provides the standard Python numerical methods
1751 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
1752 attributes and methods listed in the `ABCPolyBase` documentation.
1753
1754 Parameters
1755 ----------
1756 coef : array_like
1757 Laguerre coefficients in order of increasing degree, i.e,
1758 ``(1, 2, 3)`` gives ``1*He_0(x) + 2*He_1(X) + 3*He_2(x)``.
1759 domain : (2,) array_like, optional
1760 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
1761 to the interval ``[window[0], window[1]]`` by shifting and scaling.
1762 The default value is [-1, 1].
1763 window : (2,) array_like, optional
1764 Window, see `domain` for its use. The default value is [-1, 1].
1765
1766 .. versionadded:: 1.6.0
1767
1768 """
1769 # Virtual Functions
1770 _add = staticmethod(hermeadd)
1771 _sub = staticmethod(hermesub)
1772 _mul = staticmethod(hermemul)
1773 _div = staticmethod(hermediv)
1774 _pow = staticmethod(hermepow)
1775 _val = staticmethod(hermeval)
1776 _int = staticmethod(hermeint)
1777 _der = staticmethod(hermeder)
1778 _fit = staticmethod(hermefit)
1779 _line = staticmethod(hermeline)
1780 _roots = staticmethod(hermeroots)
1781 _fromroots = staticmethod(hermefromroots)
1782
1783 # Virtual properties
1784 nickname = 'herme'
1785 domain = np.array(hermedomain)
1786 window = np.array(hermedomain)