
%matplotlib notebook

from pathlib import Path
from brainiak.utils import fmrisim
import nibabel
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as ndimage
import scipy.spatial.distance as sp_distance
import sklearn.manifold as manifold
import scipy.stats as stats
import sklearn.model_selection
import sklearn.svm


home = str(Path.home())
nii = nibabel.load(home + '/Corr_MVPA/Participant_01_rest_run01.nii')
volume = nii.get_data()


dim = volume.shape  # What is the size of the volume
dimsize = nii.header.get_zooms()  # Get voxel dimensions from the nifti header
tr = dimsize[3]
if tr > 100:  # If high then these values are likely in ms and so fix it
    tr /= 1000
print('Volume dimensions:', dim)
print('TR duration: %0.2fs' % tr)

Volume dimensions: (64, 64, 27, 294)
TR duration: 1.50s


mask, template = fmrisim.mask_brain(volume=volume, 
                                    mask_self=True,
                                    )


# Calculate the noise parameters from the data. Set it up to be matched.
noise_dict = {'voxel_size': [dimsize[0], dimsize[1], dimsize[2]], 'matched': 1}
noise_dict = fmrisim.calc_noise(volume=volume,
                                mask=mask,
                                template=template,
                                noise_dict=noise_dict,
                                )


print('Noise parameters of the data were estimated as follows:')
print('SNR: ' + str(noise_dict['snr']))
print('SFNR: ' + str(noise_dict['sfnr']))
print('FWHM: ' + str(noise_dict['fwhm']))

Noise parameters of the data were estimated as follows:
SNR: 23.175648200023815
SFNR: 70.7171164884859
FWHM: 5.661164924881941


# Calculate the noise given the parameters
noise = fmrisim.generate_noise(dimensions=dim[0:3],
                               tr_duration=int(tr),
                               stimfunction_tr=[0] * dim[3], 
                               mask=mask,
                               template=template,
                               noise_dict=noise_dict,
                               )

/Users/cellis/anaconda/envs/brainiak_test/lib/python3.6/site-packages/scipy/stats/stats.py:2279: RuntimeWarning: divide by zero encountered in true_divide
  np.expand_dims(sstd, axis=axis))
/Users/cellis/anaconda/envs/brainiak_test/lib/python3.6/site-packages/scipy/stats/stats.py:2279: RuntimeWarning: invalid value encountered in true_divide
  np.expand_dims(sstd, axis=axis))


# Plot a slice through the noise brain
plt.figure()
plt.title('Axial slice from template used for generating noise')
plt.imshow(noise[:, :, int(dim[2] / 2), 0], cmap=plt.cm.gray)
plt.axis('off')
txt="This template is generated by averaging the real functional data in time"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


# Plot spatial noise
low_spatial = fmrisim._generate_noise_spatial(dim[0:3],
                                              fwhm=4.0,
                                              )

high_spatial = fmrisim._generate_noise_spatial(dim[0:3],
                                               fwhm=1.0,
                                               )
plt.figure()
plt.subplot(1,2,1)
plt.title('FWHM = 4.0')
plt.imshow(low_spatial[:, :, 12])
plt.axis('off')

plt.subplot(1,2,2)
plt.title('FWHM = 1.0')
plt.imshow(high_spatial[:, :, 12])
plt.axis('off')

txt="Slices through the volume for different amounts of spatial noise. FWHM stands for full width half maximum, referencing the width of the Gaussian distribution used to simulate the data"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


# Create the different types of noise
total_time = 500
timepoints = list(range(0, total_time, int(tr)))

drift = fmrisim._generate_noise_temporal_drift(total_time,
                                               int(tr),
                                               )

mini_dim = np.array([2, 2, 2])
autoreg = fmrisim._generate_noise_temporal_autoregression(timepoints,
                                                          noise_dict,
                                                          mini_dim,
                                                          np.ones(mini_dim),
                                                          )
            
phys = fmrisim._generate_noise_temporal_phys(timepoints,
                                            )

stimfunc = np.zeros((int(total_time / tr), 1))
stimfunc[np.random.randint(0, int(total_time / tr), 50)] = 1
task = fmrisim._generate_noise_temporal_task(stimfunc,
                                            )


# Plot the different noise types
plt.figure()
plt.title('Noise types')

