# sim_tools.py - Routines for running a model simulation
#
# Author: Stefan Fuertinger [stefan.fuertinger@esi-frankfurt.de]
# Created: June 23 2014
# Last modified: <2017-10-19 16:40:42>
from __future__ import division
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import ipdb; import sys
import imp
import h5py
import os
import re
import psutil
from datetime import datetime
from texttable import Texttable
from scipy import ndimage
from nipy.modalities.fmri import hrf, utils
import shutil
import tempfile
from recipes import moveit
try:
from the_model import par, solve_model
except:
print "\n\tWARNING: Could not import the model - `run_model` will not work!"
print "\tTry running `make all` in a terminal first"
##########################################################################################
[docs]def run_model(V0, Z0, DA0, task, outfile, \
seed=None, paramfile='parameters.py', symsyn=True, verbose=True, ram_use=0.2,\
**kwargs):
"""
Run a simulation using the neural population model
Parameters
----------
V0 : NumPy 1darray
Initial conditions for excitatory neurons. If `N` regions are simulated then `V0` has
to have length `N`.
Z0 : NumPy 1darray
Initial conditions for inhibitory neurons. If `N` regions are simulated then `Z0` has
to have length `N`.
DA0 : NumPy 1darray
Initial conditions for dopamine levels. If `N` regions are simulated then `DA0` has
to have length `N`.
task : str
Specify which task should be simulated. Currently, only 'rest' and 'speech' are supported.
outfile : str
File-name (including path if not in working directory) of HDF5 container that will be created to
save simulation results. See Notes for the structure of the generated container. Any existing
file will be renamed. The user has to have writing permissions for the given location.
seed : int
Random number generator seed. To make meaningful comparisons between successive simulation
runs, the random number seed should be fixed so that the solver uses the same Wiener process
realizations. Also, if synaptic coupling strengths are sampled from a probability distribution,
simulation results will vary from run to run unless the seed is fixed.
paramfile : str
Parameter file-name (including path if not in working directory) that should be used for simulation.
The parameter file has to be a Python file (.py extension). For more details refer to `Examples`
below. You should have received a sample parameter file (`parameters.py`) with this copy
of `sim_tools.py`.
symsyn : bool
Boolean switch determining whether synaptic coupling strengths should be symmetrized between
hemispheres (`symsyn=True`) or not.
verbose : bool
If `True` the code will print a summary of the most important parameters and all used
keyword arguments (see below) in the simulation together with a progress bar to (roughly)
estimate run time (requires the `progressbar` module).
ram_use : float
Fraction of memory to use for caching simulation results before writing to disk
(0 < `ram_use` < 1). More available memory means fewer disk-writes and thus better performance,
i.e, the larger `ram_use` the faster this code runs. However, if too much RAM is allocated by
this routine it may stall the executing computer. By default, `ram_use = 0.2`, i.e., around 20%
of available memory is used.
kwargs : additional keyword arguments
Instead of relying solely on a static file to define parameter values, it is also possible
to pass on parameters to the code using keyword arguments (see `Examples` below). Note: parameters
given as keyword arguments have higher priority than values set in `paramfile`, i.e., if
`p1 = 1` is defined in `paramfile` but `p1 = 2` is a keyword argument, the code will use
`p1 = 2` in the simulation. This behavior was intentionally implemented to enable the use
of this function within a parameter identification framework.
Returns
-------
Nothing : None
Simulation results are saved in the HDF5 container specified by `outfile`. See `Notes` for details.
Notes
-----
Due to the (usually) high temporal resolution of simulations, results are not kept in memory (and
thus returned as variable in the caller's work-space) but saved directly to disk using the HDF5
container `outfile`. The code uses the HDF library's data chunking feature to save entire
segments on disk while running. By default the code will allocate around 20% of available
memory to cache simulation results. Hence, more memory leads to fewer disk-writes during
run-time and thus faster performance.
The structure of the generated output container is as follows: all state variables and the dopaminergic
gain `Beta` are stored at the top-level of the file. Additionally, the employed coupling
matrix `C` and dopamine connection matrix `D` are also saved in the top level group. All used parameters
are saved in the subgroup `params`.
Examples
--------
Let `V0`, `Z0`, and `DA0` (NumPy 1darrays of length `N`) be initial conditions of the
model. Assuming that a valid parameter file (called `parameters.py`) is located in the current
working directory, the following call will run a resting state simulation and save the output
in the HDF5 container `sim_rest.h5`
>>> run_model(V0,Z0,DA0,'rest','sim_rest.h5')
Assume another parameter file, say, `par_patho.py` hold parameter settings simulating a
certain pathology. Then the command
>>> run_model(V0,Z0,DA0,'rest','patho/sim_rest_patho.h5',paramfile='par_patho.py')
runs a resting state simulation with the same initial conditions and saves the result in
the container `sim_rest_patho.h5` in the sub-directory `patho` (which must already exist, otherwise
an error is raised).
If only one or two parameters should be changed from their values found in a given parameter file,
it is probably more handy to change the value of these parameters from the command line, rather
than to write a separate parameter file (that is identical to the original one except for two
values). Thus, assume the values of `VK` and `VL` should be -0.4 and -0.9 respectively, i.e.,
different than those found in (the otherwise fine) `par_patho.py`. Then the command
>>> run_model(V0,Z0,DA0,'rest','patho/sim_rest_patho.h5',paramfile='par_patho.py',VK=-0.4,VL=-0.9)
runs the same resting state simulation as above but with `VK=-0.4` and `VL=-0.9`. This feature
can also be used to efficiently embed `run_model` in a parameter identification framework.
See also
--------
plot_sim : plot simulations generated by run_model
References
----------
.. [1] S. Fuertinger, J. C. Zinn, and K. Simonyan. A Neural Population Model Incorporating
Dopaminergic Neurotransmission during Complex Voluntary Behaviors. PLoS Computational Biology,
10(11), 2014.
"""
# Sanity checks for initial conditions
n = np.zeros((3,))
vnames = ['V0','Z0','DA0']
for i, vec in enumerate([V0, Z0, DA0]):
arrcheck(vec,'vector',vnames[i])
n[i] = vec.size
if np.unique(n).size > 1:
raise ValueError('The initial conditions for `V`, `Z`, and `DA` have to have the same length!')
# Check the `task` string
if not isinstance(task,(str,unicode)):
raise TypeError('Task has to be specified as string, not '+type(task).__name__+'!')
task = str(task)
if task != 'rest' and task != 'speech':
raise ValueError("The value of `task` has to be either 'rest' or 'speech'!")
# The path to the output file should be a valid
if not isinstance(outfile,(str,unicode)):
raise TypeError('Output filename has to be a string!')
outfile = str(outfile)
if outfile.find("~") == 0:
outfile = os.path.expanduser('~') + outfile[1:]
slash = outfile.rfind(os.sep)
if slash >= 0 and not os.path.isdir(outfile[:outfile.rfind(os.sep)]):
raise ValueError('Invalid path for output file: '+outfile+'!')
# Set or get random number generator seed
if seed is not None:
scalarcheck(seed,'seed',kind='int')
else:
seed = np.random.get_state()[1][0]
seed = int(seed)
# Make sure `paramfile` is a valid path
if not isinstance(paramfile,(str,unicode)):
raise TypeError('Parameter file has to be specified using a string!')
paramfile = str(paramfile)
if paramfile.find("~") == 0:
paramfile = os.path.expanduser('~') + paramfile[1:]
if not os.path.isfile(paramfile):
raise ValueError('Parameter file: '+paramfile+' does not exist!')
# Make sure `symsyn` and `verbose` are Boolean
if not isinstance(symsyn,bool):
raise TypeError("The switch `symsyn` has to be Boolean!")
if not isinstance(verbose,bool):
raise TypeError("The switch `verbose` has to be Boolean!")
# Finally, check `ram_use`
scalarcheck(ram_use,'ram_use',bounds=[0,1])
# Append '.h5' extension to `outfile` if necessary
if outfile[-3:] != '.h5':
outfile = outfile + '.h5'
# Check if `paramfile` has an extension, if yes, rip it off
if paramfile[-3:] == '.py':
paramfile = paramfile[0:-3]
# Divide `paramfile` into file-name and path
slash = paramfile.rfind(os.sep)
if slash < 0:
pth = '.'
fname = paramfile
else:
pth = paramfile[0:slash+1]
fname = paramfile[slash+1:]
# Import parameters module and initialize corresponding dictionary (remove `__file__`, etc)
param_py = imp.load_module(fname,*imp.find_module(fname,[pth]))
p_dict = {}
for key, value in param_py.__dict__.items():
if key[0:2] != "__":
p_dict[key] = value
# Try to load coupling and dopamine pathway matrices
mfile = "None"
vnames = ['C','D']
for mat_str in vnames:
if kwargs.has_key(mat_str):
p_dict[mat_str] = kwargs[mat_str]
else:
try:
p_dict[mat_str] = h5py.File(param_py.matrices,'r')[mat_str].value
mfile = p_dict['matrices']
except:
raise ValueError("Error reading `"+param_py.matrices+"`!")
arrcheck(p_dict[mat_str],'matrix',mat_str)
# Try to load ROI names
try:
names = h5py.File(param_py.matrices,'r')['names'].value
mfile = p_dict['matrices']
except:
try:
names = kwargs['names']
except:
raise ValueError("A NumPy 1darray or Python list of ROI names has to be either specified "+\
"in a matrix container or provided as keyword argument!")
p_dict['names'] = names
# See if we have an (optional) list/array of ROI-shorthand labels
try:
p_dict['labels'] = h5py.File(param_py.matrices,'r')['labels'].value
mfile = p_dict['matrices']
except:
if kwargs.has_key('labels'):
p_dict['labels'] = kwargs['labels']
# Put ones on the diagonal of the coupling matrix to ensure compatibility with the code
np.fill_diagonal(p_dict['C'],1.0)
# Get dimension of matrix and check correspondence
N = p_dict['C'].shape[0]
if N != p_dict['D'].shape[0]:
raise ValueError("Dopamine and coupling matrices don't have the same dimension!")
if len(names) != N:
raise ValueError("Matrix is "+str(N)+"-by-"+str(N)+" but `names` has length "+str(len(names))+"!")
for nm in names:
if not isinstance(nm,(str,unicode)):
raise ValueError("Names have to be provided as Python list/NumPy array of strings!")
# If user provided some additional parameters as keyword arguments, copy them to `p_dict`
for key, value in kwargs.items():
p_dict[key] = value
# Get synaptic couplings (and set seed of random number generator)
np.random.seed(seed)
if kwargs.has_key('aei'):
aei = kwargs['aei']
else:
aei = eval(param_py.aei)
if kwargs.has_key('aie'):
aie = kwargs['aie']
else:
aie = eval(param_py.aie)
if kwargs.has_key('ani'):
ani = kwargs['ani']
else:
ani = eval(param_py.ani)
if kwargs.has_key('ane'):
ane = kwargs['ane']
else:
ane = eval(param_py.ane)
# If wanted, make sure left/right hemispheres have balanced coupling strengths
if symsyn:
# Get indices of left-hemispheric regions and throw a warning if left/right don't match up
regex = re.compile("[Ll]_*")
match = np.vectorize(lambda x:bool(regex.match(x)))(names)
l_ind = np.where(match == True)[0]
r_ind = np.where(match == False)[0]
if l_ind.size != r_ind.size:
print "WARNING: Number of left-side regions = "+str(l_ind.size)+\
" not equal to number of right-side regions = "+str(r_ind.size)
# Equalize coupling strengths
aei[l_ind] = aei[r_ind]
aie[l_ind] = aie[r_ind]
ani[l_ind] = ani[r_ind]
ane[l_ind] = ane[r_ind]
# Save updated coupling strengths and random number generator seed in dictionary
p_dict['aei'] = aei
p_dict['aie'] = aie
p_dict['ani'] = ani
p_dict['ane'] = ane
p_dict['seed'] = seed
# If a resting state simulation is done, make sure dopamine doesn't kick in, i.e., enforce `rmax == rmin`
if task == 'rest':
rmax = np.ones((N,))*p_dict['rmin']
else:
if not kwargs.has_key('rmax'):
p_dict['rmax'] = eval(param_py.rmax)
# Save given task in dictionary
p_dict['task'] = task
# Get ion channel parameters
if not kwargs.has_key('TCa'):
p_dict['TCa'] = eval(param_py.TCa)
# Compute length for simulation and speech on-/offset times
len_cycle = p_dict['stimulus'] + p_dict['production'] + p_dict['acquisition']
speechon = p_dict['stimulus']
speechoff = p_dict['stimulus'] + p_dict['production']
# Save that stuff
p_dict['len_cycle'] = len_cycle
p_dict['speechon'] = speechon
p_dict['speechoff'] = speechoff
# Set/get initial time for simulation
if p_dict.has_key('tstart'):
tstart = p_dict['tstart'] # Use `p_dict` here, since `tstart` could be a kwarg!
if verbose: print "WARNING: Using custom initial time of "+str(tstart)+" (has to be in ms)!"
else:
tstart = 0
# Set/get step-size for simulation
if p_dict.has_key('dt'):
dt = p_dict['dt']
if verbose: print "WARNING: Using custom step-size of "+str(dt)+" (has to be in ms)!"
else:
dt = 1e-1
# Get sampling step size (in ms) and check if "original" step-size makes sense
ds = 1/p_dict['s_rate']*1000
if dt > ds:
print "WARNING: Step-size dt = "+str(dt)+\
" larger than chosen sampling frequency of "+str(s_rate)+"Hz."+\
" Using dt = "+str(ds)+"ms instead. "
dt = ds
# Compute sampling rate (w.r.t `dt`)
s_step = int(np.round(ds/dt))
# Save step-size and sampling rate in dictionary for later reference
p_dict['dt'] = dt
p_dict['s_step'] = s_step
# Compute end time for simulation (in ms) and allocate time-step array
tend = tstart + len_cycle*p_dict['n_cycles']*1000
tsteps = np.arange(tstart,tend,dt)
# Get the size of the time-array
tsize = tsteps.size
# Before laying out output HDF5 container, rename existing files to not accidentally overwrite 'em
moveit(outfile)
# Chunk outifle depending on available memory (eat up ~ 100*`ram_use`% of RAM)
datype = np.dtype('float64')
meminfo = psutil.virtual_memory()
maxmem = int(meminfo.available*ram_use/(5*N)/datype.itemsize)
maxmem += s_step - np.mod(maxmem,s_step)
# If the whole array fits into memory load it once, otherwise chunk it up
if tsteps.size <= maxmem:
blocksize = [tsize]
csize = int(np.ceil(tsize/s_step))
chunksize = [csize]
chunks = True
else:
bsize = int(tsize//maxmem)
rest = int(np.mod(tsize,maxmem))
blocksize = [maxmem]*bsize
if rest > 0: blocksize = blocksize + [rest]
numblocks = len(blocksize)
csize = int(np.ceil(maxmem/s_step))
restc = int(np.ceil(blocksize[-1]/s_step))
chunksize = [csize]*(numblocks - 1) + [restc]
chunks = (N,csize)
# Convert blocksize and chunksize to NumPy arrays
blocksize = np.array(blocksize)
chunksize = np.array(chunksize)
# Get the number of elements that will be actually saved
n_elems = chunksize.sum()
# Create output HDF5 container
f = h5py.File(outfile)
# Create datasets for numeric variables
f.create_dataset('C',data=p_dict['C'],dtype=datype)
f.create_dataset('D',data=p_dict['D'],dtype=datype)
f.create_dataset('V',shape=(N,n_elems),chunks=chunks,dtype=datype)
f.create_dataset('Z',shape=(N,n_elems),chunks=chunks,dtype=datype)
f.create_dataset('DA',shape=(N,n_elems),chunks=chunks,dtype=datype)
f.create_dataset('QV',shape=(N,n_elems),chunks=chunks,dtype=datype)
f.create_dataset('Beta',shape=(N,n_elems),chunks=chunks,dtype=datype)
f.create_dataset('t',data=np.linspace(tstart,tend,n_elems),dtype=datype)
# Save parameters (but exclude stuff imported in the parameter file)
pg = f.create_group('params')
for key,value in p_dict.items():
valuetype = type(value).__name__
if valuetype != 'instance' and valuetype != 'module' and valuetype != 'function':
pg.create_dataset(key,data=value)
# Close container and write to disk
f.close()
# Initialize parameter C-class (struct) for the model
params = par(p_dict)
# Concatenate initial conditions for the "calibration" run
VZD0 = np.hstack([V0.squeeze(),Z0.squeeze(),DA0.squeeze()])
# Set up parameters for an initial `len_init` (in ms) long resting state simulation to "calibrate" the model
len_init = 100
dt = 0.1
s_step = 10
rmax = np.zeros((N,))
tinit = np.arange(0,len_init,dt)
tsize = tinit.size
csize = int(np.ceil(tsize/s_step))
# Update `p_dict` (we don't use it anymore, so just overwrite stuff)
p_dict['dt'] = dt
p_dict['s_step'] = s_step
p_dict['rmax'] = rmax
parinit = par(p_dict)
# Create a temporary container for the simulation
tmpname = tempfile.mktemp() + '.h5'
tmpfile = h5py.File(tmpname)
tmpfile.create_dataset('V',shape=(N,csize),dtype=datype)
tmpfile.create_dataset('Z',shape=(N,csize),dtype=datype)
tmpfile.create_dataset('DA',shape=(N,csize),dtype=datype)
tmpfile.create_dataset('QV',shape=(N,csize),dtype=datype)
tmpfile.create_dataset('Beta',shape=(N,csize),dtype=datype)
tmpfile.flush()
# Run 100ms of resting state to get model to a "steady state" for the initial conditions
solve_model(VZD0,tinit,parinit,np.array([tsize]),np.array([csize]),seed,0,str(tmpfile.filename))
# Use final values of `V`, `Z` and `DA` as initial conditions for the "real" simulation
V0 = tmpfile['V'][:,-1]
Z0 = tmpfile['Z'][:,-1]
DA0 = tmpfile['DA'][:,-1]
VZD0 = np.hstack([V0.squeeze(),Z0.squeeze(),DA0.squeeze()])
# Close and delete the temporary container
tmpfile.close()
os.remove(tmpname)
# Let the user know what's going to happen...
pstr = "--"
if len(kwargs) > 0:
pstr = str(kwargs.keys())
pstr = pstr.replace("[","")
pstr = pstr.replace("]","")
pstr = pstr.replace("'","")
table = Texttable()
table.set_deco(Texttable.HEADER)
table.set_cols_align(["l", "l"])
table.add_rows([["Simulating ",task.upper()],\
["#cycles: ",str(p_dict['n_cycles'])],\
["parameter file:",paramfile+".py"],\
["keyword args:",pstr],\
["matrix file:",mfile],\
["output:",outfile]])
if verbose: print "\n"+table.draw()+"\n"
# Finally... run the actual simulation
solve_model(VZD0,tsteps,params,blocksize,chunksize,seed,int(verbose),outfile)
# Done!
if verbose: print "\nDone\n"
##########################################################################################
[docs]def plot_sim(fname,names="all",raw=True,bold=False,figname=None):
"""
Plot a simulation generated by `run_model`
Parameters
----------
fname : str
File-name (including path if not in working directory) of HDF5 container that was generated
by `run_model`.
names : str or Python list/NumPy 1darray
Specify regions to plot. Either use the region's name as found in the `params` group of
the HDF5 container given by `fname` (e.g., `names='L_IFG'`) or its index in the `names` list
(e.g., `names = 3`). Use a list or NumPy 1darray to specify more than one region
(e.g., `names = ['L_IFG','R_IFG']` or `names = [3,15]`). By default, all regions are plotted.
raw : bool
If `True` then the raw model output will be plotted. Depending on the setting of `names` (see
above) the simulation length, and the model dimension (i.e., the number of modeled regions)
this may result in a very 'busy' plot.
bold : bool
If True then the previously converted simulated BOLD signals will be plotted (if no BOLD
signal is found in the input container specified by `fname`, an error is raised).
figname : str
String to be used as window title for generated figures.
Returns
-------
Nothing : None
Notes
-----
None
See also
--------
run_model : used to run a simulation
make_bold : convert raw simulation output to a BOLD signal
"""
# Make sure `fname` is a valid file-path
f = checkhdf(fname,peek=True)
# Make sure `raw` and `bold` are either `True` or `False`
if not isinstance(raw,bool):
raise TypeError("The flag `raw` has to be Boolean!")
if not isinstance(bold,bool):
raise TypeError("The flag `bold` has to be Boolean!")
# Try to access given HDF5 container
try:
rois_infile = f['params']['names'].value.tolist()
f.close()
except:
raise ValueError("HDF5 file "+fname+" does not have the required fields!")
# Check if list of ROI-names to plot was provided. If yes, make sure they make sense
if not isinstance(names,(str,unicode)):
doleg = True
try:
names = list(names)
except:
raise TypeError("Regions to plot have to be provided as Python list or NumPy 1darray!")
idx = []
for name in names:
try:
idx.append(rois_infile.index(name))
except:
raise ValueError("Region "+name+"not found in file!")
else:
if names == "all":
idx = range(len(rois_infile))
doleg = False
else:
try:
idx = rois_infile.index(names)
except:
raise ValueError("Region "+names+"not found in file!")
doleg = True
# Check if `figname` is actually printable
if figname != None:
if not isinstance(figname,(str,unicode)):
raise TypeError("Figure name has to be a string!"+type(figname).__name__+"!")
figname = str(figname)
# Fix sorting of idx so that smart indexing below works
if doleg:
idx = np.array(idx)
names = np.array(names)
sorted = idx.argsort()
names = names[sorted]
idx = idx[sorted]
# After all the error checking, reopen the file
f = h5py.File(fname,'r')
# Turn on interactive plotting
plt.ion()
# Plot raw model output
if (raw):
# Compute plotting step size s.t. we plot every 100ms (=10Hz) (if possible)
s_rate = f['params']['s_rate'].value
if s_rate > 10:
p_rate = int(s_rate//10)
else:
p_rate = s_rate
# Get quantities for plotting
t = f['t'][::p_rate]
V = f['V'][idx,::p_rate].T
QV = f['QV'][idx,::p_rate].T
Beta = f['Beta'][idx,::p_rate].T
DA = f['DA'][idx,::p_rate].T
tmin = t.min() - t[1]
tmax = t.max() + t[1]
# Prepare window and plot stuff
fig = plt.figure(figsize=(10,7.5))
if figname != None: fig.canvas.set_window_title(figname)
plt.suptitle("Raw Model Output from "+fname,fontsize=18)
sp = plt.subplot(4,1,1)
plt.plot(t,V)
if doleg:
plt.legend(names)
plt.ylabel('mV',fontsize=16)
[ymin,ymax] = sp.get_ylim()
plt.yticks([ymin,0,ymax],fontsize=10)
sp.set_xlim(left=tmin,right=tmax)
plt.xticks([],fontsize=8)
plt.title("V",fontsize=16)
plt.draw()
sp = plt.subplot(4,1,2)
plt.plot(t,QV)
plt.ylabel('Firing',fontsize=16)
[ymin,ymax] = sp.get_ylim()
plt.yticks([0,ymax],fontsize=10)
sp.set_xlim(left=tmin,right=tmax)
plt.xticks([],fontsize=8)
plt.title("QV",fontsize=16)
plt.draw()
sp = plt.subplot(4,1,3)
plt.plot(t,Beta)
plt.ylabel('Gain Factor',fontsize=16)
[ymin,ymax] = sp.get_ylim()
plt.yticks([ymin,ymax],fontsize=10)
sp.set_xlim(left=tmin,right=tmax)
plt.xticks([],fontsize=8)
plt.title(r"$\beta$",fontsize=16)
plt.draw()
sp = plt.subplot(4,1,4)
plt.plot(t,DA)
plt.ylabel('mM',fontsize=16)
[ymin,ymax] = sp.get_ylim()
plt.yticks([ymin,ymax],fontsize=10)
sp.set_xlim(left=tmin,right=tmax)
plt.xticks(fontsize=10)
plt.xlabel("ms")
plt.title("DA",fontsize=16)
plt.draw()
# Plot raw model output
if (bold):
# Try to load the BOLD data from file
try:
BOLD = f['BOLD'][idx,:].T
except:
f.close()
raise ValueError("No BOLD data found in file "+fname+"!")
# Get x-extent of data and create x-ticks vector
xmax = BOLD.shape[0] + 1
xtv = np.arange(-1,xmax)
# Prepare window and plot stuff
fig = plt.figure(figsize=(10,7.5))
if figname != None: fig.canvas.set_window_title(figname)
plt.title("BOLD data "+fname,fontsize=18)
sp = plt.subplot(111)
plt.plot(BOLD)
if doleg:
plt.legend(names)
plt.xticks(xtv,fontsize=10)
[ymin,ymax] = sp.get_ylim()
plt.yticks([ymin,0,ymax],fontsize=10)
# Close file and return
f.close()
##########################################################################################
[docs]def make_D(target,source,names,values=None):
"""
Create matrix of afferent/efferent dopamine regions in the model
Parameters
----------
target : Python list or NumPy 1darray
Python list or NumPy 1darray of region names that are affected by dopamine release
(has to be the same length as `source`).
source : Python list or NumPy 1darray
Python list or NumPy 1darray of region names that steer dopamine release.
(has to be the same length as `target`).
names : Python list or NumPy 1darray
Names of all regions in the model. If `names` has length `N`, the resulting
dopamine connection matrix will by `N`-by-`N`.
values : Python list or NumPy 1darray
By default, dopamine links are binary, i.e., all links have unit weight. By passing
a `values` array, certain links can be emphasized (weight > 1) or weakened (weight < 1).
Entries of the `values` array have to be calibrated based on the values of `b_hi`, `b_lo`,
and `a`.
Returns
-------
D : NumPy 2darray
A `N`-by-`N` array. Every row that has a non-zero entry signifies a dopamine target,
and every non-zero column corresponds to a dopamine source.
Notes
-----
None
Examples
--------
For the sake of simplicity consider a brain "parcellation" consisting of three (bilateral) regions,
called `A`, `B`, and `C`. Thus, we define the following `names` array
>>> names = ['L_A','L_B','L_C','R_A','R_B','R_C']
Assume that in the left hemisphere dopamine is mainly concentrated in region `B`, its release is
steered by neural firing in region `A`. In the right hemisphere, dopamine release in region `C` is
depending on neural activity in area `B`. Then the `target` and `source` arrays are given by
>>> target = ['L_B','R_C']
>>> source = ['L_A','R_B']
Then the call
>>> D = make_D(target,source,names)
>>> D
array([[ 0., 0., 0., 0., 0., 0.],
[ 1., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 1., 0.]])
generates a 6-by-6 NumPy array `D`, that has non-zero entries in rows 2 and 6 (because
'L_B' is the second, and 'R_C' the sixth element of the list `names`), and columns 1 and 5
(because 'L_A' is the first, and 'R_B' the 5th element of the list `names`). Thus, the row/column
index of each non-zero entry in `D` has the format target-by-source.
See also
--------
None
"""
# Sanity checks
vnames = ["target","source","names"]
for tk, tsn in enumerate([target,source,names]):
try:
tsn = list(tsn)
except:
msg = "Inputs target, source and names have to be NumPy 1darrays or Python lists, not "+\
type(tsn).__name__
TypeError(msg)
for el in tsn:
if not isinstance(el,(str,unicode)):
raise ValueError("All elements of `"+vnames[tk]+"` have to be strings!")
if len(source) != len(target):
raise ValueError("Length of source and target lists/arrays does not match up!")
for tk, tsn in enumerate([target,source]):
for el in tsn:
if el not in names:
raise ValueError("Element `"+el+"` not found in `"+vnames[tk]+"`!")
if values != None:
values = np.array(values)
arrcheck(values,'vector','values')
if len(values) != len(target):
raise ValueError("Length of `values` list/array does not match up!")
else:
values = np.ones((len(target),))
# Convert (if we're having a NumPy array) `names` to a Python list
names = list(names)
# Get dimension we're dealing with here
N = len(names)
# Create dopamine matrix
D = np.zeros((N,N))
# Fill the matrix
for i in xrange(len(source)):
# Get row and column indices
row = names.index(target[i])
col = names.index(source[i])
# Matrix is targets-by-sources
D[row,col] = values[i]
return D
##########################################################################################
[docs]def make_bold(fname, stim_onset=None):
"""
Convert raw model output to BOLD signal
Parameters
----------
fname : str
File-name (including path if not in working directory) of HDF5 container that was generated
by `run_model`.
stim_onset : float
Time (in seconds) of stimulus onset. By default, onset/offset timings of
stimuli are stored in the HDF5 container generated by `run_model`. Only override
this setting, if you know what you are doing.
Returns
-------
Nothing : None
The computed BOLD signal is stored as dataset `BOLD` at the top level of the HDF5 container
specified by `fname`.
Notes
-----
Regional neural voltages are converted to BOLD time-series using the linear hemodynamic response
function proposed by Glover [1]_. For details consult the supporting information of our
paper [2]_.
References
----------
.. [1] Glover G. Deconvolution of Impulse Response in Event-Related BOLD FMRI.
NeuroImage 9: 416-429, 1999.
.. [2] S. Fuertinger, J. C. Zinn, and K. Simonyan. A Neural Population Model Incorporating
Dopaminergic Neurotransmission during Complex Voluntary Behaviors. PLoS Computational Biology,
10(11), 2014.
"""
# Make sure container exists and is valid
f = checkhdf(fname,peek=True)
try:
V = f['V'].value
except:
f.close()
raise ValueError("HDF5 file "+fname+" does not have the required fields!")
# Compute cycle length based on the sampling rate used to generate the file
N = f['params']['names'].size
s_rate = f['params']['s_rate'].value
n_cycles = f['params']['n_cycles'].value
len_cycle = f['params']['len_cycle'].value
cycle_idx = int(np.round(s_rate*len_cycle))
# Make sure `stim_onset` makes sense
if stim_onset != None:
scalarcheck(stim_onset,'stim_onset',bounds=[0,len_cycle])
# Get task from file to start sub-sampling procedure
task = f['params']['task'].value
# Compute cycle length based on the sampling rate used to generate the file
N = f['params']['names'].size
n_cycles = f['params']['n_cycles'].value
len_cycle = f['params']['len_cycle'].value
cycle_idx = int(np.round(s_rate*len_cycle))
# Compute step size and (if not provided by the user) compute stimulus onset time
dt = 1/s_rate
if stim_onset == None:
stim_onset = f['params']['stimulus'].value
# Use Glover's Hemodynamic response function as convolution kernel (with default length 32)
hrft = utils.lambdify_t(hrf.glover(utils.T))
hrf_kernel = np.hstack((np.zeros((int(s_rate*stim_onset),)),hrft(np.arange(0,32,dt))))
# Convolve the de-meaned model time-series with the kernel
convV = ndimage.filters.convolve1d((V.T - V.mean(axis=1)).T,hrf_kernel,mode='constant')
# Allocate space for BOLD signal
BOLD = np.zeros((N,n_cycles))
# Sub-sample convoluted data depending on task to get BOLD signal
if task == 'speech':
# Get interval to be considered for boldification
start = int(np.round(f['params']['speechoff'].value*s_rate))
stop = start + int(np.round(f['params']['acquisition'].value*s_rate))
elif task == 'rest':
# Get interval to be considered for boldification
start = 0
stop = int(np.round(f['params']['stimulus'].value*s_rate))
else:
raise ValueError("Don't know what to do for task "+task)
# Compute BOLD signal for all time points
for j in xrange(n_cycles):
BOLD[:,j] = convV[:,start:stop].mean(axis=1)
start += cycle_idx
stop += cycle_idx
# Re-scale the the signal
BOLD = BOLD*0.02
# Save it to the file
try:
f.create_dataset('BOLD',data=BOLD)
except:
f['BOLD'].write_direct(BOLD)
f.close()
##########################################################################################
[docs]def show_params(fname):
"""
Pretty-print all parameters used in a simulation
Parameters
----------
fname : str
Filename (including path if not in working directory) of HDF5 container that was generated
by `run_model`.
Returns
-------
Nothing : None
Notes
-----
None
See also
--------
None
"""
# Make sure container exists and is valid
f = checkhdf(fname,peek=True)
try:
par_grp = f['params']
f.close()
except:
raise ValueError("HDF5 file "+fname+" does not have the required fields!")
# Create a list for Texttable
tlist = []
for key in f['params'].keys():
tlist.append([key,f['params'][key].value])
# Close file
f.close()
# Print table
print"\nShowing parameters of file "+fname
table = Texttable()
table.set_deco(Texttable.HEADER)
table.set_cols_align(["l", "l"])
table.add_rows([["Parameter","Value"]],header=True)
table.add_rows(tlist,header=False)
print "\n"+table.draw()+"\n"
##########################################################################################
def arrcheck(arr,kind,varname,bounds=None):
"""
Local helper function performing sanity checks on arrays (1d/2d/3d)
"""
if not isinstance(arr,np.ndarray):
raise TypeError('Input `'+varname+'` must be a NumPy array, not '+type(arr).__name__+'!')
sha = arr.shape
if kind == 'tensor':
if len(sha) != 3:
raise ValueError('Input `'+varname+'` must be a `N`-by-`N`-by-`k` NumPy array')
if (min(sha[0],sha[1])==1) or (sha[0]!=sha[1]):
raise ValueError('Input `'+varname+'` must be a `N`-by-`N`-by-`k` NumPy array!')
dim_msg = '`N`-by-`N`-by-`k`'
elif kind == 'matrix':
if len(sha) != 2:
raise ValueError('Input `'+varname+'` must be a `N`-by-`N` NumPy array')
if (min(sha)==1) or (sha[0]!=sha[1]):
raise ValueError('Input `'+varname+'` must be a `N`-by-`N` NumPy array!')
dim_msg = '`N`-by-`N`'
elif kind == 'vector':
sha = arr.squeeze().shape
if len(sha) != 1:
raise ValueError('Input `'+varname+'` must be a NumPy 1darray')
if sha[0] <= 1:
raise ValueError('Input `'+varname+'` must be a NumPy 1darray of length `N`!')
dim_msg = ''
else:
print "Error checking could not be performed - something's wrong here..."
if not np.issubdtype(arr.dtype, np.number) or not np.isreal(arr).all():
raise ValueError('Input `'+varname+'` must be a real-valued '+dim_msg+' NumPy array!')
if np.isfinite(arr).min() == False:
raise ValueError('Input `'+varname+'` must be a real-valued NumPy array without Infs or NaNs!')
if bounds is not None:
if arr.min() < bounds[0] or arr.max() > bounds[1]:
raise ValueError("Values of input array `"+varname+"` must be between "+str(bounds[0])+\
" and "+str(bounds[1])+"!")
##########################################################################################
def scalarcheck(val,varname,kind=None,bounds=None):
"""
Local helper function performing sanity checks on scalars
"""
if not np.isscalar(val) or not plt.is_numlike(val):
raise TypeError("Input `"+varname+"` must be a scalar!")
if not np.isfinite(val) or not np.isreal(val):
raise ValueError("Input `"+varname+"` must be real and finite!")
if kind == 'int':
if (round(val) != val):
raise ValueError("Input `"+varname+"` must be an integer!")
if bounds is not None:
if val < bounds[0] or val > bounds[1]:
raise ValueError("Input scalar `"+varname+"` must be between "+str(bounds[0])+" and "+str(bounds[1])+"!")
##########################################################################################
def checkhdf(fname,peek=False):
"""
Local helper function performing sanity checks on file-names
"""
# Make sure `fname` is a valid file-path
if not isinstance(fname,(str,unicode)):
raise TypeError("Name of HDF5 file has to be a string!")
fname = str(fname)
if fname.find("~") == 0:
fname = os.path.expanduser('~') + fname[1:]
if not os.path.isfile(fname):
raise ValueError('File: '+fname+' does not exist!')
# Additionally, try opening the container if wanted
if peek:
try:
f = h5py.File(fname)
except:
raise ValueError("Cannot open "+fname+"!")
return f