Commit ed9e7f41 authored by Paul McCarthy's avatar Paul McCarthy 🚵
Browse files

Nearly there with numpy prac

parent 27317639
......@@ -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
#!/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)
#!/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))
% 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
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
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
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
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
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
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
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
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
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
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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment