Tools to calculate model using HRF and design


In [70]:
%matplotlib inline

import pandas as pd
from moss import glm
import numpy as np
from IPython.html.widgets import interact
import random
import scipy as sp

import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt

sns.set(context="notebook", style="ticks", font="Arial")
pd.set_option('display.precision', 3)

Calculate model based on event onsets


In [6]:
def calc_fromonsets(onset_list, tr=2, ax=None, 
                    hrf_model='GammaDifferenceHRF', temporal_deriv=False):
    
    """Plot the convolved response
    Parameters
    ----------
    onset_list : 1d array
        List of onset times relative to start (seconds)
    tr: int
        Length of TR (seconds)

    Prints
    -------
    Plot
    """
    if ax is None:
        ax = plt.gca()
    
    # Determine HRF
    hrf = getattr(glm, hrf_model)
    hrf = hrf(temporal_deriv, tr)

    # Make dataframe and convolve w/hrf
    onsets = []
    for onset in onset_list:
        onsets.append({'condition': 'trial', 'onset': onset})
    trial_onsets = pd.DataFrame(onsets)

    num_trs = int(np.ceil((trial_onsets.onset.max() + 10)/tr))
    model_ex = glm.DesignMatrix(trial_onsets, hrf, num_trs, tr=tr)

    # Plot
    ax.vlines(x=trial_onsets.onset/tr, ymin=-0.02, ymax=-0.017, 
           colors='darkgray', lw=4)
    ax.plot(model_ex.design_matrix.trial, color='dodgerblue', lw=2)
    ax.set_ylim(-0.02, model_ex.design_matrix.trial.max()+.001)
    ax.set_xlim(0, num_trs)
    ax.set_xlabel('TR')
    
    return ax

In [13]:
onset_list = [3, 12]
tr = 2
calc_fromonsets(onset_list)


Out[13]:
<matplotlib.axes._subplots.AxesSubplot at 0x10fa4c150>

In [15]:
trial_onsets


Out[15]:
condition onset
0 cond_1 3
1 cond_2 3
2 cond_1 12
3 cond_2 12

In [16]:
onsets = []
for onset, cond_num in zip(onset_list, range(1, len(onset_list)+1)):
    onsets.append({'condition': 'event_' + str(cond_num), 'onset': onset})
trial_onsets = pd.DataFrame(onsets)

num_trs = int(np.ceil((trial_onsets.onset.max() + 10)/tr))
model_ex = glm.DesignMatrix(trial_onsets, hrf, num_trs, tr=tr)

X = model_ex.design_matrix


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-16-8f33f2d3c9a4> in <module>()
      5 
      6 num_trs = int(np.ceil((trial_onsets.onset.max() + 10)/tr))
----> 7 model_ex = glm.DesignMatrix(trial_onsets, hrf, num_trs, tr=tr)
      8 
      9 X = model_ex.design_matrix

NameError: name 'hrf' is not defined

In [451]:
X.max().max()


Out[451]:
0.023849570082594689

In [443]:
sns.heatmap(X.T.dot(X), square=True)


Out[443]:
<matplotlib.axes._subplots.AxesSubplot at 0x11f9c6710>

In [444]:
sns.corrplot(X)


Out[444]:
<matplotlib.axes._subplots.AxesSubplot at 0x120bc4a90>

In [17]:
def calc_fromonsets_mult(onset_list, tr=2, num_cond=2,
                         hrf_model='GammaDifferenceHRF', 
                         temporal_deriv=False):
    
    """Plot the convolved response
    Parameters
    ----------
    onset_list : 1d array
        List of onset times relative to start (seconds)
    tr: int
        Length of TR (seconds)

    Prints
    -------
    Plot
    """
    f, axes = plt.subplots(nrows=2, figsize=(8,10))
    
    # Determine HRF
    hrf = getattr(glm, hrf_model)
    hrf = hrf(temporal_deriv, tr)
    
    # Make dataframe and convolve w/hrf
    onsets = []
    for onset, cond_num in zip(onset_list, range(1, len(onset_list)+1)):
        onsets.append({'condition': 'event_' + str(cond_num), 'onset': onset})
    trial_onsets = pd.DataFrame(onsets)

    num_trs = int(np.ceil((trial_onsets.onset.max() + 10)/tr))
    model_ex = glm.DesignMatrix(trial_onsets, hrf, num_trs, tr=tr)
    
    X = model_ex.design_matrix
    
