From ed9e7f41a816ed372c43551f780a97742458bf8a Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Thu, 1 Feb 2018 12:37:07 +0000
Subject: [PATCH] Nearly there with numpy prac
---
getting_started/04_numpy.md | 190 ++++++++++++++++--
.../04_numpy/.solutions/calc_column_mean.py | 10 +
.../04_numpy/.solutions/concat_affines.py | 50 +++++
getting_started/04_numpy/04_numpy.ipynb | 0
getting_started/04_numpy/2d_array.txt | 13 ++
getting_started/04_numpy/comma_separated.txt | 10 +
getting_started/04_numpy/space_separated.txt | 20 ++
.../tsetest/subj1/example_func2highres.mat | 4 +
.../tsetest/subj1/highres2standard.mat | 4 +
.../tsetest/subj2/example_func2highres.mat | 4 +
.../tsetest/subj2/highres2standard.mat | 4 +
.../xfms/subj1/example_func2highres.mat | 4 +
.../04_numpy/xfms/subj1/highres2standard.mat | 4 +
.../xfms/subj2/example_func2highres.mat | 4 +
.../04_numpy/xfms/subj2/highres2standard.mat | 4 +
15 files changed, 307 insertions(+), 18 deletions(-)
create mode 100644 getting_started/04_numpy/.solutions/calc_column_mean.py
create mode 100644 getting_started/04_numpy/.solutions/concat_affines.py
create mode 100644 getting_started/04_numpy/04_numpy.ipynb
create mode 100644 getting_started/04_numpy/2d_array.txt
create mode 100644 getting_started/04_numpy/comma_separated.txt
create mode 100644 getting_started/04_numpy/space_separated.txt
create mode 100755 getting_started/04_numpy/tsetest/subj1/example_func2highres.mat
create mode 100755 getting_started/04_numpy/tsetest/subj1/highres2standard.mat
create mode 100755 getting_started/04_numpy/tsetest/subj2/example_func2highres.mat
create mode 100755 getting_started/04_numpy/tsetest/subj2/highres2standard.mat
create mode 100755 getting_started/04_numpy/xfms/subj1/example_func2highres.mat
create mode 100755 getting_started/04_numpy/xfms/subj1/highres2standard.mat
create mode 100755 getting_started/04_numpy/xfms/subj2/example_func2highres.mat
create mode 100755 getting_started/04_numpy/xfms/subj2/highres2standard.mat
diff --git a/getting_started/04_numpy.md b/getting_started/04_numpy.md
index 4133268..0e128dd 100644
--- a/getting_started/04_numpy.md
+++ b/getting_started/04_numpy.md
@@ -13,12 +13,18 @@ important Python libraries, and it (along with its partners
alternative to Matlab as a scientific computing platform.
+The `fslpython` environment currently includes [Numpy
+1.11.1](https://docs.scipy.org/doc/numpy-1.11.0/index.html), which is a little
+out of date, but we will update it for the next release of FSL.
+
+
## Contents
* [The Python list versus the Numpy array](#the-python-list-versus-the-numpy-array)
* [Numpy basics](#numpy-basics)
* [Creating arrays](#creating-arrays)
+ * [Loading text files](#loading-text-files)
* [Array properties](#array-properties)
* [Descriptive statistics](#descriptive-statistics)
* [Reshaping and rearranging arrays](#reshaping-and-rearranging-arrays)
@@ -84,7 +90,7 @@ are trying to write efficient code. It is _crucial_ to be able to distinguish
between a Python list and a Numpy array.
-___Python list == Matlab cell array:___ A list in python is akin to a cell
+___Python list == Matlab cell array:___ A list in Python is akin to a cell
array in Matlab - they can store anything, but are extremely inefficient, and
unwieldy when you have more than a couple of dimensions.
@@ -196,6 +202,52 @@ print(o)
> second to columns - just like in Matlab.
+<a class="anchor" id="loading-text-files"></a>
+### Loading text files
+
+
+The `numpy.loadtxt` function is capable of loading numerical data from
+plain-text files. By default it expects space-separated data:
+
+
+```
+data = np.loadtxt('04_numpy/space_separated.txt')
+print('data in 04_numpy/space_separated.txt:')
+print(data)
+```
+
+
+But you can also specify the delimiter to expect<sup>1</sup>:
+
+
+```
+data = np.loadtxt('04_numpy/comma_separated.txt', delimiter=',')
+print('data in 04_numpy/comma_separated.txt:')
+print(data)
+```
+
+
+> <sup>1</sup> And many other things such as file headers, footers, comments,
+> and newline characters - see the
+> [docs](https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html)
+> for more information.
+
+
+ Of course you can also save data out to a text file just as easily, with
+[`numpy.savetxt`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html):
+
+
+```
+data = np.random.randint(1, 10, (10, 10))
+np.savetxt('mydata.txt', data, delimiter=',', fmt='%i')
+
+with open('mydata.txt', 'rt') as f:
+ for line in f:
+ print(line.strip())
+```
+
+
+
<a class="anchor" id="array-properties"></a>
### Array properties
@@ -227,8 +279,8 @@ print('Length of first dimension: ', len(z))
### Descriptive statistics
-Similarly, a Numpy array has a set of methods<sup>1</sup> which allow you to
-calculate basic descriptive statisics on an array:
+Similarly, a Numpy array has a set of methods<sup>2</sup> which allow you to
+calculate basic descriptive statistics on an array:
```
@@ -246,11 +298,12 @@ print('prod: ', a.prod())
```
-> <sup>1</sup> Python, being an object-oriented language, distinguishes
-> between _functions_ and _methods_. _Method_ is simply the term used to refer
-> to a function that is associated with a specific object. Similarly, the term
-> _attribute_ is used to refer to some piece of information that is attached
-> to an object, such as `z.shape`, or `z.dtype`.
+> <sup>2</sup> Python, being an object-oriented language, distinguishes
+> between _functions_ and _methods_. Hopefully we all know what a function is
+> - a _method_ is simply the term used to refer to a function that is
+> associated with a specific object. Similarly, the term _attribute_ is used
+> to refer to some piece of information that is attached to an object, such as
+> `z.shape`, or `z.dtype`.
<a class="anchor" id="reshaping-and-rearranging-arrays"></a>
@@ -489,14 +542,14 @@ print(b.dot(a))
### Broadcasting
-One of the coolest features of Numpy is _broadcasting_<sup>2</sup>.
+One of the coolest features of Numpy is _broadcasting_<sup>3</sup>.
Broadcasting allows you to perform element-wise operations on arrays which
have a different shape. For each axis in the two arrays, Numpy will implicitly
expand the shape of the smaller axis to match the shape of the larger one. You
never need to use `repmat` ever again!
-> <sup>2</sup>Mathworks have shamelessly stolen Numpy's broadcasting behaviour
+> <sup>3</sup>Mathworks have shamelessly stolen Numpy's broadcasting behaviour
> and included it in Matlab versions from 2016b onwards, referring to it as
> _implicit expansion_.
@@ -640,7 +693,8 @@ print('a:', a)
Multi-dimensional array indexing works in much the same way as one-dimensional
-indexing but with, well, more dimensions:
+indexing but with, well, more dimensions. Use commas within the square
+brackets to separate the slices for each dimension:
```
@@ -785,19 +839,118 @@ a = np.arange(16).reshape((4, 4))
print('a:')
print(a)
-evenx, eveny = np.where(a % 2 == 0)
+evenrows, evencols = np.where(a % 2 == 0)
-print('even X coordinates:', evenx)
-print('even Y coordinates:', eveny)
+print('even row coordinates:', evenx)
+print('even col coordinates:', eveny)
-print(a[evenx, eveny])
+print(a[evenrows, evencols])
```
+
+<a class="anchor" id="exercises"></a>
+## 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)!
+
+
+<a class="anchor" id="load-an-array-from-a-file-and-do-stuff-with-it"></a>
+### Load an array from a file and do stuff with it
+
+
+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.
+
+
+> 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>.
+
+
+TODO Test coordinates
+
+
+
+> <sup>4</sup> Even though these are FLIRT transforms, this is just a toy
+> example. Look
+> [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)
+> and
+> [here](https://git.fmrib.ox.ac.uk/fsl/fslpy/blob/1.6.2/fsl/utils/transform.py#L537)
+> 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
+Numpy's
+[`numpy.random`](https://docs.scipy.org/doc/numpy/reference/routines.random.html)
+module is where you should go if you want to introduce a little randomness
+into your code. You have already seen a couple of functions for generating
+uniformly distributed real or integer data:
+
+
+```
+import numpy.random as npr
+
+print('Random floats between 0.0 (inclusive) and 1.0 (exclusive):')
+print(npr.random((3, 3)))
+
+print('Random integers in a specified range:')
+print(npr.randint(1, 100, (3, 3)))
+```
+
+
+You can also draw random data from other distributions - here are just a few
+examples:
+
+
+```
+print('Gaussian (mean: 0, stddev: 1):')
+print(npr.normal(0, 1, (3, 3)))
+
+print('Gamma (shape: 1, scale: 1):')
+print(npr.normal(1, 1, (3, 3)))
+
+print('Chi-square (dof: 10):')
+print(npr.chisquare(10, (3, 3)))
+```
+
+
+The `numpy.random` module also has a couple of other handy functions for
+random sampling of existing data:
+
+
+```
+data = np.arange(5)
+
+print('data: ', data)
+print('two random values: ', npr.choice(data, 2))
+print('random permutation: ', npr.permutation(data))
+
+# The numpy.random.shuffle function
+# will shuffle an array *in-place*.
+npr.shuffle(data)
+print('randomly shuffled: ', data)
+```
+
+
<a class="anchor" id="appendix-importing-numpy"></a>
## Appendix B: Importing Numpy
@@ -894,7 +1047,8 @@ print(np.atleast_2d(r).T)
* [The Numpy manual](https://docs.scipy.org/doc/numpy/)
-* [Linear algebra with `numpy.linalg`](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html)
-* [Numpy broadcasting](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html)
-* [Numpy indexing](https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html)
+* [Linear algebra in `numpy.linalg`](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html)
+* [Broadcasting in Numpy](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html)
+* [Indexing in Numpy](https://docs.scipy.org/doc/numpy/reference/arrays.indexing.html)
+* [Random sampling in `numpy.random`](https://docs.scipy.org/doc/numpy/reference/routines.random.html)
* [Python slicing](https://www.pythoncentral.io/how-to-slice-listsarrays-and-tuples-in-python/)
\ No newline at end of file
diff --git a/getting_started/04_numpy/.solutions/calc_column_mean.py b/getting_started/04_numpy/.solutions/calc_column_mean.py
new file mode 100644
index 0000000..91abf3c
--- /dev/null
+++ b/getting_started/04_numpy/.solutions/calc_column_mean.py
@@ -0,0 +1,10 @@
+#!/usr/bin/env python
+import numpy as np
+
+data = np.loadtxt('04_numpy/2d_array.txt', comments='%')
+colmeans = data.mean(axis=0)
+msg = ''.join([chr(int(round(m))) for m in colmeans])
+
+print('Column means')
+print('\n'.join(['{}: {}'.format(i, m) for i, m in enumerate(colmeans)]))
+print('Secret message:', msg)
diff --git a/getting_started/04_numpy/.solutions/concat_affines.py b/getting_started/04_numpy/.solutions/concat_affines.py
new file mode 100644
index 0000000..92d8a49
--- /dev/null
+++ b/getting_started/04_numpy/.solutions/concat_affines.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+
+import numpy as np
+import numpy.linalg as npla
+
+def concat(*xforms):
+ """Combines the given matrices (returns the dot product)."""
+
+ result = xforms[0]
+
+ for i in range(1, len(xforms)):
+ result = np.dot(result, xforms[i])
+
+ return result
+
+
+def transform(xform, coord):
+ """Transform the given coordinates with the given affine. """
+ return np.dot(xform[:3, :3], coord) + xform[:3, 3]
+
+
+s1func2struc = np.loadtxt('04_numpy/xfms/subj1/example_func2highres.mat')
+s1struc2std = np.loadtxt('04_numpy/xfms/subj1/highres2standard.mat')
+s2func2struc = np.loadtxt('04_numpy/xfms/subj2/example_func2highres.mat')
+s2struc2std = np.loadtxt('04_numpy/xfms/subj2/highres2standard.mat')
+
+s1func2std = concat(s1struc2std, s1func2struc)
+s2func2std = concat(s2struc2std, s2func2struc)
+s2std2func = npla.inv(s2func2std)
+s1func2s2func = concat(s2std2func, s1func2std)
+
+print('Subject 1 functional -> Subject 2 functional affine:')
+print(s1func2s2func)
+
+testcoords = np.array([[ 0, 0, 0],
+ [-5, -20, 10],
+ [20, 25, 60]], dtype=np.float32)
+
+
+def dist(p1, p2):
+ return np.sqrt(np.sum((p1 - p2) ** 2))
+
+
+
+for c in testcoords:
+ xc = transform(s1func2s2func, c)
+ d = dist(c, xc)
+ c = '{:6.2f} {:6.2f} {:6.2f}'.format(*c)
+ xc = '{:6.2f} {:6.2f} {:6.2f}'.format(*xc)
+ print('Transform: [{}] -> [{}] (distance: {})'.format(c, xc, d))
diff --git a/getting_started/04_numpy/04_numpy.ipynb b/getting_started/04_numpy/04_numpy.ipynb
new file mode 100644
index 0000000..e69de29
diff --git a/getting_started/04_numpy/2d_array.txt b/getting_started/04_numpy/2d_array.txt
new file mode 100644
index 0000000..923f35e
--- /dev/null
+++ b/getting_started/04_numpy/2d_array.txt
@@ -0,0 +1,13 @@
+% The purpose of this comment is purely
+% to make life difficult for you. Sorry
+% about that. I'm not really sorry.
+128.3 100.8 67.2 120.1 150.2 53.0 64.2 139.3 46.7 118.1 125.8 153.1 83.2 115.9 3.4 126.3 92.2 104.7 131.2 29.3
+ 89.3 118.8 119.2 80.1 121.2 35.0 66.2 153.3 43.7 102.1 147.8 160.1 94.2 140.9 70.4 124.3 124.2 93.7 100.2 32.3
+ 91.3 146.8 137.2 129.1 157.2 -22.0 75.2 90.3 33.7 86.1 74.8 90.1 114.2 81.9 72.4 135.3 85.2 102.7 58.2 55.3
+ 75.3 87.8 84.2 112.1 131.2 57.0 150.2 73.3 41.7 110.1 78.8 102.1 116.2 134.9 34.4 74.3 150.2 52.7 89.2 21.3
+ 79.3 83.8 78.2 103.1 96.2 26.0 151.2 122.3 46.7 113.1 143.8 80.1 61.2 136.9 8.4 94.3 76.2 123.7 150.2 10.3
+ 93.3 157.8 154.2 79.1 155.2 8.0 63.2 120.3 12.7 130.1 105.8 111.1 66.2 99.9 54.4 120.3 113.2 128.7 110.2 46.3
+111.3 126.8 65.2 90.1 101.2 36.0 116.2 139.3 36.7 91.1 88.8 94.1 90.2 92.9 -0.6 46.3 80.2 105.7 88.2 -11.7
+149.3 119.8 134.2 124.1 89.2 21.0 133.2 113.3 28.7 139.1 136.8 82.1 141.2 90.9 11.4 97.3 109.2 137.7 120.2 32.3
+141.3 149.8 115.2 115.1 105.2 51.0 77.2 130.3 -4.3 124.1 141.8 162.1 95.2 130.9 -10.6 114.3 132.2 134.7 79.2 69.3
+141.3 77.8 135.2 167.1 103.2 55.0 153.2 68.3 33.7 136.1 125.8 85.1 148.2 114.9 76.4 57.3 147.2 125.7 153.2 45.3
diff --git a/getting_started/04_numpy/comma_separated.txt b/getting_started/04_numpy/comma_separated.txt
new file mode 100644
index 0000000..9aa1ff6
--- /dev/null
+++ b/getting_started/04_numpy/comma_separated.txt
@@ -0,0 +1,10 @@
+13,7,3,18,15,3,12,11,9,9
+13,1,7,17,16,13,18,9,18,6
+9,19,16,3,18,3,19,12,9,6
+2,11,6,12,2,11,15,9,3,9
+2,12,1,7,4,3,6,6,2,4
+10,8,14,1,17,19,8,19,2,9
+18,14,2,11,17,14,6,16,14,18
+6,8,13,16,11,17,16,5,16,15
+6,5,4,18,6,14,19,8,4,15
+15,17,12,10,17,12,5,9,18,6
diff --git a/getting_started/04_numpy/space_separated.txt b/getting_started/04_numpy/space_separated.txt
new file mode 100644
index 0000000..f5fce0d
--- /dev/null
+++ b/getting_started/04_numpy/space_separated.txt
@@ -0,0 +1,20 @@
+1 3
+14 3
+4 7
+12 5
+15 5
+14 4
+15 19
+7 17
+7 15
+16 3
+15 14
+17 9
+7 8
+4 5
+15 4
+11 5
+17 10
+3 19
+14 4
+18 3
diff --git a/getting_started/04_numpy/tsetest/subj1/example_func2highres.mat b/getting_started/04_numpy/tsetest/subj1/example_func2highres.mat
new file mode 100755
index 0000000..c5f9b62
--- /dev/null
+++ b/getting_started/04_numpy/tsetest/subj1/example_func2highres.mat
@@ -0,0 +1,4 @@
+1.027343449 -0.01854307437 0.02702291658 -3.429567172
+0.01572791398 0.9487250481 -0.30721057 17.65058627
+-0.02892982284 0.2868726435 1.077910798 -3.247641029
+0 0 0 1
diff --git a/getting_started/04_numpy/tsetest/subj1/highres2standard.mat b/getting_started/04_numpy/tsetest/subj1/highres2standard.mat
new file mode 100755
index 0000000..e11d282
--- /dev/null
+++ b/getting_started/04_numpy/tsetest/subj1/highres2standard.mat
@@ -0,0 +1,4 @@
+0.9998006214 -0.007506170473 -0.00206080048 -37.86830353
+0.01249141095 0.863260387 0.4957794812 -33.77969415
+0.03530937349 -0.4932384691 1.042606617 27.59239864
+0 0 0 1
diff --git a/getting_started/04_numpy/tsetest/subj2/example_func2highres.mat b/getting_started/04_numpy/tsetest/subj2/example_func2highres.mat
new file mode 100755
index 0000000..101bf79
--- /dev/null
+++ b/getting_started/04_numpy/tsetest/subj2/example_func2highres.mat
@@ -0,0 +1,4 @@
+1.0271367 -0.01315117882 -0.005127727479 -3.084298759
+0.002836833979 0.9879611849 -0.002646076731 -0.4778179727
+0.004268282598 -0.0002074879811 1.067805712 25.30182119
+0 0 0 1
diff --git a/getting_started/04_numpy/tsetest/subj2/highres2standard.mat b/getting_started/04_numpy/tsetest/subj2/highres2standard.mat
new file mode 100755
index 0000000..3f7c3cd
--- /dev/null
+++ b/getting_started/04_numpy/tsetest/subj2/highres2standard.mat
@@ -0,0 +1,4 @@
+1.077742594 -0.04288767681 -0.08886849541 -40.1774847
+0.03613742149 1.017656652 0.2363464233 -27.41273922
+0.07295070007 -0.3129031707 1.14704109 -0.4162223149
+0 0 0 1
diff --git a/getting_started/04_numpy/xfms/subj1/example_func2highres.mat b/getting_started/04_numpy/xfms/subj1/example_func2highres.mat
new file mode 100755
index 0000000..c5f9b62
--- /dev/null
+++ b/getting_started/04_numpy/xfms/subj1/example_func2highres.mat
@@ -0,0 +1,4 @@
+1.027343449 -0.01854307437 0.02702291658 -3.429567172
+0.01572791398 0.9487250481 -0.30721057 17.65058627
+-0.02892982284 0.2868726435 1.077910798 -3.247641029
+0 0 0 1
diff --git a/getting_started/04_numpy/xfms/subj1/highres2standard.mat b/getting_started/04_numpy/xfms/subj1/highres2standard.mat
new file mode 100755
index 0000000..e11d282
--- /dev/null
+++ b/getting_started/04_numpy/xfms/subj1/highres2standard.mat
@@ -0,0 +1,4 @@
+0.9998006214 -0.007506170473 -0.00206080048 -37.86830353
+0.01249141095 0.863260387 0.4957794812 -33.77969415
+0.03530937349 -0.4932384691 1.042606617 27.59239864
+0 0 0 1
diff --git a/getting_started/04_numpy/xfms/subj2/example_func2highres.mat b/getting_started/04_numpy/xfms/subj2/example_func2highres.mat
new file mode 100755
index 0000000..101bf79
--- /dev/null
+++ b/getting_started/04_numpy/xfms/subj2/example_func2highres.mat
@@ -0,0 +1,4 @@
+1.0271367 -0.01315117882 -0.005127727479 -3.084298759
+0.002836833979 0.9879611849 -0.002646076731 -0.4778179727
+0.004268282598 -0.0002074879811 1.067805712 25.30182119
+0 0 0 1
diff --git a/getting_started/04_numpy/xfms/subj2/highres2standard.mat b/getting_started/04_numpy/xfms/subj2/highres2standard.mat
new file mode 100755
index 0000000..3f7c3cd
--- /dev/null
+++ b/getting_started/04_numpy/xfms/subj2/highres2standard.mat
@@ -0,0 +1,4 @@
+1.077742594 -0.04288767681 -0.08886849541 -40.1774847
+0.03613742149 1.017656652 0.2363464233 -27.41273922
+0.07295070007 -0.3129031707 1.14704109 -0.4162223149
+0 0 0 1
--
GitLab