def clean_axis(ax):
    # Remove the borders and ticks of the plots (but different than plt.axis('off') since you can set a label)
    ax.spines["top"].set_visible(False)
    ax.spines["bottom"].set_visible(False)
    ax.spines["left"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.xticks([])
    plt.yticks([])

ax = plt.subplot(4, 1, 1)
plt.plot(drift)
clean_axis(ax)
plt.ylabel('Drift')

ax = plt.subplot(4, 1, 2)
plt.plot(autoreg[0, 0, 0, :])
clean_axis(ax)
plt.ylabel('AR')

ax = plt.subplot(4, 1, 3)
plt.plot(phys)
clean_axis(ax)
plt.ylabel('Physio')

ax = plt.subplot(4, 1, 4)
plt.plot(task)
clean_axis(ax)
plt.ylabel('Task')

txt="Time course of different noise properties that are usable in fmrisim."
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


# Dilate the mask so as to only take voxels far from the brain (performed in calc_noise)
mask_dilated = ndimage.morphology.binary_dilation(mask, iterations=10)

# Remove all non brain voxels
system_all = volume[mask_dilated == 0]  # Pull out all the non brain voxels in the first TR
system_baseline = volume - (template.reshape(dim[0], dim[1], dim[2], 1) * noise_dict['max_activity'])  # Subtract the baseline before masking
system_baseline = system_baseline[mask_dilated == 0]

# Plot the distribution of voxels
plt.figure(figsize=(10,8))
plt.subplot(1, 3, 1)
plt.hist(system_all[:,0].flatten(),100)
plt.title('Non-brain distribution')
plt.xlabel('Activity')
plt.ylabel('Frequency')

# Identify a subset of voxels to plot
idxs = list(range(system_all.shape[0]))
np.random.shuffle(idxs)

temporal = system_all[idxs[:100], :100]
plt.subplot(1, 3, 2)
plt.imshow(temporal)
plt.xticks([], [])
plt.yticks([], [])
plt.ylabel('voxel ID')
plt.xlabel('time')
plt.title('Voxel x time')

# Plot the difference
ax=plt.subplot(1, 3, 3)
plt.hist(system_baseline[:,0].flatten(),100)
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
plt.title('Demeaned non-brain distribution')
plt.xlabel('Activity difference')

txt="Histogram of non-brain voxel intensity. Left shows the raw intensity histogram, middle shows the voxel by time matrix (showing that there is variance in mean but not much in time) and then right shows the demeaned voxel intensity for the first time point."
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=10);


# Compute the noise parameters for the simulated noise
noise_dict_sim = {'voxel_size': [dimsize[0], dimsize[1], dimsize[2]], 'matched': 1}
noise_dict_sim = fmrisim.calc_noise(volume=noise,
                                    mask=mask,
                                    template=template,
                                    noise_dict=noise_dict_sim,
                                    )


print('Compare noise parameters for the real and simulated noise:')
print('SNR: %0.2f vs %0.2f' % (noise_dict['snr'], noise_dict_sim['snr']))
print('SFNR: %0.2f vs %0.2f' % (noise_dict['sfnr'], noise_dict_sim['sfnr']))
print('FWHM: %0.2f vs %0.2f' % (noise_dict['fwhm'], noise_dict_sim['fwhm']))
print('AR: %0.2f vs %0.2f' % (noise_dict['auto_reg_rho'][0], noise_dict_sim['auto_reg_rho'][0]))

Compare noise parameters for the real and simulated noise:
SNR: 23.18 vs 22.53
SFNR: 70.72 vs 69.19
FWHM: 5.66 vs 5.69
AR: 0.86 vs 0.84


# Create the region of activity where signal will appear
coordinates = np.array([[21, 21, 21]])  # Where in the brain is the signal
feature_size = 3  # How big, in voxels, is the size of the ROI
signal_volume = fmrisim.generate_signal(dimensions=dim[0:3],
                                        feature_type=['cube'],
                                        feature_coordinates=coordinates,
                                        feature_size=[feature_size],
                                        signal_magnitude=[1],
                                        )


plt.figure()
plt.imshow(signal_volume[:, :, 21], cmap=plt.cm.gray)
plt.imshow(mask[:, :, 21], cmap=plt.cm.gray, alpha=.5)
plt.axis('off')

