def gen_design(stimtime_files, scan_duration, TR, style='FSL',
               temp_res=0.01,
               hrf_para={'response_delay': 6, 'undershoot_delay': 12,
                         'response_dispersion': 0.9,
                         'undershoot_dispersion': 0.9,
                         'undershoot_scale': 0.035}):
    """ Generate design matrix based on a list of names of stimulus
        timing files. The function will read each file, and generate
        a numpy array of size [time_points \\* condition], where
        time_points equals duration / TR, and condition is the size of
        stimtime_filenames. Each column is the hypothetical fMRI response
        based on the stimulus timing in the corresponding file
        of stimtime_files.
        This function uses generate_stimfunction and double_gamma_hrf
        of brainiak.utils.fmrisim.

    Parameters
    ----------

    stimtime_files: a string or a list of string.
        Each string is the name of the file storing
        the stimulus timing information of one task condition.
        The contents in the files will be interpretated
        based on the style parameter.
        Details are explained under the style parameter.

    scan_duration: float or a list (or a 1D numpy array) of numbers.
        Total duration of each fMRI scan, in unit of seconds.
        If there are multiple runs, the duration should be
        a list (or 1-d numpy array) of numbers.
        If it is a list, then each number in the list
        represents the duration of the corresponding scan
        in the stimtime_files.
        If only a number is provided, it is assumed that
        there is only one fMRI scan lasting for scan_duration.

    TR: float.
        The sampling period of fMRI, in unit of seconds.

    style: string, default: 'FSL'
        Acceptable inputs: 'FSL', 'AFNI'
        The formating style of the stimtime_files.
        'FSL' style has one line for each event of the same condition.
        Each line contains three numbers. The first number is the onset
        of the event relative to the onset of the first scan,
        in units of seconds.
        (Multiple scans should be treated as a concatenated long scan
        for the purpose of calculating onsets.
        However, the design matrix from one scan won't leak into the next).
        The second number is the duration of the event,
        in unit of seconds.
        The third number is the amplitude modulation (or weight)
        of the response.
        It is acceptable to not provide the weight,
        or not provide both duration and weight.
        In such cases, these parameters will default to 1.0.
        This code will accept timing files with only 1 or 2 columns for
        convenience but please note that the FSL package does not allow this

        'AFNI' style has one line for each scan (run).
        Each line has a few triplets in the format of
        stim_onsets*weight:duration
        (or simpler, see below), separated by spaces.
        For example, 3.2\\*2.0:1.5 means that one event starts at 3.2s,
        modulated by weight of 2.0 and lasts for 1.5s.
        If some run does not include a single event
        of a condition (stimulus type), then you can put \\*,
        or a negative number, or a very large number in that line.
        Either duration or weight can be neglected. In such
        cases, they will default to 1.0.
        For example, 3.0, 3.0\\*1.0, 3.0:1.0 and 3.0\\*1.0:1.0 all
        means an event starting at 3.0s, lasting for 1.0s, with
        amplitude modulation of 1.0.

    temp_res: float, default: 0.01
        Temporal resolution of fMRI, in second.

    hrf_para: dictionary
        The parameters of the double-Gamma hemodynamic response function.
        To set different parameters, supply a dictionary with
        the same set of keys as the default, and replace the corresponding
        values with the new values.

    Returns
    -------

    design: 2D numpy array
        design matrix. Each time row represents one TR
        (fMRI sampling time point) and each column represents
        one experiment condition, in the order in stimtime_files

    """
    if np.ndim(scan_duration) == 0:
        scan_duration = [scan_duration]
    scan_duration = np.array(scan_duration)
    assert np.all(scan_duration > TR), \
        'scan duration should be longer than a TR'
    if type(stimtime_files) is str:
        stimtime_files = [stimtime_files]
    assert TR > 0, 'TR should be positive'
    assert style == 'FSL' or style == 'AFNI', 'style can only be FSL or AFNI'
    n_C = len(stimtime_files)  # number of conditions
    n_S = np.size(scan_duration)  # number of scans
    if n_S > 1:
        design = [np.empty([int(np.round(duration / TR)), n_C])
                  for duration in scan_duration]
    else:
        design = [np.empty([int(np.round(scan_duration / TR)), n_C])]
    scan_onoff = np.insert(np.cumsum(scan_duration), 0, 0)
    if style == 'FSL':
        design_info = _read_stimtime_FSL(stimtime_files, n_C, n_S, scan_onoff)
    elif style == 'AFNI':
        design_info = _read_stimtime_AFNI(stimtime_files, n_C, n_S, scan_onoff)

    response_delay = hrf_para['response_delay']
    undershoot_delay = hrf_para['undershoot_delay']
    response_disp = hrf_para['response_dispersion']
    undershoot_disp = hrf_para['undershoot_dispersion']
    undershoot_scale = hrf_para['undershoot_scale']
    # generate design matrix
    for i_s in range(n_S):
        for i_c in range(n_C):
            if len(design_info[i_s][i_c]['onset']) > 0:
                stimfunction = generate_stimfunction(
                    onsets=design_info[i_s][i_c]['onset'],
                    event_durations=design_info[i_s][i_c]['duration'],
                    total_time=scan_duration[i_s],
                    weights=design_info[i_s][i_c]['weight'],
                    temporal_resolution=1.0/temp_res)
                hrf = _double_gamma_hrf(response_delay=response_delay,
                                        undershoot_delay=undershoot_delay,
                                        response_dispersion=response_disp,
                                        undershoot_dispersion=undershoot_disp,
                                        undershoot_scale=undershoot_scale,
                                        temporal_resolution=1.0/temp_res)
                design[i_s][:, i_c] = convolve_hrf(
                    stimfunction, TR, hrf_type=hrf, scale_function=False,
                    temporal_resolution=1.0 / temp_res).transpose() * temp_res
            else:
                design[i_s][:, i_c] = 0.0
            # We multiply the resulting design matrix with
            # the temporal resolution to normalize it.
            # We do not use the internal normalization
            # in double_gamma_hrf because it does not guarantee
            # normalizing with the same constant.
    return np.concatenate(design, axis=0)
