%matplotlib inline
Download
This notebook can be downloaded as 01_fundamentals_of_pynapple.ipynb. See the button at the top right to download as markdown or pdf.
Jupyter Lab tip
Newer versions of Jupyter Lab have addressed an issue with skipping around the notebook while scrolling. To make sure this fix is enabled, in the Jupyter Lab GUI, navigate to Settings > Settings Editor > Notebook and scroll down to the Windowing mode setting and make sure it is set to contentVisibility.
Also reminder to presenter: Go to View > Appearance, select Simple Interface and turn off everything else to hide as many bars as possible. And maybe activate Presentation Mode.
And turn on View > Render side-by-side (shortcut Shift+R).
Learning the fundamentals of pynapple#
Learning objectives#
Instantiate pynapple objects
Make pynapple objects interact
Use numpy with pynapple
Slice pynapple objects
Add metadata to pynapple objects
Apply core functions of pynapple
Resources:
Let’s start by importing the pynapple package and matplotlib to see if everything is correctly installed.
If an import fails, you can do !pip install pynapple matplotlib in a cell to fix it.
import workshop_utils
import pynapple as nap
import matplotlib.pyplot as plt
import numpy as np
WARNING:2026-02-05 17:43:18,713:jax._src.xla_bridge:876: An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
/home/agent/workspace/rorse_ccn-software-feb-2026_main/.venv/lib/python3.12/site-packages/nemos/_documentation_utils/plotting.py:39: UserWarning: plotting functions contained within `_documentation_utils` are intended for nemos's documentation. Feel free to use them, but they will probably not work as intended with other datasets / in other contexts.
warnings.warn(
For this notebook we will work with fake data. The following cells generate a set of variables that we will use to create the different pynapple objects.
var1 = np.random.randn(100) # Variable 1
tsp1 = np.arange(100) # The timesteps of variable 1
var2 = np.random.randn(100, 3) # Variable 2
tsp2 = np.arange(0, 100, 1) # The timesteps of variable 2
col2 = ['pineapple', 'banana', 'tomato'] # The name of each columns of var2
var3 = np.random.randn(1000, 4, 5) # Variable 3
tsp3 = np.arange(0, 100, 0.1) # The timesteps of variable 3
random_times_1 = np.array([3.14, 37.0, 42.0])
random_times_2 = np.array([10, 25, 50, 70])
random_times_3 = np.sort(np.random.uniform(10, 80, 100))
starts_1 = np.array([10000, 60000, 90000]) # starts of an epoch in `ms`
ends_1 = np.array([20000, 80000, 95000]) # ends in `ms`
Instantiate pynapple objects#
This is a lot of variables to carry around. pynapple can help reduce the size of the workspace. Here we will instantiate all the different pynapple objects with the variables created above.
Let’s start with the simple ones.
Question: Can you instantiate the right pynapple objects for var1, var2 and var3? Objects should be named respectively tsd1, tsd2 and tsd3. Don’t forget the column name for var2.
tsd1 = nap.Tsd(t=tsp1, d=var1)
tsd2 = nap.TsdFrame(t=tsp2, d=var2, columns = col2)
tsd3 = nap.TsdTensor(t=tsp3, d=var3)
Question: Can you print tsd1?
print(tsd1)
Time (s)
---------- ----------
0.0 1.51965
1.0 -1.14084
2.0 -0.303191
3.0 -0.547789
4.0 -0.114832
5.0 1.04467
6.0 -0.185713
...
93.0 0.941601
94.0 1.27355
95.0 -1.20719
96.0 -1.92255
97.0 -0.0836754
98.0 1.309
99.0 -1.29158
dtype: float64, shape: (100,)
Question: Can you print tsd2?
print(tsd2)
Time (s) pineapple banana tomato
---------- ----------- ---------- ---------
0.0 0.470571 0.078878 1.02145
1.0 0.679294 1.06242 -0.439902
2.0 2.41966 1.35998 -0.188659
3.0 1.34656 -0.986796 -0.163483
4.0 -1.11433 -1.70699 2.05377
5.0 0.108866 -0.900795 0.451592
6.0 -0.190908 0.138056 -0.376092
...
93.0 -0.918117 -0.551183 0.308826
94.0 0.484801 -0.850241 -0.196055
95.0 -0.993834 0.78459 0.189073
96.0 0.465457 0.548933 -0.227605
97.0 0.390903 -0.360996 0.213092
98.0 -0.671159 0.0408639 -1.0541
99.0 0.398195 0.556957 0.497806
dtype: float64, shape: (100, 3)
Question: Can you print tsd3?
print(tsd3)
Time (s)
---------- -------------------------------
0.0 [[-0.14112 ... -0.570291] ...]
0.1 [[-0.351042 ... 0.006963] ...]
0.2 [[-1.742792 ... -0.677973] ...]
0.3 [[-0.866435 ... 0.489275] ...]
0.4 [[ 1.707552 ... -2.451565] ...]
0.5 [[0.628804 ... 1.188248] ...]
0.6 [[0.8164 ... 1.469656] ...]
...
99.3 [[-1.626194 ... 0.435116] ...]
99.4 [[-1.095492 ... 2.128051] ...]
99.5 [[-2.507073 ... -0.315374] ...]
99.6 [[-1.058535 ... 0.089936] ...]
99.7 [[ 0.748943 ... -0.577088] ...]
99.8 [[-0.460184 ... -1.511363] ...]
99.9 [[-1.396615 ... 1.208751] ...]
dtype: float64, shape: (1000, 4, 5)
Question: Can you create an IntervalSet called ep out of starts_1 and ends_1 and print it? Be careful, times given above are in ms.
ep = nap.IntervalSet(
start=starts_1,
end=ends_1, time_units='ms')
print(ep)
index start end
0 10 20
1 60 80
2 90 95
shape: (3, 2), time unit: sec.
The experiment generated a set of timestamps from 3 different channels.
Question: Can you instantiate the corresponding pynapple object (ts1, ts2, ts3) for each one of them?
ts1 = nap.Ts(t=random_times_1)
ts2 = nap.Ts(t=random_times_2)
ts3 = nap.Ts(t=random_times_3)
This is a lot of timestamps to carry around as well.
Question: Can you instantiate the right pynapple object (call it tsgroup) to group them together?
tsgroup = nap.TsGroup({0:ts1, 1:ts2, 2:ts3})
Question: Can you print the result?
print(tsgroup)
Index rate
------- -------
0 0.0391
1 0.05213
2 1.30326
Interaction between pynapple objects#
We reduced 12 variables in our workspace to 5 using pynapple. Now we can see how the objects interact.
Question: Can you print the time_support of tsgroup?
print(tsgroup.time_support)
index start end
0 3.14 79.8704
shape: (1, 2), time unit: sec.
The experiment ran from 0 to 100 seconds and as you can see, the TsGroup object shows the rate. But the rate is not accurate as it was computed over the default time_support.
Question: can you recreate the tsgroup object passing the right time_support during initialisation?
tsgroup = nap.TsGroup({0:ts1, 1:tsd2, 2:ts3}, time_support = nap.IntervalSet(0, 100))
Question: Can you print the time_support and rate to see how they changed?
print(tsgroup.time_support)
print(tsgroup.rate)
index start end
0 0 100
shape: (1, 2), time unit: sec.
0 0.03
1 1.00
2 1.00
Name: rate, dtype: float64
Now you notice that variable tsd1 has some noise. The good signal is between 10 and 30 seconds and 50 and 100.
Question: Can you create an IntervalSet object called ep_signal and use it to restrict the variable tsd1?
ep_signal = nap.IntervalSet(start=[10, 50], end=[30, 100])
tsd1 = tsd1.restrict(ep_signal)
print(tsd1)
print(tsd1.time_support)
Time (s)
---------- ----------
10.0 0.831463
11.0 -0.289732
12.0 0.827301
13.0 1.43589
14.0 -1.21685
15.0 0.781709
16.0 -1.58538
...
93.0 0.941601
94.0 1.27355
95.0 -1.20719
96.0 -1.92255
97.0 -0.0836754
98.0 1.309
99.0 -1.29158
dtype: float64, shape: (71,)
index start end
0 10 30
1 50 100
shape: (2, 2), time unit: sec.
ep_tmp = nap.IntervalSet(np.sort(np.random.uniform(0, 100, 20)))
print(ep_tmp)
index start end
0 9.95717 12.2562
1 16.9324 17.5525
2 18.3674 19.6771
3 24.4852 25.0099
4 37.0834 41.181
5 44.634 55.0222
6 56.5475 65.6439
7 69.8724 70.0106
8 86.1069 87.1082
9 95.5764 97.8659
shape: (10, 2), time unit: sec.
Question: Can you do the intersection of ep_signal and ep_tmp?
print(ep_signal.intersect(ep_tmp))
index start end
0 10 12.2562
1 16.9324 17.5525
2 18.3674 19.6771
3 24.4852 25.0099
4 50 55.0222
5 56.5475 65.6439
6 69.8724 70.0106
7 86.1069 87.1082
8 95.5764 97.8659
shape: (9, 2), time unit: sec.
workshop_utils.visualize_intervals([ep_signal, ep_tmp, ep_signal.intersect(ep_tmp)])
Question: Can you do the union of ep_signal and ep_tmp?
print(ep_signal.union(ep_tmp))
index start end
0 9.95717 30
1 37.0834 41.181
2 44.634 100
shape: (3, 2), time unit: sec.
Question: Can you visualize it?
workshop_utils.visualize_intervals([ep_signal, ep_tmp, ep_signal.union(ep_tmp)])
Question: Can you do the difference of ep_signal and ep_tmp?
print(ep_signal.set_diff(ep_tmp))
index start end
0 12.2562 16.9324
1 17.5525 18.3674
2 19.6771 24.4852
3 25.0099 30
4 55.0222 56.5475
5 65.6439 69.8724
6 70.0106 86.1069
7 87.1082 95.5764
8 97.8659 100
shape: (9, 2), time unit: sec.
Question: Can you visualize it?
workshop_utils.visualize_intervals([ep_signal, ep_tmp, ep_signal.set_diff(ep_tmp)])
Numpy & pynapple#
Pynapple objects behave similarly to numpy arrays. They can be sliced with the following syntax :
tsd[0:10] # First 10 elements
Arithmetical operations are available as well :
tsd = tsd + 1
Finally numpy functions works directly. Let’s imagine tsd3 is a movie with frame size (4,5).
Question: Can you compute the average frame along the time axis using np.mean and print the result?
print(np.mean(tsd3, 0))
[[-0.06125967 -0.01240235 -0.01773691 0.00693313 -0.01709398]
[ 0.01596417 -0.01752083 -0.00498526 -0.0195211 -0.00851743]
[ 0.0509874 0.03321083 0.02193723 -0.01876114 0.00176951]
[-0.00925244 0.0396285 0.00580824 0.01368578 -0.03701121]]
Question: can you compute the average of tsd2 for each timestamps and print it?
print(np.mean(tsd2, 1))
Time (s)
---------- -----------
0.0 0.523634
1.0 0.433937
2.0 1.19699
3.0 0.0654257
4.0 -0.255849
5.0 -0.113446
6.0 -0.142981
...
93.0 -0.386825
94.0 -0.187165
95.0 -0.00672393
96.0 0.262262
97.0 0.0809995
98.0 -0.561465
99.0 0.484319
dtype: float64, shape: (100,)
Notice how the output in the second case is still a pynapple object. In most cases, applying a numpy function will return a pynapple object if the time index is preserved.
Slicing pynapple objects#
Multiple methods exists to slice pynapple object. This parts reviews them.
IntervalSet also behaves like numpy array.
Question: Can you extract the first and last epoch of ep in a new IntervalSet?
print(ep[[0,2]])
index start end
0 10 20
1 90 95
shape: (2, 2), time unit: sec.
Sometimes you want to get a data point as close as possible in time to another timestamps.
Question: Using the get method, can you get the data point from tsd3 as close as possible to the time 50.1 seconds?
print(tsd3.get(50.1))
[[-1.21637847 -0.7562785 0.94851104 -0.96255766 1.4022758 ]
[-0.4085867 0.47283963 -0.49703136 -1.14131376 -0.48089351]
[ 0.04213654 -1.96039741 1.20745791 0.13981411 -0.80148362]
[ 0.14189759 -1.99680842 -2.73893293 -0.90623652 1.61415433]]
Metadata#
Metadata allow you to attach labels and additional information to your data objects. They are ubiquitous in neuroscience. They can be added to 3 pynapple objects :
TsGroup: to label neurons in electrophysiologyIntervalSet: to label intervalsTsdFrame: to label neurons in calcium imaging
Question: Can you run the following command tsgroup['planet'] = ['mars', 'venus', 'saturn']
tsgroup['planet'] = ['mars', 'venus', 'saturn']
Question: Can you print the result?
print(tsgroup)
Index rate planet
------- ------ --------
0 0.03 mars
1 1 venus
2 1 saturn
The object ep has 3 epochs labelled ['left', 'right', 'left'].
Question: Can you add them as a metadata column called direction?
ep['direction'] = ['left', 'right', 'left']
print(ep)
index start end direction
0 10 20 left
1 60 80 right
2 90 95 left
shape: (3, 2), time unit: sec.
The object tsd2 has 3 columns. Each column correspond to the rgb colors [(0,0,1), (0.5, 0.5, 1), (0.1, 0.2, 0.3)].
Question: Can you add them as metadata of tsd2?
tsd2['colors'] = [(0,0,1), (0.5, 0.5, 1), (0.1, 0.2, 0.3)]
print(tsd2)
Time (s) pineapple banana tomato
---------- -------------------- -------------------- --------------------
0.0 0.4705710372926658 0.07887795400250472 1.0214532743127482
1.0 0.6792940110274641 1.062419834076445 -0.43990158913769356
2.0 2.419661863397433 1.3599758706633793 -0.18865932633358853
3.0 1.346555980876706 -0.9867957297689638 -0.16348312872745002
4.0 -1.1143303626517014 -1.7069890309533662 2.0537725170452874
5.0 0.1088663575587133 -0.9007945053681811 0.45159152726650004
6.0 -0.19090750267120077 0.1380557489471097 -0.3760918062886235
...
93.0 -0.9181165680756305 -0.5511833377429957 0.30882551124771535
94.0 0.48480119834648827 -0.8502411767864009 -0.19605532370560635
95.0 -0.9938341709781316 0.7845897076805289 0.18907267519060916
96.0 0.46545707673611264 0.5489330462366683 -0.22760549699367366
97.0 0.39090254204606306 -0.360996173904464 0.21309223010969022
98.0 -0.6711585297632471 0.040863860805026644 -1.0541005731035211
99.0 0.3981946473791555 0.5569565601529481 0.4978055467219967
Metadata
colors [0. 0. 1.] [0.5 0.5 1. ] [0.1 0.2 0.3]
dtype: float64, shape: (100, 3)
You can also add metadata at initialization as a dictionary using the keyword argument metadata:
tsgroup = nap.TsGroup({0:ts1, 1:ts2, 2:ts3}, metadata={'planet':['mars','venus', 'saturn']})
print(tsgroup)
Index rate planet
------- ------- --------
0 0.0391 mars
1 0.05213 venus
2 1.30326 saturn
Metadata are accessible either as attributes (i.e. tsgroup.planet) or as dictionary-like keys (i.e. ep['direction']).
They can be used to slice objects.
Question: Can you select only the elements of tsgroup with rate below 1Hz?
print(tsgroup[tsgroup.rate<1.0])
print(tsgroup[tsgroup['rate']<1.0])
print(tsgroup.getby_threshold("rate", 1, "<"))
Index rate planet
------- ------- --------
0 0.0391 mars
1 0.05213 venus
Index rate planet
------- ------- --------
0 0.0391 mars
1 0.05213 venus
Index rate planet
------- ------- --------
0 0.0391 mars
1 0.05213 venus
Question: Can you select the intervals in ep labelled as 'left'?
print(ep[ep.direction=='left'])
index start end direction
0 10 20 left
1 90 95 left
shape: (2, 2), time unit: sec.
Special case of slicing : TsdFrame#
tsdframe = nap.TsdFrame(t=np.arange(4), d=np.random.randn(4,3),
columns = [12, 0, 1], metadata={'alpha':[2,1,0]})
print(tsdframe)
Time (s) 12 0 1
---------- ------------------- ------------------- -------------------
0.0 0.05644196234886669 -1.540174265647538 -0.5566673329617481
1.0 -1.059300574337659 0.8256134528334312 0.15514417788088422
2.0 -1.8499242285283148 -0.3701121735630015 -0.4438701925345612
3.0 -0.8469457607106837 0.5276399630306888 -0.7987527470943198
Metadata
alpha 2 1 0
dtype: float64, shape: (4, 3)
Question: What happen when you do tsdframe[0] vs tsdframe[:,0] vs tsdframe[[12,1]]
print(tsdframe[0])
[ 0.05644196 -1.54017427 -0.55666733]
print(tsdframe[:,0])
Time (s)
---------- ---------
0 0.056442
1 -1.0593
2 -1.84992
3 -0.846946
dtype: float64, shape: (4,)
try:
print(tsdframe[[12,1]])
except Exception as e:
print(e)
index 12 is out of bounds for axis 0 with size 4
Question: What happen when you do tsdframe.loc[0] and tsdframe.loc[[0,1]]
print(tsdframe.loc[0])
Time (s)
---------- ---------
0 -1.54017
1 0.825613
2 -0.370112
3 0.52764
dtype: float64, shape: (4,)
print(tsdframe.loc[[0,1]])
Time (s) 0 1
---------- ------------------- -------------------
0.0 -1.540174265647538 -0.5566673329617481
1.0 0.8256134528334312 0.15514417788088422
2.0 -0.3701121735630015 -0.4438701925345612
3.0 0.5276399630306888 -0.7987527470943198
Metadata
alpha 1 0
dtype: float64, shape: (4, 2)
Question: What happen when you do tsdframe[:,tsdframe.alpha==2]
print(tsdframe[:,tsdframe.alpha==2])
Time (s) 12
---------- -------------------
0.0 0.05644196234886669
1.0 -1.059300574337659
2.0 -1.8499242285283148
3.0 -0.8469457607106837
Metadata
alpha 2
dtype: float64, shape: (4, 1)
Core functions of pynapple#
This part focuses on the most important core functions of pynapple.
Question: Using the count function, can you count the number of events within 1 second bins for tsgroup over the ep_signal intervals?
count = tsgroup.count(1, ep_signal)
print(count)
Time (s) 0 1 2
---------- ---- ----- ------
10.5 0 1 0
11.5 0 0 2
12.5 0 0 2
13.5 0 0 2
14.5 0 0 4
15.5 0 0 2
16.5 0 0 0
...
93.5 0 0 0
94.5 0 0 0
95.5 0 0 0
96.5 0 0 0
97.5 0 0 0
98.5 0 0 0
99.5 0 0 0
Metadata
planet mars venus saturn
dtype: int64, shape: (70, 3)
Pynapple works directly with matplotlib. Passing a time series object to plt.plot will display the figure with the correct time axis.
Question: In two subplots, can you show the count and events over time?
plt.figure()
ax = plt.subplot(211)
plt.plot(count, 'o-')
plt.subplot(212, sharex=ax)
plt.plot(tsgroup.restrict(ep_signal).to_tsd(), '|')
[<matplotlib.lines.Line2D at 0x7fd8e1af8200>]
From a set of timestamps, you want to assign them a set of values with the closest point in time of another time series.
Question: Using the function value_from, can you assign values to ts2 from the tsd1 time series and call the output new_tsd?
new_tsd = ts2.value_from(tsd1)
Question: Can you plot together tsd1, ts2 and new_tsd?
plt.figure()
plt.plot(tsd1)
plt.plot(new_tsd, 'o-')
plt.plot(ts2.fillna(0), 'o')
[<matplotlib.lines.Line2D at 0x7fd8e199bfb0>]
One important aspect of data analysis is to bring data to the same size. Pynapple provides the bin_average function to downsample data.
Question: Can you downsample tsd2 to one time point every 5 seconds?
new_tsd2 = tsd2.bin_average(5.0)
Question: Can you plot the tomato column from tsd2 as well as the downsampled version?
plt.figure()
plt.plot(tsd2['tomato'])
plt.plot(new_tsd2['tomato'], 'o-')
[<matplotlib.lines.Line2D at 0x7fd8e186c170>]
For tsd1, you want to find all the epochs for which the value is above 0.0. Pynapple provides the function threshold to get 1 dimensional time series above or below a certain value.
Question: Can you print the epochs for which tsd1 is above 0.0?
ep_above = tsd1.threshold(0.0).time_support
print(ep_above)
index start end
0 10.0 10.5
1 11.5 13.5
2 14.5 15.5
3 16.5 18.5
4 20.5 21.5
5 22.5 23.5
6 25.5 29.5
... ... ...
11 69.5 71.5
12 75.5 76.5
13 79.5 82.5
14 86.5 89.5
15 90.5 91.5
16 92.5 94.5
17 97.5 98.5
shape: (18, 2), time unit: sec.
Question: can you plot tsd1 as well as the epochs for which tsd1 is above 0.0?
plt.figure()
plt.plot(tsd1)
plt.plot(tsd1.threshold(0.0), 'o-')
[plt.axvspan(s, e, alpha=0.2) for s,e in ep_above.values]
[<matplotlib.patches.Rectangle at 0x7fd8e12bb8f0>,
<matplotlib.patches.Rectangle at 0x7fd8e162b3b0>,
<matplotlib.patches.Rectangle at 0x7fd8e12ba030>,
<matplotlib.patches.Rectangle at 0x7fd8e1aaef60>,
<matplotlib.patches.Rectangle at 0x7fd8e12b9850>,
<matplotlib.patches.Rectangle at 0x7fd8e12b97f0>,
<matplotlib.patches.Rectangle at 0x7fd8e1a2d190>,
<matplotlib.patches.Rectangle at 0x7fd8e17d1910>,
<matplotlib.patches.Rectangle at 0x7fd8e17d1fa0>,
<matplotlib.patches.Rectangle at 0x7fd8e165cce0>,
<matplotlib.patches.Rectangle at 0x7fd8e12b8fe0>,
<matplotlib.patches.Rectangle at 0x7fd8e12b8b60>,
<matplotlib.patches.Rectangle at 0x7fd8e12b8a10>,
<matplotlib.patches.Rectangle at 0x7fd8e12b87d0>,
<matplotlib.patches.Rectangle at 0x7fd8e12ba8d0>,
<matplotlib.patches.Rectangle at 0x7fd8e12b8080>,
<matplotlib.patches.Rectangle at 0x7fd8e12b8560>,
<matplotlib.patches.Rectangle at 0x7fd8e17daba0>]
First high level function : compute_tuning_curves#
Pynapple provides functions for standard analysis in systems neuroscience. The first function we will try is compute_tuning_curves that calculate the response of a cell to a particular feature.
A good practice when using a function for the first time is to check the docstrings to learn how to pass the argument.
Question: can you examine the docstring of nap.compute_tuning_curves?
print(nap.compute_tuning_curves.__doc__)
Computes n-dimensional tuning curves relative to n features.
Parameters
----------
data : TsGroup, TsdFrame, Ts, Tsd
The data for which the tuning curves will be computed. This usually corresponds to the activity of the
neurons, either as spike times (TsGroup or Ts) or continuous values (TsdFrame or Tsd).
features : Tsd, TsdFrame
The features (i.e. one column per feature). This usually corresponds to behavioral variables such as
position, head direction, speed, etc.
bins : sequence or int
The bin specification:
* A sequence of arrays describing the monotonically increasing bin
edges along each dimension.
* The number of bins for each dimension (nx, ny, ... =bins)
* The number of bins for all dimensions (nx=ny=...=bins).
range : sequence, optional
A sequence of entries per feature, each an optional (lower, upper) tuple giving
the outer bin edges to be used if the edges are not given explicitly in
`bins`.
An entry of None in the sequence results in the minimum and maximum
values being used for the corresponding dimension.
The default, None, is equivalent to passing a tuple of D None values.
epochs : IntervalSet, optional
The epochs on which tuning curves are computed.
If None, the epochs are the time support of the features.
fs : float, optional
The exact sampling frequency of the features used to normalise the tuning curves.
Unit should match that of the features. If not passed, it is estimated.
feature_names : list, optional
A list of feature names. If not passed, the column names in `features` are used.
return_pandas : bool, optional
If True, the function returns a pandas.DataFrame instead of an xarray.DataArray.
Note that this will not work if the features are not 1D and that occupancy and bin edges
will not be stored as attributes.
return_counts : bool, optional
If True, does not divide the spike counts by occupancy, but returns the counts directly.
The occupancy is stored in the xarray attributes, so the division can be performed after any
particular processing steps.
If the input is a TsdFrame, this does not do anything.
Returns
-------
xarray.DataArray
A tensor containing the tuning curves with labeled bin centres.
The bin edges and occupancy are stored as attributes.
Examples
--------
In the simplest case, we can pass a group of spikes per neuron and a single feature:
>>> import pynapple as nap
>>> import numpy as np; np.random.seed(42)
>>> group = nap.TsGroup({
... 1: nap.Ts(np.arange(0, 100, 0.1)),
... 2: nap.Ts(np.arange(0, 100, 0.2))
... })
>>> feature = nap.Tsd(d=np.arange(0, 100, 0.1) % 1, t=np.arange(0, 100, 0.1))
>>> tcs = nap.compute_tuning_curves(group, feature, bins=10)
>>> tcs
<xarray.DataArray (unit: 2, 0: 10)> Size: 160B
array([[10., 10., 10., 10., 10., 10., 10., 10., 10., 10.],
[10., 0., 10., 0., 10., 0., 10., 0., 10., 0.]])
Coordinates:
* unit (unit) int64 16B 1 2
* 0 (0) float64 80B 0.045 0.135 0.225 0.315 ... 0.585 0.675 0.765 0.855
Attributes:
occupancy: [100. 100. 100. 100. 100. 100. 100. 100. 100. 100.]
bin_edges: [array([0. , 0.09, 0.18, 0.27, 0.36, 0.45, 0.54, 0.63, 0.72,...
fs: 10.0
rates: [10.01001001 5.00500501]
The function can also take multiple features, in which case it computes n-dimensional tuning curves.
We can specify the number of bins for each feature:
>>> features = nap.TsdFrame(
... d=np.stack(
... [
... np.arange(0, 100, 0.1) % 1,
... np.arange(0, 100, 0.1) % 2
... ],
... axis=1
... ),
... t=np.arange(0, 100, 0.1)
... )
>>> tcs = nap.compute_tuning_curves(group, features, bins=[5, 3])
>>> tcs
<xarray.DataArray (unit: 2, 0: 5, 1: 3)> Size: 240B
array([[[10., 10., nan],
[10., 10., 10.],
[10., nan, 10.],
[10., 10., 10.],
[nan, 10., 10.]],
...
[[ 5., 5., nan],
[ 5., 10., 0.],
[ 5., nan, 5.],
[10., 0., 5.],
[nan, 5., 5.]]])
Coordinates:
* unit (unit) int64 16B 1 2
* 0 (0) float64 40B 0.09 0.27 0.45 0.63 0.81
* 1 (1) float64 24B 0.3167 0.95 1.583
Attributes:
occupancy: [[100. 100. nan]\n [100. 50. 50.]\n [100. nan 100.]\n [ 5...
bin_edges: [array([0. , 0.18, 0.36, 0.54, 0.72, 0.9 ]), array([0. ...
fs: 10.0
rates: [10.01001001 5.00500501]
Or even specify the bin edges directly:
>>> tcs = nap.compute_tuning_curves(
... group,
... features,
... bins=[np.linspace(0, 1, 5), np.linspace(0, 2, 3)]
... )
>>> tcs
<xarray.DataArray (unit: 2, 0: 4, 1: 2)> Size: 128B
array([[[10. , 10. ],
[10. , 10. ],
[10. , 10. ],
[10. , 10. ]],
...
[[ 6.66666667, 6.66666667],
[ 5. , 5. ],
[ 3.33333333, 3.33333333],
[ 5. , 5. ]]])
Coordinates:
* unit (unit) int64 16B 1 2
* 0 (0) float64 32B 0.125 0.375 0.625 0.875
* 1 (1) float64 16B 0.5 1.5
Attributes:
occupancy: [[150. 150.]\n [100. 100.]\n [150. 150.]\n [100. 100.]]
bin_edges: [array([0. , 0.25, 0.5 , 0.75, 1. ]), array([0., 1., 2.])]
fs: 10.0
rates: [10.01001001 5.00500501]
In all of these cases, it is also possible to pass continuous values instead of spikes (e.g. calcium imaging data), in that case the mean response is computed:
>>> frame = nap.TsdFrame(d=np.random.rand(2000, 3), t=np.arange(0, 100, 0.05))
>>> tcs = nap.compute_tuning_curves(frame, feature, bins=10)
>>> tcs
<xarray.DataArray (unit: 3, 0: 10)> Size: 240B
array([[0.49147343, 0.50190395, 0.50971339, 0.50128013, 0.54332711,
0.49712328, 0.49594611, 0.5110517 , 0.52247351, 0.52057658],
[0.51132036, 0.46410557, 0.47732505, 0.49830908, 0.53523019,
0.53099429, 0.48668499, 0.44198555, 0.49222208, 0.47453398],
[0.46591801, 0.50662914, 0.46875882, 0.48734997, 0.51836574,
0.50722266, 0.48943577, 0.49730095, 0.47944075, 0.48623693]])
Coordinates:
* unit (unit) int64 24B 0 1 2
* 0 (0) float64 80B 0.045 0.135 0.225 0.315 ... 0.585 0.675 0.765 0.855
Attributes:
occupancy: [100. 100. 100. 100. 100. 100. 100. 100. 100. 100.]
bin_edges: [array([0. , 0.09, 0.18, 0.27, 0.36, 0.45, 0.54, 0.63, 0.72,...
Question: Can you compute the response (i.e. firing rate) of the units in tsgroup as function of the feature tsd1 using the function nap.compute_tuning_curves?
tc = nap.compute_tuning_curves(tsgroup, tsd1, bins=5, feature_names=["feat1"])
tc
<xarray.DataArray (unit: 3, feat1: 5)> Size: 120B
array([[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0.08 , 0.08333333, 0. ],
[0.5 , 0.83333333, 1.04 , 1.16666667, 1.375 ]])
Coordinates:
* unit (unit) int64 24B 0 1 2
* feat1 (feat1) float64 40B -2.207 -1.303 -0.3989 0.5049 1.409
Attributes:
occupancy: [ 2. 12. 25. 24. 8.]
bin_edges: [array([-2.65847473, -1.75464293, -0.85081112, 0.05302068, ...
fs: 1.0
rates: [ nan 0.05714286 1.08571429]The output is an xarray object. It is a wrapper of numpy array with extra attributes. It allows to give coordinates to each dimensions as well as attaching attributes. We can make the output look better by labelling the feature we used.
The coordinates can be accessed with the coords attribute. The feature position (i.e. center of the bin) can be accessed with the attribute.
Question: Can you print the underlying the units number, bin center and bin edges of the tuning curve xarray object?
print(tc.unit.values)
print(tc.feat1.values)
print(tc.occupancy)
print(tc.bin_edges)
print(tc.fs)
[0 1 2]
[-2.20655883 -1.30272702 -0.39889522 0.50493658 1.40876839]
[ 2. 12. 25. 24. 8.]
[array([-2.65847473, -1.75464293, -0.85081112, 0.05302068, 0.95685249,
1.86068429])]
1.0
Question: Can you plot the tuning curves for all units?
# tc.plot()
# tc.plot(row="unit")
# tc.plot(col="unit")
# tc[1].plot()
# plt.plot(tc[1].feat1, tc[1].values)
plt.plot(tc.feat1, tc.values.T)
[<matplotlib.lines.Line2D at 0x7fd8e124d460>,
<matplotlib.lines.Line2D at 0x7fd8e18df620>,
<matplotlib.lines.Line2D at 0x7fd8e1741a90>]
Verify Your Setup#
Question: Does the following data download work correctly? If not, please ask a TA.
import workshop_utils
path = workshop_utils.fetch_data("Mouse32-140822.nwb")
print(path)
/home/agent/workspace/rorse_ccn-software-feb-2026_main/data/Mouse32-140822.nwb