diff --git a/talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m b/talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m
index 81fdfe00fa0a86a1bb2e3549c3c3de20217638af..e14287e6f6aa5aaf5adeeadc5e8396db015cbd08 100644
--- a/talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m
+++ b/talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m
@@ -1,19 +1,18 @@
-
-
-% MIGP - MELODIC's Incremental Group-PCA
+%% MIGP - MELODIC's Incremental Group-PCA
% Steve Smith, FMRIB, 2012-2013
% not yet finally tweaked/evaluated, or published - please do not pass on.
% adapted by Sean Fitz, 2020
-%%%%%%%%%%%%%%%%%%%%%%%%%% USER OPTIONS
+%% User Options (equivalent to melodic defaults)
-INlist=dir('data/cobre/fmri_*.nii.gz'); % list of input 4D NIFTI standard space datasets
-INmask='data/brain_mask.nii.gz'; % 3D NIFTI mask image
-GO='matMIGP'; % output filename
-dPCAint=550; % internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)
-dPCAout=100; % number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality
+INlist=dir('~/nilearn_data/cobre/fmri_*_smooth.nii.gz'); % list of input 4D NIFTI standard space datasets
+INmask='~/nilearn_data/brain_mask.nii.gz'; % 3D NIFTI mask image
+GO='matMIGP'; % output filename
+dPCAint=299; % internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)
+dPCAout=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
-%%%%%%%%%%%%%%%%%%%%%%%%%% END OF USER OPTIONS
+%% Run MIGP
Nsub=length(INlist); [~,r]=sort(rand(Nsub,1)); % will process subjects in random order
mask=read_avw(INmask);
@@ -27,26 +26,33 @@ for i = 1:Nsub
% read data
filename=[INlist(r(i)).folder, '/', INlist(r(i)).name];
+ fprintf('Reading data file %s\n', filename);
grot=double(read_avw(filename));
grot=reshape(grot,prod(masksize),size(grot,4));
+
+ % demean
+ fprintf('\tRemoving mean image\n');
grot=demean(grot(mask~=0,:)');
- % var-norm
- [uu,ss,vv]=ss_svds(grot,30);
- vv(abs(vv)<2.3*std(vv(:)))=0;
- stddevs=max(std(grot-uu*ss*vv'),0.001);
- grot=grot./repmat(stddevs,size(grot,1),1);
+ % var-norm separately for each subject
+ if sep_vn==true
+ fprintf('\tNormalising by voxel-wise variance\r');
+ [uu,ss,vv]=ss_svds(grot,30);
+ vv(abs(vv)<2.3*std(vv(:)))=0;
+ stddevs=max(std(grot-uu*ss*vv'),0.001);
+ grot=grot./repmat(stddevs,size(grot,1),1);
+ end
if (i==1)
W=demean(grot); clear grot;
else
-
% concat
W=[W; demean(grot)]; clear grot;
% reduce W to dPCAint eigenvectors
if size(W,1)-10 > dPCAint
+ fprintf('\tReducing data matrix to a %i dimensional subspace\n', dPCAint);
[uu,dd]=eigs(W*W',dPCAint);
W=uu'*W;
clear uu;
@@ -56,8 +62,11 @@ for i = 1:Nsub
end
+%% Reshape and save
+
grot=zeros(prod(masksize),dPCAout);
grot(mask~=0,:)=W(1:dPCAout,:)'; grot=reshape(grot,[masksize ,dPCAout]);
+fprintf('Save to %s.nii.gz\n', GO)
save_avw(grot,GO,'f',[1 1 1 1]);
system(sprintf('fslcpgeom %s %s -d',filename,GO));
diff --git a/talks/matlab_vs_python/migp/fetch_data.ipynb b/talks/matlab_vs_python/migp/fetch_data.ipynb
index 33f4f27bb97a47483b9f2f0308b95c6c460369e2..44ae65f3ebbc9bafd5da9239817550537d1dc1a0 100644
--- a/talks/matlab_vs_python/migp/fetch_data.ipynb
+++ b/talks/matlab_vs_python/migp/fetch_data.ipynb
@@ -16,7 +16,7 @@
"* [Clean the data](#clean-the-data)\n",
"* [Run `melodic`](#run-melodic)\n",
"\n",
- "Firstly we will import the necessary packages for this notebook: "
+ "Firstly we will import the necessary packages for this notebook: "
]
},
{
@@ -69,7 +69,9 @@
"<a class=\"anchor\" id=\"download-the-data\"></a>\n",
"## Download the data\n",
"\n",
- "We use a method from [`nilearn`](https://nilearn.github.io/index.html) called `fetch_cobre` to download the fMRI data:"
+ "Create a directory in the users home directory to store the downloaded data:\n",
+ "\n",
+ "`expanduser` will expand the `~` to the be users home directory:"
]
},
{
@@ -78,13 +80,25 @@
"metadata": {},
"outputs": [],
"source": [
- "# Create a directory in the users home directory to store the downloaded data\n",
"data_dir = op.expanduser('~/nilearn_data')\n",
"\n",
"if not op.exists(data_dir):\n",
- " os.makedirs(data_dir)\n",
- " \n",
- "# Download the data (if not already downloaded)\n",
+ " os.makedirs(data_dir)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Download the data (if not already downloaded). We use a method from [`nilearn`](https://nilearn.github.io/index.html) called `fetch_cobre` to download the fMRI data:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
"d = datasets.fetch_cobre(data_dir=data_dir) "
]
},
@@ -158,7 +172,7 @@
"outputs": [],
"source": [
"# generate melodic command line string\n",
- "melodic_cmd = f\"melodic -i {','.join(smooth)} --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --Oall -o cobre.gica \"\n",
+ "melodic_cmd = f\"melodic -i {','.join(smooth)} --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v -o cobre.gica \"\n",
"print(melodic_cmd)\n",
"\n",
"# run melodic\n",
@@ -183,6 +197,13 @@
"ics = nb.load('cobre.gica/melodic_IC.nii.gz')\n",
"fig = map_plot(ics)"
]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
diff --git a/talks/matlab_vs_python/migp/matlab_MIGP.ipynb b/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
index 830c173746af46ee6a603d0e700c26dbd6c7e19e..779e4648870b3fba904c68beb0caa5836fd4c142 100644
--- a/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
+++ b/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
@@ -1,5 +1,19 @@
{
"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Matlab MIGP\n",
+ "\n",
+ "This notebook will load the dimension reduced data from Matlab MIGP, run group ICA, and then plot the group ICs.\n",
+ "\n",
+ "* [Run `melodic`](#run-melodic)\n",
+ "* [Plot group ICs](#plot-group-ics)\n",
+ "\n",
+ "Firstly we will import the necessary packages for this notebook: "
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -10,7 +24,17 @@
"from nilearn import image\n",
"import nibabel as nb\n",
"import matplotlib.pyplot as plt\n",
- "import numpy as np"
+ "import numpy as np\n",
+ "import os.path as op"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "It will be necessary to know the location where the data was stored so that we can load the brainmask. \n",
+ "\n",
+ "`expanduser` will expand the `~` to the be users home directory:"
]
},
{
@@ -19,16 +43,28 @@
"metadata": {},
"outputs": [],
"source": [
- "%%bash\n",
+ "data_dir = op.expanduser('~/nilearn_data')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"run-melodic\"></a>\n",
+ "### Run ```melodic```\n",
"\n",
- "melodic -i matMIGP.nii.gz \\\n",
- "\t--mask=data/brain_mask.nii.gz \\\n",
- "\t-d 20 \\\n",
- "\t-v \\\n",
- "\t--nobet \\\n",
- "\t--disableMigp \\\n",
- "\t--varnorm \\\n",
- "\t-o matMIGP_dim20.ica"
+ "Generate a command line string and run group ```melodic``` on the Matlab MIGP dimension reduced data with a dimension of 10 components:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# generate melodic command line string\n",
+ "melodic_cmd = f\"melodic -i matMIGP.nii.gz --mask={op.join(data_dir, 'brain_mask.nii.gz')} -d 10 -v --nobet --disableMigp -o matmigp.gica\"\n",
+ "print(melodic_cmd)"
]
},
{
@@ -37,15 +73,56 @@
"metadata": {},
"outputs": [],
"source": [
- "ics = nb.load('matMIGP_dim20.ica/melodic_IC.nii.gz')\n",
+ "# run melodic\n",
+ "! {melodic_cmd}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"plot-group-ics\"></a>\n",
+ "### Plot group ICs\n",
+ "\n",
+ "Now we can load and plot the group ICs generated by ```melodic```.\n",
+ "\n",
+ "This function will be used to plot ICs:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def map_plot(d):\n",
"\n",
- "N = ics.shape[-1]\n",
+ " N = d.shape[-1]\n",
"\n",
- "fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))\n",
+ " fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))\n",
"\n",
- "for img, ax0 in zip(image.iter_img(ics), ax.ravel()):\n",
- " coord = plotting.find_xyz_cut_coords(img, activation_threshold=2.3)\n",
- " plotting.plot_stat_map(img, cut_coords=coord, vmax=5, axes=ax0)"
+ " for img, ax0 in zip(image.iter_img(d), ax.ravel()):\n",
+ " coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)\n",
+ " plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)\n",
+ " \n",
+ " return fig"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Hopefully you can see some familiar looking RSN spatial patterns:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ics = nb.load('matmigp.gica/melodic_IC.nii.gz')\n",
+ "fig = map_plot(ics)"
]
},
{
diff --git a/talks/matlab_vs_python/migp/python_MIGP.ipynb b/talks/matlab_vs_python/migp/python_MIGP.ipynb
index f029476f0aacb22f53922c473df75a78788d27df..739457281bde1ff8dafa88cf0cf46d9d32eb1212 100644
--- a/talks/matlab_vs_python/migp/python_MIGP.ipynb
+++ b/talks/matlab_vs_python/migp/python_MIGP.ipynb
@@ -1,5 +1,20 @@
{
"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Python MIGP\n",
+ "\n",
+ "This notebook will perform *python* MIGP dimension reduction, run group ICA, and then plot the group ICs.\n",
+ "\n",
+ "* [Run python `MIGP`](#run-python-migp)\n",
+ "* [Run `melodic`](#run-melodic)\n",
+ "* [Plot group ICs](#plot-group-ics)\n",
+ "\n",
+ "Firstly we will import the necessary packages for this notebook: "
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -13,7 +28,17 @@
"from scipy.sparse.linalg import svds, eigs\n",
"import matplotlib.pyplot as plt\n",
"from nilearn import plotting\n",
- "from nilearn import image"
+ "from nilearn import image\n",
+ "import os.path as op"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "It will be necessary to know the location where the data was stored so that we can load the brainmask. \n",
+ "\n",
+ "`expanduser` will expand the `~` to the be users home directory:"
]
},
{
@@ -22,11 +47,17 @@
"metadata": {},
"outputs": [],
"source": [
- "in_list = [nb.load(f) for f in glob.glob('data/cobre/fmri_*.nii.gz')]\n",
- "in_mask = nb.load('data/brain_mask.nii.gz')\n",
- "go = '3DMIGP'\n",
- "dPCA_int = 550\n",
- "dPCA_out = 100"
+ "data_dir = op.expanduser('~/nilearn_data')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"run-python-migp\"></a>\n",
+ "### Run python `MIGP`\n",
+ "\n",
+ "Firstly we need to set the MIGP parameters:"
]
},
{
@@ -35,7 +66,15 @@
"metadata": {},
"outputs": [],
"source": [
- "random.shuffle(in_list)"
+ "# create lists of (nibabel) image objects\n",
+ "in_list = [nb.load(f) for f in glob.glob(f'{data_dir}/cobre/fmri_*_smooth.nii.gz')]\n",
+ "in_mask = nb.load(f'{data_dir}/brain_mask.nii.gz')\n",
+ "\n",
+ "# set user parameters (equivalent to melodic defaults)\n",
+ "GO = 'pyMIGP.nii.gz' # output filename\n",
+ "dPCA_int = 299 # internal number of components - typically 2-4 times number of timepoints in each run (if you have enough RAM for that)\n",
+ "dPCA_out = 299 # number of eigenvectors to output - should be less than dPCAint and more than the final ICA dimensionality\n",
+ "sep_vn = False # switch on separate variance nomalisation for each input dataset"
]
},
{
@@ -44,26 +83,36 @@
"metadata": {},
"outputs": [],
"source": [
+ "# randomise the subject order\n",
+ "random.shuffle(in_list)\n",
+ "\n",
+ "# load and unravel brainmask\n",
"mask = in_mask.get_fdata().ravel()\n",
"\n",
+ "# function to demean the data\n",
"def demean(x):\n",
" return x - np.mean(x, axis=0)\n",
"\n",
+ "# loop through input files/subjects\n",
"for i, f in enumerate(in_list):\n",
" \n",
- " print(i, end=' ')\n",
- " \n",
" # read data\n",
+ " print(f'Reading data file {f.get_filename()}')\n",
" grot = f.get_fdata()\n",
" grot = np.reshape(grot, [-1, grot.shape[-1]])\n",
" grot = grot[mask!=0, :].T\n",
+ " \n",
+ " # demean\n",
+ " print(f'\\tRemoving mean image')\n",
" grot = demean(grot)\n",
" \n",
" # var-norm\n",
- " [uu, ss, vt] = svds(grot, k=30)\n",
- " vt[np.abs(vt) < (2.3 * np.std(vt))] = 0;\n",
- " stddevs = np.maximum(np.std(grot - (uu @ np.diag(ss) @ vt), axis=0), 0.001)\n",
- " grot = grot/stddevs\n",
+ " if sep_vn:\n",
+ " print(f'\\tNormalising by voxel-wise variance')\n",
+ " [uu, ss, vt] = svds(grot, k=30)\n",
+ " vt[np.abs(vt) < (2.3 * np.std(vt))] = 0;\n",
+ " stddevs = np.maximum(np.std(grot - (uu @ np.diag(ss) @ vt), axis=0), 0.001)\n",
+ " grot = grot/stddevs\n",
" \n",
" if i == 0:\n",
" W = demean(grot)\n",
@@ -73,23 +122,28 @@
" \n",
" # reduce W to dPCA_int eigenvectors\n",
" if W.shape[0]-10 > dPCA_int:\n",
+ " print(f'\\tReducing data matrix to a {dPCA_int} dimensional subspace')\n",
" uu = eigs(W@W.T, dPCA_int)[1]\n",
" uu = np.real(uu)\n",
" W = uu.T @ W\n",
" \n",
- " f.uncache()\n",
- " \n"
+ "# reshape and save\n",
+ "grot = np.zeros([mask.shape[0], dPCA_out])\n",
+ "grot[mask!=0, :] = W[:dPCA_out, :].T\n",
+ "grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))\n",
+ "\n",
+ "print(f'Save to {GO}')\n",
+ "nb.Nifti1Image(grot, affine=in_list[0].affine).to_filename(GO)\n"
]
},
{
- "cell_type": "code",
- "execution_count": null,
+ "cell_type": "markdown",
"metadata": {},
- "outputs": [],
"source": [
- "grot = np.zeros([mask.shape[0], dPCA_out])\n",
- "grot[mask!=0, :] = W[:dPCA_out, :].T\n",
- "grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))"
+ "<a class=\"anchor\" id=\"run-melodic\"></a>\n",
+ "### Run ```melodic```\n",
+ "\n",
+ "Generate a command line string and run group ```melodic``` on the Python MIGP dimension reduced data with a dimension of 10 components:"
]
},
{
@@ -98,7 +152,9 @@
"metadata": {},
"outputs": [],
"source": [
- "nb.Nifti1Image(grot, affine=in_list[0].affine, header=in_list[0].header).to_filename('pyMIGP.nii.gz')"
+ "# generate melodic command line string\n",
+ "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\"\n",
+ "print(melodic_cmd)"
]
},
{
@@ -107,16 +163,20 @@
"metadata": {},
"outputs": [],
"source": [
- "%%bash\n",
- "\n",
- "melodic -i pyMIGP.nii.gz \\\n",
- "\t--mask=data/brain_mask.nii.gz \\\n",
- "\t-d 20 \\\n",
- "\t-v \\\n",
- "\t--nobet \\\n",
- "\t--disableMigp \\\n",
- "\t--varnorm \\\n",
- "\t-o pymigp_dim20.ica"
+ "# run melodic\n",
+ "! {melodic_cmd}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<a class=\"anchor\" id=\"plot-group-ics\"></a>\n",
+ "### Plot group ICs\n",
+ "\n",
+ "Now we can load and plot the group ICs generated by ```melodic```.\n",
+ "\n",
+ "This function will be used to plot ICs:"
]
},
{
@@ -125,15 +185,34 @@
"metadata": {},
"outputs": [],
"source": [
- "ics = nb.load('pymigp_dim20.ica/melodic_IC.nii.gz')\n",
+ "def map_plot(d):\n",
"\n",
- "N = ics.shape[-1]\n",
+ " N = d.shape[-1]\n",
"\n",
- "fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))\n",
+ " fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))\n",
"\n",
- "for img, ax0 in zip(image.iter_img(ics), ax.ravel()):\n",
- " coord = plotting.find_xyz_cut_coords(img, activation_threshold=1.5)\n",
- " plotting.plot_stat_map(img, cut_coords=coord, vmax=5, axes=ax0)"
+ " for img, ax0 in zip(image.iter_img(d), ax.ravel()):\n",
+ " coord = plotting.find_xyz_cut_coords(img, activation_threshold=3.5)\n",
+ " plotting.plot_stat_map(img, cut_coords=coord, vmax=10, axes=ax0)\n",
+ " \n",
+ " return fig"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Hopefully you can see some familiar looking RSN spatial patterns:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ics = nb.load('pymigp.gica/melodic_IC.nii.gz')\n",
+ "fig = map_plot(ics)"
]
},
{