Merge pull request #482 from pynapple-org/476-template-matching

476 template matching

Author: gviejo

.github/workflows/main.yml

@@ -90,7 +90,13 @@ jobs:
         with:
           directory: "doc/_build/html"
           # The directory to scan
-          arguments: --checks Links,Scripts --ignore-urls "https://fonts.gstatic.com,https://mkdocs-gallery.github.io,./doc/_build/html/_static/,https://www.nature.com/articles/s41593-022-01020-w" --assume-extension --check-external-hash --ignore-status-codes 403 --ignore-files "/.+\/html\/_static\/.+/"
+          arguments: 
+            --checks Links,Scripts 
+            --ignore-urls "https://fonts.gstatic.com,https://mkdocs-gallery.github.io,./doc/_build/html/_static/,https://www.nature.com/articles/s41593-022-01020-w,https://elifesciences.org/reviewed-preprints/85786" 
+            --assume-extension 
+            --check-external-hash 
+            --ignore-status-codes 403 
+            --ignore-files "/.+\/html\/_static\/.+/"
           # The arguments to pass to HTMLProofer
 
   check:

doc/examples.md

@@ -24,7 +24,7 @@ Streaming data from Dandi <examples/tutorial_pynapple_dandi>
 :::{card} 
 ```{toctree}
 :maxdepth: 3
-Computing calcium imaging tuning curves <examples/tutorial_calcium_imaging>
+Analyzing calcium imaging data <examples/tutorial_calcium_imaging>
 ```
 :::
 

doc/examples/tutorial_HD_dataset.md

