Alternative versions: LabNotebook, Nbviewer

Format: straight-up

Overview¶

I've been doing a lot of work recently with cortical surface laplace-beltrami operator (LBO) eigenmodes. The LBO modes are a set of spherical harmonic-like basis functions that can be used to describe both the shape and the spatiotemporal dynamics of neural activity on the cortex.

There are nice tools from moo chung and martin reuter + arno klein, in matlab and python respectively for computing the LBO from arbitrary surface meshes, both of which work very well for human cortical surface data. Strangely, however, it appears that both of these tools fail on monkey (rhesus macaque) cortical surface meshes.

Initially the following notes were simply a report on this problem and essentially a plea for help. However I'm happy to report I've now found a relatively simple regulariation-based solution.

Notebook Setup¶

Define some variables

# put some system-specific variables in the namespace ('le' dict)
%run ~/set_localenv_vars.py

nb_name = 'macaque_lbo_modes_problem'

# this is where I look for things
outdir = '%s/notebooks/%s' %(le['data_dir'],nb_name)


# this is where things actually live
outdir_loc = '%s/%s' %('/mnt/titania-hecuba/mcintosh_lab/john/Data/notebooks', nb_name)
#outdir_loc = '%s/%s' %('/hecuba/mcintosh_lab/john/Data/notebooks', nb_name)
import os,sys
if not os.path.isdir(outdir_loc): os.makedirs(outdir_loc)
if not os.path.isdir(outdir): os.system('ln -s %s %s' %(outdir_loc,outdir))
    
!mkdir $outdir/pics
    
# freesurfer surfaces
fsav5_dir = '/software/freesurfer/subjects/fsaverage5'
fsav5_lhp_file = fsav5_dir + '/surf/lh.pial'
fsav5_lhc_file = fsav5_dir + '/surf/lh.curv'


# macaque datasets: 

# 1. cmtk library example data
m1_f = '/mnt/rotmangrid-scratch2/mcintosh_lab/jgriffiths/Data/cmtk_data/macaca_mulatta_01.cff'

# 2. f99 macaque data (from gleb bezgin)
m2_f = '/mnt/rotmangrid-scratch2/mcintosh_lab/jgriffiths/Data/gleb_data/f99_vertices_faces.mat'

# 3. NIMH macaque template https://github.com/jms290/NMT/tree/master/NMT_v1.2
m3_f = '/alexandra/mcintosh_lab/john/Code/libraries_of_others/github/NMT/NMT_v1.2/lh.GM.gii'



mindboggle_dir = le['code_dir'] + '/libraries_of_others/github/mindboggle'
sys.path.append(mindboggle_dir)


# stuff for workdocs-cloudfiles

aws_key = 'drowssaperucesyreva'
aws_secret = '?teytidettopsuoyevah'

Importage

# Generic imports

import os,sys,glob,numpy as np,pandas 
from scipy.sparse.linalg import eigsh,lobpcg
from scipy.io import loadmat


# Neuroimaging stuff

import nibabel as nib
from mindboggle.shapes.laplace_beltrami import computeAB
import cfflib

# Visualization stuff

%matplotlib inline
from matplotlib import pyplot as plt
from IPython.display import Image,display as d,clear_output
from nilearn.plotting import plot_surf_stat_map


# aws api and workdocs-cloudfiles stuff

sys.path.append(le['ipynb_workdocs_dir'])
from new_utils import nb_fig,cloudfiles_nb

Calico document extensions

%%javascript
IPython.load_extensions('calico-spell-check', 
                        'calico-document-tools', 'calico-cell-tools');

Initialize aws api and workdocs-cloudfiles folder

cnb = cloudfiles_nb('aws', [aws_key,aws_secret])
cnb.initialize_folder(nb_name)

Go to output folder

os.chdir(outdir)

Ok, let's get cracking

Load data and look at surfaces¶

fsaverage 5 human data

fsav5_lhp_vtx,fsav5_lhp_tri = nib.freesurfer.read_geometry(fsav5_lhp_file)

macaque dataset 1

