code: def protocol(recarr, design_type, *hrfs):
      """Create an object that can evaluate the FIAC
  
      Subclass of formulae.Formula, but not necessary.
  
      Parameters
      ----------
      recarr : (N,) structured array
         with fields 'time' and 'event'
      design_type : str
         one of ['event', 'block'].  Handles how the 'begin' term is
         handled.  For 'block', the first event of each block is put in
         this group. For the 'event', only the first event is put in this
         group. The 'begin' events are convolved with hrf.glover.
      hrfs: symbolic HRFs
         Each event type ('SSt_SSp','SSt_DSp','DSt_SSp','DSt_DSp') is
         convolved with each of these HRFs in order.
  
      Returns
      -------
      f: Formula
         Formula for constructing design matrices.
      contrasts : dict
         Dictionary of the contrasts of the experiment.
      """
      event_types = np.unique(recarr["event"])
      N = recarr.size
      if design_type == "block":
          keep = np.not_equal((np.arange(N)) % 6, 0)
      else:
          keep = np.greater(np.arange(N), 0)
      # This first frame was used to model out a potentially
      # 'bad' first frame....
      _begin = recarr["time"][~keep]
  
      termdict = {}
      termdict["begin"] = utils.define(
          "begin", utils.events(_begin, f=hrf.glover)
      )
      drift = formulae.natural_spline(
          utils.T, knots=[N_ROWS / 2.0 + 1.25], intercept=True
      )
      for i, t in enumerate(drift.terms):
          termdict["drift%d" % i] = t
      # After removing the first frame, keep the remaining
      # events and times
      times = recarr["time"][keep]
      events = recarr["event"][keep]
  
      # Now, specify the experimental conditions.  This creates expressions named
      # SSt_SSp0, SSt_SSp1, etc.  with one expression for each (eventtype, hrf)
      # pair
      for v in event_types:
          k = np.array([events[i] == v for i in range(times.shape[0])])
          for l, h in enumerate(hrfs):
              # Make sure event type is a string (not byte string)
              term_name = "%s%d" % (_to_str(v), l)
              termdict[term_name] = utils.define(
                  term_name, utils.events(times[k], f=h)
              )
      f = formulae.Formula(list(termdict.values()))
      Tcontrasts = {}
      Tcontrasts["average"] = (
          termdict["SSt_SSp0"]
          + termdict["SSt_DSp0"]
          + termdict["DSt_SSp0"]
          + termdict["DSt_DSp0"]
      ) / 4.0
      Tcontrasts["speaker"] = (
          termdict["SSt_DSp0"]
          - termdict["SSt_SSp0"]
          + termdict["DSt_DSp0"]
          - termdict["DSt_SSp0"]
      ) * 0.5
      Tcontrasts["sentence"] = (
          termdict["DSt_DSp0"]
          + termdict["DSt_SSp0"]
          - termdict["SSt_DSp0"]
          - termdict["SSt_SSp0"]
      ) * 0.5
      Tcontrasts["interaction"] = (
          termdict["SSt_SSp0"]
          - termdict["SSt_DSp0"]
          - termdict["DSt_SSp0"]
          + termdict["DSt_DSp0"]
      )
      # Ftest
      Fcontrasts = {}
      Fcontrasts["overall1"] = formulae.Formula(
          list(Tcontrasts.values())
      )
  
      return f, Tcontrasts, Fcontrasts
  
file_path: nipy/nipy/modalities/fmri/fmristat/tests/test_FIAC.py
imports: import numpy
  from nipy.algorithms.statistics.formula import formulae
  from nipy.algorithms.statistics.models.regression import OLSModel
  from nipy.testing import assert_almost_equal
  from ...None import design, hrf, utils
  from ..None import hrf
  from .FIACdesigns import N_ROWS, altdescr, descriptions, fmristat, time_vector
label_desc: Use image node_id 2 for calling a global function with example usage: protocol(recarr, design_type) and returns: f,  Tcontrasts,  Fcontrasts