MNT,DOC clean up notebooks and improve FIR notebook

Author: mvdoc

tutorials/Makefile

@@ -12,8 +12,9 @@ help:
 	@echo "  merge-notebooks Merge notebook files into single files for Colab"
 	@echo "  build           Build the jupyter book"
 	@echo "  push-pages      Push built book to GitHub Pages"
+	@echo "  preview         Preview the book locally"
 
-.PHONY: help Makefile clean merge-notebooks build push-pages
+.PHONY: help Makefile clean merge-notebooks build push-pages preview
 
 clean:
 	rm -rf $(BUILDDIR)/

tutorials/notebooks/shortclips/05_fit_wordnet_model.ipynb

@@ -122,14 +122,14 @@
       "outputs": [],
       "source": [
         "from scipy.stats import zscore\n",
+        "from voxelwise_tutorials.utils import zscore_runs\n",
         "\n",
         "# indice of first sample of each run\n",
         "run_onsets = load_hdf5_array(file_name, key=\"run_onsets\")\n",
         "print(run_onsets)\n",
         "\n",
         "# zscore each training run separately\n",
-        "Y_train = np.split(Y_train, run_onsets[1:])\n",
-        "Y_train = np.concatenate([zscore(run, axis=0) for run in Y_train], axis=0)\n",
+        "Y_train = zscore_runs(Y_train, run_onsets)\n",
         "# zscore each test run separately\n",
         "Y_test = zscore(Y_test, axis=1)"
       ]
@@ -156,7 +156,6 @@
       "source": [
         "Y_test = Y_test.mean(0)\n",
         "# We need to zscore the test data again, because we took the mean across repetitions.\n",
-        "# This averaging step makes the standard deviation approximately equal to 1/sqrt(n_repeats)\n",
         "Y_test = zscore(Y_test, axis=0)\n",
         "\n",
         "print(\"(n_samples_test, n_voxels) =\", Y_test.shape)"

tutorials/notebooks/shortclips/06_visualize_hemodynamic_response.ipynb

@@ -84,6 +84,7 @@
         "import numpy as np\n",
         "from scipy.stats import zscore\n",
         "from voxelwise_tutorials.io import load_hdf5_array\n",
+        "from voxelwise_tutorials.utils import zscore_runs\n",
         "\n",
         "file_name = os.path.join(directory, \"responses\", f\"{subject}_responses.hdf\")\n",
         "Y_train = load_hdf5_array(file_name, key=\"Y_train\")\n",
@@ -96,8 +97,7 @@
         "run_onsets = load_hdf5_array(file_name, key=\"run_onsets\")\n",
         "\n",
         "# zscore each training run separately\n",
-        "Y_train = np.split(Y_train, run_onsets[1:])\n",
-        "Y_train = np.concatenate([zscore(run, axis=0) for run in Y_train], axis=0)\n",
+        "Y_train = zscore_runs(Y_train, run_onsets)\n",
         "# zscore each test run separately\n",
         "Y_test = zscore(Y_test, axis=1)"
       ]
@@ -287,14 +287,16 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Intermission: understanding delays\n",
+        "## Understanding delays\n",
         "\n",
         "To have an intuitive understanding of what we accomplish by delaying the\n",
         "features before model fitting, we will simulate one voxel and a single\n",
         "feature. We will then create a ``Delayer`` object (which was used in the\n",
