Merge pull request #9 from gallantlab/doc/improveflow

DOC improve flow with ngrok and cortex, add example delays

Author: mvdoc

tutorials/movies_3T/00_setup_colab.py

@@ -8,7 +8,6 @@
 skip it if you run the tutorials on your machine.
 """
 # sphinx_gallery_thumbnail_path = "static/colab.png"
-#
 ###############################################################################
 # Change runtime to use a GPU
 # ---------------------------

tutorials/movies_3T/01_plot_explainable_variance.py

@@ -37,8 +37,6 @@
 references [1]_ [2]_ [3]_.
 """
 # sphinx_gallery_thumbnail_number = 1
-
-
 ###############################################################################
 # Path of the data directory
 # --------------------------
@@ -172,19 +170,11 @@
 #
 # The second mapper we provide maps the voxel data to a Freesurfer
 # average surface ("fsaverage"), that can be used in ``pycortex``.
-# First, let's download the "fsaverage" surface.
-
-import cortex
-
-surface = "fsaverage"
-
-if not hasattr(cortex.db, surface):
-    cortex.utils.download_subject(subject_id=surface)
-
-###############################################################################
+#
 # If you are running the notebook on Colab, you might need to update the
 # pycortex filestore as following:
 
+import cortex
 try:
     import google.colab  # noqa
     in_colab = True
@@ -201,11 +191,20 @@
     cortex.quickflat.utils.db = cortex.database.db
     cortex.quickflat.composite.db = cortex.database.db
 
+###############################################################################
+# Now, let's download the "fsaverage" surface.
+
+surface = "fsaverage"
+
+if not hasattr(cortex.db, surface):
+    cortex.utils.download_subject(subject_id=surface)
+
 ###############################################################################
 # Then, we load the "fsaverage" mapper. The mapper is a matrix of shape
 # (n_vertices, n_voxels), which maps each voxel to some vertices in the
 # fsaverage surface. It is stored as a sparse CSR matrix. The mapper is applied
 # with a dot product ``@`` (equivalent to ``np.dot``).
+
 from voxelwise_tutorials.io import load_hdf5_sparse_array
 voxel_to_fsaverage = load_hdf5_sparse_array(mapper_file,
                                             key='voxel_to_fsaverage')
@@ -220,14 +219,9 @@
 vertex = cortex.Vertex(ev_projected, surface, vmin=0, vmax=0.7, cmap='viridis')
 
 ###############################################################################
-# To start an interactive 3D viewer in the browser, use the ``webshow``
-# function.
-
-if False:
-    cortex.webshow(vertex, open_browser=False, port=8050)
-
-###############################################################################
-# If you are running the notebook on Colab, you need to tunnel the pycortex
+# To start an interactive 3D viewer in the browser, we can use the ``webshow``
+# function in pycortex.
+# If you are running the notebook on Colab, you first need to tunnel the pycortex
 # application out of Colab. To do so, use the following cell to start a tunnel
 # with ``ngrok`` and to get an address where the pycortex viewer will be made
 # accessible.
@@ -246,6 +240,17 @@
           f"{result}\n"
           "and not the one proposed by pycortex ('Open viewer: ...')\n")
 
+
+###############################################################################
+# Now you can start an interactive 3D viewer by changing ``run_webshow`` to
+# ``True`` and running the following cell.  If you are using Colab, remember to
+# use the address returned by ngrok in the cell above rather than the address
+# returned by this cell.
+
+run_webshow = False
+if run_webshow:
+    cortex.webshow(vertex, open_browser=False, port=8050)
+
 ###############################################################################
 # Alternatively, to plot a flatmap in a ``matplotlib`` figure, use the
 # `quickshow` function.
@@ -255,7 +260,6 @@
 
 from cortex.testing_utils import has_installed
 
-
 fig = cortex.quickshow(vertex, colorbar_location='right',
                        with_rois=has_installed("inkscape"))
 plt.show()