res_lh,res_rh = cfflib.load(m1_f).get_connectome_surface()
res_lh.load()
m1_vtx = res_lh.data.darrays[0].data
m1_tri = res_lh.data.darrays[1].data

Loading file. Created temporary file: /tmp/macaque_connectivity_right_hemisphere514709/Gifti/macaca01.gii
Succeed.

macaque dataset 2

f2_mat = loadmat(m2_f, struct_as_record=False,squeeze_me=True)
m2_vtx = f2_mat['vertices_original']
m2_tri = f2_mat['faces']
m2_tri-=1

macaque dataset 3

img = nib.gifti.giftiio.read(m3_f)
m3_vtx = img.darrays[0].data
m3_tri = img.darrays[1].data

So we have one human cortical surface and three different macaque cortical surfaces. Let's take a look at these together:

(0_0)

f = outdir + '/pics/all_surfs.png'


fig = plot_surf([fsav5_lhp_vtx,fsav5_lhp_tri],view='lateral');
plt.savefig('test_pic1.png', bbox_inches='tight', transparent=True)
plt.close()

fig = plot_surf([fsav5_lhp_vtx,fsav5_lhp_tri],view='medial');
plt.savefig('test_pic2.png', bbox_inches='tight', transparent=True)
plt.close()

fig = plot_surf([m1_vtx,m1_tri],view='lateral');
plt.savefig('test_pic3.png', bbox_inches='tight', transparent=True)
plt.close()

fig = plot_surf([m1_vtx,m1_tri],view='medial');
plt.savefig('test_pic4.png', bbox_inches='tight', transparent=True)
plt.close()


fig = plot_surf([m2_vtx,m2_tri],view='lateral');
plt.savefig('test_pic5.png', bbox_inches='tight', transparent=True)
plt.close()

fig = plot_surf([m2_vtx,m2_tri],view='medial');
plt.savefig('test_pic6.png', bbox_inches='tight', transparent=True)
plt.close()


fig = plot_surf([m3_vtx,m3_tri],view='lateral');
plt.savefig('test_pic7.png', bbox_inches='tight', transparent=True)
plt.close()

fig = plot_surf([m3_vtx,m3_tri],view='medial');
plt.savefig('test_pic8.png', bbox_inches='tight', transparent=True)
plt.close()




fs = ['test_pic%s.png' %t for t in range(1,9)]

fig = plt.figure(figsize=(12,12))

ax1 = fig.add_axes([0,   1-0.2,0.25,0.25])
ax1.imshow(plt.imread(fs[0])[50:-50,80:-70])

ax2 = fig.add_axes([0.25,1-0.2,0.25,0.25])
ax2.imshow(plt.imread(fs[1])[50:-50,80:-70])

ax3 = fig.add_axes([0,   1-0.35,0.25,0.25])
ax3.imshow(plt.imread(fs[2])[50:-50,80:-70])

ax4 = fig.add_axes([0.25,1-0.35,0.25,0.25])
ax4.imshow(plt.imread(fs[3])[50:-50,80:-70])

ax5 = fig.add_axes([0,1-0.50,0.25,0.25])
ax5.imshow(plt.imread(fs[4])[50:-50,80:-70])

ax6 = fig.add_axes([0.25,1-0.50,0.25,0.25])
ax6.imshow(plt.imread(fs[5])[50:-50,80:-70])

ax7 = fig.add_axes([0,1-0.65,0.25,0.25])
ax7.imshow(plt.imread(fs[6])[50:-50,80:-70])

ax8 = fig.add_axes([0.25,1-0.65,0.25,0.25])
ax8.imshow(plt.imread(fs[7])[50:-50,80:-70])


for a in [ax1,ax2,ax3,ax4,ax5,ax6,ax7,ax8]: a.axis('off')

plt.savefig(f, bbox_inches='tight', transparent=True)
plt.close()



# Embed in LN
cap = ''; label = ''
im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500));
d(im)

.

Reference point: human brain LBO modes¶

Here is how the human brain cortical surface LBO modes look. This is what I've trying to get from the macaque data. Showing here modes 2-4, whic suffice to give the general picture of what we're looking for. The first mode is uniform and not very informative to look at, so we skip that one.

