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)
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)
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)
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)
# <!-- 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)
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)
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)
# <!-- 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)
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)
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)
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)
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)
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)
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)
# <!-- 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)
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
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']
0.66686113330923602
toc/60
62.16640124718348
res['x']
array([ 0.0307907 , 6.69387755, 0.10224276])
No improvement there then...
scipy.optimize.brute¶
>>> 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...