-        "previous pipeline) and visualize its effect on our single feature. Let's\n",
-        "start by simulating the data.\n",
-        "\n"
+        "previous pipeline) and visualize its effect on our single feature. \n",
+        "\n",
+        "Let's start by simulating the data. We assume a simple scenario in which an event in\n",
+        "our experiment occurred at $t = 20$ seconds and lasted for 10 seconds. For each timepoint, our simulated feature\n",
+        "is a simple variable that indicates whether the event occurred or not."
       ]
     },
     {
@@ -305,71 +307,83 @@
       },
       "outputs": [],
       "source": [
-        "# number of total trs\n",
-        "n_trs = 50\n",
-        "# repetition time for the simulated data\n",
-        "TR = 2.0\n",
-        "rng = np.random.RandomState(42)\n",
-        "y = rng.randn(n_trs)\n",
-        "x = np.zeros(n_trs)\n",
-        "# add some arbitrary value to our feature\n",
-        "x[15:20] = 0.5\n",
-        "x += rng.randn(n_trs) * 0.1  # add some noise\n",
+        "from voxelwise_tutorials.delays_toy import create_voxel_data\n",
         "\n",
-        "# create a delayer object and delay the features\n",
-        "delayer = Delayer(delays=[0, 1, 2, 3, 4])\n",
-        "x_delayed = delayer.fit_transform(x[:, None])"
+        "# simulate an activation pulse at 20 s for 10 s of duration\n",
+        "simulated_X, simulated_Y, times = create_voxel_data(onset=20, duration=10)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We next plot the simulated data. In this toy example, we assumed a \"canonical\" \n",
+        "hemodynamic response function (HRF) (a double gamma function). This is an idealized\n",
+        "HRF that is often used in the literature to model the BOLD response. In practice, \n",
+        "however, the HRF can vary significantly across brain areas.\n",
+        "\n",
+        "Because of the HRF, notice that even though the event occurred at $t = 20$ seconds, \n",
+        "the BOLD response is delayed in time. "
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "import matplotlib.pyplot as plt\n",
+        "from voxelwise_tutorials.delays_toy import plot_delays_toy\n",
+        "\n",
+        "plot_delays_toy(simulated_X, simulated_Y, times)\n",
+        "plt.show()"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "In the next cell we are plotting six lines. The subplot at the top shows the\n",
-        "simulated BOLD response, while the other subplots show the simulated feature\n",
-        "at different delays. The effect of the delayer is clear: it creates multiple\n",
+        "We next create a `Delayer` object and use it to delay the simulated feature. \n",
+        "The effect of the delayer is clear: it creates multiple\n",
         "copies of the original feature shifted forward in time by how many samples we\n",
         "requested (in this case, from 0 to 4 samples, which correspond to 0, 2, 4, 6,\n",
         "and 8 s in time with a 2 s TR).\n",
         "\n",
         "When these delayed features are used to fit a voxelwise encoding model, the\n",
         "brain response $y$ at time $t$ is simultaneously modeled by the\n",
-        "feature $x$ at times $t-0, t-2, t-4, t-6, t-8$. In the remaining\n",
-        "of this example we will see that this method improves model prediction\n",
-        "accuracy and it allows to account for the underlying shape of the hemodynamic\n",
-        "response function.\n",
-        "\n"
+        "feature $x$ at times $t-0, t-2, t-4, t-6, t-8$. For example, the time sample highlighted\n",
+        "in the plot below ($t = 30$ seconds) is modeled by the features at \n",
+        "$t = 30, 28, 26, 24, 22$ seconds."
       ]
     },
     {
       "cell_type": "code",
       "execution_count": null,
-      "metadata": {
-        "collapsed": false
-      },
+      "metadata": {},
       "outputs": [],
       "source": [
-        "import matplotlib.pyplot as plt\n",
+        "# create a delayer object and delay the features\n",
+        "delayer = Delayer(delays=[0, 1, 2, 3, 4])\n",
+        "simulated_X_delayed = delayer.fit_transform(simulated_X[:, None])\n",
         "\n",
-        "fig, axs = plt.subplots(6, 1, figsize=(6, 6), constrained_layout=True, sharex=True)\n",
-        "times = np.arange(n_trs) * TR\n",
-        "\n",
-        "axs[0].plot(times, y, color=\"r\")\n",
-        "axs[0].set_title(\"BOLD response\")\n",
-        "for i, (ax, xx) in enumerate(zip(axs.flat[1:], x_delayed.T)):\n",
-        "    ax.plot(times, xx, color=\"k\")\n",
-        "    ax.set_title(\n",
-        "        \"$x(t - {0:.0f})$ (feature delayed by {1} sample{2})\".format(\n",
-        "            i * TR, i, \"\" if i == 1 else \"s\"\n",
-        "        )\n",
-        "    )\n",
-        "for ax in axs.flat:\n",
-        "    ax.axvline(40, color=\"gray\")\n",
-        "    ax.set_yticks([])\n",
-        "_ = axs[-1].set_xlabel(\"Time [s]\")\n",
+        "# plot the simulated data and highlight t = 30\n",
+        "plot_delays_toy(simulated_X_delayed, simulated_Y, times, highlight=30)\n",
         "plt.show()"
       ]
     },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "This simple example shows how the delayed features take into account of the HRF. \n",
