neurolib-dev / neurolib Goto Github PK
View Code? Open in Web Editor NEWEasy whole-brain modeling for computational neuroscientists ๐ง ๐ป๐ฉ๐ฟโ๐ฌ
License: MIT License
Easy whole-brain modeling for computational neuroscientists ๐ง ๐ป๐ฉ๐ฟโ๐ฌ
License: MIT License
+
operator&
operatorThe BOLD model gets input from various different sources. For now, models which don't ouput a firing rate in Hz
can specify normalize_bold_input = True
to normalize their output (via normalize_bold_input_max = 50
) with a very simple predefined transformation that is baked into the BOLD model. This leads to some errors, specifically for models with negative outputs and non-oscillating states where there is no real amplitude to scale.
Therefore, it might be a better option to define a model-specific transformation function of the outputs such that the BOLD model can work with them. I am thinking something along the lines of
model.boldTransorm = (lambda x: (x + 50)*2.5)
and initialize the BOLD model in neurolib/models/model.py
with
bold.BOLDModel(self.params["N"], self.params["dt"], self.boldTransorm)
In run()
in model/bold/model.py
, the input would then be passed through this transformation function before passed to the bold timeintegration.
This way, every model would have the freedom of transforming their output to whatever they want the BOLD model to get as an input.
Thoughts?
from neurolib.utils.parameterSpace import ParameterSpace
shape = 4
pars = [np.random.rand(shape, shape) for i in range(10)]
parameters = ParameterSpace({"p": pars})
parameters.getRandom()
Returns
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-32-1b6c4906dd1f> in <module>
----> 1 parameters.getRandom()
~/Documents/PhD/projects/neurolib/neurolib/utils/parameterSpace.py in getRandom(self, safe)
97 else:
98 for key, value in self.parameters.items():
---> 99 randomPar[key] = np.random.choice(value)
100 return randomPar
101
mtrand.pyx in numpy.random.mtrand.RandomState.choice()
ValueError: a must be 1-dimensional
Needs fix.
WIP, maybe netcdf?
Comments were written without black
style guide in mind. This introduces ugly line breaks in the variable definitions and should be fixed. All code files should be cleaned up with black
in mind.
Ugly line breaks found in
neurolib/models/hopf/timeIntegration.py
neurolib/models/aln/model.py
neurolib/models/aln/timeIntegration.py
I was thinking about how to support sampling_dt
parameter. Reasons:
dt
should be 0.01) or ALN-thalamus mini-network (I am using 0.005 rn for best results)Currently, I am leaning towards:
sampling_dt
parameter for all models in loadDefaultParams.py
(e.g. default value at 10ms)?sampling_dt = None
would be valid and no subsampling would be doneModel.setOutput
functionthis way, I believe all neurolib
features should work (with exploration / evolution saving subsampled timeseries)
@caglorithm what do you think? I can implement it no problemo, just wanted to be sure it won't break anything and you know the best neurolib
's internals...
...
The default parameters of the WC model have this weird difference in E vs I. I'm not sure why this is in there, maybe @ChristophMetzner can tell, but I think this needs to change to be zero for both.
# values of the external inputs
params.exc_ext = 1.0 * np.ones((params.N,))
params.inh_ext = np.zeros((params.N,))
Use case: you want to explore the phase space over connectivity. Since the connectivity matrix is really a matrix, this won't work with the ParameterSpace
class, since it only allows scalars as parameters to search over.
Optimization can be done since you need to provide a cost function, and there you can do anything, e.g. build a connectivity matrix and pass it to the model simply with model.params[*connectivity|local]=<new matrix>
inside a cost function and it will work...
However, BoxSearch is not usable... BoxSearch simply calls model.run()
with each new parameter set, but matrices cannot be part of ParameterSpace
, i.e., parameter set during exploration.
So we need to think of a way how to allow this.
A run
, an exploration
and an optimization
should log all details about the run, including the version of neurolib
used, the parameters of the models and the specific run parameters, such as duration
, the hdf filename, the trajectory name, duration of exploration, number of cores used, machine name, ...
Working on this on branch fix/seed
.
by calling
for i, val in enumerate(inputs):
if i==0:
wc.run()
else:
wc.run(continue_run=True)
the following error is returned as i=2:
Traceback (most recent call last):
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/multiprocessing/process.py", line 297, in _bootstrap
self.run()
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/multiprocessing/process.py", line 99, in run
self._target(*self._args, **self._kwargs)
File "/mnt/antares_raid/home/ronjastroms/Bifurcation.py", line 139, in continuedRun_wc
self.model.run(continue_run=True)
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/site-packages/neurolib/models/model.py", line 227, in run
self.integrate(append_outputs=append, simulate_bold=bold)
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/site-packages/neurolib/models/model.py", line 267, in integrate
t, *variables = self.integration(self.params)
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/site-packages/neurolib/models/wc/timeIntegration.py", line 164, in timeIntegration
tau_adap,
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/site-packages/numba/dispatcher.py", line 401, in _compile_for_args
error_rewrite(e, 'typing')
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/site-packages/numba/dispatcher.py", line 344, in error_rewrite
reraise(type(e), e, None)
File "/mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/site-packages/numba/six.py", line 668, in reraise
raise value.with_traceback(tb)
numba.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of Function(<built-in function setitem>) with argument(s) of type(s): (array(float64, 2d, C), UniTuple(int64 x 2), array(float64, 1d, C))
* parameterized
In definition 0:
All templates rejected with literals.
In definition 1:
All templates rejected without literals.
In definition 2:
All templates rejected with literals.
In definition 3:
All templates rejected without literals.
In definition 4:
All templates rejected with literals.
In definition 5:
All templates rejected without literals.
In definition 6:
All templates rejected with literals.
In definition 7:
All templates rejected without literals.
In definition 8:
All templates rejected with literals.
In definition 9:
All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of setitem at /mnt/antares_raid/home/ronjastroms/anaconda3/envs/trial/lib/python3.7/site-packages/neurolib/models/wc/timeIntegration.py (270)
File "anaconda3/envs/trial/lib/python3.7/site-packages/neurolib/models/wc/timeIntegration.py", line 270:
def S_I(x):
<source elided>
# Euler integration
excs[no, i] = excs[no, i - 1] + dt * exc_rhs
^
There should be a way to list all available models. Similar to how TVB handles it would be nice: http://docs.thevirtualbrain.com/_modules/tvb/simulator/models.html
Dear all,
I've tried to reproduce the following example:
from neurolib.models.aln import ALNModel
from neurolib.optimize.exploration import BoxSearch
from neurolib.utils.parameterSpace import ParameterSpace
model = ALNModel()
parameters = ParameterSpace({"mue_ext_mean": np.linspace(0, 3, 21),
"mui_ext_mean": np.linspace(0, 3, 21)})
search = BoxSearch(model, parameters)
search.run()
search.loadResults()
in the line search.run()
I have the following error:
Exception has occurred: PicklingError
Can't pickle <class 'neurolib.utils.parameterSpace.ParameterSpace'>: it's not the same object as neurolib.utils.parameterSpace.ParameterSpace
File "/neurolib/optimize/exploration/exploration.py", line 281, in run self.env.run(self.evalFunction)
After some trials and errors, I've changed multiproc=True
to multiproc=False
here and everything works fine.
Looking forward to your comments!
Valerii
Hello,
I tried to run the following code which was mostly taken from here with following output. My Problem is that it even though it seems like it runs the BoxSearch
correctly, the resulting .hdf file is missing the results directory. This issues becomes visible when running search.loadResults
.
"""
Some code taken from https://github.com/caglarcakan/sleeping_brain
Copyright (c) 2020 Caglar Cakan
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""
import logging
import warnings
import matplotlib.pyplot as plt
import numpy as np
import neurolib.optimize.exploration.explorationUtils as eu
import neurolib.utils.functions as func
from neurolib.models.aln import ALNModel
from neurolib.optimize.exploration import BoxSearch
from neurolib.utils.loadData import Dataset
from neurolib.utils.parameterSpace import ParameterSpace
from neurolib.utils.stimulus import construct_stimulus
logger = logging.getLogger()
warnings.filterwarnings("ignore")
logger.setLevel(logging.INFO)
plt.rcParams['text.usetex'] = False
plt.rcParams['svg.fonttype'] = 'none'
plt.style.reload_library()
plt.style.use("seaborn-white")
plt.rcParams['image.cmap'] = 'plasma'
def evaluateSimulation(traj):
# get the model from the trajectory using `search.getModelFromTraj(traj)`
model = search.getModelFromTraj(traj)
# initiate the model with random initial contitions
model.randomICs()
defaultDuration = model.params['duration']
# -------- stage wise simulation --------
# Stage 3: full and final simulation
# ---------------------------------------
model.params['duration'] = defaultDuration
rect_stimulus = construct_stimulus(stim="rect",
duration=model.params.duration,
dt=model.params.dt)
model.params['ext_exc_current'] = rect_stimulus * 5.0
model.run()
# up down difference
state_length = 2000
# last_state = (model.t > defaultDuration - state_length)
# time period in ms where we expect the down-state
down_window = ((defaultDuration/2 - state_length < model.t)
& (model.t < defaultDuration/2))
# and up state
up_window = ((defaultDuration - state_length < model.t)
& (model.t < defaultDuration))
up_state_rate = np.mean(model.output[:, up_window], axis=1)
down_state_rate = np.mean(model.output[:, down_window], axis=1)
up_down_difference = np.max(up_state_rate - down_state_rate)
# check rates!
max_amp_output = np.max(
np.max(model.output[:, up_window], axis=1)
- np.min(model.output[:, up_window], axis=1)
)
max_output = np.max(model.output[:, up_window])
model_frs, model_pwrs = func.getMeanPowerSpectrum(model.output,
dt=model.params.dt,
maxfr=40,
spectrum_windowsize=10)
model_frs, model_pwrs = func.getMeanPowerSpectrum(
model.output[:, up_window], dt=model.params.dt,
maxfr=40, spectrum_windowsize=5)
domfr = model_frs[np.argmax(model_pwrs)]
result = {
"max_output": max_output,
"max_amp_output": max_amp_output,
"domfr": domfr,
"up_down_difference": up_down_difference
}
search.saveToPypet(result, traj)
return
ds = Dataset("gw")
model = ALNModel(Cmat=ds.Cmats[0], Dmat=ds.Dmats[0])
model.params['dt'] = 0.1
model.params['duration'] = 20 * 1000 # ms
# add custom parameter for downsampling results
# 10 ms sampling steps for saving data, should be multiple of dt
model.params['save_dt'] = 10.0
model.params["tauA"] = 600.0
model.params["sigma_ou"] = 0.0
model.params["b"] = 20.0
model.params["Ke_gl"] = 300.0
model.params["signalV"] = 80.0
parameters = ParameterSpace({"mue_ext_mean": np.linspace(0.0, 4, 2),
"mui_ext_mean": np.linspace(0.0, 4, 2),
"b": [20.0]
})
search = BoxSearch(evalFunction=evaluateSimulation, model=model,
parameterSpace=parameters,
filename='exploration-8.0-brain.hdf')
if __name__ == '__main__':
search.run()
fname = "C:\\Users\\HannoChan\\Documents\\uni\\neural_information_processing_project\\NI-project\\data\\hdf\\exploration-8.0-brain.hdf"
search.loadResults(filename=fname, all=False)
print(search.dfResults)
plot_key_label = "Max. $r_E$ [Hz]"
eu.plotExplorationResults(search.dfResults,
par1=['mue_ext_mean', 'Input to E [nA]'],
par2=['mui_ext_mean', 'Input to I [nA]'],
by=['b'],
plot_key='max_output',
plot_clim=[0.0, 80.0],
nan_to_zero=False,
plot_key_label=plot_key_label,
one_figure=True,
multiply_axis=0.2,
contour=["max_amp_output", "up_down_difference"],
contour_color=[['white'], ['springgreen']],
contour_levels=[[10], [10]],
contour_alpha=[1.0, 1.0],
contour_kwargs={0: {"linewidths": (5,)},
1: {"linestyles": "--",
"linewidths": (5,)}},
# alpha_mask="relative_amplitude_BOLD",
mask_threshold=0.1,
mask_alpha=0.2,
savename="gw_brain_nap_001.pdf")
INFO:root:Loading dataset gw from C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\utils\..\data\datasets\gw. INFO:root:Dataset gw loaded. INFO:root:aln: Model initialized. INFO:root:Number of processes: 12 MainProcess pypet.storageservice.HDF5StorageService INFO I will use the hdf5 file `./data/hdf/exploration-8.0-brain.hdf`. MainProcess pypet.environment.Environment INFO Environment initialized. MainProcess root INFO Number of parameter configurations: 4 MainProcess root INFO BoxSearch: Environment initialized. MainProcess pypet.environment.Environment INFO I am preparing the Trajectory for the experiment and initialise the store. MainProcess pypet.environment.Environment INFO Initialising the storage for the trajectory. MainProcess pypet.storageservice.HDF5StorageService INFO Initialising storage or updating meta data of Trajectory `results-2021-05-27-09H-25M-35S`. MainProcess pypet.storageservice.HDF5StorageService INFO Finished init or meta data update for `results-2021-05-27-09H-25M-35S`. MainProcess pypet.environment.Environment INFO ************************************************************ STARTING runs of trajectory `results-2021-05-27-09H-25M-35S`. ************************************************************I am running a Windows 10 computer with Python 3.7.9. We also tried running the code on a MacBook on which it worked flawlessly.MainProcess pypet.storageservice.HDF5StorageService INFO Initialising storage or updating meta data of Trajectory `results-2021-05-27-09H-25M-35S`. MainProcess pypet.storageservice.HDF5StorageService INFO Finished init or meta data update for `results-2021-05-27-09H-25M-35S`. MainProcess pypet.environment.Environment INFO Starting multiprocessing with at most 12 processes running at the same time. MainProcess pypet INFO PROGRESS: Finished 0/4 runs [ ] 0.0% INFO:root:Loading dataset gw from C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\utils\..\data\datasets\gw. INFO:root:Dataset gw loaded. INFO:root:aln: Model initialized. INFO:root:Number of processes: 12 Process-1 pypet.storageservice.HDF5StorageService INFO I will use the hdf5 file `./data/hdf/exploration-8.0-brain.hdf`. Process-1 pypet.environment.Environment INFO Environment initialized. Process-1 root INFO Number of parameter configurations: 4 Process-1 root INFO BoxSearch: Environment initialized. INFO:root:Loading dataset gw from C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\utils\..\data\datasets\gw. INFO:root:Dataset gw loaded. INFO:root:aln: Model initialized. INFO:root:Number of processes: 12 Process-2 pypet.storageservice.HDF5StorageService INFO I will use the hdf5 file `./data/hdf/exploration-8.0-brain.hdf`. Process-2 pypet.environment.Environment INFO Environment initialized. Process-2 root INFO Number of parameter configurations: 4 Process-2 root INFO BoxSearch: Environment initialized. INFO:root:Loading dataset gw from C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\utils\..\data\datasets\gw. INFO:root:Dataset gw loaded. INFO:root:aln: Model initialized. INFO:root:Number of processes: 12 Process-3 pypet.storageservice.HDF5StorageService INFO I will use the hdf5 file `./data/hdf/exploration-8.0-brain.hdf`. Process-3 pypet.environment.Environment INFO Environment initialized. Process-3 root INFO Number of parameter configurations: 4 Process-3 root INFO BoxSearch: Environment initialized. INFO:root:Loading dataset gw from C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\utils\..\data\datasets\gw. INFO:root:Dataset gw loaded. INFO:root:aln: Model initialized. INFO:root:Number of processes: 12 Process-4 pypet.storageservice.HDF5StorageService INFO I will use the hdf5 file `./data/hdf/exploration-8.0-brain.hdf`. Process-4 pypet.environment.Environment INFO Environment initialized. Process-4 root INFO Number of parameter configurations: 4 Process-4 root INFO BoxSearch: Environment initialized. MainProcess pypet INFO PROGRESS: Finished 1/4 runs [===== ] 25.0%, remaining: 0:00:42 MainProcess pypet INFO PROGRESS: Finished 2/4 runs [========== ] 50.0%, remaining: 0:00:17 MainProcess pypet INFO PROGRESS: Finished 3/4 runs [=============== ] 75.0%, remaining: 0:00:06 MainProcess pypet INFO PROGRESS: Finished 4/4 runs [====================]100.0% MainProcess pypet.storageservice.HDF5StorageService INFO Initialising storage or updating meta data of Trajectory `results-2021-05-27-09H-25M-35S`. MainProcess pypet.storageservice.HDF5StorageService INFO Finished init or meta data update for `results-2021-05-27-09H-25M-35S`. MainProcess pypet.environment.Environment INFO ************************************************************ FINISHED all runs of trajectory `results-2021-05-27-09H-25M-35S`. ************************************************************ MainProcess pypet.environment.Environment INFO ************************************************************ STARTING FINAL STORING of trajectory `results-2021-05-27-09H-25M-35S` ************************************************************ MainProcess pypet.storageservice.HDF5StorageService INFO Start storing Trajectory `results-2021-05-27-09H-25M-35S`. MainProcess pypet.storageservice.HDF5StorageService INFO Storing branch `config`. MainProcess pypet.storageservice.HDF5StorageService INFO Storing branch `parameters`. MainProcess pypet.storageservice.HDF5StorageService INFO Finished storing Trajectory `results-2021-05-27-09H-25M-35S`. MainProcess pypet.environment.Environment INFO ************************************************************ FINISHED FINAL STORING of trajectory `results-2021-05-27-09H-25M-35S`. ************************************************************ MainProcess pypet.environment.Environment INFO All runs of trajectory `results-2021-05-27-09H-25M-35S` were completed successfully. MainProcess root INFO Loading results from C:\Users\HannoChan\Documents\uni\neural_information_processing_project\NI-project\data\hdf\exploration-8.0-brain.hdf MainProcess root INFO Analyzing trajectory results-2021-05-27-09H-25M-35S MainProcess pypet.storageservice.HDF5StorageService INFO I will use the hdf5 file `C:\Users\HannoChan\Documents\uni\neural_information_processing_project\NI-project\data\hdf\exploration-8.0-brain.hdf`. MainProcess pypet.storageservice.HDF5StorageService INFO Loading trajectory `results-2021-05-27-09H-25M-35S`. MainProcess pypet.storageservice.HDF5StorageService INFO Loading branch `config` in mode `2`. MainProcess pypet.storageservice.HDF5StorageService INFO Loading branch `parameters` in mode `2`. MainProcess root INFO Creating `dfResults` dataframe ... MainProcess root INFO Aggregating results to `dfResults` ... 0%| | 0/4 [00:00<?, ?it/s]MainProcess pypet.storageservice.HDF5StorageService ERROR Cannot find `results` the hdf5 node `results` does not exist! MainProcess pypet.storageservice.HDF5StorageService ERROR Failed loading `ResultGroup results: ` MainProcess pypet.naturalnaming.NaturalNamingInterface ERROR Error while auto-loading `run_00000000` under `results`. 0%| | 0/4 [00:00<?, ?it/s] Traceback (most recent call last): File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\pypet\storageservice.py", line 948, in load self._tree_load_sub_branch(stuff_to_load, *args, **kwargs) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\pypet\storageservice.py", line 2115, in _tree_load_sub_branch name=hdf5_group_name) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\tables\file.py", line 1635, in get_node node = where._v_file._get_node(nodepath) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\tables\file.py", line 1590, in _get_node node = self._node_manager.get_node(nodepath) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\tables\file.py", line 432, in get_node node = self.node_factory(key) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\tables\group.py", line 1178, in _g_load_child node_type = self._g_check_has_child(childname) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\tables\group.py", line 395, in _g_check_has_child % (self._v_pathname, name)) tables.exceptions.NoSuchNodeError: group ``/`` does not have a child named ``/results-2021-05-27-09H-25M-35S/results`` During handling of the above exception, another exception occurred: Traceback (most recent call last): File "RESULT-bifurcation-diagrams.py", line 136, in <module> search.loadResults(filename=fname, all=False) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\optimize\exploration\exploration.py", line 340, in loadResults self.aggregateResultsToDfResults() File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\optimize\exploration\exploration.py", line 369, in aggregateResultsToDfResults result = self.getRun(runId) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\optimize\exploration\exploration.py", line 543, in getRun return pu.getRun(runId, pypetTrajectory, pypetShortNames=pypetShortNames) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\neurolib\utils\pypetUtils.py", line 65, in getRun pypetTrajectory.results[runId].f_load() File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\pypet\naturalnaming.py", line 2984, in __getitem__ with_links=self.v_root.v_with_links) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\pypet\naturalnaming.py", line 2329, in _get try_auto_load_directly2) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\pypet\naturalnaming.py", line 2422, in _perform_get load_data=pypetconstants.LOAD_DATA) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\pypet\naturalnaming.py", line 3376, in f_load_child max_depth=max_depth) File "C:\Users\HannoChan\Documents\uni\neural_information_processing_project\.venv\lib\site-packages\pypet\storageservice.py", line 958, in load raise pex.DataNotInStorageError(repr(exc)) pypet.pypetexceptions.DataNotInStorageError: NoSuchNodeError('group ``/`` does not have a child named ``/results-2021-05-27-09H-25M-35S/results``')
As I could not find a solution to my problem, I was wondering if somebody else had an idea.
Looking forward to your comments!
Best Regards
HannoBOLD output is always appended
Output to
BOLDModel
is always appended, even when theModel
is rerun withmodel.run(append_outputs=False)
. This leas to a very strange behaviour. Workaround is to reinitialize theModel
for every new run.Crash if delays change but model is not reinitialized
If the model it not reinitialized, the maxDelay is not recalculated but assumed to be the same value as before. This leads to a crash if the delays decrease (because state_variables become shorter than the initial_values). Code that reproduces this:
from neurolib.models.aln import ALNModel from neurolib.utils.loadData import Dataset ds = Dataset("hcp") aln = ALNModel(Cmat=ds.Cmat, Dmat=ds.Dmat) aln.params['signalV'] = 10.0 aln.run(chunkwise=True, chunksize=10000) aln.params['signalV'] = 1.0 aln.run(chunkwise=True, chunksize=10000)Error:
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-81-e10453abeca9> in <module> 2 aln.run(chunkwise=True, chunksize=10000) 3 aln.params['signalV'] = 1.0 ----> 4 aln.run(chunkwise=True, chunksize=10000) ~/Documents/PhD/projects/neurolib/neurolib/models/model.py in run(self, inputs, onedt, chunkwise, chunksize, bold, append_outputs) 145 logging.warn(f"{self.name}: BOLD model not initialized, not simulating BOLD. Use `run(bold=True)`") 146 bold = False --> 147 self.integrate_chunkwise(chunksize=chunksize, bold=bold, append_outputs=append_outputs) 148 return 149 ~/Documents/PhD/projects/neurolib/neurolib/models/model.py in integrate_chunkwise(self, chunksize, bold, append_outputs) 187 currentChunkSize += self.max_delay + 1 188 --> 189 self.autochunk(duration=currentChunkSize * dt, append_outputs=append_outputs) 190 191 if bold and self.bold_initialized: ~/Documents/PhD/projects/neurolib/neurolib/models/model.py in autochunk(self, inputs, duration, append_outputs) 219 220 # run integration --> 221 self.integrate(append_outputs=append_outputs) 222 223 # reset initial conditions to last state ~/Documents/PhD/projects/neurolib/neurolib/models/model.py in integrate(self, append_outputs) 155 """ 156 # run integration --> 157 t, *variables = self.integration(self.params) 158 159 # save time array ~/Documents/PhD/projects/neurolib/neurolib/models/aln/timeIntegration.py in timeIntegration(params) 209 rates_inh[:, startind:] = np.random.standard_normal((N, len(t) - startind)) 210 --> 211 rates_exc[:, :startind] = rates_exc_init 212 rates_inh[:, :startind] = rates_inh_init 213 ValueError: could not broadcast input array from shape (80,249) into shape (80,2484)Exploraiton with random ICs
TODO:
Exploration.useRandomICs
not implemented yetTodo: Implement Kuramoto model
bug in removing self connections
the piece of code in
loadDefaultParams.py
for various models which is supposed to get rid of self connections:params.Cmat = Cmat.copy() # coupling matrix np.fill_diagonal(Cmat, 0) # no self connectionsis not working, since the second line fills diagonal in the original
Cmat
variable, and not insideparams.Cmat
Discussion: Should models be in derivative form?
Hey peoples,
I've been thinking about this, specifically in the context of @jajcayn's recent findings that the thalamic model is really sensitive todt
in the Euler integration scheme. It would be great if we could just run models in other schemes as well, to be more flexible. But the current "philosophy" is to have an integrator for each model, i.e. every model deals with it for itself.Implementing each model as in an
dx/dt = f(x)
kind of form would allow us to build different integrators around it. Obviously this means a lot of work: who integrates the noise? How do we handle delays?This would be a big effort but maybe it's the safest way to go... thoughts?
Create `MultiModel` examples
The plan is to have three notebooks showcasing
MultiModel
.Basic run + connected model (ALN + thalamus)
- done ?
- basic
MultiModel
run withneurolib
wrapper, show star parameters and how they work- more advanced heterogeneous model (2 nodes -> brain net), use example with thalamus + ALN
Creating a new
MultiModel
model
- done ?
- build a Jansen-Rit model
- show and build
NeuralMasses
for JR model- build a single node for the model
- showcase full network - JR network with few nodes, and-or JR for cortex + thalamus
Advanced topics
- done ?
- showcase both backends and why there are two backends and which to which
- show exploration with
MultiModel
- show evolution with
MultiModel
- show how already implemented models can be easily adapted by subclassing
NeuralMasses
and/orNodes
Easier loading of datasets
- The file
neurolib/utils/loadData.py
should be called something likedataset.py
because it only includes that class.- Data loading should be easier. Right now we have
ds = Dataset("gw")
and thenmodel = ALNModel(Cmat=ds.Cmat, Dmat=ds.Dmat)
but this could also bemodel = ALNModel(ds=ds)
Submission to Conda-forge
Hi,
I was wondering if you are planning to put the library on Conda-forge? I'm trying to use the library in Julia instead of Python but it's not on the default channel so it's a bit of headache.
Thanks,
Mahtaaln: add example notebook with external stimulation
c/p form Issue #10
- Add stimulation functionality (i.e. you specify tDCS or tACS stimulation parameters and under the hood this is automatically translated to model parameter changes for the specific model the user chose)
This is actually already possible! At least with the
aln
model, you can set the external input parameters for stimulation before each run and they will be integrated over time.They are called
ext_exc_current
/ext_inh_current
for current input innA
, andext_exc_rate
/ext_inh_rate
from other spike rate sources inkHz
.I actually wanted to know how fast and easy I could implement this with
neurolib
! Let's see what happens if we ramp up a current from0
to0.1 nA
:from neurolib.models.aln import ALNModel import matplotlib.pyplot as plt import numpy as np aln = ALNModel() # comes with default parameters aln.run() plt.plot(aln.t, aln.rates_exc.T, c='k', label='no stim', lw=3) aln.params['ext_exc_current'] = np.linspace(0, 0.1, aln.params['duration']/aln.params['dt']) aln.run() plt.plot(aln.t, aln.rates_exc.T, c='r', alpha=0.5, label='ramp stim', lw=2) plt.legend()
๐คฏ
PS: I'm gonna add this to the examples!
Imports don't work in python 3.6
Chunkwise integration if chunksize is set
model.run(chunksize=10000)
should implychunkwise=True
.Feature request: progress bar reports ETA of simulation run
tqdm
offers a very nice wrapper around for loops for easy progress bars. This could be implemented easily inModel.integrate_chunkwise()
.Wilson Cowan handles x_ext parameter differently than the rest
It expects a list or an array (node-wise) whereas other models also accept a scalar. Should have either support for both or be the same as all other models.
Naming conventions for `MultiModel`
We would need to agree on naming conventions for anything relevant to neuroscience in
neurolib
in general, mainly for the newMultiModel
framework.This includes:
- how to call individual masses (Freeman KI set, i.e. neural population of the same type, typically excitatory or inhibitory) - now they are called
Mass
(but might bePopulation
, or anything similar)- how to call a collection of masses (Freeman KII set) - now they are
Node
- how to call brain network -
Network
is used now- consolidate parameter names between different
MultiModel
models (in particular if the parameter stands for the same parameter or refers to the same biophysical process) and betweenMultiModel
model andneurolib
native - in particular for ALN- possibly anything else that could cause confusion...
I would address this in a separate PR as then it'll be clear from the git history what exactly are the naming conventions
Support for physical units?
neurolib.utils.signal
has cute and simple support for units. Example:class RatesSignal(Signal): name = "Population firing rate" label = "q" signal_type = "rate" unit = "Hz"Model inputs & outputs could also have (optional) physical units. This could allow for translation of one value to another.
Some errors running the tutorials
Hi Caglar,
What an amazing tool! I've managed to simulate FC using your notebooks, but I've ran into some problems just running the notebooks as they are using the example code and data. Maybe it is a problem with my installation or virtual environment (conda environment attached, python 3.8)?
[neurolib2_packages.txt] (https://github.com/neurolib-dev/neurolib/files/5050707/neurolib2_packages.txt)With example-1-aln-paramter-exploration I get the following error:
parameters = ParameterSpace({"mue_ext_mean": np.linspace(0, 3, 21), "mui_ext_mean": np.linspace(0, 3, 21)})
search = BoxSearch(aln, parameters, filename="example-1.hdf")
ValueError Traceback (most recent call last)
in
1 parameters = ParameterSpace({"mue_ext_mean": np.linspace(0, 3, 21), "mui_ext_mean": np.linspace(0, 3, 21)})
2 # info: chose np.linspace(0, 3, 21) or more, values here are low for testing
----> 3 search = BoxSearch(aln, parameters)~/.conda/envs/neurolib2/lib/python3.8/site-packages/neurolib/optimize/exploration/exploration.py in init(self, model, parameterSpace, evalFunction, filename, saveAllModelOutputs)
71 # bool to check whether pypet was initialized properly
72 self.initialized = False
---> 73 self.initializeExploration(self.filename)
74
75 def initializeExploration(self, filename="exploration.hdf"):~/.conda/envs/neurolib2/lib/python3.8/site-packages/neurolib/optimize/exploration/exploration.py in initializeExploration(self, filename)
93
94 # set up the pypet environment
---> 95 env = pypet.Environment(
96 trajectory=trajectoryName,
97 filename=trajectoryfilename,~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/configparsing.py in new_func(env, *args, **kwargs)
18 # Pass the config data to the kwargs
19 new_kwargs = config_interpreter.interpret()
---> 20 init_func(env, *args, **new_kwargs)
21 # Add parameters and config data from the.ini
file
22 config_interpreter.add_parameters(env.traj)~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/decorators.py in new_func(*args, **kwargs)
161 kwargs[new_name] = value
162
--> 163 return func(*args, **kwargs)
164
165 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/decorators.py in new_func(*args, **kwargs)
161 kwargs[new_name] = value
162
--> 163 return func(*args, **kwargs)
164
165 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/decorators.py in new_func(*args, **kwargs)
161 kwargs[new_name] = value
162
--> 163 return func(*args, **kwargs)
164
165 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/decorators.py in new_func(*args, **kwargs)
161 kwargs[new_name] = value
162
--> 163 return func(*args, **kwargs)
164
165 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/decorators.py in new_func(*args, **kwargs)
161 kwargs[new_name] = value
162
--> 163 return func(*args, **kwargs)
164
165 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/decorators.py in new_func(*args, **kwargs)
161 kwargs[new_name] = value
162
--> 163 return func(*args, **kwargs)
164
165 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/utils/decorators.py in new_func(*args, **kwargs)
161 kwargs[new_name] = value
162
--> 163 return func(*args, **kwargs)
164
165 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/pypetlogging.py in new_func(self, *args, **kwargs)
170 _change_logging_kwargs(kwargs)
171
--> 172 return func(self, *args, **kwargs)
173
174 return new_func~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/environment.py in init(self, trajectory, add_time, comment, dynamic_imports, wildcard_functions, automatic_storing, log_config, log_stdout, report_progress, multiproc, ncores, use_scoop, use_pool, freeze_input, timeout, cpu_cap, memory_cap, swap_cap, niceness, wrap_mode, queue_maxsize, port, gc_interval, clean_up_runs, immediate_postproc, resumable, resume_folder, delete_resume, storage_service, git_repository, git_message, git_fail, sumatra_project, sumatra_reason, sumatra_label, do_single_runs, graceful_exit, lazy_debug, **kwargs)
1143 log_stdout=log_stdout,
1144 report_progress=report_progress)
-> 1145 self._logging_manager.check_log_config()
1146 self._logging_manager.add_null_handler()
1147 self._set_logger()~/.conda/envs/neurolib2/lib/python3.8/site-packages/pypet/pypetlogging.py in check_log_config(self)
521 if isinstance(self.log_config, str):
522 if not os.path.isfile(self.log_config):
--> 523 raise ValueError('Could not find the logger init file '
524 '%s
.' % self.log_config)
525 parser = NoInterpolationParser()ValueError: Could not find the logger init file
./neurolib/utils/pypet_logging.ini
.Any help would be appreciated!
Many thanks,
Joffchunkwise and normal runs differ in their output if noise is enabled
This is because every chunkwise run calls
integreate()
which generates independent noise for each run.Multimodel: Model input / stimulus shape
As per #171 all inputs should be
(space, time)
, even forMultiModel
. Make sure, this is the case and make sure also all stimuli fromModelInput
classes respect this.ParameterSpace can't initialize if n_pars=2 and values are not 1d (assumes boundaries)
from neurolib.utils.parameterSpace import ParameterSpace shape = 2 n_pars = 2 pars = [np.random.rand(shape, shape) for i in range(n_pars)] parameters = ParameterSpace({"p": pars})Returns
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-34-8e5c938267d2> in <module> ----> 1 parameters = ParameterSpace({"p": pars}) ~/Documents/PhD/projects/neurolib/neurolib/utils/parameterSpace.py in __init__(self, parameters, parameterValues, kind, allow_star_notation) 36 parameters, dict 37 ), "Parameters must be a dict, if no values are given in `parameterValues`" ---> 38 processedParameters = self._processParameterDict(parameters) 39 else: 40 # check if all names are strings ~/Documents/PhD/projects/neurolib/neurolib/utils/parameterSpace.py in _processParameterDict(self, parameters) 180 if self.kind == "bound": 181 # check the boundaries --> 182 self._validate_param_bounds(list(parameters.values())) 183 184 # set all parameters as attributes for easy access ~/Documents/PhD/projects/neurolib/neurolib/utils/parameterSpace.py in _validate_param_bounds(self, param_bounds) 142 # check every single parameter bound 143 for single_bound in param_bounds: --> 144 self._validate_single_bound(single_bound) 145 146 def _processParameterDict(self, parameters): ~/Documents/PhD/projects/neurolib/neurolib/utils/parameterSpace.py in _validate_single_bound(single_bound) 129 ), "An error occured while validating the ParameterSpace of kind 'bound': Only two bounds (min and max) are allowed" 130 assert ( --> 131 single_bound[1] > single_bound[0] 132 ), "An error occured while validating the ParameterSpace of kind 'bound': Minimum parameter value can't be larger than the maximum!" 133 ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
Works for
n_pars = 3
.ParameterSpace
assumes boundaries which it shouldn'tWhat should neurolib be?
Apart from the main purpose of this lib (that being the meso-/large- scale modelling) what other functionality we want to address?
optimisation - ?
- basic grid search
- evolutionary algorithm (DEAP)
- Gaussian Processes Surrogate
- deep density estimator [??]
plotting tools - ?
basic analysis tools?
- e.g. time-frequency stuff - ?
- cross-frequency algos - ?
- network analysis algos (small-world, clustering, avg. path length, etc...) - ?
kittens - ?
Simplify `construct_stimulus` and document stimulus in `example-0.5-aln-external-stimulus`
A couple of things:
- delete
construct_stimulus
function fromutils.stimulus
; it has three modes and all of them can be achieved withModelInput
class:
ac
is the same asSinusoidalInput
dc
is the same asStepInput
rect
can be easily done with a combination ofExponentialInput
andStepInput
with concatenation.. I will createRectifiedInput
directly in the file- create a shortcut for creating stimulus to core neurolib models... before it was something like
construct_stimulus("ac", duration=model.params.duration, dt=model.params.dt, stim_amp=1.0, stim_freq=1)
; I proposeSinusoidalInput(amplitude=1.0, period=1000.0).to_model(model)
, such that duration and dt of stimulation will be directly read frommodel.params
- revamp example-0.5-aln-external-stimulus:
- swap all calls to
construct_stimulus
with new functions- create a section about
ModelInput
, what is implemented, how it looks like, what are the options (to_model
,as_array
,as_cubic_splines
)- showcase
SummedInput
andConcatenatedInput
(need to wait for #162)Implement base class for neural masses
It would serve as a base building block for all models!
I.e. should be as general as possible and easy-to-read, so new models can be implemented with ease.See also #8 discussion for notes on how to implement.
How to save results of single simulations and explorations?
Right now, each model just saves its results (meaning the integrated state variables) into their own attributes. Clearly, this is not a good approach since different types of models can have very different number and types of variables. A better approach could be to create a
Result
class. A compromise (to get going) could be to attach anxarr
object to each model such asmodel.results
and let the user access this xarr as they wish.Similarly, explorations could be stored in this manner. An
exploration
could haveexploration.results
and the user can again decide how to access that.Thoughts?
neurolib docker
What do you think of creating a docker file with neurolib preinstalled?
There would be a couple of advantages: very easy install, can be run anywhere (Windows without problems), etc...
I can imagine some people not so comfortable with setting up the whole python environment, but ok enough with docker.Base class for integration / whole-brain network?
How exactly should we implement "collection of nodes"?
I'd go with
Circuit(nodes=[list_of_nodes], SC=SC, delays=delays)
and thencircuit = Circuit(...); circuit.run(t=..., dt=..., sampling_dt=..., ...)
the
run()
method would actually run as follows:
- run would call some python library for SDEs, e.g.
sdeint
- each node within
Circuit
would need to implement its own dy/dt update- the circuit class would first resolve coupling variables, i.e. variables from nodes that are present as input in other nodes within the network
- then gather all dy/dt updates for all the nodes
Model base class for nodes
one instance for excitatory and inhibitory? As in ALN? or just one actual population, so the E/I node would be represented using two instances?
what attributes are useful?
what methods should a node implement?
other comments?
Speed up kuramato() function
Dear all,
I am planning to use
neurolib.utils.functions.kuromato
function to estimate Kuramoto order parameter but I faced the problem that running this function takes too long.Here is the time profiler report (it takes almost 10min to calculate Kuramoto order parameter for timeseries of 8 nodes and 20 seconds):
Timer unit: 1e-06 s Total time: 588.338 s File: /opt/anaconda3/envs/neurolib/lib/python3.7/site-packages/neurolib/utils/functions.py Function: kuramoto at line 9 Line # Hits Time Per Hit % Time Line Contents ============================================================== 9 def kuramoto(traces, dt=0.1, smoothing=0.0, peakrange=[0.1, 0.2]): 10 """ 11 Computes the Kuramoto order parameter of a timeseries which is a measure for synchrony. 12 Can smooth timeseries if there is noise. 13 Peaks are then detected using a peakfinder. From these peaks a phase is derived and then 14 the amount of phase synchrony (the Kuramoto order parameter) is computed. 15 16 :param traces: Multidimensional timeseries array 17 :type traces: numpy.ndarray 18 :param dt: Integration time step 19 :type dt: float 20 :param smoothing: Gaussian smoothing strength 21 :type smoothing: float, optional 22 :param peakrange: Width range of peaks for peak detection with `scipy.signal.find_peaks_cwt` 23 :type peakrange: list[float], length 2 24 25 :return: Timeseries of Kuramoto order paramter 26 :rtype: numpy.ndarray 27 """ 28 1 39.0 39.0 0.0 phases = [] 29 1 3.0 3.0 0.0 nTraces = len(traces) 30 9 12.0 1.3 0.0 for n in range(nTraces): 31 8 209572.0 26196.5 0.0 tList = np.dot(range(len(traces[n])), dt / 1000) 32 8 102.0 12.8 0.0 a = traces[n] 33 34 # find peaks 35 8 27.0 3.4 0.0 if smoothing > 0: 36 a = scipy.ndimage.filters.gaussian_filter(traces[n], smoothing) # smooth data 37 8 579281716.0 72410214.5 98.5 maximalist = scipy.signal.find_peaks_cwt(a, np.arange(peakrange[0], peakrange[1])) 38 8 243.0 30.4 0.0 maximalist = np.append(maximalist, len(traces[n]) - 1).astype(int) 39 40 8 9.0 1.1 0.0 if len(maximalist) > 1: 41 8 12.0 1.5 0.0 phases.append([]) 42 8 6.0 0.8 0.0 lastMax = 0 43 4384 3695.0 0.8 0.0 for m in maximalist: 44 1604368 1019132.0 0.6 0.2 for t in range(lastMax, m): 45 # compute instantaneous phase 46 1599992 1871774.0 1.2 0.3 phi = 2 * np.pi * float(t - lastMax) / float(m - lastMax) 47 1599992 1174814.0 0.7 0.2 phases[n].append(phi) 48 4376 2871.0 0.7 0.0 lastMax = m 49 8 6.0 0.8 0.0 phases[n].append(2 * np.pi) 50 else: 51 logging.warning("Kuramoto: No peaks found, returning 0.") 52 return 0 53 54 # determine kuramoto order paramter 55 1 1.0 1.0 0.0 kuramoto = [] 56 200001 139197.0 0.7 0.0 for t in range(len(tList)): 57 200000 130176.0 0.7 0.0 R = 1j * 0 58 1800000 1258227.0 0.7 0.2 for n in range(nTraces): 59 1600000 2636652.0 1.6 0.4 R += np.exp(1j * phases[n][t]) 60 200000 184074.0 0.9 0.0 R /= nTraces 61 200000 425820.0 2.1 0.1 kuramoto.append(np.absolute(R)) 62 63 1 1.0 1.0 0.0 return kuramotoHere is the full code to reproduce the line profiler output:
from neurolib.utils.loadData import Dataset from neurolib.models.aln import ALNModel from neurolib.utils.functions import kuramoto from line_profiler import LineProfiler ds = Dataset("gw") model = ALNModel(Cmat=ds.Cmat, Dmat=ds.Dmat) model.params['dt'] = 0.1 model.params['duration'] = 20 * 1000 # ms # add custom parameter for downsampling results # 10 ms sampling steps for saving data, should be multiple of dt model.params['save_dt'] = 10.0 model.params["tauA"] = 600.0 model.params["sigma_ou"] = 0.0 model.params["b"] = 20.0 model.params["Ke_gl"] = 300.0 model.params["signalV"] = 80.0 model.params["mue_ext_mean"] = 1.5 model.params["mui_ext_mean"] = 0.2 model.run() lp = LineProfiler() lp_wrapper = lp(kuramoto) kur = lp_wrapper(model.rates_exc[::10]) lp.print_stats()I am wondering whether using
scipy.signal.find_peaks_cwt
is essential. With the basic scipy functionscipy.signal.find_peaks
it works 63 times faster:Timer unit: 1e-06 s Total time: 9.29577 s File: <ipython-input-11-6b8850e688bf> Function: fast_kuramoto at line 2 Line # Hits Time Per Hit % Time Line Contents ============================================================== 2 def fast_kuramoto(traces, dt=0.1, smoothing=0.0): 3 """ 4 Computes the Kuramoto order parameter of a timeseries which is a measure for synchrony. 5 Can smooth timeseries if there is noise. 6 Peaks are then detected using a peakfinder. From these peaks a phase is derived and then 7 the amount of phase synchrony (the Kuramoto order parameter) is computed. 8 9 :param traces: Multidimensional timeseries array 10 :type traces: numpy.ndarray 11 :param dt: Integration time step 12 :type dt: float 13 :param smoothing: Gaussian smoothing strength 14 :type smoothing: float, optional 15 :param peakrange: Width range of peaks for peak detection with `scipy.signal.find_peaks_cwt` 16 :type peakrange: list[float], length 2 17 18 :return: Timeseries of Kuramoto order paramter 19 :rtype: numpy.ndarray 20 """ 21 1 3.0 3.0 0.0 phases = [] 22 1 2.0 2.0 0.0 nTraces = len(traces) 23 9 12.0 1.3 0.0 for n in range(nTraces): 24 8 189800.0 23725.0 2.0 tList = np.dot(range(len(traces[n])), dt / 1000) 25 8 99.0 12.4 0.0 a = traces[n] 26 27 # find peaks 28 8 27.0 3.4 0.0 if smoothing > 0: 29 a = scipy.ndimage.filters.gaussian_filter(traces[n], smoothing) # smooth data 30 8 27252.0 3406.5 0.3 maximalist = scipy.signal.find_peaks(a, distance=10, prominence=5)[0] 31 8 239.0 29.9 0.0 maximalist = np.append(maximalist, len(traces[n]) - 1).astype(int) 32 33 8 8.0 1.0 0.0 if len(maximalist) > 1: 34 8 11.0 1.4 0.0 phases.append([]) 35 8 6.0 0.8 0.0 lastMax = 0 36 4345 3584.0 0.8 0.0 for m in maximalist: 37 1604329 1052751.0 0.7 11.3 for t in range(lastMax, m): 38 # compute instantaneous phase 39 1599992 1912686.0 1.2 20.6 phi = 2 * np.pi * float(t - lastMax) / float(m - lastMax) 40 1599992 1197383.0 0.7 12.9 phases[n].append(phi) 41 4337 3018.0 0.7 0.0 lastMax = m 42 8 7.0 0.9 0.0 phases[n].append(2 * np.pi) 43 else: 44 logging.warning("Kuramoto: No peaks found, returning 0.") 45 return 0 46 47 # determine kuramoto order paramter 48 1 1.0 1.0 0.0 kuramoto = [] 49 200001 143749.0 0.7 1.5 for t in range(len(tList)): 50 200000 134629.0 0.7 1.4 R = 1j * 0 51 1800000 1307135.0 0.7 14.1 for n in range(nTraces): 52 1600000 2686423.0 1.7 28.9 R += np.exp(1j * phases[n][t]) 53 200000 190310.0 1.0 2.0 R /= nTraces 54 200000 446639.0 2.2 4.8 kuramoto.append(np.absolute(R)) 55 56 1 1.0 1.0 0.0 return kuramotoThe correlation between the results of the
fast_kuramoto
andkuramoto
is around 0.99. So the output is almost identical.Actually, 9 seconds is also too much, so I used
numba
to speed up the code a bit more. Here is the final version:from numba import jit def fast_numba_kuramoto(traces, smoothing=5.0, dt=0.1): nTraces, nTimes = traces.shape phases = np.empty_like(traces) for n in range(nTraces): a = traces[n] # find peaks if smoothing > 0: # smooth data a = scipy.ndimage.filters.gaussian_filter(traces[n], smoothing) maximalist = scipy.signal.find_peaks(a, distance=10, prominence=5)[0] maximalist = np.append(maximalist, len(traces[n])-1).astype(int) if len(maximalist) > 1: phases[n, :] = _estimate_phase(maximalist, nTimes) else: return 0 # determine kuramoto order paramter kuramoto = _estimate_r(nTraces, nTimes, phases) return kuramoto @jit(nopython=True) def _estimate_phase(maximalist, n_times): lastMax = 0 phases = np.empty((n_times), dtype=np.float64) n = 0 for m in maximalist: for t in range(lastMax, m): phi = 2 * np.pi * float(t - lastMax) / float(m - lastMax) phases[n] = phi n += 1 lastMax = m phases[-1] = 2 * np.pi return phases @jit(nopython=True) def _estimate_r(ntraces, times, phases): kuramoto = np.empty((times), dtype=np.float64) for t in range(times): R = 1j*0 for n in range(ntraces): R += np.exp(1j * phases[n, t]) R /= ntraces kuramoto[t] = np.absolute(R) return kuramotoThe final version takes around 63ms.
Looking forward to your comments.
Valerii
Fix Wilson-Cowan variable names
The
WCModel
's variables are calledx
andy
instead ofe
andi
.`aln` timeIntegration kHz --> Hz
aln
'stimeIntegration()
multiplies the rates with1e3
at every time step.For performance increase, this could be done once after
timeIntegration()
, using something like apostprocess
function.Add Neurolib to Open Neuroscience
Hello!
We are reaching out because we would love to have your project listed on Open Neuroscience, and also share information about this project:
Open Neuroscience is a community run project, where we are curating and highlighting open source projects related to neurosciences!
Briefly, we have a website where short descritptions about projects are listed, with links to the projects themselves, their authors, together with images and other links.
Once a new entry is made, we make a quick check for spam, and publish it.
Once published, we make people aware of the new entry by Twitter and a Facebook group.
To add information about their project, developers only need to fill out this form
In the form, people can add subfields and tags to their entries, so that projects are filterable and searchable on the website!
The reason why we have the form system is that it makes it open for everyone to contribute to the website and allows developers themselves to describe their projects!
Also, there are so many amazing projects coming out in Neurosciences that it would be impossible for us to keep track and log them all!
Open Neuroscience tech stack leverages open source tools as much as possible:
- The website is based on HUGO + Academic Theme
- Everything is hosted on github here
- We use plausible.io to see visit stats on the website. It respects visitors privacy, doesn't install cookies on your computer
- You can check our visitor stats here
Please get in touch if you have any questions or would like to collaborate!
Reference for the empirical data
Many thanks for sharing this,
Where can I find the reference of empirical data for FC? detail fo the experiment.Fix: MultiModel dot_dict params fails when connectivity or delay matrices are integer arrays
No idea why, I found it very randomly.. however, fix is easy: explicitly cast connectivity and delays into floats using numpy's
.astype(float)
func.kuramoto() is shorter than timeseries
kuramoto()
inneurolib.utils.functions
returns an array that is shorter than its input. This is probably because no padding is done before calculating the value. Needs to be fixed.Recommend Projects
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow
An Open Source Machine Learning Framework for Everyone
Django
The Web framework for perfectionists with deadlines.
Laravel
A PHP framework for web artisans
D3
Bring data to life with SVG, Canvas and HTML. ๐๐๐
Recommend Topics
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web
Some thing interesting about web. New door for the world.
server
A server is a program made to process requests and deliver data to clients.
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization
Some thing interesting about visualization, use data art
Game
Some thing interesting about game, make everyone happy.
Recommend Org
We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba
Alibaba Open Source for everyone
D3
Data-Driven Documents codes.
Tencent
China tencent open source team.