JDG Lab Notebook

TVB SC-FC Model Optimization

Overview

Fitting TVB simulations to empirical data of various kinds is a major area of current development, in concert with related issues such as improving simulation speed. The situation with both of these is likely to change considerably over the next few years, but at present there is basically no single canonical approach to model fitting using the core TVB code library.

The following notes assess a couple of relatively simple working solutions using available python tools and the standard python TVB code.

A major consideration point is the factor of execution time, and getting a solution that runs in a reasonable amount of time. 'Reasonable' here is very much a moving goalpost, as it depends on model objectives, available resources, etc. One thing that is certain is that fully global optimization is not an option, even with hardcore computational beef. So what we are looking for are heuristics, compromises, and good-enoughs. Parallelization helps a lot here: whilst the core TVB code is not multi-threaded, but we can make use of paralellization over multiple cores and/or multiple CPUs to speed things up significantly. Generally speaking, however, optimization techniques such as gradient descent are very serial operations. Possibly the magic combination for us may be tools such as the playdoh library or spearmint library that support distributed (over cores and CPUs) optimization.

A second major consideration that I will just flag up now is that I am not attempting a proper technical treatment of the problem. Model optimization is a highly developed field of scientific computing and applied mathematics, and there are a number of core components of optimization solutions (e.g. properly evaluating partial derivatives, jacobians, tolerances, choosing most appropriate algorithms, objective functions, etc.) that I am entirely glossing over here. That said, the scipy algorithms do take care of a lot of these factors automatically to some extent, and also allow a lot of additional customization along these lines if and when such evaluations can be made.

To start with, however, this is just a first pass flyover, to assess the lay of the land, and see how well we can do with simple solutions.

Notebook Setup

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

nb_name = 'tvb_sc_fc_model_optimization'

outdir = '%s/notebooks/%s' %(le['data_dir'], nb_name)
!mkdir -p $outdir/pics

# root path to tvb-library, tvb-data, etc. libraries. 
tvb_folder = le['code_dir'] + '/libraries_of_others/github/tvb-pack'

# stuff for workdocs-cloudfiles

aws_key = 'drowssaperucesyreva'
aws_secret = '?teytidettopsuoyevah'  
aws_key = 'AKIAIQG63YDHSA4UJ74A'
aws_secret = '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'
# Generic imports

import os,sys,glob,h5py,itertools,multiprocessing,\
       numpy as np,pandas as pd
from datetime import datetime
from time import time
from scipy import optimize

# Visualization stuff

%matplotlib inline
from matplotlib import pyplot as plt
from IPython.display import Image,display as d
import seaborn as sns


# aws api and workdocs-cloudfiles stuff

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


# TVB stuff

nbso,nbse = sys.stdout,sys.stderr # hack part 1/2 to keep output printing properly
sys.path.append(tvb_folder +'/tvb-library')
sys.path.append(tvb_folder + '/tvb-data')
from tvb.simulator.lab import *
sys.stdout,sys.stderr = nbso,nbse  # ...hack part 2/2


from copy import copy # seems that this needs to be done last otherwise it gets 
                      # replaced with numpy copy function

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

Overview of the models

At present the need for an optimization solution arises most often in the context of fitting models of resting state functional connectivity (FC). In this simplest case this involves tuning TVB model parameters such to maximize the correlation between a simulated FC matrix and an empirical FC matrix. Other scenarios we shall consider later are fitting M/EEG phenomena such as power spectra and ERPs.

However to simplify things, for present purposes we can drop the empirical data component, and instead try to optimize the match between simulated FC and structural connectivity (SC). This is demonstrated very nicely in an example demo by the Duke, and the following is based off of that example. Whilst it is a little circular, matching simulated FC to SC is actually a useful approach, because empirical FC matrices normally have a fairly high correlation with empirical SC matrices, and this is not the case for large domains of parameter space. So, minimally, matching simulated FC to SC gets our models into physiologically plausible operating regimes, and is enough of an optimization problem to test out some methods.

It's also important to stress that the following FC matrices are computed from M/EEG/LFP time scale data, so they are not models of fMRI functional connectivity per se. Simulating and fitting a typical 5 minute resting state fMRI scan takes at least an order of magnitude longer than the 500-2000ms examples below, and so isn't really an option, and so isn't really an option for the optimization strategies explored here.

That said, in my experience you get a very good idea of what long-duration covariance patterns will look like from shorter simulation runs. Quantifying that relationship is an interesting and potentially extremely practical problem that we will leave for future work.