#     sns.heatmap(np.cov(X.T), ax=axes[0], square=True, linewidths=0)
#     axes[0].set_title('Covariance Matrix')
    axes[0].set_title('Correlation Plot')
    sns.heatmap(np.corrcoef(X.T), ax=axes[0], square=True, linewidths=0)

    
    # Plot
    axes[1].vlines(x=trial_onsets.onset/tr, ymin=-0.02, ymax=-0.017, 
           colors='darkgray', lw=4)
    
    color_list = sns.color_palette("husl", len(X.columns))
    for ev, color in zip(X.columns, color_list):
        axes[1].plot(X[ev], color=color, lw=2)
        
    axes[1].set_ylim(-0.02, X.max().max()+.001)
    axes[1].set_xlim(0, num_trs)
    axes[1].set_xlabel('TR')
    axes[1].set_title('Convolved timeseries')

    f.tight_layout()

In [123]:
distributions = {'normal': np.random.normal(base_iti, 
                                            sd_iti, 
                                            onset_num-1),
                 'poisson': np.random.poisson(lam=base_iti,
                                              size=onset_num-1),
                 'exponential': np.random.exponential(scale=1,
                                                      size=41)+5}

iti_list = distributions['exponential']

onset_list = iti_list + 4
onset_list = np.cumsum(np.insert(onset_list, 0, 2))

tr=2
num_cond=2
hrf_model='GammaDifferenceHRF'
temporal_deriv=False

# Determine HRF
hrf = getattr(glm, hrf_model)
hrf = hrf(temporal_deriv, tr)

# Make dataframe and convolve w/hrf
onsets = []
for onset, cond_num in zip(onset_list, range(1, len(onset_list)+1)):
    onsets.append({'condition': 'event_' + str(cond_num), 'onset': onset})
trial_onsets = pd.DataFrame(onsets)

num_trs = int(np.ceil((trial_onsets.onset.max() + 10)/tr))
model_ex = glm.DesignMatrix(trial_onsets, hrf, num_trs, tr=tr)

X = model_ex.design_matrix

print X.columns
c = np.matrix(np.identity(42)); print c
1/sp.trace(c.T * np.mat(sp.linalg.pinv(np.dot(np.array(X).T, np.array(X))))* c)

sns.corrplot(X)

X.to_csv('/Users/steph-backup/Desktop/dmat.csv')


Index([u'event_1', u'event_10', u'event_11', u'event_12', u'event_13', u'event_14', u'event_15', u'event_16', u'event_17', u'event_18', u'event_19', u'event_2', u'event_20', u'event_21', u'event_22', u'event_23', u'event_24', u'event_25', u'event_26', u'event_27', u'event_28', u'event_29', u'event_3', u'event_30', u'event_31', u'event_32', u'event_33', u'event_34', u'event_35', u'event_36', u'event_37', u'event_38', u'event_39', u'event_4', u'event_40', u'event_41', u'event_42', u'event_5', u'event_6', u'event_7', u'event_8', u'event_9'], dtype='object')
[[ 1.  0.  0. ...,  0.  0.  0.]
 [ 0.  1.  0. ...,  0.  0.  0.]
 [ 0.  0.  1. ...,  0.  0.  0.]
 ..., 
 [ 0.  0.  0. ...,  1.  0.  0.]
 [ 0.  0.  0. ...,  0.  1.  0.]
 [ 0.  0.  0. ...,  0.  0.  1.]]

In [112]:



