added migp demo

talks/matlab_vs_python/migp/demo_MIGP_NIFTI.m

+
+
+% 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));
+

talks/matlab_vs_python/migp/matlab_MIGP.ipynb

+{
+ "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
+}
+%% Cell type:code id: tags:
+
+``` python
+from nilearn import plotting
+from nilearn import image
+import nibabel as nb
+import matplotlib.pyplot as plt
+import numpy as np
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%%bash
+
+melodic -i matMIGP.nii.gz \
+	--mask=data/brain_mask.nii.gz \
+	-d 20 \
+	-v \
+	--nobet \
+	--disableMigp \
+	--varnorm \
+	-o matMIGP_dim20.ica
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ics = nb.load('matMIGP_dim20.ica/melodic_IC.nii.gz')
+
+N = ics.shape[-1]
+
+fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
+
+for img, ax0 in zip(image.iter_img(ics), ax.ravel()):
+    coord = plotting.find_xyz_cut_coords(img, activation_threshold=2.3)
+    plotting.plot_stat_map(img, cut_coords=coord, vmax=5, axes=ax0)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```

talks/matlab_vs_python/migp/python_MIGP.ipynb

+{
+ "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
+}
+%% Cell type:code id: tags:
+
+``` python
+import glob
+import random
+import nibabel as nb
+import numpy as np
+from scipy.sparse.linalg import svds, eigs
+import matplotlib.pyplot as plt
+from nilearn import plotting
+from nilearn import image
+```
+
+%% Cell type:code id: tags:
+
+``` python
+in_list = [nb.load(f) for f in glob.glob('data/cobre/fmri_*.nii.gz')]
+in_mask = nb.load('data/brain_mask.nii.gz')
+go = '3DMIGP'
+dPCA_int = 550
+dPCA_out = 100
+```
+
+%% Cell type:code id: tags:
+
+``` python
+random.shuffle(in_list)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+mask = in_mask.get_fdata().ravel()
+
+def demean(x):
+    return x - np.mean(x, axis=0)
+
+for i, f in enumerate(in_list):
+
+    print(i, end=' ')
+
+    # read data
+    grot = f.get_fdata()
+    grot = np.reshape(grot, [-1, grot.shape[-1]])
+    grot = grot[mask!=0, :].T
+    grot = demean(grot)
+
+    # var-norm
+    [uu, ss, vt] = svds(grot, k=30)
+    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)
+    grot = grot/stddevs
+
+    if i == 0:
+        W = demean(grot)
+    else:
+        # concat
+        W = np.concatenate((W, demean(grot)), axis=0)
+
+        # reduce W to dPCA_int eigenvectors
+        if W.shape[0]-10 > dPCA_int:
+            uu = eigs(W@W.T, dPCA_int)[1]
+            uu = np.real(uu)
+            W = uu.T @ W
+
+    f.uncache()
+
+```
+
+%% Cell type:code id: tags:
+
+``` python
+grot = np.zeros([mask.shape[0], dPCA_out])
+grot[mask!=0, :] = W[:dPCA_out, :].T
+grot = np.reshape(grot, in_list[0].shape[:3] + (dPCA_out,))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+nb.Nifti1Image(grot, affine=in_list[0].affine, header=in_list[0].header).to_filename('pyMIGP.nii.gz')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+%%bash
+
+melodic -i pyMIGP.nii.gz \
+	--mask=data/brain_mask.nii.gz \
+	-d 20 \
+	-v \
+	--nobet \
+	--disableMigp \
+	--varnorm \
+	-o pymigp_dim20.ica
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ics = nb.load('pymigp_dim20.ica/melodic_IC.nii.gz')
+
+N = ics.shape[-1]
+
+fig, ax = plt.subplots(int(np.ceil((N/2))),2, figsize=(12, N))
+
+for img, ax0 in zip(image.iter_img(ics), ax.ravel()):
+    coord = plotting.find_xyz_cut_coords(img, activation_threshold=1.5)
+    plotting.plot_stat_map(img, cut_coords=coord, vmax=5, axes=ax0)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```