From 2c905e2cb694ffc30bf2410c2696acee81e1ebae Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Thu, 1 Feb 2018 14:40:32 +0000
Subject: [PATCH] I think numpy is finished
---
getting_started/04_numpy.ipynb | 559 ++++++++++++++++++++++++++++-----
getting_started/04_numpy.md | 26 +-
2 files changed, 497 insertions(+), 88 deletions(-)
diff --git a/getting_started/04_numpy.ipynb b/getting_started/04_numpy.ipynb
index 545ae75..e877863 100644
--- a/getting_started/04_numpy.ipynb
+++ b/getting_started/04_numpy.ipynb
@@ -19,28 +19,40 @@
"alternative to Matlab as a scientific computing platform.\n",
"\n",
"\n",
+ "The `fslpython` environment currently includes [Numpy\n",
+ "1.11.1](https://docs.scipy.org/doc/numpy-1.11.0/index.html), which is a little\n",
+ "out of date, but we will update it for the next release of FSL.\n",
+ "\n",
+ "\n",
"## Contents\n",
"\n",
"\n",
"* [The Python list versus the Numpy array](#the-python-list-versus-the-numpy-array)\n",
"* [Numpy basics](#numpy-basics)\n",
" * [Creating arrays](#creating-arrays)\n",
- " * [Operating on arrays](#operating-on-arrays)\n",
+ " * [Loading text files](#loading-text-files)\n",
" * [Array properties](#array-properties)\n",
" * [Descriptive statistics](#descriptive-statistics)\n",
" * [Reshaping and rearranging arrays](#reshaping-and-rearranging-arrays)\n",
- "* [Multi-variate operations](#multi-variate-operations)\n",
+ "* [Operating on arrays](#operating-on-arrays)\n",
+ " * [Scalar operations](#scalar-operations)\n",
+ " * [Multi-variate operations](#multi-variate-operations)\n",
" * [Matrix multplication](#matrix-multiplication)\n",
" * [Broadcasting](#broadcasting)\n",
+ " * [Linear algebra](#linear-algebra)\n",
"* [Array indexing](#array-indexing)\n",
" * [Indexing multi-dimensional arrays](#indexing-multi-dimensional-arrays)\n",
" * [Boolean indexing](#boolean-indexing)\n",
" * [Coordinate array indexing](#coordinate-array-indexing)\n",
- "* [Generating random numbers](#generating-random-numbers)\n",
+ "* [Exercises](#exercises)\n",
+ " * [Load an array from a file and do stuff with it](#load-an-array-from-a-file-and-do-stuff-with-it)\n",
+ " * [Concatenate affine transforms](#concatenate-affine-transforms)\n",
"\n",
- "* [Appendix: Importing Numpy](#appendix-importing-numpy)\n",
- "* [Appendix: Vectors in Numpy](#appendix-vectors-in-numpy)\n",
+ "* [Appendix A: Generating random numbers](#appendix-generating-random-numbers)\n",
+ "* [Appendix B: Importing Numpy](#appendix-importing-numpy)\n",
+ "* [Appendix C: Vectors in Numpy](#appendix-vectors-in-numpy)\n",
"\n",
+ "* [Useful references](#useful-references)\n",
"\n",
"\n",
"<a class=\"anchor\" id=\"the-python-list-versus-the-numpy-array\"></a>\n",
@@ -98,13 +110,12 @@
"computing!\n",
"\n",
"\n",
- "This is a major source of confusion for those poor souls who have spent their\n",
- "lives working in Matlab, but have finally seen the light and switched to\n",
- "Python. It is _crucial_ to be able to distinguish between a Python list and a\n",
- "Numpy array.\n",
+ "This is a major source of confusion for people who are learning Python, and\n",
+ "are trying to write efficient code. It is _crucial_ to be able to distinguish\n",
+ "between a Python list and a Numpy array.\n",
"\n",
"\n",
- "___Python list == Matlab cell array:___ A list in python is akin to a cell\n",
+ "___Python list == Matlab cell array:___ A list in Python is akin to a cell\n",
"array in Matlab - they can store anything, but are extremely inefficient, and\n",
"unwieldy when you have more than a couple of dimensions.\n",
"\n",
@@ -216,7 +227,7 @@
"print('np.ones gives us ones: ', np.ones(5))\n",
"print('np.arange gives us a range: ', np.arange(5))\n",
"print('np.linspace gives us N linearly spaced numbers:', np.linspace(0, 1, 5))\n",
- "print('np.random.random gives us random numbers: ', np.random.random(5))\n",
+ "print('np.random.random gives us random numbers [0-1]:', np.random.random(5))\n",
"print('np.random.randint gives us random integers: ', np.random.randint(1, 10, 5))\n",
"print('np.eye gives us an identity matrix:')\n",
"print(np.eye(4))\n",
@@ -228,7 +239,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "> There will be more on random numbers [below](#generating-random-numbers).\n",
+ "> See the [appendix](#appendix-generating-random-numbers) for more information\n",
+ "> on generating random numbers in Numpy.\n",
"\n",
"\n",
"The `zeros` and `ones` functions can also be used to generate N-dimensional\n",
@@ -255,12 +267,12 @@
"> second to columns - just like in Matlab.\n",
"\n",
"\n",
- "<a class=\"anchor\" id=\"operating-on-arrays\"></a>\n",
- "### Operating on arrays\n",
+ "<a class=\"anchor\" id=\"loading-text-files\"></a>\n",
+ "### Loading text files\n",
"\n",
"\n",
- "All of the mathematical operators you know and love can be applied to Numpy\n",
- "arrays:"
+ "The `numpy.loadtxt` function is capable of loading numerical data from\n",
+ "plain-text files. By default it expects space-separated data:"
]
},
{
@@ -269,15 +281,55 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (3, 3))\n",
- "print('a:')\n",
- "print(a)\n",
- "print('a + 2:')\n",
- "print( a + 2)\n",
- "print('a * 3:')\n",
- "print( a * 3)\n",
- "print('a % 2:')\n",
- "print( a % 2)"
+ "data = np.loadtxt('04_numpy/space_separated.txt')\n",
+ "print('data in 04_numpy/space_separated.txt:')\n",
+ "print(data)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "But you can also specify the delimiter to expect<sup>1</sup>:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "data = np.loadtxt('04_numpy/comma_separated.txt', delimiter=',')\n",
+ "print('data in 04_numpy/comma_separated.txt:')\n",
+ "print(data)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> <sup>1</sup> And many other things such as file headers, footers, comments,\n",
+ "> and newline characters - see the\n",
+ "> [docs](https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html)\n",
+ "> for more information.\n",
+ "\n",
+ "\n",
+ " Of course you can also save data out to a text file just as easily, with\n",
+ "[`numpy.savetxt`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "data = np.random.randint(1, 10, (10, 10))\n",
+ "np.savetxt('mydata.txt', data, delimiter=',', fmt='%i')\n",
+ "\n",
+ "with open('mydata.txt', 'rt') as f:\n",
+ " for line in f:\n",
+ " print(line.strip())"
]
},
{
@@ -323,8 +375,8 @@
"### Descriptive statistics\n",
"\n",
"\n",
- "Similarly, a Numpy array has a set of methods<sup>1</sup> which allow you to\n",
- "calculate basic descriptive statisics on an array:"
+ "Similarly, a Numpy array has a set of methods<sup>2</sup> which allow you to\n",
+ "calculate basic descriptive statistics on an array:"
]
},
{
@@ -350,11 +402,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "> <sup>1</sup> Python, being an object-oriented language, distinguishes\n",
- "> between _functions_ and _methods_. _Method_ is simply the term used to refer\n",
- "> to a function that is associated with a specific object. Similarly, the term\n",
- "> _attribute_ is used to refer to some piece of information that is attached\n",
- "> to an object, such as `z.shape`, or `z.dtype`.\n",
+ "> <sup>2</sup> Python, being an object-oriented language, distinguishes\n",
+ "> between _functions_ and _methods_. Hopefully we all know what a function is\n",
+ "> - a _method_ is simply the term used to refer to a function that is\n",
+ "> associated with a specific object. Similarly, the term _attribute_ is used\n",
+ "> to refer to some piece of information that is attached to an object, such as\n",
+ "> `z.shape`, or `z.dtype`.\n",
"\n",
"\n",
"<a class=\"anchor\" id=\"reshaping-and-rearranging-arrays\"></a>\n",
@@ -370,7 +423,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (4, 4))\n",
+ "a = np.arange(16).reshape(4, 4)\n",
"b = a.reshape((2, 8))\n",
"print('a:')\n",
"print(a)\n",
@@ -414,7 +467,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (4, 4))\n",
+ "a = np.arange(16).reshape((4, 4))\n",
"b = np.array(a.reshape(2, 8))\n",
"a[3, 3] = 12345\n",
"b[0, 7] = 54321\n",
@@ -437,7 +490,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (4, 4))\n",
+ "a = np.arange(16).reshape((4, 4))\n",
"print(a)\n",
"print(a.T)"
]
@@ -447,7 +500,7 @@
"metadata": {},
"source": [
"The `transpose` method allows you to obtain more complicated rearrangements\n",
- "of an array's axes:"
+ "of an N-dimensional array's axes:"
]
},
{
@@ -456,7 +509,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (2, 3, 4))\n",
+ "a = np.arange(24).reshape((2, 3, 4))\n",
"b = a.transpose((2, 0, 1))\n",
"print('a: ', a.shape)\n",
"print(a)\n",
@@ -533,15 +586,53 @@
"print( dstacked)"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"operating-on-arrays\"></a>\n",
+ "## Operating on arrays\n",
+ "\n",
+ "\n",
+ "If you are coming from Matlab, you should read this section as, while many\n",
+ "Numpy operations behave similarly to Matlab, there are a few important\n",
+ "behaviours which differ from what you might expect.\n",
+ "\n",
+ "\n",
+ "<a class=\"anchor\" id=\"scalar-operations\"></a>\n",
+ "### Scalar operations\n",
+ "\n",
+ "\n",
+ "All of the mathematical operators you know and love can be applied to Numpy\n",
+ "arrays:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = np.arange(1, 10).reshape((3, 3))\n",
+ "print('a:')\n",
+ "print(a)\n",
+ "print('a + 2:')\n",
+ "print( a + 2)\n",
+ "print('a * 3:')\n",
+ "print( a * 3)\n",
+ "print('a % 2:')\n",
+ "print( a % 2)"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a class=\"anchor\" id=\"multi-variate-operations\"></a>\n",
- "## Multi-variate operations\n",
+ "### Multi-variate operations\n",
"\n",
"\n",
- "Many operations in Numpy operate on an elementwise basis. For example:"
+ "Many operations in Numpy operate on an element-wise basis. For example:"
]
},
{
@@ -550,8 +641,8 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (5))\n",
- "b = np.random.randint(1, 10, (5))\n",
+ "a = np.ones(5)\n",
+ "b = np.random.randint(1, 10, 5)\n",
"\n",
"print('a: ', a)\n",
"print('b: ', b)\n",
@@ -572,8 +663,8 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (4, 4))\n",
- "b = np.random.randint(1, 10, (4, 4))\n",
+ "a = np.ones((4, 4))\n",
+ "b = np.arange(16).reshape((4, 4))\n",
"\n",
"print('a:')\n",
"print(a)\n",
@@ -593,15 +684,17 @@
"Wait ... what's that you say? Oh, I couldn't understand because of all the\n",
"froth coming out of your mouth. I guess you're angry that `a * b` didn't give\n",
"you the matrix product, like it would have in Matlab. Well all I can say is\n",
- "that Python is not Matlab. Get over it. Take a calmative.\n",
+ "that Numpy is not Matlab. Matlab operations are typically consistent with\n",
+ "linear algebra notation. This is not the case in Numpy. Get over it. Take a\n",
+ "calmative.\n",
"\n",
"\n",
"<a class=\"anchor\" id=\"matrix-multiplication\"></a>\n",
- "*## Matrix multiplication\n",
+ "### Matrix multiplication\n",
"\n",
"\n",
"When your heart rate has returned to its normal caffeine-induced state, you\n",
- "can use the `dot` method, or the `@` operator, to perform matrix\n",
+ "can use the `@` operator or the `dot` method, to perform matrix\n",
"multiplication:"
]
},
@@ -611,8 +704,8 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (4, 4))\n",
- "b = np.random.randint(1, 10, (4, 4))\n",
+ "a = np.arange(1, 5).reshape((2, 2))\n",
+ "b = a.T\n",
"\n",
"print('a:')\n",
"print(a)\n",
@@ -633,24 +726,136 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "> The `@` matrix multiplication operator is a relatively recent addition\n",
- "> to Python and Numpy, so you might not see it all that often in existing\n",
- "> code. But it's here to stay, so go ahead and use it!\n",
+ "> The `@` matrix multiplication operator is a relatively recent addition to\n",
+ "> Python and Numpy, so you might not see it all that often in existing\n",
+ "> code. But it's here to stay, so if you don't need to worry about\n",
+ "> backwards-compatibility, go ahead and use it!\n",
"\n",
"\n",
"<a class=\"anchor\" id=\"broadcasting\"></a>\n",
"### Broadcasting\n",
"\n",
"\n",
- "One of the coolest (and possibly confusing) features of Numpy is its\n",
- "[_broadcasting_\n",
- "rules](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html).\n",
+ "One of the coolest features of Numpy is _broadcasting_<sup>3</sup>.\n",
+ "Broadcasting allows you to perform element-wise operations on arrays which\n",
+ "have a different shape. For each axis in the two arrays, Numpy will implicitly\n",
+ "expand the shape of the smaller axis to match the shape of the larger one. You\n",
+ "never need to use `repmat` ever again!\n",
"\n",
"\n",
+ "> <sup>3</sup>Mathworks have shamelessly stolen Numpy's broadcasting behaviour\n",
+ "> and included it in Matlab versions from 2016b onwards, referring to it as\n",
+ "> _implicit expansion_.\n",
+ "\n",
+ "\n",
+ "Broadcasting allows you to, for example, add the elements of a 1D vector to\n",
+ "all of the rows or columns of a 2D array:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = np.arange(9).reshape((3, 3))\n",
+ "b = np.arange(1, 4)\n",
+ "print('a:')\n",
+ "print(a)\n",
+ "print('b: ', b)\n",
+ "print('a * b (row-wise broadcasting):')\n",
+ "print(a * b)\n",
+ "print('a * b.T (column-wise broadcasting):')\n",
+ "print(a * b.reshape(-1, 1))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> Here we used a handy feature of the `reshape` method - if you pass `-1` for\n",
+ "> the size of one dimension, it will automatically determine the size to use\n",
+ "> for that dimension. Take a look at [the\n",
+ "> appendix](#appendix-vectors-in-numpy) for a discussion on vectors in Numpy.\n",
+ "\n",
+ "\n",
+ "Here is a more useful example, where we use broadcasting to de-mean the rows\n",
+ "or columns of an array:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = np.arange(9).reshape((3, 3))\n",
+ "print('a:')\n",
+ "print(a)\n",
+ "\n",
+ "print('a (cols demeaned):')\n",
+ "print(a - a.mean(axis=0))\n",
+ "\n",
+ "print('a (rows demeaned):')\n",
+ "print(a - a.mean(axis=1).reshape(-1, 1))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> As demonstrated above, many functions in Numpy accept an `axis` parameter,\n",
+ "> allowing you to apply the function along a specific axis. Omitting the\n",
+ "> `axis` parameter will apply the function to the whole array.\n",
+ "\n",
+ "\n",
+ "Broadcasting can sometimes be confusing, but the rules which Numpy follows to\n",
+ "align arrays of different sizes, and hence determine how the broadcasting\n",
+ "should be applied, are pretty straightforward. If something is not working,\n",
+ "and you can't figure out why refer to the [official\n",
+ "documentation](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html).\n",
+ "\n",
+ "\n",
+ "<a class=\"anchor\" id=\"linear-algebra\"></a>\n",
+ "### Linear algebra\n",
+ "\n",
+ "\n",
+ "Numpy is first and foremost a library for general-purpose numerical computing.\n",
+ "But it does have a range of linear algebra functionality, hidden away in the\n",
+ "[`numpy.linalg`](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html)\n",
+ "module. Here are a couple of quick examples:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy.linalg as npla\n",
+ "\n",
+ "a = np.array([[1, 2, 3, 4],\n",
+ " [0, 5, 6, 7],\n",
+ " [0, 0, 8, 9],\n",
+ " [0, 0, 0, 10]])\n",
"\n",
+ "print('a:')\n",
+ "print(a)\n",
"\n",
+ "print('inv(a)')\n",
+ "print(npla.inv(a))\n",
"\n",
+ "eigvals, eigvecs = npla.eig(a)\n",
"\n",
+ "print('eigenvalues and vectors of a:')\n",
+ "for val, vec in zip(eigvals, eigvecs):\n",
+ " print('{:2.0f} - {}'.format(val, vec))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
"<a class=\"anchor\" id=\"array-indexing\"></a>\n",
"## Array indexing\n",
"\n",
@@ -666,8 +871,10 @@
"\n",
"\n",
"Let's whet our appetites with some basic 1D array slicing. Numpy supports the\n",
- "standard Python __slice__ notation for indexing, where you can specify the\n",
- "start and end indices, and the step size, via the `start:stop:step` syntax:"
+ "standard Python\n",
+ "[__slice__](https://www.pythoncentral.io/how-to-slice-listsarrays-and-tuples-in-python/)\n",
+ "notation for indexing, where you can specify the start and end indices, and\n",
+ "the step size, via the `start:stop:step` syntax:"
]
},
{
@@ -676,7 +883,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, 10)\n",
+ "a = np.arange(10)\n",
"\n",
"print('a: ', a)\n",
"print('first element: ', a[0])\n",
@@ -685,8 +892,8 @@
"print('last element again: ', a[-1])\n",
"print('last two elements: ', a[-2:])\n",
"print('middle four elements: ', a[3:7])\n",
- "print('Every second element: ', a[::2])\n",
- "print('Every second element, reversed: ', a[1::-2])"
+ "print('Every second element: ', a[1::2])\n",
+ "print('Every second element, reversed: ', a[-1::-2])"
]
},
{
@@ -703,7 +910,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, 10)\n",
+ "a = np.arange(10)\n",
"print('a:', a)\n",
"every2nd = a[::2]\n",
"print('every 2nd:', every2nd)\n",
@@ -720,7 +927,8 @@
"\n",
"\n",
"Multi-dimensional array indexing works in much the same way as one-dimensional\n",
- "indexing but with, well, more dimensions:"
+ "indexing but with, well, more dimensions. Use commas within the square\n",
+ "brackets to separate the slices for each dimension:"
]
},
{
@@ -729,7 +937,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (5, 5))\n",
+ "a = np.arange(25).reshape((5, 5))\n",
"print('a:')\n",
"print(a)\n",
"print(' First row: ', a[ 0, :])\n",
@@ -754,7 +962,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (3, 3, 3))\n",
+ "a = np.arange(27).reshape((3, 3, 3))\n",
"print('a:')\n",
"print(a)\n",
"print('All elements at x=0:')\n",
@@ -783,7 +991,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, 10)\n",
+ "a = np.arange(10)\n",
"\n",
"print('a: ', a)\n",
"print('a > 5: ', a > 4)\n",
@@ -804,7 +1012,7 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, 10)\n",
+ "a = np.arange(10)\n",
"b = a[a > 5]\n",
"print('a: ', a)\n",
"print('b: ', b)\n",
@@ -837,19 +1045,43 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, 10)\n",
+ "a = np.arange(10)\n",
"gt5 = a > 5\n",
"even = a % 2 == 0\n",
"\n",
"print('a: ', a)\n",
- "print('elements in a which are > 5: ', a[gt5])\n",
+ "print('elements in a which are > 5: ', a[ gt5])\n",
"print('elements in a which are <= 5: ', a[~gt5])\n",
- "print('elements in a which are even: ', a[even])\n",
+ "print('elements in a which are even: ', a[ even])\n",
"print('elements in a which are odd: ', a[~even])\n",
"print('elements in a which are > 5 and even: ', a[gt5 & even])\n",
"print('elements in a which are > 5 or odd: ', a[gt5 | ~even])"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Numpy also has two handy functions, `all` and `any`, which respectively allow\n",
+ "you to perform boolean `and` and `or` operations along the axes of an array:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = np.arange(9).reshape((3, 3))\n",
+ "\n",
+ "print('a:')\n",
+ "print(a)\n",
+ "print('rows with any element divisible by 3: ', np.any(a % 3 == 0, axis=1))\n",
+ "print('cols with any element divisible by 3: ', np.any(a % 3 == 0, axis=0))\n",
+ "print('rows with all elements divisible by 3:', np.all(a % 3 == 0, axis=1))\n",
+ "print('cols with all elements divisible by 3:', np.all(a % 3 == 0, axis=0))"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
@@ -861,7 +1093,7 @@
"You can index a numpy array using another array containing coordinates into\n",
"the first array. As with boolean indexing, this will result in a copy of the\n",
"data. Generally, you will need to have a separate array, or list, of\n",
- "coordinates into each data axis:"
+ "coordinates for each axis of the array you wish to index:"
]
},
{
@@ -870,7 +1102,8 @@
"metadata": {},
"outputs": [],
"source": [
- "a = np.random.randint(1, 10, (4, 4))\n",
+ "a = np.arange(16).reshape((4, 4))\n",
+ "print('a:')\n",
"print(a)\n",
"\n",
"rows = [0, 2, 3]\n",
@@ -885,12 +1118,167 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "<a class=\"anchor\" id=\"generating-random-numbers\"></a>\n",
- "## Generating random numbers\n",
+ "The `numpy.where` function can be combined with boolean arrays to easily\n",
+ "generate of coordinate arrays for values which meet some condition:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = np.arange(16).reshape((4, 4))\n",
+ "print('a:')\n",
+ "print(a)\n",
+ "\n",
+ "evenrows, evencols = np.where(a % 2 == 0)\n",
+ "\n",
+ "print('even row coordinates:', evenx)\n",
+ "print('even col coordinates:', eveny)\n",
+ "\n",
+ "print(a[evenrows, evencols])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"exercises\"></a>\n",
+ "## Exercises\n",
+ "\n",
+ "\n",
+ "The challenge for each of these exercises is to complete them in as few lines\n",
+ "of code as possible!\n",
+ "\n",
+ "\n",
+ "> You can find example answers to the exercises in `04_numpy/.solutions`.\n",
+ "\n",
+ "\n",
+ "<a class=\"anchor\" id=\"load-an-array-from-a-file-and-do-stuff-with-it\"></a>\n",
+ "### Load an array from a file and do stuff with it\n",
+ "\n",
+ "\n",
+ "Load the file `04_numpy/2d_array.txt`, and calculate and print the mean for\n",
+ "each column. If your code doesn't work, you might want to __LOOK AT YOUR\n",
+ "DATA__, as you will have learnt if you have ever attended the FSL course.\n",
+ "\n",
+ "\n",
+ "> Bonus: Find the hidden message (hint:\n",
+ "> [chr](https://docs.python.org/3/library/functions.html#chr))\n",
+ "\n",
+ "\n",
+ "<a class=\"anchor\" id=\"concatenate-affine-transforms\"></a>\n",
+ "### Concatenate affine transforms\n",
+ "\n",
+ "\n",
+ "Given all of the files in `04_numpy/xfms/`, create a transformation matrix\n",
+ "which can transform coordinates from subject 1 functional space to subject 2\n",
+ "functional space<sup>4</sup>.\n",
+ "\n",
+ "Write some code to transform the following coordinates from subject 1\n",
+ "functional space to subject 2 functional space, to test that your matrix is\n",
+ "correct:\n",
+ "\n",
+ "\n",
+ "| __Subject 1 coordinates__ | __Subject 2 coordinates__ |\n",
+ "|:-------------------------:|:-------------------------:|\n",
+ "| `[ 0.0, 0.0, 0.0]` | `[ 3.21, 4.15, -9.89]` |\n",
+ "| `[-5.0, -20.0, 10.0]` | `[-0.87, -14.36, -1.13]` |\n",
+ "| `[20.0, 25.0, 60.0]` | `[29.58, 27.61, 53.37]` |\n",
+ "\n",
+ "\n",
+ "> <sup>4</sup> Even though these are FLIRT transforms, this is just a toy\n",
+ "> example. Look\n",
+ "> [here](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the_format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to_the_transformation_parameters.3F)\n",
+ "> and\n",
+ "> [here](https://git.fmrib.ox.ac.uk/fsl/fslpy/blob/1.6.2/fsl/utils/transform.py#L537)\n",
+ "> if you actually need to work with FLIRT transforms.\n",
+ "\n",
+ "\n",
+ "\n",
+ "<a class=\"anchor\" id=\"appendix-generating-random-numbers\"></a>\n",
+ "## Appendix A: Generating random numbers\n",
+ "\n",
+ "\n",
+ "Numpy's\n",
+ "[`numpy.random`](https://docs.scipy.org/doc/numpy/reference/routines.random.html)\n",
+ "module is where you should go if you want to introduce a little randomness\n",
+ "into your code. You have already seen a couple of functions for generating\n",
+ "uniformly distributed real or integer data:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy.random as npr\n",
+ "\n",
+ "print('Random floats between 0.0 (inclusive) and 1.0 (exclusive):')\n",
+ "print(npr.random((3, 3)))\n",
+ "\n",
+ "print('Random integers in a specified range:')\n",
+ "print(npr.randint(1, 100, (3, 3)))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You can also draw random data from other distributions - here are just a few\n",
+ "examples:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print('Gaussian (mean: 0, stddev: 1):')\n",
+ "print(npr.normal(0, 1, (3, 3)))\n",
"\n",
+ "print('Gamma (shape: 1, scale: 1):')\n",
+ "print(npr.normal(1, 1, (3, 3)))\n",
"\n",
+ "print('Chi-square (dof: 10):')\n",
+ "print(npr.chisquare(10, (3, 3)))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `numpy.random` module also has a couple of other handy functions for\n",
+ "random sampling of existing data:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "data = np.arange(5)\n",
+ "\n",
+ "print('data: ', data)\n",
+ "print('two random values: ', npr.choice(data, 2))\n",
+ "print('random permutation: ', npr.permutation(data))\n",
+ "\n",
+ "# The numpy.random.shuffle function\n",
+ "# will shuffle an array *in-place*.\n",
+ "npr.shuffle(data)\n",
+ "print('randomly shuffled: ', data)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
"<a class=\"anchor\" id=\"appendix-importing-numpy\"></a>\n",
- "## Appendix: Importing Numpy\n",
+ "## Appendix B: Importing Numpy\n",
"\n",
"\n",
"For interactive exploration/experimentation, you might want to import\n",
@@ -934,7 +1322,7 @@
"metadata": {},
"source": [
"But if you are writing a script or application using Numpy, I implore you to\n",
- "Numpy like this instead:"
+ "Numpy (and its commonly used sub-modules) like this instead:"
]
},
{
@@ -943,7 +1331,9 @@
"metadata": {},
"outputs": [],
"source": [
- "import numpy as np"
+ "import numpy as np\n",
+ "import numpy.random as npr\n",
+ "import numpy.linalg as npla"
]
},
{
@@ -963,10 +1353,12 @@
"e = np.array([1, 2, 3, 4, 5])\n",
"z = np.zeros((100, 100))\n",
"d = np.diag([2, 3, 4, 5])\n",
+ "r = npr.random(5)\n",
"\n",
"print(e)\n",
"print(z)\n",
- "print(d)"
+ "print(d)\n",
+ "print(r)"
]
},
{
@@ -980,7 +1372,7 @@
"\n",
"\n",
"<a class=\"anchor\" id=\"appendix-vectors-in-numpy\"></a>\n",
- "## Appendix: Vectors in Numpy\n",
+ "## Appendix C: Vectors in Numpy\n",
"\n",
"\n",
"One aspect of Numpy which might trip you up, and which can be quite\n",
@@ -1016,9 +1408,16 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "> Here we used a handy feature of the `reshape` method - if you pass `-1` for\n",
- "> the size of one dimension, it will automatically determine the size to use\n",
- "> for that dimension."
+ "<a class=\"anchor\" id=\"useful-references\"></a>\n",
+ "## Useful references\n",
+ "\n",
+ "\n",
+ "* [The Numpy manual](https://docs.scipy.org/doc/numpy/)\n",
+ "* [Linear algebra in `numpy.linalg`](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html)\n",
+ "* [Broadcasting in Numpy](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html)\n",
+ "* [Indexing in Numpy](https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html)\n",
+ "* [Random sampling in `numpy.random`](https://docs.scipy.org/doc/numpy/reference/routines.random.html)\n",
+ "* [Python slicing](https://www.pythoncentral.io/how-to-slice-listsarrays-and-tuples-in-python/)"
]
}
],
diff --git a/getting_started/04_numpy.md b/getting_started/04_numpy.md
index 0e128dd..dce3532 100644
--- a/getting_started/04_numpy.md
+++ b/getting_started/04_numpy.md
@@ -38,6 +38,9 @@ out of date, but we will update it for the next release of FSL.
* [Indexing multi-dimensional arrays](#indexing-multi-dimensional-arrays)
* [Boolean indexing](#boolean-indexing)
* [Coordinate array indexing](#coordinate-array-indexing)
+* [Exercises](#exercises)
+ * [Load an array from a file and do stuff with it](#load-an-array-from-a-file-and-do-stuff-with-it)
+ * [Concatenate affine transforms](#concatenate-affine-transforms)
* [Appendix A: Generating random numbers](#appendix-generating-random-numbers)
* [Appendix B: Importing Numpy](#appendix-importing-numpy)
@@ -853,8 +856,11 @@ print(a[evenrows, evencols])
## Exercises
-The challenge for each of these exercises is to complete them with as few
-lines of code as possible (whitespace and comments don't count)!
+The challenge for each of these exercises is to complete them in as few lines
+of code as possible!
+
+
+> You can find example answers to the exercises in `04_numpy/.solutions`.
<a class="anchor" id="load-an-array-from-a-file-and-do-stuff-with-it"></a>
@@ -863,25 +869,31 @@ lines of code as possible (whitespace and comments don't count)!
Load the file `04_numpy/2d_array.txt`, and calculate and print the mean for
each column. If your code doesn't work, you might want to __LOOK AT YOUR
-DATA__, as you will have learnt during the FSL course.
+DATA__, as you will have learnt if you have ever attended the FSL course.
> Bonus: Find the hidden message (hint:
> [chr](https://docs.python.org/3/library/functions.html#chr))
-
<a class="anchor" id="concatenate-affine-transforms"></a>
### Concatenate affine transforms
Given all of the files in `04_numpy/xfms/`, create a transformation matrix
which can transform coordinates from subject 1 functional space to subject 2
-functional space<sup>4></sup>.
+functional space<sup>4</sup>.
+Write some code to transform the following coordinates from subject 1
+functional space to subject 2 functional space, to test that your matrix is
+correct:
-TODO Test coordinates
+| __Subject 1 coordinates__ | __Subject 2 coordinates__ |
+|:-------------------------:|:-------------------------:|
+| `[ 0.0, 0.0, 0.0]` | `[ 3.21, 4.15, -9.89]` |
+| `[-5.0, -20.0, 10.0]` | `[-0.87, -14.36, -1.13]` |
+| `[20.0, 25.0, 60.0]` | `[29.58, 27.61, 53.37]` |
> <sup>4</sup> Even though these are FLIRT transforms, this is just a toy
@@ -892,8 +904,6 @@ TODO Test coordinates
> if you actually need to work with FLIRT transforms.
-> You can find example answers to the exercises in `04_numpy/.solutions`.
-
<a class="anchor" id="appendix-generating-random-numbers"></a>
## Appendix A: Generating random numbers
--
GitLab