+        "This approach is often referred to as a \"finite impulse response\" (FIR) model.\n",
+        "By delaying the features, the regression model learns the weights for each voxel \n",
+        "separately. Therefore, the FIR approach is able to adapt to the shape of the HRF in each \n",
+        "voxel, without assuming a fixed canonical HRF shape. \n",
+        "As we will see in the remaining of this notebook, this approach improves model \n",
+        "prediction accuracy significantly."
+      ]
+    },
     {
       "cell_type": "markdown",
       "metadata": {},

tutorials/notebooks/shortclips/08_fit_motion_energy_model.ipynb

@@ -91,6 +91,7 @@
         "import numpy as np\n",
         "from scipy.stats import zscore\n",
         "from voxelwise_tutorials.io import load_hdf5_array\n",
+        "from voxelwise_tutorials.utils import zscore_runs\n",
         "\n",
         "file_name = os.path.join(directory, \"responses\", f\"{subject}_responses.hdf\")\n",
         "Y_train = load_hdf5_array(file_name, key=\"Y_train\")\n",
@@ -103,8 +104,7 @@
         "run_onsets = load_hdf5_array(file_name, key=\"run_onsets\")\n",
         "\n",
         "# zscore each training run separately\n",
-        "Y_train = np.split(Y_train, run_onsets[1:])\n",
-        "Y_train = np.concatenate([zscore(run, axis=0) for run in Y_train], axis=0)\n",
+        "Y_train = zscore_runs(Y_train, run_onsets)\n",
         "# zscore each test run separately\n",
         "Y_test = zscore(Y_test, axis=1)"
       ]
@@ -486,7 +486,7 @@
         "semantic information.\n",
         "\n",
         "To better disentangle the two feature spaces, we developed a joint model\n",
-        "called `banded ridge regression` {cite}`nunez2019,dupre2022`, which fits multiple feature spaces\n",
+        "called **banded ridge regression** {cite}`nunez2019,dupre2022`, which fits multiple feature spaces\n",
         "simultaneously with optimal regularization for each feature space. This model\n",
         "is described in the next example.\n",
         "\n"

tutorials/notebooks/shortclips/09_fit_banded_ridge_model.ipynb

@@ -77,6 +77,7 @@
         "import numpy as np\n",
         "from scipy.stats import zscore\n",
         "from voxelwise_tutorials.io import load_hdf5_array\n",
+        "from voxelwise_tutorials.utils import zscore_runs\n",
         "\n",
         "file_name = os.path.join(directory, \"responses\", f\"{subject}_responses.hdf\")\n",
         "Y_train = load_hdf5_array(file_name, key=\"Y_train\")\n",
@@ -89,8 +90,7 @@
         "run_onsets = load_hdf5_array(file_name, key=\"run_onsets\")\n",
         "\n",
         "# zscore each training run separately\n",
-        "Y_train = np.split(Y_train, run_onsets[1:])\n",
-        "Y_train = np.concatenate([zscore(run, axis=0) for run in Y_train], axis=0)\n",
+        "Y_train = zscore_runs(Y_train, run_onsets)\n",
         "# zscore each test run separately\n",
         "Y_test = zscore(Y_test, axis=1)"
       ]

voxelwise_tutorials/delays_toy.py