The following functions are based closely on the examples in Marmaduke's demo, extended for a little more flexibility in model and parameter specification:

# Simulation workhorse function

def run_sim(gca=0.3,spd=3.,nsig=1,dt=0.1,                # default global params
            mt = 'G2DO', a=0., 
            alpha= 1.0,beta=1.0,
            k11=0.5,k12=0.5,k21=0.5,                     # default model params
            sim_len=500, returntime = True):
    
  tic = time()
    
  #if mt == 'G2DO':  mod = models.Generic2dOscillator()        
  if mt == 'G2DO':  mod = models.Generic2dOscillator(a=a,alpha=alpha,beta=beta)      
  elif mt == 'SJ3D': mod = models.ReducedSetHindmarshRose(K11=k11,K12=k12,K21=k21)
    
  conn   = connectivity.Connectivity(load_default=True, speed=spd)
  cpl    = coupling.Linear(a=gca)
  solver = integrators.HeunStochastic(dt=dt, noise=noise.Additive(nsig=array([nsig])))
  mons   = (monitors.TemporalAverage(period=1.0)) # 200 Hz

  sim = simulator.Simulator(model=mod,connectivity=conn,coupling=cpl,
                            integrator=solver, monitors=mons)
  sim.configure()
    
  ts,ys = sim.run(simulation_length=sim_len)[0]
    
  if mt == 'SJ3D': ys = np.squeeze(ys.sum(axis=3))        
  elif mt == 'G2DO': ys = np.squeeze(ys)
            
  toc = time() - tic
        
  if returntime: 
    return ts,ys,toc
  else: 
    return ts,ys


# Function for distributing runs across multiple cores (very handy!)

def multiprocfunc(func,args):
  pool = multiprocessing.Pool()
  tic = time()
  results = pool.map(func, args)
  toc = time() - tic
  return results,toc


# Function to compute FC from the outputs of the `run_sim`

def compFC(ys,ts):
  t_mask = ts > 150.0    
  return corrcoef(ys[t_mask].T)

# Function to compare outputs of `compSC` to structural connectivity

def corrSCFC(FC):
  return corrcoef(FC.ravel(), SC.ravel())[0, 1]


# Another that flips the above function (for use where we need to 
# minimize rather than maximize the objective function)

def oneminuscorrSCFC(FC):
  return 1-corrSCFC(FC)

# Wrapping the objective function (SC-FC correlation) 
# and simulation function up together:

def corrSCFC_4params(params):
  ts,ys,toc = run_sim(**params)
  fc = compFC(ys,ts)
  r = corrSCFC(fc)
  return r

def oneminuscorrSCFC_4params(params):
  ts,ys,toc = run_sim(**params)
  fc = compFC(ys,ts)
  r = oneminuscorrSCFC(fc)
  return r



# Plotting function; also displays SC-FC correlation for good measure

def plot_FCSC(FC): 
  
  fig, ax = plt.subplots(ncols=2, figsize=(12,3))
  sns.heatmap(FC-np.eye(FC.shape[0]), xticklabels='', 
              yticklabels='', ax=ax[0], 
              cmap='coolwarm') #, vmin=-1, vmax=1)
    
  sns.heatmap(SC/SC.max(), xticklabels='', yticklabels='', 
              ax=ax[1], cmap='coolwarm', vmin=-1,vmax=1)
    
  r = corrSCFC(FC)
    
  ax[0].set_title('simulated FC. \n(SC-FC r = %1.4s )'  %r)
  ax[1].set_title('SC')
    

Also at this point we will define the initial guess values and search ranges for all parameters that will be of interest

gca_init,  gca_bounds  = 0.004,(0.001,0.5)  # global coupling 'a' param 
nsig_init, nsig_bounds = 0.01, (0.0001,0.1) # noise stdev
spd_init,  spd_bounds  = 3.,   (1, 15)      # conduction velocity
a_init,    a_bounds    = 0,    (-5,5)       # G2DO 'a' parameter
k11_init,  k11_bounds  = 0.5,  (0.05,0.95)  # SJ3D K11
k12_init,  k12_bounds  = 0.5,  (0.05,0.95)  # SJ3D K12
k21_init,  k21_bounds  = 0.5,  (0.05,0.95)  # SJ3D K21

And the set the default connectivity weights matrix as a global variable, which we shall be using extensively

