def export_epoch_file(stimfunction,
                      filename,
                      tr_duration,
                      temporal_resolution=100.0
                      ):
    """ Output an epoch file, necessary for some inputs into brainiak

    This takes in the time course of stimulus events and outputs the epoch
    file used in Brainiak. The epoch file is a way to structure the timing
    information in fMRI that allows you to flexibly input different stimulus
    sequences. This is a list with each entry a 3d matrix corresponding to a
    participant. The dimensions of the 3d matrix are condition by epoch by
    time. For the i-th condition, if its k-th epoch spans time points t_m to
    t_n-1, then [i, k, t_m:t_n] are 1 in the epoch file.

    Parameters
    ----------

    stimfunction : list of timepoint by condition arrays
        The stimulus function describing the time course of events. Each
        list entry is from a different participant, each row is a different
        timepoint (with the given temporal precision), each column is a
        different condition. export_epoch_file is looking for differences in
        the value of stimfunction to identify the start and end of an
        epoch. If epochs in stimfunction are coded with the same weight and
        there is no time between blocks then export_epoch_file won't be able to
        label them as different epochs

    filename : str
        The name of the epoch file to be output

    tr_duration : float
        How long is each TR in seconds

    temporal_resolution : float
        How many elements per second are you modeling with the
        stimfunction?

    """

    # Cycle through the participants, different entries in the list
    epoch_file = [0] * len(stimfunction)
    for ppt_counter in range(len(stimfunction)):

        # What is the time course for the participant (binarized)
        stimfunction_ppt = np.abs(stimfunction[ppt_counter]) > 0

        # Down sample the stim function
        stride = tr_duration * temporal_resolution
        stimfunction_downsampled = stimfunction_ppt[::int(stride), :]

        # Calculates the number of event onsets. This uses changes in value
        # to reflect different epochs. This might be false in some cases (the
        # weight is non-uniform over an epoch or there is no break between
        # identically weighted epochs).
        epochs = 0  # Preset
        conditions = stimfunction_ppt.shape[1]
        for condition_counter in range(conditions):

            weight_change = (np.diff(stimfunction_downsampled[:,
                                     condition_counter], 1, 0) != 0)

            # If the first or last events are 'on' then make these
            # represent a epoch change
            if stimfunction_downsampled[0, condition_counter] == 1:
                weight_change[0] = True
            if stimfunction_downsampled[-1, condition_counter] == 1:
                weight_change[-1] = True

            epochs += int(np.max(np.sum(weight_change, 0)) / 2)

        # Get other information
        trs = stimfunction_downsampled.shape[0]

        # Make a timing file for this participant
        epoch_file[ppt_counter] = np.zeros((conditions, epochs, trs))

        # Cycle through conditions
        epoch_counter = 0  # Reset and count across conditions
        tr_counter = 0
        while tr_counter < stimfunction_downsampled.shape[0]:

            for condition_counter in range(conditions):

                # Is it an event?
                if tr_counter < stimfunction_downsampled.shape[0] and \
                                stimfunction_downsampled[
                                    tr_counter, condition_counter] == 1:
                    # Add a one for this TR
                    epoch_file[ppt_counter][condition_counter,
                                            epoch_counter, tr_counter] = 1

                    # Find the next non event value
                    end_idx = np.where(stimfunction_downsampled[tr_counter:,
                                       condition_counter] == 0)[
                        0][0]
                    tr_idxs = list(range(tr_counter, tr_counter + end_idx))

                    # Add ones to all the trs within this event time frame
                    epoch_file[ppt_counter][condition_counter,
                                            epoch_counter, tr_idxs] = 1

                    # Start from this index
                    tr_counter += end_idx

                    # Increment
                    epoch_counter += 1

                # Increment the counter
                tr_counter += 1

        # Convert to boolean
        epoch_file[ppt_counter] = epoch_file[ppt_counter].astype('bool')

    # Save the file
    np.save(filename, epoch_file)
Output an epoch file, necessary for some inputs into brainiak

    This takes in the time course of stimulus events and outputs the epoch
    file used in Brainiak. The epoch file is a way to structure the timing
    information in fMRI that allows you to flexibly input different stimulus
    sequences. This is a list with each entry a 3d matrix corresponding to a
    participant. The dimensions of the 3d matrix are condition by epoch by
    time. For the i-th condition, if its k-th epoch spans time points t_m to
    t_n-1, then [i, k, t_m:t_n] are 1 in the epoch file.

    Parameters
    ----------

    stimfunction : list of timepoint by condition arrays
        The stimulus function describing the time course of events. Each
        list entry is from a different participant, each row is a different
        timepoint (with the given temporal precision), each column is a
        different condition. export_epoch_file is looking for differences in
        the value of stimfunction to identify the start and end of an
        epoch. If epochs in stimfunction are coded with the same weight and
        there is no time between blocks then export_epoch_file won't be able to
        label them as different epochs

    filename : str
        The name of the epoch file to be output

    tr_duration : float
        How long is each TR in seconds

    temporal_resolution : float
        How many elements per second are you modeling with the
        stimfunction?