Generate design matrix based on a list of names of stimulus
        timing files. The function will read each file, and generate
        a numpy array of size [time_points \\* condition], where
        time_points equals duration / TR, and condition is the size of
        stimtime_filenames. Each column is the hypothetical fMRI response
        based on the stimulus timing in the corresponding file
        of stimtime_files.
        This function uses generate_stimfunction and double_gamma_hrf
        of brainiak.utils.fmrisim.

    Parameters
    ----------

    stimtime_files: a string or a list of string.
        Each string is the name of the file storing
        the stimulus timing information of one task condition.
        The contents in the files will be interpretated
        based on the style parameter.
        Details are explained under the style parameter.

    scan_duration: float or a list (or a 1D numpy array) of numbers.
        Total duration of each fMRI scan, in unit of seconds.
        If there are multiple runs, the duration should be
        a list (or 1-d numpy array) of numbers.
        If it is a list, then each number in the list
        represents the duration of the corresponding scan
        in the stimtime_files.
        If only a number is provided, it is assumed that
        there is only one fMRI scan lasting for scan_duration.

    TR: float.
        The sampling period of fMRI, in unit of seconds.

    style: string, default: 'FSL'
        Acceptable inputs: 'FSL', 'AFNI'
        The formating style of the stimtime_files.
        'FSL' style has one line for each event of the same condition.
        Each line contains three numbers. The first number is the onset
        of the event relative to the onset of the first scan,
        in units of seconds.
        (Multiple scans should be treated as a concatenated long scan
        for the purpose of calculating onsets.
        However, the design matrix from one scan won't leak into the next).
        The second number is the duration of the event,
        in unit of seconds.
        The third number is the amplitude modulation (or weight)
        of the response.
        It is acceptable to not provide the weight,
        or not provide both duration and weight.
        In such cases, these parameters will default to 1.0.
        This code will accept timing files with only 1 or 2 columns for
        convenience but please note that the FSL package does not allow this

        'AFNI' style has one line for each scan (run).
        Each line has a few triplets in the format of
        stim_onsets*weight:duration
        (or simpler, see below), separated by spaces.
        For example, 3.2\\*2.0:1.5 means that one event starts at 3.2s,
        modulated by weight of 2.0 and lasts for 1.5s.
        If some run does not include a single event
        of a condition (stimulus type), then you can put \\*,
        or a negative number, or a very large number in that line.
        Either duration or weight can be neglected. In such
        cases, they will default to 1.0.
        For example, 3.0, 3.0\\*1.0, 3.0:1.0 and 3.0\\*1.0:1.0 all
        means an event starting at 3.0s, lasting for 1.0s, with
        amplitude modulation of 1.0.

    temp_res: float, default: 0.01
        Temporal resolution of fMRI, in second.

    hrf_para: dictionary
        The parameters of the double-Gamma hemodynamic response function.
        To set different parameters, supply a dictionary with
        the same set of keys as the default, and replace the corresponding
        values with the new values.

    Returns
    -------

    design: 2D numpy array
        design matrix. Each time row represents one TR
        (fMRI sampling time point) and each column represents
        one experiment condition, in the order in stimtime_files