# - compatibility with Python 3
from __future__ import print_function  # print('me') instead of print 'me'
from __future__ import division  # 1/2 == 0.5, not 0

# - show figures inside the notebook
%matplotlib inline

# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package

# - import numpy.linalg with a shorter name
import numpy.linalg as npl

# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'

Load the image ds114_sub009_t2r1.nii with nibabel. Get the data:

# Load the image 'ds114_sub009_t2r1.nii' with nibabel
# Get the data array from the image

We can think of the shape of the data as two parts - the first 3 values are the 3D shape of the individual volumes, and the last value is the number of volumes.

Put the 3D shape into a variable vol_shape and the number of volumes into a variable n_vols:

# Make variables:
# 'vol_shape' for shape of volumes
# 'n_vols' for number of volumes

We are going to start by looking at only the first two "variables" (time points).

Slice the data array to make a new array that contains only the first two volumes:

# Slice the image data array to give array with only first two volumes

How many voxels are there in one volume? Call this number variable N:

# Set N to be the number of voxels in a volume

Reshape the new two-volume data array to have a first dimension length N and second dimension length 2. So, each column corresponds to the voxels for one volume.

# Reshape to 2D array with first dimension length N

Take the transpose of this array to get a 2 by N array, ready for the PCA:

# Transpose to 2 by N array

Calculate the mean across columns (row means):

# Calculate the mean across columns

Make the row_means vector into a column vector, using reshape. (This gets the vector ready for broadcasting):

# Make row means into a column vector

Subtract the row_means column vector from the 2 by N data matrix and put the result into a variable X. Remember, broadcasting will expand the last 1-length dimension of row_means to length N, by repeating the data.

Print the means across columns (row means) to check they are now very close to 0:

# Subtract the means for each row, put the result into X
# Show the means over the columns, after the subtraction

Plot the two rows against each other to get a feel for the variation. Remember that each row in X is a volume, so you are plotting the signal from the first volume against the signal for the second volume.

# Plot the signal in the first row against the signal in the second

Time for the PCA.

We are going to use the unscaled covariance. This is given by the matrix multiplication of X with its transpose. Calculate this:

# Calculate unscaled covariance matrix for X

Use SVD to get the U, S and V matrices from the unscaled covariance:

# Use SVD to return U, S, V matrices from unscaled covariance

Confirm that the row vectors in V are both unit vectors. A unit vector has vector length (vector norm of 1):

# Show that the rows in V each have vector length 1

Confirm that the first row in V is orthogonal to the second:

# Confirm orthogonality of rows in V

Show the total sum of squares in X. Confirm that the total sum of squares in X is the same as the sum of the singular values in the S vector from the SVD:

# Show the total sum of squares in X
# Is this (nearly) the same as the sum of the values in S?

Plot the first row in X against the second row in X again. This time add a line corresponding to the first principal component. You'll have to scale this line (a unit vector) to make it long enough to see on the axes of the plot. If the first principal component is flipped ($\vec{r}$ defines same line as $-\vec{r}$) you may need to scale by a negative number for it to look nice on the plot:

# Plot the signal in the first row against the signal in the second
# Plot line corresponding to a scaled version of the first principal component
# (Scaling may need to be negative)

Remember the projection formula $c = \hat{u} \cdot \vec{v}$.

We now need to calculate the projection coefficients $c$ for each component $\hat{u}$ and each voxel (each $\vec{v}$).

This will give us a new output matrix of coefficients $C$ which is shape (2, N), where the rows give the coefficients for one component, and the columns give the coefficients for one voxel.

For example, C[0, 0] will be the result of V[0].dot(X[:, 0]), C[0, 1] will be the result of V[0].dot(X[:, 1]), and C[1, 0] will be the result of V[1].dot(X[:, 0]).

With that background, see if you can use matrix multiplication to calculate the projection coefficients C for projecting the data X onto the vectors in V:

# Calculate the projection coefficients for projecting X onto the vectors in V
# Put the result into a new array C.

Remember that C - the projection coefficients - has one column per voxel. We can think of each row as corresponding to a volume where the volumes are contained in: C[0] (first row of C) - projection coefficients for first principal component; C[1] (second row of C) - projection coefficients for second principal component.

Take the transpose of C and reshape the resulting first dimension (length N) back to vol_shape - the original shape of the 3D volumes in the data.

# Transpose C
# Reshape the first dimension of C to have the 3D shape of the original data volumes.

Break this 4D array up into two volumes (volume for first component, volume for second component) using slicing:

# Break 4D array into two 3D volumes

Show the middle slice (slice over the third dimension) from the volume of coefficients for the first component:

# Show middle slice (over third dimension) from volume of coefficients for first component

Show the middle slice (slice over the third dimension) from the volume of coefficients for the second component:

# Show middle slice (over third dimension) from volume of coefficients for second component

Now we are ready to do the same thing for all the time points, instead of just the first two.

Take the original array data matrix with n_vols volumes. Reshape to be shape (N, nvols). Take the transpose to get an array shape (n_vols, N):

# Reshape first dimension of whole image data array to N, and take transpose

Calculate the mean across columns (mean per volume). Make this 1D vector (shape (n_vols,)) vector into a column vector (shape (n_vols, 1)). Subtract this array from the (n_vols, N) data array to remove the mean across columns (mean per volume) (broadcasting will extend the single mean column to N repeats for the subtraction). Call this mean-corrected variable X:

# Calculate mean across columns
# Reshape mean to a column vector
# Subtract from data array to remove mean over columns (row means)
# Put result into array X

Get the unscaled covariance matrix of X:

# Calculate unscaled covariance matrix of X

Get U, S, V outputs from the SVD of the unscaled covariance

# Calculate U, S, V with SVD on unscaled covariance matrix

Use plt.subplots to make a column of 10 axes. Iterate over these axes, plotting one principal component vector per axis. You should then have a plot of each of the first 10 principal component vectors:

# Use subplots to make axes to plot first 10 principle component vectors
# Plot one component vector per sub-plot.

Use V and matrix multiplication to calculate the projection coefficients for projecting the data X onto the principal components V. Put the result into C:

# Calculate projection coefficients for projecting X onto V
# Put results into array C.

Remember, each row of C is a full volume of projection coefficients, one row per principal component.

Reconstruct these rows as volumes by taking the transpose of C and reshaping the first dimension length N to the original three dimensions of the original data volumes.

# Transpose C
# Reshape the first dimension of C to have the 3D shape of the original data volumes.

Take the first volume (corresponding to the first principal component) and display the middle slice (slicing over the third dimension):

# Show middle slice (over third dimension) of first principal component volume

What does this first component show us?

Remember that the projections coefficients tell us how much of the component vector there is in the data, at each voxel. If the voxel has a high value it means there is a large amplitude of the component time-course at this voxel, and if the component is dark it means there is a large negative amplitude of the component time course.

Have a look again at the first component time course. How would a large positive or negative amplitude of the time course come about?

As a hint while you are thinking, get the mean over time from the image data (mean over the last axis), and show the middle slice (slicing over the third axis):

# Make the mean volume (mean over the last axis)
# Show the middle slice (slicing over the third axis)

Display the middle slice (slicing over the third dimension) for the second principal component volume.

Looking at the principal component plot - what kind of changes over time does this principal component represent?

# Show middle slice (over third dimension) of second principal component volume

Do the same for the third principal component volume:

# Show middle slice (over third dimension) of third principal component volume