---
jupytext:
formats: ipynb,md:myst
text_representation:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.11.5
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
```{code-cell} ipython3
:tags: [remove-cell]
from ndslib.config import jupyter_startup
jupyter_startup()
```
(numpy)=
# The scientific Python ecosystem
Python has an excellent eco-system of scientific computing tools. What do we
mean when we say "scientific computing"? There's no one definition, but we
roughly mean those parts of computing that allow us to work with scientific
data. Data that represents images or time-series from neuroimaging experiments
definitely under that definition. It's also often called "numerical computing",
because these data are often stored as numbers. A common structure for numerical
computing is that of an array. An array is an ordered collection of items that
are stored such that a user can access one or more items in the array. The
[`numpy`](https://numpy.org/) Python library provides an implementation of an
array data structure that is both powerful and light-weight.
What do we mean by "eco-system"? This means that many different individuals and
organizations develop software that interoperates to perform the various
scientific computing tasks that they need to perform. Tools compete for users
and developers' attention and mind-share, but more often they might coordinate.
No one entity manages this coordination. In that sense, the collection of
software tools that exists in Python is more like an organically evolving
ecosystem than a centrally-managed organization. One place in which the
ecosystem comes together as a community is an annual conference called
["Scientific Computing in Python"](https://conference.scipy.org/), which takes
place in Austin, Texas every summer.
```{eval-rst}
.. index::
single: Numpy
```
+++
## Numerical computing in Python
### From naive loops to n-dimensional arrays
In many disciplines, data analysis consists largely of the analysis of tabular
data. This means two-dimensional tables where data are structured into rows
and columns, with each observation typically taking up a row, and each column
representing a single variable. Neuroimaging data analysis also often involves
tabular data, and we'll spend a lot of time in {numref}`pandas` discussing the
tools Python provides for working efficiently with tabular data.
However, one of the notable features of neuroscience experiments is that
datasets are often large, and naturally tend to have many dimensions. This means
that data from neuroscience experiments is often stored in big tables of data
called arrays. Much like Python lists or tuples, the items in an array follow a
particular order. But in contrast to a list or a tuple, they can have multiple
dimensions: arrays can have one dimension (for example a time-series from a
single channel of an EEG recording), two dimensions (for example all the
time-series from multiple EEG channels recorded simultaneously or a single slice
from an MRI image), three dimensions (for example, an anatomical scan of a human
brain), four dimensions (for example, an fMRI dataset with the three spatial
dimensions and time-points arranged as a series of volumes along the last
dimension) or even more.
In all of these cases, the arrays are *continuous*. That means that there are no
holes in an array, where nothing is stored. They are also *homogeneous*, which
means that all of the items in an array are of the same numerical data type.
That is, if one number stored in the array is a floating point number, all
numbers in the array will be represented in their floating point representation.
It also means that the dimensions of the array are predictable. If the first
channel of the EEG time-series has 20,000 time-points in it, the second channel
would also have 20,000 time-points in it. If the first slice of an MRI image has
64-by-64 voxels, other slices in the MRI image would also have 64-by-64 voxels.
![](./figures/multid_arrays.png)
+++
To understand why this is challenging, let's consider a typical fMRI study, in
which participants lie in the scanner, while their brain responses are measured.
Suppose we have 20 participants, each scanned for roughly 30 minutes, with a
repetition time (TR) -— i.e., the duration of acquisition of each fMRI volume -—
of 1 second. If the data are acquired at an isotropic spatial resolution of 2mm
(i.e., each brain "voxel", or 3-dimensional pixel is 2 mm along each dimension),
then the resulting dataset might have approximately 20 x 1800 x 100 x 100 x 100
= 36 billion observations. That's a lot of data! Moreover, each subject's data
has a clear 4 dimensional structure — the 3 spatial dimensions, plus time. If we
wanted to, we could also potentially represent subjects as the 5th dimension,
though that involves some complications, since at least initially, different
subjects' brains won't be aligned with one another —- we'd need to spatially
register them for that (we'll learn more about that in {numref}`registration`).
This means that a single subject's data would look a little bit like the
4-dimensional example that we looked at just a moment ago.
We will typically want to access this data in pretty specific ways. That is,
rather than applying an operation to every single voxel in the brain, at every
single point in time, we usually want to pull out specific slices of the data,
and only apply an operation to those slices. Say for example we're interested in
a voxel in the amygdala. How would we access only that voxel, at every time
point?
A very naive approach that we could implement in pure Python would be to store
all our data as a nested series of Python lists: each element in the first list
would be a list containing data for one time point; each element within the list
for each time point would itself be a list containing the 2d slices at each
x-coordinate; and so on. Then we could write a series of 4 nested for-loops to
sequentially access every data point (or voxel) in our array. For each voxel we
inspect, we could then determine whether the voxel is one we want to work with,
and if so, apply some operation to it.
Basically, we'd have something like this (note that this is just pseudocode, not
valid Python code, and you wouldn't be able to execute this snippet!):
```python
for t in time:
for x in t:
for y in x:
for z in y:
if z is in amygdala:
apply_my_function(z)
```
Something like the above would probably work just fine. But it would be very
inefficient, and might take a while to execute. Even if we only want to access a
single amygdala voxel at each timepoint (i.e., 1800 data points in total, if we
stick with our resting state example from above), our naive looping approach
still requires us to inspect every single one of the 1.8 billion voxels in our
dataset —- even if all we're doing for the vast majority of voxels is deciding
that we don't actually need to use that voxel in our analysis! This is, to put
it mildly, A Bad Thing.
Computer scientists have developed a notation for more precisely describing how
much of a bad thing it is: Big O notation, which gives us a formal way of
representing the time complexity of an algorithm. In Big O notation, the code
above runs in $O(n^4)$ time. Meaning, for any given dimension size $n$ (we
assume for simplicity that $n$ is the same for all dimensions), we need to
perform on the order of $n^4$ computations. To put this in perspective, just
going from $n = 2$ to $n = 4$ means that we go from $2^4 = 16$ to $4^4 = 256$
operations. You can see how this might start to become problematic as $n$ gets
large.
The lesson we should take away from this is that writing naive Python loops just
to access the values at a specific index in our dataset is a bad idea.
### Towards a specialized array structure
What about list indexing? Fortunately, since we've already covered list indexing
in Python, we have a better way to access the amygdala data we want. We know
that the outer dimension in our nested list represents time, and we also know
that the spatial index of the amygdala voxel we want will be constant across
time points. This means that instead of writing a 4-level nested for-loop, we
can do the job in just a single loop. If we know that the index of the amygdala
voxel we want in our 3-d image is [20, 18, 32], then the following code suffices
(again note that this code can't be executed):
```python
amygdala_values = []
for t in dataset:
amygdala_values.append(t[20][18][32])
```
That's a big improvement, in terms of time complexity. We've gone from an
$O(n^4)$ to an $O(n)$ algorithm. The time our new algorithm takes to run is now
linear in the number of timepoints in our dataset, rather than growing as a
polynomial function of it.
Unfortunately, our new approach still isn't great. One problem is that
performance still won't be so good, because Python isn't naturally optimized for
the kind of operation we just saw -— indexing repeatedly into nested lists. For
one thing, Python is a dynamically typed language, so there will be some
overhead involved in accessing each individual element (because the interpreter
has to examine the object and determine its type), and this can quickly add up
if we have millions or billions of elements.
The other big limitation is that, even if lists themselves were super fast to
work with, the Python standard library provides limited functionality for doing
numerical analysis on them. There are a lot of basic numerical operations we
might want to apply to our data that would be really annoying to write using
only the built-in functionality of Python. For example, Python lists don't
natively support matrix operations, so if we wanted to multiply one matrix
(represented as a list of lists) by another, we'd probably need to write out a
series of summations and multiplications ourselves inside for-loops. And that
would be completely impractical given how basic and common an operation matrix
multiplication is (not to important for our point now, but if you're wondering
what a matrix multiplication is, you can see a formal definition in
{numref}`matmul`).
The bottom line is that, as soon as we start working with large, or highly
structured datasets —- as is true of most datasets in neuroimaging —- we're
going to have to look beyond what Python provides us in its standard library.
Lists just aren't going to cut it; we need some other kind of data structure —-
one that's optimized for numerical analysis on large, N-dimensional datasets.
This is where numpy comes into play.
## Introducing `Numpy`
Numpy is the backbone of the numerical and scientific computing stack in Python;
many of the libraries we'll cover in this book (pandas, scikit-learn, etc.)
depend on it internally. Numpy provides many data structures optimized for
efficient representation and manipulation of different kinds of high-dimensional
data, as well as an enormous range of numerical tools that help us work with
those structures.
### Importing numpy
Recall that the default Python namespace contains only a small number of
built-in functions. To use any other functionality, we need to explicitly import
it into our namespace. Let's do that for Numpy. The library can be imported
using the `import` key word. We use `import numpy as np` to create a short name
for the library that we can refer to. Everything that is in `numpy` will now be
accessible through the short name `np`. Importing the library specifically as
`np` is not a requirement, but it is a very strongly-held convention and when
you read code that others have written, you will find that this the form that is
very often used.
```{code-cell} ipython3
import numpy as np
```
### Data is represented in arrays
The core data structure in Numpy is the n-dimensional array (or `ndarray`). As
the name suggests, an `ndarray` is an array with an arbitrary number of
dimensions. Unlike Python lists, Numpy arrays are homogeneously typed -—
meaning, every element in the array has to have the same data type. You can have
an array of floats, or an array of integers, but you can't have an array that
mixes floats and integers. This is exactly what we described above: a
homogeneous dense table that holds items in it.
### Creating `ndarray`s
Like any other Python object, we need to initialize an `ndarray` before we can
do anything with it. Numpy provides us with several ways to create new arrays.
Let's explore a couple.
#### Initializing an array from an existing list
Let's start by constructing an array from existing data. Assume we have some
values already stored in a native Python iterable object (typically, a list;
"iterable" means that it is an object that we can iterate over with a for loop,
see {numref}`forloop`). If we want to convert that object to an ndarray so that
we can perform more efficient numerical operations on it, we can pass the object
directly to the `np.narray()` function. In a relatively simple case, the input
to the `np.array` function is a sequence of items, a list that holds some
integers.
```{code-cell} ipython3
list_of_numbers = [1, 1, 2, 3, 5]
arr1 = np.array(list_of_numbers)
print(arr1)
```
The result still looks a lot like a list of numbers, but let's look at some of
the attributes that this variable has that are not available in lists.
```{code-cell} ipython3
print("The shape of the array is", arr1.shape)
print("The dtype of the array is", arr1.dtype)
print("The size of each item is", arr1.itemsize)
print("The strides of the array is", arr1.strides)
```
The shape of an array is its dimensionality. In this case, the array is
one-dimensional (like the single channel of EEG data in the example above). All
of the items in this array are integers, so we can use an integer representation
as the dtype of the array. This means that each element of the array takes up 8
bytes of memory. Using 64-bit integers, we can represent the numbers from
-9223372036854775808 to 9223372036854775807. And each one would take up 8 bytes
of memory in our computer (which are 64 bits). This also explains the strides.
Strides are the number of bytes in memory that we would have to skip to get to
the next value in the array. Because the array is packed densely in a contiguous
segment of memory, to get from the value `1` to the value `2`, for example, we
have to skip 8 bytes.
![](./figures/numpy_arrays1.png)
+++
Remember that the data we would like to represent in arrays can have multiple
dimensions. Numpy is pretty good about figuring out what kind of array we want
based on the input data; for example, if we pass in a list of lists, where each
of the inner lists has the same number of elements, Numpy will infer that we
want to construct a 2-dimensional array
```{code-cell} ipython3
list_of_lists = [[1, 1, 2, 3, 5], [8, 13, 21, 34, 55]]
arr2 = np.array(list_of_lists)
print(arr2)
print("The shape of this array is", arr2.shape)
print("The dtype of this array is", arr2.dtype)
print("The size of each item is", arr2.itemsize)
print("The strides of this array are", arr2.strides)
```
Now, the shape of the array is 2 items: the first is the number of lists in the
array (2) and the second is the number of items in each list (5). You can also
think of that as rows (2 rows) and columns (5 columns), because each item in the
first row of the array has an item matching it in terms of its column in the
second row. For example, the number `2` in the first row is equivalent to the
number `21` in the second row in that they are both the third item in their
particular rows. Or, in other words, they share the third column in the array.
The items are still each 8 bytes in size, but now there are two items in the
strides: `(40, 8)`. The first item tells us how many bytes we have to move to
get from one row to the next row (for example, from the `2` in the first row to
the `21` in the second row). That's because the part of the computer's memory in
which the entire array is stored is contiguous and linear. Here, we've drawn
this contiguous bit of memory, such that the part that represents the first row
of the array is green and the part that represents the second row of the array
is blue.
![](./figures/numpy_arrays2.png)
+++
The first time you see a printed representation of a Numpy array, as in the
above output, it might look a little confusing. But the visual representation
closely resembles what we'd see if we were to print the nested list in Python.
In the 2-dimensional case above, the elements of the outer list are the values
along the first dimension, and each inner list gives the values along the second
dimension. This ends up giving us a nice tabular representation, where, at least
for small arrays, we can just read off the values from the output. For example,
the first row of the array contains the values [1, 1, 2, 3, 5]; the second
column of the array contains the values [1, 13]; and so on.
(numpy_ex1)=
#### Exercise
Extend the above principle into 3 dimensions: create a list of list of lists,
initialize an array from the list, and print into the screen. Make sure you
understand how the values displayed visually map onto the dimensions of the
array.
#### Initializing an empty array
```{eval-rst}
.. index::
single: zeros
```
Alternatively, we can construct a new array from scratch and fill it with some
predetermined value -- most commonly zero. We can do this easily in Numpy using
the conveniently - named `np.zeros` function. In many practical applications, we
can think of this as an "empty" array (though technically we could also create a
truly empty array that has no assigned values using the `np.empty` function).
The `zeros` function takes a mandatory shape tuple as its first argument; this
specifies the dimensions of the desired array.
```{code-cell} ipython3
arr_2d = np.zeros((5, 10))
print("The shape of this array is", arr_2d.shape)
print("The dtype of this array is", arr_2d.dtype)
print("The size of each item is", arr_2d.itemsize)
print("The strides of this array are", arr_2d.strides)
arr_3d = np.zeros((2, 4, 8))
print("The shape of this array is", arr_3d.shape)
print("The dtype of this array is", arr_3d.dtype)
print("The size of each item is", arr_3d.itemsize)
print("The strides of this array are", arr_3d.strides)
```
(numpy_ex2)=
#### Exercise
Can you explain the outputs of the previous code cell? What does the first item in the `strides` represent in each case?
### Neuroscience data is stored in arrays
```{eval-rst}
.. index::
single: Strides
```
As we described above, many different kinds of neuroscience data are naturally
stored in arrays. To demonstrate this, let's look at some neuroscience data that
we've stored in an array. We've stored some data from an fMRI experiment in a
Numpy array, you can load it using our the `load_data` utility function from the
`nsdlib` software library that we have implemented:
```{code-cell} ipython3
from ndslib import load_data
```
This `load_data` function takes as arguments the name of the dataset that we
would like to download and alternatively also the name of a file that we would
like to save the data into (for reuse). The `npy` file format is a format that
was developed specifically to store Numpy arrays and that's all that it knows
how to store. When you load one of these files, you can assign into a variable
the contents of the array that was stored in the file.
```{code-cell} ipython3
bold = load_data("bold_numpy", fname="bold.npy")
print("The shape of the data is", bold.shape)
print("The dtype of the data is", bold.dtype)
print("The itemsize of the data is", bold.itemsize)
print("The strides of the data is", bold.strides)
```
This data is a little bit more complicated than the arrays you have seen so far.
It has 4 dimensions: the first two dimensions are the inplane dimensions of each
slice: 64 by 64 voxels. The next dimension is the slice dimension: 25 slices.
Finally, the last dimension is the time-points that were measured in the
experiment. The dtype of this array is `float64`. These floating point numbers
are what computers use to represent real numbers. Each one of these items also
uses up 64 bits, or 8 bytes, so that's the item size here as well. Just like the
shape, there are 4 strides: the first item in the strides tells how far we need
to go to get from one voxel to another voxel in the same row, in the same slice,
in the same time-point. The second stride (512) tells us how far we would have
to go to find the same voxel in the same slice and time-point, but one row over.
That's 64 times 8. The next one is about finding the same voxel in a neighboring
slice. Finally the large one at the end is about finding the same slice in the
next time step. That large number is equal to 64 by 64 by 25 (the number of
voxels in each time point) by 8 (the number of bytes per element in the array).
+++
### Indexing is used to access elements of an array
```{eval-rst}
.. index::
single: Indexing (arrays)
```
We've seen how we can create arrays and describe them; now let's talk about how
we can get data in and out of arrays. We already know how to index Python lists,
and array indexing in Python will look quite similar. But Numpy indexing adds
considerably more flexibility and power, and developing array indexing facility
is a critical step towards acquiring general proficiency with the package. Let's
start simple: accessing individual items in a one-dimensional array looks a lot
like indexing into a list (see {numref}`listindexing`). For example, the first
item in the one-dimensional array we saw earlier, you would do the following:
```{code-cell} ipython3
arr1[0]
```
or if you want the second element, you would do something like this (don't forget that Python uses 0-based indexing!):
```{code-cell} ipython3
arr1[1]
```
Numpy arrays also support two other important syntactic conventions that are
also found in Python lists: indexing from the end of the array, and "slicing"
the array to extract multiple contiguous values.
To index from the end, we use the minus sign (-):
```{code-cell} ipython3
arr1[-2]
```
This gives us the 2nd-from-last value in the array. To slice an array, we use
the colon (:) operator, passing in the positions we want to start and end at:
```{eval-rst}
.. index::
single: Slicing (arrays)
```
```{code-cell} ipython3
arr1[2:5]
```
As in lists, the start index is inclusive and the end index is exclusive (i.e.,
in the above example, the resulting array includes the value at position 2, but
excludes the one at position 6). We can also omit the start or end indexes, in
which case, Numpy will return all positions up to, or starting from, the provided
index. For example, to get the first 4 elements:
```{code-cell} ipython3
arr1[:4]
```
### Indexing in multiple dimensions
Once we start working with arrays with more than one dimension, indexing gets a
little more complicated. Both because the syntax is a little different, and
because there's a lot more we can do with multi-dimensional arrays.
For example, if you try one of these operations with a two-dimensional array,
something funny will happen:
```{code-cell} ipython3
arr2[0]
```
It looks like the first item of a two-dimensional array is the entire first row
of that array. To access the first element in the first row, we need to use two
indices, separated by a comma:
```{code-cell} ipython3
arr2[0, 0]
```
One way to think about this is that the first index points to the row and the
second index points to the column, so to get the very first item in the array,
we need to explicitly ask for the first column in the first row. Consider what
happens when you ask for third item in the second row:
```{code-cell} ipython3
arr2[1, 2]
```
Under the hood, Python asks how many bytes into the array it would have to go to
get that item. It looks at the strides and computes: 1 times the first stride =
40 plus 2 times the second stride = 40 + 16 = 56. Indeed, the third element of
the second row of the array is 7 times 8 bytes into the array:
![](./figures/numpy_arrays3.png)
+++
How would you access data in the BOLD array? As an example, let's access the
voxel in the very first time-point, in the very first slice, and in this slice
in the first row and column:
```{code-cell} ipython3
bold[0, 0, 0, 0]
```
This element of the array is the very corner of the volume, so it is way outside
of the part of the data that contains the brain. That is why it contains a value
of 0. Instead, we can look at the voxel at the very center of the center slice:
```{code-cell} ipython3
bold[32, 32, 12, 0]
```
This voxel has a value that corresponds to the intensity of the MR image in this
voxel at the first time point. If we drop the last index, that would be like
asking for the data in all of the time-points in that particular voxel:
```{code-cell} ipython3
bold[32, 32, 12]
```
This is now a time-series of data, corresponding to the intensity of the BOLD
response in the voxel at the very center of the central slice, across all time
points in the measurement.
+++
### Slicing a dimension with "`:`"
As you saw before, we can use slicing to choose a sub-array from within our
array. For example, you can slice out the time-series in the central voxel, as
we just did:
```{code-cell} ipython3
bold[32, 32, 12, :]
```
This is equivalent to dropping the last index. If we want to slice through all
of the slices, taking the values of the central voxel in each slice, in the very
first time-point, we would instead slice on the third (slice) dimension of the
array:
```{code-cell} ipython3
bold[32, 32, :, 0]
```
These 25 values correspond to the values in the center of each slice in the
first time point. We can also slice on multiple dimensions of the array. For
example, if we would like to get the values of all the voxels in the central
slice in the very first time-point, we can do:
```{code-cell} ipython3
central_slice_t0 = bold[:, :, 12, 0]
```
This is a two-dimensional array, with the dimensions of the slice:
```{code-cell} ipython3
central_slice_t0.shape
```
Which we can verify by printing out the values inside the array:
```{code-cell} ipython3
print(central_slice_t0)
```
```{code-cell} ipython3
bold.shape
```
(boolean_indexing_numpy)=
### Indexing with conditionals
One more way of indexing is to to use logical operations, or conditionals. This
allows us to choose values from an array, only if they fulfill certain
conditions. For example, if we'd like to exclude measurements with a value of 0
or less, we might perform a logical operation like the following:
```{code-cell} ipython3
larger_than_0 = bold>0
```
This creates a new array that has the same shape as the original array, but
contains Boolean (`True`/`False`) values. Wherever we had a number larger than 0
in the original array, this new array will have the value `True` and wherever
the numbers in the original array were 0 or smaller, this new array will have
the value `False`. This is already helpful in some regards, but the next thing
that we can do with this is that we can use this array to index into the
original array. What this does is that it pulls out of the original array only
the values for which `larger_than_0` is `True` and reorganizes these values as a
long one-dimensional array:
```{code-cell} ipython3
bold[larger_than_0]
```
This can be useful, for example, if you would like to calculate statistics of
the original array, but only above some threshold value. For example, if you
would like to calculate the average of the array (using the `mean` method of the
array) with and without values that are 0 or less:
```{code-cell} ipython3
print(bold.mean())
print(bold[larger_than_0].mean())
```
(numpy_ex3)=
#### Exercise
Arrays have a variety of methods attached to them to do calculations. Use the
array method `std` to calculate the standard deviation of the `bold` array, and
the standard deviation of just values of the array that are between 1 and 100.
### Arithmetic with arrays
In addition to storing data from a measurement, and providing access to the this
data through indexing, arrays also support a variety of mathematical operations
with arrays. For example, you can perform arithmetic operations between an array
and a number. Consider for example what happens when we add a single number to
an entire array:
```{code-cell} ipython3
bold[32, 32, 12] + 1000
```
The result of this operation is that the number is added to every item in the
array. Similarly:
```{code-cell} ipython3
bold[32, 32, 12] / 2
```
In this case, the number is used to separately divide each one of the elements
of the array.
We can also do arithmetic between arrays. For example, we can add the values in
this voxel to the items in the central voxel in the neighboring slice:
```{code-cell} ipython3
bold[32, 32, 12] + bold[32, 32, 11]
```
```{eval-rst}
.. index::
single: Array-scalar arithmetic
```
```{eval-rst}
.. index::
single: Array-array arithmetic
```
When we perform arithmetic operations between an array and a number (or
array-scalar operations), the number is separately added to each element. Here,
when we perform arithmetic operations between arrays (array-array operations),
the operation is instead performed element-by-element. In both cases the result
is an array that has the same number of items as the original array (or arrays).
```{code-cell} ipython3
tsnr = bold.mean(-1) / bold.std(-1)
print(tsnr[tsnr.shape[0]//2, tsnr.shape[1]//2, tsnr.shape[2]//2])
```
(numpy_ex4)=
#### Exercises
1. Creating an array: create an array that contains the sequence of numbers from 1 to 99 in steps of 2: [1, 3, 5, ..., 99]
2. Creating an array of 1's: create an array of shape (3, 5) containing only the number 1
3. Find all the items in an array containing numbers that have values that are larger than 0 and smaller or equal to 100
4. Find all the items in an array that can be divided by 3 with no remainder (hint: the function `np.mod` executes the modulo operation on every item in the array).
5. Generate an array with 100 randomly distributed numbers between 0 and 10 (hint: look at the `np.random` sub-module).
### The Scipy library
The Numpy array serves as the basis for numerical computing in Python. But it
does only so much. Scientific computing involves a variety of different
algorithms that do more interesting things with numbers. For example, we'd like
to use fancy image processing and signal processing algorithms with our data. Or
rely on various facts about statistical distributions. Over the last 100 years,
scientists and engineers have developed many different algorithms for scientific
computing. In the Python scientific computing ecosystem, many of the fundamental
algorithms are collected together in one library, called ``Scipy``. It provides
everything but the kitchen sink, in terms of basic algorithms for science.
Because it's such a large collection of things, it is composed of several
different sub-libraries that perform various scientific computing tasks, such as
statistics, image processing, clustering, linear algebra, optimization, etc. We
will revisit some of these in later chapters. For now, as an example, we will
focus on the signal processing sub-library ``scipy.signal``. We import the
sub-library and name it `sps`:
```{code-cell} ipython3
import scipy.signal as sps
```
Remember that the BOLD dataset that we were looking at before had 180
time-points. That means that the shape of this data was `(64, 64, 25, 180)`. One
basic signal processing operation is to efficiently and accurately resample a
time-series into a different sampling rate. For example, let's imagine that we'd
like to resample the data to a slower sampling rate. The original data was
collected with a sampling interval or repetition time (TR) of 2.0 seconds. If we
would like to resample this to half the sampling rate (TR=4.0 seconds), we would
halve the number of samples on the last dimension using the `sps.resample`
function:
```{code-cell} ipython3
resampled_bold = sps.resample(bold, 90, axis=3)
print(resampled_bold.shape)
```
If, instead, we'd like to resample this to twice the sampling rate (TR=1.0), we
would double the number of samples.
```{code-cell} ipython3
resampled_bold = sps.resample(bold, 360, axis=3)
print(resampled_bold.shape)
```
Although an entire book could be written about Scipy, we will not go into much more detail here. One reason is that Scipy is *under the hood* of many of the applications that you will see later in the book. For example, the "high-level" libraries that we will use for image-processing ({numref}`skimage`) and for machine-learning ({numref}`sklearn`) rely on implementations of fundamental computational methods in Scipy. But, unless you are building such tools from scractch, you may not need to use it directly. Nevertheless, it is a wonderful resource to know about, when you do want to build something new that is not already implemented in the existing libraries.
(numpy_ex5)=
#### Exercise
Use the `scipy.interpolate` sub-library to interpolate the signal in one voxel (take the central voxel, `bold[32, 32, 12]`) to even higher resolution (TR=0.5 seconds, or 720 samples) using a cubic spline. Hint: Look at the documentation of the `scipy.interpolate.interp1d` function to choose the method of interpolation and for example usage.
### Additional resources
For a description of the scientific Python ecosystem, see Fernando Perez, Brian
Granger and John Hunter's paper from 2011 {cite}`Perez2010-nn`. This paper
describes in some detail how Python has grown from a general-purpose programming
language into a flexible and powerful tool that covers many of the uses that
scientists have, through a distributed, eco-system of development. More recent
papers describe the tools at the core of the ecosystem: The Numpy library is
described in one recent paper {cite}`Harris2020-so` and the Scipy library is
described in another {cite}`Virtanen2020-jl`.