code: def taylor_approx(hrf2decompose, time=None, delta=None):
      """A Taylor series approximation of an HRF shifted by times `delta`
  
      Returns original HRF and gradient of HRF
  
      Parameters
      ----------
      hrf2decompose : sympy expression
          An expression that can be lambdified as a function of 't'. This
          is the HRF to be expanded in PCA
      time : None or np.ndarray, optional
          None gives default value of np.linspace(-15,50,3251) chosen to
          match fMRIstat implementation.  This corresponds to a time
          interval of 0.02.  Presumed to be equally spaced.
      delta : None or np.ndarray, optional
          None results in default value of np.arange(-4.5, 4.6, 0.1)
          chosen to match fMRIstat implementation.
  
      Returns
      -------
      hrf : [sympy expressions]
          Sequence length 2 comprising (`hrf2decompose`, ``dhrf``) where
          ``dhrf`` is the first derivative of `hrf2decompose`.
      approx :
          TODO
  
      References
      ----------
      Liao, C.H., Worsley, K.J., Poline, J-B., Aston, J.A.D., Duncan, G.H.,
      Evans, A.C. (2002). \'Estimating the delay of the response in fMRI
      data.\' NeuroImage, 16:593-606.
      """
      if time is None:
          time = np.linspace(-15, 50, 3251)
      dt = time[1] - time[0]
      if delta is None:
          delta = np.arange(-4.5, 4.6, 0.1)
      # make numerical implementation from hrf function and symbol t.
      # hrft returns function values when called with values for time as
      # input.
      hrft = lambdify_t(hrf2decompose(T))
      # interpolator for negative gradient of hrf
      dhrft = Interp1dNumeric(
          time,
          -np.gradient(hrft(time), dt),
          bounds_error=False,
          fill_value=0.0,
      )
      dhrft.y *= 2
      # Create stack of time-shifted HRFs.  Time varies over row, delta
      # over column.
      ts_hrf_vals = np.array([hrft(time - d) for d in delta]).T
      # hrf, dhrf
      W = np.array([hrft(time), dhrft(time)]).T
      # regress hrf, dhrf at times against stack of time-shifted hrfs
      WH = np.dot(npl.pinv(W), ts_hrf_vals)
      # put these into interpolators to get estimated coefficients for any
      # value of delta
      coef = [
          Interp1dNumeric(delta, w, bounds_error=False, fill_value=0.0)
          for w in WH
      ]
  
      def approx(time, delta):
          value = coef[0](delta) * hrft(time) + coef[1](delta) * dhrft(
              time
          )
          return value
  
      approx.coef = coef
      approx.components = [hrft, dhrft]
      (
          approx.theta,
          approx.inverse,
          approx.dinverse,
          approx.forward,
          approx.dforward,
      ) = invertR(delta, approx.coef)
      dhrf = implemented_function(f"d{str(hrf2decompose)}", dhrft)
      return [hrf2decompose, dhrf], approx
  
file_path: nipy/nipy/modalities/fmri/fmristat/hrf.py
imports: import numpy
  import numpy.linalg
  from sympy.utilities.lambdify import implemented_function
  from ..None import hrf
  from ..utils import Interp1dNumeric, T, lambdify_t
  from .invert import invertR
label_desc: Use image node_id 2 for calling a global function with example usage: taylor_approx(hrf2decompose, time, delta) and returns: list,  approx, value