code: def cumtrapz(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.
          Typically this value should be 0.  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.
  
      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
      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 matplotlib.pyplot as plt
  
      >>> x = np.linspace(-2, 2, num=20)
      >>> y = x
      >>> y_int = integrate.cumtrapz(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 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 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
  
imports: from __future__ import division, print_function, absolute_import
  import numpy
  import math
  import warnings
  from numpy import trapz
  from scipy.special import roots_legendre
  from scipy.special import gammaln
  from scipy._lib.six import xrange
label_desc: Use image node_id 6 for calling a global function with example usage: cumtrapz(y, x, dx, axis, initial) and returns: res