txt="Mask of volume showing all voxels in gray. In white is the cube where signal is inserted"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


# Create a pattern for each voxel in our signal ROI
voxels = feature_size ** 3

# Pull the conical voxel activity from a uniform distribution
pattern_A = np.random.rand(voxels).reshape((voxels, 1))  
pattern_B = np.random.rand(voxels).reshape((voxels, 1))


# Plot pattern of activity for each condition
plt.figure()
plt.subplot(1,2,1)
plt.imshow(pattern_A)
plt.ylabel('Voxels')
plt.tick_params(which='both', left='off', labelleft='off', bottom='off', labelbottom='off')
plt.xticks([])
plt.xlabel('Condition A')

plt.subplot(1,2,2)
plt.imshow(pattern_B)
plt.tick_params(which='both', left='off', labelleft='off', bottom='off', labelbottom='off')
plt.xticks([])
plt.xlabel('Condition B')

txt="Randomly generated pattern that is used for simulating the response evoked by two different events"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


# Set up stimulus event time course parameters
event_duration = 2  # How long is each event
isi = 7  # What is the time between each event
burn_in = 1  # How long before the first event

total_time = int(dim[3] * tr) + burn_in  # How long is the total event time course
events = int((total_time - ((event_duration + isi) * 2))  / ((event_duration + isi) * 2)) * 2  # How many events are there?
onsets_all = np.linspace(burn_in, events * (event_duration + isi), events)  # Space the events out
np.random.shuffle(onsets_all)  # Shuffle their order
onsets_A = onsets_all[:int(events / 2)]  # Assign the first half of shuffled events to condition A
onsets_B = onsets_all[int(events / 2):]  # Assign the second half of shuffled events to condition B
temporal_res = 10.0 # How many timepoints per second of the stim function are to be generated?


# Create a time course of events 
stimfunc_A = fmrisim.generate_stimfunction(onsets=onsets_A,
                                           event_durations=[event_duration],
                                           total_time=total_time,
                                           temporal_resolution=temporal_res,
                                           )

stimfunc_B = fmrisim.generate_stimfunction(onsets=onsets_B,
                                           event_durations=[event_duration],
                                           total_time=total_time,
                                           temporal_resolution=temporal_res,
                                           )


fmrisim.export_epoch_file(stimfunction=[np.hstack((stimfunc_A, stimfunc_B))],
                          filename=home + '/epoch_file.npy',
                          tr_duration=tr,
                          temporal_resolution=temporal_res,
                          )

fmrisim.export_3_column(stimfunction=stimfunc_A,
                        filename=home + '/Condition_A.txt',
                        temporal_resolution=temporal_res,
                        )

fmrisim.export_3_column(stimfunction=stimfunc_B,
                        filename=home + '/Condition_B.txt',
                        temporal_resolution=temporal_res,
                        )


# Multiply each pattern by each voxel time course
weights_A = np.matlib.repmat(stimfunc_A, 1, voxels).transpose() * pattern_A
weights_B = np.matlib.repmat(stimfunc_B, 1, voxels).transpose() * pattern_B

# Sum these time courses together
stimfunc_weighted = weights_A + weights_B
stimfunc_weighted = stimfunc_weighted.transpose()


plt.figure(figsize=(6,7))
plt.plot(stimfunc_weighted[:, 0])
plt.title('Example voxel response time course')
plt.xlabel('Upsampled time course')
plt.ylabel('Response size')

txt="A signal voxel's simulated response over time, before convolution. The two heights of the bars refer to different event types"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


signal_func = fmrisim.convolve_hrf(stimfunction=stimfunc_weighted,
                                   tr_duration=tr,
                                   temporal_resolution=temporal_res,
                                   scale_function=1,
                                   )


# Prepare the data to be plotted
response = signal_func[0:100,0] * 2
downsample_A = stimfunc_A[0:int(100*temporal_res * tr):int(temporal_res * tr), 0]
downsample_B = stimfunc_B[0:int(100*temporal_res * tr):int(temporal_res * tr), 0]

