diff --git a/talks/matlab_vs_python/migp/fetch_data.ipynb b/talks/matlab_vs_python/migp/MIGP.ipynb
similarity index 56%
rename from talks/matlab_vs_python/migp/fetch_data.ipynb
rename to talks/matlab_vs_python/migp/MIGP.ipynb
index 6778e404ee079e8b4f132865f8ff52a2b7307a33..1a0baf9f839a6b10ddc443deb4d46048084d1e6d 100644
--- a/talks/matlab_vs_python/migp/fetch_data.ipynb
+++ b/talks/matlab_vs_python/migp/MIGP.ipynb
@@ -4,15 +4,29 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Fetch Data\n",
+ "# MIGP\n",
"\n",
- "This notebook will download an open fMRI dataset (~50MB) for use in the MIGP demo. It also regresses confounds from the data and performs spatial smoothing with 10mm FWHM.\n",
+ "For group ICA, `melodic` uses multi-session temporal concatenation. This will perform a single 2D ICA run on the concatenated data matrix (obtained by stacking all 2D data matrices of every single data set on top of each other).\n",
"\n",
- "This data is a derivative from the COBRE sample found in the International Neuroimaging Data-sharing Initiative (http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html), originally released under Creative Commons - Attribution Non-Commercial.\n",
+ "\n",
"\n",
- "It comprises 10 preprocessed resting-state fMRI selected from 72 patients diagnosed with schizophrenia and 74 healthy controls (6mm isotropic, TR=2s, 150 volumes).\n",
+ "Resulting in **high dimension** datasets!\n",
+ "\n",
+ "Furthermore, with ICA we are typically only interested in a comparitively low dimension decomposition so that we can capture spatially extended networks.\n",
+ "\n",
+ "Therefore the first step is to reduce the dimensionality of the data. This can be achieved in a number of ways, but `melodic`, by default, uses `MIGP`.\n",
+ "\n",
+ "> MIGP is an incremental approach that aims to provide a very close approximation to full temporal concatenation followed by PCA, but without the large memory requirements *(Smith et al., 2014)*.\n",
+ "\n",
+ "Essentially, MIGP stacks the datasets incrementally in the temporal dimension, and whenever the temporal dimension exceeds a specified size, a PCA-based temporal reduction is performed.\n",
"\n",
- "* [Download the data](#download-the-data)\n",
+ "> MIGP does not increase at all in memory requirement with increasing numbers of subjects, no large matrices are ever formed, and the computation time scales linearly with the number of subjects. It is easily parallelisable, simply by applying the approach in parallel to subsets of subjects, and then combining across these with the same “concatenate and reduce†approach described above *(Smith et al., 2014)*.\n",
+ "\n",
+ "## This notebook\n",
+ "\n",
+ "This notebook will download an open fMRI dataset (~50MB) for use in the MIGP demo, regresses confounds from the data, performs spatial smoothing with 10mm FWHM, and then runs group `melodic` with `MIGP`.\n",
+ "\n",
+ "* [Fetch the data](#download-the-data)\n",
"* [Clean the data](#clean-the-data)\n",
"* [Run `melodic`](#run-melodic)\n",
"* [Plot group ICs](#plot-group-ics)\n",
@@ -42,11 +56,15 @@
"metadata": {},
"source": [
"<a class=\"anchor\" id=\"download-the-data\"></a>\n",
- "## Download the data\n",
+ "## Fetch the data\n",
+ "\n",
+ "This data is a derivative from the [COBRE](http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html) sample found in the International Neuroimaging Data-sharing Initiative, originally released under Creative Commons - Attribution Non-Commercial.\n",
+ "\n",
+ "It comprises 10 preprocessed resting-state fMRI selected from 72 patients diagnosed with schizophrenia and 74 healthy controls (6mm isotropic, TR=2s, 150 volumes).\n",
"\n",
"Create a directory in the users home directory to store the downloaded data:\n",
"\n",
- "> **NOTE:** `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:"
]
},
{
@@ -67,7 +85,7 @@
"source": [
"Download the data (if not already downloaded):\n",
"\n",
- "> **Note:** We use a method from [`nilearn`](https://nilearn.github.io/index.html) called `fetch_cobre` to download the fMRI data"
+ "> **Note:** We use a method from [`nilearn`](https://nilearn.github.io/index.html) called [`fetch_cobre`](https://nilearn.github.io/modules/generated/nilearn.datasets.fetch_cobre.html) to download the fMRI data"
]
},
{
@@ -89,9 +107,9 @@
"Regress confounds from the data and to spatially smooth the data with a gaussian filter of 10mm FWHM.\n",
"\n",
"> **Note:**\n",
- "> 1. We use `clean_img` from the [`nilearn`](https://nilearn.github.io/index.html) package to regress confounds from the data\n",
- "> 2. We use `smooth_img` from the [`nilearn`](https://nilearn.github.io/index.html) package to spatially smooth the data\n",
- "> 3. `zip` takes iterables and aggregates them in a tuple. Here it is used to iterate through four lists simultaneously\n",
+ "> 1. We use [`clean_img`](https://nilearn.github.io/modules/generated/nilearn.image.clean_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to regress confounds from the data\n",
+ "> 2. We use [`smooth_img`](https://nilearn.github.io/modules/generated/nilearn.image.smooth_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to spatially smooth the data\n",
+ "> 3. [`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 four lists simultaneously\n",
"> 4. We use list comprehension to loop through all the filenames and append suffixes\n"
]
},
@@ -120,10 +138,10 @@
"To run ```melodic``` we will need a brain mask in MNI152 space at the same resolution as the fMRI. \n",
"\n",
"> **Note:**\n",
- "> 1. We use `load_mni152_brain_mask` from the [`nilearn`](https://nilearn.github.io/index.html) package to load the MNI152 mask\n",
- "> 2. We use `resample_to_img` from the [`nilearn`](https://nilearn.github.io/index.html) package to resample the mask to the resolution of the fMRI \n",
- "> 3. We use `math_img` from the [`nilearn`](https://nilearn.github.io/index.html) package to binarize the resample mask\n",
- "> 4. The mask is plotted using `plot_anat` from the [`nilearn`](https://nilearn.github.io/index.html) package"
+ "> 1. We use [`load_mni152_brain_mask`](https://nilearn.github.io/modules/generated/nilearn.datasets.load_mni152_brain_mask.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to load the MNI152 mask\n",
+ "> 2. We use [`resample_to_img`](https://nilearn.github.io/modules/generated/nilearn.image.resample_to_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to resample the mask to the resolution of the fMRI \n",
+ "> 3. We use [`math_img`](https://nilearn.github.io/modules/generated/nilearn.image.math_img.html) from the [`nilearn`](https://nilearn.github.io/index.html) package to binarize the resample mask\n",
+ "> 4. The mask is plotted using [`plot_anat`](https://nilearn.github.io/modules/generated/nilearn.plotting.plot_anat.html) from the [`nilearn`](https://nilearn.github.io/index.html) package"
]
},
{
@@ -157,8 +175,8 @@
"\n",
"> **Note**: \n",
"> 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\n",
- "> 2. `op.join` will join path strings using the platform-specific directory separator\n",
- "> 3. `','.join(smooth)` will create a comma seprated string of all the items in the list `smooth`"
+ "> 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\n",
+ "> 3. [`','.join(smooth)`](https://docs.python.org/3/library/stdtypes.html#str.join) will create a comma seprated string of all the items in the list `smooth`"
]
},
{
@@ -203,11 +221,11 @@
"This function will be used to plot ICs:\n",
"\n",
"> **NOTE:**\n",
- "> 1. Here we use `plot_stat_map` from the `nilearn` package to plot the orthographic images\n",
- "> 2. `subplots` from `matplotlib.pyplot` creates a figure and multiple subplots\n",
- "> 3. `find_xyz_cut_coords` from the `nilearn` package will find the image coordinates of the center of the largest activation connected component\n",
- "> 4. `zip` takes iterables and aggregates them in a tuple. Here it is used to iterate through two lists simultaneously\n",
- "> 5. `iter_img` from the `nilearn` package creates an iterator from an image that steps through each volume/time-point of the image"
+ "> 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\n",
+ "> 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots\n",
+ "> 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\n",
+ "> 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\n",
+ "> 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"
]
},
{
@@ -248,13 +266,6 @@
"# plot\n",
"fig = map_plot(ics)"
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {
@@ -273,7 +284,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
+ "version": "3.7.6"
}
},
"nbformat": 4,
diff --git a/talks/matlab_vs_python/migp/concat_diag.png b/talks/matlab_vs_python/migp/concat_diag.png
new file mode 100644
index 0000000000000000000000000000000000000000..252b34034c426ffa797762ae538f090312dfb208
Binary files /dev/null and b/talks/matlab_vs_python/migp/concat_diag.png differ
diff --git a/talks/matlab_vs_python/migp/matlab_MIGP.ipynb b/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
index ca0089a2ca37a9db15c420cebb6b9196bce93963..b017170e77e4e04a6d8d2e69c930323f9c18566b 100644
--- a/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
+++ b/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
@@ -34,7 +34,7 @@
"source": [
"It will be necessary to know the location where the data was stored so that we can load the brainmask:\n",
"\n",
- "> **Note**: `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"
]
},
{
@@ -57,7 +57,7 @@
"\n",
"> **Note**: \n",
"> 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\n",
- "> 2. `op.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"
]
},
{
@@ -102,11 +102,11 @@
"This function will be used to plot ICs:\n",
"\n",
"> **NOTE:**\n",
- "> 1. Here we use `plot_stat_map` from the `nilearn` package to plot the orthographic images\n",
- "> 2. `subplots` from `matplotlib.pyplot` creates a figure and multiple subplots\n",
- "> 3. `find_xyz_cut_coords` from the `nilearn` package will find the image coordinates of the center of the largest activation connected component\n",
- "> 4. `zip` takes iterables and aggregates them in a tuple. Here it is used to iterate through two lists simultaneously\n",
- "> 5. `iter_img` from the `nilearn` package creates an iterator from an image that steps through each volume/time-point of the image"
+ "> 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\n",
+ "> 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots\n",
+ "> 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\n",
+ "> 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\n",
+ "> 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"
]
},
{
@@ -144,13 +144,6 @@
"ics = nb.load('matmigp.gica/melodic_IC.nii.gz')\n",
"fig = map_plot(ics)"
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {
@@ -169,7 +162,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
+ "version": "3.7.6"
}
},
"nbformat": 4,
diff --git a/talks/matlab_vs_python/migp/python_MIGP.ipynb b/talks/matlab_vs_python/migp/python_MIGP.ipynb
index 7c3188f27adce0a5dd8dd182713242dbd5494e57..c447f9162a0f30400cbdfef98070c14940de13ac 100644
--- a/talks/matlab_vs_python/migp/python_MIGP.ipynb
+++ b/talks/matlab_vs_python/migp/python_MIGP.ipynb
@@ -38,7 +38,7 @@
"source": [
"It will be necessary to know the location where the data was stored so that we can load the brainmask. \n",
"\n",
- "> **Note:** `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"
]
},
{
@@ -60,9 +60,9 @@
"Firstly we need to set the MIGP parameters:\n",
"\n",
"> **Note:**\n",
- "> 1. `glob.glob` will create a list of filenames that match the glob/wildcard pattern\n",
- "> 2. `nb.load` from the `nibabel` package will load the image into `nibabel.Nifti1Image` object. This will not load the actual data though.\n",
- "> 3. We use a list comprehension to loop through all the filenames and load them with `nibabel`"
+ "> 1. [`glob.glob`](https://docs.python.org/3/library/glob.html) will create a list of filenames that match the glob/wildcard pattern\n",
+ "> 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.\n",
+ "> 3. We use a list comprehension to loop through all the filenames and load them with [`nibabel`](https://nipy.org/nibabel/index.html)"
]
},
{
@@ -87,13 +87,13 @@
"metadata": {},
"source": [
"> **Note:**\n",
- "> 1. `random.shuffle` will shuffle a list, in this instance it shuffles the list of `nibabel.Nifti1Image` objects\n",
- "> 2. `ravel` will unfold a n-d array into vector. Similar to the `:` operator in Matlab\n",
- "> 3. `reshape` works similarly to reshape in Matlab, but be careful becase the default order is different from Matlab.\n",
- "> 4. `.T` does a transpose in `numpy`\n",
+ "> 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\n",
+ "> 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\n",
+ "> 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.\n",
+ "> 4. [`.T`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.T.html) does a transpose in `numpy`\n",
"> 5. The final element of an array is indexed with `-1` in `numpy`, as opposed to `end` in Matlab\n",
- "> 6. `svds` and `eigs` come from the `scipy.sparse.linalg` package\n",
- "> 7. `svds` and `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$\n",
+ "> 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\n",
+ "> 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$\n",
"> 8. We index into the output of `eigs(W@W.T, dPCA_int)[1]` to only return the 2nd output (index 1)"
]
},
@@ -107,7 +107,7 @@
"random.shuffle(in_list)\n",
"\n",
"# load and unravel brainmask\n",
- "mask = in_mask.get_fdata().ravel()\n",
+ "mask = in_mask.get_fdata(caching='unchanged') != 0.0\n",
"\n",
"# function to demean the data\n",
"def demean(x):\n",
@@ -118,9 +118,7 @@
" \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",
+ " grot = f.get_fdata(caching='unchanged')[mask].T\n",
" \n",
" # demean\n",
" print(f'\\tRemoving mean image')\n",
@@ -143,13 +141,13 @@
" # 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 = eigs(np.dot(W, W.T), dPCA_int)[1]\n",
" uu = np.real(uu)\n",
- " W = uu.T @ W\n",
+ " W = np.dot(uu.T, W)\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.zeros([mask.size, dPCA_out])\n",
+ "grot[mask.ravel(), :] = W[:dPCA_out, :].T\n",
"grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))\n",
"\n",
"print(f'Save to {GO}')\n",
@@ -163,7 +161,11 @@
"<a class=\"anchor\" id=\"run-python-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:"
+ "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.\n",
+ "\n",
+ "> **Note**: \n",
+ "> 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\n",
+ "> 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"
]
},
{
@@ -177,6 +179,15 @@
"print(melodic_cmd)"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> **Note:** \n",
+ "> 1. Here we use the `!` operator to execute the command in the shell\n",
+ "> 2. The `{}` will expand the contained python variable in the shell"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -196,7 +207,14 @@
"\n",
"Now we can load and plot the group ICs generated by ```melodic```.\n",
"\n",
- "This function will be used to plot ICs:"
+ "This function will be used to plot ICs:\n",
+ "\n",
+ "> **NOTE:**\n",
+ "> 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\n",
+ "> 2. [`subplots`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html) from `matplotlib.pyplot` creates a figure and multiple subplots\n",
+ "> 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\n",
+ "> 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\n",
+ "> 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"
]
},
{
@@ -259,7 +277,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
+ "version": "3.7.6"
}
},
"nbformat": 4,