Case Study: Creating an Integrator

This case study illustrates how to convert an existing numerical analysis routine, scipy.integrate.quad(), to work with gvar.GVars.

The Problem

We want a Python code that can evaluate one dimensional integrals such as

where any of the integration limits or are gvar.GVars and is an arbitrary function coded as a Python function.

One approach is to implement an integration function directly in Python, as then it is likely to work just as well for gvar.GVars as for floats. For example, the code

>>> import gvar as gv
>>> import numpy as np
>>>
>>> def trap_integral(f, interval, n=100):
...     """ Estimate integral of f(x) on interval=(a,b) using the Trapezoidal Rule. """
...     a, b = interval
...     x = a + (b - a) * np.linspace(0, 1., n+1)
...     fx = np.array([f(xi) for xi in x])
...     I =  np.sum(fx[:-1] + fx[1:]) * (b - a) / (2. * n)
...     return I
...
>>> A = gv.gvar(2, 0.1)
>>> K = gv.gvar(1, 0.11)
>>> D = gv.gvar(1., 0.4)
>>>
>>> def f(x):
...     return A * np.cos(K * x**2 + D) ** 2
...
>>> a = gv.gvar(0, 0.1)
>>> b = gv.gvar(4, 0.1)
>>> Itrap = trap_integral(f, (a, b), n=100)
>>> print(f'Itrap = {Itrap:#P}')
Itrap = 3.45 ± 0.32

estimates the integral of f(x) over the interval between 0 ± 0.1 and 4 ± 0.1 using the Trapezoidal Rule.

This code is simple because we are using one of the simplest numerical estimates of the integral. A general purpose integrators needs a much more robust algorithm. For example, trap_integral fails badly when applied to a much more singular function:

>>> def g(x):
...    return A * x /(K * x**2 + 1e-6)
...
>>> Itrap_g = trap_integral(g, (a, b), n=100)
>>> print(f'Itrap_g = {Itrap_g:#P}')
Itrap_g = 10.3633 ± 4.0e+03

The correct answer is 16.6 ± 1.9. We need a much larger number of integrand samples n (100x larger) to get reasonable results.

Leveraging Existing Code

Coding a more robust integrator is complicated and time consuming. A better strategy is, if possible, to build on existing libraries. Here we will use integrators from the scipy.integrate module.

The integral is a function of its endpoints and of any parameters buried in the definition of the function : where and for are the parameters implicit in the integrand (e.g., A and K in the examples above). We want an integrator that works when any of these parameters is replaced by a gvar.GVar.

We can do this using gvar.gvar_function(p, I, dI_dp) where p is an array of the gvar.GVar-valued parameters, I is the integral evaluated with these parameters replaced by their mean values, and dI_dp is the array of derivatives of I with respect to each of these parameters — — again evaluated with their mean values.

The integral I (with the parameters replaced by their mean values) can be evaluated using standard routines as no gvar.GVars are involved. The derivatives with respect to the endpoints are also easily evaluated:

The derivatives with respect to the function parameters involve different integrals, which again can be evaluated using standard routines:

In the following code we use the integrators quad(...) and quad_vec(...) from scipy.integrate to evaluate the integrals needed to calculate I and elements of dI_dp, respectively:

import scipy.integrate

def integral(f, interval, tol=1e-8):
    """ GVar-compatible integrator """
    a, b = interval

    # collect GVar-valued parameters
    p = []
    dI_dp = []
    if isinstance(a, gv.GVar):
        p += [a]
        dI_dp += [-f(a).mean]
        a = a.mean
    if isinstance(b, gv.GVar):
        p += [b]
        dI_dp += [f(b).mean]
        b = b.mean

    # evaluate integral I of f(x).mean
    sum_fx = [0]
    def fmean(x):
        fx = f(x)
        if isinstance(fx, gv.GVar):
            sum_fx[0] += f(x)
            return f(x).mean
        else:
            return fx
    I = scipy.integrate.quad(fmean, a, b, epsrel=tol)[0]

    # parameters from the integrand
    pf = gv.dependencies(sum_fx[0], all=True)

    # evaluate dI/dpf
    if len(pf) > 0:
        # vector-valued integrand returns df(x)/dpf
        def df_dpf(x):
            fx = f(x)
            if isinstance(fx, gv.GVar):
                return fx.deriv(pf)
            else:
                return np.array(len(pf) * [0.0])

        # integrate df/dpf to obtain dI/dpf
        dI_dpf = scipy.integrate.quad_vec(df_dpf, a, b, epsrel=tol)[0]

        # combine with other parameters, if any
        p += list(pf)
        dI_dp += list(dI_dpf)

    return gv.gvar_function(p, I, dI_dp) if len(p) > 0 else I

A key ingredient of this code is the use of gvar.dependencies() to obtain an array pf of the gvar.GVar-valued parameters implicit in the integrand f(x). This is done without knowing anything about f(x) beyond the sum sum_fx[0] of its values at all the integration points used to calculate I. Given parameters pf[i], the derivatives of f(x) with respect to those parameters are obtained using f(x).deriv(pf) (see the documentation for gvar.GVar.deriv()).

This new integrator works well with the first example above and gives the same result:

>>> I = integral(f, (a, b))
>>> print(f'I = {I:#P}')
I = 3.45 ± 0.32

It also works well with the much more singular integrand g(x):

>>> I_g = integral(g, (a, b))
>>> print(f'I_g = {I_g:#P}')
I_g = 16.6 ± 1.9

gvar comes with a different integrator, gvar.ode.integral(), that gives the same results with similar performance: for example,

>>> Iode = gv.ode.integral(f, (a, b))
>>> print(f'Iode = {Iode:#P}')
Iode = 3.45 ± 0.32
>>> Iode_g = gv.ode.integral(g, (a, b))
>>> print(f'Iode_g = {Iode_g:#P}')
Iode_g = 16.6 ± 1.9

We can generate error budgets for each of the integral estimates to see where the final uncertainties come from:

>>> inputs = dict(a=a, b=b, A=A, K=K, D=D)
>>> outputs = dict(I=I, Iode=Iode, Itrap=Itrap)
>>> print(gv.fmt_errorbudget(inputs=inputs, outputs=outputs))
Partial % Errors:
                   I      Iode     Itrap
----------------------------------------
        a:      1.69      1.69      1.69
        b:      0.44      0.44      0.52
        A:      5.00      5.00      5.00
        K:      4.53      4.53      4.35
        D:      6.29      6.29      6.25
----------------------------------------
    total:      9.39      9.39      9.28

As expected the different methods are in good agreement (the Trapezoidal Rule gives slightly different results because n is a bit too small).