# Display signal
plt.figure(figsize=(5,6))
plt.title('Example event time course and voxel response')
Event_A = plt.plot(downsample_A, 'r', label='Event_A')
Event_B = plt.plot(downsample_B, 'g', label='Event_B')
Response = plt.plot(response, 'b', label='Response')
plt.legend(loc=1)
plt.yticks([],'')
plt.xlabel('nth TR')

txt="A signal voxel's response convolved with a double-gamma HRF. Event types are also shown."
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=8);


# Specify the parameters for signal
signal_method = 'CNR_Amp/Noise-SD'
signal_magnitude = [0.5]

# Where in the brain are there stimulus evoked voxels
signal_idxs = np.where(signal_volume == 1)

# Pull out the voxels corresponding to the noise volume
noise_func = noise[signal_idxs[0], signal_idxs[1], signal_idxs[2], :].T


# Compute the signal appropriate scaled
signal_func_scaled = fmrisim.compute_signal_change(signal_func,
                                                  noise_func,
                                                  noise_dict,
                                                  magnitude=signal_magnitude,
                                                  method=signal_method,
                                                  )


signal = fmrisim.apply_signal(signal_func_scaled,
                              signal_volume,
                              )


brain = signal + noise


hrf_lag = 4  # Assumed time from stimulus onset to HRF peak

# Get the lower and upper bounds of the ROI
lb = (coordinates - ((feature_size - 1) / 2)).astype('int')[0]
ub = (coordinates + ((feature_size - 1) / 2) + 1).astype('int')[0]

# Pull out voxels in the ROI for the specified timepoints
trials_A = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], ((onsets_A + hrf_lag) / tr).astype('int')]
trials_B = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], ((onsets_B + hrf_lag) / tr).astype('int')]

# Reshape data for easy handling
trials_A = trials_A.reshape((voxels, trials_A.shape[3]))
trials_B = trials_B.reshape((voxels, trials_B.shape[3]))


# Plot the pattern of activity for our signal voxels at each timepoint
plt.figure()
plt.subplot(1,2,1)
plt.imshow(trials_A)
plt.ylabel('Voxels')
plt.xlabel('Trials')
plt.subplot(1,2,2)
plt.imshow(trials_B)
plt.xlabel('Trials')

txt="Time course of activity for each signal voxel when combining signal and noise"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


# Calculate the distance matrix between trial types
distance_matrix = sp_distance.squareform(sp_distance.pdist(np.vstack([trials_A.transpose(), trials_B.transpose()])))

mds = manifold.MDS(n_components=2, dissimilarity='precomputed')  # Fit the mds object
coords = mds.fit(distance_matrix).embedding_  # Find the mds coordinates

# Plot the data
plt.figure()
plt.scatter(coords[:, 0], coords[:, 1], c=['red'] * trials_A.shape[1] + ['green'] * trials_B.shape[1])
plt.axis('off')
plt.title('Low Dimensional Representation of conditions A and B')

txt="Low dimensional representation of the two conditions. Different colors refer to different conditions"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment='center', fontsize=12);


mean_difference = (np.mean(trials_A,0) - np.mean(trials_B,0))
ttest = stats.ttest_1samp(mean_difference, 0)

print('Mean difference between condition A and B: %0.2f\np value: %0.3f' % (mean_difference.mean(), ttest.pvalue))

Mean difference between condition A and B: 3.14
p value: 0.229


# Get the inputs to the SVM
input_mat = np.vstack([trials_A.transpose(), trials_B.transpose()])
input_labels = trials_A.shape[1] * [1] + trials_B.shape[1] * [0]

# Set up the classifier
X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(
    input_mat, input_labels, test_size=0.2, random_state=0)

clf = sklearn.svm.SVC(kernel='linear', C=1, gamma='auto').fit(X_train, y_train)

score = clf.score(X_test, y_test)
print('Classification accuracy between condition A and B: %0.3f' % score)

Classification accuracy between condition A and B: 0.700

fmrisim demo script¶

Overview¶

Annotated bibliography¶

Table of Contents¶

1. Set parameters. ¶

2. Generate noise¶

3. Generate signal¶

4. Analyse data¶

Summary ¶

References¶

fmrisim demo script¶

Overview¶

Annotated bibliography¶

Table of Contents¶

1. Set parameters. ¶

2. Generate noise¶

3. Generate signal¶

4. Analyse data¶

Summary¶

References¶

Summary ¶