Skip to content

Commit 8f7fb72

Browse files
jonathfjonathf
committed
including genz-keister back again
1 parent 9050080 commit 8f7fb72

File tree

5 files changed

+321
-5
lines changed

5 files changed

+321
-5
lines changed

‎chaospy/distributions/approximation.py

+2-3
Original file line numberDiff line numberDiff line change
@@ -156,7 +156,7 @@ def approximate_inverse(
156156
def approximate_moment(
157157
distribution,
158158
k_loc,
159-
order=None,
159+
order=100000,
160160
rule="clenshaw_curtis",
161161
**kwargs
162162
):
@@ -171,8 +171,7 @@ def approximate_moment(
171171
k_loc (Sequence[int, ...]):
172172
The exponents of the moments of interest with ``shape == (dim,)``.
173173
order (int):
174-
The quadrature order used in approximation. If omitted, defaults to
175-
``1e5**(1./len(distribution))``.
174+
The quadrature order used in approximation.
176175
rule (str):
177176
Quadrature rule for integrating moments.
178177
kwargs:

‎chaospy/quadrature/__init__.py

+7
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@
1212
from .fejer_1 import fejer_1
1313
from .fejer_2 import fejer_2
1414
from .gaussian import gaussian
15+
from .genz_keister import (
16+
genz_keister_16, genz_keister_18, genz_keister_22, genz_keister_24)
1517
from .gegenbauer import gegenbauer
1618
from .grid import grid
1719
from .hermite import hermite
@@ -34,6 +36,10 @@
3436
"fejer_1": fejer_1,
3537
"fejer_2": fejer_2,
3638
"gaussian": gaussian,
39+
"genz_keister_16": genz_keister_16,
40+
"genz_keister_18": genz_keister_18,
41+
"genz_keister_22": genz_keister_22,
42+
"genz_keister_24": genz_keister_24,
3743
"grid": grid,
3844
"kronrod": kronrod,
3945
"legendre": legendre_proxy,
@@ -73,3 +79,4 @@ def wrapped(*args, **kwargs):
7379
quadrature_deprecation_warning("gauss_lobatto", lobatto)
7480
quadrature_deprecation_warning("gauss_patterson", patterson)
7581
quadrature_deprecation_warning("gauss_radau", radau)
82+
quadrature_deprecation_warning("genz_keister", genz_keister_24)

‎chaospy/quadrature/frontend.py

+10
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,10 @@
2020
"n": "newton_cotes", "newton_cotes": "newton_cotes",
2121
"d": "discrete", "discrete": "discrete",
2222
"i": "grid", "grid": "grid",
23+
"z16": "genz_keister_16", "genz_keister_16": "genz_keister_16",
24+
"z18": "genz_keister_18", "genz_keister_18": "genz_keister_18",
25+
"z22": "genz_keister_22", "genz_keister_22": "genz_keister_22",
26+
"z24": "genz_keister_24", "genz_keister_24": "genz_keister_24",
2327
}
2428
DEPRECATED_SHORT_NAMES = {
2529
"f": "f2",
@@ -29,6 +33,8 @@
2933
"gauss_patterson": "patterson",
3034
"gauss_radau": "radau",
3135
"gauss_legendre": "legendre",
36+
"z": "genz_keister_24",
37+
"genz_keister": "genz_keister_24",
3238
}
3339

3440

@@ -113,6 +119,10 @@ def generate_quadrature(
113119
:func:`chaospy.quadrature.newton_cotes`
114120
:func:`chaospy.quadrature.discrete`
115121
:func:`chaospy.quadrature.grid`
122+
:func:`chaospy.quadrature.genz_keister_16`
123+
:func:`chaospy.quadrature.genz_keister_18`
124+
:func:`chaospy.quadrature.genz_keister_22`
125+
:func:`chaospy.quadrature.genz_keister_24`
116126
117127
"""
118128
if not rule:

‎chaospy/quadrature/gaussian.py

+2-2
Original file line numberDiff line numberDiff line change
@@ -25,10 +25,10 @@ def gaussian(
2525
algorithms exists.
2626
2727
Args:
28-
dist (chaospy.Distribution):
29-
The distribution which density will be used as weight function.
3028
order (int):
3129
The order of the quadrature.
30+
dist (chaospy.Distribution):
31+
The distribution which density will be used as weight function.
3232
recurrence_algorithm (str):
3333
Name of the algorithm used to generate abscissas and weights.
3434
rule (str):

‎chaospy/quadrature/genz_keister.py

+300
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,300 @@
1+
"""
2+
Hermite Genz-Keister quadrature rules
3+
4+
Adapted from John Burkardt's implementation in Matlab
5+
"""
6+
import numpy
7+
import scipy
8+
9+
from .utils import combine_quadrature
10+
11+
GENZ_KEISTER_STORE = {
12+
1: ((0.0000000000000000e+00,), (1.7724538509055159,)),
13+
3: ((0.0000000000000000e+00, 1.2247448713915889),
14+
(1.1816359006036772, 0.29540897515091930)),
15+
7: ((0.0000000000000000, 0.52403354748695763,
16+
1.2247448713915889, 2.9592107790638380),
17+
(0.81310410832613500, 0.23286251787386100,
18+
0.24557928535031393, 0.0012330680655153448)),
19+
9: ((0.0000000000000000, 0.52403354748695763, 1.2247448713915889,
20+
2.0232301911005157, 2.9592107790638380),
21+
(0.45014700975378197, 0.47869428549114124, 0.16811892894767771,
22+
0.014173117873979098, 1.6708826306882348e-4)),
23+
17: ((0.0000000000000000, 0.52403354748695763, 0.87004089535290285,
24+
1.2247448713915889, 1.8357079751751868, 2.0232301911005157,
25+
2.9592107790638380, 3.6677742159463378, 4.4995993983103881),
26+
(0.47310733504965385, 0.45119803602358544, 0.025155825701712934,
27+
0.15718298376652240, 0.0034840719346803800, 0.012466519132805918,
28+
1.8723818949278350e-04, -1.4542843387069391e-06, 3.7463469943051758e-08)),
29+
19: ((0.0000000000000000, 0.52403354748695763, 0.87004089535290285,
30+
1.2247448713915889, 1.8357079751751868, 2.0232301911005157,
31+
2.2665132620567876, 2.9592107790638380,
32+
3.6677742159463378, 4.4995993983103881),
33+
(0.53788160700510168, 0.36924643368920851, 0.10838861955003017,
34+
0.11360729895748269, 0.032055243099445879, -0.011232438489069229,
35+
5.1133174390883855e-03, 1.0656589772852267e-04,
36+
1.0802767206624762e-06, 1.5295717705322357e-09)),
37+
31: ((0.0000000000000000, 0.17606414208200893, 0.52403354748695763,
38+
0.87004089535290285, 1.2247448713915889, 1.5794121348467671,
39+
1.8357079751751868, 2.0232301911005157, 2.2665132620567876,
40+
2.5705583765842968, 2.9592107790638380, 3.6677742159463378,
41+
4.4995993983103881, 5.0360899444730940,
42+
5.6432578578857449, 6.3759392709822356),
43+
(0.45888839636756751, 0.049855761893293160, 0.35393889029580544,
44+
0.11594930984853116, 0.10939325071860877, 0.0031210210352682834,
45+
0.029409427580350787, -0.0098566270434610019, 0.0048385208205502612,
46+
2.6665159778939428e-05, 1.0541662394746661e-04, 1.0889219692128120e-06,
47+
1.4055252024722478e-09, 9.0675288231679823e-12,
48+
-2.6304696458548942e-13, 2.2365645607044459e-15)),
49+
33: ((0.0000000000000000, 0.17606414208200893, 0.52403354748695763,
50+
0.87004089535290285, 1.2247448713915889, 1.5794121348467671,
51+
1.8357079751751868, 2.0232301911005157, 2.2665132620567876,
52+
2.5705583765842968, 2.9592107790638380, 3.6677742159463378,
53+
4.0292201405043713, 4.4995993983103881, 5.0360899444730940,
54+
5.6432578578857449, 6.3759392709822356),
55+
(2.4656644932829619e-01, 1.8411696047725790e-01, 3.1208656194697448e-01,
56+
1.3726521191567551e-01, 9.6913444944583621e-02, 1.3032872699027960e-02,
57+
2.0435058359107205e-02, -4.9118576123877555e-03, 3.7580026604304793e-03,
58+
1.4753204901862772e-04, 9.8710009197409173e-05, 1.2245220967158438e-06,
59+
-2.3903343382803510e-08, 2.7547825138935901e-09, -3.4281570530349562e-11,
60+
4.7219278666417693e-13, -1.7602932805372496e-15)),
61+
35: ((0.0000000000000000e+00, 1.7606414208200893e-01, 5.2403354748695763e-01,
62+
8.7004089535290285e-01, 1.2247448713915889e+00, 1.5794121348467671e+00,
63+
1.8357079751751868e+00, 2.0232301911005157e+00, 2.2665132620567876e+00,
64+
2.5705583765842968e+00, 2.9592107790638380e+00, 3.3491639537131945e+00,
65+
3.6677742159463378e+00, 4.0292201405043713e+00, 4.4995993983103881e+00,
66+
5.0360899444730940e+00, 5.6432578578857449e+00, 6.3759392709822356e+00),
67+
(9.1262675363737921e-04, 3.3988595585585218e-01, 2.6244871488784277e-01,
68+
1.6371221555735804e-01, 8.0245518147390893e-02, 2.7780508908535097e-02,
69+
5.5928828911469180e-03, 4.0967527720344047e-03, 1.4515580425155904e-03,
70+
4.8785399304443770e-04, 6.3328620805617891e-05, 4.8462799737020461e-06,
71+
4.3737818040926989e-07, 3.7920222392319532e-08, 8.1553721816916897e-10,
72+
5.4896836948499462e-12, 9.6599466278563243e-15, 1.8684014894510604e-18)),
73+
37: ((0.000000000000000, 0.214618180588171, 0.524033547486958,
74+
0.870040895352903, 1.224744871391589, 1.561553427651873,
75+
1.835707975175187, 2.023230191100516, 2.266513262056788,
76+
2.597288631188366, 2.959210779063838, 3.315584617593290,
77+
3.667774215946338, 4.057956316089741, 4.499599398310388,
78+
4.986551454150765, 5.521865209868350, 6.124527854622158,
79+
6.853200069757519),
80+
(0.968824552928425499e-01, 0.147655710402686249e+00,
81+
0.143099302896833389e+00, 0.937208280655245902e-01,
82+
0.442116442189845444e-01, 0.15513109874859354e-01,
83+
0.43334988122723492e-02, 0.176802225818295443e-02,
84+
0.985827582996483824e-03, 0.234940366465975222e-03,
85+
0.32265185983739747e-04, 0.330975870979203419e-05,
86+
0.295907520230744049e-06, 0.16595448809389819e-07,
87+
0.422525843963111041e-09, 0.45661763676186859e-11,
88+
0.182242751549129356e-13, 0.187781893143728947e-16,
89+
0.19030350940130498e-20)),
90+
41: ((0.0000000000000000, 0.195324784415805, 0.52403354748695763,
91+
0.87004089535290285, 1.2247448713915889, 1.585873011819188,
92+
1.8357079751751868, 2.0232301911005157, 2.043834754429505,
93+
2.2665132620567876, 2.630415236459871, 2.9592107790638380,
94+
3.296114596212218, 3.6677742159463378, 4.070919267883068,
95+
4.4995993983103881, 4.953574342912980, 5.437443360177798,
96+
5.961461043404500, 6.547083258397540, 7.251792998192644),
97+
(0.562793426043218877e-01, 0.165639740400529554e+00,
98+
0.145966293895926429e+00, 0.928338228510111845e-01,
99+
0.45109010335859128e-01, 0.165445526705860772e-01,
100+
0.705471110122962612e-03, 0.178852543033699732e-01,
101+
- 0.144528422206988237e-01, 0.140697424065246825e-02,
102+
0.189010909805097887e-03, 0.288976780274478689e-04,
103+
0.381182791749177506e-05, 0.315372265852264871e-06,
104+
0.149158210417831408e-07, 0.400784141604834759e-09,
105+
0.581803393170320419e-11, 0.408820161202505983e-13,
106+
0.1140700785308509e-15, 0.860427172512207236e-19,
107+
0.664195893812757801e-23)),
108+
43: ((0.0000000000000000, 0.196029453662011, 0.52403354748695763,
109+
0.87004089535290285, 1.2247448713915889, 1.583643465293944,
110+
1.8357079751751868, 2.0232301911005157, 2.089340389294661,
111+
2.2665132620567876, 2.633356763661946, 2.9592107790638380,
112+
3.295265921534226, 3.6677742159463378, 4.071335874253583,
113+
4.4995993983103881, 4.952329763008589, 5.434053000365068,
114+
5.954781975039809, 6.535398426382995, 7.231746029072501,
115+
10.167574994881873),
116+
(0.579595986101181095e-01, 0.164880913687436689e+00,
117+
0.145863292632147353e+00, 0.928711584442575456e-01,
118+
0.450612329041864976e-01, 0.163616873493832402e-01,
119+
0.139966252291568061e-02, 0.67354758901013295e-02,
120+
-0.38799558623877157e-02, 0.150909333211638847e-02,
121+
0.184789465688357423e-03, 0.286802318064777813e-04,
122+
0.383880761947398577e-05, 0.316018363221289247e-06,
123+
0.148653643571796457e-07, 0.400030575425776948e-09,
124+
0.586915885251734856e-11, 0.421921851448196032e-13,
125+
0.122619614947864357e-15, 0.992619971560149097e-19,
126+
0.87544909871323873e-23, 0.546191947478318097e-37)),
127+
}
128+
129+
RULES = {
130+
16: [1, 3, 7, 9, 17, 19, 31],
131+
18: [1, 3, 9, 19, 37],
132+
22: [1, 3, 9, 19, 41],
133+
24: [1, 3, 9, 19, 43],
134+
}
135+
136+
def genz_keister_16(order, dist=None):
137+
"""
138+
Create Genz-Keister variant 16 quadrature nodes and weights.
139+
140+
Args:
141+
order (int, Sequence[int]):
142+
The order of the quadrature.
143+
dist (Optional[chaospy.Distribution]):
144+
The distribution which density will be used as weight function.
145+
If omitted, standard Gaussian is assumed.
146+
147+
Returns:
148+
(numpy.ndarray, numpy.ndarray):
149+
Genz-Keister quadrature abscissas and weights.
150+
151+
Examples:
152+
>>> nodes, weights = genz_keister_16(4)
153+
>>> nodes.round(2)
154+
array([[-6.36, -5.19, -4.18, -2.86, -2.6 , -1.73, -1.23, -0.74, 0. ,
155+
0.74, 1.23, 1.73, 2.6 , 2.86, 4.18, 5.19, 6.36]])
156+
>>> weights.round(8)
157+
array([ 2.0000000e-08, -8.2000000e-07, 1.0564000e-04, 7.0334800e-03,
158+
1.9656800e-03, 8.8681000e-02, 1.4192650e-02, 2.5456123e-01,
159+
2.6692223e-01, 2.5456123e-01, 1.4192650e-02, 8.8681000e-02,
160+
1.9656800e-03, 7.0334800e-03, 1.0564000e-04, -8.2000000e-07,
161+
2.0000000e-08])
162+
163+
"""
164+
return genz_keister(order, dist, rule=16)
165+
166+
167+
def genz_keister_18(order, dist=None):
168+
"""
169+
Create Genz-Keister variant 18 quadrature nodes and weights.
170+
171+
Args:
172+
order (int, Sequence[int]):
173+
The order of the quadrature.
174+
dist (Optional[chaospy.Distribution]):
175+
The distribution which density will be used as weight function.
176+
If omitted, standard Gaussian is assumed.
177+
178+
Returns:
179+
(numpy.ndarray, numpy.ndarray):
180+
Genz-Keister quadrature abscissas and weights.
181+
182+
Examples:
183+
>>> nodes, weights = genz_keister_18(2)
184+
>>> nodes.round(2)
185+
array([[-4.18, -2.86, -1.73, -0.74, 0. , 0.74, 1.73, 2.86, 4.18]])
186+
>>> weights.round(8)
187+
array([9.4270000e-05, 7.9963300e-03, 9.4850950e-02, 2.7007433e-01,
188+
2.5396825e-01, 2.7007433e-01, 9.4850950e-02, 7.9963300e-03,
189+
9.4270000e-05])
190+
191+
"""
192+
return genz_keister(order, dist, rule=18)
193+
194+
195+
def genz_keister_22(order, dist=None):
196+
"""
197+
Create Genz-Keister variant 22 quadrature nodes and weights.
198+
199+
Args:
200+
order (int, Sequence[int]):
201+
The order of the quadrature.
202+
dist (Optional[chaospy.Distribution]):
203+
The distribution which density will be used as weight function.
204+
If omitted, standard Gaussian is assumed.
205+
206+
Returns:
207+
(numpy.ndarray, numpy.ndarray):
208+
Genz-Keister quadrature abscissas and weights.
209+
210+
Examples:
211+
>>> nodes, weights = genz_keister_22(2)
212+
>>> nodes.round(2)
213+
array([[-4.18, -2.86, -1.73, -0.74, 0. , 0.74, 1.73, 2.86, 4.18]])
214+
>>> weights.round(8)
215+
array([9.4270000e-05, 7.9963300e-03, 9.4850950e-02, 2.7007433e-01,
216+
2.5396825e-01, 2.7007433e-01, 9.4850950e-02, 7.9963300e-03,
217+
9.4270000e-05])
218+
219+
"""
220+
return genz_keister(order, dist, rule=22)
221+
222+
223+
def genz_keister_24(order, dist=None):
224+
"""
225+
Create Genz-Keister variant 24 quadrature nodes and weights.
226+
227+
Args:
228+
order (int, Sequence[int]):
229+
The order of the quadrature.
230+
dist (Optional[chaospy.Distribution]):
231+
The distribution which density will be used as weight function.
232+
If omitted, standard Gaussian is assumed.
233+
234+
Returns:
235+
(numpy.ndarray, numpy.ndarray):
236+
Genz-Keister quadrature abscissas and weights.
237+
238+
Examples:
239+
>>> nodes, weights = genz_keister_24(2)
240+
>>> nodes.round(2)
241+
array([[-4.18, -2.86, -1.73, -0.74, 0. , 0.74, 1.73, 2.86, 4.18]])
242+
>>> weights.round(8)
243+
array([9.4270000e-05, 7.9963300e-03, 9.4850950e-02, 2.7007433e-01,
244+
2.5396825e-01, 2.7007433e-01, 9.4850950e-02, 7.9963300e-03,
245+
9.4270000e-05])
246+
247+
"""
248+
return genz_keister(order, dist, rule=24)
249+
250+
251+
def genz_keister(order, dist=None, rule=24):
252+
"""
253+
Create Genz-Keister quadrature nodes and weights.
254+
255+
Args:
256+
order (int, Sequence[int]):
257+
The order of the quadrature.
258+
dist (Optional[chaospy.Distribution]):
259+
The distribution which density will be used as weight function.
260+
If omitted, standard Gaussian is assumed.
261+
rule (int, Sequence[int]):
262+
The Genz-Keister rule name. Supported rules are 16, 18, 22 and 24.
263+
264+
Returns:
265+
(numpy.ndarray, numpy.ndarray):
266+
Genz-Keister quadrature abscissas and weights.
267+
268+
Examples:
269+
>>> genz_keister(0)
270+
(array([[0.]]), array([1.]))
271+
>>> genz_keister(1) # doctest: +NORMALIZE_WHITESPACE
272+
(array([[-1.73205081, 0. , 1.73205081]]),
273+
array([0.16666667, 0.66666667, 0.16666667]))
274+
275+
"""
276+
shape = (1,) if dist is None else (len(dist),)
277+
order = numpy.broadcast_to(order, shape)
278+
rule = numpy.broadcast_to(rule, shape)
279+
nodes, weights = zip(*[_genz_keister(order_, rule_)
280+
for order_, rule_ in zip(order, rule)])
281+
nodes, weights = combine_quadrature(nodes, weights)
282+
if dist is not None:
283+
nodes = dist.inv(scipy.special.ndtr(nodes))
284+
return nodes, weights
285+
286+
287+
def _genz_keister(order, rule):
288+
assert rule in RULES, "rule %d not in known rules: %s" % (rule, list(RULES))
289+
assert order <= len(RULES[rule]), (
290+
"rule genz_keister_%d limited at order %d" % (rule, order))
291+
order = RULES[rule][order]
292+
nodes, weights = GENZ_KEISTER_STORE[order]
293+
length = len(nodes)
294+
nodes = numpy.array(nodes[::-1]+nodes[1:])
295+
nodes[:length-1] *= -1
296+
nodes *= numpy.sqrt(2)
297+
weights = numpy.array(weights[::-1]+weights[1:])
298+
weights /= numpy.sum(weights)
299+
300+
return nodes, weights

0 commit comments

Comments
 (0)