diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 7c6d0819a..41c5d73dd 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -66,25 +66,27 @@ use an MNI grey matter atlas but any image can be supplied to create an estimate. -plot_brain -Display the brain, timepoint by timepoint, with above threshold voxels -highlighted against the outline of the brain. - +compute_signal_change +Convert the signal function into useful metric units according to metrics +used by others (Welvaert & Rosseel, 2013) Authors: - Cameron Ellis (Princeton) 2016-2017 + Cameron Ellis (Princeton & Yale) 2016-2018 Chris Baldassano (Princeton) 2016-2017 + Mingbo Cai (Princeton) 2017 """ import logging from itertools import product -import nitime.algorithms.autoregressive as ar +from statsmodels.tsa.arima_model import ARMA import math import numpy as np +from numpy.linalg import LinAlgError from pkg_resources import resource_stream from scipy import stats from scipy import signal import scipy.ndimage as ndimage +import copy __all__ = [ "generate_signal", @@ -96,7 +98,7 @@ "calc_noise", "generate_noise", "mask_brain", - "plot_brain", + "compute_signal_change", ] logger = logging.getLogger(__name__) @@ -487,7 +489,7 @@ def generate_stimfunction(onsets, 'Consider increasing' \ ' the temporal ' \ 'resolution.' - logging.warning(warning) + logger.warning(warning) onsets.append(float(onset)) event_durations.append(float(duration)) @@ -608,7 +610,9 @@ def export_epoch_file(stimfunction, 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 + 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 ---------- @@ -624,7 +628,7 @@ def export_epoch_file(stimfunction, label them as different epochs filename : str - The name of the three column text file to be output + The name of the epoch file to be output tr_duration : float How long is each TR in seconds @@ -637,58 +641,65 @@ def export_epoch_file(stimfunction, # Cycle through the participants, different entries in the list epoch_file = [0] * len(stimfunction) - for participant_counter in range(len(stimfunction)): + for ppt_counter in range(len(stimfunction)): # What is the time course for the participant (binarized) - stimfunction_ppt = np.abs(stimfunction[participant_counter]) > 0 + stimfunction_ppt = np.abs(stimfunction[ppt_counter]) > 0 - # Cycle through conditions + # 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): - # Down sample the stim function - stride = tr_duration * temporal_resolution - stimfunction_temp = stimfunction_ppt[:, condition_counter] - stimfunction_temp = stimfunction_temp[::int(stride)] + weight_change = (np.diff(stimfunction_downsampled[:, + condition_counter], 1, 0) != 0) - if condition_counter == 0: - # Calculates the number of event onsets (max of all - # conditions). This uses changes in value to reflect - # different epochs. This might be false in some cases (the - # weight is supposed to unfold over an epoch or there is no - # break between identically weighted epochs). In such cases - # this will not work - weight_change = (np.diff(stimfunction_temp, 1, 0) != 0) - epochs = int(np.max(np.sum(weight_change, 0)) / 2) + # 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 - # Get other information - trs = stimfunction_temp.shape[0] + epochs += int(np.max(np.sum(weight_change, 0)) / 2) - # Make a timing file for this participant - epoch_file[participant_counter] = np.zeros((conditions, - epochs, trs)) + # Get other information + trs = stimfunction_downsampled.shape[0] - epoch_counter = 0 - tr_counter = 0 - while tr_counter < stimfunction_temp.shape[0]: + # Make a timing file for this participant + epoch_file[ppt_counter] = np.zeros((conditions, epochs, trs)) - # Is it an event? - if stimfunction_temp[tr_counter] == 1: + # 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[participant_counter][condition_counter, - epoch_counter, - tr_counter] = 1 + epoch_file[ppt_counter][condition_counter, + epoch_counter, tr_counter] = 1 # Find the next non event value - end_idx = np.where(stimfunction_temp[tr_counter:] == 0)[ + 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[participant_counter][condition_counter, - epoch_counter, - tr_idxs] = 1 + epoch_file[ppt_counter][condition_counter, + epoch_counter, tr_idxs] = 1 # Start from this index tr_counter += end_idx @@ -699,6 +710,9 @@ def export_epoch_file(stimfunction, # 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) @@ -808,7 +822,7 @@ def convolve_hrf(stimfunction, Parameters ---------- - stimfunction : timepoint by timecourse array + stimfunction : timepoint by feature array What is the time course of events to be modelled in this experiment. This can specify one or more timecourses of events. The events can be weighted or binary @@ -836,6 +850,11 @@ def convolve_hrf(stimfunction, columns in this array. """ + + # Check if it is timepoint by feature + if stimfunction.shape[0] < stimfunction.shape[1]: + logger.warning('Stimfunction may be the wrong shape') + # How will stimfunction be resized stride = int(temporal_resolution * tr_duration) duration = int(stimfunction.shape[0] / stride) @@ -1106,16 +1125,17 @@ def _calc_sfnr(volume, def _calc_snr(volume, mask, - tr=None, + dilation=5, + reference_tr=None, ): """ Calculate the the SNR of a volume Calculates the Signal to Noise Ratio, the mean of brain voxels divided by the standard deviation across non-brain voxels. Specify a TR value to calculate the mean and standard deviation for that TR. To - calculate the standard deviation this subtracts any baseline structure - in the non-brain voxels, hence getting at deviations due to the system - noise and not something like high baseline values in non-brain parts of - the body. + calculate the standard deviation of non-brain voxels we can subtract + any baseline structure away first, hence getting at deviations due to the + system noise and not something like high baseline values in non-brain + parts of the body. Parameters ---------- @@ -1126,8 +1146,15 @@ def _calc_snr(volume, mask : 3d array, binary A binary mask the same size as the volume - tr : int - Integer specifying TR to calculate the SNR for + dilation : int + How many binary dilations do you want to perform on the mask to + determine the non-brain voxels. If you increase this the SNR + increases and the non-brain voxels (after baseline subtraction) more + closely resemble a gaussian + + reference_tr : int or list + Specifies the TR to calculate the SNR for. If multiple are supplied + then it will use the average of them. Returns ------- @@ -1137,86 +1164,128 @@ def _calc_snr(volume, """ - # If no TR is specified then take the middle one - if tr is None: - tr = int(np.ceil(volume.shape[3] / 2)) + # If no TR is specified then take all of them + if reference_tr is None: + reference_tr = list(range(volume.shape[3])) - # Make a matrix of brain and non_brain voxels by time - brain_voxels = volume[mask > 0] - nonbrain_voxels = volume[:, :, :, tr][mask == 0] + # Dilate the mask in order to ensure that non-brain voxels are far from + # the brain + if dilation > 0: + mask_dilated = ndimage.morphology.binary_dilation(mask, + iterations=dilation) + else: + mask_dilated = mask + + # Make a matrix of brain and non_brain voxels, selecting the timepoint/s + brain_voxels = volume[mask > 0][:, reference_tr] + nonbrain_voxels = (volume[:, :, :, reference_tr]).astype('float64') - # Find the mean of the non_brain voxels (deals with structure that may - # exist outside of the mask) - nonbrain_voxels_mean = np.mean(volume[mask == 0], 1) + # If you have multiple TRs + if len(brain_voxels.shape) > 1: + brain_voxels = np.mean(brain_voxels, 1) + nonbrain_voxels = np.mean(nonbrain_voxels, 3) + + nonbrain_voxels = nonbrain_voxels[mask_dilated == 0] # Take the means of each voxel over time mean_voxels = np.nanmean(brain_voxels) - std_voxels = np.nanstd(nonbrain_voxels - nonbrain_voxels_mean) + + # Find the standard deviation of the voxels + std_voxels = np.nanstd(nonbrain_voxels) # Return the snr return mean_voxels / std_voxels -def _calc_temporal_noise(volume, - mask, - auto_reg_order=1, - ): - """ Calculate the the temporal noise of a volume - This calculates the variability of the volume over time and the - proportion of variance over time that is due to autoregression and how - much is due to scanner drift. +def _calc_ARMA_noise(volume, + mask, + auto_reg_order=1, + ma_order=1, + sample_num=100, + ): + """ Calculate the the ARMA noise of a volume + This calculates the autoregressive and moving average noise of the volume + over time by sampling brain voxels and averaging them. Parameters ---------- - volume : 4d array, float + volume : 4d array or 1d array, float Take a volume time series to extract the middle slice from the - middle TR + middle TR. Can also accept a one dimensional time course (mask input + is then ignored). mask : 3d array, binary A binary mask the same size as the volume auto_reg_order : int - What order of the autoregression do you want to pull out + What order of the autoregression do you want to estimate + sample_num : int + How many voxels would you like to sample to calculate the AR values. + The AR distribution of real data is approximately exponential maxing + at 1. From analyses across a number of participants, to get less + than 3% standard deviation of error from the true mean it is + necessary to sample at least 100 voxels. Returns ------- + auto_reg_rho : list of floats + Rho of a specific order for the autoregression noise in the data - sfnr : float
 - The SFNR of the volume (mean brain activity divided by
 temporal - variability in the averaged non brain voxels)

 + na_rho : list of floats + Moving average of a specific order for the data - auto_reg_sigma : float - A sigma of the autoregression in the data + """ - drift_sigma : float - Sigma of the drift in the data + # Pull out the non masked voxels + if len(volume.shape) > 1: + brain_timecourse = volume[mask > 0] + else: + # If a 1 dimensional input is supplied then reshape it to make the + # timecourse + brain_timecourse = volume.reshape(1, len(volume)) - """ + # Identify some brain voxels to assess + voxel_idxs = list(range(brain_timecourse.shape[0])) + np.random.shuffle(voxel_idxs) + + # If there are more samples than voxels, take all of the voxels + if len(voxel_idxs) < sample_num: + sample_num = len(voxel_idxs) - # Calculate sfnr and convert from memmap - sfnr = _calc_sfnr(volume, - mask, - ) + auto_reg_rho_all = np.zeros((sample_num, auto_reg_order)) + ma_all = np.zeros((sample_num, ma_order)) + for voxel_counter in range(sample_num): - # Calculate the time course of voxels within the brain - timecourse = np.mean(volume[mask > 0], 0) - demeaned_timecourse = timecourse-timecourse.mean() + # Get the timecourse and demean it + timecourse = brain_timecourse[voxel_idxs[voxel_counter], :] + demeaned_timecourse = timecourse - timecourse.mean() - # Pull out the AR values (depends on order) - auto_reg_sigma = ar.AR_est_YW(demeaned_timecourse, auto_reg_order) - auto_reg_sigma = np.sqrt(auto_reg_sigma[1]) + # Pull out the ARMA values (depends on order) + try: + model = ARMA(demeaned_timecourse, [auto_reg_order, ma_order]) + model_fit = model.fit(disp=False) + params = model_fit.params + except (ValueError, LinAlgError): + params = np.ones(auto_reg_order + ma_order + 1) * np.nan - # What is the size of the change in the time course - drift_sigma = timecourse.std().tolist() + # Add to the list + auto_reg_rho_all[voxel_counter, :] = params[1:auto_reg_order + 1] + ma_all[voxel_counter, :] = params[auto_reg_order + 1:] - return sfnr, auto_reg_sigma, drift_sigma + # Average all of the values and then convert them to a list + auto_reg_rho = np.nanmean(auto_reg_rho_all, 0).tolist() + ma_rho = np.nanmean(ma_all, 0).tolist() + + # Return the coefficients + return auto_reg_rho, ma_rho def calc_noise(volume, - mask=None, + mask, + template, noise_dict=None, ): """ Calculates the noise properties of the volume supplied. @@ -1234,7 +1303,13 @@ def calc_noise(volume, mask : 3d numpy array, binary A binary mask of the brain, the same size as the volume + template : 3d array, float + A continuous (0 -> 1) volume describing the likelihood a voxel is in + the brain. This can be used to contrast the brain and non brain. + noise_dict : dict + The initialized dictionary of the calculated noise parameters of the + provided dataset (usually it is only the voxel size) Returns ------- @@ -1244,9 +1319,13 @@ def calc_noise(volume, """ + # Check the inputs + if template.max() > 1.1: + raise ValueError('Template out of range') + # Create the mask if not supplied and set the mask size if mask is None: - mask = np.ones(volume.shape[:-1]) + raise ValueError('Mask not supplied') # Update noise dict if it is not yet created if noise_dict is None: @@ -1258,12 +1337,24 @@ def calc_noise(volume, # convert between the mask and the mean of the brain volume) noise_dict['max_activity'] = np.nanmax(np.mean(volume, 3)) - # Since you are deriving the 'true' values then you want your noise to - # be set to that level - # Calculate the temporal variability of the volume - sfnr, auto_reg, drift = _calc_temporal_noise(volume, mask) - noise_dict['sfnr'] = sfnr + noise_dict['auto_reg_rho'], noise_dict['ma_rho'] = _calc_ARMA_noise( + volume, mask) + + # Set it such that all of the temporal variability will be accounted for + # by the AR component + noise_dict['auto_reg_sigma'] = 1 + + # Preset these values to be zero, as in you are not attempting to + # simulate them + noise_dict['physiological_sigma'] = 0 + noise_dict['task_sigma'] = 0 + noise_dict['drift_sigma'] = 0 + + # Calculate the sfnr + noise_dict['sfnr'] = _calc_sfnr(volume, + mask, + ) # Calculate the fwhm on a subset of volumes if volume.shape[3] > 100: @@ -1286,13 +1377,6 @@ def calc_noise(volume, mask, ) - # Total temporal noise, since these values only make sense relatively - total_temporal_noise = auto_reg + drift - - # What proportion of noise is accounted for by these variables? - noise_dict['auto_reg_sigma'] = auto_reg / total_temporal_noise - noise_dict['drift_sigma'] = drift / total_temporal_noise - # Return the noise dictionary return noise_dict @@ -1300,15 +1384,23 @@ def calc_noise(volume, def _generate_noise_system(dimensions_tr, spatial_sd, temporal_sd, - spatial_noise_type='exponential', - temporal_noise_type='exponential', + spatial_noise_type='gaussian', + temporal_noise_type='gaussian', ): """Generate the scanner noise - Generate system noise, either rician or exponential, for the scanner. - Low SNR scans tend to have rician noise whereas high SNR scans (>30) are - better modelled by exponential noise. Generates a distribution with a SD - of 1. + Generate system noise, either rician, gaussian or exponential, for the + scanner. Generates a distribution with a SD of 1. If you look at the + distribution of non-brain voxel intensity in modern scans you will see + it is rician. However, depending on how you have calculated the SNR and + whether the template is being used you will want to use this function + differently: the voxels outside the brain tend to be stable over time and + usually reflect structure in the MR signal (e.g. the + baseline MR of the head coil or skull). Hence the template captures this + rician noise structure. If you are adding the machine noise to the + template, as is done in generate_noise, then you are likely doubling up + on the addition of machine noise. In such cases, machine noise seems to + be better modelled by gaussian noise on top of this rician structure. Parameters ---------- @@ -1317,9 +1409,20 @@ def _generate_noise_system(dimensions_tr, What are the dimensions of the volume you wish to insert noise into. This can be a volume of any size + spatial_sd : float + What is the standard deviation in space of the noise volume to be + generated + + temporal_sd : float + What is the standard deviation in time of the noise volume to be + generated + noise_type : str - String specifying the noise type. Rician is appropriate when the SNR is - low but is insufficiently skewed to appropriately model high SNR data. + String specifying the noise type. If you aren't specifying the noise + template then Rician is the appropriate model of noise. However, + if you are subtracting the template, as is default, then you should + use gaussian. (If the dilation parameter of _calc_snr is <10 then + gaussian is only an approximation) Returns ---------- @@ -1328,9 +1431,9 @@ def _generate_noise_system(dimensions_tr, Create a volume with system noise """ - def generate_noise_volume(dimensions, - noise_type, - ): + def noise_volume(dimensions, + noise_type, + ): if noise_type == 'rician': # Generate the Rician noise (has an SD of 1) @@ -1351,38 +1454,24 @@ def generate_noise_volume(dimensions, 1]) # Generate noise - spatial_noise = generate_noise_volume(dimensions, spatial_noise_type) - temporal_noise = generate_noise_volume(dimensions_tr, temporal_noise_type) - - # Since you are combining spatial and temporal noise, you need to - # subtract the variance of the two to get the spatial sd - if spatial_sd > temporal_sd: - spatial_sd = np.sqrt(spatial_sd ** 2 - temporal_sd ** 2) - else: - # If this is below zero then all the noise will be temporal - spatial_sd = 0 - - # # Mean centre, while preserving the SD - # spatial_noise = spatial_noise - spatial_noise.mean() + spatial_noise = noise_volume(dimensions, spatial_noise_type) + temporal_noise = noise_volume(dimensions_tr, temporal_noise_type) # Make the system noise have a specific spatial variability spatial_noise *= spatial_sd - # # Mean centre, while preserving the SD - # temporal_noise = temporal_noise - temporal_noise.mean() - # Set the size of the noise temporal_noise *= temporal_sd # The mean in time of system noise needs to be zero, so subtract the - # means of the temporal noise in time and spatial noise + # means of the temporal noise in time temporal_noise_mean = np.mean(temporal_noise, 3).reshape(dimensions[0], dimensions[1], dimensions[2], 1) - temporal_noise = temporal_noise - (temporal_noise_mean - spatial_noise) + temporal_noise = temporal_noise - temporal_noise_mean - # Save the size of the noise + # Save the combination system_noise = spatial_noise + temporal_noise return system_noise @@ -1426,7 +1515,7 @@ def _generate_noise_temporal_task(stimfunction_tr, noise_task = stimfunction_tr + noise # Normalize - noise_task = stats.zscore(noise_task) + noise_task = stats.zscore(noise_task).flatten() return noise_task @@ -1515,13 +1604,14 @@ def _generate_noise_temporal_drift(trs, def _generate_noise_temporal_autoregression(timepoints, - auto_reg_order=1, - auto_reg_rho=[0.5], + noise_dict, + dimensions, + mask, ): """Generate the autoregression noise Make a slowly drifting timecourse with the given autoregression - parameters. The output should have an autoregression coefficient of 1 + parameters. This can take in both AR and MA components Parameters ---------- @@ -1529,15 +1619,21 @@ def _generate_noise_temporal_autoregression(timepoints, timepoints : 1 Dimensional array What time points are sampled by a TR - auto_reg_order : float - How many timepoints ought to be taken into consideration for the - autoregression function + noise_dict : dict + A dictionary specifying the types of noise in this experiment. The + noise types interact in important ways. First, all noise types + ending with sigma (e.g. motion sigma) are mixed together in + _generate_temporal_noise. The sigma values describe the proportion of + mixing of these elements. However critically, SFNR is the + parameter that describes how much noise these components contribute + to the brain. If you set the noise dict to matched then it will fit + the parameters to match the participant as best as possible. - auto_reg_rho : float - What is the scaling factor on the predictiveness of the previous - time point. This value is below 1 to avoid brownian motion (and - growing variance). Values near or greater than one may produce drift or - other unwanted trends. + dimensions : 3 length array, int + What is the shape of the volume to be generated + + mask : 3 dimensional array, binary + The masked brain, thresholded to distinguish brain and non-brain Returns ---------- @@ -1546,34 +1642,72 @@ def _generate_noise_temporal_autoregression(timepoints, """ - if len(auto_reg_rho) == 1: - auto_reg_rho = auto_reg_rho * auto_reg_order # Duplicate this so that - # there is one - # for each value + # Pull out the relevant noise parameters + auto_reg_rho = noise_dict['auto_reg_rho'] + ma_rho = noise_dict['ma_rho'] + + # Specify the order based on the number of rho supplied + auto_reg_order = len(auto_reg_rho) + ma_order = len(ma_rho) + + # This code assumes that the AR order is higher than the MA order + if ma_order > auto_reg_order: + msg = 'MA order (%d) is greater than AR order (%d). Cannot run.' % ( + ma_order, auto_reg_order) + raise ValueError(msg) # Generate a random variable at each time point that is a decayed value # of the previous time points - noise_autoregression = [] + noise_autoregression = np.zeros((dimensions[0], dimensions[1], + dimensions[2], len(timepoints))) + err_vols = np.zeros((dimensions[0], dimensions[1], dimensions[2], + len(timepoints))) for tr_counter in range(len(timepoints)): + # Create a brain shaped volume with appropriate smoothing properties + noise = _generate_noise_spatial(dimensions=dimensions, + mask=mask, + fwhm=noise_dict['fwhm'], + ) + + # Store all of the noise volumes + err_vols[:, :, :, tr_counter] = noise + if tr_counter == 0: - noise_autoregression.append(np.random.normal(0, 1)) + noise_autoregression[:, :, :, tr_counter] = noise else: - temp = [] + # Preset the volume to collect the AR estimated process + AR_vol = np.zeros((dimensions[0], dimensions[1], dimensions[2])) + + # Iterate through both the AR and MA values for pCounter in list(range(1, auto_reg_order + 1)): + past_TR = int(tr_counter - pCounter) + if tr_counter - pCounter >= 0: - past_trs = noise_autoregression[int(tr_counter - pCounter)] - past_reg = auto_reg_rho[pCounter - 1] - temp.append(past_trs * past_reg) - random = np.random.normal(0, 1) - noise_autoregression.append(np.sum(temp) + random) + # Pull out a previous TR + past_vols = noise_autoregression[:, :, :, past_TR] - # N.B. You don't want to normalize. Although that may make the sigma of - # this timecourse 1, it will change the autoregression coefficient to be - # much lower. + # Add the discounted previous volume + AR_vol += past_vols * auto_reg_rho[pCounter - 1] + + # If the MA order has at least this many coefficients + # then consider the error terms + if ma_order >= pCounter: + + # Pull out a previous TR + past_noise = err_vols[:, :, :, past_TR] + + # Add the discounted previous noise + AR_vol += past_noise * ma_rho[pCounter - 1] + + noise_autoregression[:, :, :, tr_counter] = AR_vol + noise + + # Z score the data so that all of the standard deviations of the voxels + # are one (but the ARMA coefs are unchanged) + noise_autoregression = stats.zscore(noise_autoregression, 3) return noise_autoregression @@ -1594,10 +1728,10 @@ def _generate_noise_temporal_phys(timepoints, What time points, in seconds, are sampled by a TR resp_freq : float - What is the frequency of respiration + What is the frequency of respiration (in Hz) heart_freq : float - What is the frequency of heart beat + What is the frequency of heart beat (in Hz) Returns ---------- @@ -1606,18 +1740,20 @@ def _generate_noise_temporal_phys(timepoints, """ - noise_phys = [] # Preset resp_phase = (np.random.rand(1) * 2 * np.pi)[0] heart_phase = (np.random.rand(1) * 2 * np.pi)[0] - for tr_counter in timepoints: - # Calculate the radians for each variable at this - # given TR - resp_radians = resp_freq * tr_counter * 2 * np.pi + resp_phase - heart_radians = heart_freq * tr_counter * 2 * np.pi + heart_phase + # Find the rate for each timepoint + resp_rate = (resp_freq * 2 * np.pi) + heart_rate = (heart_freq * 2 * np.pi) + + # Calculate the radians for each variable at this + # given TR + resp_radians = np.multiply(timepoints, resp_rate) + resp_phase + heart_radians = np.multiply(timepoints, heart_rate) + heart_phase - # Combine the two types of noise and append - noise_phys.append(np.cos(resp_radians) + np.sin(heart_radians)) + # Combine the two types of noise and append + noise_phys = np.cos(resp_radians) + np.sin(heart_radians) # Normalize noise_phys = stats.zscore(noise_phys) @@ -1626,7 +1762,6 @@ def _generate_noise_temporal_phys(timepoints, def _generate_noise_spatial(dimensions, - template=None, mask=None, fwhm=4.0, ): @@ -1678,10 +1813,12 @@ def _generate_noise_spatial(dimensions, """ + # Check the input is correct if len(dimensions) == 4: - return + logger.warning('4 dimensions have been supplied, only using 3') + dimensions = dimensions[0:3] - def logfunc(x, a, b, c): + def _logfunc(x, a, b, c): """Solve for y given x for log function. Parameters @@ -1707,44 +1844,88 @@ def logfunc(x, a, b, c): """ return (np.log(x + a) / np.log(b)) + c + def _fftIndgen(n): + """# Specify the fft coefficents + + Parameters + ---------- + + n : int + Dim size to estimate over + + Returns + ---------- + + array of ints + fft indexes + """ + + # Pull out the ascending and descending indexes + ascending = np.linspace(0, int(n / 2), int(n / 2 + 1)) + elements = int(np.ceil(n / 2 - 1)) # Round up so that len(output)==n + descending = np.linspace(-elements, -1, elements) + + return np.concatenate((ascending, descending)) + + def _Pk2(idxs, sigma): + """# Specify the amplitude given the fft coefficents + + Parameters + ---------- + + idxs : 3 by voxel array int + fft indexes + + sigma : float + spatial sigma + + Returns + ---------- + + amplitude : 3 by voxel array + amplitude of the fft coefficients + """ + + # The first set of idxs ought to be zero so make the first value + # zero to avoid a divide by zero error + amp_start = np.array((0)) + + # Compute the amplitude of the function for a series of indices + amp_end = np.sqrt(np.sqrt(np.sum(idxs[:, 1:] ** 2, 0)) ** (-1 * sigma)) + amplitude = np.append(amp_start, amp_end) + + # Return the output + return amplitude + # Convert from fwhm to sigma (relationship discovered empirical, only an # approximation up to sigma = 0 -> 5 which corresponds to fwhm = 0 -> 8, # relies on an assumption of brain size). - spatial_sigma = logfunc(fwhm, -0.36778719, 2.10601011, 2.15439247) + spatial_sigma = _logfunc(fwhm, -0.36778719, 2.10601011, 2.15439247) - # Set up the input to the fast fourier transform - def fftIndgen(n): - a = list(range(0, int(n / 2 + 1))) - b = list(range(1, int(n / 2))) - b.reverse() - b = [-i for i in b] - return a + b + noise = np.fft.fftn(np.random.normal(size=dimensions)) - # Take in an array of fft values and determine the amplitude for those - # values - def Pk2(idxs): + # Create a meshgrid of the object + fft_vol = np.meshgrid(_fftIndgen(dimensions[0]), _fftIndgen(dimensions[1]), + _fftIndgen(dimensions[2])) - # If all the indexes are zero then set the out put to zero - if np.all(idxs == 0): - return 0.0 - return np.sqrt(np.sqrt(np.sum(idxs ** 2)) ** (-1 * spatial_sigma)) + # Reshape the data into a vector + fft_vec = np.asarray((fft_vol[0].flatten(), fft_vol[1].flatten(), fft_vol[ + 2].flatten())) - noise = np.fft.fftn(np.random.normal(size=dimensions)) - amplitude = np.zeros(dimensions) + # Compute the amplitude for each element in the grid + amp_vec = _Pk2(fft_vec, spatial_sigma) - for x, fft_x in enumerate(fftIndgen(dimensions[0])): - for y, fft_y in enumerate(fftIndgen(dimensions[1])): - for z, fft_z in enumerate(fftIndgen(dimensions[2])): - amplitude[x, y, z] = Pk2(np.array([fft_x, fft_y, fft_z])) + # Reshape to be a brain volume + amplitude = amp_vec.reshape(dimensions) - # The output + # Inverse FFT of the noise plus amplitude noise_spatial = np.fft.ifftn(noise * amplitude) # Mask or not, then z score - if mask is not None and template is not None: + if mask is not None: # Mask the output - noise_spatial = noise_spatial.real * template + noise_spatial = noise_spatial.real * mask # Z score the specific to the brain noise_spatial[mask > 0] = stats.zscore(noise_spatial[mask > 0]) @@ -1759,11 +1940,7 @@ def _generate_noise_temporal(stimfunction_tr, dimensions, template, mask, - fwhm, - motion_sigma, - drift_sigma, - auto_reg_sigma, - physiological_sigma, + noise_dict ): """Generate the temporal noise Generate the time course of the average brain voxel. To change the @@ -1790,27 +1967,15 @@ def _generate_noise_temporal(stimfunction_tr, mask : 3 dimensional array, binary The masked brain, thresholded to distinguish brain and non-brain - fwhm : float - What is the full width half max of the gaussian fields being created - to model spatial noise. - - motion_sigma : float - This is noise that only occurs for the task events, potentially - representing something like noise due to motion - - drift_sigma : float - - What is the sigma on the size of the sine wave - - auto_reg_sigma : float, list - How large is the sigma on the autocorrelation. Higher means more - variable over time. If there are multiple entries then this is - inferred as higher orders of the autoregression - - physiological_sigma : float - - How variable is the signal as a result of physiology, - like heart beat and breathing + noise_dict : dict + A dictionary specifying the types of noise in this experiment. The + noise types interact in important ways. First, all noise types + ending with sigma (e.g. motion sigma) are mixed together in + _generate_temporal_noise. The sigma values describe the proportion of + mixing of these elements. However critically, SFNR is the + parameter that describes how much noise these components contribute + to the brain. If you set the noise dict to matched then it will fit + the parameters to match the participant as best as possible. Returns ---------- @@ -1827,79 +1992,90 @@ def _generate_noise_temporal(stimfunction_tr, # What time points are sampled by a TR? timepoints = list(np.linspace(0, (trs - 1) * tr_duration, trs)) - noise_drift = _generate_noise_temporal_drift(trs, - tr_duration, - ) + # Preset the volume + noise_volume = np.zeros((dimensions[0], dimensions[1], dimensions[2], trs)) - noise_phys = _generate_noise_temporal_phys(timepoints, + # Generate the drift noise + if noise_dict['drift_sigma'] != 0: + # Calculate the drift time course + noise = _generate_noise_temporal_drift(trs, + tr_duration, ) + # Create a volume with the drift properties + volume = np.ones(dimensions) - noise_autoregression = _generate_noise_temporal_autoregression(timepoints, - ) + # Combine the volume and noise + noise_volume += np.multiply.outer(volume, noise) * noise_dict[ + 'drift_sigma'] - # Generate the volumes that will differ depending on the type of noise - # that it will be used for. For drift you want the volume to not have - # the shape of the brain, for the other types of noise you want them to - # have brain shapes - volume_drift = np.ones(dimensions) + # Generate the physiological noise + if noise_dict['physiological_sigma'] != 0: - volume_phys = _generate_noise_spatial(dimensions=dimensions, - template=template, - mask=mask, - fwhm=fwhm, - ) + # Calculate the physiological time course + noise = _generate_noise_temporal_phys(timepoints, + ) - volume_autoreg = _generate_noise_spatial(dimensions=dimensions, - template=template, - mask=mask, - fwhm=fwhm, - ) - - # Multiply the noise by the spatial volume - noise_drift_volume = np.multiply.outer(volume_drift, noise_drift) - noise_phys_volume = np.multiply.outer(volume_phys, noise_phys) - noise_autoregression_volume = np.multiply.outer(volume_autoreg, - noise_autoregression) - - # Sum the noise (it would have been nice to just add all of them in a - # single line but this was causing formatting problems) - noise_temporal = noise_drift_volume * drift_sigma - noise_temporal = noise_temporal + (noise_phys_volume * physiological_sigma) - noise_temporal = noise_temporal + (noise_autoregression_volume * - auto_reg_sigma) - - # Only do this if you are making motion variance - if motion_sigma != 0 and np.sum(stimfunction_tr) > 0: - # Make each noise type - noise_task = _generate_noise_temporal_task(stimfunction_tr, - ) - volume_task = _generate_noise_spatial(dimensions=dimensions, - template=template, - mask=mask, - fwhm=fwhm, + # Create a brain shaped volume with similar smoothing properties + volume = _generate_noise_spatial(dimensions=dimensions, + mask=mask, + fwhm=noise_dict['fwhm'], + ) + + # Combine the volume and noise + noise_volume += np.multiply.outer(volume, noise) * noise_dict[ + 'physiological_sigma'] + + # Generate the AR noise + if noise_dict['auto_reg_sigma'] != 0: + + # Calculate the AR time course volume + noise = _generate_noise_temporal_autoregression(timepoints, + noise_dict, + dimensions, + mask, + ) + + # Combine the volume and noise + noise_volume += noise * noise_dict['auto_reg_sigma'] + + # Generate the task related noise + if noise_dict['task_sigma'] != 0 and np.sum(stimfunction_tr) > 0: + + # Calculate the task based noise time course + noise = _generate_noise_temporal_task(stimfunction_tr, ) - noise_task_volume = np.multiply.outer(volume_task, noise_task) - noise_temporal = noise_temporal + (noise_task_volume * motion_sigma) + + # Create a brain shaped volume with similar smoothing properties + volume = _generate_noise_spatial(dimensions=dimensions, + mask=mask, + fwhm=noise_dict['fwhm'], + ) + # Combine the volume and noise + noise_volume += np.multiply.outer(volume, noise) * noise_dict[ + 'task_sigma'] # Finally, z score each voxel so things mix nicely - noise_temporal = stats.zscore(noise_temporal, 3) + noise_volume = stats.zscore(noise_volume, 3) # If it is a nan it is because you just divided by zero (since some # voxels are zeros in the template) - noise_temporal[np.isnan(noise_temporal)] = 0 + noise_volume[np.isnan(noise_volume)] = 0 - return noise_temporal + return noise_volume def mask_brain(volume, template_name=None, mask_threshold=None, - mask_self=0, + mask_self=True, ): """ Mask the simulated volume - This creates a mask specifying the likelihood (kind of) a voxel is + This creates a mask specifying the approximate likelihood that a voxel is part of the brain. All values are bounded to the range of 0 to 1. An - appropriate threshold to isolate brain voxels is >0.2 + appropriate threshold to isolate brain voxels is >0.2. Critically, + the data that should be used to create a template shouldn't already be + masked/skull stripped. If it is then it will give in accurate estimates + of non-brain noise and corrupt estimations of SNR. Parameters ---------- @@ -1910,8 +2086,8 @@ def mask_brain(volume, template_name : str What is the path to the template to be loaded? If empty then it - defaults to an MNI152 grey matter mask. This is ignored if mask_self is - True. + defaults to an MNI152 grey matter mask. This is ignored if mask_self + is True. mask_threshold : float What is the threshold (0 -> 1) for including a voxel in the mask? If @@ -1920,9 +2096,10 @@ def mask_brain(volume, the minima before that peak as the threshold. Won't work when the data is not bimodal. - mask_self : bool + mask_self : bool or None If set to true then it makes a mask from the volume supplied (by - averaging across time points and changing the range). + averaging across time points and changing the range). If it is set + to false then it will use the template_name as an input. Returns ---------- @@ -1949,7 +2126,7 @@ def mask_brain(volume, else: mask_raw = np.load(template_name) - # Make the masks 3d + # Make the masks 3dremove_baseline if len(mask_raw.shape) == 3: mask_raw = np.array(mask_raw) elif len(mask_raw.shape) == 4 and mask_raw.shape[3] == 1: @@ -2037,7 +2214,34 @@ def _noise_dict_update(noise_dict): _generate_temporal_noise. These values describe the proportion of mixing of these elements. However critically, SFNR is the parameter that describes how much noise these components contribute - to the brain. + to the brain. If you set the noise dict to matched then it will fit + the parameters to match the participant as best as possible. + The noise variables are as follows: + + snr [float]: Ratio of MR signal to the spatial noise + sfnr [float]: Ratio of the MR signal to the temporal noise. This is the + total variability that the following sigmas 'sum' to: + + task_sigma [float]: Size of the variance of task specific noise + drift_sigma [float]: Size of the variance of drift noise + auto_reg_sigma [float]: Size of the variance of autoregressive + noise. This is an ARMA process where the AR and MA components can be + separately specified + physiological_sigma [float]: Size of the variance of physiological + noise + + auto_reg_rho [list]: The coefficients of the autoregressive + components you are modeling + ma_rho [list]:The coefficients of the moving average components you + are modeling + max_activity [float]: The max value of the averaged brain in order + to reference the template + voxel_size [list]: The mm size of the voxels + fwhm [float]: The gaussian smoothing kernel size (mm) + matched [bool]: Specify whether you are fitting the noise parameters + + The volumes of brain noise that are generated have smoothness + specified by 'fwhm' Returns ------- @@ -2046,38 +2250,345 @@ def _noise_dict_update(noise_dict): Updated dictionary """ + # Create the default dictionary + default_dict = {'task_sigma': 0, 'drift_sigma': 0, 'auto_reg_sigma': 1, + 'auto_reg_rho': [0.5], 'ma_rho': [0.0], + 'physiological_sigma': 0, 'sfnr': 90, 'snr': 50, + 'max_activity': 1000, 'voxel_size': [1.0, 1.0, 1.0], + 'fwhm': 4, 'matched': 1} # Check what noise is in the dictionary and add if necessary. Numbers # determine relative proportion of noise - - if 'motion_sigma' not in noise_dict: - noise_dict['motion_sigma'] = 0 - if 'drift_sigma' not in noise_dict: - noise_dict['drift_sigma'] = 0.45 - if 'auto_reg_sigma' not in noise_dict: - noise_dict['auto_reg_sigma'] = 0.45 - if 'physiological_sigma' not in noise_dict: - noise_dict['physiological_sigma'] = 0.1 - if 'sfnr' not in noise_dict: - noise_dict['sfnr'] = 30 - if 'snr' not in noise_dict: - noise_dict['snr'] = 30 - if 'max_activity' not in noise_dict: - noise_dict['max_activity'] = 1000 - if 'voxel_size' not in noise_dict: - noise_dict['voxel_size'] = [1.0, 1.0, 1.0] - if 'fwhm' not in noise_dict: - noise_dict['fwhm'] = 4 + for default_key in default_dict: + if default_key not in noise_dict: + noise_dict[default_key] = default_dict[default_key] return noise_dict +def _fit_spatial(noise, + noise_temporal, + mask, + template, + spatial_sd, + temporal_sd, + noise_dict, + fit_thresh, + fit_delta, + iterations, + ): + """ + Fit the noise model to match the SNR of the data + + Parameters + ---------- + + noise : multidimensional array, float + Initial estimate of the noise + + noise_temporal : multidimensional array, float + The temporal noise that was generated by _generate_temporal_noise + + tr_duration : float + What is the duration, in seconds, of each TR? + + template : 3d array, float + A continuous (0 -> 1) volume describing the likelihood a voxel + is in the brain. This can be used to contrast the brain and non + brain. + + mask : 3d array, binary + The mask of the brain volume, distinguishing brain from non-brain + + spatial_sd : float + What is the standard deviation in space of the noise volume to be + generated + + temporal_sd : float + What is the standard deviation in time of the noise volume to be + generated + + noise_dict : dict + A dictionary specifying the types of noise in this experiment. The + noise types interact in important ways. First, all noise types + ending with sigma (e.g. motion sigma) are mixed together in + _generate_temporal_noise. These values describe the proportion of + mixing of these elements. However critically, SFNR is the + parameter that describes how much noise these components contribute + to the brain. If you set the noise dict to matched then it will + fit the parameters to match the participant as best as possible. + + fit_thresh : float + What proportion of the target parameter value is sufficient + error to warrant finishing fit search. + + fit_delta : float + How much are the parameters attenuated during the fitting process, + in terms of the proportion of difference between the target + parameter and the actual parameter + + iterations : int + The first element is how many steps of fitting the SFNR and SNR + values will be performed. Usually converges after < 5. The + second element is the number of iterations for the AR fitting. + This is much more time consuming (has to make a new timecourse + on each iteration) so be careful about setting this appropriately. + + Returns + ------- + + noise : multidimensional array, float + Generates the noise volume given these parameters + + """ + + # Pull out information that is needed + dim_tr = noise.shape + base = template * noise_dict['max_activity'] + base = base.reshape(dim_tr[0], dim_tr[1], dim_tr[2], 1) + mean_signal = (base[mask > 0]).mean() + target_snr = noise_dict['snr'] + + # Iterate through different parameters to fit SNR and SFNR + spat_sd_orig = np.copy(spatial_sd) + iteration = 0 + for iteration in list(range(iterations)): + + # Calculate the new metrics + new_snr = _calc_snr(noise, mask) + + # Calculate the difference between the real and simulated data + diff_snr = abs(new_snr - target_snr) / target_snr + + # If the AR is sufficiently close then break the loop + if diff_snr < fit_thresh: + logger.info('Terminated SNR fit after ' + str( + iteration) + ' iterations.') + break + + # Convert the SFNR and SNR + spat_sd_new = mean_signal / new_snr + + # Update the variable + spatial_sd -= ((spat_sd_new - spat_sd_orig) * fit_delta) + + # Prevent these going out of range + if spatial_sd < 0 or np.isnan(spatial_sd): + spatial_sd = 10e-3 + + # Set up the machine noise + noise_system = _generate_noise_system(dimensions_tr=dim_tr, + spatial_sd=spatial_sd, + temporal_sd=temporal_sd, + ) + + # Sum up the noise of the brain + noise = base + (noise_temporal * temporal_sd) + noise_system + + # Reject negative values (only happens outside of the brain) + noise[noise < 0] = 0 + + # Failed to converge + if iterations == 0: + logger.info('No fitting iterations were run') + elif iteration == iterations: + logger.warning('SNR failed to converge.') + + # Return the updated noise + return noise, spatial_sd + + +def _fit_temporal(noise, + mask, + template, + stimfunction_tr, + tr_duration, + spatial_sd, + temporal_proportion, + temporal_sd, + noise_dict, + fit_thresh, + fit_delta, + iterations, + ): + """ + Fit the noise model to match the SFNR and AR of the data + + Parameters + ---------- + + noise : multidimensional array, float + Initial estimate of the noise + + mask : 3d array, binary + The mask of the brain volume, distinguishing brain from non-brain + + template : 3d array, float + A continuous (0 -> 1) volume describing the likelihood a voxel + is in the brain. This can be used to contrast the brain and non + brain. + + stimfunction_tr : Iterable, list + When do the stimuli events occur. Each element is a TR + + tr_duration : float + What is the duration, in seconds, of each TR? + + spatial_sd : float + What is the standard deviation in space of the noise volume to be + generated + + temporal_proportion, float + What is the proportion of the temporal variance (as specified by + the SFNR noise parameter) that is accounted for by the system + noise. If this number is high then all of the temporal + variability is due to system noise, if it is low then all of the + temporal variability is due to brain variability. + + temporal_sd : float + What is the standard deviation in time of the noise volume to be + generated + + noise_dict : dict + A dictionary specifying the types of noise in this experiment. The + noise types interact in important ways. First, all noise types + ending with sigma (e.g. motion sigma) are mixed together in + _generate_temporal_noise. These values describe the proportion of + mixing of these elements. However critically, SFNR is the + parameter that describes how much noise these components contribute + to the brain. If you set the noise dict to matched then it will + fit the parameters to match the participant as best as possible. + + fit_thresh : float + What proportion of the target parameter value is sufficient + error to warrant finishing fit search. + + fit_delta : float + How much are the parameters attenuated during the fitting process, + in terms of the proportion of difference between the target + parameter and the actual parameter + + iterations : list, int + The first element is how many steps of fitting the SFNR and SNR + values will be performed. Usually converges after < 5. The + second element is the number of iterations for the AR fitting. + This is much more time consuming (has to make a new timecourse + on each iteration) so be careful about setting this appropriately. + + Returns + ------- + + noise : multidimensional array, float + Generates the noise volume given these parameters + + """ + + # Pull out the + dim_tr = noise.shape + dim = dim_tr[0:3] + base = template * noise_dict['max_activity'] + base = base.reshape(dim[0], dim[1], dim[2], 1) + mean_signal = (base[mask > 0]).mean() + + # Iterate through different parameters to fit SNR and SFNR + temp_sd_orig = np.copy(temporal_sd) + + # Make a copy of the dictionary so it can be modified + new_nd = copy.deepcopy(noise_dict) + + # What SFNR do you want + target_sfnr = noise_dict['sfnr'] + + # What AR do you want? + target_ar = noise_dict['auto_reg_rho'][0] + + # Iterate through different MA parameters to fit AR + for iteration in list(range(iterations)): + + # If there are iterations left to perform then recalculate the + # metrics and try again + + # Calculate the new SFNR + new_sfnr = _calc_sfnr(noise, mask) + + # Calculate the AR + new_ar, _ = _calc_ARMA_noise(noise, + mask, + len(noise_dict['auto_reg_rho']), + len(noise_dict['ma_rho']), + ) + + # Calculate the difference between the real and simulated data + sfnr_diff = abs(new_sfnr - target_sfnr) / target_sfnr + + # Calculate the difference in the first AR component + ar_diff = new_ar[0] - target_ar + + # If the SFNR and AR is sufficiently close then break the loop + if (abs(ar_diff) / target_ar) < fit_thresh and sfnr_diff < fit_thresh: + msg = 'Terminated AR fit after ' + str(iteration) + ' iterations.' + logger.info(msg) + break + + # Otherwise update the noise metrics. Get the new temporal noise value + temp_sd_new = mean_signal / new_sfnr + temporal_sd -= ((temp_sd_new - temp_sd_orig) * fit_delta) + + # Prevent these going out of range + if temporal_sd < 0 or np.isnan(temporal_sd): + temporal_sd = 10e-3 + + # Set the new system noise + temp_sd_system_new = np.sqrt((temporal_sd ** 2) * temporal_proportion) + + # Get the new AR value + new_nd['auto_reg_rho'][0] -= (ar_diff * fit_delta) + + # Don't let the AR coefficient exceed 1 + if new_nd['auto_reg_rho'][0] >= 1: + new_nd['auto_reg_rho'][0] = 0.99 + + # Generate the noise. The appropriate + noise_temporal = _generate_noise_temporal(stimfunction_tr, + tr_duration, + dim, + template, + mask, + new_nd, + ) + + # Set up the machine noise + noise_system = _generate_noise_system(dimensions_tr=dim_tr, + spatial_sd=spatial_sd, + temporal_sd=temp_sd_system_new, + ) + + # Sum up the noise of the brain + noise = base + (noise_temporal * temporal_sd) + noise_system + + # Reject negative values (only happens outside of the brain) + noise[noise < 0] = 0 + + # Failed to converge + if iterations == 0: + logger.info('No fitting iterations were run') + elif iteration == iterations: + logger.warning('AR failed to converge.') + + # Return the updated noise + return noise + + def generate_noise(dimensions, stimfunction_tr, tr_duration, template, mask=None, noise_dict=None, + temporal_proportion=0.5, + iterations=None, + fit_thresh=0.05, + fit_delta=0.5, ): """ Generate the noise to be added to the signal. Default noise parameters will create a noise volume with a standard @@ -2106,8 +2617,57 @@ def generate_noise(dimensions, noise_dict : dictionary, float This is a dictionary which describes the noise parameters of the - data. If there are no other variables provided then it will use default - values + data. If there are no other variables provided then it will use + default values. The noise variables are as follows: + + snr [float]: Ratio of MR signal to the spatial noise + sfnr [float]: Ratio of the MR signal to the temporal noise. This is the + total variability that the following sigmas 'sum' to: + + task_sigma [float]: Size of the variance of task specific noise + drift_sigma [float]: Size of the variance of drift noise + auto_reg_sigma [float]: Size of the variance of autoregressive + noise. This is an ARMA process where the AR and MA components can be + separately specified + physiological_sigma [float]: Size of the variance of physiological + noise + + auto_reg_rho [list]: The coefficients of the autoregressive + components you are modeling + ma_rho [list]:The coefficients of the moving average components you + are modeling + max_activity [float]: The max value of the averaged brain in order + to reference the template + voxel_size [list]: The mm size of the voxels + fwhm [float]: The gaussian smoothing kernel size (mm) + matched [bool]: Specify whether you are fitting the noise parameters + + The volumes of brain noise that are generated have smoothness + specified by 'fwhm' + + temporal_proportion, float + What is the proportion of the temporal variance (as specified by the + SFNR noise parameter) that is accounted for by the system noise. If + this number is high then all of the temporal variability is due to + system noise, if it is low then all of the temporal variability is + due to brain variability. + + iterations : list, int + The first element is how many steps of fitting the SFNR and SNR values + will be performed. Usually converges after < 5. The second element + is the number of iterations for the AR fitting. This is much more + time consuming (has to make a new timecourse on each iteration) so + be careful about setting this appropriately. + + fit_thresh : float + What proportion of the target parameter value is sufficient error to + warrant finishing fit search. + + fit_delta : float + How much are the parameters attenuated during the fitting process, + in terms of the proportion of difference between the target + parameter and the actual parameter + Returns ---------- @@ -2117,13 +2677,29 @@ def generate_noise(dimensions, """ + # Check the input data + if template.max() > 1.1: + raise ValueError('Template out of range') + # Change to be an empty dictionary if it is None if noise_dict is None: noise_dict = {} - # Take in the noise dictionary and determine whether + # Take in the noise dictionary and add any missing information noise_dict = _noise_dict_update(noise_dict) + # How many iterations will you perform? If unspecified it will set + # values based on whether you are trying to match noise specifically to + # this participant or just get in the ball park + if iterations is None: + if noise_dict['matched'] == 1: + iterations = [20, 20] + else: + iterations = [0, 0] + + if abs(noise_dict['auto_reg_rho'][0]) - abs(noise_dict['ma_rho'][0]) < 0.1: + logger.warning('ARMA coefs are close, may have trouble fitting') + # What are the dimensions of the volume, including time dimensions_tr = (dimensions[0], dimensions[1], @@ -2134,126 +2710,238 @@ def generate_noise(dimensions, if mask is None: mask = np.ones(dimensions) + # Create the base (this inverts the process to make the template) + base = template * noise_dict['max_activity'] + + # Reshape the base (to be the same size as the volume to be created) + base = base.reshape(dimensions[0], dimensions[1], dimensions[2], 1) + base = np.ones(dimensions_tr) * base + + # What is the mean signal of the non masked voxels in this template? + mean_signal = (base[mask > 0]).mean() + # Generate the noise noise_temporal = _generate_noise_temporal(stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, dimensions=dimensions, template=template, mask=mask, - fwhm=noise_dict[ - 'fwhm'], - motion_sigma=noise_dict[ - 'motion_sigma'], - drift_sigma=noise_dict[ - 'drift_sigma'], - auto_reg_sigma=noise_dict[ - 'auto_reg_sigma'], - physiological_sigma=noise_dict[ - 'physiological_sigma'], + noise_dict=noise_dict, ) - # Create the base (this inverts the process to make the template) - base = template * noise_dict['max_activity'] - - # What is the mean signal of the non masked voxels in this template? - mean_signal = (base[mask > 0]).mean() - # Convert SFNR into the size of the standard deviation of temporal # variability temporal_sd = (mean_signal / noise_dict['sfnr']) - # Calculate the sd that is necessary to be combined with itself in order - # to generate the temporal_sd - temporal_sd_element = np.sqrt(temporal_sd ** 2 / 2) + # Calculate the temporal sd of the system noise (as opposed to the noise + # attributed to the functional variability). + temporal_sd_system = np.sqrt((temporal_sd ** 2) * temporal_proportion) # What is the standard deviation of the background activity - spatial_sd = mean_signal / noise_dict['snr'] + spat_sd = mean_signal / noise_dict['snr'] + spatial_sd = np.sqrt((spat_sd ** 2) * (1 - temporal_proportion)) # Set up the machine noise noise_system = _generate_noise_system(dimensions_tr=dimensions_tr, spatial_sd=spatial_sd, - temporal_sd=temporal_sd_element, + temporal_sd=temporal_sd_system, ) - # Reshape the base (to be the same size as the volume to be created) - base = base.reshape(dimensions[0], dimensions[1], dimensions[2], 1) - base = np.ones(dimensions_tr) * base - # Sum up the noise of the brain - noise = base + (noise_temporal * temporal_sd_element) + noise_system + noise = base + (noise_temporal * temporal_sd) + noise_system # Reject negative values (only happens outside of the brain) noise[noise < 0] = 0 + # Fit the SNR + noise, spatial_sd = _fit_spatial(noise, + noise_temporal, + mask, + template, + spatial_sd, + temporal_sd_system, + noise_dict, + fit_thresh, + fit_delta, + iterations[0], + ) + + # Fit the SFNR and AR noise + noise = _fit_temporal(noise, + mask, + template, + stimfunction_tr, + tr_duration, + spatial_sd, + temporal_proportion, + temporal_sd, + noise_dict, + fit_thresh, + fit_delta, + iterations[1], + ) + + # Return the noise return noise -def plot_brain(fig, - brain, - mask=None, - percentile=99, - ): - """ Display the brain that has been generated with a given threshold - Will display the voxels above the given percentile and then a shadow of - all voxels in the mask +def compute_signal_change(signal_function, + noise_function, + noise_dict, + magnitude, + method='PSC', + ): + """ Rescale the signal to be a given magnitude, based on a specified + metric (e.g. percent signal change). Metrics are heavily inspired by + Welvaert & Rosseel (2013). The rescaling is based on the maximal + activity in the timecourse. Importantly, all values within the + signal_function are scaled to have a min of -1 or max of 1, meaning that + the voxel value will be the same as the magnitude. Parameters ---------- - fig : matplotlib object - The figure to be displayed, generated from matplotlib. import - matplotlib.pyplot as plt; fig = plt.figure() - brain : 3d array - This is a 3d array with the neural data + signal_function : timepoint by voxel array + The signal time course to be altered. This can have + multiple time courses specified as different columns in this + array. Conceivably you could use the output of + generate_stimfunction as the input but the temporal variance + will be incorrect. Critically, different values across voxels are + considered relative to each other, not independently. E.g., if the + voxel has a peak signal twice as high as another voxel's, then this + means that the signal after these transformations will still be + twice as high (according to the metric) in the first voxel relative + to the second + + noise_function : timepoint by voxel numpy array + The time course of noise (a voxel created from generate_noise) + for each voxel specified in signal_function. This is necessary + for computing the mean evoked activity and the noise variability - mask : 3d array - A binary mask describing the location that you want to specify as + noise_dict : dict + A dictionary specifying the types of noise in this experiment. The + noise types interact in important ways. First, all noise types + ending with sigma (e.g. motion sigma) are mixed together in + _generate_temporal_noise. The sigma values describe the proportion of + mixing of these elements. However critically, SFNR is the + parameter that describes how much noise these components contribute + to the brain. If you set the noise dict to matched then it will + fit the parameters to match the participant as best as possible. + + magnitude : list of floats + This specifies the size, in terms of the metric choosen below, + of the signal being generated. This can be a single number, + and thus apply to all signal timecourses, or it can be array and + thus different for each voxel. + + method : str + Select the procedure used to calculate the signal magnitude, + some of which are based on the definitions outlined in Welvaert & + Rosseel (2013): + - 'SFNR': Change proportional to the temporal variability, + as represented by the (desired) SFNR + - 'CNR_Amp/Noise-SD': Signal magnitude relative to the temporal + noise + - 'CNR_Amp2/Noise-Var_dB': Same as above but converted to decibels + - 'CNR_Signal-SD/Noise-SD': Standard deviation in signal + relative to standard deviation in noise + - 'CNR_Signal-Var/Noise-Var_dB': Same as above but converted to + decibels + - 'PSC': Calculate the percent signal change based on the + average activity of the noise (mean / 100 * magnitude) - percentile : float - What percentage of voxels will be included? Based on the values - supplied Returns ---------- - ax : matplotlib object - Object with the information to be plotted + signal_function_scaled : 4d numpy array + The new signal volume with the appropriately set signal change """ - ax = fig.add_subplot(111, projection='3d') + # If you have only one magnitude value, duplicate the magnitude for each + # timecourse you have + assert type(magnitude) is list, '"magnitude" should be a list of floats' + if len(magnitude) == 1: + magnitude *= signal_function.shape[1] - # Threshold the data - threshold = np.percentile(brain.reshape(np.prod(brain.shape[0:3])), - percentile) + # Scale all signals that to have a range of -1 to 1. This is + # so that any values less than this will be scaled appropriately + signal_function /= np.max(np.abs(signal_function)) - # How many voxels exceed a threshold - brain_threshold = np.where(np.abs(brain) > threshold) + # Iterate through the timecourses and calculate the metric + signal_function_scaled = np.zeros(signal_function.shape) + for voxel_counter in range(signal_function.shape[1]): - # Clear the way - ax.clear() + # Pull out the values for this voxel + sig_voxel = signal_function[:, voxel_counter] + noise_voxel = noise_function[:, voxel_counter] + magnitude_voxel = magnitude[voxel_counter] - ax.set_xlim(0, brain.shape[0]) - ax.set_ylim(0, brain.shape[1]) - ax.set_zlim(0, brain.shape[2]) + # Calculate the maximum signal amplitude (likely to be 1, + # but not necessarily) + max_amp = np.max(np.abs(sig_voxel)) - # If a mask is provided then plot this - if mask is not None: - mask_threshold = np.where(np.abs(mask) > 0) - ax.scatter(mask_threshold[0], - mask_threshold[1], - mask_threshold[2], - zdir='z', - c='black', - s=10, - alpha=0.01) - - # Plot the volume - ax.scatter(brain_threshold[0], - brain_threshold[1], - brain_threshold[2], - zdir='z', - c='red', - s=20) - - return ax + # Calculate the scaled time course using the specified method + if method == 'SFNR': + + # How much temporal variation is there, relative to the mean + # activity + temporal_var = noise_voxel.mean() / noise_dict['sfnr'] + + # Multiply the timecourse by the variability metric + new_sig = sig_voxel * (temporal_var * magnitude_voxel) + + elif method == 'CNR_Amp/Noise-SD': + + # What is the standard deviation of the noise + noise_std = np.std(noise_voxel) + + # Multiply the signal timecourse by the the CNR and noise ( + # rearranging eq.) + new_sig = sig_voxel * (magnitude_voxel * noise_std) + + elif method == 'CNR_Amp2/Noise-Var_dB': + + # What is the standard deviation of the noise + noise_std = np.std(noise_voxel) + + # Rearrange the equation to compute the size of signal change in + # decibels + scale = (10 ** (magnitude_voxel / 20)) * noise_std / max_amp + + new_sig = sig_voxel * scale + + elif method == 'CNR_Signal-SD/Noise-SD': + + # What is the standard deviation of the signal and noise + sig_std = np.std(sig_voxel) + noise_std = np.std(noise_voxel) + + # Multiply the signal timecourse by the the CNR and noise ( + # rearranging eq.) + new_sig = sig_voxel * ((magnitude_voxel / max_amp) * noise_std + / sig_std) + + elif method == 'CNR_Signal-Var/Noise-Var_dB': + # What is the standard deviation of the signal and noise + sig_std = np.std(sig_voxel) + noise_std = np.std(noise_voxel) + + # Rearrange the equation to compute the size of signal change in + # decibels + scale = (10 ** (magnitude_voxel / 20)) * noise_std / (max_amp * + sig_std) + + new_sig = sig_voxel * scale + + elif method == 'PSC': + + # What is the average activity divided by percentage + scale = ((noise_voxel.mean() / 100) * magnitude_voxel) + new_sig = sig_voxel * scale + + signal_function_scaled[:, voxel_counter] = new_sig + + # Return the scaled time course + return signal_function_scaled diff --git a/examples/utils/fmrisim_example.py b/examples/utils/fmrisim_example.py deleted file mode 100644 index 1c4bb194c..000000000 --- a/examples/utils/fmrisim_example.py +++ /dev/null @@ -1,191 +0,0 @@ -# Copyright 2016 Intel Corporation -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -"""fMRI Simulator example script - -Example script to generate a run of a participant's data. This generates -data representing a pair of conditions that are then combined - - Authors: Cameron Ellis (Princeton) 2016 -""" -import logging -import numpy as np -from brainiak.utils import fmrisim as sim -import matplotlib.pyplot as plt -from mpl_toolkits.mplot3d import Axes3D # noqa: F401 -import nibabel - -logger = logging.getLogger(__name__) - -# Inputs for generate_signal -dimensions = np.array([64, 64, 36]) # What is the size of the brain -feature_size = [9, 4, 9, 9] -feature_type = ['loop', 'cube', 'cavity', 'sphere'] -coordinates_A = np.array( - [[32, 32, 18], [26, 32, 18], [32, 26, 18], [32, 32, 12]]) -coordinates_B = np.array( - [[32, 32, 18], [38, 32, 18], [32, 38, 18], [32, 32, 24]]) -signal_magnitude = [1, 0.5, 0.25, -1] # In percent signal change - -# Inputs for generate_stimfunction -onsets_A = [10, 30, 50, 70, 90] -onsets_B = [0, 20, 40, 60, 80] -event_durations = [6] -tr_duration = 2 -temporal_res = 1000.0 # How many elements per second are there -duration = 100 - -# Specify a name to save this generated volume. -savename = 'examples/utils/example.nii' - -# Generate a volume representing the location and quality of the signal -volume_signal_A = sim.generate_signal(dimensions=dimensions, - feature_coordinates=coordinates_A, - feature_type=feature_type, - feature_size=feature_size, - signal_magnitude=signal_magnitude, - ) - -volume_signal_B = sim.generate_signal(dimensions=dimensions, - feature_coordinates=coordinates_B, - feature_type=feature_type, - feature_size=feature_size, - signal_magnitude=signal_magnitude, - ) - -# Visualize the signal that was generated for condition A -fig = plt.figure() -sim.plot_brain(fig, - volume_signal_A) -plt.show() - -# Create the time course for the signal to be generated -stimfunction_A = sim.generate_stimfunction(onsets=onsets_A, - event_durations=event_durations, - total_time=duration, - temporal_resolution=temporal_res, - ) - -stimfunction_B = sim.generate_stimfunction(onsets=onsets_B, - event_durations=event_durations, - total_time=duration, - temporal_resolution=temporal_res, - ) - -# Convolve the HRF with the stimulus sequence -signal_function_A = sim.convolve_hrf(stimfunction=stimfunction_A, - tr_duration=tr_duration, - temporal_resolution=temporal_res, - ) - -signal_function_B = sim.convolve_hrf(stimfunction=stimfunction_B, - tr_duration=tr_duration, - temporal_resolution=temporal_res, - ) - -# Multiply the HRF timecourse with the signal -signal_A = sim.apply_signal(signal_function=signal_function_A, - volume_signal=volume_signal_A, - ) - -signal_B = sim.apply_signal(signal_function=signal_function_B, - volume_signal=volume_signal_B, - ) - -# Combine the signals from the two conditions -signal = signal_A + signal_B - -# Combine the stim functions -stimfunction = list(np.add(stimfunction_A, stimfunction_B)) -stimfunction_tr = stimfunction[::int(tr_duration * temporal_res)] - -# Generate the mask of the signal -mask, template = sim.mask_brain(signal, mask_threshold=0.2) - -# Mask the signal to the shape of a brain (attenuates signal according to grey -# matter likelihood) -signal *= mask.reshape(dimensions[0], dimensions[1], dimensions[2], 1) - -# Generate original noise dict for comparison later -orig_noise_dict = sim._noise_dict_update({}) - -# Create the noise volumes (using the default parameters -noise = sim.generate_noise(dimensions=dimensions, - stimfunction_tr=stimfunction_tr, - tr_duration=tr_duration, - mask=mask, - template=template, - noise_dict=orig_noise_dict, - ) - -# Standardize the signal activity to make it percent signal change -mean_act = (mask * orig_noise_dict['max_activity']).sum() / (mask > 0).sum() -signal = signal * mean_act / 100 - -# Combine the signal and the noise -brain = signal + noise - -# Display the brain -fig = plt.figure() -for tr_counter in list(range(0, brain.shape[3])): - - # Get the axis to be plotted - ax = sim.plot_brain(fig, - brain[:, :, :, tr_counter], - mask=mask, - percentile=99.9) - - # Wait for an input - logging.info(tr_counter) - plt.pause(0.5) - -# Save the volume -affine_matrix = np.diag([-1, 1, 1, 1]) # LR gets flipped -brain_nifti = nibabel.Nifti1Image(brain, affine_matrix) # Create a nifti brain -nibabel.save(brain_nifti, savename) - -# Load in the test dataset and generate a random volume based on it - -# Pull out the data and associated data -volume = nibabel.load(savename).get_data() -dimensions = volume.shape[0:3] -total_time = volume.shape[3] * tr_duration -stimfunction = sim.generate_stimfunction(onsets=[], - event_durations=[0], - total_time=total_time, - ) -stimfunction_tr = stimfunction[::int(tr_duration * temporal_res)] - -# Calculate the mask -mask, template = sim.mask_brain(volume=volume, - mask_self=True, - ) - -# Calculate the noise parameters -noise_dict = sim.calc_noise(volume=volume, - mask=mask, - ) - -# Create the noise volumes (using the default parameters -noise = sim.generate_noise(dimensions=dimensions, - tr_duration=tr_duration, - stimfunction_tr=stimfunction_tr, - template=template, - mask=mask, - noise_dict=noise_dict, - ) - -# Create a nifti brain -brain_noise = nibabel.Nifti1Image(noise, affine_matrix) -nibabel.save(brain_noise, 'examples/utils/example2.nii') # Save diff --git a/examples/utils/fmrisim_multivariate_example.ipynb b/examples/utils/fmrisim_multivariate_example.ipynb index e4ca5d68f..176ad9923 100644 --- a/examples/utils/fmrisim_multivariate_example.ipynb +++ b/examples/utils/fmrisim_multivariate_example.ipynb @@ -24,11 +24,11 @@ "# fMRI Simulator example script for multivariate analyses\n", "\n", "Example script to demonstrate fmrisim functionality. This generates\n", - "data for a two condition, event related design in which each condition\n", + "data for a two condition, event-related design in which each condition\n", "evokes different activity within the same voxels. It then runs simple \n", "univariate and multivariate analyses on the data\n", "\n", - "Authors: Cameron Ellis (Yale) 2017\n" + "Authors: Cameron Ellis (Yale) 2018\n" ] }, { @@ -49,11 +49,18 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n", + " from pandas.core import datetools\n" + ] + } + ], "source": [ "%matplotlib notebook\n", "\n", @@ -62,6 +69,7 @@ "import nibabel\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import scipy.ndimage as ndimage\n", "import scipy.spatial.distance as sp_distance\n", "import sklearn.manifold as manifold\n", "import scipy.stats as stats\n", @@ -80,7 +88,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": { "collapsed": true }, @@ -97,19 +105,25 @@ "source": [ "*1.2\tSpecify participant dimensions and resolution*\n", "\n", - "It is possible to manually specify all parameters necessary for fmrisim. However, it is also possible to instead provide an fMRI dataset as input and extract the necessary parameters from that dataset. Such an example is described below and will be followed throughout. Here the size of the volume and the resolution of the voxels within it are determined." + "The size of the volume and the resolution of the voxels must be specified (or extracted from the real data as is the case below)." ] }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(64, 64, 27, 294)\n" + ] + } + ], "source": [ - "dim = volume.shape\n", - "dimsize = nii.header.get_zooms()\n", + "dim = volume.shape # What is the size of the volume\n", + "dimsize = nii.header.get_zooms() # Get voxel dimensions from the nifti header\n", "tr = dimsize[3]\n", "if tr > 100: # If high then these values are likely in ms\n", " tr /= 1000\n", @@ -122,20 +136,20 @@ "source": [ "*1.3 Generate an activity template and a mask*\n", "\n", - "Functions in fmrisim require a continuous map that describes the appropriate average MR value for each voxel in the brain and a mask which specifies voxels in the brain versus voxels outside of the brain. One way to generate both of these volumes is the mask_brain function. At a minimum, this takes as an input the fMRI volume to be simulated. To create the template this volume is averaged over time and bounded to a range from 0 to 1. In other words, voxels with a high value in the template have high activity over time. To create a mask, the template is thresholded. This threshold can be set manually or instead an appropriate value can be determined by looking for the minima between the two first peaks in the histogram of voxel values.\n" + "Functions in fmrisim require a continuous map that describes the appropriate average MR value for each voxel in the brain and a mask which specifies voxels in the brain versus voxels outside of the brain. One way to generate both of these volumes is the mask_brain function. At a minimum, this takes as an input the fMRI volume to be simulated. To create the template this volume is averaged over time and bounded to a range from 0 to 1. In other words, voxels with a high value in the template have high activity over time. To create a mask, the template is thresholded. This threshold can be set manually or instead an appropriate value can be determined by looking for the minima between the two first peaks in the histogram of voxel values. If you would prefer, you could use the [compute_epi_mask](http://nilearn.github.io/modules/generated/nilearn.masking.compute_epi_mask.html) function in nilearn which uses a similar method." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "mask, template = fmrisim.mask_brain(volume=volume, \n", - " mask_self=True,\n", - " )" + " mask_self=True,\n", + " )" ] }, { @@ -144,33 +158,73 @@ "source": [ "*1.4 Determine noise parameters*\n", "\n", - "A critical step in the fmrisim toolbox is determining the noise parameters of the volume to be created. Many noise parameters are available for specification and if any are not set then they will default to reasonable values. As before, it is instead possible to provide raw fMRI data that will be used to estimate these noise parameters. The goal of the noise estimation is to calculate general descriptive statistics about the noise in the brain that are thought to be important. The simulations are thought to be useful for understanding how signals will survive analyses when embedded in realistic neural noise. \n", + "A critical step in the fmrisim toolbox is determining the noise parameters of the volume to be created. Many noise parameters are available for specification and if any are not set then they will default to reasonable values. As mentioned before, it is instead possible to provide raw fMRI data that will be used to estimate these noise parameters. The goal of the noise estimation is to calculate general descriptive statistics about the noise in the brain that are thought to be important. The simulations are then useful for understanding how signals will survive analyses when embedded in realistic neural noise. \n", + "\n", + "Now the disclaimers: the values here are only an estimate and will depend on noise properties combining in the ways assumed. In addition, because of the non-linearity and stochasticity of this simulation, this estimation is not fully invertible: if you generate a dataset with a set of noise parameters it will have similar but not the same noise parameters as a result. Moreover, complex interactions between brain regions that likely better describe brain noise are not modelled here: this toolbox pays no attention to regions of the brain or their interactions. Finally, for best results use raw fMRI because if the data has been preprocessed then assumptions this algorithm makes are likely to be erroneous. For instance, if the brain has been masked then this will eliminate variance in non-brain voxels which will mean that calculations of noise dependent on those voxels as a reference will fail.\n", "\n", - "Now the disclaimers: the values here are only an estimate and will depend on noise properties combining in the ways specified. In addition, because of the non-linearity and stochasticity of this simulation, this estimation is not fully invertible: if you generate a dataset with a set of noise parameters it will have similar but not the same noise parameters as a result. Moreover, complex interactions between brain regions that likely better describe brain noise are not modelled here: this toolbox pays no attention to regions of the brain or their interactions. Finally, for best results use raw fMRI because if the data has been preprocessed then assumptions this algorithm makes are likely to be erroneous. For instance, if the brain has been masked then this will eliminate variance in non-brain voxels which will mean that calculations of noise dependent on those voxels as a reference will fail.\n", + "To ameliorate some of these concerns, it is possible to fit the spatial and temporal noise properties of the data. This iterates over the noise generation process and tunes parameters in order to match those that are provided. This is time consuming (especially for fitting the temporal noise) but is helpful in matching the specified noise properties. \n", "\n", - "This toolbox separates noise in two: spatial noise and temporal noise. To estimate spatial noise both the smoothness and the amount of non-brain noise of the data must be quantified. For smoothness, the Full Width Half Max (FWHM) of the volume is averaged for the X, Y and Z dimension and then averaged across a sample of time points. To calculate the Signal to Noise Ratio the mean activity in brain voxels for the middle time point is divided by the standard deviation in activity across non-brain voxels for that time point. For temporal noise the drift, temporal autocorrelation, and functional variability is estimated. The drift is estimated by averaging all non-brain voxels and looking at the variance in this average across time. This time course is also used to estimate the temporal autoregression by taking the first slope coefficient of an autoregression estimation function from the Nitime package . The Signal to Fluctuation Noise Ratio is calculated by dividing the average activity of voxels in the brain with that voxel’s noise (Friedman & Glover, 2006). That noise is calculated by taking the standard deviation of that voxel over time after it has been detrended with a second order polynomial. Other types of noise can be generated, such as physiological noise, but are not estimated by this function.\n" + "This toolbox separates noise in two: spatial noise and temporal noise. To estimate spatial noise both the smoothness and the amount of non-brain noise of the data must be quantified. For smoothness, the Full Width Half Max (FWHM) of the volume is averaged for the X, Y and Z dimension and then averaged across a sample of time points. To calculate the Signal to Noise Ratio (SNR) the mean activity in brain voxels for the middle time point is divided by the standard deviation in activity across non-brain voxels for that time point. For temporal noise an auto-regressive and moving average (ARMA) process is estimated, along with the overall size of temporal variability. A sample of brain voxels is used to estimate the first AR component and the first MA component of each voxel's activity over time using the statsmodels package. The Signal to Fluctuation Noise Ratio (SFNR) is calculated by dividing the average activity of voxels in the brain with that voxel’s noise (Friedman & Glover, 2006). That noise is calculated by taking the standard deviation of that voxel over time after it has been detrended with a second order polynomial. The SFNR then controls the amount of functional variability. Other types of noise can be generated, such as physiological noise, but are not estimated by this function.\n" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n", + " \"Check mle_retvals\", ConvergenceWarning)\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/tsa/kalmanf/kalmanfilter.py:649: RuntimeWarning: divide by zero encountered in true_divide\n", + " R_mat, T_mat)\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/tsa/tsatools.py:584: RuntimeWarning: overflow encountered in exp\n", + " newparams = ((1-np.exp(-params))/\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/tsa/tsatools.py:585: RuntimeWarning: overflow encountered in exp\n", + " (1+np.exp(-params))).copy()\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/tsa/tsatools.py:585: RuntimeWarning: invalid value encountered in true_divide\n", + " (1+np.exp(-params))).copy()\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/tsa/tsatools.py:586: RuntimeWarning: overflow encountered in exp\n", + " tmp = ((1-np.exp(-params))/\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/tsa/tsatools.py:587: RuntimeWarning: overflow encountered in exp\n", + " (1+np.exp(-params))).copy()\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/tsa/tsatools.py:587: RuntimeWarning: invalid value encountered in true_divide\n", + " (1+np.exp(-params))).copy()\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/base/model.py:473: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available\n", + " 'available', HessianInversionWarning)\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n", + " \"Check mle_retvals\", ConvergenceWarning)\n" + ] + } + ], "source": [ - "noise_dict = {'voxel_size': [dimsize[0], dimsize[1], dimsize[2]]}\n", + "# Calculate the noise parameters from the data. Set it up to be matched.\n", + "noise_dict = {'voxel_size': [dimsize[0], dimsize[1], dimsize[2]], 'matched': 1}\n", "noise_dict = fmrisim.calc_noise(volume=volume,\n", " mask=mask,\n", + " template=template,\n", " noise_dict=noise_dict,\n", " )" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Noise parameters of the data were estimated as follows:\n", + "SNR: 23.1756482\n", + "SFNR: 70.7171164885\n", + "FWHM: 5.65994469633\n" + ] + } + ], "source": [ "print('Noise parameters of the data were estimated as follows:')\n", "print('SNR: ' + str(noise_dict['snr']))\n", @@ -182,114 +236,5186 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### **2. Generate signal**\n", - "\n", - "fmrisim can be used to generate signal in a number of different ways depending on the type of effect being simulated. Several tools are supplied to help with different types of signal that may be required; however, custom scripts may be necessary for unique effects. Below an experiment will be simulated in which two conditions, A and B, evoke different patterns of activity in the same set of voxels in the brain. This pattern does not manifest as a univariate change in voxel activity across voxels but instead each condition evokes a consistent pattern across voxels. These conditions are randomly intermixed trial by trial. This code could be easily changed to instead compare univariate changes evoked by stimuli in different brain regions. " + "### **2. Generate noise**\n", + "fmrisim can generate realistic fMRI noise when supplied with the appropriate inputs. A single function receives these inputs and deals with generating the noise. The necessary inputs are described below; however, the steps performed by this function are also described in detail for clarity." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/cellis/anaconda/lib/python3.6/site-packages/scipy/stats/stats.py:2246: RuntimeWarning: invalid value encountered in true_divide\n", + " np.expand_dims(sstd, axis=axis))\n", + "/Users/cellis/anaconda/lib/python3.6/site-packages/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n", + " \"Check mle_retvals\", ConvergenceWarning)\n" + ] + } + ], + "source": [ + "# Calculate the noise given the parameters\n", + "noise = fmrisim.generate_noise(dimensions=dim[0:3],\n", + " tr_duration=int(tr),\n", + " stimfunction_tr=[0] * dim[3], \n", + " mask=mask,\n", + " template=template,\n", + " noise_dict=noise_dict,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "(-0.5, 63.5, 63.5, -0.5)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Plot spatial noise\n", + "low_spatial = fmrisim._generate_noise_spatial(dim[0:3],\n", + " fwhm=4.0,\n", + " )\n", + "\n", + "high_spatial = fmrisim._generate_noise_spatial(dim[0:3],\n", + " fwhm=1.0,\n", + " )\n", + "plt.figure()\n", + "plt.subplot(1,2,1)\n", + "plt.title('FWHM = 4.0')\n", + "plt.imshow(low_spatial[:, :, 12])\n", + "plt.axis('off')\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('FWHM = 1.0')\n", + "plt.imshow(high_spatial[:, :, 12])\n", + "plt.axis('off')" + ] + }, + { + "cell_type": "code", + "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ - "coordinates = np.array([[24, 24, 24]])\n", - "feature_size = 3\n", - "signal_volume = fmrisim.generate_signal(dimensions=dim[0:3],\n", - " feature_type=['cube'],\n", - " feature_coordinates=coordinates,\n", - " feature_size=[feature_size],\n", - " signal_magnitude=[1],\n", - " )" + "# Create the different types of noise\n", + "total_time = 500\n", + "timepoints = list(range(0, total_time, int(tr)))\n", + "\n", + "drift = fmrisim._generate_noise_temporal_drift(total_time,\n", + " int(tr),\n", + " )\n", + "\n", + "mini_dim = np.array([2, 2, 2])\n", + "autoreg = fmrisim._generate_noise_temporal_autoregression(timepoints,\n", + " noise_dict,\n", + " mini_dim,\n", + " np.ones(mini_dim),\n", + " np.ones(mini_dim),\n", + " )\n", + " \n", + "phys = fmrisim._generate_noise_temporal_phys(timepoints,\n", + " )\n", + "\n", + "stimfunc = np.zeros((int(total_time / tr), 1))\n", + "stimfunc[np.random.randint(0, int(total_time / tr), 50)] = 1\n", + "task = fmrisim._generate_noise_temporal_task(stimfunc,\n", + " )" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Text(0.5,0,'Activity difference')" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Dilate the mask so as to only take voxels far from the brain (performed in calc_noise)\n", + "mask_dilated = ndimage.morphology.binary_dilation(mask, iterations=10)\n", + "\n", + "# Remove all non brain voxels\n", + "system_all = volume[mask_dilated == 0] # Pull out all the non brain voxels in the first TR\n", + "system_baseline = volume - (template.reshape(dim[0], dim[1], dim[2], 1) * noise_dict['max_activity']) # Subtract the baseline before masking\n", + "system_baseline = system_baseline[mask_dilated == 0]\n", + "\n", + "# Plot the distribution of voxels\n", + "plt.figure()\n", + "plt.subplot(1, 3, 1)\n", + "plt.hist(system_all[:,0].flatten(),100)\n", + "plt.title('Non-brain distribution')\n", + "plt.xlabel('Activity')\n", + "plt.ylabel('Frequency')\n", + "\n", + "# Identify a subset of voxels to plot\n", + "idxs = list(range(system_all.shape[0]))\n", + "np.random.shuffle(idxs)\n", + "\n", + "temporal = system_all[idxs[:100], :100]\n", + "plt.subplot(1, 3, 2)\n", + "plt.imshow(temporal)\n", + "plt.xticks([], [])\n", + "plt.yticks([], [])\n", + "plt.ylabel('voxel ID')\n", + "plt.xlabel('time')\n", + "plt.title('Voxel x time')\n", + "\n", + "# Plot the difference\n", + "ax=plt.subplot(1, 3, 3)\n", + "plt.hist(system_baseline[:,0].flatten(),100)\n", + "ax.yaxis.tick_right()\n", + "ax.yaxis.set_label_position(\"right\")\n", + "plt.title('Demeaned non-brain distribution')\n", + "plt.xlabel('Activity difference')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*2.2 Establish effect size*\n", + "*2.3 Combine noise and template*\n", + " \n", + "The template volume is used to estimate the appropriate baseline distribution of MR values. This estimate is then combined with the temporal noise and the system noise to make an estimate of the noise. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*2.4 Fit the data to the noise parameters*\n", + "\n", + "The generate_noise function does its best to estimate the appropriate noise parameters using assumptions about noise sources; however, because of the complexity of these different noise types, it is often wrong. To compensate, fitting is performed in which parameters involved in the noise generation process are changed and the noise metrics are recalculated to see whether those changes helped the fit. Due to their importance, the parameters that can be fit are SNR, SFNR and AR.\n", + "\n", + "The fitting of SNR/SFNR involves reweighting spatial and temporal metrics of noise. This analysis is relatively quick because this reweighting does not require that any timecourses are recreated, only that they are reweighted. At least 10 iterations are recommended because the initial guesses tend to underestimate SFNR and SNR (although the size of this error depends on the data). In the case of fitting the AR, the MA rho is adjusted until the AR is appropriate and in doing so the timecourse needs to be recreated for each iteration. In the noise_dict, one of the keys is 'matched' which is a binary value determining whether any fitting will be done\n", "\n", - "When specifying the signal we must determine the amount of activity change each voxel undergoes. A useful metric for this is the SFNR value determined from noise calculations because it can be used to estimate the variability in the average voxel. This is set up so that the evoked activity is proportional to the mean activity of a voxel. For a univariate effect, to estimate activity with a Cohen’s d of 1, the size of the change must be equivalent to one standard deviation. For multivariate effects the effect size depends on multiple factors including the number of voxels and conditions. Different measures for effect size could also be calculated, such as percent signal change. Note that that this signal change is based on the average voxel. Instead it might be preferable to model signal change based on the mean of each voxel (i.e. the template value)." + "In terms of timing, for a medium size dataset (64x64x27x300 voxels) it takes approximately 2 minutes to generate the data when fitting on a Mac 2014 laptop. " ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 26, "metadata": { "collapsed": true }, "outputs": [], "source": [ - "effect_size = 1\n", - "temporal_sd = template * noise_dict['max_activity'] / noise_dict['sfnr']\n", + "# Compute the noise parameters for the simulated noise\n", + "noise_dict_sim = {'voxel_size': [dimsize[0], dimsize[1], dimsize[2]], 'matched': 1}\n", + "noise_dict_sim = fmrisim.calc_noise(volume=noise,\n", + " mask=mask,\n", + " template=template,\n", + " noise_dict=noise_dict_sim,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Compare noise parameters for the real and simulated noise:\n", + "SNR: 23.18 vs 23.26\n", + "SFNR: 70.72 vs 68.53\n", + "FWHM: 5.66 vs 5.76\n", + "AR: 0.82 vs 0.84\n" + ] + } + ], + "source": [ + "print('Compare noise parameters for the real and simulated noise:')\n", + "print('SNR: %0.2f vs %0.2f' % (noise_dict['snr'], noise_dict_sim['snr']))\n", + "print('SFNR: %0.2f vs %0.2f' % (noise_dict['sfnr'], noise_dict_sim['sfnr']))\n", + "print('FWHM: %0.2f vs %0.2f' % (noise_dict['fwhm'], noise_dict_sim['fwhm']))\n", + "print('AR: %0.2f vs %0.2f' % (noise_dict['auto_reg_rho'][0], noise_dict_sim['auto_reg_rho'][0]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **3. Generate signal**\n", "\n", - "signal_idxs = np.where(signal_volume == 1)\n", + "fmrisim can be used to generate signal in a number of different ways depending on the type of effect being simulated. Several tools are supplied to help with different types of signal that may be required; however, custom scripts may be necessary for unique effects. Below an experiment will be simulated in which two conditions, A and B, evoke different patterns of activity in the same set of voxels in the brain. This pattern does not manifest as a uniform change in activity across voxels but instead each condition evokes a consistent pattern across voxels. These conditions are randomly intermixed trial by trial. This code could be easily changed to instead compare univariate changes evoked by stimuli in different brain regions. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*3.1 Specify which voxels in the brain contain signal*\n", "\n", - "signal_change = np.zeros((int(signal_volume.sum()), 1))\n", - "for idx_counter in list(range(0, int(signal_volume.sum()))):\n", - " x = signal_idxs[0][idx_counter]\n", - " y = signal_idxs[1][idx_counter]\n", - " z = signal_idxs[2][idx_counter]\n", - " signal_change[idx_counter] = effect_size * temporal_sd[x, y, z]" + "fmrisim provides tools to specify certain voxels in the brain that contain signal. The generate_signal function can produce regions of activity in a brain of different shapes, such as cubes, loops and spheres. Alternatively a volume could be loaded in that specifies the signal voxels (e.g. for ROIs from nilearn). The value of each voxel can be specified here, or set to be a random value." + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Create the region of activity where signal will appear\n", + "coordinates = np.array([[21, 21, 21]]) # Where in the brain is the signal\n", + "feature_size = 3 # How big, in voxels, is the size of the ROI\n", + "signal_volume = fmrisim.generate_signal(dimensions=dim[0:3],\n", + " feature_type=['cube'],\n", + " feature_coordinates=coordinates,\n", + " feature_size=[feature_size],\n", + " signal_magnitude=[1],\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/cellis/anaconda/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.\n", + " warnings.warn(message, mplDeprecation, stacklevel=1)\n" + ] + }, + { + "data": { + "text/plain": [ + "Text(0.5,0,'Condition B')" + ] + }, + "execution_count": 42, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ + "# Plot pattern of activity for each condition\n", "plt.figure()\n", "plt.subplot(1,2,1)\n", "plt.imshow(pattern_A)\n", "plt.ylabel('Voxels')\n", "plt.tick_params(which='both', left='off', labelleft='off', bottom='off', labelbottom='off')\n", "plt.xlabel('Condition A')\n", + "\n", "plt.subplot(1,2,2)\n", "plt.imshow(pattern_B)\n", "plt.tick_params(which='both', left='off', labelleft='off', bottom='off', labelbottom='off')\n", @@ -300,29 +5426,42 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.4 Generate event time course*\n", + "*3.3 Generate event time course*\n", "\n", "generate_stimfunction can be used to specify the time points at which task stimulus events occur. The timing of events can be specified by describing the onset and duration of each event. Alternatively, it is possible to provide a path to a 3 column timing file, used by fMRI software packages like FSL, which specifies event onset, duration and weight. \n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ - "event_duration = 2\n", - "isi = 3\n", - "burn_in = 3\n", - "total_time = int(dim[3] * tr) + burn_in\n", - "events = int((total_time - ((event_duration + isi) * 2)) / ((event_duration + isi) * 2)) * 2\n", - "onsets_all = np.linspace(burn_in, events * (event_duration + isi), events) \n", - "np.random.shuffle(onsets_all)\n", - "onsets_A = onsets_all[:int(events / 2)]\n", - "onsets_B = onsets_all[int(events / 2):]\n", - "temporal_res = 1.0 # How many timepoints per second of the stim function are to be generated?\n", + "# Set up stimulus event time course parameters\n", + "event_duration = 2 # How long is each event\n", + "isi = 7 # What is the time between each event\n", + "burn_in = 1 # How long before the first event\n", + "\n", + "total_time = int(dim[3] * tr) + burn_in # How long is the total event time course\n", + "events = int((total_time - ((event_duration + isi) * 2)) / ((event_duration + isi) * 2)) * 2 # How many events are there?\n", + "onsets_all = np.linspace(burn_in, events * (event_duration + isi), events) # Space the events out\n", + "np.random.shuffle(onsets_all) # Shuffle their order\n", + "onsets_A = onsets_all[:int(events / 2)] # Assign the first half of shuffled events to condition A\n", + "onsets_B = onsets_all[int(events / 2):] # Assign the second half of shuffled events to condition B\n", + "temporal_res = 10.0 # How many timepoints per second of the stim function are to be generated?" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Create a time course of events \n", "stimfunc_A = fmrisim.generate_stimfunction(onsets=onsets_A,\n", " event_durations=[event_duration],\n", " total_time=total_time,\n", @@ -340,14 +5479,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.5 Export stimulus time course for analysis*\n", + "*3.4 Export stimulus time course for analysis*\n", "\n", "If a time course of events is generated, as is the case here, it may be useful to store this in a certain format for future analyses. The export_3_column function can be used to export the time course to be a three column (event onset, duration and weight) timing file that might readable to FSL. Alternatively, the export_epoch_file function can be used to export numpy files that are necessary inputs for MVPA and FCMA in BrainIAK.\n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 36, "metadata": { "collapsed": true }, @@ -374,267 +5513,1772 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.6 Estimate the voxel weight for each event*\n", + "*3.5 Estimate the voxel weight for each event*\n", "\n", - "According to the logic of this example, each signal voxel will respond a different amount for condition A and B, but this amount will also differ across voxels. To simulate this we multiply a voxel’s response to each condition by the time course of events and then combine these conditions time courses to make a single time course. This time course describes each voxel’s response to stimuli over time." + "According to the logic of this example, each voxel carrying signal will respond a different amount for condition A and B. To simulate this we multiply a voxel’s response to each condition by the time course of events and then combine these to make a single time course. This time course describes each voxel’s response to signal over time." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 37, "metadata": { "collapsed": true }, "outputs": [], "source": [ + "# Multiply each pattern by each voxel time course\n", "weights_A = np.matlib.repmat(stimfunc_A, 1, voxels).transpose() * pattern_A\n", "weights_B = np.matlib.repmat(stimfunc_B, 1, voxels).transpose() * pattern_B\n", - "weights_all = weights_A + weights_B\n", - "weights_all = weights_all.transpose()" + "\n", + "# Sum these time courses together\n", + "stimfunc_weighted = weights_A + weights_B\n", + "stimfunc_weighted = stimfunc_weighted.transpose()" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Text(0.5,0,'nth TR')" + ] + }, + "execution_count": 45, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ + "# Prepare the data to be plotted\n", + "response = signal_func[0:100,0] * 2\n", + "downsample_A = stimfunc_A[0:int(100*temporal_res * tr):int(temporal_res * tr), 0]\n", + "downsample_B = stimfunc_B[0:int(100*temporal_res * tr):int(temporal_res * tr), 0]\n", + "\n", "# Display signal\n", "plt.figure()\n", - "\n", - "response = stats.zscore(signal_func[0:100,0])\n", "plt.title('Example event time course and voxel response')\n", - "Event_A = plt.plot(stimfunc_A[0:100, 0], 'r', label='Event_A')\n", - "Event_B = plt.plot(stimfunc_B[0:100, 0], 'g', label='Event_B')\n", + "Event_A = plt.plot(downsample_A, 'r', label='Event_A')\n", + "Event_B = plt.plot(downsample_B, 'g', label='Event_B')\n", "Response = plt.plot(response, 'b', label='Response')\n", - "plt.legend(loc=1)" + "plt.legend(loc=1)\n", + "plt.yticks([],'')\n", + "plt.xlabel('nth TR')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*2.8 Multiply the convolved response with the signal voxels*\n", + "*3.7 Establish signal magnitude*\n", "\n", - "If you have a time course of simulated response for one or more voxels and a three dimensional volume representing voxels that ought to respond to these events then apply_signal will combine these appropriately. This function multiplies each signal voxel in the brain by the convolved event time course. \n" + "When specifying the signal we must determine the amount of activity change each voxel undergoes. fmrisim contains a tool to allow you to choose between a variety of different metrics that you could use to scale the signal. For instance, we can calculate percent signal change (referred to as PSC) by taking the average activity of a voxel in the noise volume and multiplying the maximal activation of the signal by a percentage of this number. This metric doesn't take into account the variance in the noise but other metrics in this tool do. One metric that does take account of variance, and is used below, is the signal amplitude divided by the temporal variability. The choices that are available for computing the signal scale are based on Welvaert and Rosseel (2013)." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 46, "metadata": { "collapsed": true }, "outputs": [], "source": [ - "signal = fmrisim.apply_signal(signal_func,\n", - " signal_volume,\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### **3. Add noise to signal**\n", - "fmrisim can generate realistic fMRI noise when supplied with the appropriate inputs. A single function receives these inputs and deals with generating the noise. The necessary inputs are described below; however, the steps performed by this function are also described in detail for clarity.\n" + "# Specify the parameters for signal\n", + "signal_method = 'CNR_Amp/Noise-SD'\n", + "signal_magnitude = [0.5]\n", + "\n", + "# Where in the brain are there stimulus evoked voxels\n", + "signal_idxs = np.where(signal_volume == 1)\n", + "\n", + "# Pull out the voxels corresponding to the noise volume\n", + "noise_func = noise[signal_idxs[0], signal_idxs[1], signal_idxs[2], :]" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 47, "metadata": { "collapsed": true }, "outputs": [], "source": [ - "noise = fmrisim.generate_noise(dimensions=dim[0:3],\n", - " tr_duration=int(tr),\n", - " stimfunction_tr=weights_all[::int(tr)], \n", - " mask=mask,\n", - " template=template,\n", - " noise_dict=noise_dict,\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure()\n", - "plt.imshow(noise[:, :, 24, 0], cmap=plt.cm.gray)\n", - "plt.axis('off')" + "# Compute the signal appropriate scaled\n", + "signal_func_scaled = fmrisim.compute_signal_change(signal_func,\n", + " noise_func,\n", + " noise_dict,\n", + " magnitude=signal_magnitude,\n", + " method=signal_method,\n", + " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*3.1 Create temporal noise*\n", - "\n", - "The temporal noise of fMRI data is comprised of multiple components: drift, autoregression, task related motion and physiological noise. To estimate drift, a sine wave with a default period of 300s, is used. To estimate drift, cosine basis functions are combined, with longer runs being comprised of more basis functions (Welvaert, et al., 2011). This drift is then multiplied by a three-dimensional volume of Gaussian random fields of a specific FWHM. Autoregression noise is estimated by creating a time course of Gaussian noise values that are weighted by previous values of the time course. This autoregressive time course is multiplied by a brain shaped volume of Gaussian random fields. Physiological noise is modeled by sine waves comprised of heart rate (1.17Hz) and respiration rate (0.2Hz) (Biswal, et al., 1996) with random phase. This time course is also multiplied by brain shaped spatial noise. Finally, task related noise is simulated by adding Gaussian or Rician noise to time points where there are events (according to the event time course) and in turn this is multiplied by a brain shaped spatial noise volume. These four noise components are then mixed together in proportion to the size of their corresponding noise values. This aggregated volume is then Z scored and the SFNR is used to estimate the appropriate standard deviation of these values across time. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Plot spatial noise\n", - "low_spatial = fmrisim._generate_noise_spatial(dim[0:3],\n", - " fwhm=4.0,\n", - " )\n", - "\n", - "high_spatial = fmrisim._generate_noise_spatial(dim[0:3],\n", - " fwhm=1.0,\n", - " )\n", - "plt.figure()\n", - "plt.subplot(1,2,1)\n", - "plt.title('Low noise')\n", - "plt.imshow(low_spatial[:, :, 12])\n", - "plt.axis('off')\n", - "\n", - "plt.subplot(1,2,2)\n", - "plt.title('High noise')\n", - "plt.imshow(high_spatial[:, :, 12])\n", - "plt.axis('off')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Create the different types of noise\n", - "timepoints = list(range(0, total_time, int(tr)))\n", - "\n", - "drift = fmrisim._generate_noise_temporal_drift(total_time,\n", - " int(tr),\n", - " )\n", + "*3.8 Multiply the convolved response with the signal voxels*\n", "\n", - "autoreg = fmrisim._generate_noise_temporal_autoregression(timepoints,\n", - " )\n", - " \n", - "phys = fmrisim._generate_noise_temporal_phys(timepoints,\n", - " )\n", - "\n", - "task = fmrisim._generate_noise_temporal_task(abs(stimfunc_A[::int(tr)]),\n", - " )\n", - "\n", - "# Plot the different noise types\n", - "plt.figure()\n", - "plt.title('Noise types')\n", - "\n", - "plt.subplot(4, 1, 1)\n", - "plt.plot(drift)\n", - "plt.axis('off')\n", - "plt.xlabel('Drift')\n", - "\n", - "plt.subplot(4, 1, 2)\n", - "plt.plot(autoreg)\n", - "plt.axis('off')\n", - "plt.xlabel('Autoregression')\n", - "\n", - "plt.subplot(4, 1, 3)\n", - "plt.plot(phys)\n", - "plt.axis('off')\n", - "plt.xlabel('Physiological')\n", - "\n", - "plt.subplot(4, 1, 4)\n", - "plt.plot(task)\n", - "plt.axis('off')\n", - "plt.xlabel('Task')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*3.2 Create system noise*\n", - " \n", - "In addition to temporal noise from fluctuations in the scanner there is also machine noise that causes fluctuations in all voxels. When SNR is low, Rician noise is a good estimate of background noise data (Gudbjartsson, & Patz, 1995). From our testing, when SNR is higher than 30 then noise with an exponential distribution better describes the data. The SNR value that is supplied determines the standard deviation of this machine noise.\t" + "If you have a time course of simulated response for one or more voxels and a three dimensional volume representing voxels that ought to respond to these events then apply_signal will combine these appropriately. This function multiplies each signal voxel in the brain by the convolved event time course. " ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, + "execution_count": 48, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ - "system = fmrisim._generate_noise_system(dimensions_tr=dim,\n", - " spatial_sd=1.5,\n", - " temporal_sd=1,\n", - " )\n", - "\n", - "plt.figure()\n", - "plt.subplot(1, 3, 1)\n", - "plt.hist(system.flatten())\n", - "plt.title('Activity distribution')\n", - "plt.xlabel('Activity')\n", - "plt.ylabel('Frequency')\n", - "\n", - "plt.subplot(1, 3, 2)\n", - "plt.imshow(system[:, :, 0, 0])\n", - "plt.axis('off')\n", - "plt.title('Spatial plane')\n", - "plt.clim([system.min(), np.percentile(system, 95)])\n", - "\n", - "plt.subplot(1, 3, 3)\n", - "plt.imshow(system[0, :64, 0, :64])\n", - "plt.axis('off')\n", - "plt.title('Temporal plane')\n", - "plt.clim([system.min(), np.percentile(system, 95)])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*3.3 Combine noise and template*\n", - " \n", - "The template volume is used to estimate the appropriate baseline distribution of MR values. This estimate is then combined with the temporal noise and the system noise to make an estimate of the noise. " + "signal = fmrisim.apply_signal(signal_func_scaled,\n", + " signal_volume,\n", + " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*3.4 Combine signal and noise*\n", + "*3.9 Combine signal and noise*\n", "\n", "Since the brain signal is expected to be small and sparse relative to the noise, it is assumed sufficient to simply add the volume containing signal with the volume modeling noise to make the simulated brain. " ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 49, "metadata": { "collapsed": true }, @@ -649,7 +7293,7 @@ "source": [ "### **4. Analyse data**\n", "\n", - "Several tools are available for multivariate analysis in BrainIAK. These greatly speed up computation and are critical in some cases, such as a whole brain searchlight. However, for this example data we will only look at data in the ROI that we know contains signal." + "Several tools are available for multivariate analysis in BrainIAK. These greatly speed up computation and are critical in some cases, such as a whole brain searchlight. However, for this example data we will only look at data in the ROI that we know contains signal and so do not need these advanced tools optimized for whole-brain analyses." ] }, { @@ -663,7 +7307,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 50, "metadata": { "collapsed": true }, @@ -675,19 +7319,825 @@ "lb = (coordinates - ((feature_size - 1) / 2)).astype('int')[0]\n", "ub = (coordinates + ((feature_size - 1) / 2) + 1).astype('int')[0]\n", "\n", - "trials_A = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], (onsets_A + hrf_lag / tr).astype('int')]\n", - "trials_B = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], (onsets_B + hrf_lag / tr).astype('int')]\n", + "# Pull out voxels in the ROI for the specified timepoints\n", + "trials_A = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], ((onsets_A + hrf_lag) / tr).astype('int')]\n", + "trials_B = brain[lb[0]:ub[0], lb[1]:ub[1], lb[2]:ub[2], ((onsets_B + hrf_lag) / tr).astype('int')]\n", "\n", + "# Reshape data for easy handling\n", "trials_A = trials_A.reshape((voxels, trials_A.shape[3]))\n", "trials_B = trials_B.reshape((voxels, trials_B.shape[3]))" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 51, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Text(0.5,1,'Low Dimensional Representation of conditions A and B')" + ] + }, + "execution_count": 52, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ + "# Calculate the distance matrix between trial types\n", "distance_matrix = sp_distance.squareform(sp_distance.pdist(np.vstack([trials_A.transpose(), trials_B.transpose()])))\n", "\n", "mds = manifold.MDS(n_components=2, dissimilarity='precomputed') # Fit the mds object\n", "coords = mds.fit(distance_matrix).embedding_ # Find the mds coordinates\n", + "\n", + "# Plot the data\n", "plt.figure()\n", "plt.scatter(coords[:, 0], coords[:, 1], c=['red'] * trials_A.shape[1] + ['green'] * trials_B.shape[1])\n", "plt.axis('off')\n", @@ -734,15 +8990,23 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 53, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean difference between condition A and B: 1.57\n", + "p value: 0.489\n" + ] + } + ], "source": [ "mean_difference = (np.mean(trials_A,0) - np.mean(trials_B,0))\n", "ttest = stats.ttest_1samp(mean_difference, 0)\n", "\n", - "print('Mean difference between condition A and B: ' + str(mean_difference.mean())[0:5])\n", - "print('pvalue: '+ str(ttest.pvalue)[0:5])" + "print('Mean difference between condition A and B: %0.2f\\np value: %0.3f' % (mean_difference.mean(), ttest.pvalue))" ] }, { @@ -756,14 +9020,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 54, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Classification accuracy between condition A and B: 1.000\n" + ] + } + ], "source": [ + "# Get the inputs to the SVM\n", "input_mat = np.vstack([trials_A.transpose(), trials_B.transpose()])\n", "input_labels = trials_A.shape[1] * [1] + trials_B.shape[1] * [0]\n", "\n", - "\n", "# Set up the classifier\n", "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(\n", " input_mat, input_labels, test_size=0.2, random_state=0)\n", @@ -771,7 +9043,7 @@ "clf = sklearn.svm.SVC(kernel='linear', C=1).fit(X_train, y_train)\n", "\n", "score = clf.score(X_test, y_test)\n", - "print('Classification accuracy between condition A and B: ' + str(score)[0:5])" + "print('Classification accuracy between condition A and B: %0.3f' % score)" ] }, { @@ -789,17 +9061,10 @@ "\n", "Gudbjartsson, H. and Patz, S. (1995) The Rician distribution of noisy MRI data. Magnetic resonance in medicine 34, 910-914\n", "\n", - "Welvaert, M., et al. (2011) neuRosim: An R package for generating fMRI data. Journal of Statistical Software 44, 1-18\n" + "Welvaert, M., et al. (2011) neuRosim: An R package for generating fMRI data. Journal of Statistical Software 44, 1-18\n", + "\n", + "Welvaert, M., & Rosseel, Y. (2013). On the definition of signal-to-noise ratio and contrast-to-noise ratio for fMRI data. PloS one, 8(11), e77089.\n" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/setup.py b/setup.py index a0d883ec9..7cb7e972c 100644 --- a/setup.py +++ b/setup.py @@ -127,6 +127,7 @@ def finalize_options(self): 'numpy', 'scikit-learn[alldeps]>=0.18', 'scipy!=1.0.0', # See https://github.com/scipy/scipy/pull/8082 + 'statsmodels', 'pymanopt', 'theano', 'pybind11>=1.7', diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index a1b0f8278..e7254dc23 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -21,6 +21,7 @@ import numpy as np import math from brainiak.utils import fmrisim as sim +import pytest def test_generate_signal(): @@ -29,8 +30,7 @@ def test_generate_signal(): dimensions = np.array([10, 10, 10]) # What is the size of the brain feature_size = [3] feature_type = ['cube'] - feature_coordinates = np.array( - [[5, 5, 5]]) + feature_coordinates = np.array([[5, 5, 5]]) signal_magnitude = [30] # Generate a volume representing the location and quality of the signal @@ -51,6 +51,7 @@ def test_generate_signal(): feature_coordinates = np.array( [[5, 5, 5], [3, 3, 3], [7, 7, 7]]) + # Check feature size is correct volume = sim.generate_signal(dimensions=dimensions, feature_coordinates=feature_coordinates, feature_type=['loop', 'cavity', 'sphere'], @@ -60,6 +61,35 @@ def test_generate_signal(): assert volume[3, 3, 3] == 0, "Cavity is empty" assert volume[7, 7, 7] != 0, "Sphere is not empty" + # Check feature size manipulation + volume = sim.generate_signal(dimensions=dimensions, + feature_coordinates=feature_coordinates, + feature_type=['loop', 'cavity', 'sphere'], + feature_size=[1], + signal_magnitude=signal_magnitude) + assert volume[5, 6, 6] == 0, "Loop is too big" + assert volume[3, 5, 5] == 0, "Cavity is too big" + assert volume[7, 9, 9] == 0, "Sphere is too big" + + # Check that out of bounds feature coordinates are corrected + feature_coordinates = np.array([0, 2, dimensions[2]]) + x, y, z = sim._insert_idxs(feature_coordinates, feature_size[0], + dimensions) + assert x[1] - x[0] == 2, "x min not corrected" + assert y[1] - y[0] == 3, "y was corrected when it shouldn't be" + assert z[1] - z[0] == 1, "z max not corrected" + + # Check that signal patterns are created + feature_coordinates = np.array([[5, 5, 5]]) + volume = sim.generate_signal(dimensions=dimensions, + feature_coordinates=feature_coordinates, + feature_type=feature_type, + feature_size=feature_size, + signal_magnitude=signal_magnitude, + signal_constant=0, + ) + assert volume[4:7, 4:7, 4:7].std() > 0, "Signal is constant" + def test_generate_stimfunction(): @@ -104,6 +134,69 @@ def test_generate_stimfunction(): assert 25 < max_response <= 30, "HRF has the incorrect length" assert np.sum(signal_function < 0) > 0, "No values below zero" + # Export a stimfunction + sim.export_3_column(stimfunction, + 'temp.txt', + ) + + # Load in the stimfunction + stimfunc_new = sim.generate_stimfunction(onsets=None, + event_durations=None, + total_time=duration, + timing_file='temp.txt', + ) + + assert np.all(stimfunc_new == stimfunction), "Export/import failed" + + # Break the timing precision of the generation + stimfunc_new = sim.generate_stimfunction(onsets=None, + event_durations=None, + total_time=duration, + timing_file='temp.txt', + temporal_resolution=0.5, + ) + + assert stimfunc_new.sum() == 0, "Temporal resolution not working right" + + # Set the duration to be too short so you should get an error + onsets = [10, 30, 50, 70, 90] + event_durations = [5] + with pytest.raises(ValueError): + sim.generate_stimfunction(onsets=onsets, + event_durations=event_durations, + total_time=89, + ) + + # Clip the event offset + stimfunc_new = sim.generate_stimfunction(onsets=onsets, + event_durations=event_durations, + total_time=95, + ) + assert stimfunc_new[-1] == 1, 'Event offset was not clipped' + + # Test exporting a group of participants to an epoch file + cond_a = sim.generate_stimfunction(onsets=onsets, + event_durations=event_durations, + total_time=110, + ) + + cond_b = sim.generate_stimfunction(onsets=[x + 5 for x in onsets], + event_durations=event_durations, + total_time=110, + ) + + stimfunction_group = [np.hstack((cond_a, cond_b))] * 2 + sim.export_epoch_file(stimfunction_group, + 'temp.txt', + tr_duration, + ) + + # Check that convolve throws a warning when the shape is wrong + sim.convolve_hrf(stimfunction=np.hstack((cond_a, cond_b)).T, + tr_duration=tr_duration, + temporal_resolution=1, + ) + def test_apply_signal(): @@ -138,6 +231,140 @@ def test_apply_signal(): tr_duration=tr_duration, ) + # Check that you can compute signal change appropriately + # Preset a bunch of things + stimfunction_tr = stimfunction[::int(tr_duration * 100)] + mask, template = sim.mask_brain(dimensions, mask_self=False) + noise_dict = sim._noise_dict_update({}) + noise = sim.generate_noise(dimensions=dimensions, + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict=noise_dict, + iterations=[0, 0] + ) + coords = feature_coordinates[0] + noise_function_a = noise[coords[0], coords[1], coords[2], :] + noise_function_a = noise_function_a.reshape(duration // tr_duration, 1) + + noise_function_b = noise[coords[0] + 1, coords[1], coords[2], :] + noise_function_b = noise_function_b.reshape(duration // tr_duration, 1) + + # Create the calibrated signal with PSC + method = 'PSC' + sig_a = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [0.5], + method, + ) + sig_b = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [1.0], + method, + ) + + assert sig_b.max() / sig_a.max() == 2, 'PSC modulation failed' + + # Create the calibrated signal with SFNR + method = 'SFNR' + sig_a = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [0.5], + method, + ) + scaled_a = sig_a / (noise_function_a.mean() / noise_dict['sfnr']) + sig_b = sim.compute_signal_change(signal_function, + noise_function_b, + noise_dict, + [1.0], + method, + ) + scaled_b = sig_b / (noise_function_b.mean() / noise_dict['sfnr']) + + assert scaled_b.max() / scaled_a.max() == 2, 'SFNR modulation failed' + + # Create the calibrated signal with CNR_Amp/Noise-SD + method = 'CNR_Amp/Noise-SD' + sig_a = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [0.5], + method, + ) + scaled_a = sig_a / noise_function_a.std() + sig_b = sim.compute_signal_change(signal_function, + noise_function_b, + noise_dict, + [1.0], + method, + ) + scaled_b = sig_b / noise_function_b.std() + + assert scaled_b.max() / scaled_a.max() == 2, 'CNR_Amp modulation failed' + + # Create the calibrated signal with CNR_Amp/Noise-Var_dB + method = 'CNR_Amp2/Noise-Var_dB' + sig_a = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [0.5], + method, + ) + scaled_a = np.log(sig_a.max() / noise_function_a.std()) + sig_b = sim.compute_signal_change(signal_function, + noise_function_b, + noise_dict, + [1.0], + method, + ) + scaled_b = np.log(sig_b.max() / noise_function_b.std()) + + assert np.round(scaled_b / scaled_a) == 2, 'CNR_Amp dB modulation failed' + + # Create the calibrated signal with CNR_Signal-SD/Noise-SD + method = 'CNR_Signal-SD/Noise-SD' + sig_a = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [0.5], + method, + ) + scaled_a = sig_a.std() / noise_function_a.std() + sig_b = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [1.0], + method, + ) + scaled_b = sig_b.std() / noise_function_a.std() + + assert (scaled_b / scaled_a) == 2, 'CNR signal modulation failed' + + # Create the calibrated signal with CNR_Amp/Noise-Var_dB + method = 'CNR_Signal-Var/Noise-Var_dB' + sig_a = sim.compute_signal_change(signal_function, + noise_function_a, + noise_dict, + [0.5], + method, + ) + + scaled_a = np.log(sig_a.std() / noise_function_a.std()) + sig_b = sim.compute_signal_change(signal_function, + noise_function_b, + noise_dict, + [1.0], + method, + ) + scaled_b = np.log(sig_b.std() / noise_function_b.std()) + + assert np.round(scaled_b / scaled_a) == 2, 'CNR signal dB modulation ' \ + 'failed' + # Convolve the HRF with the stimulus sequence signal = sim.apply_signal(signal_function=signal_function, volume_signal=volume, @@ -153,6 +380,15 @@ def test_apply_signal(): assert np.any(signal == signal_magnitude), "The stimfunction is not binary" + # Check that there is an error if the number of signal voxels doesn't + # match the number of non zero brain voxels + with pytest.raises(IndexError): + sig_vox = (volume > 0).sum() + vox_pattern = np.tile(stimfunction, (1, sig_vox - 1)) + sim.apply_signal(signal_function=vox_pattern, + volume_signal=volume, + ) + def test_generate_noise(): @@ -193,36 +429,167 @@ def test_generate_noise(): ) # Generate the mask of the signal - mask, template = sim.mask_brain(signal, mask_threshold=0.1) + mask, template = sim.mask_brain(signal, + mask_self=None) assert min(mask[mask > 0]) > 0.1, "Mask thresholding did not work" assert len(np.unique(template) > 2), "Template creation did not work" stimfunction_tr = stimfunction[::int(tr_duration * 100)] + # Create the noise volumes (using the default parameters) noise = sim.generate_noise(dimensions=dimensions, stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, template=template, mask=mask, + iterations=[1, 0], ) assert signal.shape == noise.shape, "The dimensions of signal and noise " \ "the same" - assert np.std(signal) < np.std(noise), "Noise was not created" - - noise = sim.generate_noise(dimensions=dimensions, - stimfunction_tr=stimfunction_tr, - tr_duration=tr_duration, - template=template, - mask=mask, - noise_dict={'sfnr': 10000, 'snr': 10000}, - ) - - system_noise = np.std(noise[mask > 0], 1).mean() + noise_high = sim.generate_noise(dimensions=dimensions, + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict={'sfnr': 50, 'snr': 25}, + iterations=[1, 0], + ) + + noise_low = sim.generate_noise(dimensions=dimensions, + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict={'sfnr': 100, 'snr': 25}, + iterations=[1, 0], + ) + + system_high = np.std(noise_high[mask > 0], 1).mean() + system_low = np.std(noise_low[mask > 0], 1).mean() + + assert system_low < system_high, "SFNR noise could not be manipulated" + + # Check that you check for the appropriate template values + with pytest.raises(ValueError): + sim.generate_noise(dimensions=dimensions, + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template * 2, + mask=mask, + noise_dict={}, + ) + + # Check that iterations does what it should + sim.generate_noise(dimensions=dimensions, + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict={}, + iterations=[0, 0], + ) + + sim.generate_noise(dimensions=dimensions, + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict={}, + iterations=None, + ) + + # Test drift noise + trs = 1000 + period = 100 + drift = sim._generate_noise_temporal_drift(trs, + tr_duration, + 'sine', + period, + ) + + # Check that the max frequency is the appropriate frequency + power = abs(np.fft.fft(drift))[1:trs // 2] + freq = np.linspace(1, trs // 2 - 1, trs // 2 - 1) / trs + period_freq = np.where(freq == 1 / (period // tr_duration)) + max_freq = np.argmax(power) + + assert period_freq == max_freq, 'Max frequency is not where it should be' + + # Do the same but now with cosine basis functions, answer should be close + drift = sim._generate_noise_temporal_drift(trs, + tr_duration, + 'discrete_cos', + period, + ) + + # Check that the appropriate frequency is peaky (may not be the max) + power = abs(np.fft.fft(drift))[1:trs // 2] + freq = np.linspace(1, trs // 2 - 1, trs // 2 - 1) / trs + period_freq = np.where(freq == 1 / (period // tr_duration))[0][0] + + assert power[period_freq] > power[period_freq + 1], 'Power is low' + assert power[period_freq] > power[period_freq - 1], 'Power is low' + + # Check it gives a warning if the duration is too short + drift = sim._generate_noise_temporal_drift(50, + tr_duration, + 'discrete_cos', + period, + ) + + # Test physiological noise (using unrealistic parameters so that it's easy) + timepoints = list(np.linspace(0, (trs - 1) * tr_duration, trs)) + resp_freq = 0.2 + heart_freq = 1.17 + phys = sim._generate_noise_temporal_phys(timepoints, + resp_freq, + heart_freq, + ) - assert system_noise <= 0.1, "Noise strength could not be manipulated" + # Check that the max frequency is the appropriate frequency + power = abs(np.fft.fft(phys))[1:trs // 2] + freq = np.linspace(1, trs // 2 - 1, trs // 2 - 1) / (trs * tr_duration) + peaks = (power > (power.mean() + power.std())) # Where are the peaks + peak_freqs = freq[peaks] + + assert np.any(resp_freq == peak_freqs), 'Resp frequency not found' + assert len(peak_freqs) == 2, 'Two peaks not found' + + # Test task noise + sim._generate_noise_temporal_task(stimfunction_tr, + motion_noise='gaussian', + ) + sim._generate_noise_temporal_task(stimfunction_tr, + motion_noise='rician', + ) + + # Test ARMA noise + with pytest.raises(ValueError): + noise_dict = {'fwhm': 4, 'auto_reg_rho': [1], 'ma_rho': [1, 1]} + sim._generate_noise_temporal_autoregression(stimfunction_tr, + noise_dict, + dimensions, + mask, + ) + + # Generate spatial noise + vol = sim._generate_noise_spatial(np.array([10, 10, 10, trs])) + assert len(vol.shape) == 3, 'Volume was not reshaped to ignore TRs' + + # Switch some of the noise types on + noise_dict = dict(physiological_sigma=1, drift_sigma=1, task_sigma=1, + auto_reg_sigma=0) + sim.generate_noise(dimensions=dimensions, + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict=noise_dict, + iterations=[0, 0], + ) def test_mask_brain(): @@ -244,7 +611,7 @@ def test_mask_brain(): ) # Mask the volume to be the same shape as a brain - mask, _ = sim.mask_brain(volume) + mask, _ = sim.mask_brain(dimensions, mask_self=None,) brain = volume * mask assert np.sum(brain != 0) == np.sum(volume != 0), "Masking did not work" @@ -262,11 +629,25 @@ def test_mask_brain(): ) # Mask the volume to be the same shape as a brain - mask, _ = sim.mask_brain(volume) + mask, _ = sim.mask_brain(dimensions, mask_self=None, ) brain = volume * mask assert np.sum(brain != 0) < np.sum(volume != 0), "Masking did not work" + # Test that you can load the default + dimensions = np.array([100, 100, 100]) + mask, template = sim.mask_brain(dimensions, mask_self=False) + + assert mask[20, 80, 50] == 0, 'Masking didn''t work' + assert mask[25, 80, 50] == 1, 'Masking didn''t work' + assert int(template[25, 80, 50] * 100) == 57, 'Template not correct' + + # Check that you can mask self + mask_self, template_self = sim.mask_brain(template, mask_self=True) + + assert (template_self - template).sum() < 1e2, 'Mask self error' + assert (mask_self - mask).sum() == 0, 'Mask self error' + def test_calc_noise(): @@ -275,27 +656,25 @@ def test_calc_noise(): event_durations = [6] tr_duration = 2 duration = 200 + temporal_res = 100 tr_number = int(np.floor(duration / tr_duration)) dimensions_tr = np.array([10, 10, 10, tr_number]) # Preset the noise dict - nd_orig = {'auto_reg_sigma': 0.6, - 'drift_sigma': 0.4, - 'snr': 30, - 'sfnr': 30, - 'max_activity': 1000, - 'fwhm': 4, - } + nd_orig = sim._noise_dict_update({}) # Create the time course for the signal to be generated stimfunction = sim.generate_stimfunction(onsets=onsets, event_durations=event_durations, total_time=duration, + temporal_resolution=temporal_res, ) # Mask the volume to be the same shape as a brain - mask, template = sim.mask_brain(dimensions_tr, mask_threshold=0.2) - stimfunction_tr = stimfunction[::int(tr_duration * 100)] + mask, template = sim.mask_brain(dimensions_tr, mask_self=None) + stimfunction_tr = stimfunction[::int(tr_duration * temporal_res)] + + nd_orig['matched'] = 0 noise = sim.generate_noise(dimensions=dimensions_tr[0:3], stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, @@ -304,26 +683,89 @@ def test_calc_noise(): noise_dict=nd_orig, ) - # Check that noise_system is being calculated correctly - spatial_sd = 5 - temporal_sd = 5 - noise_system = sim._generate_noise_system(dimensions_tr, - spatial_sd, - temporal_sd) + # Check the spatial noise match + nd_orig['matched'] = 1 + noise_matched = sim.generate_noise(dimensions=dimensions_tr[0:3], + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict=nd_orig, + iterations=[50, 0] + ) - precision = abs(noise_system[0, 0, 0, :].std() - spatial_sd) - assert precision < spatial_sd, 'noise_system calculated incorrectly' + # Calculate the noise parameters from this newly generated volume + nd_new = sim.calc_noise(noise, mask, template) + nd_matched = sim.calc_noise(noise_matched, mask, template) + + # Check the values are reasonable" + assert nd_new['snr'] > 0, 'snr out of range' + assert nd_new['sfnr'] > 0, 'sfnr out of range' + assert nd_new['auto_reg_rho'][0] > 0, 'ar out of range' + + # Check that the dilation increases SNR + no_dilation_snr = sim._calc_snr(noise_matched, + mask, + dilation=0, + reference_tr=tr_duration, + ) + + assert nd_new['snr'] > no_dilation_snr, "Dilation did not increase SNR" + + # Check that template size is in bounds + with pytest.raises(ValueError): + sim.calc_noise(noise, mask, template * 2) + + # Check that Mask is set is checked + with pytest.raises(ValueError): + sim.calc_noise(noise, None, template) + + # Check that it can deal with missing noise parameters + temp_nd = sim.calc_noise(noise, mask, template, noise_dict={}) + assert temp_nd['voxel_size'][0] == 1, 'Default voxel size not set' + + temp_nd = sim.calc_noise(noise, mask, template, noise_dict=None) + assert temp_nd['voxel_size'][0] == 1, 'Default voxel size not set' + + # Check that the fitting worked + snr_diff = abs(nd_orig['snr'] - nd_new['snr']) + snr_diff_match = abs(nd_orig['snr'] - nd_matched['snr']) + assert snr_diff > snr_diff_match, 'snr fit incorrectly' + + # Test that you can generate rician and exponential noise + sim._generate_noise_system(dimensions_tr, + 1, + 1, + spatial_noise_type='exponential', + temporal_noise_type='rician', + ) + + # Check the temporal noise match + nd_orig['matched'] = 1 + noise_matched = sim.generate_noise(dimensions=dimensions_tr[0:3], + stimfunction_tr=stimfunction_tr, + tr_duration=tr_duration, + template=template, + mask=mask, + noise_dict=nd_orig, + iterations=[0, 50] + ) - precision = abs(noise_system[:, :, :, 0].std() - temporal_sd) - assert precision < spatial_sd, 'noise_system calculated incorrectly' + nd_matched = sim.calc_noise(noise_matched, mask, template) - # Calculate the noise - nd_calc = sim.calc_noise(volume=noise, - mask=mask) + sfnr_diff = abs(nd_orig['sfnr'] - nd_new['sfnr']) + sfnr_diff_match = abs(nd_orig['sfnr'] - nd_matched['sfnr']) + assert sfnr_diff > sfnr_diff_match, 'sfnr fit incorrectly' - # How precise are these estimates - precision = abs(nd_calc['snr'] - nd_orig['snr']) - assert precision < nd_orig['snr'], 'snr calculated incorrectly' + ar1_diff = abs(nd_orig['auto_reg_rho'][0] - nd_new['auto_reg_rho'][0]) + ar1_diff_match = abs(nd_orig['auto_reg_rho'][0] - nd_matched[ + 'auto_reg_rho'][0]) + assert ar1_diff > ar1_diff_match, 'AR1 fit incorrectly' - precision = abs(nd_calc['sfnr'] - nd_orig['sfnr']) - assert precision < nd_orig['sfnr'], 'sfnr calculated incorrectly' + # Check that you can calculate ARMA for a single voxel + vox = noise[5, 5, 5, :] + arma = sim._calc_ARMA_noise(vox, + None, + sample_num=2, + ) + assert len(arma) == 2, "Two outputs not given by ARMA"