Skip to content

Commit e7a177d

Browse files
JesseLivezeyJesseLivezey
committed
added tutorial notebook
1 parent d41f291 commit e7a177d

1 file changed

Lines changed: 334 additions & 0 deletions

File tree

‎notebooks/tutorial.ipynb‎

Lines changed: 334 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,334 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import os, urllib\n",
10+
"import numpy as np\n",
11+
"import matplotlib.pyplot as plt\n",
12+
"from sklearn.decomposition import PCA\n",
13+
"from sklearn.linear_model import Ridge as RR\n",
14+
"from sklearn.metrics import r2_score\n",
15+
"\n",
16+
"from dca.data_util import load_sabes_data\n",
17+
"from dca.dca import DynamicalComponentsAnalysis as DCA"
18+
]
19+
},
20+
{
21+
"cell_type": "markdown",
22+
"metadata": {},
23+
"source": [
24+
"# Dimensionality and neural data?\n",
25+
"Due to some combination of simple experimental paradigms and more fundamental unknown principles of cortical processing, the neural data we recorded is often much lower dimensional than the recording dimensionality. This tutorial explores a few different ways of thinking about neural dimensionality. It contrasts DCA with the commonly used PCA.\n",
26+
"\n",
27+
"# Run once to download the M1 reaching data\n",
28+
"The file is 1.1 GB, so it may take a few minutes to download. This data is from the Sabes lab and is recordings from M1 while a monkey is reaching to different locations on a grid. This is the same data used in the DCA paper. In this tutorial, we won't cross validate the results as was done in the paper to keep things simple.\n",
29+
"\n",
30+
"More information and data can be found here:\n",
31+
"\n",
32+
"O'Doherty, Joseph E., Cardoso, Mariana M. B., Makin, Joseph G., & Sabes, Philip N. (2017). Nonhuman Primate Reaching with Multichannel Sensorimotor Cortex Electrophysiology [Data set]. Zenodo. http://doi.org/10.5281/zenodo.583331"
33+
]
34+
},
35+
{
36+
"cell_type": "code",
37+
"execution_count": null,
38+
"metadata": {},
39+
"outputs": [],
40+
"source": [
41+
"fname = 'indy_20160627_01.mat'\n",
42+
"if not os.path.isfile(fname): # check if file was already downloaded\n",
43+
" urllib.request.urlretrieve('https://zenodo.org/record/583331/files/indy_20160627_01.mat?download=1', fname)"
44+
]
45+
},
46+
{
47+
"cell_type": "markdown",
48+
"metadata": {},
49+
"source": [
50+
"# Let's load and visualize some of the data\n",
51+
"We'll use 50ms bins and preprocess the data by removing neurons with very low firing rates, square-root transforming the spike counts, and high-pass filtering the data to remove slow nonstationarity (30s timescale).\n",
52+
"\n",
53+
"This should load a dictionary that contains the (preprocessed) spike counts for 109 neurons along with the cursor location sampled at the same rate. We'll visualize the spike raster and cursor location for 1 minute of data."
54+
]
55+
},
56+
{
57+
"cell_type": "code",
58+
"execution_count": null,
59+
"metadata": {},
60+
"outputs": [],
61+
"source": [
62+
"data = load_sabes_data(fname, bin_width_s=.05, preprocess=True)"
63+
]
64+
},
65+
{
66+
"cell_type": "code",
67+
"execution_count": null,
68+
"metadata": {},
69+
"outputs": [],
70+
"source": [
71+
"keys = data.keys()\n",
72+
"print(data.keys())\n",
73+
"print(*[(key, data[key].shape) for key in keys])"
74+
]
75+
},
76+
{
77+
"cell_type": "code",
78+
"execution_count": null,
79+
"metadata": {},
80+
"outputs": [],
81+
"source": [
82+
"X = data['M1']\n",
83+
"Xn = X / X.std(axis=0, keepdims=True) # normalized version will be used later\n",
84+
"Y = data['cursor']"
85+
]
86+
},
87+
{
88+
"cell_type": "code",
89+
"execution_count": null,
90+
"metadata": {},
91+
"outputs": [],
92+
"source": [
93+
"plt.figure(figsize=(10, 5))\n",
94+
"plt.imshow(X[:1200].T, extent=[0, 1199*.05, 0, 108], cmap='gray_r', aspect='auto')\n",
95+
"plt.xlabel('Time (s)')\n",
96+
"plt.ylabel('Neuron')\n",
97+
"\n",
98+
"plt.figure(figsize=(5, 5))\n",
99+
"plt.plot(*Y[:1200].T, c='k')\n",
100+
"plt.xlabel('cursor x')\n",
101+
"plt.ylabel('cursor y')"
102+
]
103+
},
104+
{
105+
"cell_type": "markdown",
106+
"metadata": {},
107+
"source": [
108+
"# What is the dimensionality of the neural data?\n",
109+
"There are many ways of defining dimensionality. PCA organizes the data along projections of decreasing variance explained. If the variance of your dataset is confined to a lower dimensional subspace in the original measurement space, PCA will find this manifold. DCA, on the other hand, looks for subspaces with highly predictive dynamics, as measured by Predictive Information (PI).\n",
110+
"\n",
111+
"We'll focus on the first 30 dimensions, but you could extend the analysis out to 109. We'll look at the objective of PCA (variance explained) and DCA (PI) as a function of projection dimensionality. This is a purely unsupervised analysis of dimensionality. We can also ask how well the projections (found in an unsupervised manner) can be used to predict behavior for each method.\n",
112+
"\n",
113+
"One weakness of PCA, which motivated the development of DCA, is that it cannot distinguish high variance dynamics from high variance noise. Let's first look at the variance explained by PCA projections and their $R^2$ in predicting the behavioral data."
114+
]
115+
},
116+
{
117+
"cell_type": "code",
118+
"execution_count": null,
119+
"metadata": {},
120+
"outputs": [],
121+
"source": [
122+
"max_dim = 30\n",
123+
"lag = 4 # 200ms lag for the neural data for predicting behavior"
124+
]
125+
},
126+
{
127+
"cell_type": "code",
128+
"execution_count": null,
129+
"metadata": {},
130+
"outputs": [],
131+
"source": [
132+
"ds = np.arange(1, max_dim+1)\n",
133+
"pca_model = PCA()\n",
134+
"pca_model.fit(X)\n",
135+
"var = np.sum(pca_model.explained_variance_)\n",
136+
"\n",
137+
"pca1_scores = np.zeros(ds.size)\n",
138+
"for ii, d in enumerate(ds):\n",
139+
" Xd = pca_model.transform(X)[:, :d]\n",
140+
" rr_model = RR(alpha=1e-6)\n",
141+
" rr_model.fit(Xd[:-lag], Y[lag:])\n",
142+
" pca1_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))\n",
143+
"rr_model = RR(alpha=1e-6)\n",
144+
"rr_model.fit(X[:-lag], Y[lag:])\n",
145+
"max_score = r2_score(Y[lag:], rr_model.predict(X[:-lag]))\n",
146+
"\n",
147+
"plt.plot(ds, np.cumsum(pca_model.explained_variance_)[:ds.size] / var, label='Var. explained')\n",
148+
"plt.ylim(0, 1)\n",
149+
"plt.plot(ds, pca1_scores / max_score, label=r'$R^2')\n",
150+
"plt.legend(loc='best')\n",
151+
"plt.xlabel('Projected dimensionality')\n",
152+
"plt.ylabel('0-1 normalized metric')"
153+
]
154+
},
155+
{
156+
"cell_type": "markdown",
157+
"metadata": {},
158+
"source": [
159+
"To make visualization easier, we've normalized the y-axis so that 1 is the value of the metric at the full dimensionality of the dataset (109). So, all plots would eventually go to 1 at $d=109$.\n",
160+
"\n",
161+
"### Questions:\n",
162+
"- Does the variance explained look low dimensional?\n",
163+
"- Does the $R^2$ look low dimensional?\n",
164+
"\n",
165+
"We can run the equivalent analysis for DCA. Rather than explained variance, we'll look at PI. $R^2$ is computed in the same way across methods. This analysis will be somewhat slower since the projections need to be re-fit for each dimensionality."
166+
]
167+
},
168+
{
169+
"cell_type": "code",
170+
"execution_count": null,
171+
"metadata": {},
172+
"outputs": [],
173+
"source": [
174+
"pi = np.zeros(ds.size)\n",
175+
"dca_scores = np.zeros(ds.size)\n",
176+
"dca_model = DCA(T=10, d=109)\n",
177+
"dca_model.estimate_data_statistics(X) # only need to estimate this once\n",
178+
"max_pi = dca_model.score() # PI of data with no dimensionality reduction"
179+
]
180+
},
181+
{
182+
"cell_type": "code",
183+
"execution_count": null,
184+
"metadata": {
185+
"scrolled": true
186+
},
187+
"outputs": [],
188+
"source": [
189+
"for ii, d in enumerate(ds):\n",
190+
" print(d)\n",
191+
" dca_model.fit_projection(d=d)\n",
192+
" pi[ii] = dca_model.score()\n",
193+
" Xd = dca_model.transform(X)\n",
194+
" rr_model = RR(alpha=1e-6)\n",
195+
" rr_model.fit(Xd[:-lag], Y[lag:])\n",
196+
" dca_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))"
197+
]
198+
},
199+
{
200+
"cell_type": "code",
201+
"execution_count": null,
202+
"metadata": {},
203+
"outputs": [],
204+
"source": [
205+
"plt.plot(ds, pi / max_pi, label='PI')\n",
206+
"plt.plot(ds, dca_scores / max_score, label=r'$R^2')\n",
207+
"plt.ylim(0, 1)\n",
208+
"plt.legend(loc='best')\n",
209+
"plt.xlabel('Projected dimensionality')\n",
210+
"plt.ylabel('0-1 normalized metric')"
211+
]
212+
},
213+
{
214+
"cell_type": "markdown",
215+
"metadata": {},
216+
"source": [
217+
"### Questions:\n",
218+
"- Does the PI look low dimensional?\n",
219+
"- Does the $R^2$ look low dimensional?\n",
220+
"\n",
221+
"# Preprocessing\n",
222+
"Data analysis methods are often sensitive to certain types of preprocessing. Spiking data is sometimes variance normalized per neuron before PCA is performed. Does this change the picture of the data PCA gives?\n",
223+
"\n",
224+
"What about DCA? What does this say about the invariances encoded in the choice of metrics: PI versus variance?"
225+
]
226+
},
227+
{
228+
"cell_type": "code",
229+
"execution_count": null,
230+
"metadata": {},
231+
"outputs": [],
232+
"source": [
233+
"pca_model = PCA()\n",
234+
"pca_model.fit(Xn) # Xn rather than X\n",
235+
"var = np.sum(pca_model.explained_variance_)\n",
236+
"\n",
237+
"pca2_scores = np.zeros(ds.size)\n",
238+
"for ii, d in enumerate(ds):\n",
239+
" Xd = pca_model.transform(Xn)[:, :d]\n",
240+
" rr_model = RR(alpha=1e-6)\n",
241+
" rr_model.fit(Xd[:-lag], Y[lag:])\n",
242+
" pca2_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))\n",
243+
"\n",
244+
"plt.plot(ds, np.cumsum(pca_model.explained_variance_)[:ds.size] / var, label='Var. explained')\n",
245+
"plt.ylim(0, 1)\n",
246+
"plt.plot(ds, pca2_scores / max_score, label=r'$R^2')\n",
247+
"plt.legend(loc='best')\n",
248+
"plt.xlabel('Projected dimensionality')\n",
249+
"plt.ylabel('0-1 normalized metric')"
250+
]
251+
},
252+
{
253+
"cell_type": "code",
254+
"execution_count": null,
255+
"metadata": {},
256+
"outputs": [],
257+
"source": [
258+
"pi2 = np.zeros(ds.size)\n",
259+
"dca_scores2 = np.zeros(ds.size)\n",
260+
"dca_model = DCA(T=10, d=109)\n",
261+
"dca_model.estimate_data_statistics(Xn) # only need to estimate this once"
262+
]
263+
},
264+
{
265+
"cell_type": "code",
266+
"execution_count": null,
267+
"metadata": {},
268+
"outputs": [],
269+
"source": [
270+
"for ii, d in enumerate(ds):\n",
271+
" print(d)\n",
272+
" dca_model.fit_projection(d=d)\n",
273+
" pi2[ii] = dca_model.score()\n",
274+
" Xd = dca_model.transform(Xn)\n",
275+
" rr_model = RR(alpha=1e-6)\n",
276+
" rr_model.fit(Xd[:-lag], Y[lag:])\n",
277+
" dca_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))"
278+
]
279+
},
280+
{
281+
"cell_type": "code",
282+
"execution_count": null,
283+
"metadata": {},
284+
"outputs": [],
285+
"source": [
286+
"plt.plot(ds, pi / max_pi, label='PI')\n",
287+
"plt.plot(ds, dca_scores / max_score, label=r'$R^2')\n",
288+
"plt.ylim(0, 1)\n",
289+
"plt.legend(loc='best')\n",
290+
"plt.xlabel('Projected dimensionality')\n",
291+
"plt.ylabel('0-1 normalized metric')"
292+
]
293+
},
294+
{
295+
"cell_type": "markdown",
296+
"metadata": {},
297+
"source": [
298+
"# Ideas for other simple visualization and analysis\n",
299+
"\n",
300+
"- You could visualize the predicted cursor locations versus the true cursor locations. What features does it predict well? What does it seem to miss?\n",
301+
"- You could visualized low-dimensional projections of the neural data for PCA versus DCA. Do you see any qualitative differences?\n",
302+
"- The cursor velocities are sometimes included as variables to predict in addition or as an alternative to the location. How does this change the results?"
303+
]
304+
},
305+
{
306+
"cell_type": "code",
307+
"execution_count": null,
308+
"metadata": {},
309+
"outputs": [],
310+
"source": []
311+
}
312+
],
313+
"metadata": {
314+
"kernelspec": {
315+
"display_name": "Python 3",
316+
"language": "python",
317+
"name": "python3"
318+
},
319+
"language_info": {
320+
"codemirror_mode": {
321+
"name": "ipython",
322+
"version": 3
323+
},
324+
"file_extension": ".py",
325+
"mimetype": "text/x-python",
326+
"name": "python",
327+
"nbconvert_exporter": "python",
328+
"pygments_lexer": "ipython3",
329+
"version": "3.7.7"
330+
}
331+
},
332+
"nbformat": 4,
333+
"nbformat_minor": 4
334+
}

0 commit comments

Comments
 (0)