Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/integrate/_quad_vec.py: 13%
264 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1import sys
2import copy
3import heapq
4import collections
5import functools
7import numpy as np
9from scipy._lib._util import MapWrapper, _FunctionWrapper
12class LRUDict(collections.OrderedDict):
13 def __init__(self, max_size):
14 self.__max_size = max_size
16 def __setitem__(self, key, value):
17 existing_key = (key in self)
18 super().__setitem__(key, value)
19 if existing_key:
20 self.move_to_end(key)
21 elif len(self) > self.__max_size:
22 self.popitem(last=False)
24 def update(self, other):
25 # Not needed below
26 raise NotImplementedError()
29class SemiInfiniteFunc:
30 """
31 Argument transform from (start, +-oo) to (0, 1)
32 """
33 def __init__(self, func, start, infty):
34 self._func = func
35 self._start = start
36 self._sgn = -1 if infty < 0 else 1
38 # Overflow threshold for the 1/t**2 factor
39 self._tmin = sys.float_info.min**0.5
41 def get_t(self, x):
42 z = self._sgn * (x - self._start) + 1
43 if z == 0:
44 # Can happen only if point not in range
45 return np.inf
46 return 1 / z
48 def __call__(self, t):
49 if t < self._tmin:
50 return 0.0
51 else:
52 x = self._start + self._sgn * (1 - t) / t
53 f = self._func(x)
54 return self._sgn * (f / t) / t
57class DoubleInfiniteFunc:
58 """
59 Argument transform from (-oo, oo) to (-1, 1)
60 """
61 def __init__(self, func):
62 self._func = func
64 # Overflow threshold for the 1/t**2 factor
65 self._tmin = sys.float_info.min**0.5
67 def get_t(self, x):
68 s = -1 if x < 0 else 1
69 return s / (abs(x) + 1)
71 def __call__(self, t):
72 if abs(t) < self._tmin:
73 return 0.0
74 else:
75 x = (1 - abs(t)) / t
76 f = self._func(x)
77 return (f / t) / t
80def _max_norm(x):
81 return np.amax(abs(x))
84def _get_sizeof(obj):
85 try:
86 return sys.getsizeof(obj)
87 except TypeError:
88 # occurs on pypy
89 if hasattr(obj, '__sizeof__'):
90 return int(obj.__sizeof__())
91 return 64
94class _Bunch:
95 def __init__(self, **kwargs):
96 self.__keys = kwargs.keys()
97 self.__dict__.update(**kwargs)
99 def __repr__(self):
100 return "_Bunch({})".format(", ".join("{}={}".format(k, repr(self.__dict__[k]))
101 for k in self.__keys))
104def quad_vec(f, a, b, epsabs=1e-200, epsrel=1e-8, norm='2', cache_size=100e6, limit=10000,
105 workers=1, points=None, quadrature=None, full_output=False,
106 *, args=()):
107 r"""Adaptive integration of a vector-valued function.
109 Parameters
110 ----------
111 f : callable
112 Vector-valued function f(x) to integrate.
113 a : float
114 Initial point.
115 b : float
116 Final point.
117 epsabs : float, optional
118 Absolute tolerance.
119 epsrel : float, optional
120 Relative tolerance.
121 norm : {'max', '2'}, optional
122 Vector norm to use for error estimation.
123 cache_size : int, optional
124 Number of bytes to use for memoization.
125 limit : float or int, optional
126 An upper bound on the number of subintervals used in the adaptive
127 algorithm.
128 workers : int or map-like callable, optional
129 If `workers` is an integer, part of the computation is done in
130 parallel subdivided to this many tasks (using
131 :class:`python:multiprocessing.pool.Pool`).
132 Supply `-1` to use all cores available to the Process.
133 Alternatively, supply a map-like callable, such as
134 :meth:`python:multiprocessing.pool.Pool.map` for evaluating the
135 population in parallel.
136 This evaluation is carried out as ``workers(func, iterable)``.
137 points : list, optional
138 List of additional breakpoints.
139 quadrature : {'gk21', 'gk15', 'trapezoid'}, optional
140 Quadrature rule to use on subintervals.
141 Options: 'gk21' (Gauss-Kronrod 21-point rule),
142 'gk15' (Gauss-Kronrod 15-point rule),
143 'trapezoid' (composite trapezoid rule).
144 Default: 'gk21' for finite intervals and 'gk15' for (semi-)infinite
145 full_output : bool, optional
146 Return an additional ``info`` dictionary.
147 args : tuple, optional
148 Extra arguments to pass to function, if any.
150 .. versionadded:: 1.8.0
152 Returns
153 -------
154 res : {float, array-like}
155 Estimate for the result
156 err : float
157 Error estimate for the result in the given norm
158 info : dict
159 Returned only when ``full_output=True``.
160 Info dictionary. Is an object with the attributes:
162 success : bool
163 Whether integration reached target precision.
164 status : int
165 Indicator for convergence, success (0),
166 failure (1), and failure due to rounding error (2).
167 neval : int
168 Number of function evaluations.
169 intervals : ndarray, shape (num_intervals, 2)
170 Start and end points of subdivision intervals.
171 integrals : ndarray, shape (num_intervals, ...)
172 Integral for each interval.
173 Note that at most ``cache_size`` values are recorded,
174 and the array may contains *nan* for missing items.
175 errors : ndarray, shape (num_intervals,)
176 Estimated integration error for each interval.
178 Notes
179 -----
180 The algorithm mainly follows the implementation of QUADPACK's
181 DQAG* algorithms, implementing global error control and adaptive
182 subdivision.
184 The algorithm here has some differences to the QUADPACK approach:
186 Instead of subdividing one interval at a time, the algorithm
187 subdivides N intervals with largest errors at once. This enables
188 (partial) parallelization of the integration.
190 The logic of subdividing "next largest" intervals first is then
191 not implemented, and we rely on the above extension to avoid
192 concentrating on "small" intervals only.
194 The Wynn epsilon table extrapolation is not used (QUADPACK uses it
195 for infinite intervals). This is because the algorithm here is
196 supposed to work on vector-valued functions, in an user-specified
197 norm, and the extension of the epsilon algorithm to this case does
198 not appear to be widely agreed. For max-norm, using elementwise
199 Wynn epsilon could be possible, but we do not do this here with
200 the hope that the epsilon extrapolation is mainly useful in
201 special cases.
203 References
204 ----------
205 [1] R. Piessens, E. de Doncker, QUADPACK (1983).
207 Examples
208 --------
209 We can compute integrations of a vector-valued function:
211 >>> from scipy.integrate import quad_vec
212 >>> import numpy as np
213 >>> import matplotlib.pyplot as plt
214 >>> alpha = np.linspace(0.0, 2.0, num=30)
215 >>> f = lambda x: x**alpha
216 >>> x0, x1 = 0, 2
217 >>> y, err = quad_vec(f, x0, x1)
218 >>> plt.plot(alpha, y)
219 >>> plt.xlabel(r"$\alpha$")
220 >>> plt.ylabel(r"$\int_{0}^{2} x^\alpha dx$")
221 >>> plt.show()
223 """
224 a = float(a)
225 b = float(b)
227 if args:
228 if not isinstance(args, tuple):
229 args = (args,)
231 # create a wrapped function to allow the use of map and Pool.map
232 f = _FunctionWrapper(f, args)
234 # Use simple transformations to deal with integrals over infinite
235 # intervals.
236 kwargs = dict(epsabs=epsabs,
237 epsrel=epsrel,
238 norm=norm,
239 cache_size=cache_size,
240 limit=limit,
241 workers=workers,
242 points=points,
243 quadrature='gk15' if quadrature is None else quadrature,
244 full_output=full_output)
245 if np.isfinite(a) and np.isinf(b):
246 f2 = SemiInfiniteFunc(f, start=a, infty=b)
247 if points is not None:
248 kwargs['points'] = tuple(f2.get_t(xp) for xp in points)
249 return quad_vec(f2, 0, 1, **kwargs)
250 elif np.isfinite(b) and np.isinf(a):
251 f2 = SemiInfiniteFunc(f, start=b, infty=a)
252 if points is not None:
253 kwargs['points'] = tuple(f2.get_t(xp) for xp in points)
254 res = quad_vec(f2, 0, 1, **kwargs)
255 return (-res[0],) + res[1:]
256 elif np.isinf(a) and np.isinf(b):
257 sgn = -1 if b < a else 1
259 # NB. explicitly split integral at t=0, which separates
260 # the positive and negative sides
261 f2 = DoubleInfiniteFunc(f)
262 if points is not None:
263 kwargs['points'] = (0,) + tuple(f2.get_t(xp) for xp in points)
264 else:
265 kwargs['points'] = (0,)
267 if a != b:
268 res = quad_vec(f2, -1, 1, **kwargs)
269 else:
270 res = quad_vec(f2, 1, 1, **kwargs)
272 return (res[0]*sgn,) + res[1:]
273 elif not (np.isfinite(a) and np.isfinite(b)):
274 raise ValueError("invalid integration bounds a={}, b={}".format(a, b))
276 norm_funcs = {
277 None: _max_norm,
278 'max': _max_norm,
279 '2': np.linalg.norm
280 }
281 if callable(norm):
282 norm_func = norm
283 else:
284 norm_func = norm_funcs[norm]
286 parallel_count = 128
287 min_intervals = 2
289 try:
290 _quadrature = {None: _quadrature_gk21,
291 'gk21': _quadrature_gk21,
292 'gk15': _quadrature_gk15,
293 'trapz': _quadrature_trapezoid, # alias for backcompat
294 'trapezoid': _quadrature_trapezoid}[quadrature]
295 except KeyError as e:
296 raise ValueError("unknown quadrature {!r}".format(quadrature)) from e
298 # Initial interval set
299 if points is None:
300 initial_intervals = [(a, b)]
301 else:
302 prev = a
303 initial_intervals = []
304 for p in sorted(points):
305 p = float(p)
306 if not (a < p < b) or p == prev:
307 continue
308 initial_intervals.append((prev, p))
309 prev = p
310 initial_intervals.append((prev, b))
312 global_integral = None
313 global_error = None
314 rounding_error = None
315 interval_cache = None
316 intervals = []
317 neval = 0
319 for x1, x2 in initial_intervals:
320 ig, err, rnd = _quadrature(x1, x2, f, norm_func)
321 neval += _quadrature.num_eval
323 if global_integral is None:
324 if isinstance(ig, (float, complex)):
325 # Specialize for scalars
326 if norm_func in (_max_norm, np.linalg.norm):
327 norm_func = abs
329 global_integral = ig
330 global_error = float(err)
331 rounding_error = float(rnd)
333 cache_count = cache_size // _get_sizeof(ig)
334 interval_cache = LRUDict(cache_count)
335 else:
336 global_integral += ig
337 global_error += err
338 rounding_error += rnd
340 interval_cache[(x1, x2)] = copy.copy(ig)
341 intervals.append((-err, x1, x2))
343 heapq.heapify(intervals)
345 CONVERGED = 0
346 NOT_CONVERGED = 1
347 ROUNDING_ERROR = 2
348 NOT_A_NUMBER = 3
350 status_msg = {
351 CONVERGED: "Target precision reached.",
352 NOT_CONVERGED: "Target precision not reached.",
353 ROUNDING_ERROR: "Target precision could not be reached due to rounding error.",
354 NOT_A_NUMBER: "Non-finite values encountered."
355 }
357 # Process intervals
358 with MapWrapper(workers) as mapwrapper:
359 ier = NOT_CONVERGED
361 while intervals and len(intervals) < limit:
362 # Select intervals with largest errors for subdivision
363 tol = max(epsabs, epsrel*norm_func(global_integral))
365 to_process = []
366 err_sum = 0
368 for j in range(parallel_count):
369 if not intervals:
370 break
372 if j > 0 and err_sum > global_error - tol/8:
373 # avoid unnecessary parallel splitting
374 break
376 interval = heapq.heappop(intervals)
378 neg_old_err, a, b = interval
379 old_int = interval_cache.pop((a, b), None)
380 to_process.append(((-neg_old_err, a, b, old_int), f, norm_func, _quadrature))
381 err_sum += -neg_old_err
383 # Subdivide intervals
384 for dint, derr, dround_err, subint, dneval in mapwrapper(_subdivide_interval, to_process):
385 neval += dneval
386 global_integral += dint
387 global_error += derr
388 rounding_error += dround_err
389 for x in subint:
390 x1, x2, ig, err = x
391 interval_cache[(x1, x2)] = ig
392 heapq.heappush(intervals, (-err, x1, x2))
394 # Termination check
395 if len(intervals) >= min_intervals:
396 tol = max(epsabs, epsrel*norm_func(global_integral))
397 if global_error < tol/8:
398 ier = CONVERGED
399 break
400 if global_error < rounding_error:
401 ier = ROUNDING_ERROR
402 break
404 if not (np.isfinite(global_error) and np.isfinite(rounding_error)):
405 ier = NOT_A_NUMBER
406 break
408 res = global_integral
409 err = global_error + rounding_error
411 if full_output:
412 res_arr = np.asarray(res)
413 dummy = np.full(res_arr.shape, np.nan, dtype=res_arr.dtype)
414 integrals = np.array([interval_cache.get((z[1], z[2]), dummy)
415 for z in intervals], dtype=res_arr.dtype)
416 errors = np.array([-z[0] for z in intervals])
417 intervals = np.array([[z[1], z[2]] for z in intervals])
419 info = _Bunch(neval=neval,
420 success=(ier == CONVERGED),
421 status=ier,
422 message=status_msg[ier],
423 intervals=intervals,
424 integrals=integrals,
425 errors=errors)
426 return (res, err, info)
427 else:
428 return (res, err)
431def _subdivide_interval(args):
432 interval, f, norm_func, _quadrature = args
433 old_err, a, b, old_int = interval
435 c = 0.5 * (a + b)
437 # Left-hand side
438 if getattr(_quadrature, 'cache_size', 0) > 0:
439 f = functools.lru_cache(_quadrature.cache_size)(f)
441 s1, err1, round1 = _quadrature(a, c, f, norm_func)
442 dneval = _quadrature.num_eval
443 s2, err2, round2 = _quadrature(c, b, f, norm_func)
444 dneval += _quadrature.num_eval
445 if old_int is None:
446 old_int, _, _ = _quadrature(a, b, f, norm_func)
447 dneval += _quadrature.num_eval
449 if getattr(_quadrature, 'cache_size', 0) > 0:
450 dneval = f.cache_info().misses
452 dint = s1 + s2 - old_int
453 derr = err1 + err2 - old_err
454 dround_err = round1 + round2
456 subintervals = ((a, c, s1, err1), (c, b, s2, err2))
457 return dint, derr, dround_err, subintervals, dneval
460def _quadrature_trapezoid(x1, x2, f, norm_func):
461 """
462 Composite trapezoid quadrature
463 """
464 x3 = 0.5*(x1 + x2)
465 f1 = f(x1)
466 f2 = f(x2)
467 f3 = f(x3)
469 s2 = 0.25 * (x2 - x1) * (f1 + 2*f3 + f2)
471 round_err = 0.25 * abs(x2 - x1) * (float(norm_func(f1))
472 + 2*float(norm_func(f3))
473 + float(norm_func(f2))) * 2e-16
475 s1 = 0.5 * (x2 - x1) * (f1 + f2)
476 err = 1/3 * float(norm_func(s1 - s2))
477 return s2, err, round_err
480_quadrature_trapezoid.cache_size = 3 * 3
481_quadrature_trapezoid.num_eval = 3
484def _quadrature_gk(a, b, f, norm_func, x, w, v):
485 """
486 Generic Gauss-Kronrod quadrature
487 """
489 fv = [0.0]*len(x)
491 c = 0.5 * (a + b)
492 h = 0.5 * (b - a)
494 # Gauss-Kronrod
495 s_k = 0.0
496 s_k_abs = 0.0
497 for i in range(len(x)):
498 ff = f(c + h*x[i])
499 fv[i] = ff
501 vv = v[i]
503 # \int f(x)
504 s_k += vv * ff
505 # \int |f(x)|
506 s_k_abs += vv * abs(ff)
508 # Gauss
509 s_g = 0.0
510 for i in range(len(w)):
511 s_g += w[i] * fv[2*i + 1]
513 # Quadrature of abs-deviation from average
514 s_k_dabs = 0.0
515 y0 = s_k / 2.0
516 for i in range(len(x)):
517 # \int |f(x) - y0|
518 s_k_dabs += v[i] * abs(fv[i] - y0)
520 # Use similar error estimation as quadpack
521 err = float(norm_func((s_k - s_g) * h))
522 dabs = float(norm_func(s_k_dabs * h))
523 if dabs != 0 and err != 0:
524 err = dabs * min(1.0, (200 * err / dabs)**1.5)
526 eps = sys.float_info.epsilon
527 round_err = float(norm_func(50 * eps * h * s_k_abs))
529 if round_err > sys.float_info.min:
530 err = max(err, round_err)
532 return h * s_k, err, round_err
535def _quadrature_gk21(a, b, f, norm_func):
536 """
537 Gauss-Kronrod 21 quadrature with error estimate
538 """
539 # Gauss-Kronrod points
540 x = (0.995657163025808080735527280689003,
541 0.973906528517171720077964012084452,
542 0.930157491355708226001207180059508,
543 0.865063366688984510732096688423493,
544 0.780817726586416897063717578345042,
545 0.679409568299024406234327365114874,
546 0.562757134668604683339000099272694,
547 0.433395394129247190799265943165784,
548 0.294392862701460198131126603103866,
549 0.148874338981631210884826001129720,
550 0,
551 -0.148874338981631210884826001129720,
552 -0.294392862701460198131126603103866,
553 -0.433395394129247190799265943165784,
554 -0.562757134668604683339000099272694,
555 -0.679409568299024406234327365114874,
556 -0.780817726586416897063717578345042,
557 -0.865063366688984510732096688423493,
558 -0.930157491355708226001207180059508,
559 -0.973906528517171720077964012084452,
560 -0.995657163025808080735527280689003)
562 # 10-point weights
563 w = (0.066671344308688137593568809893332,
564 0.149451349150580593145776339657697,
565 0.219086362515982043995534934228163,
566 0.269266719309996355091226921569469,
567 0.295524224714752870173892994651338,
568 0.295524224714752870173892994651338,
569 0.269266719309996355091226921569469,
570 0.219086362515982043995534934228163,
571 0.149451349150580593145776339657697,
572 0.066671344308688137593568809893332)
574 # 21-point weights
575 v = (0.011694638867371874278064396062192,
576 0.032558162307964727478818972459390,
577 0.054755896574351996031381300244580,
578 0.075039674810919952767043140916190,
579 0.093125454583697605535065465083366,
580 0.109387158802297641899210590325805,
581 0.123491976262065851077958109831074,
582 0.134709217311473325928054001771707,
583 0.142775938577060080797094273138717,
584 0.147739104901338491374841515972068,
585 0.149445554002916905664936468389821,
586 0.147739104901338491374841515972068,
587 0.142775938577060080797094273138717,
588 0.134709217311473325928054001771707,
589 0.123491976262065851077958109831074,
590 0.109387158802297641899210590325805,
591 0.093125454583697605535065465083366,
592 0.075039674810919952767043140916190,
593 0.054755896574351996031381300244580,
594 0.032558162307964727478818972459390,
595 0.011694638867371874278064396062192)
597 return _quadrature_gk(a, b, f, norm_func, x, w, v)
600_quadrature_gk21.num_eval = 21
603def _quadrature_gk15(a, b, f, norm_func):
604 """
605 Gauss-Kronrod 15 quadrature with error estimate
606 """
607 # Gauss-Kronrod points
608 x = (0.991455371120812639206854697526329,
609 0.949107912342758524526189684047851,
610 0.864864423359769072789712788640926,
611 0.741531185599394439863864773280788,
612 0.586087235467691130294144838258730,
613 0.405845151377397166906606412076961,
614 0.207784955007898467600689403773245,
615 0.000000000000000000000000000000000,
616 -0.207784955007898467600689403773245,
617 -0.405845151377397166906606412076961,
618 -0.586087235467691130294144838258730,
619 -0.741531185599394439863864773280788,
620 -0.864864423359769072789712788640926,
621 -0.949107912342758524526189684047851,
622 -0.991455371120812639206854697526329)
624 # 7-point weights
625 w = (0.129484966168869693270611432679082,
626 0.279705391489276667901467771423780,
627 0.381830050505118944950369775488975,
628 0.417959183673469387755102040816327,
629 0.381830050505118944950369775488975,
630 0.279705391489276667901467771423780,
631 0.129484966168869693270611432679082)
633 # 15-point weights
634 v = (0.022935322010529224963732008058970,
635 0.063092092629978553290700663189204,
636 0.104790010322250183839876322541518,
637 0.140653259715525918745189590510238,
638 0.169004726639267902826583426598550,
639 0.190350578064785409913256402421014,
640 0.204432940075298892414161999234649,
641 0.209482141084727828012999174891714,
642 0.204432940075298892414161999234649,
643 0.190350578064785409913256402421014,
644 0.169004726639267902826583426598550,
645 0.140653259715525918745189590510238,
646 0.104790010322250183839876322541518,
647 0.063092092629978553290700663189204,
648 0.022935322010529224963732008058970)
650 return _quadrature_gk(a, b, f, norm_func, x, w, v)
653_quadrature_gk15.num_eval = 15