@@ -139,17 +139,12 @@ color = xr.DataArray(
 To make the tuning curves look nice, we will smooth them before plotting:
 
 ```{code-cell} ipython3
-from scipy.ndimage import gaussian_filter1d
-
-tmp = np.concatenate(
-    [
-        tuning_curves.values, 
-        tuning_curves.values, 
-        tuning_curves.values
-    ], 
-    axis=1)
-tmp = gaussian_filter1d(tmp, sigma=3, axis=1)
-tuning_curves.values = tmp[:, tuning_curves.shape[1]:2*tuning_curves.shape[1]]
+tuning_curves.values = scipy.ndimage.gaussian_filter1d(
+    tuning_curves.values, 
+    sigma=3, 
+    axis=1, 
+    mode="wrap" # important for circular variables!
+)
 ```
 
 What does this look like? Let's plot them!
@@ -180,7 +175,7 @@ Now that we have HD tuning curves, we can go one step further. Using only the po
 We will then compare this to the real head-direction of the animal, and discover that population activity in the ADn indeed codes for HD.
 
 To decode the population activity, we will be using a bayesian decoder as implemented in Pynapple.
-Again, just a single line of code!
+Again, just a single line of code:
 
 ```{code-cell} ipython3
 decoded, proba_feature = nap.decode_bayes(
@@ -197,7 +192,7 @@ What does this look like?
 print(decoded)
 ```
 
-The variable 'decoded' contains the most probable angle, and 'proba_feature' contains the probability of a given angular bin at a given time point:
+The variable ``decoded`` contains the most probable angle, and ``proba_feature`` contains the probability of a given angular bin at a given time point:
 
 ```{code-cell} ipython3
 print(proba_feature)
@@ -227,15 +222,15 @@ plt.ylabel("Neurons")
 plt.show()
 ```
 
-From this plot, we can see that the decoder is able to estimate the head-direction based on the population activity in ADn. Amazing!
+From this plot, we can see that the decoder is able to estimate the head-direction based on the population activity in ADn.
 
-What does the probability distribution in this example event look like?
-Ideally, the bins with the highest probability will correspond to the bins having the most spikes. Let's plot the probability matrix to visualize this.
+We can also visualize the probability distribution.
+Ideally, the bins with the highest probability correspond to the bins with the most spikes. 
 
 ```{code-cell} ipython3
 smoothed = scipy.ndimage.gaussian_filter(
     proba_feature, 1
-)  # Smoothening the probability distribution
+)  # Smoothing the probability distribution
 
 # Create a DataFrame with the smoothed distribution
 p_feature = pd.DataFrame(
@@ -270,8 +265,7 @@ plt.show()
 ```
 
 <!-- #region -->
-From this probability distribution, we observe that the decoded HD closely matches the actual HD.
-Hence, the population activity in ADn is a reliable estimate of the heading direction of the animal.
+The decoded HD (dashed grey line) closely matches the actual HD (solid white line), and thus the population activity in ADn is a reliable estimate of the heading direction of the animal.
 
 I hope this tutorial was helpful. If you have any questions, comments or suggestions, please feel free to reach out to the Pynapple Team!
 

doc/examples/tutorial_calcium_imaging.md

@@ -23,6 +23,7 @@ The NWB file for the example is hosted on [OSF](https://osf.io/sbnaw). We show b
 
 ```{code-cell} ipython3
 :tags: [hide-output]
+import numpy as np
 import pynapple as nap
 import matplotlib.pyplot as plt
 import seaborn as sns
@@ -91,7 +92,6 @@ As you can see, we have a longer recording for our tracking of the animal's head
 
 ```{code-cell} ipython3
 transients.time_support
-angle.time_support
 ```
 
 ***
@@ -133,25 +133,180 @@ We start by finding the midpoint of the recording, using the function [`get_inte
 center = transients.time_support.get_intervals_center()
 
 halves = nap.IntervalSet(
-	start = [transients.time_support.start[0], center.t[0]],
+    start = [transients.time_support.start[0], center.t[0]],
     end = [center.t[0], transients.time_support.end[0]]
-    )
+)
 ```
 
 Now, we can compute the tuning curves for each half of the recording and plot the tuning curves again.
 
 ```{code-cell} ipython3
-half1 = nap.compute_tuning_curves(transients, angle, bins = 120, epochs = halves.loc[[0]])
-half2 = nap.compute_tuning_curves(transients, angle, bins = 120, epochs = halves.loc[[1]])
+tuning_curves_half1 = nap.compute_tuning_curves(transients, angle, bins = 120, epochs = halves.loc[[0]])
+tuning_curves_half2 = nap.compute_tuning_curves(transients, angle, bins = 120, epochs = halves.loc[[1]])
 
 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
-set_metadata(half1[4]).plot(ax=ax1)
+set_metadata(tuning_curves_half1[4]).plot(ax=ax1)
 ax1.set_title("First half")
-set_metadata(half2[4]).plot(ax=ax2)
+set_metadata(tuning_curves_half2[4]).plot(ax=ax2)
 ax2.set_title("Second half")
 plt.show()
 ```
 
+***
+Calcium decoding
+---------------------
+
+Given some tuning curves, we can also try to decode head direction from the population activity.
+For calcium imaging data, Pynapple has `decode_template`, which implements a template matching algorithm.
+
+```{code-cell} ipython3
+epochs = nap.IntervalSet([50, 150])
+decoded, dist = nap.decode_template(
+    tuning_curves=tuning_curves,
+    data=transients,
+    epochs=epochs,
+    bin_size=0.1,
+    metric="correlation",
+)
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+# normalize distance for better visualization
+dist_norm = (dist - np.min(dist.values, axis=1, keepdims=True)) / np.ptp(
+    dist.values, axis=1, keepdims=True
+)
+
+fig, (ax1, ax2, ax3) = plt.subplots(figsize=(8, 8), nrows=3, ncols=1, sharex=True)
+ax1.plot(angle.restrict(epochs), label="True")
+ax1.scatter(decoded.times(), decoded.values, label="Decoded", c="orange")
+ax1.legend(frameon=False, bbox_to_anchor=(1.0, 1.0))
+ax1.set_ylabel("Angle [rad]")
+
+im = ax2.imshow(
+    dist.values.T, 
+    aspect="auto", 
+    origin="lower", 
+    cmap="inferno_r", 
+    extent=(epochs.start[0], epochs.end[0], 0.0, 2*np.pi)
+)
+ax2.set_ylabel("Angle [rad]")
+cbar_ax2 = fig.add_axes([0.95, ax2.get_position().y0, 0.015, ax2.get_position().height])
+fig.colorbar(im, cax=cbar_ax2, label="Distance")
+
+im = ax3.imshow(
+    dist_norm.values.T, 
+    aspect="auto", 
+    origin="lower", 
+    cmap="inferno_r", 
+    extent=(epochs.start[0], epochs.end[0], 0.0, 2*np.pi)
+)
+cbar_ax3 = fig.add_axes([0.95, ax3.get_position().y0, 0.015, ax3.get_position().height])
+fig.colorbar(im, cax=cbar_ax3, label="Norm. distance")
+ax3.set_xlabel("Time (s)")
+ax3.set_ylabel("Angle [rad]")
+plt.show()
+```
+
+The distance metric you choose can influence how well we decode.
+Internally, ``decode_template`` uses `scipy.spatial.distance.cdist` to compute the distance matrix; 
+you can take a look at [its documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html) 
+to see which metrics are supported. Here are a couple examples:
+
+```{code-cell} ipython3
+:tags: [hide-input]
+metrics = [
+    "chebyshev",
+    "dice",
+    "canberra",
+    "sqeuclidean",
+    "minkowski",
+    "euclidean",
+    "cityblock",
+    "mahalanobis",
+    "correlation",
+    "cosine",
+    "seuclidean",
+    "braycurtis",
+    "jensenshannon",
+]
+
+fig, axs = plt.subplots(5, 1, figsize=(8,12), sharex=True, sharey=True)
+for metric, ax in zip(metrics[-5:], axs.flatten()):
+    decoded, dist = nap.decode_template(
+        tuning_curves=tuning_curves,
+        data=transients,
+        bin_size=0.1,
+        metric=metric,
+        epochs=epochs,
+    )
+    # normalize distance for better visualization
+    dist_norm = (dist - np.min(dist.values, axis=1, keepdims=True)) / np.ptp(
+        dist.values, axis=1, keepdims=True
+    )
+    ax.plot(angle.restrict(epochs), label="True")
+    im = ax.imshow(
+        dist_norm.values.T, 
+        aspect="auto", 
+        origin="lower", 
+        cmap="inferno_r", 
+        extent=(epochs.start[0], epochs.end[0], 0.0, 2*np.pi)
+    )
+    if metric != metrics[-1]:
+        ax.spines['bottom'].set_visible(False)
+        ax.tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False)
+        ax.set_yticks([])
+        ax.spines['left'].set_visible(False)
+    ax.set_ylabel(metric)
+cbar_ax = fig.add_axes([0.92, ax.get_position().y0, 0.015, ax.get_position().height])
+cbar=fig.colorbar(im, cax=cbar_ax)
+cbar.set_label("Norm. distance")
+ax.set_xlabel("Time (s)")
+plt.show()
+```
+
+We recommend trying a bunch to see which one works best for you.
+In the case of head direction, we can quantify how well we decode using the absolute angular error.
+To get a fair estimate of error, we will compute the tuning curves on the first half of the data 
+and compute the error for predictions of the second half.
+
+```{code-cell} ipython3
+def absolute_angular_error(x, y):
+    return np.abs(np.angle(np.exp(1j * (x - y))))
+
+# Compute errors
+errors = {}
+for metric in metrics:
+    decoded, dist = nap.decode_template(
+        tuning_curves=tuning_curves_half1,
+        data=transients,
+        bin_size=0.1,
+        metric=metric,
+        epochs=halves.loc[[1]],
+    )
+    errors[metric] = absolute_angular_error(
+        angle.interpolate(decoded).values, decoded.values
+    )
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+sorted_items = sorted(errors.items(), key=lambda item: np.median(item[1]))
+sorted_labels, sorted_values = zip(*sorted_items)
+
+fig, ax = plt.subplots(figsize=(8, 8))
+bp = ax.boxplot(
+    x=sorted_values,
+    tick_labels=sorted_labels,
+    vert=False,
+    showfliers=False
+)
+ax.set_xlabel("Angular error [rad]")
+plt.show()
+```
+
+In this case, `jensenshannon` yields the lowest angular error.
+
 :::{card}
 Authors
 ^^^