Out[112]:
array([[ 1123.90292564,    43.25373958,    26.60649883, ...,
          127.41593722,    92.24679286,    63.48396943],
       [   43.25373958,   889.60006584,   220.9078943 , ...,
          149.15300552,   175.80065566,   217.04097656],
       [   26.60649883,   220.9078943 ,   887.88405615, ...,
          134.18675988,   149.57178468,   170.71850146],
       ..., 
       [  127.41593722,   149.15300552,   134.18675988, ...,
          914.32324819,   247.45561171,   172.02925334],
       [   92.24679286,   175.80065566,   149.57178468, ...,
          247.45561171,   907.90006756,   234.2275804 ],
       [   63.48396943,   217.04097656,   170.71850146, ...,
          172.02925334,   234.2275804 ,   893.39788803]])

In [88]:
# onset_list = [2, 10, 20, 33, 40, 36, 43]
onset_list[5] = 20
print onset_list[1:10]
model_ex = calc_fromonsets_mult(onset_list[1:10])


[ 15.53131793  24.66133677  35.43892717  45.6526862   20.          63.93373834
  74.62810612  83.79677847  93.54839177]

Calculate models based on different ITIs


In [19]:
def calc_fromITI(iti_list, stim_dur, onset_num=2, tr=2, 
                 hrf_model='GammaDifferenceHRF', temporal_deriv=False):

    """Plot the convolved response
    Parameters
    ----------
    iti_list : 1d array
        List of ITIs (seconds)
    stim_dur : int
        Duration of stimulus
    onset_num : int
        The number of stimulus onsets
    tr: int
        Length of TR (seconds)

    Prints
    -------
    Plot
    """
    
    color_list = sns.color_palette("husl", len(iti_list))
    max_tr = int(np.ceil(((max(iti_list)+stim_dur)*onset_num + 2)/tr))
    
    for iti, color in zip(iti_list, color_list): 
        # Make stim list for this iti
        stim_list = [stim_dur] * (onset_num-1) # make list of dur for stim after 1st
        stim_list = np.array(stim_list) + iti # add the ITI to those values
        stim_list = np.insert(stim_list,0,2) # add in the first event
        onset_list = np.cumsum(stim_list) # find onset relative to last event

        # Determine HRF
        hrf = getattr(glm, hrf_model)
        hrf = hrf(temporal_deriv, tr)

        # Make dataframe and convolve w/hrf
        onsets = []
        for onset in onset_list:
            onsets.append({'condition': 'trial', 'onset': onset})
        trial_onsets = pd.DataFrame(onsets)

        model_ex = glm.DesignMatrix(trial_onsets, hrf, max_tr, tr=tr)

        # Plot
        plt.plot(model_ex.design_matrix.trial, color=color, lw=2)
        plt.ylim(-0.02, 0.05)
        plt.xlim(0, max_tr)
        plt.xlabel('TR')

    plt.legend(iti_list, title='ITI (s)')

In [20]:
iti_list = [6, 8, 10, 12]
stim_dur = 4
calc_fromITI(iti_list, stim_dur)



In [21]:
iti_list = [4, 10]
stim_dur = 4
calc_fromITI(iti_list, stim_dur, onset_num=5)



In [21]:

Calculate model from jittered ITIs


In [41]:
def calc_jitterITI(base_iti, dev_iti, stim_dur, 
                   onset_num=10, dist_iti='normal', 
                   hrf_model='GammaDifferenceHRF', 
                   plot=True):
    
    distributions = {'normal': np.random.normal(base_iti, 
                                                sd_iti, 
                                                onset_num-1),
                     'poisson': np.random.poisson(lam=base_iti,
                                                  size=onset_num-1),
                     'exponential': np.random.exponential(scale=dev_iti,
                                                          size=onset_num-1)+base_iti}

    iti_list = distributions[dist_iti]
    
    onset_list = iti_list + stim_dur
    onset_list = np.cumsum(np.insert(onset_list, 0, 2))
    
    if plot:
        print 'Total time (min): ' + str((onset_list[-1] + 22)/60) # 24 sec lead-in/out time


        f, axes = plt.subplots(nrows=2)
        sns.distplot(iti_list, ax=axes[0], rug=True, color='purple')
        axes[0].set_title('ITI distribution')
        calc_fromonsets(onset_list, ax=axes[1], hrf_model=hrf_model)
        axes[1].set_title('Convolved Timeseries')

        plt.tight_layout()
    else:
        total_seconds = onset_list[-1] + 22
        return total_seconds
    
    return iti_list