fsav5_lhp_A,fsav5_lhp_B = computeAB(fsav5_lhp_vtx,fsav5_lhp_tri)
fsav5_lhp_lbo_evals,fsav5_lhp_lbo_evecs = eigsh(fsav5_lhp_A,M=fsav5_lhp_B,sigma=0,k=10)

(0_0)

.

The problem: doesn't work for macaque¶

Ok so this is the error that we get when trying to compute LBO modes as above from macaque surfaces. The mindboggle computeAB function seems to run ok, but the sparse eigendecomposition with eigsh fails because A is singular.

macaque dataset 1

m1_A,m1_B = computeAB(m1_vtx,m1_tri)
m1_lbo_evals,m1_lbo_evecs = eigsh(m1_A,M=m1_B,sigma=0,k=10)


---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-204-9938091651e0> in <module>()
      1 m1_A,m1_B = computeAB(m1_vtx,m1_tri)
----> 2 m1_lbo_evals,m1_lbo_evecs = eigsh(m1_A,M=m1_B,sigma=0,k=10)

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in eigsh(A, k, M, sigma, which, v0, ncv, maxiter, tol, return_eigenvectors, Minv, OPinv, mode)
   1543             if OPinv is None:
   1544                 Minv_matvec = get_OPinv_matvec(A, M, sigma,
-> 1545                                                symmetric=True, tol=tol)
   1546             else:
   1547                 OPinv = _aslinearoperator_with_dtype(OPinv)

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in get_OPinv_matvec(A, M, sigma, symmetric, tol)
   1020 def get_OPinv_matvec(A, M, sigma, symmetric=False, tol=0):
   1021     if sigma == 0:
-> 1022         return get_inv_matvec(A, symmetric=symmetric, tol=tol)
   1023 
   1024     if M is None:

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in get_inv_matvec(M, symmetric, tol)
   1013         if isspmatrix_csr(M) and symmetric:
   1014             M = M.T
-> 1015         return SpLuInv(M).matvec
   1016     else:
   1017         return IterInv(M, tol=tol).matvec

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/interface.pyc in __new__(cls, *args, **kwargs)
    140                                 " at least one of _matvec and _matmat.")
    141 
--> 142             obj.__init__(*args, **kwargs)
    143             return obj
    144 

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in __init__(self, M)
    897     """
    898     def __init__(self, M):
--> 899         self.M_lu = splu(M)
    900         self.shape = M.shape
    901         self.dtype = M.dtype

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.pyc in splu(A, permc_spec, diag_pivot_thresh, drop_tol, relax, panel_size, options)
    255         _options.update(options)
    256     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
--> 257                           ilu=False, options=_options)
    258 
    259 

RuntimeError: Factor is exactly singular

macaque dataset 2

m2_A,m2_B = computeAB(m2_vtx,m2_tri)
m2_lbo_evals,m2_lbo_evecs = eigsh(m2_A,M=m2_B,sigma=0,k=10)


---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-208-d917db00a7e1> in <module>()
      1 m2_A,m2_B = computeAB(m2_vtx,m2_tri)
----> 2 m2_lbo_evals,m2_lbo_evecs = eigsh(m2_A,M=m2_B,sigma=0,k=10)

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in eigsh(A, k, M, sigma, which, v0, ncv, maxiter, tol, return_eigenvectors, Minv, OPinv, mode)
   1543             if OPinv is None:
   1544                 Minv_matvec = get_OPinv_matvec(A, M, sigma,
-> 1545                                                symmetric=True, tol=tol)
   1546             else:
   1547                 OPinv = _aslinearoperator_with_dtype(OPinv)

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in get_OPinv_matvec(A, M, sigma, symmetric, tol)
   1020 def get_OPinv_matvec(A, M, sigma, symmetric=False, tol=0):
   1021     if sigma == 0:
-> 1022         return get_inv_matvec(A, symmetric=symmetric, tol=tol)
   1023 
   1024     if M is None:

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in get_inv_matvec(M, symmetric, tol)
   1013         if isspmatrix_csr(M) and symmetric:
   1014             M = M.T
-> 1015         return SpLuInv(M).matvec
   1016     else:
   1017         return IterInv(M, tol=tol).matvec

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/interface.pyc in __new__(cls, *args, **kwargs)
    140                                 " at least one of _matvec and _matmat.")
    141 
--> 142             obj.__init__(*args, **kwargs)
    143             return obj
    144 

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in __init__(self, M)
    897     """
    898     def __init__(self, M):
