Mercurial > hg > vamp-build-and-test
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) |