From d21b3862e1a065939a0056c2e7d491c8ac639cf3 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Sun, 10 Dec 2017 17:26:41 -0500 Subject: [PATCH 01/16] Changed the order of downsampling in convolve_hrf to accelerate computation. Changed drift to default to cosine functions --- brainiak/utils/fmrisim.py | 109 ++++++++++++++++++++++++++++++-------- 1 file changed, 87 insertions(+), 22 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 1e9d07833..32d8fd436 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -469,6 +469,16 @@ def generate_stimfunction(onsets, # Pull out the onsets, weights and durations, set as a float for line in text: onset, duration, weight = line.strip().split() + + # Check if the onset is more precise than the temporal resolution + dp = len(onsets[onsets.find('.') + 1:]) + + if 10 ** (dp - 1) > temporal_resolution: + raise ValueError('Temporal resolution is lower than the ' + 'decimal place precision of the timing ' + 'file. This can mean that events are ' + 'missed. Aborting') + onsets.append(float(onset)) event_durations.append(float(duration)) weights.append(float(weight)) @@ -771,7 +781,19 @@ def convolve_hrf(stimfunction, scale_function=True, temporal_resolution=1000.0, ): - """ Convolve the specified hrf with the timecourse + """ Convolve the specified hrf with the timecourse. In order + to accelerate the (slow) convolution, stimfunction is downsampled (mean + filtering) by the specified temporal resolution (e.g. a stimfunction of + 100000 elements becomes 100 elements if the temporal_resolution is + 1000). If temporal_resolution is 1 then the output will be the same + length as stimfunction. + Be aware that if scaling is on and events are very short + (> temporal_resolution * tr_duration) then the hrf may or may not come + out as anticipated. This is because very short events would evoke a small + absolute response after convolution but if there are only short events + and you scale then this will look identical to a convolution with longer + events. In general scaling is useful, which is why it is the default, + but be aware of this edge case Parameters ---------- @@ -804,10 +826,13 @@ def convolve_hrf(stimfunction, columns in this array. """ + # How will stimfunction be resized + stride = int(temporal_resolution * tr_duration) + duration = int(stimfunction.shape[0] / stride) # Generate the hrf to use in the convolution if hrf_type == 'double_gamma': - hrf = _double_gamma_hrf(temporal_resolution=temporal_resolution) + hrf = _double_gamma_hrf(temporal_resolution=1 / tr_duration) elif isinstance(hrf_type, list): hrf = hrf_type @@ -817,19 +842,21 @@ def convolve_hrf(stimfunction, # Create signal functions for each list in the stimfunction for list_counter in range(list_num): - # Take the stim function - stimfunction_temp = stimfunction[:, list_counter] + # Down sample the stim function so that it only has one element per + # TR. This accelerates the convolution greatly + stimfunction_temp = np.zeros((duration,)) + for sample in list(range(duration)): + idx_start = stride * sample + idx_end = stride * (sample + 1) + stimfunction_idx = stimfunction[idx_start : idx_end, list_counter] + stimfunction_temp[sample] = np.mean(stimfunction_idx) + # Perform the convolution signal_function_temp = np.convolve(stimfunction_temp, hrf) - # Decimate the signal function so that it only has one element per TR - decimate_interval = int(tr_duration * temporal_resolution) - signal_function_temp = signal_function_temp[0::decimate_interval] - - # Cut off the HRF - last_timepoint = stimfunction_temp.shape[0] / tr_duration - last_timepoint /= temporal_resolution - signal_function_temp = signal_function_temp[0:int(last_timepoint)] + # Shorten the output if the convolution made it grow + if len(signal_function_temp) > duration: + signal_function_temp = signal_function_temp[0:duration,] # Scale the function so that the peak response is 1 if scale_function: @@ -1403,13 +1430,14 @@ def _generate_noise_temporal_task(stimfunction_tr, def _generate_noise_temporal_drift(trs, tr_duration, - period=300, + basis="discrete_cos", + period=150, ): """Generate the drift noise - Create a sinewave, of a given period and random phase, to represent the - drift of the signal over time + Create a trend (either sine or discrete_cos), of a given period and random + phase, to represent the drift of the signal over time Parameters ---------- @@ -1420,6 +1448,11 @@ def _generate_noise_temporal_drift(trs, tr_duration : float How long in seconds is each volume acqusition + basis : str + What is the basis function for the drift. Could be made of discrete + cosines (number of bases scale with duration) or a sine wave with + the temporal order + period : int How many seconds is the period of oscillation of the drift @@ -1430,14 +1463,46 @@ def _generate_noise_temporal_drift(trs, """ - # Calculate the cycles of the drift for a given function. - cycles = trs * tr_duration / period + # Calculate drift differently depending on the basis function + if basis == 'discrete_cos': + + # Specify each tr in terms of its phase with the given period + timepoints = np.linspace(0, trs - 1, trs) + timepoints = ((timepoints * tr_duration) / period) * 2 * np.pi + + # Specify the other timing information + duration = trs * tr_duration + basis_funcs = int(np.floor(duration / period)) # How bases do you have + + if basis_funcs == 0: + logger.warning('Too few timepoints (' + str(trs) + ') to ' + 'accurately ' + 'model drift') + basis_funcs = 1 + + noise_drift = np.zeros((timepoints.shape[0], basis_funcs)) + for basis_counter in list(range(1, basis_funcs + 1)): + + # What steps do you want to take for this basis function + timepoints_basis = (timepoints/basis_counter) + (np.random.rand() + * np.pi * 2) + + # Store the drift from this basis func + noise_drift[:, basis_counter - 1]= np.cos(timepoints_basis) + + # Average the drift + noise_drift = np.mean(noise_drift, 1) + + elif basis == 'sine': + + # Calculate the cycles of the drift for a given function. + cycles = trs * tr_duration / period - # Create a sine wave with a given number of cycles and random phase - timepoints = np.linspace(0, trs - 1, trs) - phaseshift = np.pi * 2 * np.random.random() - phase = (timepoints / (trs - 1) * cycles * 2 * np.pi) + phaseshift - noise_drift = np.sin(phase) + # Create a sine wave with a given number of cycles and random phase + timepoints = np.linspace(0, trs - 1, trs) + phaseshift = np.pi * 2 * np.random.random() + phase = (timepoints / (trs - 1) * cycles * 2 * np.pi) + phaseshift + noise_drift = np.sin(phase) # Normalize so the sigma is 1 noise_drift = stats.zscore(noise_drift) From fcfe125629992b2fd0ee996cae9c1dcf74e7b19f Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Sun, 10 Dec 2017 17:32:38 -0500 Subject: [PATCH 02/16] Extended testing to protect against the HRF being miscalculated/convolved incorrectly --- tests/utils/test_fmrisim.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index 9aa3834f8..71d045aa5 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -88,7 +88,10 @@ def test_generate_stimfunction(): stim_dur = stimfunction.shape[0] / (tr_duration * 1000) assert signal_function.shape[0] == stim_dur, "The length did not change" - onsets = [10] + # Test + onsets = [0] + tr_duration = 1 + event_durations = [1] stimfunction = sim.generate_stimfunction(onsets=onsets, event_durations=event_durations, total_time=duration, @@ -97,6 +100,9 @@ def test_generate_stimfunction(): signal_function = sim.convolve_hrf(stimfunction=stimfunction, tr_duration=tr_duration, ) + + max_response = np.where(signal_function != 0)[0].max() + assert 25 < max_response < 30, "HRF is incorrect length" assert np.sum(signal_function < 0) > 0, "No values below zero" @@ -170,7 +176,7 @@ def test_generate_noise(): onsets = [10, 30, 50, 70, 90] event_durations = [6] tr_duration = 2 - duration = 100 + duration = 200 # Create the time course for the signal to be generated stimfunction = sim.generate_stimfunction(onsets=onsets, @@ -269,7 +275,7 @@ def test_calc_noise(): onsets = [10, 30, 50, 70, 90] event_durations = [6] tr_duration = 2 - duration = 100 + duration = 200 tr_number = int(np.floor(duration / tr_duration)) dimensions_tr = np.array([10, 10, 10, tr_number]) From 86044bfbb38a7ed052e42e6ee4643933059fd2a3 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Sun, 10 Dec 2017 21:09:40 -0500 Subject: [PATCH 03/16] PEP errors --- brainiak/utils/fmrisim.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 32d8fd436..2d88c631e 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -848,7 +848,7 @@ def convolve_hrf(stimfunction, for sample in list(range(duration)): idx_start = stride * sample idx_end = stride * (sample + 1) - stimfunction_idx = stimfunction[idx_start : idx_end, list_counter] + stimfunction_idx = stimfunction[idx_start:idx_end, list_counter] stimfunction_temp[sample] = np.mean(stimfunction_idx) # Perform the convolution @@ -856,7 +856,7 @@ def convolve_hrf(stimfunction, # Shorten the output if the convolution made it grow if len(signal_function_temp) > duration: - signal_function_temp = signal_function_temp[0:duration,] + signal_function_temp = signal_function_temp[0:duration, ] # Scale the function so that the peak response is 1 if scale_function: @@ -1475,9 +1475,9 @@ def _generate_noise_temporal_drift(trs, basis_funcs = int(np.floor(duration / period)) # How bases do you have if basis_funcs == 0: - logger.warning('Too few timepoints (' + str(trs) + ') to ' - 'accurately ' - 'model drift') + err_msg = 'Too few timepoints (' + str(trs) + ') to accurately ' \ + 'model drift' + logger.warning(err_msg) basis_funcs = 1 noise_drift = np.zeros((timepoints.shape[0], basis_funcs)) @@ -1488,7 +1488,7 @@ def _generate_noise_temporal_drift(trs, * np.pi * 2) # Store the drift from this basis func - noise_drift[:, basis_counter - 1]= np.cos(timepoints_basis) + noise_drift[:, basis_counter - 1] = np.cos(timepoints_basis) # Average the drift noise_drift = np.mean(noise_drift, 1) From cc778dac1d95cfa2537cbadbfb18745c12559b29 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Tue, 12 Dec 2017 11:26:49 -0500 Subject: [PATCH 04/16] Error --- brainiak/utils/fmrisim.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 2d88c631e..a0e6e44cd 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -471,7 +471,7 @@ def generate_stimfunction(onsets, onset, duration, weight = line.strip().split() # Check if the onset is more precise than the temporal resolution - dp = len(onsets[onsets.find('.') + 1:]) + dp = len(onset[onset.find('.') + 1:]) if 10 ** (dp - 1) > temporal_resolution: raise ValueError('Temporal resolution is lower than the ' From 6fabc57f414d4090baaf02aaf42fee0203223f70 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Thu, 21 Dec 2017 20:26:31 -0500 Subject: [PATCH 05/16] In response to Mingbo's review: now downsample after convolution (as before) in convolve_hrf, changed some defaults and updated the test. Also added some more plots to the example --- brainiak/utils/fmrisim.py | 76 +- .../utils/fmrisim_multivariate_example.ipynb | 1206 +++++++++++++++-- tests/utils/test_fmrisim.py | 8 +- 3 files changed, 1132 insertions(+), 158 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index a0e6e44cd..3eab03f06 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -408,7 +408,7 @@ def generate_stimfunction(onsets, total_time, weights=[1], timing_file=None, - temporal_resolution=1000.0, + temporal_resolution=100.0, ): """Return the function for the timecourse events @@ -471,9 +471,9 @@ def generate_stimfunction(onsets, onset, duration, weight = line.strip().split() # Check if the onset is more precise than the temporal resolution - dp = len(onset[onset.find('.') + 1:]) + upsampled_onset = float(onset) * temporal_resolution - if 10 ** (dp - 1) > temporal_resolution: + if upsampled_onset - np.round(upsampled_onset) != 0: raise ValueError('Temporal resolution is lower than the ' 'decimal place precision of the timing ' 'file. This can mean that events are ' @@ -521,7 +521,7 @@ def generate_stimfunction(onsets, def export_3_column(stimfunction, filename, - temporal_resolution=1000.0 + temporal_resolution=100.0 ): """ Output a tab separated three column timing file @@ -591,7 +591,7 @@ def export_3_column(stimfunction, def export_epoch_file(stimfunction, filename, tr_duration, - temporal_resolution=1000.0 + temporal_resolution=100.0 ): """ Output an epoch file, necessary for some inputs into brainiak @@ -700,7 +700,7 @@ def _double_gamma_hrf(response_delay=6, undershoot_dispersion=0.9, response_scale=1, undershoot_scale=0.035, - temporal_resolution=1000.0, + temporal_resolution=100.0, ): """Create the double gamma HRF with the timecourse evoked activity. Default values are based on Glover, 1999 and Walvaert, Durnez, @@ -779,21 +779,20 @@ def convolve_hrf(stimfunction, tr_duration, hrf_type='double_gamma', scale_function=True, - temporal_resolution=1000.0, + temporal_resolution=100.0, ): - """ Convolve the specified hrf with the timecourse. In order - to accelerate the (slow) convolution, stimfunction is downsampled (mean - filtering) by the specified temporal resolution (e.g. a stimfunction of - 100000 elements becomes 100 elements if the temporal_resolution is - 1000). If temporal_resolution is 1 then the output will be the same - length as stimfunction. - Be aware that if scaling is on and events are very short - (> temporal_resolution * tr_duration) then the hrf may or may not come - out as anticipated. This is because very short events would evoke a small - absolute response after convolution but if there are only short events - and you scale then this will look identical to a convolution with longer - events. In general scaling is useful, which is why it is the default, - but be aware of this edge case + """ Convolve the specified hrf with the timecourse. + The output of this is a downsampled convolution of the stimfunction and + the HRF function. If temporal_resolution is 1 / tr_duration then the + output will be the same length as stimfunction. + + Be aware that if scaling is on and event durations are less than the + duration of a TR then the hrf may or may not come out as anticipated. + This is because very short events would evoke a small absolute response + after convolution but if there are only short events and you scale then + this will look similar to a convolution with longer events. In general + scaling is useful, which is why it is the default, but be aware of this + edge case and if it is a concern, set the scale_function to false. Parameters ---------- @@ -832,7 +831,7 @@ def convolve_hrf(stimfunction, # Generate the hrf to use in the convolution if hrf_type == 'double_gamma': - hrf = _double_gamma_hrf(temporal_resolution=1 / tr_duration) + hrf = _double_gamma_hrf(temporal_resolution=temporal_resolution) elif isinstance(hrf_type, list): hrf = hrf_type @@ -842,32 +841,27 @@ def convolve_hrf(stimfunction, # Create signal functions for each list in the stimfunction for list_counter in range(list_num): + # Perform the convolution + signal_temp = np.convolve(stimfunction[:, list_counter], hrf) + # Down sample the stim function so that it only has one element per # TR. This accelerates the convolution greatly - stimfunction_temp = np.zeros((duration,)) + signal_vox = np.zeros((duration,)) for sample in list(range(duration)): idx_start = stride * sample idx_end = stride * (sample + 1) - stimfunction_idx = stimfunction[idx_start:idx_end, list_counter] - stimfunction_temp[sample] = np.mean(stimfunction_idx) - - # Perform the convolution - signal_function_temp = np.convolve(stimfunction_temp, hrf) - - # Shorten the output if the convolution made it grow - if len(signal_function_temp) > duration: - signal_function_temp = signal_function_temp[0:duration, ] + idxs = signal_temp[idx_start:idx_end] + signal_vox[sample] = np.mean(idxs) # Scale the function so that the peak response is 1 if scale_function: - signal_function_temp = signal_function_temp / np.max( - signal_function_temp) + signal_vox = signal_vox / np.max(signal_vox) # Add this function to the stack if list_counter == 0: - signal_function = np.zeros((len(signal_function_temp), list_num)) + signal_function = np.zeros((len(signal_vox), list_num)) - signal_function[:, list_counter] = signal_function_temp + signal_function[:, list_counter] = signal_vox return signal_function @@ -1414,11 +1408,11 @@ def _generate_noise_temporal_task(stimfunction_tr, # Make the noise to be added stimfunction_tr = stimfunction_tr != 0 if motion_noise == 'gaussian': - noise = stimfunction_tr * np.random.normal(0, 1, size=len( - stimfunction_tr)) + noise = stimfunction_tr * np.random.normal(0, 1, + size=stimfunction_tr.shape) elif motion_noise == 'rician': - noise = stimfunction_tr * stats.rice.rvs(0, 1, size=len( - stimfunction_tr)) + noise = stimfunction_tr * stats.rice.rvs(0, 1, + size=stimfunction_tr.shape) noise_task = stimfunction_tr + noise @@ -1450,8 +1444,8 @@ def _generate_noise_temporal_drift(trs, basis : str What is the basis function for the drift. Could be made of discrete - cosines (number of bases scale with duration) or a sine wave with - the temporal order + cosines (for longer run durations, more basis functions are + created) or a sine wave. period : int How many seconds is the period of oscillation of the drift diff --git a/examples/utils/fmrisim_multivariate_example.ipynb b/examples/utils/fmrisim_multivariate_example.ipynb index 02bece51e..e8115ef21 100644 --- a/examples/utils/fmrisim_multivariate_example.ipynb +++ b/examples/utils/fmrisim_multivariate_example.ipynb @@ -47,12 +47,20 @@ "*1.1 Import necessary Python packages*" ] }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.chdir('/Users/cellis/Documents/MATLAB/Analysis_BrainIAK')" + ] + }, { "cell_type": "code", "execution_count": 2, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook\n", @@ -177,7 +185,7 @@ "Noise parameters of the data were estimated as follows:\n", "SNR: 69.5116700015\n", "SFNR: 70.7171164885\n", - "FWHM: 5.65860178162\n" + "FWHM: 5.65977419548\n" ] } ], @@ -230,7 +238,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 59, "metadata": { "collapsed": true }, @@ -244,7 +252,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 60, "metadata": {}, "outputs": [ { @@ -1027,7 +1035,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -1039,10 +1047,10 @@ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 10, + "execution_count": 60, "metadata": {}, "output_type": "execute_result" } @@ -1071,7 +1079,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 61, "metadata": { "collapsed": true }, @@ -1079,9 +1087,10 @@ "source": [ "event_duration = 2\n", "isi = 3\n", - "total_time = int(dim[3] * tr)\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(0, events * (event_duration + isi), events) \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", @@ -1094,7 +1103,7 @@ "\n", "stimfunc_B = fmrisim.generate_stimfunction(onsets=onsets_B,\n", " event_durations=[event_duration],\n", - " total_time=int(dim[3] * tr),\n", + " total_time=total_time,\n", " temporal_resolution=temporal_res,\n", " )" ] @@ -1110,10 +1119,8 @@ }, { "cell_type": "code", - "execution_count": 12, - "metadata": { - "collapsed": true - }, + "execution_count": 62, + "metadata": {}, "outputs": [], "source": [ "fmrisim.export_epoch_file(stimfunction=[np.hstack((stimfunc_A, stimfunc_B))],\n", @@ -1144,7 +1151,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 63, "metadata": { "collapsed": true }, @@ -1167,7 +1174,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 64, "metadata": { "collapsed": true }, @@ -1182,7 +1189,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 65, "metadata": {}, "outputs": [ { @@ -1965,7 +1972,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -1977,10 +1984,10 @@ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 15, + "execution_count": 65, "metadata": {}, "output_type": "execute_result" } @@ -2025,7 +2032,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 57, "metadata": {}, "outputs": [ { @@ -2808,7 +2815,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -2823,7 +2830,7 @@ "(-0.5, 63.5, 63.5, -0.5)" ] }, - "execution_count": 17, + "execution_count": 57, "metadata": {}, "output_type": "execute_result" } @@ -2884,7 +2891,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 56, "metadata": {}, "outputs": [ { @@ -3667,7 +3674,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -3682,7 +3689,7 @@ "(-0.5, 63.5, 63.5, -0.5)" ] }, - "execution_count": 20, + "execution_count": 56, "metadata": {}, "output_type": "execute_result" } @@ -3699,74 +3706,12 @@ "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. Although other simulators use a combination of discrete cosines for estimating drift (Welvaert, et al., 2011), our testing suggests that this provides poor estimates in long scans (>200s). 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. \n", - "\t\n", - "*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\n", - "\n", - "*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. \n", - "\n", - "*3.4 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. \n" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "brain = signal + noise" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*4.1 Pull out data for each trial*\n", - "\n", - "Identify which voxels are in the signal ROI by using the coordinates provided earlier. To identify the relevant timepoints, assume that the peak of the neural response occurs 4 - 6s after each event onset. Take the TR corresponding to this peak response as the TR for that trial. In longer event/block designs you might instead average over each event." - ] - }, - { - "cell_type": "code", - "execution_count": 22, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "hrf_lag = 4 # Assumed time from stimulus onset to HRF peak\n", - "\n", - "# Get the lower and upper bounds of the ROI\n", - "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", - "\n", - "trials_A = trials_A.reshape((voxels, trials_A.shape[3]))\n", - "trials_B = trials_B.reshape((voxels, trials_B.shape[3]))" + "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. Although other simulators use a combination of discrete cosines for estimating drift (Welvaert, et al., 2011), our testing suggests that this provides poor estimates in long scans (>200s). 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": 23, + "execution_count": 31, "metadata": {}, "outputs": [ { @@ -4549,7 +4494,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -4561,37 +4506,38 @@ { "data": { "text/plain": [ - "" + "(-0.5, 63.5, 63.5, -0.5)" ] }, - "execution_count": 23, + "execution_count": 31, "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.imshow(trials_A)\n", - "plt.ylabel('Voxels')\n", - "plt.xlabel('Trials')\n", - "plt.subplot(1,2,2)\n", - "plt.imshow(trials_B)\n", - "plt.xlabel('Trials')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*4.2 Represent the data*\n", + "plt.title('Low noise')\n", + "plt.imshow(low_spatial[:, :, 12])\n", + "plt.axis('off')\n", "\n", - "Treat each voxel as a dimension and each trial as a point in this voxel space. It is then possible to display the different conditions and determine whether these are separable in this lower dimensionality (note that the conditions may be separable in higher dimensionality but unsupervised techniques like Multidimensional Scaling used below, might not show such a difference)" + "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": 24, + "execution_count": 30, "metadata": {}, "outputs": [ { @@ -5374,7 +5320,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -5383,6 +5329,1040 @@ "metadata": {}, "output_type": "display_data" }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "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", + "\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" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "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", + " this.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", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\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": [ - "" - ] - }, - "execution_count": 60, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "plt.figure()\n", "plt.subplot(1,2,1)\n", @@ -1079,7 +275,7 @@ }, { "cell_type": "code", - "execution_count": 61, + "execution_count": null, "metadata": { "collapsed": true }, @@ -1119,8 +315,10 @@ }, { "cell_type": "code", - "execution_count": 62, - "metadata": {}, + "execution_count": null, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "fmrisim.export_epoch_file(stimfunction=[np.hstack((stimfunc_A, stimfunc_B))],\n", @@ -1151,7 +349,7 @@ }, { "cell_type": "code", - "execution_count": 63, + "execution_count": null, "metadata": { "collapsed": true }, @@ -1174,7 +372,7 @@ }, { "cell_type": "code", - "execution_count": 64, + "execution_count": null, "metadata": { "collapsed": true }, @@ -1189,809 +387,11 @@ }, { "cell_type": "code", - "execution_count": 65, - "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", - " this.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", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\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": 57, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "plt.figure()\n", "plt.imshow(signal_volume[:, :, 24], cmap=plt.cm.gray)\n", @@ -2853,7 +455,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": null, "metadata": { "collapsed": true }, @@ -2874,7 +476,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": null, "metadata": { "collapsed": true }, @@ -2891,809 +493,11 @@ }, { "cell_type": "code", - "execution_count": 56, - "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", - " this.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", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\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": 31, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "# Plot spatial noise\n", "low_spatial = fmrisim._generate_noise_spatial(dim[0:3],\n", @@ -4537,809 +543,11 @@ }, { "cell_type": "code", - "execution_count": 30, - "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", - " this.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", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\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" - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "system = fmrisim._generate_noise_system(dimensions_tr=dim,\n", " spatial_sd=1.5,\n", @@ -6232,7 +652,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": null, "metadata": { "collapsed": true }, @@ -6261,7 +681,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, "metadata": { "collapsed": true }, @@ -6282,41 +702,11 @@ }, { "cell_type": "code", - "execution_count": 23, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 23, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "plt.figure()\n", "plt.subplot(1,2,1)\n", @@ -6339,41 +729,11 @@ }, { "cell_type": "code", - "execution_count": 24, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 24, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "distance_matrix = sp_distance.squareform(sp_distance.pdist(np.vstack([trials_A.transpose(), trials_B.transpose()])))\n", "\n", @@ -6396,18 +756,11 @@ }, { "cell_type": "code", - "execution_count": 25, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Mean difference between condition A and B: -0.52\n", - "pvalue: 0.677\n" - ] - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "mean_difference = (np.mean(trials_A,0) - np.mean(trials_B,0))\n", "ttest = stats.ttest_1samp(mean_difference, 0)\n", @@ -6427,17 +780,11 @@ }, { "cell_type": "code", - "execution_count": 26, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Classification accuracy between condition A and B: 0.833\n" - ] - } - ], + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], "source": [ "input_mat = np.vstack([trials_A.transpose(), trials_B.transpose()])\n", "input_labels = trials_A.shape[1] * [1] + trials_B.shape[1] * [0]\n", diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index 0d6bf0f00..a1b0f8278 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -101,7 +101,7 @@ def test_generate_stimfunction(): ) max_response = np.where(signal_function != 0)[0].max() - assert 25 < max_response <= 30, "HRF is incorrect length" + assert 25 < max_response <= 30, "HRF has the incorrect length" assert np.sum(signal_function < 0) > 0, "No values below zero" From b00aa129a4ee163b6611bb6c4fa30ae30a75f86c Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Mon, 8 Jan 2018 13:52:41 -0500 Subject: [PATCH 09/16] Changed raise to warning for temporal precision error --- brainiak/utils/fmrisim.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index a568b710c..acd2f5af3 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -473,13 +473,21 @@ def generate_stimfunction(onsets, # Check if the onset is more precise than the temporal resolution upsampled_onset = float(onset) * temporal_resolution - # Because of float precision, there can be issues. E.g. - # float('1.001') * 1000 = 1000.99 + # Because of float precision, the upsampled values might + # not round as expected . + # E.g. float('1.001') * 1000 = 1000.99 if np.allclose(upsampled_onset, np.round(upsampled_onset)): - raise ValueError('Temporal resolution is lower than the ' - 'decimal place precision of the timing ' - 'file. This can mean that events are ' - 'missed. Aborting') + warning = 'Your onset: ' + str(onset) + ' has more decimal ' \ + 'points than the ' \ + 'specified temporal ' \ + 'resolution can ' \ + 'resolve. This means' \ + ' that events might' \ + ' be missed. ' \ + 'Consider increasing' \ + ' the temporal ' \ + 'resolution.' + logging.warning(warning) onsets.append(float(onset)) event_durations.append(float(duration)) From c810e9779f053563cc2a8dc2da9541c0084124d8 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Mon, 8 Jan 2018 20:49:26 -0500 Subject: [PATCH 10/16] Made the evoked activity proportional to the mean activity of the voxel --- .../utils/fmrisim_multivariate_example.ipynb | 168 ++++++++---------- 1 file changed, 71 insertions(+), 97 deletions(-) diff --git a/examples/utils/fmrisim_multivariate_example.ipynb b/examples/utils/fmrisim_multivariate_example.ipynb index ada8d7f1c..e4ca5d68f 100644 --- a/examples/utils/fmrisim_multivariate_example.ipynb +++ b/examples/utils/fmrisim_multivariate_example.ipynb @@ -112,18 +112,8 @@ "dimsize = nii.header.get_zooms()\n", "tr = dimsize[3]\n", "if tr > 100: # If high then these values are likely in ms\n", - " tr /= 1000" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "dim" + " tr /= 1000\n", + "print(dim)" ] }, { @@ -179,9 +169,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "print('Noise parameters of the data were estimated as follows:')\n", @@ -203,9 +191,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.1 Establish effect size*\n", + "*2.1 Specify which voxels in the brain contain signal*\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. 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)." + "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 ROI analyses). The value of each voxel can be specified here, or set to be a random value." ] }, { @@ -216,18 +204,35 @@ }, "outputs": [], "source": [ - "effect_size = 1\n", - "temporal_sd = (template[mask > 0]).mean() * noise_dict['max_activity'] / noise_dict['sfnr']\n", - "effect_signal_change = effect_size * temporal_sd" + "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", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.imshow(signal_volume[:, :, 24], cmap=plt.cm.gray)\n", + "plt.imshow(mask[:, :, 24], cmap=plt.cm.gray, alpha=.5)\n", + "plt.axis('off')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*2.2 Characterize signal for voxels*\n", + "*2.2 Establish effect size*\n", "\n", - "Specify the pattern of activity across a given number of voxels that characterizes each condition. This pattern can simply be random, as is done here, or can be structured, like the position of voxels in high dimensional representation space." + "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)." ] }, { @@ -238,10 +243,26 @@ }, "outputs": [], "source": [ - "feature_size = 3\n", - "voxels = feature_size ** 3\n", - "pattern_A = np.random.randn(voxels).reshape((voxels, 1)) * effect_signal_change\n", - "pattern_B = np.random.randn(voxels).reshape((voxels, 1)) * effect_signal_change" + "effect_size = 1\n", + "temporal_sd = template * noise_dict['max_activity'] / noise_dict['sfnr']\n", + "\n", + "signal_idxs = np.where(signal_volume == 1)\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]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*2.3 Characterize signal for voxels*\n", + "\n", + "Specify the pattern of activity across a given number of voxels that characterizes each condition. This pattern can simply be random, as is done here, or can be structured, like the position of voxels in high dimensional representation space." ] }, { @@ -251,6 +272,17 @@ "collapsed": true }, "outputs": [], + "source": [ + "voxels = feature_size ** 3\n", + "pattern_A = np.random.randn(voxels).reshape((voxels, 1)) * signal_change\n", + "pattern_B = np.random.randn(voxels).reshape((voxels, 1)) * signal_change" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "plt.figure()\n", "plt.subplot(1,2,1)\n", @@ -268,7 +300,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.3 Generate event time course*\n", + "*2.4 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" ] @@ -308,7 +340,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.4 Export stimulus time course for analysis*\n", + "*2.5 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" ] @@ -342,7 +374,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.5 Estimate the voxel weight for each event*\n", + "*2.6 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." ] @@ -365,7 +397,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*2.6 Convolve each voxel’s time course with the Hemodynamic Response Function*\n", + "*2.7 Convolve each voxel’s time course with the Hemodynamic Response Function*\n", "\n", "With the time course of stimulus events it is necessary to estimate the brain’s response to those events, which can be estimated by convolving it with using a Hemodynamic Response Function (HRF). By default, convolve_hrf assumes a double gamma HRF appropriately models a brain’s response to events, as modeled by fMRI (Friston, et al., 1998). To do this, each voxel’s time course is convolved to make a function of the signal activity. Hence this produces an estimate of the voxel’s activity, after considering the temporal blurring of the HRF. This can take a single vector of events or multiple time courses." ] @@ -388,9 +420,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# Display signal\n", @@ -404,46 +434,6 @@ "plt.legend(loc=1)" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*2.7 Specify which voxels in the brain contain signal*\n", - "\n", - "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 ROI analyses). The value of each voxel can be specified here, or set to be a random value." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "coordinates = np.array([[24, 24, 24]])\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": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "plt.figure()\n", - "plt.imshow(signal_volume[:, :, 24], cmap=plt.cm.gray)\n", - "plt.imshow(mask[:, :, 24], cmap=plt.cm.gray, alpha=.5)\n", - "plt.axis('off')" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -494,9 +484,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", @@ -516,9 +504,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# Plot spatial noise\n", @@ -544,9 +530,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# Create the different types of noise\n", @@ -602,9 +586,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "system = fmrisim._generate_noise_system(dimensions_tr=dim,\n", @@ -703,9 +685,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", @@ -730,9 +710,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "distance_matrix = sp_distance.squareform(sp_distance.pdist(np.vstack([trials_A.transpose(), trials_B.transpose()])))\n", @@ -757,9 +735,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "mean_difference = (np.mean(trials_A,0) - np.mean(trials_B,0))\n", @@ -781,9 +757,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "input_mat = np.vstack([trials_A.transpose(), trials_B.transpose()])\n", From ad70aa1a70fb756de68049bc2b05bf8736832af2 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Tue, 9 Jan 2018 21:22:19 -0500 Subject: [PATCH 11/16] Downsample the HRF rather average, as we are doing slice time correction --- brainiak/utils/fmrisim.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index acd2f5af3..e6d7f5a61 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -794,7 +794,9 @@ def convolve_hrf(stimfunction, """ Convolve the specified hrf with the timecourse. The output of this is a downsampled convolution of the stimfunction and the HRF function. If temporal_resolution is 1 / tr_duration then the - output will be the same length as stimfunction. + output will be the same length as stimfunction. This time course assumes + that slice time correction has occurred and all slices have been aligned + to the middle time point in the TR. Be aware that if scaling is on and event durations are less than the duration of a TR then the hrf may or may not come out as anticipated. @@ -854,10 +856,12 @@ def convolve_hrf(stimfunction, # Perform the convolution signal_temp = np.convolve(stimfunction[:, list_counter], hrf) - # Down sample the stim function so that it only has one element per - # TR. This accelerates the convolution greatly + # Down sample the signal function so that it only has one element per + # TR. This assumes that all slices are collected at the same time, + # which is often the result of slice time correction. In other + # words, the output assumes slice time correction signal_temp = signal_temp[:duration * stride] - signal_vox = np.mean(signal_temp.reshape(-1, stride), 1) + signal_vox = signal_temp[int(stride / 2)::stride] # Scale the function so that the peak response is 1 if scale_function: From 38b2fdb18ec08e5be12efeb0b4fee1aae5838df8 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 10 Jan 2018 00:36:08 -0500 Subject: [PATCH 12/16] Updated gen_design to solve an issue resulting from changes to convolve_hrf --- brainiak/utils/utils.py | 2 ++ tests/utils/test_utils.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/brainiak/utils/utils.py b/brainiak/utils/utils.py index 4b60497e7..33a0a29dd 100644 --- a/brainiak/utils/utils.py +++ b/brainiak/utils/utils.py @@ -379,6 +379,8 @@ def gen_design(stimtime_files, scan_duration, TR, style='FSL', It is acceptable to not provide the weight, or not provide both duration and weight. In such cases, these parameters will default to 1.0. + This code will accept timing files with only 1 or 2 columns for + convenience but please note that the FSL package does not allow this 'AFNI' style has one line for each scan (run). Each line has a few triplets in the format of diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index f58a2768c..2decfd6b4 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -123,7 +123,7 @@ def test_gen_design(): 'gen_design does not treat missing values correctly') design5 = gen_design(stimtime_files=[files['FSL2']], scan_duration=[48, 20], TR=1) - assert np.all(np.isclose(design4, design5[::2])), ( + assert np.allclose(design4, design5[::2], rtol=0.01), ( 'design matrices sampled at different frequency do not match' ' at corresponding time points') design6 = gen_design(stimtime_files=[files['AFNI1']], From 63795d2226ded4b498039d2611621a4e81f5f108 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 10 Jan 2018 00:47:22 -0500 Subject: [PATCH 13/16] Updated for test --- tests/utils/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index 2decfd6b4..820f8ae79 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -123,7 +123,7 @@ def test_gen_design(): 'gen_design does not treat missing values correctly') design5 = gen_design(stimtime_files=[files['FSL2']], scan_duration=[48, 20], TR=1) - assert np.allclose(design4, design5[::2], rtol=0.01), ( + assert np.allclose(design4, design5[::2], rtol=0.1), ( 'design matrices sampled at different frequency do not match' ' at corresponding time points') design6 = gen_design(stimtime_files=[files['AFNI1']], From 56d949b76f690efaefadad45b2b038eee699c149 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 10 Jan 2018 09:59:30 -0500 Subject: [PATCH 14/16] New changes to convolve hrf add a bit more slop to some tests for gen_design which need to be relaxed --- tests/utils/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index 820f8ae79..385a22be5 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -123,7 +123,7 @@ def test_gen_design(): 'gen_design does not treat missing values correctly') design5 = gen_design(stimtime_files=[files['FSL2']], scan_duration=[48, 20], TR=1) - assert np.allclose(design4, design5[::2], rtol=0.1), ( + assert (design4 - design5[::2]).mean() < 0.1, ( 'design matrices sampled at different frequency do not match' ' at corresponding time points') design6 = gen_design(stimtime_files=[files['AFNI1']], From 435a5d998839ab6317cb19e8c39cb8b8ce77f87f Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 10 Jan 2018 11:41:23 -0500 Subject: [PATCH 15/16] Error in the directionality of the warning for generate_stimfunction --- brainiak/utils/fmrisim.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index e6d7f5a61..69cd91de8 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -476,7 +476,7 @@ def generate_stimfunction(onsets, # Because of float precision, the upsampled values might # not round as expected . # E.g. float('1.001') * 1000 = 1000.99 - if np.allclose(upsampled_onset, np.round(upsampled_onset)): + if np.allclose(upsampled_onset, np.round(upsampled_onset)) == 0: warning = 'Your onset: ' + str(onset) + ' has more decimal ' \ 'points than the ' \ 'specified temporal ' \ From 45e7b932164de64c677162865e79077f6865b1ad Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Fri, 12 Jan 2018 14:42:50 -0500 Subject: [PATCH 16/16] Updated test --- tests/utils/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index 385a22be5..fb8f7e24c 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -123,7 +123,7 @@ def test_gen_design(): 'gen_design does not treat missing values correctly') design5 = gen_design(stimtime_files=[files['FSL2']], scan_duration=[48, 20], TR=1) - assert (design4 - design5[::2]).mean() < 0.1, ( + assert (np.abs(design4 - design5[::2])).mean() < 0.1, ( 'design matrices sampled at different frequency do not match' ' at corresponding time points') design6 = gen_design(stimtime_files=[files['AFNI1']],