--> 899         self.M_lu = splu(M)
    900         self.shape = M.shape
    901         self.dtype = M.dtype

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.pyc in splu(A, permc_spec, diag_pivot_thresh, drop_tol, relax, panel_size, options)
    255         _options.update(options)
    256     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
--> 257                           ilu=False, options=_options)
    258 
    259 

RuntimeError: Factor is exactly singular

macaque dataset 3

m3_A,m3_B = computeAB(m3_vtx,m3_tri)
m3_lbo_evals,m3_lbo_evecs = eigsh(m3_A,M=m3_B,sigma=0,k=10)


---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-210-246be7d601b0> in <module>()
      1 m3_A,m3_B = computeAB(m3_vtx,m3_tri)
----> 2 m3_lbo_evals,m3_lbo_evecs = eigsh(m3_A,M=m3_B,sigma=0,k=10)

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in eigsh(A, k, M, sigma, which, v0, ncv, maxiter, tol, return_eigenvectors, Minv, OPinv, mode)
   1543             if OPinv is None:
   1544                 Minv_matvec = get_OPinv_matvec(A, M, sigma,
-> 1545                                                symmetric=True, tol=tol)
   1546             else:
   1547                 OPinv = _aslinearoperator_with_dtype(OPinv)

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in get_OPinv_matvec(A, M, sigma, symmetric, tol)
   1020 def get_OPinv_matvec(A, M, sigma, symmetric=False, tol=0):
   1021     if sigma == 0:
-> 1022         return get_inv_matvec(A, symmetric=symmetric, tol=tol)
   1023 
   1024     if M is None:

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in get_inv_matvec(M, symmetric, tol)
   1013         if isspmatrix_csr(M) and symmetric:
   1014             M = M.T
-> 1015         return SpLuInv(M).matvec
   1016     else:
   1017         return IterInv(M, tol=tol).matvec

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/interface.pyc in __new__(cls, *args, **kwargs)
    140                                 " at least one of _matvec and _matmat.")
    141 
