{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "9c086566",
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"from ndslib.config import jupyter_startup\n",
"jupyter_startup()"
]
},
{
"cell_type": "markdown",
"id": "32455d1b",
"metadata": {},
"source": [
"(numpy)=\n",
"# The scientific Python ecosystem\n",
"\n",
"Python has an excellent eco-system of scientific computing tools. What do we\n",
"mean when we say \"scientific computing\"? There's no one definition, but we\n",
"roughly mean those parts of computing that allow us to work with scientific\n",
"data. Data that represents images or time-series from neuroimaging experiments\n",
"definitely under that definition. It's also often called \"numerical computing\",\n",
"because these data are often stored as numbers. A common structure for numerical\n",
"computing is that of an array. An array is an ordered collection of items that\n",
"are stored such that a user can access one or more items in the array. The\n",
"[`numpy`](https://numpy.org/) Python library provides an implementation of an\n",
"array data structure that is both powerful and light-weight.\n",
"\n",
"What do we mean by \"eco-system\"? This means that many different individuals and\n",
"organizations develop software that interoperates to perform the various\n",
"scientific computing tasks that they need to perform. Tools compete for users\n",
"and developers' attention and mind-share, but more often they might coordinate.\n",
"No one entity manages this coordination. In that sense, the collection of\n",
"software tools that exists in Python is more like an organically evolving\n",
"ecosystem than a centrally-managed organization. One place in which the\n",
"ecosystem comes together as a community is an annual conference called\n",
"[\"Scientific Computing in Python\"](https://conference.scipy.org/), which takes\n",
"place in Austin, Texas every summer.\n",
"\n",
"```{eval-rst}\n",
".. index::\n",
" single: Numpy\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "d1ba0a5a",
"metadata": {},
"source": [
"## Numerical computing in Python\n",
"\n",
"### From naive loops to n-dimensional arrays\n",
"\n",
"In many disciplines, data analysis consists largely of the analysis of tabular\n",
"data. This means two-dimensional tables where data are structured into rows\n",
"and columns, with each observation typically taking up a row, and each column\n",
"representing a single variable. Neuroimaging data analysis also often involves\n",
"tabular data, and we'll spend a lot of time in {numref}`pandas` discussing the\n",
"tools Python provides for working efficiently with tabular data.\n",
"\n",
"However, one of the notable features of neuroscience experiments is that\n",
"datasets are often large, and naturally tend to have many dimensions. This means\n",
"that data from neuroscience experiments is often stored in big tables of data\n",
"called arrays. Much like Python lists or tuples, the items in an array follow a\n",
"particular order. But in contrast to a list or a tuple, they can have multiple\n",
"dimensions: arrays can have one dimension (for example a time-series from a\n",
"single channel of an EEG recording), two dimensions (for example all the\n",
"time-series from multiple EEG channels recorded simultaneously or a single slice\n",
"from an MRI image), three dimensions (for example, an anatomical scan of a human\n",
"brain), four dimensions (for example, an fMRI dataset with the three spatial\n",
"dimensions and time-points arranged as a series of volumes along the last\n",
"dimension) or even more.\n",
"\n",
"In all of these cases, the arrays are *continuous*. That means that there are no\n",
"holes in an array, where nothing is stored. They are also *homogeneous*, which\n",
"means that all of the items in an array are of the same numerical data type.\n",
"That is, if one number stored in the array is a floating point number, all\n",
"numbers in the array will be represented in their floating point representation.\n",
"It also means that the dimensions of the array are predictable. If the first\n",
"channel of the EEG time-series has 20,000 time-points in it, the second channel\n",
"would also have 20,000 time-points in it. If the first slice of an MRI image has\n",
"64-by-64 voxels, other slices in the MRI image would also have 64-by-64 voxels.\n",
"\n",
"![](./figures/multid_arrays.png)"
]
},
{
"cell_type": "markdown",
"id": "a00cba56",
"metadata": {},
"source": [
"To understand why this is challenging, let's consider a typical fMRI study, in\n",
"which participants lie in the scanner, while their brain responses are measured.\n",
"Suppose we have 20 participants, each scanned for roughly 30 minutes, with a\n",
"repetition time (TR) -— i.e., the duration of acquisition of each fMRI volume -—\n",
"of 1 second. If the data are acquired at an isotropic spatial resolution of 2mm\n",
"(i.e., each brain \"voxel\", or 3-dimensional pixel is 2 mm along each dimension),\n",
"then the resulting dataset might have approximately 20 x 1800 x 100 x 100 x 100\n",
"= 36 billion observations. That's a lot of data! Moreover, each subject's data\n",
"has a clear 4 dimensional structure — the 3 spatial dimensions, plus time. If we\n",
"wanted to, we could also potentially represent subjects as the 5th dimension,\n",
"though that involves some complications, since at least initially, different\n",
"subjects' brains won't be aligned with one another —- we'd need to spatially\n",
"register them for that (we'll learn more about that in {numref}`registration`).\n",
"\n",
"This means that a single subject's data would look a little bit like the\n",
"4-dimensional example that we looked at just a moment ago.\n",
"\n",
"We will typically want to access this data in pretty specific ways. That is,\n",
"rather than applying an operation to every single voxel in the brain, at every\n",
"single point in time, we usually want to pull out specific slices of the data,\n",
"and only apply an operation to those slices. Say for example we're interested in\n",
"a voxel in the amygdala. How would we access only that voxel, at every time\n",
"point?\n",
"\n",
"A very naive approach that we could implement in pure Python would be to store\n",
"all our data as a nested series of Python lists: each element in the first list\n",
"would be a list containing data for one time point; each element within the list\n",
"for each time point would itself be a list containing the 2d slices at each\n",
"x-coordinate; and so on. Then we could write a series of 4 nested for-loops to\n",
"sequentially access every data point (or voxel) in our array. For each voxel we\n",
"inspect, we could then determine whether the voxel is one we want to work with,\n",
"and if so, apply some operation to it.\n",
"\n",
"Basically, we'd have something like this (note that this is just pseudocode, not\n",
"valid Python code, and you wouldn't be able to execute this snippet!):\n",
"\n",
"```python\n",
"for t in time:\n",
" for x in t:\n",
" for y in x:\n",
" for z in y:\n",
" if z is in amygdala:\n",
" apply_my_function(z)\n",
"```\n",
"\n",
"Something like the above would probably work just fine. But it would be very\n",
"inefficient, and might take a while to execute. Even if we only want to access a\n",
"single amygdala voxel at each timepoint (i.e., 1800 data points in total, if we\n",
"stick with our resting state example from above), our naive looping approach\n",
"still requires us to inspect every single one of the 1.8 billion voxels in our\n",
"dataset —- even if all we're doing for the vast majority of voxels is deciding\n",
"that we don't actually need to use that voxel in our analysis! This is, to put\n",
"it mildly, A Bad Thing.\n",
"\n",
"Computer scientists have developed a notation for more precisely describing how\n",
"much of a bad thing it is: Big O notation, which gives us a formal way of\n",
"representing the time complexity of an algorithm. In Big O notation, the code\n",
"above runs in $O(n^4)$ time. Meaning, for any given dimension size $n$ (we\n",
"assume for simplicity that $n$ is the same for all dimensions), we need to\n",
"perform on the order of $n^4$ computations. To put this in perspective, just\n",
"going from $n = 2$ to $n = 4$ means that we go from $2^4 = 16$ to $4^4 = 256$\n",
"operations. You can see how this might start to become problematic as $n$ gets\n",
"large.\n",
"\n",
"The lesson we should take away from this is that writing naive Python loops just\n",
"to access the values at a specific index in our dataset is a bad idea.\n",
"\n",
"### Towards a specialized array structure\n",
"\n",
"What about list indexing? Fortunately, since we've already covered list indexing\n",
"in Python, we have a better way to access the amygdala data we want. We know\n",
"that the outer dimension in our nested list represents time, and we also know\n",
"that the spatial index of the amygdala voxel we want will be constant across\n",
"time points. This means that instead of writing a 4-level nested for-loop, we\n",
"can do the job in just a single loop. If we know that the index of the amygdala\n",
"voxel we want in our 3-d image is [20, 18, 32], then the following code suffices\n",
"(again note that this code can't be executed):\n",
"\n",
"```python\n",
"\n",
"amygdala_values = []\n",
"\n",
"for t in dataset:\n",
" amygdala_values.append(t[20][18][32])\n",
"```\n",
"\n",
"That's a big improvement, in terms of time complexity. We've gone from an\n",
"$O(n^4)$ to an $O(n)$ algorithm. The time our new algorithm takes to run is now\n",
"linear in the number of timepoints in our dataset, rather than growing as a\n",
"polynomial function of it.\n",
"\n",
"Unfortunately, our new approach still isn't great. One problem is that\n",
"performance still won't be so good, because Python isn't naturally optimized for\n",
"the kind of operation we just saw -— indexing repeatedly into nested lists. For\n",
"one thing, Python is a dynamically typed language, so there will be some\n",
"overhead involved in accessing each individual element (because the interpreter\n",
"has to examine the object and determine its type), and this can quickly add up\n",
"if we have millions or billions of elements.\n",
"\n",
"The other big limitation is that, even if lists themselves were super fast to\n",
"work with, the Python standard library provides limited functionality for doing\n",
"numerical analysis on them. There are a lot of basic numerical operations we\n",
"might want to apply to our data that would be really annoying to write using\n",
"only the built-in functionality of Python. For example, Python lists don't\n",
"natively support matrix operations, so if we wanted to multiply one matrix\n",
"(represented as a list of lists) by another, we'd probably need to write out a\n",
"series of summations and multiplications ourselves inside for-loops. And that\n",
"would be completely impractical given how basic and common an operation matrix\n",
"multiplication is (not to important for our point now, but if you're wondering\n",
"what a matrix multiplication is, you can see a formal definition in\n",
"{numref}`matmul`).\n",
"\n",
"The bottom line is that, as soon as we start working with large, or highly\n",
"structured datasets —- as is true of most datasets in neuroimaging —- we're\n",
"going to have to look beyond what Python provides us in its standard library.\n",
"Lists just aren't going to cut it; we need some other kind of data structure —-\n",
"one that's optimized for numerical analysis on large, N-dimensional datasets.\n",
"This is where numpy comes into play.\n",
"\n",
"## Introducing `Numpy`\n",
"\n",
"Numpy is the backbone of the numerical and scientific computing stack in Python;\n",
"many of the libraries we'll cover in this book (pandas, scikit-learn, etc.)\n",
"depend on it internally. Numpy provides many data structures optimized for\n",
"efficient representation and manipulation of different kinds of high-dimensional\n",
"data, as well as an enormous range of numerical tools that help us work with\n",
"those structures.\n",
"\n",
"### Importing numpy\n",
"\n",
"Recall that the default Python namespace contains only a small number of\n",
"built-in functions. To use any other functionality, we need to explicitly import\n",
"it into our namespace. Let's do that for Numpy. The library can be imported\n",
"using the `import` key word. We use `import numpy as np` to create a short name\n",
"for the library that we can refer to. Everything that is in `numpy` will now be\n",
"accessible through the short name `np`. Importing the library specifically as\n",
"`np` is not a requirement, but it is a very strongly-held convention and when\n",
"you read code that others have written, you will find that this the form that is\n",
"very often used."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "b14e7b9f",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"id": "cbe98f52",
"metadata": {},
"source": [
"### Data is represented in arrays\n",
"\n",
"The core data structure in Numpy is the n-dimensional array (or `ndarray`). As\n",
"the name suggests, an `ndarray` is an array with an arbitrary number of\n",
"dimensions. Unlike Python lists, Numpy arrays are homogeneously typed -—\n",
"meaning, every element in the array has to have the same data type. You can have\n",
"an array of floats, or an array of integers, but you can't have an array that\n",
"mixes floats and integers. This is exactly what we described above: a\n",
"homogeneous dense table that holds items in it.\n",
"\n",
"### Creating `ndarray`s\n",
"\n",
"Like any other Python object, we need to initialize an `ndarray` before we can\n",
"do anything with it. Numpy provides us with several ways to create new arrays.\n",
"Let's explore a couple.\n",
"\n",
"#### Initializing an array from an existing list\n",
"\n",
"Let's start by constructing an array from existing data. Assume we have some\n",
"values already stored in a native Python iterable object (typically, a list;\n",
"\"iterable\" means that it is an object that we can iterate over with a for loop,\n",
"see {numref}`forloop`). If we want to convert that object to an ndarray so that\n",
"we can perform more efficient numerical operations on it, we can pass the object\n",
"directly to the `np.narray()` function. In a relatively simple case, the input\n",
"to the `np.array` function is a sequence of items, a list that holds some\n",
"integers."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "2926a442",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1 1 2 3 5]\n"
]
}
],
"source": [
"list_of_numbers = [1, 1, 2, 3, 5]\n",
"arr1 = np.array(list_of_numbers)\n",
"print(arr1)"
]
},
{
"cell_type": "markdown",
"id": "a2afd6eb",
"metadata": {},
"source": [
"The result still looks a lot like a list of numbers, but let's look at some of\n",
"the attributes that this variable has that are not available in lists."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "6af23e3a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The shape of the array is (5,)\n",
"The dtype of the array is int64\n",
"The size of each item is 8\n",
"The strides of the array is (8,)\n"
]
}
],
"source": [
"print(\"The shape of the array is\", arr1.shape)\n",
"print(\"The dtype of the array is\", arr1.dtype)\n",
"print(\"The size of each item is\", arr1.itemsize)\n",
"print(\"The strides of the array is\", arr1.strides)"
]
},
{
"cell_type": "markdown",
"id": "09665abb",
"metadata": {},
"source": [
"The shape of an array is its dimensionality. In this case, the array is\n",
"one-dimensional (like the single channel of EEG data in the example above). All\n",
"of the items in this array are integers, so we can use an integer representation\n",
"as the dtype of the array. This means that each element of the array takes up 8\n",
"bytes of memory. Using 64-bit integers, we can represent the numbers from\n",
"-9223372036854775808 to 9223372036854775807. And each one would take up 8 bytes\n",
"of memory in our computer (which are 64 bits). This also explains the strides.\n",
"Strides are the number of bytes in memory that we would have to skip to get to\n",
"the next value in the array. Because the array is packed densely in a contiguous\n",
"segment of memory, to get from the value `1` to the value `2`, for example, we\n",
"have to skip 8 bytes.\n",
"\n",
"![](./figures/numpy_arrays1.png)"
]
},
{
"cell_type": "markdown",
"id": "f664eb4d",
"metadata": {},
"source": [
"Remember that the data we would like to represent in arrays can have multiple\n",
"dimensions. Numpy is pretty good about figuring out what kind of array we want\n",
"based on the input data; for example, if we pass in a list of lists, where each\n",
"of the inner lists has the same number of elements, Numpy will infer that we\n",
"want to construct a 2-dimensional array"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c99d7f4a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 1 1 2 3 5]\n",
" [ 8 13 21 34 55]]\n",
"The shape of this array is (2, 5)\n",
"The dtype of this array is int64\n",
"The size of each item is 8\n",
"The strides of this array are (40, 8)\n"
]
}
],
"source": [
"list_of_lists = [[1, 1, 2, 3, 5], [8, 13, 21, 34, 55]]\n",
"arr2 = np.array(list_of_lists)\n",
"\n",
"print(arr2)\n",
"print(\"The shape of this array is\", arr2.shape)\n",
"print(\"The dtype of this array is\", arr2.dtype)\n",
"print(\"The size of each item is\", arr2.itemsize)\n",
"print(\"The strides of this array are\", arr2.strides)"
]
},
{
"cell_type": "markdown",
"id": "e9c2441f",
"metadata": {},
"source": [
"Now, the shape of the array is 2 items: the first is the number of lists in the\n",
"array (2) and the second is the number of items in each list (5). You can also\n",
"think of that as rows (2 rows) and columns (5 columns), because each item in the\n",
"first row of the array has an item matching it in terms of its column in the\n",
"second row. For example, the number `2` in the first row is equivalent to the\n",
"number `21` in the second row in that they are both the third item in their\n",
"particular rows. Or, in other words, they share the third column in the array.\n",
"The items are still each 8 bytes in size, but now there are two items in the\n",
"strides: `(40, 8)`. The first item tells us how many bytes we have to move to\n",
"get from one row to the next row (for example, from the `2` in the first row to\n",
"the `21` in the second row). That's because the part of the computer's memory in\n",
"which the entire array is stored is contiguous and linear. Here, we've drawn\n",
"this contiguous bit of memory, such that the part that represents the first row\n",
"of the array is green and the part that represents the second row of the array\n",
"is blue.\n",
"\n",
"![](./figures/numpy_arrays2.png)"
]
},
{
"cell_type": "markdown",
"id": "e5ca8d63",
"metadata": {},
"source": [
"The first time you see a printed representation of a Numpy array, as in the\n",
"above output, it might look a little confusing. But the visual representation\n",
"closely resembles what we'd see if we were to print the nested list in Python.\n",
"In the 2-dimensional case above, the elements of the outer list are the values\n",
"along the first dimension, and each inner list gives the values along the second\n",
"dimension. This ends up giving us a nice tabular representation, where, at least\n",
"for small arrays, we can just read off the values from the output. For example,\n",
"the first row of the array contains the values [1, 1, 2, 3, 5]; the second\n",
"column of the array contains the values [1, 13]; and so on.\n",
"\n",
"(numpy_ex1)=\n",
"#### Exercise\n",
"\n",
"Extend the above principle into 3 dimensions: create a list of list of lists,\n",
"initialize an array from the list, and print into the screen. Make sure you\n",
"understand how the values displayed visually map onto the dimensions of the\n",
"array.\n",
"\n",
"#### Initializing an empty array\n",
"\n",
"```{eval-rst}\n",
".. index::\n",
" single: zeros\n",
"```\n",
"\n",
"Alternatively, we can construct a new array from scratch and fill it with some\n",
"predetermined value -- most commonly zero. We can do this easily in Numpy using\n",
"the conveniently - named `np.zeros` function. In many practical applications, we\n",
"can think of this as an \"empty\" array (though technically we could also create a\n",
"truly empty array that has no assigned values using the `np.empty` function).\n",
"\n",
"The `zeros` function takes a mandatory shape tuple as its first argument; this\n",
"specifies the dimensions of the desired array."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "e07a132b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The shape of this array is (5, 10)\n",
"The dtype of this array is float64\n",
"The size of each item is 8\n",
"The strides of this array are (80, 8)\n",
"The shape of this array is (2, 4, 8)\n",
"The dtype of this array is float64\n",
"The size of each item is 8\n",
"The strides of this array are (256, 64, 8)\n"
]
}
],
"source": [
"arr_2d = np.zeros((5, 10))\n",
"print(\"The shape of this array is\", arr_2d.shape)\n",
"print(\"The dtype of this array is\", arr_2d.dtype)\n",
"print(\"The size of each item is\", arr_2d.itemsize)\n",
"print(\"The strides of this array are\", arr_2d.strides)\n",
"\n",
"arr_3d = np.zeros((2, 4, 8))\n",
"print(\"The shape of this array is\", arr_3d.shape)\n",
"print(\"The dtype of this array is\", arr_3d.dtype)\n",
"print(\"The size of each item is\", arr_3d.itemsize)\n",
"print(\"The strides of this array are\", arr_3d.strides)"
]
},
{
"cell_type": "markdown",
"id": "a03ce371",
"metadata": {},
"source": [
"(numpy_ex2)=\n",
"#### Exercise\n",
"\n",
"Can you explain the outputs of the previous code cell? What does the first item in the `strides` represent in each case?\n",
"\n",
"### Neuroscience data is stored in arrays\n",
"\n",
"```{eval-rst}\n",
".. index::\n",
" single: Strides\n",
"```\n",
"\n",
"As we described above, many different kinds of neuroscience data are naturally\n",
"stored in arrays. To demonstrate this, let's look at some neuroscience data that\n",
"we've stored in an array. We've stored some data from an fMRI experiment in a\n",
"Numpy array, you can load it using our the `load_data` utility function from the\n",
"`nsdlib` software library that we have implemented:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "1b58f84b",
"metadata": {},
"outputs": [],
"source": [
"from ndslib import load_data"
]
},
{
"cell_type": "markdown",
"id": "650cfd71",
"metadata": {},
"source": [
"This `load_data` function takes as arguments the name of the dataset that we\n",
"would like to download and alternatively also the name of a file that we would\n",
"like to save the data into (for reuse). The `npy` file format is a format that\n",
"was developed specifically to store Numpy arrays and that's all that it knows\n",
"how to store. When you load one of these files, you can assign into a variable\n",
"the contents of the array that was stored in the file."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "6e6bcf7a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The shape of the data is (64, 64, 25, 180)\n",
"The dtype of the data is float64\n",
"The itemsize of the data is 8\n",
"The strides of the data is (8, 512, 32768, 819200)\n"
]
}
],
"source": [
"bold = load_data(\"bold_numpy\", fname=\"bold.npy\")\n",
"\n",
"print(\"The shape of the data is\", bold.shape)\n",
"print(\"The dtype of the data is\", bold.dtype)\n",
"print(\"The itemsize of the data is\", bold.itemsize)\n",
"print(\"The strides of the data is\", bold.strides)"
]
},
{
"cell_type": "markdown",
"id": "34a8537f",
"metadata": {},
"source": [
"This data is a little bit more complicated than the arrays you have seen so far.\n",
"It has 4 dimensions: the first two dimensions are the inplane dimensions of each\n",
"slice: 64 by 64 voxels. The next dimension is the slice dimension: 25 slices.\n",
"Finally, the last dimension is the time-points that were measured in the\n",
"experiment. The dtype of this array is `float64`. These floating point numbers\n",
"are what computers use to represent real numbers. Each one of these items also\n",
"uses up 64 bits, or 8 bytes, so that's the item size here as well. Just like the\n",
"shape, there are 4 strides: the first item in the strides tells how far we need\n",
"to go to get from one voxel to another voxel in the same row, in the same slice,\n",
"in the same time-point. The second stride (512) tells us how far we would have\n",
"to go to find the same voxel in the same slice and time-point, but one row over.\n",
"That's 64 times 8. The next one is about finding the same voxel in a neighboring\n",
"slice. Finally the large one at the end is about finding the same slice in the\n",
"next time step. That large number is equal to 64 by 64 by 25 (the number of\n",
"voxels in each time point) by 8 (the number of bytes per element in the array)."
]
},
{
"cell_type": "markdown",
"id": "485d13d2",
"metadata": {},
"source": [
"### Indexing is used to access elements of an array\n",
"\n",
"```{eval-rst}\n",
".. index::\n",
" single: Indexing (arrays)\n",
"```\n",
"\n",
"We've seen how we can create arrays and describe them; now let's talk about how\n",
"we can get data in and out of arrays. We already know how to index Python lists,\n",
"and array indexing in Python will look quite similar. But Numpy indexing adds\n",
"considerably more flexibility and power, and developing array indexing facility\n",
"is a critical step towards acquiring general proficiency with the package. Let's\n",
"start simple: accessing individual items in a one-dimensional array looks a lot\n",
"like indexing into a list (see {numref}`listindexing`). For example, the first\n",
"item in the one-dimensional array we saw earlier, you would do the following:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "677dd999",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr1[0]"
]
},
{
"cell_type": "markdown",
"id": "e1b42276",
"metadata": {},
"source": [
"or if you want the second element, you would do something like this (don't forget that Python uses 0-based indexing!):"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "18d11802",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr1[1]"
]
},
{
"cell_type": "markdown",
"id": "aa682cfb",
"metadata": {},
"source": [
"Numpy arrays also support two other important syntactic conventions that are\n",
"also found in Python lists: indexing from the end of the array, and \"slicing\"\n",
"the array to extract multiple contiguous values.\n",
"\n",
"To index from the end, we use the minus sign (-):"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "2cec3381",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr1[-2]"
]
},
{
"cell_type": "markdown",
"id": "a0c0793a",
"metadata": {},
"source": [
"This gives us the 2nd-from-last value in the array. To slice an array, we use\n",
"the colon (:) operator, passing in the positions we want to start and end at:\n",
"\n",
"```{eval-rst}\n",
".. index::\n",
" single: Slicing (arrays)\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "7b52f62f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([2, 3, 5])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr1[2:5]"
]
},
{
"cell_type": "markdown",
"id": "b5a9cc0b",
"metadata": {},
"source": [
"As in lists, the start index is inclusive and the end index is exclusive (i.e.,\n",
"in the above example, the resulting array includes the value at position 2, but\n",
"excludes the one at position 6). We can also omit the start or end indexes, in\n",
"which case, Numpy will return all positions up to, or starting from, the provided\n",
"index. For example, to get the first 4 elements:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "e828f630",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 1, 2, 3])"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr1[:4]"
]
},
{
"cell_type": "markdown",
"id": "a464a05e",
"metadata": {},
"source": [
"### Indexing in multiple dimensions\n",
"\n",
"Once we start working with arrays with more than one dimension, indexing gets a\n",
"little more complicated. Both because the syntax is a little different, and\n",
"because there's a lot more we can do with multi-dimensional arrays.\n",
"\n",
"For example, if you try one of these operations with a two-dimensional array,\n",
"something funny will happen:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "0245f797",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 1, 2, 3, 5])"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr2[0]"
]
},
{
"cell_type": "markdown",
"id": "951e79db",
"metadata": {},
"source": [
"It looks like the first item of a two-dimensional array is the entire first row\n",
"of that array. To access the first element in the first row, we need to use two\n",
"indices, separated by a comma:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "daae0019",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr2[0, 0]"
]
},
{
"cell_type": "markdown",
"id": "965638ed",
"metadata": {},
"source": [
"One way to think about this is that the first index points to the row and the\n",
"second index points to the column, so to get the very first item in the array,\n",
"we need to explicitly ask for the first column in the first row. Consider what\n",
"happens when you ask for third item in the second row:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "ca81f289",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"21"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arr2[1, 2]"
]
},
{
"cell_type": "markdown",
"id": "5e1b65f7",
"metadata": {},
"source": [
"Under the hood, Python asks how many bytes into the array it would have to go to\n",
"get that item. It looks at the strides and computes: 1 times the first stride =\n",
"40 plus 2 times the second stride = 40 + 16 = 56. Indeed, the third element of\n",
"the second row of the array is 7 times 8 bytes into the array:\n",
"\n",
"![](./figures/numpy_arrays3.png)"
]
},
{
"cell_type": "markdown",
"id": "da36834d",
"metadata": {},
"source": [
"How would you access data in the BOLD array? As an example, let's access the\n",
"voxel in the very first time-point, in the very first slice, and in this slice\n",
"in the first row and column:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "7157ab4d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[0, 0, 0, 0]"
]
},
{
"cell_type": "markdown",
"id": "0061817e",
"metadata": {},
"source": [
"This element of the array is the very corner of the volume, so it is way outside\n",
"of the part of the data that contains the brain. That is why it contains a value\n",
"of 0. Instead, we can look at the voxel at the very center of the center slice:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "8beea2bb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1080.0"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[32, 32, 12, 0]"
]
},
{
"cell_type": "markdown",
"id": "e40cbcd6",
"metadata": {},
"source": [
"This voxel has a value that corresponds to the intensity of the MR image in this\n",
"voxel at the first time point. If we drop the last index, that would be like\n",
"asking for the data in all of the time-points in that particular voxel:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "d3de489c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1080., 1052., 1056., 1087., 1146., 1147., 1105., 1064., 1128.,\n",
" 1089., 1095., 1049., 1109., 1051., 1074., 1073., 1112., 1086.,\n",
" 1090., 1062., 1086., 1023., 1047., 1139., 1065., 1117., 1070.,\n",
" 1070., 1089., 1074., 1051., 1034., 1096., 1060., 1096., 1076.,\n",
" 1032., 1067., 1030., 1072., 1056., 1069., 1061., 1054., 1072.,\n",
" 1072., 1035., 1018., 1116., 1056., 1051., 1084., 1075., 1080.,\n",
" 1036., 1022., 1076., 1060., 1031., 1079., 1048., 1002., 1055.,\n",
" 1027., 1014., 1006., 1072., 1003., 1026., 1039., 1096., 1078.,\n",
" 1025., 1029., 1009., 1065., 1023., 1098., 1045., 1094., 1016.,\n",
" 1015., 1027., 1020., 1030., 1049., 1034., 1053., 1018., 1038.,\n",
" 1072., 1020., 1007., 1037., 1082., 1050., 1011., 1027., 972.,\n",
" 992., 984., 1021., 1018., 1029., 1063., 1070., 1033., 1048.,\n",
" 1052., 1040., 1007., 996., 1026., 1011., 1023., 958., 982.,\n",
" 1008., 1059., 1006., 1012., 1042., 1000., 1066., 1026., 1032.,\n",
" 1009., 1022., 989., 1020., 1078., 1034., 1036., 976., 1026.,\n",
" 960., 1021., 1009., 1049., 1029., 1065., 1002., 1007., 1031.,\n",
" 1015., 1044., 1012., 999., 1002., 979., 965., 1011., 1027.,\n",
" 1013., 1012., 1022., 981., 1030., 1061., 1056., 1014., 957.,\n",
" 968., 1015., 1100., 994., 1000., 973., 997., 980., 1026.,\n",
" 1004., 1008., 988., 993., 1042., 1016., 1044., 1002., 1006.])"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[32, 32, 12]"
]
},
{
"cell_type": "markdown",
"id": "0fab8bcc",
"metadata": {},
"source": [
"This is now a time-series of data, corresponding to the intensity of the BOLD\n",
"response in the voxel at the very center of the central slice, across all time\n",
"points in the measurement."
]
},
{
"cell_type": "markdown",
"id": "b48c7b06",
"metadata": {},
"source": [
"### Slicing a dimension with \"`:`\"\n",
"\n",
"As you saw before, we can use slicing to choose a sub-array from within our\n",
"array. For example, you can slice out the time-series in the central voxel, as\n",
"we just did:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "104af3eb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1080., 1052., 1056., 1087., 1146., 1147., 1105., 1064., 1128.,\n",
" 1089., 1095., 1049., 1109., 1051., 1074., 1073., 1112., 1086.,\n",
" 1090., 1062., 1086., 1023., 1047., 1139., 1065., 1117., 1070.,\n",
" 1070., 1089., 1074., 1051., 1034., 1096., 1060., 1096., 1076.,\n",
" 1032., 1067., 1030., 1072., 1056., 1069., 1061., 1054., 1072.,\n",
" 1072., 1035., 1018., 1116., 1056., 1051., 1084., 1075., 1080.,\n",
" 1036., 1022., 1076., 1060., 1031., 1079., 1048., 1002., 1055.,\n",
" 1027., 1014., 1006., 1072., 1003., 1026., 1039., 1096., 1078.,\n",
" 1025., 1029., 1009., 1065., 1023., 1098., 1045., 1094., 1016.,\n",
" 1015., 1027., 1020., 1030., 1049., 1034., 1053., 1018., 1038.,\n",
" 1072., 1020., 1007., 1037., 1082., 1050., 1011., 1027., 972.,\n",
" 992., 984., 1021., 1018., 1029., 1063., 1070., 1033., 1048.,\n",
" 1052., 1040., 1007., 996., 1026., 1011., 1023., 958., 982.,\n",
" 1008., 1059., 1006., 1012., 1042., 1000., 1066., 1026., 1032.,\n",
" 1009., 1022., 989., 1020., 1078., 1034., 1036., 976., 1026.,\n",
" 960., 1021., 1009., 1049., 1029., 1065., 1002., 1007., 1031.,\n",
" 1015., 1044., 1012., 999., 1002., 979., 965., 1011., 1027.,\n",
" 1013., 1012., 1022., 981., 1030., 1061., 1056., 1014., 957.,\n",
" 968., 1015., 1100., 994., 1000., 973., 997., 980., 1026.,\n",
" 1004., 1008., 988., 993., 1042., 1016., 1044., 1002., 1006.])"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[32, 32, 12, :]"
]
},
{
"cell_type": "markdown",
"id": "6dd79bdd",
"metadata": {},
"source": [
"This is equivalent to dropping the last index. If we want to slice through all\n",
"of the slices, taking the values of the central voxel in each slice, in the very\n",
"first time-point, we would instead slice on the third (slice) dimension of the\n",
"array:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "41312427",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 31., 52., 117., 423., 77., 312., 454., 806., 576.,\n",
" 581., 403., 1031., 1080., 482., 439., 494., 546., 502.,\n",
" 453., 362., 385., 367., 344., 169., 25.])"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[32, 32, :, 0]"
]
},
{
"cell_type": "markdown",
"id": "c71e35a4",
"metadata": {},
"source": [
"These 25 values correspond to the values in the center of each slice in the\n",
"first time point. We can also slice on multiple dimensions of the array. For\n",
"example, if we would like to get the values of all the voxels in the central\n",
"slice in the very first time-point, we can do:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "3f1e2aaf",
"metadata": {},
"outputs": [],
"source": [
"central_slice_t0 = bold[:, :, 12, 0]"
]
},
{
"cell_type": "markdown",
"id": "dbc055eb",
"metadata": {},
"source": [
"This is a two-dimensional array, with the dimensions of the slice:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "ed22e379",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(64, 64)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"central_slice_t0.shape"
]
},
{
"cell_type": "markdown",
"id": "8a2627cf",
"metadata": {},
"source": [
"Which we can verify by printing out the values inside the array:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "e2c949e9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0. 0. 0. ... 0. 0. 0.]\n",
" [13. 4. 3. ... 16. 3. 6.]\n",
" [ 6. 13. 6. ... 13. 26. 7.]\n",
" ...\n",
" [13. 11. 17. ... 16. 6. 5.]\n",
" [10. 3. 7. ... 8. 6. 8.]\n",
" [ 7. 7. 7. ... 20. 16. 4.]]\n"
]
}
],
"source": [
"print(central_slice_t0)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "57d1d419",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(64, 64, 25, 180)"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold.shape"
]
},
{
"cell_type": "markdown",
"id": "b434ad08",
"metadata": {},
"source": [
"(boolean_indexing_numpy)=\n",
"### Indexing with conditionals\n",
"\n",
"One more way of indexing is to to use logical operations, or conditionals. This\n",
"allows us to choose values from an array, only if they fulfill certain\n",
"conditions. For example, if we'd like to exclude measurements with a value of 0\n",
"or less, we might perform a logical operation like the following:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "d53f0105",
"metadata": {},
"outputs": [],
"source": [
"larger_than_0 = bold>0"
]
},
{
"cell_type": "markdown",
"id": "bb6bda9a",
"metadata": {},
"source": [
"This creates a new array that has the same shape as the original array, but\n",
"contains Boolean (`True`/`False`) values. Wherever we had a number larger than 0\n",
"in the original array, this new array will have the value `True` and wherever\n",
"the numbers in the original array were 0 or smaller, this new array will have\n",
"the value `False`. This is already helpful in some regards, but the next thing\n",
"that we can do with this is that we can use this array to index into the\n",
"original array. What this does is that it pulls out of the original array only\n",
"the values for which `larger_than_0` is `True` and reorganizes these values as a\n",
"long one-dimensional array:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "669252e8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 6., 11., 5., ..., 12., 4., 14.])"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[larger_than_0]"
]
},
{
"cell_type": "markdown",
"id": "6c854816",
"metadata": {},
"source": [
"This can be useful, for example, if you would like to calculate statistics of\n",
"the original array, but only above some threshold value. For example, if you\n",
"would like to calculate the average of the array (using the `mean` method of the\n",
"array) with and without values that are 0 or less:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "f5e8f1c7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"145.360578125\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"148.077239538181\n"
]
}
],
"source": [
"print(bold.mean())\n",
"print(bold[larger_than_0].mean())"
]
},
{
"cell_type": "markdown",
"id": "85b3582f",
"metadata": {},
"source": [
"(numpy_ex3)=\n",
"#### Exercise\n",
"\n",
"Arrays have a variety of methods attached to them to do calculations. Use the\n",
"array method `std` to calculate the standard deviation of the `bold` array, and\n",
"the standard deviation of just values of the array that are between 1 and 100.\n",
"\n",
"\n",
"### Arithmetic with arrays\n",
"\n",
"In addition to storing data from a measurement, and providing access to the this\n",
"data through indexing, arrays also support a variety of mathematical operations\n",
"with arrays. For example, you can perform arithmetic operations between an array\n",
"and a number. Consider for example what happens when we add a single number to\n",
"an entire array:"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "a338863f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([2080., 2052., 2056., 2087., 2146., 2147., 2105., 2064., 2128.,\n",
" 2089., 2095., 2049., 2109., 2051., 2074., 2073., 2112., 2086.,\n",
" 2090., 2062., 2086., 2023., 2047., 2139., 2065., 2117., 2070.,\n",
" 2070., 2089., 2074., 2051., 2034., 2096., 2060., 2096., 2076.,\n",
" 2032., 2067., 2030., 2072., 2056., 2069., 2061., 2054., 2072.,\n",
" 2072., 2035., 2018., 2116., 2056., 2051., 2084., 2075., 2080.,\n",
" 2036., 2022., 2076., 2060., 2031., 2079., 2048., 2002., 2055.,\n",
" 2027., 2014., 2006., 2072., 2003., 2026., 2039., 2096., 2078.,\n",
" 2025., 2029., 2009., 2065., 2023., 2098., 2045., 2094., 2016.,\n",
" 2015., 2027., 2020., 2030., 2049., 2034., 2053., 2018., 2038.,\n",
" 2072., 2020., 2007., 2037., 2082., 2050., 2011., 2027., 1972.,\n",
" 1992., 1984., 2021., 2018., 2029., 2063., 2070., 2033., 2048.,\n",
" 2052., 2040., 2007., 1996., 2026., 2011., 2023., 1958., 1982.,\n",
" 2008., 2059., 2006., 2012., 2042., 2000., 2066., 2026., 2032.,\n",
" 2009., 2022., 1989., 2020., 2078., 2034., 2036., 1976., 2026.,\n",
" 1960., 2021., 2009., 2049., 2029., 2065., 2002., 2007., 2031.,\n",
" 2015., 2044., 2012., 1999., 2002., 1979., 1965., 2011., 2027.,\n",
" 2013., 2012., 2022., 1981., 2030., 2061., 2056., 2014., 1957.,\n",
" 1968., 2015., 2100., 1994., 2000., 1973., 1997., 1980., 2026.,\n",
" 2004., 2008., 1988., 1993., 2042., 2016., 2044., 2002., 2006.])"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[32, 32, 12] + 1000"
]
},
{
"cell_type": "markdown",
"id": "fc888150",
"metadata": {},
"source": [
"The result of this operation is that the number is added to every item in the\n",
"array. Similarly:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "e7728594",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([540. , 526. , 528. , 543.5, 573. , 573.5, 552.5, 532. , 564. ,\n",
" 544.5, 547.5, 524.5, 554.5, 525.5, 537. , 536.5, 556. , 543. ,\n",
" 545. , 531. , 543. , 511.5, 523.5, 569.5, 532.5, 558.5, 535. ,\n",
" 535. , 544.5, 537. , 525.5, 517. , 548. , 530. , 548. , 538. ,\n",
" 516. , 533.5, 515. , 536. , 528. , 534.5, 530.5, 527. , 536. ,\n",
" 536. , 517.5, 509. , 558. , 528. , 525.5, 542. , 537.5, 540. ,\n",
" 518. , 511. , 538. , 530. , 515.5, 539.5, 524. , 501. , 527.5,\n",
" 513.5, 507. , 503. , 536. , 501.5, 513. , 519.5, 548. , 539. ,\n",
" 512.5, 514.5, 504.5, 532.5, 511.5, 549. , 522.5, 547. , 508. ,\n",
" 507.5, 513.5, 510. , 515. , 524.5, 517. , 526.5, 509. , 519. ,\n",
" 536. , 510. , 503.5, 518.5, 541. , 525. , 505.5, 513.5, 486. ,\n",
" 496. , 492. , 510.5, 509. , 514.5, 531.5, 535. , 516.5, 524. ,\n",
" 526. , 520. , 503.5, 498. , 513. , 505.5, 511.5, 479. , 491. ,\n",
" 504. , 529.5, 503. , 506. , 521. , 500. , 533. , 513. , 516. ,\n",
" 504.5, 511. , 494.5, 510. , 539. , 517. , 518. , 488. , 513. ,\n",
" 480. , 510.5, 504.5, 524.5, 514.5, 532.5, 501. , 503.5, 515.5,\n",
" 507.5, 522. , 506. , 499.5, 501. , 489.5, 482.5, 505.5, 513.5,\n",
" 506.5, 506. , 511. , 490.5, 515. , 530.5, 528. , 507. , 478.5,\n",
" 484. , 507.5, 550. , 497. , 500. , 486.5, 498.5, 490. , 513. ,\n",
" 502. , 504. , 494. , 496.5, 521. , 508. , 522. , 501. , 503. ])"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[32, 32, 12] / 2"
]
},
{
"cell_type": "markdown",
"id": "6ee1af9f",
"metadata": {},
"source": [
"In this case, the number is used to separately divide each one of the elements\n",
"of the array.\n",
"\n",
"We can also do arithmetic between arrays. For example, we can add the values in\n",
"this voxel to the items in the central voxel in the neighboring slice:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "0ae44ec5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([2111., 2079., 2026., 2048., 2009., 1999., 1991., 1926., 2019.,\n",
" 1940., 1964., 1969., 2015., 1979., 1942., 1976., 1993., 2016.,\n",
" 1999., 1996., 2054., 1962., 1973., 2017., 1940., 2041., 1966.,\n",
" 1960., 1957., 2000., 1966., 1961., 1974., 1965., 1953., 2033.,\n",
" 1973., 2019., 1970., 1946., 1976., 1994., 1974., 1992., 1964.,\n",
" 2017., 1961., 1946., 2015., 1959., 1959., 1999., 1945., 2031.,\n",
" 1943., 1957., 2024., 1954., 1932., 2002., 1933., 1935., 2007.,\n",
" 2012., 1917., 1931., 2032., 1953., 1976., 1951., 1983., 1966.,\n",
" 1981., 1940., 1944., 1978., 1935., 2001., 1963., 1974., 1965.,\n",
" 1939., 1997., 1930., 1959., 1952., 1950., 1981., 1982., 1938.,\n",
" 1967., 1962., 1928., 1906., 1963., 1941., 1927., 1959., 1998.,\n",
" 1971., 1891., 1975., 1938., 1985., 1939., 1940., 1914., 1939.,\n",
" 1951., 1927., 1909., 1928., 1980., 1893., 1934., 1894., 1953.,\n",
" 1901., 1989., 1881., 1928., 1959., 1914., 1962., 1908., 1948.,\n",
" 1929., 1941., 1917., 1911., 1956., 1943., 1958., 1899., 1961.,\n",
" 1917., 1933., 1922., 1920., 1911., 1944., 1908., 1916., 1969.,\n",
" 1913., 1922., 1890., 1937., 1925., 1913., 1876., 1870., 1892.,\n",
" 1957., 1932., 1930., 1906., 1882., 1913., 1946., 1939., 1947.,\n",
" 1822., 1861., 1953., 1909., 1911., 1898., 1969., 1832., 1945.,\n",
" 1878., 1952., 1885., 1918., 1931., 1892., 1913., 1950., 1911.])"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bold[32, 32, 12] + bold[32, 32, 11]"
]
},
{
"cell_type": "markdown",
"id": "9673b91e",
"metadata": {},
"source": [
"```{eval-rst}\n",
".. index::\n",
" single: Array-scalar arithmetic\n",
"```\n",
"\n",
"\n",
"```{eval-rst}\n",
".. index::\n",
" single: Array-array arithmetic\n",
"```\n",
"\n",
"\n",
"When we perform arithmetic operations between an array and a number (or\n",
"array-scalar operations), the number is separately added to each element. Here,\n",
"when we perform arithmetic operations between arrays (array-array operations),\n",
"the operation is instead performed element-by-element. In both cases the result\n",
"is an array that has the same number of items as the original array (or arrays)."
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "83de98da",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"27.384290220586703\n"
]
}
],
"source": [
"tsnr = bold.mean(-1) / bold.std(-1)\n",
"\n",
"print(tsnr[tsnr.shape[0]//2, tsnr.shape[1]//2, tsnr.shape[2]//2])"
]
},
{
"cell_type": "markdown",
"id": "019176d9",
"metadata": {},
"source": [
"(numpy_ex4)=\n",
"#### Exercises\n",
"\n",
"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]\n",
"\n",
"2. Creating an array of 1's: create an array of shape (3, 5) containing only the number 1\n",
"\n",
"3. Find all the items in an array containing numbers that have values that are larger than 0 and smaller or equal to 100\n",
"\n",
"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).\n",
"\n",
"5. Generate an array with 100 randomly distributed numbers between 0 and 10 (hint: look at the `np.random` sub-module).\n",
"\n",
"\n",
"### The Scipy library\n",
"\n",
"The Numpy array serves as the basis for numerical computing in Python. But it\n",
"does only so much. Scientific computing involves a variety of different\n",
"algorithms that do more interesting things with numbers. For example, we'd like\n",
"to use fancy image processing and signal processing algorithms with our data. Or\n",
"rely on various facts about statistical distributions. Over the last 100 years,\n",
"scientists and engineers have developed many different algorithms for scientific\n",
"computing. In the Python scientific computing ecosystem, many of the fundamental\n",
"algorithms are collected together in one library, called ``Scipy``. It provides\n",
"everything but the kitchen sink, in terms of basic algorithms for science.\n",
"Because it's such a large collection of things, it is composed of several\n",
"different sub-libraries that perform various scientific computing tasks, such as\n",
"statistics, image processing, clustering, linear algebra, optimization, etc. We\n",
"will revisit some of these in later chapters. For now, as an example, we will\n",
"focus on the signal processing sub-library ``scipy.signal``. We import the\n",
"sub-library and name it `sps`:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "1d5a1a54",
"metadata": {},
"outputs": [],
"source": [
"import scipy.signal as sps"
]
},
{
"cell_type": "markdown",
"id": "97d6ad64",
"metadata": {},
"source": [
"Remember that the BOLD dataset that we were looking at before had 180\n",
"time-points. That means that the shape of this data was `(64, 64, 25, 180)`. One\n",
"basic signal processing operation is to efficiently and accurately resample a\n",
"time-series into a different sampling rate. For example, let's imagine that we'd\n",
"like to resample the data to a slower sampling rate. The original data was\n",
"collected with a sampling interval or repetition time (TR) of 2.0 seconds. If we\n",
"would like to resample this to half the sampling rate (TR=4.0 seconds), we would\n",
"halve the number of samples on the last dimension using the `sps.resample`\n",
"function:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "950f1ffa",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(64, 64, 25, 90)\n"
]
}
],
"source": [
"resampled_bold = sps.resample(bold, 90, axis=3)\n",
"print(resampled_bold.shape)"
]
},
{
"cell_type": "markdown",
"id": "2843d570",
"metadata": {},
"source": [
"If, instead, we'd like to resample this to twice the sampling rate (TR=1.0), we\n",
"would double the number of samples."
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "b913911d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(64, 64, 25, 360)\n"
]
}
],
"source": [
"resampled_bold = sps.resample(bold, 360, axis=3)\n",
"print(resampled_bold.shape)"
]
},
{
"cell_type": "markdown",
"id": "2794c954",
"metadata": {},
"source": [
"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.\n",
"\n",
"(numpy_ex5)=\n",
"#### Exercise\n",
"\n",
"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.\n",
"\n",
"### Additional resources\n",
"\n",
"For a description of the scientific Python ecosystem, see Fernando Perez, Brian\n",
"Granger and John Hunter's paper from 2011 {cite}`Perez2010-nn`. This paper\n",
"describes in some detail how Python has grown from a general-purpose programming\n",
"language into a flexible and powerful tool that covers many of the uses that\n",
"scientists have, through a distributed, eco-system of development. More recent\n",
"papers describe the tools at the core of the ecosystem: The Numpy library is\n",
"described in one recent paper {cite}`Harris2020-so` and the Scipy library is\n",
"described in another {cite}`Virtanen2020-jl`."
]
}
],
"metadata": {
"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"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
},
"source_map": [
15,
20,
52,
90,
232,
234,
263,
267,
272,
277,
293,
301,
310,
330,
366,
378,
398,
400,
409,
416,
434,
452,
454,
458,
460,
468,
470,
480,
482,
490,
492,
503,
505,
511,
513,
520,
522,
531,
537,
539,
545,
547,
553,
555,
561,
569,
571,
578,
580,
587,
589,
593,
595,
599,
603,
605,
615,
617,
629,
631,
638,
641,
659,
661,
666,
668,
676,
678,
698,
702,
736,
738,
750,
753,
758,
761
]
},
"nbformat": 4,
"nbformat_minor": 5
}