code: def natural_spline(tvals, knots=None, order=3, intercept=True):
      """Design matrix with columns given by a natural spline order `order`
  
      Return design matrix with natural splines with knots `knots`, order
      `order`.  If `intercept` == True (the default), add constant column.
  
      Parameters
      ----------
      tvals : np.array
          Time values
      knots : None or sequence, optional
         Sequence of float.  Default None (same as empty list)
      order : int, optional
         Order of the spline. Defaults to a cubic (==3)
      intercept : bool, optional
         If True, include a constant function in the natural
         spline. Default is False
  
      Returns
      -------
      X : np.ndarray
  
      Examples
      --------
      >>> tvals = np.linspace(0,50,101)
      >>> drift = natural_spline(tvals, knots=[10,20,30,40])
      >>> drift.shape
      (101, 8)
      """
      tvals = make_recarray(tvals, ["t"])
      t = Term("t")
      f = formulae.natural_spline(
          t, knots=knots, order=order, intercept=intercept
      )
      return f.design(tvals, return_float=True)
  
imports: import itertools
  from functools import reduce
  from operator import mul
  import numpy
  from nipy.algorithms.statistics.formula import formulae
  from nipy.algorithms.statistics.formula.formulae import Factor, Formula, Term, make_recarray
  from nipy.algorithms.statistics.utils import combinations
  from .hrf import glover
  from .utils import T, blocks, convolve_functions, events
  from .utils import fourier_basis
label_desc: Use image node_id 2 for calling a global function with example usage: natural_spline(tvals, knots, order, intercept) and returns: f