calc_pcorr_Manyneighbours.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import nibabel as nb\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#img=nb.load('BP4T_rest1_1_denoised_z35.nii.gz')\n",
    "#img=nb.load('BP4T_rest1_1_denoised.nii.gz')\n",
    "#img=nb.load('BP4T_rest1_1_denoised_z20_49.nii.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## from: https://gist.github.com/fabianp/9396204419c7b638d38f\n",
    "\n",
    "\"\"\"\n",
    "Partial Correlation in Python (clone of Matlab's partialcorr)\n",
    "This uses the linear regression approach to compute the partial \n",
    "correlation (might be slow for a huge number of variables). The \n",
    "algorithm is detailed here:\n",
    "    http://en.wikipedia.org/wiki/Partial_correlation#Using_linear_regression\n",
    "Taking X and Y two variables of interest and Z the matrix with all the variable minus {X, Y},\n",
    "the algorithm can be summarized as\n",
    "    1) perform a normal linear least-squares regression with X as the target and Z as the predictor\n",
    "    2) calculate the residuals in Step #1\n",
    "    3) perform a normal linear least-squares regression with Y as the target and Z as the predictor\n",
    "    4) calculate the residuals in Step #3\n",
    "    5) calculate the correlation coefficient between the residuals from Steps #2 and #4; \n",
    "    The result is the partial correlation between X and Y while controlling for the effect of Z\n",
    "Date: Nov 2014\n",
    "Author: Fabian Pedregosa-Izquierdo, f@bianp.net\n",
    "Testing: Valentina Borghesani, valentinaborghesani@gmail.com\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "from scipy import stats, linalg\n",
    "\n",
    "def partial_corr(C):\n",
    "    \"\"\"\n",
    "    Returns the sample linear partial correlation coefficients between pairs of variables in C, controlling \n",
    "    for the remaining variables in C.\n",
    "    Parameters\n",
    "    ----------\n",
    "    C : array-like, shape (n, p)\n",
    "        Array with the different variables. Each column of C is taken as a variable\n",
    "    Returns\n",
    "    -------\n",
    "    P : array-like, shape (p, p)\n",
    "        P[i, j] contains the partial correlation of C[:, i] and C[:, j] controlling\n",
    "        for the remaining variables in C.\n",
    "    \"\"\"\n",
    "    \n",
    "    C = np.asarray(C)\n",
    "    p = C.shape[1]\n",
    "    P_corr = np.zeros((p, p), dtype=np.float)\n",
    "    for i in range(p):\n",
    "        P_corr[i, i] = 1\n",
    "        for j in range(i+1, p):\n",
    "            idx = np.ones(p, dtype=np.bool)\n",
    "            idx[i] = False\n",
    "            idx[j] = False\n",
    "            beta_i = linalg.lstsq(C[:, idx], C[:, j])[0]\n",
    "            beta_j = linalg.lstsq(C[:, idx], C[:, i])[0]\n",
    "\n",
    "            res_j = C[:, j] - C[:, idx].dot( beta_i)\n",
    "            res_i = C[:, i] - C[:, idx].dot(beta_j)\n",
    "            \n",
    "            corr = stats.pearsonr(res_i, res_j)[0]\n",
    "            P_corr[i, j] = corr\n",
    "            P_corr[j, i] = corr\n",
    "        \n",
    "    return P_corr\n",
    "\n",
    "def calc_vox_pcorr_dirs(d,m,nbours):\n",
    "    \n",
    "    d_shape=np.shape(d)\n",
    "    res=np.zeros(d_shape[0:3]+(np.shape(nbours)[0]*2,)) #x,y,z and the neighbourhoods as the 4th dim, mult*2 b/c +/-\n",
    "    xdim=d_shape[0]\n",
    "    ydim=d_shape[1]\n",
    "    zdim=d_shape[2]\n",
    "    for x in range(0+1,xdim-1):\n",
    "        for y in range(0+1,ydim-1):\n",
    "            for z in range(0+1,zdim-1):\n",
    "                for idx,nbour in enumerate(nbours):\n",
    "                    nidx = idx*2\n",
    "                    loc = [x,y,z] + nbour\n",
    "                    locm = [x,y,z] + nbour*-1\n",
    "#                    corrmat=np.hstack((np.ndarray.flatten(d[tuple(loc)]).reshape(1,-1).T,\n",
    "#                                   np.ndarray.flatten(d[x,y,z]).reshape(1,-1).T,\n",
    "#                                   np.ndarray.flatten(d[tuple(locm)]).reshape(1,-1).T))\n",
    "#                    res_corr=partial_corr(corrmat)\n",
    "                    \n",
    "#                    res[x,y,z,nidx]=res_corr[1,0] #pcorr btw the first two\n",
    "#                    res[x,y,z,nidx+1]=res_corr[1,2] #pcorr btw the last two\n",
    "                    \n",
    "                    print(nidx)\n",
    "    return res\n",
    "                    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "0\n",
      "1\n",
      "0\n",
      "1\n",
      "0\n",
      "1\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "for idx,a in enumerate(range(1,10)):\n",
    "    #print(idx)\n",
    "    print(idx%2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## setup neighbourhood matrices to use as loop for calculating the partial correlations\n",
    "\n",
    "lval = 1 #value of voxels to lookup (or down, sideways etc), mult by n to increase range\n",
    "\n",
    "# setup so that we can mult by -1 to give the other side of neighbourhood\n",
    "neighbours_6=np.array([[lval, 0, 0],            [0,lval,0],             [0,0,lval]])\n",
    "neighbours_12=np.array([[lval, 0, 0],           [0,lval,0],             [0,0,lval],\n",
    "                       [lval,lval,0],           [lval,0,lval],          [0,lval,lval],\n",
    "                       [lval,lval*-1,0],        [lval,0,lval*-1],       [0,lval,lval*-1]])\n",
    "#currently wrong\n",
    "neighbours_24=np.array([[lval, 0, 0],           [0,lval,0],             [0,0,lval],\n",
    "                        [lval,lval,0],          [lval,0,lval],          [0,lval,lval],\n",
    "                        [lval,lval*-1,0],       [lval,0,lval*-1],       [0,lval,lval*-1],\n",
    "                        [lval,lval*-1,lval],    [lval,lval,lval*-1],    [lval,lval,lval],\n",
    "                        [lval*-1,lval*-1,lval], [lval,lval*-1,lval*-1], [lval*-1,lval,lval*-1]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "2\n",
      "4\n",
      "0\n",
      "2\n",
      "4\n",
      "0\n",
      "2\n",
      "4\n",
      "0\n",
      "2\n",
      "4\n",
      "0\n",
      "2\n",
      "4\n",
      "0\n",
      "2\n",
      "4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[[[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]]],\n",
       "\n",
       "\n",
       "       [[[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]]],\n",
       "\n",
       "\n",
       "       [[[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]],\n",
       "\n",
       "        [[ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.],\n",
       "         [ 0.,  0.,  0.,  0.,  0.,  0.]]]])"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_vox_pcorr_dirs(np.zeros((3,5,4)),[1],neighbours_6)\n",
    "#np.zeros((3,5,4))[tuple([0,0,0])]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d=img.get_data() #this takes a while on a full image\n",
    "a=img.affine\n",
    "mask=d[:,:,:,0]\n",
    "mask[np.logical_not(mask==0)]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7322.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "128 128 1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "xdim=np.shape(d)[0]\n",
    "ydim=np.shape(d)[1]\n",
    "zdim=np.shape(d)[2] #not used here, since I only loaded a single slice to check this\n",
    "print xdim,ydim,zdim\n",
    "print \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nbours=6\n",
    "res = np.zeros(np.shape(d)[0:3]+(nbours,)) #x,y,z and the neighbourhoods\n",
    "print(np.shape(res))\n",
    "# TODO: get rid of enumerate\n",
    "print(\"x slice:\")\n",
    "for idxx,x in enumerate(range(0+1,xdim-1)):\n",
    "    print x\n",
    "    for idxy,y in enumerate(range(0+1,ydim-1)):\n",
    "        for idxz,z in enumerate(range(0+1,zdim-1)):\n",
    "        #    if x<0 or y<0 or z<0 or x>xdim-1 or y>ydim-1 or z>zdim-1:\n",
    "        #        pass\n",
    "        #    else:\n",
    "            if not(mask[x,y,z]==0):\n",
    "                #print x,y,z\n",
    "                #for now, 2 directions of correlations\n",
    "                corrmat=np.hstack((np.ndarray.flatten(d[x-1,y,z]).reshape(1,-1).T,\n",
    "                                   np.ndarray.flatten(d[x,y,z]).reshape(1,-1).T,\n",
    "                                   np.ndarray.flatten(d[x+1,y,z]).reshape(1,-1).T))\n",
    "                res_corr=partial_corr(corrmat)\n",
    "                res[x,y,z,0]=res_corr[1,0] #pcorr btw the first two\n",
    "                res[x,y,z,1]=res_corr[1,2] #pcorr btw the last two\n",
    "\n",
    "                corrmat=np.hstack((np.ndarray.flatten(d[x,y-1,z]).reshape(1,-1).T,\n",
    "                                   np.ndarray.flatten(d[x,y,z]).reshape(1,-1).T,\n",
    "                                   np.ndarray.flatten(d[x,y+1,z]).reshape(1,-1).T))\n",
    "                res_corr=partial_corr(corrmat)\n",
    "                res[x,y,z,2]=res_corr[1,0] #pcorr btw the first two\n",
    "                res[x,y,z,3]=res_corr[1,2] #pcorr btw the last two    \n",
    "\n",
    "                corrmat=np.hstack((np.ndarray.flatten(d[x,y,z-1]).reshape(1,-1).T,\n",
    "                                   np.ndarray.flatten(d[x,y,z]).reshape(1,-1).T,\n",
    "                                   np.ndarray.flatten(d[x,y,z+1]).reshape(1,-1).T))\n",
    "                res_corr=partial_corr(corrmat)\n",
    "                res[x,y,z,4]=res_corr[1,0] #pcorr btw the first two\n",
    "                res[x,y,z,5]=res_corr[1,2] #pcorr btw the last two \n",
    "print \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res[np.isnan(res)]=0\n",
    "out_img=nb.Nifti1Image(res,a,header=img.get_header())\n",
    "nb.save(out_img,'6neighbours.nii.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Try for some parallel action?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "def calc_vox_vectors(d,mask,x,y,z):\n",
    "    #calc the partial corr between each pair in the neighbourhood, controlling for the diagonal other\n",
    "    res=np.zeros((1,6))\n",
    "    if not(mask[x,y,z]==0):\n",
    "        corrmat=np.hstack((np.ndarray.flatten(d[x-1,y,z]).reshape(1,-1).T,\n",
    "                           np.ndarray.flatten(d[x,y,z]).reshape(1,-1).T,\n",
    "                           np.ndarray.flatten(d[x+1,y,z]).reshape(1,-1).T))\n",
    "        res_corr=partial_corr(corrmat)\n",
    "        res[0]=res_corr[0,1] #pcorr btw the first two\n",
    "        res[1]=res_corr[1,2] #pcorr btw the last two\n",
    "\n",
    "        corrmat=np.hstack((np.ndarray.flatten(d[x,y-1,z]).reshape(1,-1).T,\n",
    "                           np.ndarray.flatten(d[x,y,z]).reshape(1,-1).T,\n",
    "                           np.ndarray.flatten(d[x,y+1,z]).reshape(1,-1).T))\n",
    "        res_corr=partial_corr(corrmat)\n",
    "        res[2]=res_corr[0,1] #pcorr btw the first two\n",
    "        res[3]=res_corr[1,2] #pcorr btw the last two    \n",
    "\n",
    "        corrmat=np.hstack((np.ndarray.flatten(d[x,y,z-1]).reshape(1,-1).T,\n",
    "                           np.ndarray.flatten(d[x,y,z]).reshape(1,-1).T,\n",
    "                           np.ndarray.flatten(d[x,y,z+1]).reshape(1,-1).T))\n",
    "        res_corr=partial_corr(corrmat)\n",
    "        res[4]=res_corr[0,1] #pcorr btw the first two\n",
    "        res[5]=res_corr[1,2] #pcorr btw the last two\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#!/usr/bin/env python2\n",
    "import itertools\n",
    "from multiprocessing import Pool, freeze_support\n",
    "img=nb.load('BP4T_rest1_1_denoised_z20_49.nii.gz')\n",
    "d=img.get_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "pool = Pool()\n",
    "xargs=range(0+1,xdim-1)\n",
    "yargs=range(0+1,ydim-1)\n",
    "zargs=range(0+1,zdim-1)\n",
    "a= pool.map(calc_vox_vectors, itertools.izip(d,mask,xargs,yargs,zargs))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}