SC = connectivity.Connectivity(load_default=True).weights

Some examples of how the above functions are to be called, with some non-default arguments:

# <!-- collapse=False-->

f = '%s/pics/g2do_eg1_bestFC.png' %outdir
# plot SC vs. FC
ts_eg1,ys_eg1,time_eg1 = run_sim(mt='G2DO') # Run G2DO sim 
FC_eg1 = compFC(ys_eg1,ts_eg1)              # Compute FC
plot_FCSC(FC_eg1)                           # Compare to SC and plot

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

# upload to cloud and embed link in notebook 
cap = ''; label = 'Noise vs. coupling strength PSE - FC vs. SC`; G2DO example'
im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500))
d(im)
broken link Noise vs. coupling strength PSE - FC vs. SC`; G2DO example.
# <!-- collapse=False-->

f = '%s/pics/sj3d_eg2_bestFC.png' %outdir
# plot SC vs. FC
ts_eg2,ys_eg2,time_eg2 = run_sim(mt='SJ3D', dt=0.01) # Run SJ3D sim
FC_eg2 = compFC(ys_eg2,ts_eg2)                       # Compute FC
plot_FCSC(FC_eg2)                                    # Compare to SC and plot

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

# upload to cloud and embed link in notebook 
cap = ''; label = 'Noise vs. coupling strength PSE - FC vs. SC; SJ3D example'
im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500))
d(im)
broken link Noise vs. coupling strength PSE - FC vs. SC; SJ3D example.

Run through a couple of different configurations of model type, run time to get an idea of how long each model fit iteration will take, and also the effects of run times on model fits with fixed parameters.

(FYI: Unlike G2DO, SJ3D seems to need an integration step size of ~0.01 (otherwise you just get nans); so they can't really be compared directly, but certainly SJ3D sims will take always take longer than G2DO sims.)

testparams = {'g_500':   {'mt': 'G2DO', 'sim_len': 500,  'dt': 0.1},
              'g_1000':  {'mt': 'G2DO', 'sim_len': 1000, 'dt': 0.1},
              'g_10000': {'mt': 'G2DO', 'sim_len': 10000, 'dt': 0.1},
              's_500':   {'mt': 'SJ3D', 'sim_len': 500,  'dt': 0.01},
              's_1000':  {'mt': 'SJ3D', 'sim_len': 1000, 'dt': 0.01},
              's_10000': {'mt': 'SJ3D', 'sim_len': 10000, 'dt': 0.01},
             }
testsimres = {k: run_sim(**v) for k,v in testparams.items()}

Collating the results...

ks = testparams.keys()

testsimtimes = {k: v[2]/60 for k,v in testsimres.items()}
testsimFCs = {k: compFC(v[1],v[0]) for k,v in testsimres.items()}
testsimcorrs = {k: corrSCFC(v) for k,v in testsimFCs.items()}

df_rt = pd.DataFrame(testsimtimes.values(), index=ks, columns=['runtime (mins)'])
df_fc = pd.DataFrame(testsimcorrs.values(), index=ks,columns=['SCFC_corr'])
df_params = pd.DataFrame(testparams).T
dfs = pd.concat([df_params,df_fc,df_rt], axis=1)

dfs.sort(columns=['dt', 'sim_len'], inplace=True)

d(dfs)
dt mt sim_len SCFC_corr runtime (mins)
s_500 0.01 SJ3D 500 0.050131 4.027150
s_1000 0.01 SJ3D 1000 0.072436 6.898374
s_10000 0.01 SJ3D 10000 0.177292 76.380730
g_500 0.1 G2DO 500 0.183572 0.095402
g_1000 0.1 G2DO 1000 0.240530 0.192491
g_10000 0.1 G2DO 10000 0.237397 1.928348

Benchmark approaches

The main point of this investigation is to see if intelligent optimization algorithms can get speed and quality improvements over two 'benchmark' approaches.

The first is a full parameter parameter space exploration (PSE) over a pre-defined grid. Here I am again borrowing directly from the Duke's example notebook, who gave a very nice and succinct demonstration of an SC-FC matching for two 2-dimensional parameter sweeps. That example also showed how to distribute such a simulation set over multiple cores, and we will make good use of that here too.

The second 'benchmark' approach is an iterative scheme described by Falcon et al. (2015), and also used in one form or another by various papers by Deco, Jirsa, Ritter, and co. Combinatorial explosion makes a full PSE impractical over any more than two or three parameter sets, and Falcon et al. describe an iterative procedure that progressively optimizes and fixes several parameters, in order of their influence on the resultant model behaviour.

In a sense these two benchmarks are two ends of spectrum; the first is too comprehensive, the second is in danger of being too expedient, and potentially susceptible to local minima. Still, they both work reasonable well and are good starting points.

PSEs

In the Duke's example notebook he looks at two parameter sweeps with the generic 2D oscillator model. The first is global coupling vs. noise, the second is global coupling vs. speed.

The following class wraps and generalizes a bit the code from those examples

class pse_fit_2D(object):
    
    
  def __init__(self,param_ranges,model_type, fixed_params={},sim_len=500):
    
    # params is a dict of parameter ranges (must have only 2 entries)
    # assumes that there are just two params
    
    # Initialize everything
    
    otherargs = copy(fixed_params)
    otherargs['sim_len'] = sim_len
    if model_type == 'G2DO': 
        otherargs.update({'mt': 'G2DO', 'dt': 0.1})
    elif model_type == 'SJ3D': 
        otherargs.update({'mt': 'SJ3D', 'dt': 0.01})
        
    ks = param_ranges.keys()
    combs = list(itertools.product(param_ranges[ks[0]],
                                   param_ranges[ks[1]]))
    param_combs = []
    for c_it,c in enumerate(combs):
      pcs = copy(otherargs)
      for k_it,k in enumerate(ks): pcs.update({k:c[k_it]})
      param_combs.append(pcs)
    
    self.param_combs =param_combs            
    
    # (are these needed?)
    self.param_ranges = param_ranges
    self.fixed_params = fixed_params
    self.model_type = model_type
    self.sim_len = sim_len
    
  def run(self):
         
    runtic = time()
    
    self.res,_ = multiprocfunc(corrSCFC_4params, self.param_combs)
    
    idx = argmax(self.res)
    
    self.finalfit = self.res[idx]
    self.finalvals = self.param_combs[idx]

    
    # Re-run the final parameter set to get the time series and FC matrix
    
    ts,ys,_ = run_sim(**self.param_combs[idx])
    FC = compFC(ys,ts)
    self.finalres = {'ts': ts, 'ys': ys, 'FC': FC}
    
    runtoc = time() - runtic
    
    self.run_time = runtoc/60
    self.completed = True
            
     
  def plot_pse(self, pr1k,pr2k,log10x=False,log10y=False):

    pr1vs = self.param_ranges[pr1k]
    pr2vs = self.param_ranges[pr2k]
    
    resarr = array(self.res).reshape((len(pr1vs),len(pr2vs)))
 
    X,Y = meshgrid(pr1vs,pr2vs)

    xlab,ylab =  '%s' %pr1k, '%s' %pr2k
    if log10x: 
      X = log10(X)
      xlab = '$log_{10}$ ' + xlab
    if log10y: 
      Y = log10(Y)
      ylab = '$log_{10}$ ' + ylab
            
    figure(figsize=(10, 5))
    contourf(X, Y, resarr, 15, cmap='jet')
    xlabel(xlab); ylabel(ylab)
    title('FC-SC correlation')
    colorbar();
        
        
        
  def _repr_html_(self):

    
    if self.completed: 
      msg = '<div>TVB PSE Model Fit: </div>'
      msg += '<div>completed in: %2.3f mins </div>' %self.run_time
      msg += '<div>best fit:     %s </div>' %self.finalfit
      msg += '<div>best params:     %s </div>' %self.finalvals
        
    else: 
      msg =  'not yet run'
        
    #return '<td>%s</td><td></td>' % (msg)    
    return msg    
        

This is how they went:

Parameter sweep for global coupling vs. noise

# lists of all parameters' values
gcas  =  10**r_[-2.5 : -1.0 : 50j]
nsigs =  10**r_[-5.0 : -1.0 : 50j]

param_ranges = {'gca': gcas, 'nsig': nsigs}
ps1 = pse_fit_2D(param_ranges, 'G2DO', sim_len=500)
ps1.run()
d(ps1)
TVB PSE Model Fit:
completed in: 87.962 mins
best fit: 0.279623120001
best params: {'gca': 0.070297321153254724, 'mt': 'G2DO', 'sim_len': 500, 'dt': 0.1, 'nsig': 0.00042919342601287783}
# <!-- collapse=False-->

f = '%s/pics/pse_nsig_vs_gca_pseplot.png' %outdir
# pse plot
ps1.plot_pse('nsig', 'gca', log10x=True,log10y=True)
plt.savefig(f, bbox_inches='tight'); plt.close()

# upload to cloud and embed link in notebook 
cap = ''; label = 'Noise vs. coupling strength PSE'
im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500))
d(im)
broken link Noise vs. coupling strength PSE.
# <!-- collapse=False-->

f = '%s/pics/pse_nsig_vs_gca_bestFC.png' %outdir
# pse plot
plot_FCSC(ps1.finalres['FC'])
plt.savefig(f, bbox_inches='tight'); plt.close()

# upload to cloud and embed link in notebook 
cap = ''; label = 'Noise vs. coupling strength PSE - FC vs. SC'
im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500))
d(im)
broken link Noise vs. coupling strength PSE - FC vs. SC.

Parameter sweep for coupling vs. speed

# lists of all parameters' values
gcas  =  10**r_[-2.5 : -1.0 : 50j]
spds =  r_[1.0 : 10 : 50j]

param_ranges = {'gca': gcas, 'spd': spds}

fixed_params = {'nsig': 0.0004291}

ps2 = pse_fit_2D(param_ranges, 'G2DO', sim_len=3000, 
                  fixed_params=fixed_params)
                  
ps2.run()
d(ps2)
TVB PSE Model Fit:
completed in: 394.147 mins
best fit: 0.66453012485
best params: {'gca': 0.032374575428176434, 'mt': 'G2DO', 'sim_len': 3000, 'nsig': 0.0004291, 'dt': 0.1, 'spd': 6.6938775510204085}
# <!-- collapse=False-->

f = '%s/pics/pse_spd_vs_gca_pseplot.png' %outdir
# pse plot
ps2.plot_pse('spd', 'gca',log10y=True)

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

# upload to cloud and embed link in notebook 
cap = ''; label = 'Coupling strength vs. speed PSE'
im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500))
d(im)
broken link Coupling strength vs. speed PSE.
# <!-- collapse=False-->

f = '%s/pics/pse_gca_vs_spd_bestFC.png' %outdir
# pse plot
plot_FCSC(ps2.finalres['FC'])
plt.savefig(f, bbox_inches='tight'); plt.close()

# upload to cloud and embed link in notebook 
cap = ''; label = 'coupling strength vs. conduction velocity - FC vs. SC'

im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500))
d(im)
broken link coupling strength vs. conduction velocity - FC vs. SC.

Comments on the two PSEs:

  • main difference between FC matrices from pse 1 and 2 is the XH connections; they are lower in pse 2 with better fits
  • this is probably due to weaker coupling
  • essentially the increased conduction speed compensates for the decreased coupling strength in the optimal models
  • (bear in mind this matrix may not be best for these Qs; XH conns are very sparse...)
  • the optimal speed value is still in the main physiological range of CVs

Iterative method

Falcon et al. (2015) arrive at the 'best' (in their case maximize defined as maximal variance) models for different individuals by progressively optimizing parameters of the SJ3D model in the following sequence: G, conduction velocity, K12, K21, K11.

The following class implements this sequential approach in a flexible manner, allowing arbitrary choices for parameters of interest and their order:

class stepwise_fit(object):
    
  # Params is a dictionary of lists
  # Order is an ordered list of keys
  # initvals is a dict of initial values
    
  def __init__(self,param_ranges,step_order,initvals,model_type, sim_len=500):
        
    # Initialize everything
    
    
    # Make some empty attributes to fill
    
    for v in ['param_ranges', 'step_order', 'initvals', 'model_type']:
      exec 'self.%s = %s' %(v,v)
    
    self.finalvals = {}
    self.step_history = []
    self.completed = False
    
    
    # Add info to the corrSCFC_4params initialization arguments 
    # regarding model type and integration step size
   
    self.initargs = copy(initvals)
    self.initargs['sim_len'] = sim_len
    if self.model_type == 'G2DO': 
        self.initargs.update({'mt': 'G2DO', 'dt': 0.1})
    elif self.model_type == 'SJ3D': 
        self.initargs.update({'mt': 'SJ3D', 'dt': 0.01})
    
    
  def run(self):
        
    # Iterate through the sequence prescribed in self.order. 
    # At each stage, issue a set of simulation runs in parallel, 
    # varying a single parameter. The varied parameter values 
    # are defined by the ranges set in self.param_ranges. The 
    # fixed parameters are set to be either the initialization 
    # values defined in self.initvals, or, the final values in 
    # self.finalvals arrived at in previous steps. 
        
    runtic = time()
    
    
    for step in self.step_order:
      
      # Define parameter set to run for this step
      params = self.set_step_params(step)
        
      # Run sims
      res,toc = multiprocfunc(corrSCFC_4params, params)
    
      # Identify best fit and log
      idx = np.nanargmax(res)
      self.finalvals[step] = params[idx][step]
        
      # Keep some other info
      self.step_history.append({'varying': step, 'params': params, 
                                'res': res, 'toc': toc})
        
      
    self.finalfit = res[idx]
    
    
    # Re-run the final parameter set to get the time series and FC matrix
    
    ts,ys,_ = run_sim(**params[idx])
    FC = compFC(ys,ts)
    self.finalres = {'ts': ts, 'ys': ys, 'FC': FC}
    
    runtoc = time() - runtic
    self.run_time = runtoc/60
    
    self.completed = True
    
    
  def set_step_params(self,step):
            
    # Specify a list of arguments for corrSCFC_4params for the current step
    
    params = copy(self.initargs)   # Start with the initial values
    params.update(self.finalvals)  # Add in final values obtained so far
    
    paramslist = []
    for varval in self.param_ranges[step]:
      ps = copy(params)
      ps[step] = varval
      paramslist.append(ps)
      
    return paramslist

    
  def _repr_html_(self):

    
    if self.completed: 
      msg = '<div>TVB Stepwise Model Fit: </div>'
      msg += '<div>completed in: %2.3f mins </div>' %self.run_time
      msg += '<div>best fit:     %s </div>' %self.finalfit
      msg += '<div>best params:     %s </div>' %self.finalvals
        
    else: 
      msg =  'not yet run'
        
    #return '<td>%s</td><td></td>' % (msg)    
    return msg    
        
    
  def print_summary(self):
        
    if self.completed: 
      print 'completed in %.5f minutes' %self.run_time
      print 'best fit: %s' %self.finalfit
        
    else:
      print 'not yet run'
        
        

Ok, let's put this class to work:

G2dO sims

a_bounds = [-5.0,5.0]
alpha_bounds = [-5.0,5.0]
beta_bounds = [-5.0,5.0]


a_init = -2.0
alpha_init = 1.0
beta_init = 1.0
_as = np.linspace(a_bounds[0], a_bounds[1],11)
_alphas = np.linspace(alpha_bounds[0],alpha_bounds[1],11)
_betas = np.linspace(beta_bounds[0],beta_bounds[1],11)
gca = 0.03237
nsig = 0.0004291
spd = 6.69387

param_ranges = {'a': _as, 'alpha': _alphas, 'beta': _betas}
order = ['a', 'alpha', 'beta']
inits = {'a': a_init, 'alpha': alpha_init, 'beta': beta_init,
         'gca': gca, 'nsig': nsig, 'spd': spd}
sf = stepwise_fit(param_ranges,order,inits,'G@DO',sim_len=3000)
sf.run()
d(sf)
TVB Stepwise Model Fit:
completed in: 6.290 mins
best fit: 0.664576134684
best params: {'a': 0.0, 'alpha': 1.0, 'beta': 1.0}

Ok so that returned the default params as oprtimal; so no improvement there...

Try changing order?

sf = stepwise_fit(param_ranges,['a', 'beta', 'alpha'],
                  inits,'G@DO',sim_len=3000)
sf.run()
d(sf)
TVB Stepwise Model Fit:
completed in: 6.180 mins
best fit: 0.664576134684
best params: {'a': 0.0, 'alpha': 1.0, 'beta': 1.0}

Different question: how well does the stepwise method get to the pse optimum?

# this should be equiv of pse 2; i.e. nsig is fixed

param_ranges = {'gca': gcas, 'spd': spds }
order = ['gca', 'spd']
inits = {'nsig': 0.00042919, 'gca': gca_init, 'spd': spd_init}

sf = stepwise_fit(param_ranges,order,
                  inits,'G@DO',sim_len=3000)
sf.run()
d(sf)
TVB Stepwise Model Fit:
completed in: 17.192 mins
best fit: 0.511790953542
best params: {'gca': 0.056898660290182958, 'spd': 3.2040816326530615}

Nope. That was worse.

SJ3D sims

gcas = np.linspace(0.001, 0.9, 10)
spds = np.linspace(spd_bounds[0], spd_bounds[1],10)
k12s = np.linspace(k12_bounds[0], k12_bounds[1],10)
k21s = np.linspace(k21_bounds[0], k21_bounds[1],10)
k11s = np.linspace(k11_bounds[0], k11_bounds[1],10)

param_ranges = {'gca': gcas, 'spd': spds, 'k12': k12s, 'k21': k21s,'k11': k11s}
order = ['gca', 'spd', 'k12', 'k21', 'k11']
inits = {'gca':gca,'spd':spd_init,
         'k12':k12_init,'k21':k21_init,'k11':k11_init,
         'nsig': 0.0004}
# Run fit
sf = stepwise_fit(param_ranges,order,inits,'SJ3D',sim_len=1000)
sf.run()
d(sf)
TVB Stepwise Model Fit:
completed in: 103.448 mins
best fit: 0.214318890273
best params: {'gca': 0.10088888888888889, 'k12': 0.94999999999999996, 'k11': 0.14999999999999999, 'k21': 0.25, 'spd': 11.888888888888889}

urgh.

Higher noise level...

inits = {'gca':gca,'spd':spd_init,
         'k12':k12_init,'k21':k21_init,'k11':k11_init,
         'nsig': 0.001}

# Run fit
sf = stepwise_fit(param_ranges,order,inits,'SJ3D',sim_len=1000)
sf.run()
d(sf)
TVB Stepwise Model Fit:
completed in: 102.193 mins
best fit: 0.222072471038
best params: {'gca': 0.10088888888888889, 'k12': 0.94999999999999996, 'k11': 0.84999999999999998, 'k21': 0.050000000000000003, 'spd': 11.888888888888889}
# <!-- collapse=False-->

f = '%s/pics/sj3d_stepwise_v1_bestFC.png' %outdir
plot_FCSC(sf.finalres['FC'])

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

# upload to cloud and embed link in notebook 
cap = ''; label = 'SJ3D stepwise fitting v1 - FC vs. SC`'
im = nb_fig(f,label,cap,cnb,show_fignum=False,size=(800,500))
d(im)
broken link SJ3D stepwise fitting v1 - FC vs. SC`.

No better...

I still feel like we should be able to get the model fit up higher.

Scipy optimization algorithms

The above PSEs searched through thousands of parameter values and arrived at the following best fits:

fc_sc_for_params:

  • best corr: 0.27
  • best params: cs = 0.070297321153254724, D = 0.00042919342601287783

fc_sc_for_params2:

  • best corr: 0.66452569100250425
  • best params: cs = 0.032374575428176434, cv = 6.6938775510204085

Let's see if we get the same fits using an optimization algorithm.

scipy.optimize.minimize

Modified from here

class scipy_optimize_fit(object):
    
  # Params is a dictionary of lists
  # Order is an ordered list of keys
  # initvals is a dict of initial values
    
  def __init__(self,boundvals,initvals,model_type,
               fixed_params={},sim_len=500, opt_method = 'L-BFGS-B'):
        
    # Initialize everything
    
    
    # Make some empty attributes to fill
    
    for v in ['boundvals', 'initvals', 'model_type', 'opt_method']:
      exec 'self.%s = %s' %(v,v)
    
    self.finalvals = {}
    self.completed = False

    # Add info to the corrSCFC_4params initialization arguments 
    # regarding model type and integration step size
   
    self.initargs = copy(initvals)
    self.initargs['sim_len'] = sim_len
    if self.model_type == 'G2DO': 
        self.initargs.update({'mt': 'G2DO', 'dt': 0.1})
    elif self.model_type == 'SJ3D': 
        self.initargs.update({'mt': 'SJ3D', 'dt': 0.01})
    
    
  def fitfunc(self,args):
        
    args = update(args, self.initargs)
    return oneminuscorrSCFC_4params(**args)

    params = copy(self.initargs)   # Start with the initial values
    params.update(self.finalvals)  # Add in final values obtained so far
    
    paramslist = []
    for varval in self.param_ranges[step]:
      ps = copy(params)
      ps[step] = varval
      paramslist.append(ps)
      
    return paramslist



  def run(self):
        
    runtic = time()
    
    res = scipy.optimize.minimize(self.fitfunc,self.initvals,
                                  bounds=self.boundvals,method=self.opt_method)
    

    runtoc = time() - tic


    
    for step in self.step_order:
      
      # Define parameter set to run for this step
      params = self.set_step_params(step)
        
      # Run sims
      res,toc = multiprocfunc(corrSCFC_4params, params)
    
      # Identify best fit and log
      idx = np.nanargmax(res)
      self.finalvals[step] = params[idx][step]
        
      # Keep some other info
      self.step_history.append({'varying': step, 'params': params, 
                                'res': res, 'toc': toc})
        
      
    self.finalfit = res[idx]
    
    
    # Re-run the final parameter set to get the time series and FC matrix
    
    ts,ys,_ = run_sim(**params[idx])
    FC = compFC(ys,ts)
    self.finalres = {'ts': ts, 'ys': ys, 'FC': FC}
    
    runtoc = time() - runtic
    self.run_time = runtoc/60
    
    self.completed = True
    
    
  def set_step_params(self,step):
            
    # Specify a list of arguments for corrSCFC_4params for the current step
    
    params = copy(self.initargs)   # Start with the initial values
    params.update(self.finalvals)  # Add in final values obtained so far
    
    paramslist = []
    for varval in self.param_ranges[step]:
      ps = copy(params)
      ps[step] = varval
      paramslist.append(ps)
      
    return paramslist

    
  def _repr_html_(self):

    
    if self.completed: 
      msg = '<div>TVB Stepwise Model Fit: </div>'
      msg += '<div>completed in: %2.3f mins </div>' %self.run_time
      msg += '<div>best fit:     %s </div>' %self.finalfit
      msg += '<div>best params:     %s </div>' %self.finalvals
        
    else: 
      msg =  'not yet run'
        
    #return '<td>%s</td><td></td>' % (msg)    
    return msg    
        
    
  def print_summary(self):
        
    if self.completed: 
      print 'completed in %.5f minutes' %self.run_time
      print 'best fit: %s' %self.finalfit
        
    else:
      print 'not yet run'
        
        

set inits to the best values...

nsig = 0.0004291
gca_init = 0.032374575428176434
sim_len = 3000
nsig = 0.0004291
dt = 0.1
spd_init = 6.6938775510204085
a_init = 0
def fitfunc(args):
  gca,spd,a = args  
  return oneminuscorrSCFC_4params({'gca': gca, 'spd': spd, 'a': a,
                                   'nsig': nsig, 'sim_len': sim_len})

inits = [gca_init,spd_init,a_init]
bounds = [gca_bounds,spd_bounds,a_bounds]

tic = time()
res = scipy.optimize.minimize(fitfunc,inits,bounds=bounds,method='L-BFGS-B')
toc = time() - tic
res
Out[559]:
   status: 0
  success: True
     nfev: 84
 hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
      fun: 0.33313886669076398
        x: array([ 0.0307907 ,  6.69387755,  0.10224276])
  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
      jac: array([ -1.90292226e-05,   0.00000000e+00,  -2.33146835e-07])
      nit: 18
1-res['fun']
Out[561]:
0.66686113330923602
toc/60
Out[562]:
62.16640124718348
res['x']
Out[563]:
array([ 0.0307907 ,  6.69387755,  0.10224276])

No improvement there then...

scipy.optimize.brute

http://scipy.github.io/devdocs/generated/scipy.optimize.basinhopping.html#scipy.optimize.basinhopping

>>> rranges = (slice(-4, 4, 0.25), slice(-4, 4, 0.25))
>>> from scipy import optimize
>>> resbrute = optimize.brute(f, rranges, args=params, full_output=True,
...                           finish=optimize.fmin)
>>> resbrute[0]  # global minimum
array([-1.05665192,  1.80834843])
>>> resbrute[1]  # function value at global minimum
-3.4085818767

scipy.optimize.basinhopping

nsig = 0.0004291
gca_init = 0.032374575428176434
sim_len = 3000
nsig = 0.0004291
dt = 0.1
spd_init = 6.6938775510204085
a_init = 0
def fitfunc(args):
  gca,spd,a = args  
  return oneminuscorrSCFC_4params({'gca': gca, 'spd': spd, 'a': a,
                                   'nsig': nsig, 'sim_len': sim_len})

inits = [gca_init,spd_init,a_init]
bounds = [gca_bounds,spd_bounds,a_bounds]

tic = time()
res = scipy.optimize.basinhopping(fitfunc,inits,minimizer_kwargs = {'bounds': bounds})

toc = time() - tic

Playdoh library

# to do...

Conclusions

# to do...

misc TVB