@@ -0,0 +1,154 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import gamma
+
+
+def simulate_bold(stimulus_array, TR=2.0, duration=32.0):
+    """Simulate BOLD signal by convolving a stimulus time series with a canonical HRF.
+
+    Parameters:
+    -----------
+    stimulus_array : array-like
+        Binary array representing stimulus onset (1) and absence (0)
+    TR : float, optional
+        Repetition time of the scanner in seconds (default: 2.0s)
+    duration : float, optional
+        Total duration to model the HRF in seconds (default: 32.0s)
+
+    Returns:
+    --------
+    bold_signal : ndarray
+        The simulated BOLD signal after convolution
+    """
+    # Create time array based on TR
+    t = np.arange(0, duration, TR)
+
+    # Define canonical double-gamma HRF
+    # Parameters for the canonical HRF (based on SPM defaults)
+    peak_delay = 6.0  # delay of peak in seconds
+    undershoot_delay = 16.0  # delay of undershoot
+    peak_disp = 1.0  # dispersion of peak
+    undershoot_disp = 1.0  # dispersion of undershoot
+    peak_amp = 1.0  # amplitude of peak
+    undershoot_amp = 0.1  # amplitude of undershoot
+
+    # Compute positive and negative gamma functions
+    pos_gamma = peak_amp * gamma.pdf(t, peak_delay / peak_disp, scale=peak_disp)
+    neg_gamma = undershoot_amp * gamma.pdf(
+        t, undershoot_delay / undershoot_disp, scale=undershoot_disp
+    )
+
+    # Canonical HRF = positive gamma - negative gamma
+    hrf = pos_gamma - neg_gamma
+
+    # Normalize HRF to have sum = 1
+    hrf = hrf / np.sum(hrf)
+
+    # Perform convolution
+    bold_signal = np.convolve(stimulus_array, hrf, mode="full")
+
+    # Return only the part of the signal that matches the length of the input
+    return bold_signal[: len(stimulus_array)]
+
+
+def create_voxel_data(n_trs=50, TR=2.0, onset=30, duration=10, random_seed=42):
+    """Create a toy dataset with a single voxel and a single feature.
+
+    Parameters
+    ----------
+    n_trs : int
+        Number of time points (TRs).
+    TR : float
+        Repetition time in seconds.
+    activation_t : float
+        Time point of activation onset in seconds.
+    activation_duration : float
+        Duration of activation in seconds.
+    random_seed : int
+        Seed for random number generation.
+
+    Returns
+    -------
+    X : array of shape (n_trs,)
+        The generated feature.
+    Y : array of shape (n_trs,)
+        The generated voxel data.
+    times : array of shape (n_trs,)
+        The time points corresponding to the voxel data.
+    """
+    if onset > n_trs * TR:
+        raise ValueError("onset must be less than n_trs * TR")
+    if duration > n_trs * TR:
+        raise ValueError("duration must be less than n_trs * TR")
+    if duration < 0:
+        raise ValueError("duration must be greater than 0")
+    if onset < 0:
+        raise ValueError("onset must be greater than 0")
+    rng = np.random.RandomState(random_seed)
+    # figure out slices of activation
+    activation = slice(int(onset / TR), int((onset + duration) / TR))
+    X = np.zeros(n_trs)
+    # add some arbitrary value to our feature
+    X[activation] = 1
+    Y = simulate_bold(X, TR=TR)
+    # add some noise
+    Y += rng.randn(n_trs) * 0.1
+    times = np.arange(n_trs) * TR
+    return X, Y, times
+
+
+def plot_delays_toy(X_delayed, Y, times, highlight=None):
+    """Creates a figure showing a BOLD response and delayed versions of a stimulus.
+
+    Parameters
+    ----------
+    X_delayed : ndarray, shape (n_timepoints, n_delays)
+        The delayed stimulus, where each column corresponds to a different delay.
+    Y : ndarray, shape (n_timepoints,)
+        The BOLD response time series.
+    times : ndarray, shape (n_timepoints,)
+        Time points in seconds.
+    highlight : float or None, optional
+        Time point to highlight in the plot (default: 30 seconds).
+
+    Returns
+    -------
+    axs : ndarray
+        Array of matplotlib axes objects containing the plots.
+    """
+    if X_delayed.ndim == 1:
+        X_delayed = X_delayed[:, np.newaxis]
+    n_delays = X_delayed.shape[1]
+    n_rows = n_delays + 1
+    TR = times[1] - times[0]
+
+    fig, axs = plt.subplots(
+        n_rows,
+        1,
+        figsize=(6, n_rows),
+        constrained_layout=True,
+        sharex=True,
+        sharey=True,
+    )
+    axs[0].plot(times, Y, color="r")
+    axs[0].set_title("BOLD response")
+    if len(axs) == 2:
+        axs[1].set_title("Feature")
+        axs[1].plot(times, X_delayed, color="k")
+    else:
+        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"
+                )
+            )
+    if highlight is not None:
+        for ax in axs.flat:
+            ax.axvline(highlight, color="gray")
+            ax.set_yticks([])
+    _ = axs[-1].set_xlabel("Time [s]")
+    # add more margin at the top of the y axis
+    ylim = axs[0].get_ylim()
+    axs[0].set_ylim(ylim[0], ylim[1] * 1.2)
+    return axs