import math
import itertools
from functools import cached_property
from scipy._lib._array_api import array_namespace, np_compat
from scipy.integrate._rules import NestedFixedRule
class GenzMalikCubature(NestedFixedRule):
"""
Genz-Malik cubature.
Genz-Malik is only defined for integrals of dimension >= 2.
Parameters
----------
ndim : int
The spatial dimension of the integrand.
xp : array_namespace, optional
The namespace for the node and weight arrays. Default is None, where NumPy is
used.
Attributes
----------
higher : Cubature
Higher-order rule.
lower : Cubature
Lower-order rule.
References
----------
.. [1] A.C. Genz, A.A. Malik, Remarks on algorithm 006: An adaptive algorithm for
numerical integration over an N-dimensional rectangular region, Journal of
Computational and Applied Mathematics, Volume 6, Issue 4, 1980, Pages 295-302,
ISSN 0377-0427, https://doi.org/10.1016/0771-050X(80)90039-X.
Examples
--------
Evaluate a 3D integral:
>>> import numpy as np
>>> from scipy.integrate import cubature
>>> from scipy.integrate._rules import GenzMalikCubature
>>> def f(x):
... # f(x) = cos(x_1) + cos(x_2) + cos(x_3)
... return np.sum(np.cos(x), axis=-1)
>>> rule = GenzMalikCubature(3) # Use 3D Genz-Malik
>>> a, b = np.array([0, 0, 0]), np.array([1, 1, 1])
>>> rule.estimate(f, a, b) # True value 3*sin(1), approximately 2.5244
np.float64(2.5244129547230862)
>>> rule.estimate_error(f, a, b)
np.float64(1.378269656626685e-06)
"""
def __init__(self, ndim, degree=7, lower_degree=5, xp=None):
if ndim < 2:
raise ValueError("Genz-Malik cubature is only defined for ndim >= 2")
if degree != 7 or lower_degree != 5:
raise NotImplementedError("Genz-Malik cubature is currently only supported"
"for degree=7, lower_degree=5")
self.ndim = ndim
self.degree = degree
self.lower_degree = lower_degree
if xp is None:
xp = np_compat
self.xp = array_namespace(xp.empty(0))
@cached_property
def nodes_and_weights(self):
# TODO: Currently only support for degree 7 Genz-Malik cubature, should aim to
# support arbitrary degree
l_2 = math.sqrt(9/70)
l_3 = math.sqrt(9/10)
l_4 = math.sqrt(9/10)
l_5 = math.sqrt(9/19)
its = itertools.chain(
[(0,) * self.ndim],
_distinct_permutations((l_2,) + (0,) * (self.ndim - 1)),
_distinct_permutations((-l_2,) + (0,) * (self.ndim - 1)),
_distinct_permutations((l_3,) + (0,) * (self.ndim - 1)),
_distinct_permutations((-l_3,) + (0,) * (self.ndim - 1)),
_distinct_permutations((l_4, l_4) + (0,) * (self.ndim - 2)),
_distinct_permutations((l_4, -l_4) + (0,) * (self.ndim - 2)),
_distinct_permutations((-l_4, -l_4) + (0,) * (self.ndim - 2)),
itertools.product((l_5, -l_5), repeat=self.ndim),
)
nodes_size = 1 + (2 * (self.ndim + 1) * self.ndim) + 2**self.ndim
nodes = self.xp.asarray(
list(zip(*its)),
dtype=self.xp.float64,
)
nodes = self.xp.reshape(nodes, (self.ndim, nodes_size))
# It's convenient to generate the nodes as a sequence of evaluation points
# as an array of shape (npoints, ndim), but nodes needs to have shape
# (ndim, npoints)
nodes = nodes.T
w_1 = (
(2**self.ndim) * (12824 - 9120*self.ndim + (400 * self.ndim**2)) / 19683
)
w_2 = (2**self.ndim) * 980/6561
w_3 = (2**self.ndim) * (1820 - 400 * self.ndim) / 19683
w_4 = (2**self.ndim) * (200 / 19683)
w_5 = 6859 / 19683
weights = self.xp.concat([
self.xp.asarray([w_1] * 1, dtype=self.xp.float64),
self.xp.asarray([w_2] * (2 * self.ndim), dtype=self.xp.float64),
self.xp.asarray([w_3] * (2 * self.ndim), dtype=self.xp.float64),
self.xp.asarray(
[w_4] * (2 * (self.ndim - 1) * self.ndim),
dtype=self.xp.float64,
),
self.xp.asarray([w_5] * (2**self.ndim), dtype=self.xp.float64),
])
return nodes, weights
@cached_property
def lower_nodes_and_weights(self):
# TODO: Currently only support for the degree 5 lower rule, in the future it
# would be worth supporting arbitrary degree
# Nodes are almost the same as the full rule, but there are no nodes
# corresponding to l_5.
l_2 = math.sqrt(9/70)
l_3 = math.sqrt(9/10)
l_4 = math.sqrt(9/10)
its = itertools.chain(
[(0,) * self.ndim],
_distinct_permutations((l_2,) + (0,) * (self.ndim - 1)),
_distinct_permutations((-l_2,) + (0,) * (self.ndim - 1)),
_distinct_permutations((l_3,) + (0,) * (self.ndim - 1)),
_distinct_permutations((-l_3,) + (0,) * (self.ndim - 1)),
_distinct_permutations((l_4, l_4) + (0,) * (self.ndim - 2)),
_distinct_permutations((l_4, -l_4) + (0,) * (self.ndim - 2)),
_distinct_permutations((-l_4, -l_4) + (0,) * (self.ndim - 2)),
)
nodes_size = 1 + (2 * (self.ndim + 1) * self.ndim)
nodes = self.xp.asarray(list(zip(*its)), dtype=self.xp.float64)
nodes = self.xp.reshape(nodes, (self.ndim, nodes_size))
nodes = nodes.T
# Weights are different from those in the full rule.
w_1 = (2**self.ndim) * (729 - 950*self.ndim + 50*self.ndim**2) / 729
w_2 = (2**self.ndim) * (245 / 486)
w_3 = (2**self.ndim) * (265 - 100*self.ndim) / 1458
w_4 = (2**self.ndim) * (25 / 729)
weights = self.xp.concat([
self.xp.asarray([w_1] * 1, dtype=self.xp.float64),
self.xp.asarray([w_2] * (2 * self.ndim), dtype=self.xp.float64),
self.xp.asarray([w_3] * (2 * self.ndim), dtype=self.xp.float64),
self.xp.asarray(
[w_4] * (2 * (self.ndim - 1) * self.ndim),
dtype=self.xp.float64,
),
])
return nodes, weights
def _distinct_permutations(iterable):
"""
Find the number of distinct permutations of elements of `iterable`.
"""
# Algorithm: https://w.wiki/Qai
items = sorted(iterable)
size = len(items)
while True:
# Yield the permutation we have
yield tuple(items)
# Find the largest index i such that A[i] < A[i + 1]
for i in range(size - 2, -1, -1):
if items[i] < items[i + 1]:
break
# If no such index exists, this permutation is the last one
else:
return
# Find the largest index j greater than j such that A[i] < A[j]
for j in range(size - 1, i, -1):
if items[i] < items[j]:
break
# Swap the value of A[i] with that of A[j], then reverse the
# sequence from A[i + 1] to form the new permutation
items[i], items[j] = items[j], items[i]
items[i+1:] = items[:i-size:-1] # A[i + 1:][::-1]