Skip to content
Snippets Groups Projects
Commit 45de929f authored by Sean Fitzgibbon's avatar Sean Fitzgibbon
Browse files

tweaks

parent ee1184f3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
talks/matlab_vs_python/migp/python_MIGP.ipynb
+ 13
− 8
View file @ 45de929f
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Python MIGP ## Python MIGP
This notebook will perform *python* MIGP dimension reduction, run group ICA, and then plot the group ICs. This notebook will perform *python* MIGP dimension reduction, run group ICA, and then plot the group ICs.
* [Run python `MIGP`](#run-python-migp) * [Run python `MIGP`](#run-python-migp)
* [Run `melodic`](#run-python-melodic) * [Run `melodic`](#run-python-melodic)
* [Plot group ICs](#plot-python-group-ics) * [Plot group ICs](#plot-python-group-ics)
Firstly we will import the necessary packages for this notebook: Firstly we will import the necessary packages for this notebook:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import glob import glob
import random import random
import nibabel as nb import nibabel as nb
import numpy as np import numpy as np
from scipy.sparse.linalg import svds, eigs from scipy.sparse.linalg import svds, eigs
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from nilearn import plotting from nilearn import plotting
from nilearn import image from nilearn import image
import os.path as op import os.path as op
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
It will be necessary to know the location where the data was stored so that we can load the brainmask. It will be necessary to know the location where the data was stored so that we can load the brainmask.
> **Note**: [`expanduser`](https://docs.python.org/3.7/library/os.path.html#os.path.expanduser) will expand the `~` to the be users home directory > **Note**: [`expanduser`](https://docs.python.org/3.7/library/os.path.html#os.path.expanduser) will expand the `~` to the be users home directory
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
data_dir = op.expanduser('~/nilearn_data') data_dir = op.expanduser('~/nilearn_data')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<a class="anchor" id="run-python-migp"></a> <a class="anchor" id="run-python-migp"></a>
### Run python `MIGP` ### Run python `MIGP`
Firstly we need to set the MIGP parameters: Firstly we need to set the MIGP parameters:
> **Note:** > **Note:**
> 1. [`glob.glob`](https://docs.python.org/3/library/glob.html) will create a list of filenames that match the glob/wildcard pattern > 1. [`glob.glob`](https://docs.python.org/3/library/glob.html) will create a list of filenames that match the glob/wildcard pattern
> 2. [`nb.load`](https://nipy.org/nibabel/gettingstarted.html) from the [`nibabel`](https://nipy.org/nibabel/index.html) package will load the image into [`nibabel.Nifti1Image`](https://nipy.org/nibabel/reference/nibabel.nifti1.html) object. This will not load the actual data though. > 2. [`nb.load`](https://nipy.org/nibabel/gettingstarted.html) from the [`nibabel`](https://nipy.org/nibabel/index.html) package will load the image into [`nibabel.Nifti1Image`](https://nipy.org/nibabel/reference/nibabel.nifti1.html) object. This will not load the actual data though.
> 3. We use a list comprehension to loop through all the filenames and load them with [`nibabel`](https://nipy.org/nibabel/index.html) > 3. We use a list comprehension to loop through all the filenames and load them with [`nibabel`](https://nipy.org/nibabel/index.html)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# create lists of (nibabel) image objects # create lists of (nibabel) image objects
in_list = [nb.load(f) for f in glob.glob(f'{data_dir}/cobre/fmri_*_smooth.nii.gz')] in_list = [nb.load(f) for f in glob.glob(f'{data_dir}/cobre/fmri_*_smooth.nii.gz')]
in_mask = nb.load(f'{data_dir}/brain_mask.nii.gz') in_mask = nb.load(f'{data_dir}/brain_mask.nii.gz')
# set user parameters (equivalent to melodic defaults) # set user parameters (equivalent to melodic defaults)
GO = 'pyMIGP.nii.gz' # output filename GO = 'pyMIGP.nii.gz' # output filename
dPCA_int = 299 # internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that) dPCA_int = 299 # internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)
dPCA_out = 299 # number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality dPCA_out = 299 # number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality
sep_vn = False # switch on separate variance nomalisation for each input dataset sep_vn = False # switch on separate variance nomalisation for each input dataset
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
> **Note:** > **Note:**
> 1. [`random.shuffle`](https://docs.python.org/3.7/library/random.html#random.shuffle) will shuffle a list, in this instance it shuffles the list of [`nibabel.Nifti1Image`](https://nipy.org/nibabel/reference/nibabel.nifti1.html) objects > 1. [`random.shuffle`](https://docs.python.org/3.7/library/random.html#random.shuffle) will shuffle a list, in this instance it shuffles the list of [`nibabel.Nifti1Image`](https://nipy.org/nibabel/reference/nibabel.nifti1.html) objects
> 2. [`ravel`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ravel.html) will unfold a n-d array into vector. Similar to the `:` operator in Matlab > 2. [`ravel`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ravel.html) will unfold a n-d array into vector. Similar to the `:` operator in Matlab
> 3. [`reshape`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html) works similarly to reshape in Matlab, but be careful becase the default order is different from Matlab. > 3. [`reshape`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html) works similarly to reshape in Matlab, but be careful becase the default order is different from Matlab.
> 4. [`.T`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.T.html) does a transpose in `numpy` > 4. [`.T`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.T.html) does a transpose in `numpy`
> 5. The final element of an array is indexed with `-1` in `numpy`, as opposed to `end` in Matlab > 5. The final element of an array is indexed with `-1` in `numpy`, as opposed to `end` in Matlab
> 6. [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) and [`eigs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html?highlight=eigs#scipy.sparse.linalg.eigs) come from the [`scipy.sparse.linalg`](https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html) package > 6. [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) and [`eigs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html?highlight=eigs#scipy.sparse.linalg.eigs) come from the [`scipy.sparse.linalg`](https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html) package
> 7. [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) and [`eigs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html?highlight=eigs#scipy.sparse.linalg.eigs) are very similar to their Matlab counterparts, but be careful because Matlab `svds` returns $U$, $S$, and $V$, whereas python `svds` returns $U$, $S$, and $V^T$ > 7. [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) and [`eigs`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html?highlight=eigs#scipy.sparse.linalg.eigs) are very similar to their Matlab counterparts, but be careful because Matlab `svds` returns $U$, $S$, and $V$, whereas python `svds` returns $U$, $S$, and $V^T$
> 8. We index into the output of `eigs(W@W.T, dPCA_int)[1]` to only return the 2nd output (index 1) > 8. We index into the output of `eigs(W@W.T, dPCA_int)[1]` to only return the 2nd output (index 1)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# randomise the subject order # randomise the subject order
random.shuffle(in_list) random.shuffle(in_list)
# load and unravel brainmask # load and unravel brainmask
mask = in_mask.get_fdata().ravel() mask = in_mask.get_fdata(caching='unchanged') != 0.0
# function to demean the data # function to demean the data
def demean(x): def demean(x):
return x - np.mean(x, axis=0) return x - np.mean(x, axis=0)
# loop through input files/subjects # loop through input files/subjects
for i, f in enumerate(in_list): for i, f in enumerate(in_list):
# read data # read data
print(f'Reading data file {f.get_filename()}') print(f'Reading data file {f.get_filename()}')
grot = f.get_fdata() grot = f.get_fdata(caching='unchanged')[mask].T
grot = np.reshape(grot, [-1, grot.shape[-1]])
grot = grot[mask!=0, :].T
# demean # demean
print(f'\tRemoving mean image') print(f'\tRemoving mean image')
grot = demean(grot) grot = demean(grot)
# var-norm # var-norm
if sep_vn: if sep_vn:
print(f'\tNormalising by voxel-wise variance') print(f'\tNormalising by voxel-wise variance')
[uu, ss, vt] = svds(grot, k=30) [uu, ss, vt] = svds(grot, k=30)
vt[np.abs(vt) < (2.3 * np.std(vt))] = 0; vt[np.abs(vt) < (2.3 * np.std(vt))] = 0;
stddevs = np.maximum(np.std(grot - (uu @ np.diag(ss) @ vt), axis=0), 0.001) stddevs = np.maximum(np.std(grot - (uu @ np.diag(ss) @ vt), axis=0), 0.001)
grot = grot/stddevs grot = grot/stddevs
if i == 0: if i == 0:
W = demean(grot) W = demean(grot)
else: else:
# concat # concat
W = np.concatenate((W, demean(grot)), axis=0) W = np.concatenate((W, demean(grot)), axis=0)
# reduce W to dPCA_int eigenvectors # reduce W to dPCA_int eigenvectors
if W.shape[0]-10 > dPCA_int: if W.shape[0]-10 > dPCA_int:
print(f'\tReducing data matrix to a {dPCA_int} dimensional subspace') print(f'\tReducing data matrix to a {dPCA_int} dimensional subspace')
uu = eigs(W@W.T, dPCA_int)[1] uu = eigs(np.dot(W, W.T), dPCA_int)[1]
uu = np.real(uu) uu = np.real(uu)
W = uu.T @ W W = np.dot(uu.T, W)
# reshape and save # reshape and save
grot = np.zeros([mask.shape[0], dPCA_out]) grot = np.zeros([mask.size, dPCA_out])
grot[mask!=0, :] = W[:dPCA_out, :].T grot[mask.ravel(), :] = W[:dPCA_out, :].T
grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,)) grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))
print(f'Save to {GO}') print(f'Save to {GO}')
nb.Nifti1Image(grot, affine=in_list[0].affine).to_filename(GO) nb.Nifti1Image(grot, affine=in_list[0].affine).to_filename(GO)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<a class="anchor" id="run-python-melodic"></a> <a class="anchor" id="run-python-melodic"></a>
### Run ```melodic``` ### Run ```melodic```
Generate a command line string and run group ```melodic``` on the Matlab MIGP dimension reduced data with a dimension of 10 components. We disable MIGP because it was already run separately in Matlab. Generate a command line string and run group ```melodic``` on the Matlab MIGP dimension reduced data with a dimension of 10 components. We disable MIGP because it was already run separately in Matlab.
> **Note**: > **Note**:
> 1. Here we use python [f-strings](https://www.python.org/dev/peps/pep-0498/), formally known as literal string interpolation, which allow for easy formatting > 1. Here we use python [f-strings](https://www.python.org/dev/peps/pep-0498/), formally known as literal string interpolation, which allow for easy formatting
> 2. [`op.join`](https://docs.python.org/3.7/library/os.path.html#os.path.join) will join path strings using the platform-specific directory separator > 2. [`op.join`](https://docs.python.org/3.7/library/os.path.html#os.path.join) will join path strings using the platform-specific directory separator
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# generate melodic command line string # generate melodic command line string
melodic_cmd = f"melodic -i pyMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o pymigp.gica" melodic_cmd = f"melodic -i pyMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o pymigp.gica"
print(melodic_cmd) print(melodic_cmd)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
> **Note:** > **Note:**
> 1. Here we use the `!` operator to execute the command in the shell > 1. Here we use the `!` operator to execute the command in the shell
> 2. The `{}` will expand the contained python variable in the shell > 2. The `{}` will expand the contained python variable in the shell
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# run melodic # run melodic
! {melodic_cmd} ! {melodic_cmd}
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<a class="anchor" id="plot-python-group-ics"></a> <a class="anchor" id="plot-python-group-ics"></a>
### Plot group ICs ### Plot group ICs
Now we can load and plot the group ICs generated by ```melodic```. Now we can load and plot the group ICs generated by ```melodic```.
This function will be used to plot ICs: This function will be used to plot ICs:
> **NOTE:** > **NOTE:**
> 1. Here we use [`plot_stat_map`](https://nilearn.github.io/modules/generated/nilearn.plotting.plot_stat_map.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to plot the orthographic images > 1. Here we use [`plot_stat_map`](https://nilearn.github.io/modules/generated/nilearn.plotting.plot_stat_map.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to plot the orthographic images
> 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots > 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots
> 3. [`find_xyz_cut_coords`](https://nilearn.github.io/modules/generated/nilearn.plotting.find_xyz_cut_coords.html) from the [`nilearn`](https://nilearn.github.io/index.html) package will find the image coordinates of the center of the largest activation connected component > 3. [`find_xyz_cut_coords`](https://nilearn.github.io/modules/generated/nilearn.plotting.find_xyz_cut_coords.html) from the [`nilearn`](https://nilearn.github.io/index.html) package will find the image coordinates of the center of the largest activation connected component
> 4. [`zip`](https://docs.python.org/3.3/library/functions.html#zip) takes iterables and aggregates them in a tuple. Here it is used to iterate through two lists simultaneously > 4. [`zip`](https://docs.python.org/3.3/library/functions.html#zip) takes iterables and aggregates them in a tuple. Here it is used to iterate through two lists simultaneously
> 5. [`iter_img`](https://nilearn.github.io/modules/generated/nilearn.image.iter_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package creates an iterator from an image that steps through each volume/time-point of the image > 5. [`iter_img`](https://nilearn.github.io/modules/generated/nilearn.image.iter_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package creates an iterator from an image that steps through each volume/time-point of the image
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def map_plot(d): def map_plot(d):
N = d.shape[-1] N = d.shape[-1]
fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N)) fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
for img, ax0 in zip(image.iter_img(d), ax.ravel()): for img, ax0 in zip(image.iter_img(d), ax.ravel()):
coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5) coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)
plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0) plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)
return fig return fig
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Hopefully you can see some familiar looking RSN spatial patterns: Hopefully you can see some familiar looking RSN spatial patterns:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
ics = nb.load('pymigp.gica/melodic_IC.nii.gz') ics = nb.load('pymigp.gica/melodic_IC.nii.gz')
fig = map_plot(ics) fig = map_plot(ics)
``` ```
%% Cell type:code id: tags:
``` python
```
......
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