diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 1e9d07833..69cd91de8 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 @@ -469,6 +469,26 @@ 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 + upsampled_onset = float(onset) * temporal_resolution + + # 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)) == 0: + 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)) weights.append(float(weight)) @@ -511,7 +531,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 @@ -581,7 +601,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 @@ -690,7 +710,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, @@ -769,9 +789,22 @@ 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 + """ 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. 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. + 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 ---------- @@ -804,6 +837,9 @@ 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': @@ -817,30 +853,25 @@ 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] - - 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] + # Perform the convolution + signal_temp = np.convolve(stimfunction[:, list_counter], hrf) - # 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)] + # 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 = signal_temp[int(stride / 2)::stride] # 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 @@ -1387,11 +1418,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 @@ -1403,13 +1434,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 +1452,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 (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 @@ -1430,14 +1467,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: + 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)) + 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) 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/examples/utils/fmrisim_multivariate_example.ipynb b/examples/utils/fmrisim_multivariate_example.ipynb index 02bece51e..e4ca5d68f 100644 --- a/examples/utils/fmrisim_multivariate_example.ipynb +++ b/examples/utils/fmrisim_multivariate_example.ipynb @@ -49,7 +49,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": { "collapsed": true }, @@ -80,7 +80,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": { "collapsed": true }, @@ -102,7 +102,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": { "collapsed": true }, @@ -112,7 +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" + " tr /= 1000\n", + "print(dim)" ] }, { @@ -126,7 +127,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": { "collapsed": true }, @@ -152,7 +153,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": { "collapsed": true }, @@ -167,20 +168,9 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Noise parameters of the data were estimated as follows:\n", - "SNR: 69.5116700015\n", - "SFNR: 70.7171164885\n", - "FWHM: 5.65860178162\n" - ] - } - ], + "outputs": [], "source": [ "print('Noise parameters of the data were estimated as follows:')\n", "print('SNR: ' + str(noise_dict['snr']))\n", @@ -201,852 +191,98 @@ "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." + "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": 8, + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "coordinates = np.array([[24, 24, 24]])\n", + "feature_size = 3\n", + "signal_volume = fmrisim.generate_signal(dimensions=dim[0:3],\n", + " feature_type=['cube'],\n", + " feature_coordinates=coordinates,\n", + " feature_size=[feature_size],\n", + " signal_magnitude=[1],\n", + " )" + ] + }, + { + "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 Establish effect size*\n", + "\n", + "When specifying the signal we must determine the amount of activity change each voxel undergoes. A useful metric for this is the SFNR value determined from noise calculations because it can be used to estimate the variability in the average voxel. This is set up so that the evoked activity is proportional to the mean activity of a voxel. For a univariate effect, to estimate activity with a Cohen’s d of 1, the size of the change must be equivalent to one standard deviation. For multivariate effects the effect size depends on multiple factors including the number of voxels and conditions. Different measures for effect size could also be calculated, such as percent signal change. Note that that this signal change is based on the average voxel. Instead it might be preferable to model signal change based on the mean of each voxel (i.e. the template value)." + ] + }, + { + "cell_type": "code", + "execution_count": null, "metadata": { "collapsed": true }, "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" + "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.2 Characterize signal for voxels*\n", + "*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." ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": { "collapsed": true }, "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" + "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": 10, + "execution_count": null, "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": 15, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "# Display signal\n", "plt.figure()\n", @@ -1997,844 +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": 16, - "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": 17, - "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": 20, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "plt.figure()\n", "plt.imshow(noise[:, :, 24, 0], cmap=plt.cm.gray)\n", @@ -3699,24 +498,143 @@ "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", + "The temporal noise of fMRI data is comprised of multiple components: drift, autoregression, task related motion and physiological noise. To estimate drift, a sine wave with a default period of 300s, is used. To estimate drift, cosine basis functions are combined, with longer runs being comprised of more basis functions (Welvaert, et al., 2011). This drift is then multiplied by a three-dimensional volume of Gaussian random fields of a specific FWHM. Autoregression noise is estimated by creating a time course of Gaussian noise values that are weighted by previous values of the time course. This autoregressive time course is multiplied by a brain shaped volume of Gaussian random fields. Physiological noise is modeled by sine waves comprised of heart rate (1.17Hz) and respiration rate (0.2Hz) (Biswal, et al., 1996) with random phase. This time course is also multiplied by brain shaped spatial noise. Finally, task related noise is simulated by adding Gaussian or Rician noise to time points where there are events (according to the event time course) and in turn this is multiplied by a brain shaped spatial noise volume. These four noise components are then mixed together in proportion to the size of their corresponding noise values. This aggregated volume is then Z scored and the SFNR is used to estimate the appropriate standard deviation of these values across time. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot spatial noise\n", + "low_spatial = fmrisim._generate_noise_spatial(dim[0:3],\n", + " fwhm=4.0,\n", + " )\n", + "\n", + "high_spatial = fmrisim._generate_noise_spatial(dim[0:3],\n", + " fwhm=1.0,\n", + " )\n", + "plt.figure()\n", + "plt.subplot(1,2,1)\n", + "plt.title('Low noise')\n", + "plt.imshow(low_spatial[:, :, 12])\n", + "plt.axis('off')\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.title('High noise')\n", + "plt.imshow(high_spatial[:, :, 12])\n", + "plt.axis('off')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create the different types of noise\n", + "timepoints = list(range(0, total_time, int(tr)))\n", + "\n", + "drift = fmrisim._generate_noise_temporal_drift(total_time,\n", + " int(tr),\n", + " )\n", + "\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\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": null, + "metadata": {}, + "outputs": [], + "source": [ + "system = fmrisim._generate_noise_system(dimensions_tr=dim,\n", + " spatial_sd=1.5,\n", + " temporal_sd=1,\n", + " )\n", "\n", + "plt.figure()\n", + "plt.subplot(1, 3, 1)\n", + "plt.hist(system.flatten())\n", + "plt.title('Activity distribution')\n", + "plt.xlabel('Activity')\n", + "plt.ylabel('Frequency')\n", + "\n", + "plt.subplot(1, 3, 2)\n", + "plt.imshow(system[:, :, 0, 0])\n", + "plt.axis('off')\n", + "plt.title('Spatial plane')\n", + "plt.clim([system.min(), np.percentile(system, 95)])\n", + "\n", + "plt.subplot(1, 3, 3)\n", + "plt.imshow(system[0, :64, 0, :64])\n", + "plt.axis('off')\n", + "plt.title('Temporal plane')\n", + "plt.clim([system.min(), np.percentile(system, 95)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "*3.3 Combine noise and template*\n", " \n", - "The template volume is used to estimate the appropriate baseline distribution of MR values. This estimate is then combined with the temporal noise and the system noise to make an estimate of the noise. \n", - "\n", + "The template volume is used to estimate the appropriate baseline distribution of MR values. This estimate is then combined with the temporal noise and the system noise to make an estimate of the noise. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "*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" + "Since the brain signal is expected to be small and sparse relative to the noise, it is assumed sufficient to simply add the volume containing signal with the volume modeling noise to make the simulated brain. " ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": null, "metadata": { "collapsed": true }, @@ -3745,7 +663,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, "metadata": { "collapsed": true }, @@ -3766,809 +684,9 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": null, "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": 24, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "distance_matrix = sp_distance.squareform(sp_distance.pdist(np.vstack([trials_A.transpose(), trials_B.transpose()])))\n", "\n", @@ -5416,18 +734,9 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Mean difference between condition A and B: -0.52\n", - "pvalue: 0.677\n" - ] - } - ], + "outputs": [], "source": [ "mean_difference = (np.mean(trials_A,0) - np.mean(trials_B,0))\n", "ttest = stats.ttest_1samp(mean_difference, 0)\n", @@ -5447,17 +756,9 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Classification accuracy between condition A and B: 0.833\n" - ] - } - ], + "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", @@ -5521,5 +822,5 @@ } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 1 } diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index 9aa3834f8..a1b0f8278 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -75,9 +75,8 @@ def test_generate_stimfunction(): total_time=duration, ) - assert stimfunction.shape[0] == duration * 1000, "stimfunc incorrect " \ - "length" - eventNumber = np.sum(event_durations * len(onsets)) * 1000 + assert stimfunction.shape[0] == duration * 100, "stimfunc incorrect length" + eventNumber = np.sum(event_durations * len(onsets)) * 100 assert np.sum(stimfunction) == eventNumber, "Event number" # Create the signal function @@ -85,10 +84,13 @@ def test_generate_stimfunction(): tr_duration=tr_duration, ) - stim_dur = stimfunction.shape[0] / (tr_duration * 1000) + stim_dur = stimfunction.shape[0] / (tr_duration * 100) 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 +99,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 has the incorrect length" assert np.sum(signal_function < 0) > 0, "No values below zero" @@ -170,7 +175,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, @@ -193,7 +198,7 @@ def test_generate_noise(): assert min(mask[mask > 0]) > 0.1, "Mask thresholding did not work" assert len(np.unique(template) > 2), "Template creation did not work" - stimfunction_tr = stimfunction[::int(tr_duration * 1000)] + stimfunction_tr = stimfunction[::int(tr_duration * 100)] # Create the noise volumes (using the default parameters) noise = sim.generate_noise(dimensions=dimensions, stimfunction_tr=stimfunction_tr, @@ -269,7 +274,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]) @@ -290,7 +295,7 @@ def test_calc_noise(): # Mask the volume to be the same shape as a brain mask, template = sim.mask_brain(dimensions_tr, mask_threshold=0.2) - stimfunction_tr = stimfunction[::int(tr_duration * 1000)] + stimfunction_tr = stimfunction[::int(tr_duration * 100)] noise = sim.generate_noise(dimensions=dimensions_tr[0:3], stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index f58a2768c..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 np.all(np.isclose(design4, design5[::2])), ( + 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']],