'Compute a definite integral over a piecewise polynomial.
Parameters
ranges : ndim-tuple of 2-tuples float
Sequence of lower and upper bounds for each dimension,
``[(a[0], b[0]), ..., (a[ndim-1], b[ndim-1])]``
extrapolate : bool, optional
Whether to extrapolate to out-of-bounds points based on first
and last intervals, or to return NaNs.
Returns
ig : array_like
Definite integral of the piecewise polynomial over
[a[0], b[0]] x ... x [a[ndim-1], b[ndim-1]]'

def integrate(self, ranges, extrapolate=None):

    ndim = len(self.x)
    if (extrapolate is None):
        extrapolate = self.extrapolate
    else:
        extrapolate = bool(extrapolate)
    if ((not hasattr(ranges, '__len__')) or (len(ranges) != ndim)):
        raise ValueError('Range    not    a    sequence    of    correct    length')
    self._ensure_c_contiguous()
    c = self.c
    for (n, (a, b)) in enumerate(ranges):
        swap = list(range(c.ndim))
        swap.insert(1, swap[(ndim - n)])
        del swap[((ndim - n) + 1)]
        c = c.transpose(swap)
        p = PPoly.construct_fast(c, self.x[n], extrapolate=extrapolate)
        out = p.integrate(a, b, extrapolate=extrapolate)
        c = out.reshape(c.shape[2:])
    return c