tutorials/movies_3T/03_plot_hemodynamic_response.py

@@ -141,6 +141,64 @@
 scores = backend.to_numpy(scores)
 print("(n_voxels,) =", scores.shape)
 
+###############################################################################
+# Intermission: understanding delays
+# ----------------------------------
+#
+# To have an intuitive understanding of what we accomplish by delaying the
+# features before model fitting, we will simulate one voxel and a single
+# feature. We will then create a ``Delayer`` object (which was used in the
+# previous pipeline) and visualize its effect on our single feature. Let's
+# start by simulating the data.
+
+# number of total trs
+n_trs = 50
+# repetition time for the simulated data
+TR = 2.0
+rng = np.random.RandomState(42)
+y = rng.randn(n_trs)
+x = np.zeros(n_trs)
+# add some arbitrary value to our feature
+x[15:20] = .5
+x += rng.randn(n_trs) * 0.1  # add some noise
+
+# create a delayer object and delay the features
+delayer = Delayer(delays=[0, 1, 2, 3, 4])
+x_delayed = delayer.fit_transform(x[:, None])
+
+###############################################################################
+# In the next cell we are plotting six lines. The subplot at the top shows the
+# simulated BOLD response, while the other subplots show the simulated feature
+# at different delays. The effect of the delayer is clear: it creates multiple
+# copies of the original feature shifted forward in time by how many samples we
+# requested (in this case, from 0 to 4 samples, which correspond to 0, 2, 4, 6,
+# and 8 s in time with a 2 s TR).
+#
+# When these delayed features are used to fit a voxelwise encoding model, the
+# brain response :math:`y` at time :math:`t` is simultaneously modeled by the
+# feature :math:`x` at times :math:`t-0, t-2, t-4, t-6, t-8`. In the remaining
+# of this example we will see that this method improves model prediction accuracy
+# and it allows to account for the underlying shape of the hemodynamic response
+# function.
+
+import matplotlib.pyplot as plt
+fig, axs = plt.subplots(6, 1, figsize=(8, 6.5), constrained_layout=True, 
+        sharex=True)
+times = np.arange(n_trs)*TR
+
+axs[0].plot(times, y, color="r")
+axs[0].set_title("BOLD response")
+for i, (ax, xx) in enumerate(zip(axs.flat[1:], x_delayed.T)):
+  ax.plot(times, xx, color='k')
+  ax.set_title("$x(t - {0:.0f})$ (feature delayed by {1} sample{2})".format(
+      i*TR, i, "" if i == 1 else "s"))
+for ax in axs.flat:
+  ax.axvline(40, color='gray')
+  ax.set_yticks([])
+_ = axs[-1].set_xlabel("Time [s]")
+plt.show()
+
+
 ###############################################################################
 # Compare with a model without delays
 # -----------------------------------
@@ -171,7 +229,6 @@
 # diagonal corresponds to identical prediction accuracy for both models. A
 # distibution deviating from the diagonal means that one model has better
 # prediction accuracy than the other.
-import matplotlib.pyplot as plt
 from voxelwise_tutorials.viz import plot_hist2d
 
 ax = plot_hist2d(scores_no_delay, scores)

tutorials/movies_3T/04_plot_motion_energy_model.py

@@ -147,7 +147,7 @@
 ###############################################################################
 # Plot the model performances
 # ---------------------------
-# The performances are computed using the math:`R^2` scores.
+# The performances are computed using the :math:`R^2` scores.
 
 import matplotlib.pyplot as plt
 from voxelwise_tutorials.viz import plot_flatmap_from_mapper

tutorials/movies_3T/06_extract_motion_energy.py

@@ -37,7 +37,7 @@
 # Load the stimuli images
 # -----------------------
 #
-# Here the data is not loaded in memory, we only take a peak at the data shape.
+# Here the data is not loaded in memory, we only take a peek at the data shape.
 
 import h5py