In [124]:
mean_iti = 6
sd_iti = 1
stim_dur = 4
onset_num = 43

iti_list = calc_jitterITI(mean_iti, sd_iti, stim_dur, onset_num, dist_iti='exponential')

iti_list


Total time (min): 8.2187895586
Out[124]:
array([  6.49221295,   6.20328611,   6.51859324,   7.99083164,
         9.37187985,   6.94065537,   6.5553255 ,   6.07851324,
         7.90552603,   6.04325642,   6.55160202,   8.88296853,
         8.89349571,   7.07843086,   8.25362589,   7.47624936,
         6.15383892,   6.27066917,   6.74582144,   7.78964066,
         6.93932639,   7.2247866 ,   6.78475466,   6.36581551,
         7.23401463,   7.30390242,   6.63741713,   6.58502626,
        11.90181607,   8.08080055,   7.12894261,   7.57374325,
         6.05618154,   6.17142151,   6.00833701,   7.89542954,
         6.0652307 ,   8.42235242,   6.06592273,   6.09881654,
         6.2800477 ,   8.10686484])

In [60]:
iti_list.tofile('/Users/steph-backup/Desktop/iti_list_jitter_exponential_mean5_sd1.txt', sep="\n")

In [24]:
mean_iti = 6
sd_iti = 2
stim_dur = 4
onset_num = 42

calc_jitterITI(mean_iti, sd_iti, stim_dur, onset_num)


Total time (min): 7.25265710345

In [42]:
interact(calc_jitterITI, base_iti=(2,10), dev_iti=(1.0,3.0,.5),
         stim_dur=(1,5), onset_num=(3,50), 
         hrf_model={'GammaDifferenceHRF': 'GammaDifferenceHRF',
                    'FIR': 'FIR'},
        dist_iti={'normal': 'normal', 
                  'poisson': 'poisson',
                  'exponential': 'exponential'});


Total time (min): 1.91825479095

In [29]:
def calc_jitterITI_mult(base_iti, dev_iti, stim_dur, 
                   onset_num=10, dist_iti='exponential', 
                   hrf_model='GammaDifferenceHRF', plot=True):
    
    distributions = {'normal': np.random.normal(base_iti, 
                                                dev_iti, 
                                                onset_num-1),
                     'poisson': np.random.poisson(lam=base_iti,
                                                  size=onset_num-1),
                     'exponential': np.random.exponential(scale=dev_iti,
                                                          size=onset_num-1)+base_iti}

    iti_list = distributions[dist_iti]
    
    onset_list = iti_list + stim_dur
    onset_list = np.cumsum(np.insert(onset_list, 0, 2))
    
    if plot:
        print 'Total time (min): ' + str((onset_list[-1] + 24)/60) # 24 sec lead-in/out time

        sns.distplot(iti_list, rug=True, color='purple')
        calc_fromonsets_mult(onset_list, hrf_model=hrf_model)

        plt.tight_layout()
        
    else:
        total_seconds = onset_list[-1] + 24
        return total_seconds

In [61]:
base_iti = 5
dev_iti = 1
stim_dur = 4
onset_num = 43
dist_iti = 'exponential'

calc_jitterITI_mult(base_iti, dev_iti, stim_dur, onset_num, dist_iti=dist_iti)


Total time (min): 7.38164515211

Distribution of Block Durations


In [62]:
blocklen_list = []
for i in range(1000):
    blocklen = calc_jitterITI_mult(base_iti, dev_iti, stim_dur, 
                                   onset_num, dist_iti=dist_iti, 
                                   plot=False)
    blocklen_list.append(blocklen)

In [63]:
sns.distplot(blocklen_list, color='seagreen')
plt.xlabel('Seconds/Block')

print 'Max TR Difference: %5.3f' %((max(blocklen_list) - min(blocklen_list))/2)


Max TR Difference: 23.192

In [ ]: