s474139/Traktor
2
1
Fork 0
Code Issues Pull Requests 1 Packages Projects Releases Wiki Activity
14640a2b12
Traktor/myenv/Lib/site-packages/scipy/integrate/_quadrature.py

1831 lines
64 KiB
Python
Raw Normal View History

losowanie zdjec
2024-05-26 05:12:46 +02:00
from __future__ import annotations
from typing import TYPE_CHECKING, Callable, Any, cast
import numpy as np
import numpy.typing as npt
import math
import warnings
from collections import namedtuple
from scipy.special import roots_legendre
from scipy.special import gammaln, logsumexp
from scipy._lib._util import _rng_spawn
from scipy._lib.deprecation import (_NoValue, _deprecate_positional_args,
_deprecated)
__all__ = ['fixed_quad', 'quadrature', 'romberg', 'romb',
'trapezoid', 'trapz', 'simps', 'simpson',
'cumulative_trapezoid', 'cumtrapz', 'newton_cotes',
'qmc_quad', 'AccuracyWarning', 'cumulative_simpson']
def trapezoid(y, x=None, dx=1.0, axis=-1):
r"""
Integrate along the given axis using the composite trapezoidal rule.
If `x` is provided, the integration happens in sequence along its
elements - they are not sorted.
Integrate `y` (`x`) along each 1d slice on the given axis, compute
:math:`\int y(x) dx`.
When `x` is specified, this integrates along the parametric curve,
computing :math:`\int_t y(t) dt =
\int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
Parameters
----------
y : array_like
Input array to integrate.
x : array_like, optional
The sample points corresponding to the `y` values. If `x` is None,
the sample points are assumed to be evenly spaced `dx` apart. The
default is None.
dx : scalar, optional
The spacing between sample points when `x` is None. The default is 1.
axis : int, optional
The axis along which to integrate.
Returns
-------
trapezoid : float or ndarray
Definite integral of `y` = n-dimensional array as approximated along
a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
then the result is a float. If `n` is greater than 1, then the result
is an `n`-1 dimensional array.
See Also
--------
cumulative_trapezoid, simpson, romb
Notes
-----
Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
will be taken from `y` array, by default x-axis distances between
points will be 1.0, alternatively they can be provided with `x` array
or with `dx` scalar. Return value will be equal to combined area under
the red lines.
References
----------
.. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
.. [2] Illustration image:
https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
Examples
--------
Use the trapezoidal rule on evenly spaced points:
>>> import numpy as np
>>> from scipy import integrate
>>> integrate.trapezoid([1, 2, 3])
4.0
The spacing between sample points can be selected by either the
``x`` or ``dx`` arguments:
>>> integrate.trapezoid([1, 2, 3], x=[4, 6, 8])
8.0
>>> integrate.trapezoid([1, 2, 3], dx=2)
8.0
Using a decreasing ``x`` corresponds to integrating in reverse:
>>> integrate.trapezoid([1, 2, 3], x=[8, 6, 4])
-8.0
More generally ``x`` is used to integrate along a parametric curve. We can
estimate the integral :math:`\int_0^1 x^2 = 1/3` using:
>>> x = np.linspace(0, 1, num=50)
>>> y = x**2
>>> integrate.trapezoid(y, x)
0.33340274885464394
Or estimate the area of a circle, noting we repeat the sample which closes
the curve:
>>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
>>> integrate.trapezoid(np.cos(theta), x=np.sin(theta))
3.141571941375841
``trapezoid`` can be applied along a specified axis to do multiple
computations in one call:
>>> a = np.arange(6).reshape(2, 3)
>>> a
array([[0, 1, 2],
[3, 4, 5]])
>>> integrate.trapezoid(a, axis=0)
array([1.5, 2.5, 3.5])
>>> integrate.trapezoid(a, axis=1)
array([2., 8.])
"""
y = np.asanyarray(y)
if x is None:
d = dx
else:
x = np.asanyarray(x)
if x.ndim == 1:
d = np.diff(x)
# reshape to correct shape
shape = [1]*y.ndim
shape[axis] = d.shape[0]
d = d.reshape(shape)
else:
d = np.diff(x, axis=axis)
nd = y.ndim
slice1 = [slice(None)]*nd
slice2 = [slice(None)]*nd
slice1[axis] = slice(1, None)
slice2[axis] = slice(None, -1)
try:
ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
except ValueError:
# Operations didn't work, cast to ndarray
d = np.asarray(d)
y = np.asarray(y)
ret = np.add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
return ret
# Note: alias kept for backwards compatibility. Rename was done
# because trapz is a slur in colloquial English (see gh-12924).
def trapz(y, x=None, dx=1.0, axis=-1):
"""An alias of `trapezoid`.
`trapz` is kept for backwards compatibility. For new code, prefer
`trapezoid` instead.
"""
msg = ("'scipy.integrate.trapz' is deprecated in favour of "
"'scipy.integrate.trapezoid' and will be removed in SciPy 1.14.0")
warnings.warn(msg, DeprecationWarning, stacklevel=2)
return trapezoid(y, x=x, dx=dx, axis=axis)
class AccuracyWarning(Warning):
pass
if TYPE_CHECKING:
# workaround for mypy function attributes see:
# https://github.com/python/mypy/issues/2087#issuecomment-462726600
from typing import Protocol
class CacheAttributes(Protocol):
cache: dict[int, tuple[Any, Any]]
else:
CacheAttributes = Callable
def cache_decorator(func: Callable) -> CacheAttributes:
return cast(CacheAttributes, func)
@cache_decorator
def _cached_roots_legendre(n):
"""
Cache roots_legendre results to speed up calls of the fixed_quad
function.
"""
if n in _cached_roots_legendre.cache:
return _cached_roots_legendre.cache[n]
_cached_roots_legendre.cache[n] = roots_legendre(n)
return _cached_roots_legendre.cache[n]
_cached_roots_legendre.cache = dict()
def fixed_quad(func, a, b, args=(), n=5):
"""
Compute a definite integral using fixed-order Gaussian quadrature.
Integrate `func` from `a` to `b` using Gaussian quadrature of
order `n`.
Parameters
----------
func : callable
A Python function or method to integrate (must accept vector inputs).
If integrating a vector-valued function, the returned array must have
shape ``(..., len(x))``.
a : float
Lower limit of integration.
b : float
Upper limit of integration.
args : tuple, optional
Extra arguments to pass to function, if any.
n : int, optional
Order of quadrature integration. Default is 5.
Returns
-------
val : float
Gaussian quadrature approximation to the integral
none : None
Statically returned value of None
See Also
--------
quad : adaptive quadrature using QUADPACK
dblquad : double integrals
tplquad : triple integrals
romberg : adaptive Romberg quadrature
quadrature : adaptive Gaussian quadrature
romb : integrators for sampled data
simpson : integrators for sampled data
cumulative_trapezoid : cumulative integration for sampled data
ode : ODE integrator
odeint : ODE integrator
Examples
--------
>>> from scipy import integrate
>>> import numpy as np
>>> f = lambda x: x**8
>>> integrate.fixed_quad(f, 0.0, 1.0, n=4)
(0.1110884353741496, None)
>>> integrate.fixed_quad(f, 0.0, 1.0, n=5)
(0.11111111111111102, None)
>>> print(1/9.0) # analytical result
0.1111111111111111
>>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=4)
(0.9999999771971152, None)
>>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=5)
(1.000000000039565, None)
>>> np.sin(np.pi/2)-np.sin(0) # analytical result
1.0
"""
x, w = _cached_roots_legendre(n)
x = np.real(x)
if np.isinf(a) or np.isinf(b):
raise ValueError("Gaussian quadrature is only available for "
"finite limits.")
y = (b-a)*(x+1)/2.0 + a
return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None
def vectorize1(func, args=(), vec_func=False):
"""Vectorize the call to a function.
This is an internal utility function used by `romberg` and
`quadrature` to create a vectorized version of a function.
If `vec_func` is True, the function `func` is assumed to take vector
arguments.
Parameters
----------
func : callable
User defined function.
args : tuple, optional
Extra arguments for the function.
vec_func : bool, optional
True if the function func takes vector arguments.
Returns
-------
vfunc : callable
A function that will take a vector argument and return the
result.
"""
if vec_func:
def vfunc(x):
return func(x, *args)
else:
def vfunc(x):
if np.isscalar(x):
return func(x, *args)
x = np.asarray(x)
# call with first point to get output type
y0 = func(x[0], *args)
n = len(x)
dtype = getattr(y0, 'dtype', type(y0))
output = np.empty((n,), dtype=dtype)
output[0] = y0
for i in range(1, n):
output[i] = func(x[i], *args)
return output
return vfunc
@_deprecated("`scipy.integrate.quadrature` is deprecated as of SciPy 1.12.0"
"and will be removed in SciPy 1.15.0. Please use"
"`scipy.integrate.quad` instead.")
def quadrature(func, a, b, args=(), tol=1.49e-8, rtol=1.49e-8, maxiter=50,
vec_func=True, miniter=1):
"""
Compute a definite integral using fixed-tolerance Gaussian quadrature.
.. deprecated:: 1.12.0
This function is deprecated as of SciPy 1.12.0 and will be removed
in SciPy 1.15.0. Please use `scipy.integrate.quad` instead.
Integrate `func` from `a` to `b` using Gaussian quadrature
with absolute tolerance `tol`.
Parameters
----------
func : function
A Python function or method to integrate.
a : float
Lower limit of integration.
b : float
Upper limit of integration.
args : tuple, optional
Extra arguments to pass to function.
tol, rtol : float, optional
Iteration stops when error between last two iterates is less than
`tol` OR the relative change is less than `rtol`.
maxiter : int, optional
Maximum order of Gaussian quadrature.
vec_func : bool, optional
True or False if func handles arrays as arguments (is
a "vector" function). Default is True.
miniter : int, optional
Minimum order of Gaussian quadrature.
Returns
-------
val : float
Gaussian quadrature approximation (within tolerance) to integral.
err : float
Difference between last two estimates of the integral.
See Also
--------
romberg : adaptive Romberg quadrature
fixed_quad : fixed-order Gaussian quadrature
quad : adaptive quadrature using QUADPACK
dblquad : double integrals
tplquad : triple integrals
romb : integrator for sampled data
simpson : integrator for sampled data
cumulative_trapezoid : cumulative integration for sampled data
ode : ODE integrator
odeint : ODE integrator
Examples
--------
>>> from scipy import integrate
>>> import numpy as np
>>> f = lambda x: x**8
>>> integrate.quadrature(f, 0.0, 1.0)
(0.11111111111111106, 4.163336342344337e-17)
>>> print(1/9.0) # analytical result
0.1111111111111111
>>> integrate.quadrature(np.cos, 0.0, np.pi/2)
(0.9999999999999536, 3.9611425250996035e-11)
>>> np.sin(np.pi/2)-np.sin(0) # analytical result
1.0
"""
if not isinstance(args, tuple):
args = (args,)
vfunc = vectorize1(func, args, vec_func=vec_func)
val = np.inf
err = np.inf
maxiter = max(miniter+1, maxiter)
for n in range(miniter, maxiter+1):
newval = fixed_quad(vfunc, a, b, (), n)[0]
err = abs(newval-val)
val = newval
if err < tol or err < rtol*abs(val):
break
else:
warnings.warn(
"maxiter (%d) exceeded. Latest difference = %e" % (maxiter, err),
AccuracyWarning, stacklevel=2
)
return val, err
def tupleset(t, i, value):
l = list(t)
l[i] = value
return tuple(l)
# Note: alias kept for backwards compatibility. Rename was done
# because cumtrapz is a slur in colloquial English (see gh-12924).
def cumtrapz(y, x=None, dx=1.0, axis=-1, initial=None):
"""An alias of `cumulative_trapezoid`.
`cumtrapz` is kept for backwards compatibility. For new code, prefer
`cumulative_trapezoid` instead.
"""
msg = ("'scipy.integrate.cumtrapz' is deprecated in favour of "
"'scipy.integrate.cumulative_trapezoid' and will be removed "
"in SciPy 1.14.0")
warnings.warn(msg, DeprecationWarning, stacklevel=2)
return cumulative_trapezoid(y, x=x, dx=dx, axis=axis, initial=initial)
def cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None):
"""
Cumulatively integrate y(x) using the composite trapezoidal rule.
Parameters
----------
y : array_like
Values to integrate.
x : array_like, optional
The coordinate to integrate along. If None (default), use spacing `dx`
between consecutive elements in `y`.
dx : float, optional
Spacing between elements of `y`. Only used if `x` is None.
axis : int, optional
Specifies the axis to cumulate. Default is -1 (last axis).
initial : scalar, optional
If given, insert this value at the beginning of the returned result.
0 or None are the only values accepted. Default is None, which means
`res` has one element less than `y` along the axis of integration.
.. deprecated:: 1.12.0
The option for non-zero inputs for `initial` will be deprecated in
SciPy 1.15.0. After this time, a ValueError will be raised if
`initial` is not None or 0.
Returns
-------
res : ndarray
The result of cumulative integration of `y` along `axis`.
If `initial` is None, the shape is such that the axis of integration
has one less value than `y`. If `initial` is given, the shape is equal
to that of `y`.
See Also
--------
numpy.cumsum, numpy.cumprod
cumulative_simpson : cumulative integration using Simpson's 1/3 rule
quad : adaptive quadrature using QUADPACK
romberg : adaptive Romberg quadrature
quadrature : adaptive Gaussian quadrature
fixed_quad : fixed-order Gaussian quadrature
dblquad : double integrals
tplquad : triple integrals
romb : integrators for sampled data
ode : ODE integrators
odeint : ODE integrators
Examples
--------
>>> from scipy import integrate
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-2, 2, num=20)
>>> y = x
>>> y_int = integrate.cumulative_trapezoid(y, x, initial=0)
>>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
>>> plt.show()
"""
y = np.asarray(y)
if y.shape[axis] == 0:
raise ValueError("At least one point is required along `axis`.")
if x is None:
d = dx
else:
x = np.asarray(x)
if x.ndim == 1:
d = np.diff(x)
# reshape to correct shape
shape = [1] * y.ndim
shape[axis] = -1
d = d.reshape(shape)
elif len(x.shape) != len(y.shape):
raise ValueError("If given, shape of x must be 1-D or the "
"same as y.")
else:
d = np.diff(x, axis=axis)
if d.shape[axis] != y.shape[axis] - 1:
raise ValueError("If given, length of x along axis must be the "
"same as y.")
nd = len(y.shape)
slice1 = tupleset((slice(None),)*nd, axis, slice(1, None))
slice2 = tupleset((slice(None),)*nd, axis, slice(None, -1))
res = np.cumsum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)
if initial is not None:
if initial != 0:
warnings.warn(
"The option for values for `initial` other than None or 0 is "
"deprecated as of SciPy 1.12.0 and will raise a value error in"
" SciPy 1.15.0.",
DeprecationWarning, stacklevel=2
)
if not np.isscalar(initial):
raise ValueError("`initial` parameter should be a scalar.")
shape = list(res.shape)
shape[axis] = 1
res = np.concatenate([np.full(shape, initial, dtype=res.dtype), res],
axis=axis)
return res
def _basic_simpson(y, start, stop, x, dx, axis):
nd = len(y.shape)
if start is None:
start = 0
step = 2
slice_all = (slice(None),)*nd
slice0 = tupleset(slice_all, axis, slice(start, stop, step))
slice1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
slice2 = tupleset(slice_all, axis, slice(start+2, stop+2, step))
if x is None: # Even-spaced Simpson's rule.
result = np.sum(y[slice0] + 4.0*y[slice1] + y[slice2], axis=axis)
result *= dx / 3.0
else:
# Account for possibly different spacings.
# Simpson's rule changes a bit.
h = np.diff(x, axis=axis)
sl0 = tupleset(slice_all, axis, slice(start, stop, step))
sl1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
h0 = h[sl0].astype(float, copy=False)
h1 = h[sl1].astype(float, copy=False)
hsum = h0 + h1
hprod = h0 * h1
h0divh1 = np.true_divide(h0, h1, out=np.zeros_like(h0), where=h1 != 0)
tmp = hsum/6.0 * (y[slice0] *
(2.0 - np.true_divide(1.0, h0divh1,
out=np.zeros_like(h0divh1),
where=h0divh1 != 0)) +
y[slice1] * (hsum *
np.true_divide(hsum, hprod,
out=np.zeros_like(hsum),
where=hprod != 0)) +
y[slice2] * (2.0 - h0divh1))
result = np.sum(tmp, axis=axis)
return result
# Note: alias kept for backwards compatibility. simps was renamed to simpson
# because the former is a slur in colloquial English (see gh-12924).
def simps(y, x=None, dx=1.0, axis=-1, even=_NoValue):
"""An alias of `simpson`.
`simps` is kept for backwards compatibility. For new code, prefer
`simpson` instead.
"""
msg = ("'scipy.integrate.simps' is deprecated in favour of "
"'scipy.integrate.simpson' and will be removed in SciPy 1.14.0")
warnings.warn(msg, DeprecationWarning, stacklevel=2)
# we don't deprecate positional use as the wrapper is going away completely
return simpson(y, x=x, dx=dx, axis=axis, even=even)
@_deprecate_positional_args(version="1.14")
def simpson(y, *, x=None, dx=1.0, axis=-1, even=_NoValue):
"""
Integrate y(x) using samples along the given axis and the composite
Simpson's rule. If x is None, spacing of dx is assumed.
If there are an even number of samples, N, then there are an odd
number of intervals (N-1), but Simpson's rule requires an even number
of intervals. The parameter 'even' controls how this is handled.
Parameters
----------
y : array_like
Array to be integrated.
x : array_like, optional
If given, the points at which `y` is sampled.
dx : float, optional
Spacing of integration points along axis of `x`. Only used when
`x` is None. Default is 1.
axis : int, optional
Axis along which to integrate. Default is the last axis.
even : {None, 'simpson', 'avg', 'first', 'last'}, optional
'avg' : Average two results:
1) use the first N-2 intervals with
a trapezoidal rule on the last interval and
2) use the last
N-2 intervals with a trapezoidal rule on the first interval.
'first' : Use Simpson's rule for the first N-2 intervals with
a trapezoidal rule on the last interval.
'last' : Use Simpson's rule for the last N-2 intervals with a
trapezoidal rule on the first interval.
None : equivalent to 'simpson' (default)
'simpson' : Use Simpson's rule for the first N-2 intervals with the
addition of a 3-point parabolic segment for the last
interval using equations outlined by Cartwright [1]_.
If the axis to be integrated over only has two points then
the integration falls back to a trapezoidal integration.
.. versionadded:: 1.11.0
.. versionchanged:: 1.11.0
The newly added 'simpson' option is now the default as it is more
accurate in most situations.
.. deprecated:: 1.11.0
Parameter `even` is deprecated and will be removed in SciPy
1.14.0. After this time the behaviour for an even number of
points will follow that of `even='simpson'`.
Returns
-------
float
The estimated integral computed with the composite Simpson's rule.
See Also
--------
quad : adaptive quadrature using QUADPACK
romberg : adaptive Romberg quadrature
quadrature : adaptive Gaussian quadrature
fixed_quad : fixed-order Gaussian quadrature
dblquad : double integrals
tplquad : triple integrals
romb : integrators for sampled data
cumulative_trapezoid : cumulative integration for sampled data
cumulative_simpson : cumulative integration using Simpson's 1/3 rule
ode : ODE integrators
odeint : ODE integrators
Notes
-----
For an odd number of samples that are equally spaced the result is
exact if the function is a polynomial of order 3 or less. If
the samples are not equally spaced, then the result is exact only
if the function is a polynomial of order 2 or less.
References
----------
.. [1] Cartwright, Kenneth V. Simpson's Rule Cumulative Integration with
MS Excel and Irregularly-spaced Data. Journal of Mathematical
Sciences and Mathematics Education. 12 (2): 1-9
Examples
--------
>>> from scipy import integrate
>>> import numpy as np
>>> x = np.arange(0, 10)
>>> y = np.arange(0, 10)
>>> integrate.simpson(y, x=x)
40.5
>>> y = np.power(x, 3)
>>> integrate.simpson(y, x=x)
1640.5
>>> integrate.quad(lambda x: x**3, 0, 9)[0]
1640.25
>>> integrate.simpson(y, x=x, even='first')
1644.5
"""
y = np.asarray(y)
nd = len(y.shape)
N = y.shape[axis]
last_dx = dx
first_dx = dx
returnshape = 0
if x is not None:
x = np.asarray(x)
if len(x.shape) == 1:
shapex = [1] * nd
shapex[axis] = x.shape[0]
saveshape = x.shape
returnshape = 1
x = x.reshape(tuple(shapex))
elif len(x.shape) != len(y.shape):
raise ValueError("If given, shape of x must be 1-D or the "
"same as y.")
if x.shape[axis] != N:
raise ValueError("If given, length of x along axis must be the "
"same as y.")
# even keyword parameter is deprecated
if even is not _NoValue:
warnings.warn(
"The 'even' keyword is deprecated as of SciPy 1.11.0 and will be "
"removed in SciPy 1.14.0",
DeprecationWarning, stacklevel=2
)
if N % 2 == 0:
val = 0.0
result = 0.0
slice_all = (slice(None),) * nd
# default is 'simpson'
even = even if even not in (_NoValue, None) else "simpson"
if even not in ['avg', 'last', 'first', 'simpson']:
raise ValueError(
"Parameter 'even' must be 'simpson', "
"'avg', 'last', or 'first'."
)
if N == 2:
# need at least 3 points in integration axis to form parabolic
# segment. If there are two points then any of 'avg', 'first',
# 'last' should give the same result.
slice1 = tupleset(slice_all, axis, -1)
slice2 = tupleset(slice_all, axis, -2)
if x is not None:
last_dx = x[slice1] - x[slice2]
val += 0.5 * last_dx * (y[slice1] + y[slice2])
# calculation is finished. Set `even` to None to skip other
# scenarios
even = None
if even == 'simpson':
# use Simpson's rule on first intervals
result = _basic_simpson(y, 0, N-3, x, dx, axis)
slice1 = tupleset(slice_all, axis, -1)
slice2 = tupleset(slice_all, axis, -2)
slice3 = tupleset(slice_all, axis, -3)
h = np.asarray([dx, dx], dtype=np.float64)
if x is not None:
# grab the last two spacings from the appropriate axis
hm2 = tupleset(slice_all, axis, slice(-2, -1, 1))
hm1 = tupleset(slice_all, axis, slice(-1, None, 1))
diffs = np.float64(np.diff(x, axis=axis))
h = [np.squeeze(diffs[hm2], axis=axis),
np.squeeze(diffs[hm1], axis=axis)]
# This is the correction for the last interval according to
# Cartwright.
# However, I used the equations given at
# https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
# A footnote on Wikipedia says:
# Cartwright 2017, Equation 8. The equation in Cartwright is
# calculating the first interval whereas the equations in the
# Wikipedia article are adjusting for the last integral. If the
# proper algebraic substitutions are made, the equation results in
# the values shown.
num = 2 * h[1] ** 2 + 3 * h[0] * h[1]
den = 6 * (h[1] + h[0])
alpha = np.true_divide(
num,
den,
out=np.zeros_like(den),
where=den != 0
)
num = h[1] ** 2 + 3.0 * h[0] * h[1]
den = 6 * h[0]
beta = np.true_divide(
num,
den,
out=np.zeros_like(den),
where=den != 0
)
num = 1 * h[1] ** 3
den = 6 * h[0] * (h[0] + h[1])
eta = np.true_divide(
num,
den,
out=np.zeros_like(den),
where=den != 0
)
result += alpha*y[slice1] + beta*y[slice2] - eta*y[slice3]
# The following code (down to result=result+val) can be removed
# once the 'even' keyword is removed.
# Compute using Simpson's rule on first intervals
if even in ['avg', 'first']:
slice1 = tupleset(slice_all, axis, -1)
slice2 = tupleset(slice_all, axis, -2)
if x is not None:
last_dx = x[slice1] - x[slice2]
val += 0.5*last_dx*(y[slice1]+y[slice2])
result = _basic_simpson(y, 0, N-3, x, dx, axis)
# Compute using Simpson's rule on last set of intervals
if even in ['avg', 'last']:
slice1 = tupleset(slice_all, axis, 0)
slice2 = tupleset(slice_all, axis, 1)
if x is not None:
first_dx = x[tuple(slice2)] - x[tuple(slice1)]
val += 0.5*first_dx*(y[slice2]+y[slice1])
result += _basic_simpson(y, 1, N-2, x, dx, axis)
if even == 'avg':
val /= 2.0
result /= 2.0
result = result + val
else:
result = _basic_simpson(y, 0, N-2, x, dx, axis)
if returnshape:
x = x.reshape(saveshape)
return result
def _cumulatively_sum_simpson_integrals(
y: np.ndarray,
dx: np.ndarray,
integration_func: Callable[[np.ndarray, np.ndarray], np.ndarray],
) -> np.ndarray:
"""Calculate cumulative sum of Simpson integrals.
Takes as input the integration function to be used.
The integration_func is assumed to return the cumulative sum using
composite Simpson's rule. Assumes the axis of summation is -1.
"""
sub_integrals_h1 = integration_func(y, dx)
sub_integrals_h2 = integration_func(y[..., ::-1], dx[..., ::-1])[..., ::-1]
shape = list(sub_integrals_h1.shape)
shape[-1] += 1
sub_integrals = np.empty(shape)
sub_integrals[..., :-1:2] = sub_integrals_h1[..., ::2]
sub_integrals[..., 1::2] = sub_integrals_h2[..., ::2]
# Integral over last subinterval can only be calculated from
# formula for h2
sub_integrals[..., -1] = sub_integrals_h2[..., -1]
res = np.cumsum(sub_integrals, axis=-1)
return res
def _cumulative_simpson_equal_intervals(y: np.ndarray, dx: np.ndarray) -> np.ndarray:
"""Calculate the Simpson integrals for all h1 intervals assuming equal interval
widths. The function can also be used to calculate the integral for all
h2 intervals by reversing the inputs, `y` and `dx`.
"""
d = dx[..., :-1]
f1 = y[..., :-2]
f2 = y[..., 1:-1]
f3 = y[..., 2:]
# Calculate integral over the subintervals (eqn (10) of Reference [2])
return d / 3 * (5 * f1 / 4 + 2 * f2 - f3 / 4)
def _cumulative_simpson_unequal_intervals(y: np.ndarray, dx: np.ndarray) -> np.ndarray:
"""Calculate the Simpson integrals for all h1 intervals assuming unequal interval
widths. The function can also be used to calculate the integral for all
h2 intervals by reversing the inputs, `y` and `dx`.
"""
x21 = dx[..., :-1]
x32 = dx[..., 1:]
f1 = y[..., :-2]
f2 = y[..., 1:-1]
f3 = y[..., 2:]
x31 = x21 + x32
x21_x31 = x21/x31
x21_x32 = x21/x32
x21x21_x31x32 = x21_x31 * x21_x32
# Calculate integral over the subintervals (eqn (8) of Reference [2])
coeff1 = 3 - x21_x31
coeff2 = 3 + x21x21_x31x32 + x21_x31
coeff3 = -x21x21_x31x32
return x21/6 * (coeff1*f1 + coeff2*f2 + coeff3*f3)
def _ensure_float_array(arr: npt.ArrayLike) -> np.ndarray:
arr = np.asarray(arr)
if np.issubdtype(arr.dtype, np.integer):
arr = arr.astype(float, copy=False)
return arr
def cumulative_simpson(y, *, x=None, dx=1.0, axis=-1, initial=None):
r"""
Cumulatively integrate y(x) using the composite Simpson's 1/3 rule.
The integral of the samples at every point is calculated by assuming a
quadratic relationship between each point and the two adjacent points.
Parameters
----------
y : array_like
Values to integrate. Requires at least one point along `axis`. If two or fewer
points are provided along `axis`, Simpson's integration is not possible and the
result is calculated with `cumulative_trapezoid`.
x : array_like, optional
The coordinate to integrate along. Must have the same shape as `y` or
must be 1D with the same length as `y` along `axis`. `x` must also be
strictly increasing along `axis`.
If `x` is None (default), integration is performed using spacing `dx`
between consecutive elements in `y`.
dx : scalar or array_like, optional
Spacing between elements of `y`. Only used if `x` is None. Can either
be a float, or an array with the same shape as `y`, but of length one along
`axis`. Default is 1.0.
axis : int, optional
Specifies the axis to integrate along. Default is -1 (last axis).
initial : scalar or array_like, optional
If given, insert this value at the beginning of the returned result,
and add it to the rest of the result. Default is None, which means no
value at ``x[0]`` is returned and `res` has one element less than `y`
along the axis of integration. Can either be a float, or an array with
the same shape as `y`, but of length one along `axis`.
Returns
-------
res : ndarray
The result of cumulative integration of `y` along `axis`.
If `initial` is None, the shape is such that the axis of integration
has one less value than `y`. If `initial` is given, the shape is equal
to that of `y`.
See Also
--------
numpy.cumsum
cumulative_trapezoid : cumulative integration using the composite
trapezoidal rule
simpson : integrator for sampled data using the Composite Simpson's Rule
Notes
-----
.. versionadded:: 1.12.0
The composite Simpson's 1/3 method can be used to approximate the definite
integral of a sampled input function :math:`y(x)` [1]_. The method assumes
a quadratic relationship over the interval containing any three consecutive
sampled points.
Consider three consecutive points:
:math:`(x_1, y_1), (x_2, y_2), (x_3, y_3)`.
Assuming a quadratic relationship over the three points, the integral over
the subinterval between :math:`x_1` and :math:`x_2` is given by formula
(8) of [2]_:
.. math::
\int_{x_1}^{x_2} y(x) dx\ &= \frac{x_2-x_1}{6}\left[\
\left\{3-\frac{x_2-x_1}{x_3-x_1}\right\} y_1 + \
\left\{3 + \frac{(x_2-x_1)^2}{(x_3-x_2)(x_3-x_1)} + \
\frac{x_2-x_1}{x_3-x_1}\right\} y_2\\
- \frac{(x_2-x_1)^2}{(x_3-x_2)(x_3-x_1)} y_3\right]
The integral between :math:`x_2` and :math:`x_3` is given by swapping
appearances of :math:`x_1` and :math:`x_3`. The integral is estimated
separately for each subinterval and then cumulatively summed to obtain
the final result.
For samples that are equally spaced, the result is exact if the function
is a polynomial of order three or less [1]_ and the number of subintervals
is even. Otherwise, the integral is exact for polynomials of order two or
less.
References
----------
.. [1] Wikipedia page: https://en.wikipedia.org/wiki/Simpson's_rule
.. [2] Cartwright, Kenneth V. Simpson's Rule Cumulative Integration with
MS Excel and Irregularly-spaced Data. Journal of Mathematical
Sciences and Mathematics Education. 12 (2): 1-9
Examples
--------
>>> from scipy import integrate
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-2, 2, num=20)
>>> y = x**2
>>> y_int = integrate.cumulative_simpson(y, x=x, initial=0)
>>> fig, ax = plt.subplots()
>>> ax.plot(x, y_int, 'ro', x, x**3/3 - (x[0])**3/3, 'b-')
>>> ax.grid()
>>> plt.show()
The output of `cumulative_simpson` is similar to that of iteratively
calling `simpson` with successively higher upper limits of integration, but
not identical.
>>> def cumulative_simpson_reference(y, x):
... return np.asarray([integrate.simpson(y[:i], x=x[:i])
... for i in range(2, len(y) + 1)])
>>>
>>> rng = np.random.default_rng(354673834679465)
>>> x, y = rng.random(size=(2, 10))
>>> x.sort()
>>>
>>> res = integrate.cumulative_simpson(y, x=x)
>>> ref = cumulative_simpson_reference(y, x)
>>> equal = np.abs(res - ref) < 1e-15
>>> equal # not equal when `simpson` has even number of subintervals
array([False, True, False, True, False, True, False, True, True])
This is expected: because `cumulative_simpson` has access to more
information than `simpson`, it can typically produce more accurate
estimates of the underlying integral over subintervals.
"""
y = _ensure_float_array(y)
# validate `axis` and standardize to work along the last axis
original_y = y
original_shape = y.shape
try:
y = np.swapaxes(y, axis, -1)
except IndexError as e:
message = f"`axis={axis}` is not valid for `y` with `y.ndim={y.ndim}`."
raise ValueError(message) from e
if y.shape[-1] < 3:
res = cumulative_trapezoid(original_y, x, dx=dx, axis=axis, initial=None)
res = np.swapaxes(res, axis, -1)
elif x is not None:
x = _ensure_float_array(x)
message = ("If given, shape of `x` must be the same as `y` or 1-D with "
"the same length as `y` along `axis`.")
if not (x.shape == original_shape
or (x.ndim == 1 and len(x) == original_shape[axis])):
raise ValueError(message)
x = np.broadcast_to(x, y.shape) if x.ndim == 1 else np.swapaxes(x, axis, -1)
dx = np.diff(x, axis=-1)
if np.any(dx <= 0):
raise ValueError("Input x must be strictly increasing.")
res = _cumulatively_sum_simpson_integrals(
y, dx, _cumulative_simpson_unequal_intervals
)
else:
dx = _ensure_float_array(dx)
final_dx_shape = tupleset(original_shape, axis, original_shape[axis] - 1)
alt_input_dx_shape = tupleset(original_shape, axis, 1)
message = ("If provided, `dx` must either be a scalar or have the same "
"shape as `y` but with only 1 point along `axis`.")
if not (dx.ndim == 0 or dx.shape == alt_input_dx_shape):
raise ValueError(message)
dx = np.broadcast_to(dx, final_dx_shape)
dx = np.swapaxes(dx, axis, -1)
res = _cumulatively_sum_simpson_integrals(
y, dx, _cumulative_simpson_equal_intervals
)
if initial is not None:
initial = _ensure_float_array(initial)
alt_initial_input_shape = tupleset(original_shape, axis, 1)
message = ("If provided, `initial` must either be a scalar or have the "
"same shape as `y` but with only 1 point along `axis`.")
if not (initial.ndim == 0 or initial.shape == alt_initial_input_shape):
raise ValueError(message)
initial = np.broadcast_to(initial, alt_initial_input_shape)
initial = np.swapaxes(initial, axis, -1)
res += initial
res = np.concatenate((initial, res), axis=-1)
res = np.swapaxes(res, -1, axis)
return res
def romb(y, dx=1.0, axis=-1, show=False):
"""
Romberg integration using samples of a function.
Parameters
----------
y : array_like
A vector of ``2**k + 1`` equally-spaced samples of a function.
dx : float, optional
The sample spacing. Default is 1.
axis : int, optional
The axis along which to integrate. Default is -1 (last axis).
show : bool, optional
When `y` is a single 1-D array, then if this argument is True
print the table showing Richardson extrapolation from the
samples. Default is False.
Returns
-------
romb : ndarray
The integrated result for `axis`.
See Also
--------
quad : adaptive quadrature using QUADPACK
romberg : adaptive Romberg quadrature
quadrature : adaptive Gaussian quadrature
fixed_quad : fixed-order Gaussian quadrature
dblquad : double integrals
tplquad : triple integrals
simpson : integrators for sampled data
cumulative_trapezoid : cumulative integration for sampled data
ode : ODE integrators
odeint : ODE integrators
Examples
--------
>>> from scipy import integrate
>>> import numpy as np
>>> x = np.arange(10, 14.25, 0.25)
>>> y = np.arange(3, 12)
>>> integrate.romb(y)
56.0
>>> y = np.sin(np.power(x, 2.5))
>>> integrate.romb(y)
-0.742561336672229
>>> integrate.romb(y, show=True)
Richardson Extrapolation Table for Romberg Integration
======================================================
-0.81576
4.63862 6.45674
-1.10581 -3.02062 -3.65245
-2.57379 -3.06311 -3.06595 -3.05664
-1.34093 -0.92997 -0.78776 -0.75160 -0.74256
======================================================
-0.742561336672229 # may vary
"""
y = np.asarray(y)
nd = len(y.shape)
Nsamps = y.shape[axis]
Ninterv = Nsamps-1
n = 1
k = 0
while n < Ninterv:
n <<= 1
k += 1
if n != Ninterv:
raise ValueError("Number of samples must be one plus a "
"non-negative power of 2.")
R = {}
slice_all = (slice(None),) * nd
slice0 = tupleset(slice_all, axis, 0)
slicem1 = tupleset(slice_all, axis, -1)
h = Ninterv * np.asarray(dx, dtype=float)
R[(0, 0)] = (y[slice0] + y[slicem1])/2.0*h
slice_R = slice_all
start = stop = step = Ninterv
for i in range(1, k+1):
start >>= 1
slice_R = tupleset(slice_R, axis, slice(start, stop, step))
step >>= 1
R[(i, 0)] = 0.5*(R[(i-1, 0)] + h*y[slice_R].sum(axis=axis))
for j in range(1, i+1):
prev = R[(i, j-1)]
R[(i, j)] = prev + (prev-R[(i-1, j-1)]) / ((1 << (2*j))-1)
h /= 2.0
if show:
if not np.isscalar(R[(0, 0)]):
print("*** Printing table only supported for integrals" +
" of a single data set.")
else:
try:
precis = show[0]
except (TypeError, IndexError):
precis = 5
try:
width = show[1]
except (TypeError, IndexError):
width = 8
formstr = "%%%d.%df" % (width, precis)
title = "Richardson Extrapolation Table for Romberg Integration"
print(title, "=" * len(title), sep="\n", end="\n")
for i in range(k+1):
for j in range(i+1):
print(formstr % R[(i, j)], end=" ")
print()
print("=" * len(title))
return R[(k, k)]
# Romberg quadratures for numeric integration.
#
# Written by Scott M. Ransom <ransom@cfa.harvard.edu>
# last revision: 14 Nov 98
#
# Cosmetic changes by Konrad Hinsen <hinsen@cnrs-orleans.fr>
# last revision: 1999-7-21
#
# Adapted to SciPy by Travis Oliphant <oliphant.travis@ieee.org>
# last revision: Dec 2001
def _difftrap(function, interval, numtraps):
"""
Perform part of the trapezoidal rule to integrate a function.
Assume that we had called difftrap with all lower powers-of-2
starting with 1. Calling difftrap only returns the summation
of the new ordinates. It does _not_ multiply by the width
of the trapezoids. This must be performed by the caller.
'function' is the function to evaluate (must accept vector arguments).
'interval' is a sequence with lower and upper limits
of integration.
'numtraps' is the number of trapezoids to use (must be a
power-of-2).
"""
if numtraps <= 0:
raise ValueError("numtraps must be > 0 in difftrap().")
elif numtraps == 1:
return 0.5*(function(interval[0])+function(interval[1]))
else:
numtosum = numtraps/2
h = float(interval[1]-interval[0])/numtosum
lox = interval[0] + 0.5 * h
points = lox + h * np.arange(numtosum)
s = np.sum(function(points), axis=0)
return s
def _romberg_diff(b, c, k):
"""
Compute the differences for the Romberg quadrature corrections.
See Forman Acton's "Real Computing Made Real," p 143.
"""
tmp = 4.0**k
return (tmp * c - b)/(tmp - 1.0)
def _printresmat(function, interval, resmat):
# Print the Romberg result matrix.
i = j = 0
print('Romberg integration of', repr(function), end=' ')
print('from', interval)
print('')
print('%6s %9s %9s' % ('Steps', 'StepSize', 'Results'))
for i in range(len(resmat)):
print('%6d %9f' % (2**i, (interval[1]-interval[0])/(2.**i)), end=' ')
for j in range(i+1):
print('%9f' % (resmat[i][j]), end=' ')
print('')
print('')
print('The final result is', resmat[i][j], end=' ')
print('after', 2**(len(resmat)-1)+1, 'function evaluations.')
@_deprecated("`scipy.integrate.romberg` is deprecated as of SciPy 1.12.0"
"and will be removed in SciPy 1.15.0. Please use"
"`scipy.integrate.quad` instead.")
def romberg(function, a, b, args=(), tol=1.48e-8, rtol=1.48e-8, show=False,
divmax=10, vec_func=False):
"""
Romberg integration of a callable function or method.
.. deprecated:: 1.12.0
This function is deprecated as of SciPy 1.12.0 and will be removed
in SciPy 1.15.0. Please use `scipy.integrate.quad` instead.
Returns the integral of `function` (a function of one variable)
over the interval (`a`, `b`).
If `show` is 1, the triangular array of the intermediate results
will be printed. If `vec_func` is True (default is False), then
`function` is assumed to support vector arguments.
Parameters
----------
function : callable
Function to be integrated.
a : float
Lower limit of integration.
b : float
Upper limit of integration.
Returns
-------
results : float
Result of the integration.
Other Parameters
----------------
args : tuple, optional
Extra arguments to pass to function. Each element of `args` will
be passed as a single argument to `func`. Default is to pass no
extra arguments.
tol, rtol : float, optional
The desired absolute and relative tolerances. Defaults are 1.48e-8.
show : bool, optional
Whether to print the results. Default is False.
divmax : int, optional
Maximum order of extrapolation. Default is 10.
vec_func : bool, optional
Whether `func` handles arrays as arguments (i.e., whether it is a
"vector" function). Default is False.
See Also
--------
fixed_quad : Fixed-order Gaussian quadrature.
quad : Adaptive quadrature using QUADPACK.
dblquad : Double integrals.
tplquad : Triple integrals.
romb : Integrators for sampled data.
simpson : Integrators for sampled data.
cumulative_trapezoid : Cumulative integration for sampled data.
ode : ODE integrator.
odeint : ODE integrator.
References
----------
.. [1] 'Romberg's method' https://en.wikipedia.org/wiki/Romberg%27s_method
Examples
--------
Integrate a gaussian from 0 to 1 and compare to the error function.
>>> from scipy import integrate
>>> from scipy.special import erf
>>> import numpy as np
>>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
>>> result = integrate.romberg(gaussian, 0, 1, show=True)
Romberg integration of <function vfunc at ...> from [0, 1]
::
Steps StepSize Results
1 1.000000 0.385872
2 0.500000 0.412631 0.421551
4 0.250000 0.419184 0.421368 0.421356
8 0.125000 0.420810 0.421352 0.421350 0.421350
16 0.062500 0.421215 0.421350 0.421350 0.421350 0.421350
32 0.031250 0.421317 0.421350 0.421350 0.421350 0.421350 0.421350
The final result is 0.421350396475 after 33 function evaluations.
>>> print("%g %g" % (2*result, erf(1)))
0.842701 0.842701
"""
if np.isinf(a) or np.isinf(b):
raise ValueError("Romberg integration only available "
"for finite limits.")
vfunc = vectorize1(function, args, vec_func=vec_func)
n = 1
interval = [a, b]
intrange = b - a
ordsum = _difftrap(vfunc, interval, n)
result = intrange * ordsum
resmat = [[result]]
err = np.inf
last_row = resmat[0]
for i in range(1, divmax+1):
n *= 2
ordsum += _difftrap(vfunc, interval, n)
row = [intrange * ordsum / n]
for k in range(i):
row.append(_romberg_diff(last_row[k], row[k], k+1))
result = row[i]
lastresult = last_row[i-1]
if show:
resmat.append(row)
err = abs(result - lastresult)
if err < tol or err < rtol * abs(result):
break
last_row = row
else:
warnings.warn(
"divmax (%d) exceeded. Latest difference = %e" % (divmax, err),
AccuracyWarning, stacklevel=2)
if show:
_printresmat(vfunc, interval, resmat)
return result
# Coefficients for Newton-Cotes quadrature
#
# These are the points being used
# to construct the local interpolating polynomial
# a are the weights for Newton-Cotes integration
# B is the error coefficient.
# error in these coefficients grows as N gets larger.
# or as samples are closer and closer together
# You can use maxima to find these rational coefficients
# for equally spaced data using the commands
# a(i,N) := (integrate(product(r-j,j,0,i-1) * product(r-j,j,i+1,N),r,0,N)
# / ((N-i)! * i!) * (-1)^(N-i));
# Be(N) := N^(N+2)/(N+2)! * (N/(N+3) - sum((i/N)^(N+2)*a(i,N),i,0,N));
# Bo(N) := N^(N+1)/(N+1)! * (N/(N+2) - sum((i/N)^(N+1)*a(i,N),i,0,N));
# B(N) := (if (mod(N,2)=0) then Be(N) else Bo(N));
#
# pre-computed for equally-spaced weights
#
# num_a, den_a, int_a, num_B, den_B = _builtincoeffs[N]
#
# a = num_a*array(int_a)/den_a
# B = num_B*1.0 / den_B
#
# integrate(f(x),x,x_0,x_N) = dx*sum(a*f(x_i)) + B*(dx)^(2k+3) f^(2k+2)(x*)
# where k = N // 2
#
_builtincoeffs = {
1: (1,2,[1,1],-1,12),
2: (1,3,[1,4,1],-1,90),
3: (3,8,[1,3,3,1],-3,80),
4: (2,45,[7,32,12,32,7],-8,945),
5: (5,288,[19,75,50,50,75,19],-275,12096),
6: (1,140,[41,216,27,272,27,216,41],-9,1400),
7: (7,17280,[751,3577,1323,2989,2989,1323,3577,751],-8183,518400),
8: (4,14175,[989,5888,-928,10496,-4540,10496,-928,5888,989],
-2368,467775),
9: (9,89600,[2857,15741,1080,19344,5778,5778,19344,1080,
15741,2857], -4671, 394240),
10: (5,299376,[16067,106300,-48525,272400,-260550,427368,
-260550,272400,-48525,106300,16067],
-673175, 163459296),
11: (11,87091200,[2171465,13486539,-3237113, 25226685,-9595542,
15493566,15493566,-9595542,25226685,-3237113,
13486539,2171465], -2224234463, 237758976000),
12: (1, 5255250, [1364651,9903168,-7587864,35725120,-51491295,
87516288,-87797136,87516288,-51491295,35725120,
-7587864,9903168,1364651], -3012, 875875),
13: (13, 402361344000,[8181904909, 56280729661, -31268252574,
156074417954,-151659573325,206683437987,
-43111992612,-43111992612,206683437987,
-151659573325,156074417954,-31268252574,
56280729661,8181904909], -2639651053,
344881152000),
14: (7, 2501928000, [90241897,710986864,-770720657,3501442784,
-6625093363,12630121616,-16802270373,19534438464,
-16802270373,12630121616,-6625093363,3501442784,
-770720657,710986864,90241897], -3740727473,
1275983280000)
}
def newton_cotes(rn, equal=0):
r"""
Return weights and error coefficient for Newton-Cotes integration.
Suppose we have (N+1) samples of f at the positions
x_0, x_1, ..., x_N. Then an N-point Newton-Cotes formula for the
integral between x_0 and x_N is:
:math:`\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i)
+ B_N (\Delta x)^{N+2} f^{N+1} (\xi)`
where :math:`\xi \in [x_0,x_N]`
and :math:`\Delta x = \frac{x_N-x_0}{N}` is the average samples spacing.
If the samples are equally-spaced and N is even, then the error
term is :math:`B_N (\Delta x)^{N+3} f^{N+2}(\xi)`.
Parameters
----------
rn : int
The integer order for equally-spaced data or the relative positions of
the samples with the first sample at 0 and the last at N, where N+1 is
the length of `rn`. N is the order of the Newton-Cotes integration.
equal : int, optional
Set to 1 to enforce equally spaced data.
Returns
-------
an : ndarray
1-D array of weights to apply to the function at the provided sample
positions.
B : float
Error coefficient.
Notes
-----
Normally, the Newton-Cotes rules are used on smaller integration
regions and a composite rule is used to return the total integral.
Examples
--------
Compute the integral of sin(x) in [0, :math:`\pi`]:
>>> from scipy.integrate import newton_cotes
>>> import numpy as np
>>> def f(x):
... return np.sin(x)
>>> a = 0
>>> b = np.pi
>>> exact = 2
>>> for N in [2, 4, 6, 8, 10]:
... x = np.linspace(a, b, N + 1)
... an, B = newton_cotes(N, 1)
... dx = (b - a) / N
... quad = dx * np.sum(an * f(x))
... error = abs(quad - exact)
... print('{:2d} {:10.9f} {:.5e}'.format(N, quad, error))
...
2 2.094395102 9.43951e-02
4 1.998570732 1.42927e-03
6 2.000017814 1.78136e-05
8 1.999999835 1.64725e-07
10 2.000000001 1.14677e-09
"""
try:
N = len(rn)-1
if equal:
rn = np.arange(N+1)
elif np.all(np.diff(rn) == 1):
equal = 1
except Exception:
N = rn
rn = np.arange(N+1)
equal = 1
if equal and N in _builtincoeffs:
na, da, vi, nb, db = _builtincoeffs[N]
an = na * np.array(vi, dtype=float) / da
return an, float(nb)/db
if (rn[0] != 0) or (rn[-1] != N):
raise ValueError("The sample positions must start at 0"
" and end at N")
yi = rn / float(N)
ti = 2 * yi - 1
nvec = np.arange(N+1)
C = ti ** nvec[:, np.newaxis]
Cinv = np.linalg.inv(C)
# improve precision of result
for i in range(2):
Cinv = 2*Cinv - Cinv.dot(C).dot(Cinv)
vec = 2.0 / (nvec[::2]+1)
ai = Cinv[:, ::2].dot(vec) * (N / 2.)
if (N % 2 == 0) and equal:
BN = N/(N+3.)
power = N+2
else:
BN = N/(N+2.)
power = N+1
BN = BN - np.dot(yi**power, ai)
p1 = power+1
fac = power*math.log(N) - gammaln(p1)
fac = math.exp(fac)
return ai, BN*fac
def _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log):
# lazy import to avoid issues with partially-initialized submodule
if not hasattr(qmc_quad, 'qmc'):
from scipy import stats
qmc_quad.stats = stats
else:
stats = qmc_quad.stats
if not callable(func):
message = "`func` must be callable."
raise TypeError(message)
# a, b will be modified, so copy. Oh well if it's copied twice.
a = np.atleast_1d(a).copy()
b = np.atleast_1d(b).copy()
a, b = np.broadcast_arrays(a, b)
dim = a.shape[0]
try:
func((a + b) / 2)
except Exception as e:
message = ("`func` must evaluate the integrand at points within "
"the integration range; e.g. `func( (a + b) / 2)` "
"must return the integrand at the centroid of the "
"integration volume.")
raise ValueError(message) from e
try:
func(np.array([a, b]).T)
vfunc = func
except Exception as e:
message = ("Exception encountered when attempting vectorized call to "
f"`func`: {e}. For better performance, `func` should "
"accept two-dimensional array `x` with shape `(len(a), "
"n_points)` and return an array of the integrand value at "
"each of the `n_points.")
warnings.warn(message, stacklevel=3)
def vfunc(x):
return np.apply_along_axis(func, axis=-1, arr=x)
n_points_int = np.int64(n_points)
if n_points != n_points_int:
message = "`n_points` must be an integer."
raise TypeError(message)
n_estimates_int = np.int64(n_estimates)
if n_estimates != n_estimates_int:
message = "`n_estimates` must be an integer."
raise TypeError(message)
if qrng is None:
qrng = stats.qmc.Halton(dim)
elif not isinstance(qrng, stats.qmc.QMCEngine):
message = "`qrng` must be an instance of scipy.stats.qmc.QMCEngine."
raise TypeError(message)
if qrng.d != a.shape[0]:
message = ("`qrng` must be initialized with dimensionality equal to "
"the number of variables in `a`, i.e., "
"`qrng.random().shape[-1]` must equal `a.shape[0]`.")
raise ValueError(message)
rng_seed = getattr(qrng, 'rng_seed', None)
rng = stats._qmc.check_random_state(rng_seed)
if log not in {True, False}:
message = "`log` must be boolean (`True` or `False`)."
raise TypeError(message)
return (vfunc, a, b, n_points_int, n_estimates_int, qrng, rng, log, stats)
QMCQuadResult = namedtuple('QMCQuadResult', ['integral', 'standard_error'])
def qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None,
log=False):
"""
Compute an integral in N-dimensions using Quasi-Monte Carlo quadrature.
Parameters
----------
func : callable
The integrand. Must accept a single argument ``x``, an array which
specifies the point(s) at which to evaluate the scalar-valued
integrand, and return the value(s) of the integrand.
For efficiency, the function should be vectorized to accept an array of
shape ``(d, n_points)``, where ``d`` is the number of variables (i.e.
the dimensionality of the function domain) and `n_points` is the number
of quadrature points, and return an array of shape ``(n_points,)``,
the integrand at each quadrature point.
a, b : array-like
One-dimensional arrays specifying the lower and upper integration
limits, respectively, of each of the ``d`` variables.
n_estimates, n_points : int, optional
`n_estimates` (default: 8) statistically independent QMC samples, each
of `n_points` (default: 1024) points, will be generated by `qrng`.
The total number of points at which the integrand `func` will be
evaluated is ``n_points * n_estimates``. See Notes for details.
qrng : `~scipy.stats.qmc.QMCEngine`, optional
An instance of the QMCEngine from which to sample QMC points.
The QMCEngine must be initialized to a number of dimensions ``d``
corresponding with the number of variables ``x1, ..., xd`` passed to
`func`.
The provided QMCEngine is used to produce the first integral estimate.
If `n_estimates` is greater than one, additional QMCEngines are
spawned from the first (with scrambling enabled, if it is an option.)
If a QMCEngine is not provided, the default `scipy.stats.qmc.Halton`
will be initialized with the number of dimensions determine from
the length of `a`.
log : boolean, default: False
When set to True, `func` returns the log of the integrand, and
the result object contains the log of the integral.
Returns
-------
result : object
A result object with attributes:
integral : float
The estimate of the integral.
standard_error :
The error estimate. See Notes for interpretation.
Notes
-----
Values of the integrand at each of the `n_points` points of a QMC sample
are used to produce an estimate of the integral. This estimate is drawn
from a population of possible estimates of the integral, the value of
which we obtain depends on the particular points at which the integral
was evaluated. We perform this process `n_estimates` times, each time
evaluating the integrand at different scrambled QMC points, effectively
drawing i.i.d. random samples from the population of integral estimates.
The sample mean :math:`m` of these integral estimates is an
unbiased estimator of the true value of the integral, and the standard
error of the mean :math:`s` of these estimates may be used to generate
confidence intervals using the t distribution with ``n_estimates - 1``
degrees of freedom. Perhaps counter-intuitively, increasing `n_points`
while keeping the total number of function evaluation points
``n_points * n_estimates`` fixed tends to reduce the actual error, whereas
increasing `n_estimates` tends to decrease the error estimate.
Examples
--------
QMC quadrature is particularly useful for computing integrals in higher
dimensions. An example integrand is the probability density function
of a multivariate normal distribution.
>>> import numpy as np
>>> from scipy import stats
>>> dim = 8
>>> mean = np.zeros(dim)
>>> cov = np.eye(dim)
>>> def func(x):
... # `multivariate_normal` expects the _last_ axis to correspond with
... # the dimensionality of the space, so `x` must be transposed
... return stats.multivariate_normal.pdf(x.T, mean, cov)
To compute the integral over the unit hypercube:
>>> from scipy.integrate import qmc_quad
>>> a = np.zeros(dim)
>>> b = np.ones(dim)
>>> rng = np.random.default_rng()
>>> qrng = stats.qmc.Halton(d=dim, seed=rng)
>>> n_estimates = 8
>>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
>>> res.integral, res.standard_error
(0.00018429555666024108, 1.0389431116001344e-07)
A two-sided, 99% confidence interval for the integral may be estimated
as:
>>> t = stats.t(df=n_estimates-1, loc=res.integral,
... scale=res.standard_error)
>>> t.interval(0.99)
(0.0001839319802536469, 0.00018465913306683527)
Indeed, the value reported by `scipy.stats.multivariate_normal` is
within this range.
>>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
0.00018430867675187443
"""
args = _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log)
func, a, b, n_points, n_estimates, qrng, rng, log, stats = args
def sum_product(integrands, dA, log=False):
if log:
return logsumexp(integrands) + np.log(dA)
else:
return np.sum(integrands * dA)
def mean(estimates, log=False):
if log:
return logsumexp(estimates) - np.log(n_estimates)
else:
return np.mean(estimates)
def std(estimates, m=None, ddof=0, log=False):
m = m or mean(estimates, log)
if log:
estimates, m = np.broadcast_arrays(estimates, m)
temp = np.vstack((estimates, m + np.pi * 1j))
diff = logsumexp(temp, axis=0)
return np.real(0.5 * (logsumexp(2 * diff)
- np.log(n_estimates - ddof)))
else:
return np.std(estimates, ddof=ddof)
def sem(estimates, m=None, s=None, log=False):
m = m or mean(estimates, log)
s = s or std(estimates, m, ddof=1, log=log)
if log:
return s - 0.5*np.log(n_estimates)
else:
return s / np.sqrt(n_estimates)
# The sign of the integral depends on the order of the limits. Fix this by
# ensuring that lower bounds are indeed lower and setting sign of resulting
# integral manually
if np.any(a == b):
message = ("A lower limit was equal to an upper limit, so the value "
"of the integral is zero by definition.")
warnings.warn(message, stacklevel=2)
return QMCQuadResult(-np.inf if log else 0, 0)
i_swap = b < a
sign = (-1)**(i_swap.sum(axis=-1)) # odd # of swaps -> negative
a[i_swap], b[i_swap] = b[i_swap], a[i_swap]
A = np.prod(b - a)
dA = A / n_points
estimates = np.zeros(n_estimates)
rngs = _rng_spawn(qrng.rng, n_estimates)
for i in range(n_estimates):
# Generate integral estimate
sample = qrng.random(n_points)
# The rationale for transposing is that this allows users to easily
# unpack `x` into separate variables, if desired. This is consistent
# with the `xx` array passed into the `scipy.integrate.nquad` `func`.
x = stats.qmc.scale(sample, a, b).T # (n_dim, n_points)
integrands = func(x)
estimates[i] = sum_product(integrands, dA, log)
# Get a new, independently-scrambled QRNG for next time
qrng = type(qrng)(seed=rngs[i], **qrng._init_quad)
integral = mean(estimates, log)
standard_error = sem(estimates, m=integral, log=log)
integral = integral + np.pi*1j if (log and sign < 0) else integral*sign
return QMCQuadResult(integral, standard_error)
Reference in New Issue Copy Permalink