%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.
Learning the fundamentals of pynapple#
Learning objectives#
Instantiate the pynapple objects
Make the pynapple objects interact
Use numpy with pynapple
Slicing pynapple objects
Adding metadata to pynapple objects
Learn the core functions of pynapple
The pynapple documentation can be found here.
The documentation for objects and method of the core of pynapple is here.
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:2025-11-15 00:46:51,302:jax._src.xla_bridge:850: An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
/home/jenkins/workspace/rorse_ccn-software-sfn-2025_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.21918
1.0 0.0533441
2.0 1.88599
3.0 0.0484334
4.0 0.474418
5.0 0.298899
6.0 -1.22407
...
93.0 0.906162
94.0 0.74899
95.0 -0.193397
96.0 1.08983
97.0 0.132593
98.0 -0.864925
99.0 0.830074
dtype: float64, shape: (100,)
Question: Can you print tsd2?
print(tsd2)
Time (s) pineapple banana tomato
---------- ----------- ---------- ---------
0.0 -2.35606 0.175694 -1.55505
1.0 0.927279 -1.17221 0.744864
2.0 -1.61942 0.593622 -1.20612
3.0 -0.524923 -0.0822806 0.680098
4.0 2.05214 -0.301787 -0.22075
5.0 -0.32219 -0.339772 -0.854587
6.0 -0.919266 -1.48068 1.29426
...
93.0 -1.48166 0.262022 -0.269566
94.0 0.340472 -0.822813 0.347829
95.0 1.26923 0.914358 0.508982
96.0 0.528647 0.155076 0.993064
97.0 -1.70918 -0.297416 1.36131
98.0 -0.950789 0.935986 0.517214
99.0 -0.0922875 -0.524276 0.477058
dtype: float64, shape: (100, 3)
Question: Can you print tsd3?
print(tsd3)
Time (s)
---------- -------------------------------
0.0 [[ 0.971225 ... -0.398633] ...]
0.1 [[-0.278583 ... -1.126599] ...]
0.2 [[ 0.796699 ... -0.322272] ...]
0.3 [[-1.619002 ... -0.142493] ...]
0.4 [[0.271323 ... 1.81719 ] ...]
0.5 [[1.381239 ... 1.098416] ...]
0.6 [[ 0.275037 ... -2.423651] ...]
...
99.3 [[-0.325148 ... 1.649048] ...]
99.4 [[ 0.694916 ... -0.521061] ...]
99.5 [[-0.440783 ... -0.615822] ...]
99.6 [[ 0.491117 ... -1.754317] ...]
99.7 [[-1.049159 ... -0.491122] ...]
99.8 [[-1.248918 ... 0.255819] ...]
99.9 [[ 0.671458 ... -1.662642] ...]
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: … and print it?
print(tsgroup)
Index rate
------- -------
0 0.0391
1 0.05213
2 1.30321
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.8736
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 realized the 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.186128
11.0 -0.463043
12.0 1.49725
13.0 1.1733
14.0 -1.50306
15.0 -1.04424
16.0 0.0074228
...
93.0 0.906162
94.0 0.74899
95.0 -0.193397
96.0 1.08983
97.0 0.132593
98.0 -0.864925
99.0 0.830074
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 5.68007 7.4132
1 23.2985 34.6245
2 35.3759 36.7646
3 38.7375 40.9948
4 43.5923 46.9667
5 47.7217 51.8919
6 61.3927 75.7028
7 77.1481 77.9268
8 81.9966 88.8798
9 97.4802 99.6025
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 23.2985 30
1 50 51.8919
2 61.3927 75.7028
3 77.1481 77.9268
4 81.9966 88.8798
5 97.4802 99.6025
shape: (6, 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 5.68007 7.4132
1 10 34.6245
2 35.3759 36.7646
3 38.7375 40.9948
4 43.5923 46.9667
5 47.7217 100
shape: (6, 2), time unit: sec.
Question: … and 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 10 23.2985
1 51.8919 61.3927
2 75.7028 77.1481
3 77.9268 81.9966
4 88.8798 97.4802
5 99.6025 100
shape: (6, 2), time unit: sec.
Question: … and visualize it?
workshop_utils.visualize_intervals([ep_signal, ep_tmp, ep_signal.set_diff(ep_tmp)])
Numpy & pynapple#
Pynapple objects behaves very similarly like numpy array. 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.01460362 -0.04258991 0.02104224 0.00705391 0.02599554]
[ 0.06060399 -0.03064132 -0.03175223 -0.01278073 -0.04610488]
[ 0.00454288 0.00494842 0.00160614 0.03519576 -0.03957175]
[-0.02638122 0.01718459 -0.04386483 -0.02277178 -0.03056191]]
Question:: can you compute the average of tsd2 for each timestamps and print it?
print(np.mean(tsd2, 1))
Time (s)
---------- ----------
0.0 -1.24514
1.0 0.166644
2.0 -0.743973
3.0 0.0242983
4.0 0.509866
5.0 -0.505516
6.0 -0.368563
...
93.0 -0.496402
94.0 -0.0448372
95.0 0.897523
96.0 0.558929
97.0 -0.215098
98.0 0.16747
99.0 -0.0465017
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))
[[-0.34010665 0.29079866 1.1586195 0.75325926 0.04871718]
[-1.97931525 -1.14318835 -1.10872367 0.39001127 0.74280212]
[-0.19467827 0.1890571 -0.29143654 -0.6466056 1.60366266]
[-0.11479806 0.36815689 0.43145111 -0.31955783 -0.77646658]]
Metadata#
Metadata 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: … and print it?
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 -2.3560559242662644 0.17569368934995985 -1.555054713589844
1.0 0.9272787573343956 -1.1722099327734883 0.7448637201623065
2.0 -1.6194185881168184 0.59362177532366 -1.2061207493148416
3.0 -0.5249226622352758 -0.08228060704446603 0.6800980473905671
4.0 2.052135344698557 -0.30178664653559983 -0.22074950780535216
5.0 -0.3221903427905951 -0.33977161213798257 -0.8545865047730654
6.0 -0.919265764026655 -1.480683083466595 1.2942598791123627
...
93.0 -1.4816623279170766 0.2620215598420029 -0.26956564660703425
94.0 0.3404723746983879 -0.8228128654944309 0.3478288514034553
95.0 1.2692286189432118 0.9143576299931303 0.5089819122262199
96.0 0.5286473676864707 0.15507615282573572 0.9930642290412443
97.0 -1.7091828649632772 -0.2974163365341157 1.3613055649131436
98.0 -0.9507887701492719 0.9359856769884788 0.5172140443044121
99.0 -0.0922875451401194 -0.524275897995294 0.47705834269864145
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 dictionnary 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.30321 saturn
Metadata are accessible either as attributes (i.e. tsgroup.planet) or as dictionnary-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.6454909024717835 -0.4251212613642054 0.4759214732648534
1.0 -1.0164884865935153 0.7704257647350657 -0.6619299642284323
2.0 -0.9630138087757621 0.12036774271438518 0.11849319405464497
3.0 -0.6269131275749616 -0.22969203794798787 0.8988949281843825
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.6454909 -0.42512126 0.47592147]
Question: What happen when you do tsdframe.loc[0] and tsdframe.loc[[0,1]]
print(tsdframe.loc[0])
print(tsdframe.loc[[0,1]])
Time (s)
---------- ---------
0 -0.425121
1 0.770426
2 0.120368
3 -0.229692
dtype: float64, shape: (4,)
Time (s) 0 1
---------- -------------------- -------------------
0.0 -0.4251212613642054 0.4759214732648534
1.0 0.7704257647350657 -0.6619299642284323
2.0 0.12036774271438518 0.11849319405464497
3.0 -0.22969203794798787 0.8988949281843825
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.6454909024717835
1.0 -1.0164884865935153
2.0 -0.9630138087757621
3.0 -0.6269131275749616
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 4
11.5 0 0 1
12.5 0 0 2
13.5 0 0 0
14.5 0 0 0
15.5 0 0 0
16.5 0 0 2
...
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 0x7fbc470e4aa0>]
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 0x7fbc44fcee10>]
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 0x7fbc44fb5d60>]
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 11.5 13.5
1 15.5 17.5
2 19.5 21.5
3 22.5 23.5
4 24.5 25.5
5 28.5 30.0
6 50.5 51.5
... ... ...
14 76.5 79.5
15 80.5 81.5
16 85.5 86.5
17 88.5 89.5
18 92.5 94.5
19 95.5 97.5
20 98.5 99.0
shape: (21, 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 0x7fbc44e6e420>,
<matplotlib.patches.Rectangle at 0x7fbc4734abd0>,
<matplotlib.patches.Rectangle at 0x7fbc4747fad0>,
<matplotlib.patches.Rectangle at 0x7fbc473f52e0>,
<matplotlib.patches.Rectangle at 0x7fbc470e5430>,
<matplotlib.patches.Rectangle at 0x7fbc44b31250>,
<matplotlib.patches.Rectangle at 0x7fbc44afea50>,
<matplotlib.patches.Rectangle at 0x7fbc44cea450>,
<matplotlib.patches.Rectangle at 0x7fbc44ca7dd0>,
<matplotlib.patches.Rectangle at 0x7fbc44cc3020>,
<matplotlib.patches.Rectangle at 0x7fbc44c92e10>,
<matplotlib.patches.Rectangle at 0x7fbc44e15850>,
<matplotlib.patches.Rectangle at 0x7fbc477bca40>,
<matplotlib.patches.Rectangle at 0x7fbc44b13d40>,
<matplotlib.patches.Rectangle at 0x7fbc44ec63c0>,
<matplotlib.patches.Rectangle at 0x7fbc44fcc710>,
<matplotlib.patches.Rectangle at 0x7fbc44d00aa0>,
<matplotlib.patches.Rectangle at 0x7fbc473020c0>,
<matplotlib.patches.Rectangle at 0x7fbc44ea91c0>,
<matplotlib.patches.Rectangle at 0x7fbc473428a0>,
<matplotlib.patches.Rectangle at 0x7fbc44ea9e20>]
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,...
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. ...
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.])]
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.04347826, 0.08695652, 0. , 0.33333333],
[0.66666667, 0.86956522, 1.26086957, 1.1875 , 2. ]])
Coordinates:
* unit (unit) int64 24B 0 1 2
* feat1 (feat1) float64 40B -1.779 -0.7877 0.2038 1.195 2.187
Attributes:
occupancy: [ 6. 23. 23. 16. 3.]
bin_edges: [array([-2.27493437, -1.28342894, -0.29192352, 0.69958191, ...
fs: 1.0The 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]
[-1.77918165 -0.78767623 0.20382919 1.19533462 2.18684004]
[ 6. 23. 23. 16. 3.]
[array([-2.27493437, -1.28342894, -0.29192352, 0.69958191, 1.69108733,
2.68259275])]
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 0x7fbc449a0470>,
<matplotlib.lines.Line2D at 0x7fbc449a04d0>,
<matplotlib.lines.Line2D at 0x7fbc449a05c0>]
Important#
Question: Does this work? If not, please ask a TA.
import workshop_utils
path = workshop_utils.fetch_data("Mouse32-140822.nwb")
print(path)
/home/jenkins/workspace/rorse_ccn-software-sfn-2025_main/data/Mouse32-140822.nwb