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)
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]:
In [15]:
trial_onsets
Out[15]:
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
In [451]:
X.max().max()
Out[451]:
In [443]:
sns.heatmap(X.T.dot(X), square=True)
Out[443]:
In [444]:
sns.corrplot(X)
Out[444]:
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')
In [112]:
Out[112]:
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])
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]:
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
Out[124]:
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)
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'});
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)
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)
In [ ]: