diff --git a/talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m b/talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m
new file mode 100644
index 0000000000000000000000000000000000000000..81fdfe00fa0a86a1bb2e3549c3c3de20217638af
--- /dev/null
+++ b/talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m
@@ -0,0 +1,63 @@
+
+
+% 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
+
+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
+
+%%%%%%%%%%%%%%%%%%%%%%%%%% END OF USER OPTIONS
+
+Nsub=length(INlist); [~,r]=sort(rand(Nsub,1)); % will process subjects in random order
+mask=read_avw(INmask);
+masksize=size(mask);
+mask=reshape(mask,prod(masksize),1);
+
+demean = @(x) x - repmat(mean(x,1),[size(x,1), 1, 1]);
+ss_svds = @(x,n) svds(x, n);
+
+for i = 1:Nsub
+
+ % read data
+ filename=[INlist(r(i)).folder, '/', INlist(r(i)).name];
+ grot=double(read_avw(filename));
+ grot=reshape(grot,prod(masksize),size(grot,4));
+ 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);
+
+
+ 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
+ [uu,dd]=eigs(W*W',dPCAint);
+ W=uu'*W;
+ clear uu;
+ end
+
+ end
+
+end
+
+grot=zeros(prod(masksize),dPCAout);
+grot(mask~=0,:)=W(1:dPCAout,:)'; grot=reshape(grot,[masksize ,dPCAout]);
+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/matlab_MIGP.ipynb b/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..830c173746af46ee6a603d0e700c26dbd6c7e19e
--- /dev/null
+++ b/talks/matlab_vs_python/migp/matlab_MIGP.ipynb
@@ -0,0 +1,80 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from nilearn import plotting\n",
+ "from nilearn import image\n",
+ "import nibabel as nb\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%bash\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"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ics = nb.load('matMIGP_dim20.ica/melodic_IC.nii.gz')\n",
+ "\n",
+ "N = ics.shape[-1]\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)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/talks/matlab_vs_python/migp/python_MIGP.ipynb b/talks/matlab_vs_python/migp/python_MIGP.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..f029476f0aacb22f53922c473df75a78788d27df
--- /dev/null
+++ b/talks/matlab_vs_python/migp/python_MIGP.ipynb
@@ -0,0 +1,168 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import glob\n",
+ "import random\n",
+ "import nibabel as nb\n",
+ "import numpy as np\n",
+ "from scipy.sparse.linalg import svds, eigs\n",
+ "import matplotlib.pyplot as plt\n",
+ "from nilearn import plotting\n",
+ "from nilearn import image"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "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"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "random.shuffle(in_list)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mask = in_mask.get_fdata().ravel()\n",
+ "\n",
+ "def demean(x):\n",
+ " return x - np.mean(x, axis=0)\n",
+ "\n",
+ "for i, f in enumerate(in_list):\n",
+ " \n",
+ " print(i, end=' ')\n",
+ " \n",
+ " # read data\n",
+ " grot = f.get_fdata()\n",
+ " grot = np.reshape(grot, [-1, grot.shape[-1]])\n",
+ " grot = grot[mask!=0, :].T\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",
+ " \n",
+ " if i == 0:\n",
+ " W = demean(grot)\n",
+ " else:\n",
+ " # concat\n",
+ " W = np.concatenate((W, demean(grot)), axis=0)\n",
+ " \n",
+ " # reduce W to dPCA_int eigenvectors\n",
+ " if W.shape[0]-10 > dPCA_int:\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"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "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,))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nb.Nifti1Image(grot, affine=in_list[0].affine, header=in_list[0].header).to_filename('pyMIGP.nii.gz')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "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"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ics = nb.load('pymigp_dim20.ica/melodic_IC.nii.gz')\n",
+ "\n",
+ "N = ics.shape[-1]\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)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}