--> 142             obj.__init__(*args, **kwargs)
    143             return obj
    144 

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.pyc in __init__(self, M)
    897     """
    898     def __init__(self, M):
--> 899         self.M_lu = splu(M)
    900         self.shape = M.shape
    901         self.dtype = M.dtype

/alexandra/mcintosh_lab/john/Software/miniconda2/envs/_ipython2.4b/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.pyc in splu(A, permc_spec, diag_pivot_thresh, drop_tol, relax, panel_size, options)
    255         _options.update(options)
    256     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
--> 257                           ilu=False, options=_options)
    258 
    259 

RuntimeError: Factor is exactly singular

Alternative computation 1: doesn't work¶

The mindboggle fem_laplacian.py code has an alternative for when this kind of error is encountered. Here's a snippet:


# -----------------------------------------------------------
# Use the eigsh eigensolver:
# -----------------------------------------------------------
try :
  # eigs is for nonsymmetric matrices while
  # eigsh is for real-symmetric or complex-Hermitian matrices:
  eigenvalues, eigenvectors = eigsh(A, k=spectrum_size,
                                    M=B,sigma=-0.01)
  spectrum = eigenvalues.tolist()

# -----------------------------------------------------------
# Use the lobpcg eigensolver:
# -----------------------------------------------------------
except RuntimeError:     

  if verbose:
    print("eigsh() failed. Now try lobpcg.")
    print("Warning: lobpcg can produce" 
          "different results from "
          "Reuter (2006) shapeDNA-tria software.")
  # Initial eigenvector values:
  init_eigenvecs = np.random.random((A.shape[0],
                                     spectrum_size))

  # maxiter = 40 forces lobpcg to use 20 iterations.
  # Strangely, largest=false finds largest eigenvalues
  # and largest=True gives the smallest eigenvalues:
  eigenvalues, eigenvectors = lobpcg(A, init_eigenvecs,B=B,
                                   largest=True, maxiter=40)

...etc etc

Let's try this for the first macaque dataset:

spectrum_size = 5
A,B = computeAB(m1_vtx,m1_tri)
init_evecs = np.random.random((A.shape[0],
                                spectrum_size))
m1_lobpcg_lbo_evals,m1_lobpcg_lbo_evecs = lobpcg(A, init_evecs,B=B,
                     largest=True, maxiter=40)

Unfortunately the outputs of this computation are clearly not what we're after; which fortunately we can judge because we have some geometric intuitions about what LBO eigenmodes 'should' look like.

(0_0)

.

They should not look like that.

Also, using that on human brains doesn't look right either. So looks like the lobpcg computation should be avoided for this particular application.

Alternative computation 2: also doesn't work¶

Second alternative I tried is moo chung's matlab code. Won't show that in full detail here; but it's pretty simple to use. Goes something like this:


load(f);

mac.vertices = vertices_orig;
mac.faces = faces;

[A, C] =FEM(mac);

[V, D] = eigs(C,A,999,'sm');

But this also fails at the eigedecomposition stage.

Solution: regularization¶

It's clear from the error message above RuntimeError: Factor is exactly singular that there is a problem inverting the A matrix because it is singular and/or poorly conditioned. The error message from the matlab option (not shown) complains about the same thing.

The solution to this problem isn't too hard to find: the top google hit for the "python RuntimeError: Factor is exactly singular" has a comment reply that points to this stackoverflow post on improving a badly conditioned matrix.

The discussion there receommends adding a small constant c to the diagonal of the matrix.

A simple way to find the best value of c is simply to start with some small starting value that works (o.e. when added to the diagonal doesn't throw up the singular matrix error), and then suceessively reduce the value of c and until the you start getting the error again.

Here's a dumb, very simple way to do that. I find the best (/ an acceptable) regularization factor by decreasing c until the eigsh call throws the error it was doing above


blah = 500
while blah > 100:
  c = c/10.
  print 'c=%s' %c
  m1_A,m1_B = computeAB(m1_vtx,m1_tri)
  m1_A_mod = m1_A.copy()
  for a in np.arange(m1_A_mod.shape[0]): m1_A_mod[a,a] += c
  m1_lbo_evals,m1_lbo_evecs = eigsh(m1_A_mod,M=m1_B,
                                    sigma=0,k=10)

c = 1E-50
m1_A,m1_B = computeAB(m1_vtx,m1_tri)
m1_A_mod = m1_A.copy()
for a in np.arange(m1_A_mod.shape[0]): m1_A_mod[a,a] += c
m1_lbo_evals,m1_lbo_evecs = eigsh(m1_A_mod,M=m1_B,sigma=0,k=10)

(0_0)

.

c = 1E-50
m2_A,m2_B = computeAB(m2_vtx,m2_tri)
m2_A_mod = m2_A.copy()
for a in np.arange(m2_A_mod.shape[0]): m2_A_mod[a,a] += c
m2_lbo_evals,m2_lbo_evecs = eigsh(m2_A_mod,M=m2_B,sigma=0,k=10)

(0_0)

.

c = 1E-20
m3_A,m1_B = computeAB(m3_vtx,m3_tri)
m3_A_mod = m3_A.copy()
for a in np.arange(m3_A_mod.shape[0]): m3_A_mod[a,a] += c
m3_lbo_evals,m3_lbo_evecs = eigsh(m3_A_mod,M=m3_B,sigma=0,k=10)

(0_0)

.

Ok so in the third case there the modes aren't looking so good as the first two. Possibly that can be fixed with a bit more of a refine search for a good value of c.

That aside, bonzer: we have functioning LBO decompositions for macaque cortical surface data :)

JDG Lab Notebook

Macaque LBO modes problem

Contents

Overview¶

Notebook Setup¶

Load data and look at surfaces¶

Reference point: human brain LBO modes¶

The problem: doesn't work for macaque¶

Alternative computation 1: doesn't work¶

Alternative computation 2: also doesn't work¶

Solution: regularization¶