Skip to content

cubature

Cubature rules.

PositiveCubatureRule ¤

Bases: NamedTuple

Cubature rule with positive weights.

Source code in probdiffeq/solvers/strategies/components/cubature.py
 9
10
11
12
13
class PositiveCubatureRule(containers.NamedTuple):
    """Cubature rule with positive weights."""

    points: Array
    weights_sqrtm: Array

gauss_hermite(input_shape, degree=5) -> PositiveCubatureRule ¤

(Statistician's) Gauss-Hermite cubature.

The number of cubature points is prod(input_shape)**degree.

Source code in probdiffeq/solvers/strategies/components/cubature.py
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def gauss_hermite(input_shape, degree=5) -> PositiveCubatureRule:
    """(Statistician's) Gauss-Hermite cubature.

    The number of cubature points is `prod(input_shape)**degree`.
    """
    assert len(input_shape) == 1
    (dim,) = input_shape

    # Roots of the probabilist/statistician's Hermite polynomials (in Numpy...)
    _roots = special.roots_hermitenorm(n=degree, mu=True)
    pts, weights, sum_of_weights = _roots
    weights = weights / sum_of_weights

    # Transform into jax arrays and take square root of weights
    pts = np.asarray(pts)
    weights_sqrtm = np.sqrt(np.asarray(weights))

    # Build a tensor grid and return class
    tensor_pts = _tensor_points(pts, d=dim)
    tensor_weights_sqrtm = _tensor_weights(weights_sqrtm, d=dim)
    return PositiveCubatureRule(points=tensor_pts, weights_sqrtm=tensor_weights_sqrtm)

third_order_spherical(input_shape) -> PositiveCubatureRule ¤

Third-order spherical cubature integration.

Source code in probdiffeq/solvers/strategies/components/cubature.py
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def third_order_spherical(input_shape) -> PositiveCubatureRule:
    """Third-order spherical cubature integration."""
    assert len(input_shape) <= 1
    if len(input_shape) == 1:
        (d,) = input_shape
        points_mat, weights_sqrtm = _third_order_spherical_params(d=d)
        return PositiveCubatureRule(points=points_mat, weights_sqrtm=weights_sqrtm)

    # If input_shape == (), compute weights via input_shape=(1,)
    # and 'squeeze' the points.
    points_mat, weights_sqrtm = _third_order_spherical_params(d=1)
    (S, _) = points_mat.shape
    points = np.reshape(points_mat, (S,))
    return PositiveCubatureRule(points=points, weights_sqrtm=weights_sqrtm)

unscented_transform(input_shape, r=1.0) -> PositiveCubatureRule ¤

Unscented transform.

Source code in probdiffeq/solvers/strategies/components/cubature.py
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def unscented_transform(input_shape, r=1.0) -> PositiveCubatureRule:
    """Unscented transform."""
    assert len(input_shape) <= 1
    if len(input_shape) == 1:
        (d,) = input_shape
        points_mat, weights_sqrtm = _unscented_transform_params(d=d, r=r)
        return PositiveCubatureRule(points=points_mat, weights_sqrtm=weights_sqrtm)

    # If input_shape == (), compute weights via input_shape=(1,)
    # and 'squeeze' the points.
    points_mat, weights_sqrtm = _unscented_transform_params(d=1, r=r)
    (S, _) = points_mat.shape
    points = np.reshape(points_mat, (S,))
    return PositiveCubatureRule(points=points, weights_sqrtm=weights_sqrtm)