From 5123f6a26834f9e69ad446bd780c12c74b21ba86 Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Mon, 9 Sep 2019 20:25:41 +0200 Subject: [PATCH 01/10] Added constraints handling. Version 1. Not tested. --- tvb/simulator/integrators.py | 38 ++++++++++++++++++- tvb/simulator/models/base.py | 16 ++++++++ .../models/wong_wang_exc_io_inh_i.py | 36 +++++++----------- tvb/simulator/simulator.py | 14 +++++++ .../library/simulator/integrators_test.py | 24 ++++++++++++ 5 files changed, 105 insertions(+), 23 deletions(-) diff --git a/tvb/simulator/integrators.py b/tvb/simulator/integrators.py index 6c383d49b..d3e97d7aa 100644 --- a/tvb/simulator/integrators.py +++ b/tvb/simulator/integrators.py @@ -88,6 +88,17 @@ class Integrator(core.Type): it is consitent with Monitors using sample periods corresponding to powers of 2 from 128 to 4096Hz.""") + constraint_state_variable_indices = arrays.IntegerArray( + label="indices of the state variables to be constraint by the integrators " + "within the boundaries in the constraint values array", + default=None, + order=-1) + + constraint_state_variable_boundaries = arrays.FloatArray( + label = "The boundary values of the state variables within which they are constraint", + default = None, + order=-1) + clamped_state_variable_indices = arrays.IntegerArray( label = "indices of the state variables to be clamped by the integrators to the values in the clamped_values array", default = None, @@ -98,7 +109,6 @@ class Integrator(core.Type): default = None, order=-1) - def scheme(self, X, dfun, coupling, local_coupling, stimulus): """ The scheme of integrator should take a state and provide the next @@ -109,6 +119,18 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): msg = "Integrator is a base class; please use a suitable subclass." raise NotImplementedError(msg) + def constrain_state(self, X): + if self.constraint_state_variable_boundaries is not None: + for sv_ind, sv_bounds in \ + zip(self.constraint_state_variable_indices, + self.constraint_state_variable_boundaries): + temp = X[sv_ind] + if sv_bounds[0] is not None: + temp[temp < sv_bounds[0]] = sv_bounds[0] + if sv_bounds[1] is not None: + temp[temp > sv_bounds[1]] = sv_bounds[1] + X[sv_ind] = temp + def clamp_state(self, X): if self.clamped_state_variable_values is not None: X[self.clamped_state_variable_indices] = self.clamped_state_variable_values @@ -116,6 +138,7 @@ def clamp_state(self, X): def __str__(self): return simple_gen_astr(self, 'dt') + class IntegratorStochastic(Integrator): r""" The IntegratorStochastic class is a base class for the stochastic @@ -178,11 +201,13 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): #import pdb; pdb.set_trace() m_dx_tn = dfun(X, coupling, local_coupling) inter = X + self.dt * (m_dx_tn + stimulus) + self.constrain_state(inter) self.clamp_state(inter) dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0 X_next = X + dX + self.dt * stimulus + self.constrain_state(X_next) self.clamp_state(X_next) return X_next @@ -221,12 +246,15 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): noise *= noise_gfun inter = X + self.dt * m_dx_tn + noise + self.dt * stimulus + self.constrain_state(inter) self.clamp_state(inter) dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0 X_next = X + dX + noise + self.dt * stimulus + self.constrain_state(X_next) self.clamp_state(X_next) + return X_next @@ -254,6 +282,7 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): self.dX = dfun(X, coupling, local_coupling) X_next = X + self.dt * (self.dX + stimulus) + self.constrain_state(X_next) self.clamp_state(X_next) return X_next @@ -288,6 +317,7 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): dX = dfun(X, coupling, local_coupling) * self.dt noise_gfun = self.noise.gfun(X) X_next = X + dX + noise_gfun * noise + self.dt * stimulus + self.constrain_state(X_next) self.clamp_state(X_next) return X_next @@ -326,18 +356,22 @@ def scheme(self, X, dfun, coupling, local_coupling=0.0, stimulus=0.0): k1 = dfun(X, coupling, local_coupling) inter_k1 = X + dt2 * k1 + self.constrain_state(inter_k1) self.clamp_state(inter_k1) k2 = dfun(inter_k1, coupling, local_coupling) inter_k2 = X + dt2 * k2 + self.constrain_state(inter_k2) self.clamp_state(inter_k2) k3 = dfun(inter_k2, coupling, local_coupling) inter_k3 = X + dt * k3 + self.constrain_state(inter_k3) self.clamp_state(inter_k3) k4 = dfun(inter_k3, coupling, local_coupling) dX = dt6 * (k1 + 2.0 * k2 + 2.0 * k3 + k4) X_next = X + dX + self.dt * stimulus + self.constrain_state(X_next) self.clamp_state(X_next) return X_next @@ -401,6 +435,7 @@ class SciPyODE(SciPyODEBase): def scheme(self, X, dfun, coupling, local_coupling, stimulus): X_next = self._apply_ode(X, dfun, coupling, local_coupling, stimulus) + self.constrain_state(X_next) self.clamp_state(X_next) return X_next @@ -409,6 +444,7 @@ class SciPySDE(SciPyODEBase): def scheme(self, X, dfun, coupling, local_coupling, stimulus): X_next = self._apply_ode(X, dfun, coupling, local_coupling, stimulus) X_next += self.noise.gfun(X) * self.noise.generate(X.shape) + self.constrain_state(X_next) self.clamp_state(X_next) return X_next diff --git a/tvb/simulator/models/base.py b/tvb/simulator/models/base.py index f412d49bb..abcc7a222 100644 --- a/tvb/simulator/models/base.py +++ b/tvb/simulator/models/base.py @@ -60,6 +60,7 @@ class Model(core.Type): _nvar = None number_of_modes = 1 cvar = None + state_variable_constraint = None def _build_observer(self): template = ("def observe(state):\n" @@ -85,6 +86,21 @@ def configure(self): super(Model, self).configure() self.update_derived_parameters() self._build_observer() + # Make sure that if there are any state variable constraints, ... + if isinstance(self.state_variable_constraint, dict): + for sv, constraint in self. state_variable_constraint.items(): + try: + # ...the constraints correspond to model's state variables, + self.state_variables.index(sv) + except: + # TODO: Add the correct type of error and error message + raise + # and for every two sided constraint, the left boundary is lower than the right one + if constraint[0] is not None and constraint[1] is not None: + assert constraint[0] <= constraint[1] + elif self.state_variable_constraint is not None: + # TODO: Add here a warning or, even, error? + self.state_variable_constraint = None @property def nvar(self): diff --git a/tvb/simulator/models/wong_wang_exc_io_inh_i.py b/tvb/simulator/models/wong_wang_exc_io_inh_i.py index 35213444b..dd37ba531 100644 --- a/tvb/simulator/models/wong_wang_exc_io_inh_i.py +++ b/tvb/simulator/models/wong_wang_exc_io_inh_i.py @@ -44,30 +44,16 @@ def _numba_dfun(S, c, ae, be, de, ge, te, wp, we, jn, ai, bi, di, gi, ti, wi, ji cc = g[0]*jn[0]*c[0] - if S[0] < 0.0: - S_e = 0.0 # - S[0] # TODO: clarify the boundary to be reflective or saturated!!! - elif S[0] > 1.0: - S_e = 1.0 # - S[0] # TODO: clarify the boundary to be reflective or saturated!!! - else: - S_e = S[0] - - if S[1] < 0.0: - S_i = 0.0 # - S[1] TODO: clarify the boundary to be reflective or saturated!!! - elif S[1] > 1.0: - S_i = 1.0 # - S[1] TODO: clarify the boundary to be reflective or saturated!!! - else: - S_i = S[1] - - jnSe = jn[0]*S_e - x = wp[0]*jnSe - ji[0]*S_i + we[0]*io[0] + cc + jnSe = jn[0]*S[0] + x = wp[0]*jnSe - ji[0]*S[1] + we[0]*io[0] + cc x = ae[0]*x - be[0] h = x / (1 - numpy.exp(-de[0]*x)) - dx[0] = - (S_e / te[0]) + (1.0 - S_e) * h * ge[0] + dx[0] = - (S[0] / te[0]) + (1.0 - S[0]) * h * ge[0] - x = jnSe - S_i + wi[0]*io[0] + l[0]*cc + x = jnSe - S[1] + wi[0]*io[0] + l[0]*cc x = ai[0]*x - bi[0] h = x / (1 - numpy.exp(-di[0]*x)) - dx[1] = - (S_i / ti[0]) + h * gi[0] + dx[1] = - (S[1] / ti[0]) + h * gi[0] class ReducedWongWangExcIOInhI(ModelNumbaDfun): @@ -234,13 +220,21 @@ class ReducedWongWangExcIOInhI(ModelNumbaDfun): order=22 ) + # Used for phase-plane axis ranges and to bound random initial() conditions. + state_variable_constraint = basic.Dict( + label="State Variable constraints [lo, hi]", + default={"S_e": numpy.array([0.0, 1.0]), "S_i": numpy.array([0.0, 1.0])}, + doc="""The values for each state-variable should be set to encompass + the boundaries of the dynamic range of that state-variable. Set None for one-sided boundaries""", + order=23) + variables_of_interest = basic.Enumerate( label="Variables watched by Monitors", options=['S_e', 'S_i'], default=['S_e', 'S_i'], select_multiple=True, doc="""default state variables to be monitored""", - order=23) + order=24) state_variables = ['S_e', 'S_i'] _nvar = 2 @@ -266,8 +260,6 @@ def _numpy_dfun(self, state_variables, coupling, local_coupling=0.0): """ S = state_variables[:, :] - S[S < 0] = 0. - S[S > 1] = 1. c_0 = coupling[0, :] diff --git a/tvb/simulator/simulator.py b/tvb/simulator/simulator.py index a4f5e4711..72d67c200 100644 --- a/tvb/simulator/simulator.py +++ b/tvb/simulator/simulator.py @@ -233,6 +233,18 @@ def preconfigure(self): self.coupling.configure() self.model.configure() self.integrator.configure() + if isinstance(self.model.state_variable_constraint, dict): + indices = [] + boundaries = [] + for sv, sv_bounds in self.model.state_variable_constraint.items(): + indices.append(self.model.state_variables.index(sv)) + boundaries.append(sv_bounds) + sort_inds = numpy.sort(indices) + self.integrator.constraint_state_variable_indices = numpy.array(indices)[sort_inds] + self.integrator.constraint_state_variable_boundaries = numpy.array(boundaries)[sort_inds] + else: + self.integrator.constraint_state_variable_indices = None + self.integrator.constraint_state_variable_boundaries = None # monitors needs to be a list or tuple, even if there is only one... if not isinstance(self.monitors, (list, tuple)): self.monitors = [self.monitors] @@ -465,6 +477,8 @@ def _configure_history(self, initial_conditions): history[:ic_shape[0], :, :, :] = initial_conditions history = numpy.roll(history, shift, axis=0) self.current_step += ic_shape[0] - 1 + # TODO: confirm that this is doing what I think it is doing and that this is where this line should be placed! + numpy.apply_along_axis(self.integrator.constrain_state, 0, history) LOG.info('Final initial history shape is %r', history.shape) # create initial state from history self.current_state = history[self.current_step % self.horizon].copy() diff --git a/tvb/tests/library/simulator/integrators_test.py b/tvb/tests/library/simulator/integrators_test.py index 6cd7c7a53..db118584d 100644 --- a/tvb/tests/library/simulator/integrators_test.py +++ b/tvb/tests/library/simulator/integrators_test.py @@ -132,6 +132,30 @@ def test_scipy_dop853(self): def test_scipy_dopri5(self): self._scipy_scheme_tester('Dopri5') + def test_constrain(self): + vode = integrators.VODE( + constraint_state_variable_indices=numpy.r_[0, 1, 2, 3], + constraint_state_variable_values=numpy.array([[0.0, 1.0], [None, 1.0], [0.0, None], [None, None]]) + ) + x = numpy.ones((5, 4, 2)) + x[:, 0, ] = -x[:, 0, ] + x[:, 1, ] = 2 * x[:, 1, ] + x0 = numpy.array(x) + for i in range(10): + x = vode.scheme(x, self._dummy_dfun, 0.0, 0.0, 0.0) + for idx, val in zip(vode.constraint_state_variable_indices, vode.constraint_state_variable_indices): + if idx == 0: + assert numpy.all(x[idx] >= val[0]) + assert numpy.all(x[idx] <= val[1]) + elif idx == 1: + assert numpy.all(x[idx] <= val[1]) + assert numpy.allclose(x[idx, 0], x0[idx, 0]) + elif idx == 2: + assert numpy.all(x[idx] >= val[0]) + assert numpy.allclose(x[idx, 1], x0[idx, 1]) + else: + assert numpy.allclose(x[idx], x0[idx]) + def test_clamp(self): vode = integrators.VODE( clamped_state_variable_indices=numpy.r_[0, 3], From 7047ef17541c6dcf2e27fe4ab59cc96331c6d10b Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Tue, 10 Sep 2019 17:03:02 +0200 Subject: [PATCH 02/10] Added constraints checking back into the dfun of the RedWongWang to be independent of constraints checked in the integrator's scheme. --- .../models/wong_wang_exc_io_inh_i.py | 26 +++++++++++++++---- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/tvb/simulator/models/wong_wang_exc_io_inh_i.py b/tvb/simulator/models/wong_wang_exc_io_inh_i.py index dd37ba531..edd4075dc 100644 --- a/tvb/simulator/models/wong_wang_exc_io_inh_i.py +++ b/tvb/simulator/models/wong_wang_exc_io_inh_i.py @@ -44,16 +44,30 @@ def _numba_dfun(S, c, ae, be, de, ge, te, wp, we, jn, ai, bi, di, gi, ti, wi, ji cc = g[0]*jn[0]*c[0] - jnSe = jn[0]*S[0] - x = wp[0]*jnSe - ji[0]*S[1] + we[0]*io[0] + cc + if S[0] < 0.0: + S_e = 0.0 # - S[0] # TODO: clarify the boundary to be reflective or saturated!!! + elif S[0] > 1.0: + S_e = 1.0 # - S[0] # TODO: clarify the boundary to be reflective or saturated!!! + else: + S_e = S[0] + + if S[1] < 0.0: + S_i = 0.0 # - S[1] TODO: clarify the boundary to be reflective or saturated!!! + elif S[1] > 1.0: + S_i = 1.0 # - S[1] TODO: clarify the boundary to be reflective or saturated!!! + else: + S_i = S[1] + + jnSe = jn[0]*S_e + x = wp[0]*jnSe - ji[0]*S_i + we[0]*io[0] + cc x = ae[0]*x - be[0] h = x / (1 - numpy.exp(-de[0]*x)) - dx[0] = - (S[0] / te[0]) + (1.0 - S[0]) * h * ge[0] + dx[0] = - (S_e / te[0]) + (1.0 - S_e) * h * ge[0] - x = jnSe - S[1] + wi[0]*io[0] + l[0]*cc + x = jnSe - S_i + wi[0]*io[0] + l[0]*cc x = ai[0]*x - bi[0] h = x / (1 - numpy.exp(-di[0]*x)) - dx[1] = - (S[1] / ti[0]) + h * gi[0] + dx[1] = - (S_i / ti[0]) + h * gi[0] class ReducedWongWangExcIOInhI(ModelNumbaDfun): @@ -260,6 +274,8 @@ def _numpy_dfun(self, state_variables, coupling, local_coupling=0.0): """ S = state_variables[:, :] + S[S < 0] = 0. + S[S > 1] = 1. c_0 = coupling[0, :] From 7f153db4aff7ec5fa77d6118bfbb8128dcf381ac Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Tue, 10 Sep 2019 17:44:21 +0200 Subject: [PATCH 03/10] Debugged and tested changes to simulator --- tvb/simulator/simulator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tvb/simulator/simulator.py b/tvb/simulator/simulator.py index 72d67c200..7f55565e3 100644 --- a/tvb/simulator/simulator.py +++ b/tvb/simulator/simulator.py @@ -239,7 +239,7 @@ def preconfigure(self): for sv, sv_bounds in self.model.state_variable_constraint.items(): indices.append(self.model.state_variables.index(sv)) boundaries.append(sv_bounds) - sort_inds = numpy.sort(indices) + sort_inds = numpy.argsort(indices) self.integrator.constraint_state_variable_indices = numpy.array(indices)[sort_inds] self.integrator.constraint_state_variable_boundaries = numpy.array(boundaries)[sort_inds] else: @@ -478,7 +478,7 @@ def _configure_history(self, initial_conditions): history = numpy.roll(history, shift, axis=0) self.current_step += ic_shape[0] - 1 # TODO: confirm that this is doing what I think it is doing and that this is where this line should be placed! - numpy.apply_along_axis(self.integrator.constrain_state, 0, history) + self.integrator.constrain_state(numpy.swapaxes(history, 0, 1)) LOG.info('Final initial history shape is %r', history.shape) # create initial state from history self.current_state = history[self.current_step % self.horizon].copy() From b4500d934ec67db646394c1969f58b258d3858c2 Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Tue, 10 Sep 2019 17:45:40 +0200 Subject: [PATCH 04/10] Debugged and tested changes to simulator --- tvb/simulator/simulator.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tvb/simulator/simulator.py b/tvb/simulator/simulator.py index 7f55565e3..50c217f46 100644 --- a/tvb/simulator/simulator.py +++ b/tvb/simulator/simulator.py @@ -477,7 +477,6 @@ def _configure_history(self, initial_conditions): history[:ic_shape[0], :, :, :] = initial_conditions history = numpy.roll(history, shift, axis=0) self.current_step += ic_shape[0] - 1 - # TODO: confirm that this is doing what I think it is doing and that this is where this line should be placed! self.integrator.constrain_state(numpy.swapaxes(history, 0, 1)) LOG.info('Final initial history shape is %r', history.shape) # create initial state from history From c24fd7030a43cdaaae3714fef014c74aa319455c Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Tue, 10 Sep 2019 17:46:26 +0200 Subject: [PATCH 05/10] Changes in imports in tvb/basic/traits/__init__.py. --- tvb/basic/traits/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tvb/basic/traits/__init__.py b/tvb/basic/traits/__init__.py index 126741d0b..56a835e77 100644 --- a/tvb/basic/traits/__init__.py +++ b/tvb/basic/traits/__init__.py @@ -117,8 +117,8 @@ """ -import core -import traited_interface, traited_interface2 +from . import core +from . import traited_interface, traited_interface2 from tvb.basic.profile import TvbProfile # Add interfaces based on configured parameter on classes From 68cd43bda8aa89a4831442e75bfc661fbb4d1f59 Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Thu, 12 Sep 2019 20:07:23 +0200 Subject: [PATCH 06/10] Got the if statements outside state clamp and constrain methods in integrator --- tvb/simulator/integrators.py | 92 ++++++++++++++++++++++-------------- 1 file changed, 57 insertions(+), 35 deletions(-) diff --git a/tvb/simulator/integrators.py b/tvb/simulator/integrators.py index d3e97d7aa..6a4287ecc 100644 --- a/tvb/simulator/integrators.py +++ b/tvb/simulator/integrators.py @@ -120,20 +120,18 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): raise NotImplementedError(msg) def constrain_state(self, X): - if self.constraint_state_variable_boundaries is not None: - for sv_ind, sv_bounds in \ + for sv_ind, sv_bounds in \ zip(self.constraint_state_variable_indices, self.constraint_state_variable_boundaries): - temp = X[sv_ind] - if sv_bounds[0] is not None: - temp[temp < sv_bounds[0]] = sv_bounds[0] - if sv_bounds[1] is not None: - temp[temp > sv_bounds[1]] = sv_bounds[1] - X[sv_ind] = temp + temp = X[sv_ind] + if sv_bounds[0] is not None: + temp[temp < sv_bounds[0]] = sv_bounds[0] + if sv_bounds[1] is not None: + temp[temp > sv_bounds[1]] = sv_bounds[1] + X[sv_ind] = temp def clamp_state(self, X): - if self.clamped_state_variable_values is not None: - X[self.clamped_state_variable_indices] = self.clamped_state_variable_values + X[self.clamped_state_variable_indices] = self.clamped_state_variable_values def __str__(self): return simple_gen_astr(self, 'dt') @@ -200,15 +198,19 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): """ #import pdb; pdb.set_trace() m_dx_tn = dfun(X, coupling, local_coupling) - inter = X + self.dt * (m_dx_tn + stimulus) - self.constrain_state(inter) - self.clamp_state(inter) + inter = X + self.dt * (m_dx_tn + stimulus) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(inter) + if self.clamped_state_variable_values is not None: + self.clamp_state(inter) dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0 X_next = X + dX + self.dt * stimulus - self.constrain_state(X_next) - self.clamp_state(X_next) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(X_next) + if self.clamped_state_variable_values is not None: + self.clamp_state(X_next) return X_next @@ -246,14 +248,18 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): noise *= noise_gfun inter = X + self.dt * m_dx_tn + noise + self.dt * stimulus - self.constrain_state(inter) - self.clamp_state(inter) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(inter) + if self.clamped_state_variable_values is not None: + self.clamp_state(inter) dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0 X_next = X + dX + noise + self.dt * stimulus - self.constrain_state(X_next) - self.clamp_state(X_next) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(X_next) + if self.clamped_state_variable_values is not None: + self.clamp_state(X_next) return X_next @@ -282,8 +288,10 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): self.dX = dfun(X, coupling, local_coupling) X_next = X + self.dt * (self.dX + stimulus) - self.constrain_state(X_next) - self.clamp_state(X_next) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(X_next) + if self.clamped_state_variable_values is not None: + self.clamp_state(X_next) return X_next @@ -317,8 +325,10 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): dX = dfun(X, coupling, local_coupling) * self.dt noise_gfun = self.noise.gfun(X) X_next = X + dX + noise_gfun * noise + self.dt * stimulus - self.constrain_state(X_next) - self.clamp_state(X_next) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(X_next) + if self.clamped_state_variable_values is not None: + self.clamp_state(X_next) return X_next @@ -356,23 +366,31 @@ def scheme(self, X, dfun, coupling, local_coupling=0.0, stimulus=0.0): k1 = dfun(X, coupling, local_coupling) inter_k1 = X + dt2 * k1 - self.constrain_state(inter_k1) - self.clamp_state(inter_k1) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(inter_k1) + if self.clamped_state_variable_values is not None: + self.clamp_state(inter_k1) k2 = dfun(inter_k1, coupling, local_coupling) inter_k2 = X + dt2 * k2 - self.constrain_state(inter_k2) - self.clamp_state(inter_k2) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(inter_k2) + if self.clamped_state_variable_values is not None: + self.clamp_state(inter_k2) k3 = dfun(inter_k2, coupling, local_coupling) inter_k3 = X + dt * k3 - self.constrain_state(inter_k3) - self.clamp_state(inter_k3) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(inter_k3) + if self.clamped_state_variable_values is not None: + self.clamp_state(inter_k3) k4 = dfun(inter_k3, coupling, local_coupling) dX = dt6 * (k1 + 2.0 * k2 + 2.0 * k3 + k4) X_next = X + dX + self.dt * stimulus - self.constrain_state(X_next) - self.clamp_state(X_next) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(X_next) + if self.clamped_state_variable_values is not None: + self.clamp_state(X_next) return X_next @@ -435,8 +453,10 @@ class SciPyODE(SciPyODEBase): def scheme(self, X, dfun, coupling, local_coupling, stimulus): X_next = self._apply_ode(X, dfun, coupling, local_coupling, stimulus) - self.constrain_state(X_next) - self.clamp_state(X_next) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(X_next) + if self.clamped_state_variable_values is not None: + self.clamp_state(X_next) return X_next class SciPySDE(SciPyODEBase): @@ -444,8 +464,10 @@ class SciPySDE(SciPyODEBase): def scheme(self, X, dfun, coupling, local_coupling, stimulus): X_next = self._apply_ode(X, dfun, coupling, local_coupling, stimulus) X_next += self.noise.gfun(X) * self.noise.generate(X.shape) - self.constrain_state(X_next) - self.clamp_state(X_next) + if self.constraint_state_variable_boundaries is not None: + self.constrain_state(X_next) + if self.clamped_state_variable_values is not None: + self.clamp_state(X_next) return X_next class VODE(SciPyODE, Integrator): From d42999abb49aaa729310f1a5cf31ad02f83a2eca Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Thu, 12 Sep 2019 20:08:52 +0200 Subject: [PATCH 07/10] Added if statement outside state constrain of history --- tvb/simulator/simulator.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tvb/simulator/simulator.py b/tvb/simulator/simulator.py index 50c217f46..709c667bd 100644 --- a/tvb/simulator/simulator.py +++ b/tvb/simulator/simulator.py @@ -477,7 +477,8 @@ def _configure_history(self, initial_conditions): history[:ic_shape[0], :, :, :] = initial_conditions history = numpy.roll(history, shift, axis=0) self.current_step += ic_shape[0] - 1 - self.integrator.constrain_state(numpy.swapaxes(history, 0, 1)) + if self.constraint_state_variable_boundaries is not None: + self.integrator.constrain_state(numpy.swapaxes(history, 0, 1)) LOG.info('Final initial history shape is %r', history.shape) # create initial state from history self.current_state = history[self.current_step % self.horizon].copy() From 92ed4f38bcee337646bef7109c9cbbbecc9d5b49 Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Thu, 12 Sep 2019 20:30:33 +0200 Subject: [PATCH 08/10] Removed TODOs on constraints from red wong wang model --- tvb/simulator/models/wong_wang_exc_io_inh_i.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tvb/simulator/models/wong_wang_exc_io_inh_i.py b/tvb/simulator/models/wong_wang_exc_io_inh_i.py index edd4075dc..3c95458f6 100644 --- a/tvb/simulator/models/wong_wang_exc_io_inh_i.py +++ b/tvb/simulator/models/wong_wang_exc_io_inh_i.py @@ -45,16 +45,16 @@ def _numba_dfun(S, c, ae, be, de, ge, te, wp, we, jn, ai, bi, di, gi, ti, wi, ji cc = g[0]*jn[0]*c[0] if S[0] < 0.0: - S_e = 0.0 # - S[0] # TODO: clarify the boundary to be reflective or saturated!!! + S_e = 0.0 elif S[0] > 1.0: - S_e = 1.0 # - S[0] # TODO: clarify the boundary to be reflective or saturated!!! + S_e = 1.0 else: S_e = S[0] if S[1] < 0.0: - S_i = 0.0 # - S[1] TODO: clarify the boundary to be reflective or saturated!!! + S_i = 0.0 elif S[1] > 1.0: - S_i = 1.0 # - S[1] TODO: clarify the boundary to be reflective or saturated!!! + S_i = 1.0 else: S_i = S[1] From 823b8eb8de9bc52c755e19d1de65ef790f0a0443 Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Thu, 12 Sep 2019 20:31:14 +0200 Subject: [PATCH 09/10] Converted reflective to saturated constraints for wong-wang _numba_dfun --- tvb/simulator/models/wong_wang.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/tvb/simulator/models/wong_wang.py b/tvb/simulator/models/wong_wang.py index a55112f97..a5c8c8eb4 100644 --- a/tvb/simulator/models/wong_wang.py +++ b/tvb/simulator/models/wong_wang.py @@ -39,13 +39,15 @@ def _numba_dfun(S, c, a, b, d, g, ts, w, j, io, dx): "Gufunc for reduced Wong-Wang model equations." if S[0] < 0.0: - dx[0] = 0.0 - S[0] + S_constraint = 0.0 elif S[0] > 1.0: - dx[0] = 1.0 - S[0] + S_constraint = 1.0 else: - x = w[0]*j[0]*S[0] + io[0] + j[0]*c[0] - h = (a[0]*x - b[0]) / (1 - numpy.exp(-d[0]*(a[0]*x - b[0]))) - dx[0] = - (S[0] / ts[0]) + (1.0 - S[0]) * h * g[0] + S_constraint = S[0] + + x = w[0]*j[0]*S_constraint + io[0] + j[0]*c[0] + h = (a[0]*x - b[0]) / (1 - numpy.exp(-d[0]*(a[0]*x - b[0]))) + dx[0] = - (S_constraint / ts[0]) + (1.0 - S_constraint) * h * g[0] class ReducedWongWang(ModelNumbaDfun): @@ -144,13 +146,21 @@ class ReducedWongWang(ModelNumbaDfun): order=9 ) + # Used for phase-plane axis ranges and to bound random initial() conditions. + state_variable_constraint = basic.Dict( + label="State Variable constraints [lo, hi]", + default={"S": numpy.array([0.0, 1.0])}, + doc="""The values for each state-variable should be set to encompass + the boundaries of the dynamic range of that state-variable. Set None for one-sided boundaries""", + order=10) + variables_of_interest = basic.Enumerate( label="Variables watched by Monitors", options=["S"], default=["S"], select_multiple=True, doc="""default state variables to be monitored""", - order=10) + order=11) state_variables = ['S'] _nvar = 1 From 33d2fdf6865dd2dc3ce983bb86319e4abd0094b6 Mon Sep 17 00:00:00 2001 From: Dionysios Perdikis Date: Fri, 13 Sep 2019 21:30:55 +0200 Subject: [PATCH 10/10] Changed names form constraints to bound... Removed boundary check inside wong-wang models Removed temp from bound_state method. --- tvb/simulator/integrators.py | 88 +++++++++---------- tvb/simulator/models/base.py | 18 ++-- tvb/simulator/models/wong_wang.py | 19 ++-- .../models/wong_wang_exc_io_inh_i.py | 31 ++----- tvb/simulator/simulator.py | 4 +- .../library/simulator/integrators_test.py | 8 +- 6 files changed, 71 insertions(+), 97 deletions(-) diff --git a/tvb/simulator/integrators.py b/tvb/simulator/integrators.py index 6a4287ecc..1799e223e 100644 --- a/tvb/simulator/integrators.py +++ b/tvb/simulator/integrators.py @@ -75,11 +75,11 @@ class Integrator(core.Type): _base_classes = ['Integrator', 'IntegratorStochastic', 'RungeKutta4thOrderDeterministic'] dt = basic.Float( - label = "Integration-step size (ms)", - default = 0.01220703125, #0.015625, + label="Integration-step size (ms)", + default=0.01220703125, #0.015625, #range = basic.Range(lo= 0.0048828125, hi=0.244140625, step= 0.1, base=2.) - required = True, - doc = """The step size used by the integration routine in ms. This + required=True, + doc="""The step size used by the integration routine in ms. This should be chosen to be small enough for the integration to be numerically stable. It is also necessary to consider the desired sample period of the Monitors, as they are restricted to being integral @@ -88,25 +88,25 @@ class Integrator(core.Type): it is consitent with Monitors using sample periods corresponding to powers of 2 from 128 to 4096Hz.""") - constraint_state_variable_indices = arrays.IntegerArray( - label="indices of the state variables to be constraint by the integrators " - "within the boundaries in the constraint values array", + bounded_state_variable_indices = arrays.IntegerArray( + label="indices of the state variables to be bounded by the integrators " + "within the boundaries in the boundaries' values array", default=None, order=-1) - constraint_state_variable_boundaries = arrays.FloatArray( - label = "The boundary values of the state variables within which they are constraint", - default = None, + state_variable_boundaries = arrays.FloatArray( + label="The boundary values of the state variables", + default=None, order=-1) clamped_state_variable_indices = arrays.IntegerArray( - label = "indices of the state variables to be clamped by the integrators to the values in the clamped_values array", - default = None, + label="indices of the state variables to be clamped by the integrators to the values in the clamped_values array", + default=None, order=-1) clamped_state_variable_values = arrays.FloatArray( - label = "The values of the state variables which are clamped ", - default = None, + label="The values of the state variables which are clamped ", + default=None, order=-1) def scheme(self, X, dfun, coupling, local_coupling, stimulus): @@ -119,16 +119,14 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): msg = "Integrator is a base class; please use a suitable subclass." raise NotImplementedError(msg) - def constrain_state(self, X): + def bound_state(self, X): for sv_ind, sv_bounds in \ - zip(self.constraint_state_variable_indices, - self.constraint_state_variable_boundaries): - temp = X[sv_ind] + zip(self.bounded_state_variable_indices, + self.state_variable_boundaries): if sv_bounds[0] is not None: - temp[temp < sv_bounds[0]] = sv_bounds[0] + X[sv_ind][X[sv_ind] < sv_bounds[0]] = sv_bounds[0] if sv_bounds[1] is not None: - temp[temp > sv_bounds[1]] = sv_bounds[1] - X[sv_ind] = temp + X[sv_ind][X[sv_ind] > sv_bounds[1]] = sv_bounds[1] def clamp_state(self, X): X[self.clamped_state_variable_indices] = self.clamped_state_variable_values @@ -199,16 +197,16 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): #import pdb; pdb.set_trace() m_dx_tn = dfun(X, coupling, local_coupling) inter = X + self.dt * (m_dx_tn + stimulus) - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(inter) + if self.state_variable_boundaries is not None: + self.bound_state(inter) if self.clamped_state_variable_values is not None: self.clamp_state(inter) dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0 X_next = X + dX + self.dt * stimulus - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(X_next) + if self.state_variable_boundaries is not None: + self.bound_state(X_next) if self.clamped_state_variable_values is not None: self.clamp_state(X_next) return X_next @@ -248,16 +246,16 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): noise *= noise_gfun inter = X + self.dt * m_dx_tn + noise + self.dt * stimulus - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(inter) + if self.state_variable_boundaries is not None: + self.bound_state(inter) if self.clamped_state_variable_values is not None: self.clamp_state(inter) dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0 X_next = X + dX + noise + self.dt * stimulus - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(X_next) + if self.state_variable_boundaries is not None: + self.bound_state(X_next) if self.clamped_state_variable_values is not None: self.clamp_state(X_next) @@ -288,8 +286,8 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): self.dX = dfun(X, coupling, local_coupling) X_next = X + self.dt * (self.dX + stimulus) - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(X_next) + if self.state_variable_boundaries is not None: + self.bound_state(X_next) if self.clamped_state_variable_values is not None: self.clamp_state(X_next) return X_next @@ -325,8 +323,8 @@ def scheme(self, X, dfun, coupling, local_coupling, stimulus): dX = dfun(X, coupling, local_coupling) * self.dt noise_gfun = self.noise.gfun(X) X_next = X + dX + noise_gfun * noise + self.dt * stimulus - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(X_next) + if self.state_variable_boundaries is not None: + self.bound_state(X_next) if self.clamped_state_variable_values is not None: self.clamp_state(X_next) return X_next @@ -366,20 +364,20 @@ def scheme(self, X, dfun, coupling, local_coupling=0.0, stimulus=0.0): k1 = dfun(X, coupling, local_coupling) inter_k1 = X + dt2 * k1 - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(inter_k1) + if self.state_variable_boundaries is not None: + self.bound_state(inter_k1) if self.clamped_state_variable_values is not None: self.clamp_state(inter_k1) k2 = dfun(inter_k1, coupling, local_coupling) inter_k2 = X + dt2 * k2 - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(inter_k2) + if self.state_variable_boundaries is not None: + self.bound_state(inter_k2) if self.clamped_state_variable_values is not None: self.clamp_state(inter_k2) k3 = dfun(inter_k2, coupling, local_coupling) inter_k3 = X + dt * k3 - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(inter_k3) + if self.state_variable_boundaries is not None: + self.bound_state(inter_k3) if self.clamped_state_variable_values is not None: self.clamp_state(inter_k3) k4 = dfun(inter_k3, coupling, local_coupling) @@ -387,8 +385,8 @@ def scheme(self, X, dfun, coupling, local_coupling=0.0, stimulus=0.0): dX = dt6 * (k1 + 2.0 * k2 + 2.0 * k3 + k4) X_next = X + dX + self.dt * stimulus - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(X_next) + if self.state_variable_boundaries is not None: + self.bound_state(X_next) if self.clamped_state_variable_values is not None: self.clamp_state(X_next) return X_next @@ -453,8 +451,8 @@ class SciPyODE(SciPyODEBase): def scheme(self, X, dfun, coupling, local_coupling, stimulus): X_next = self._apply_ode(X, dfun, coupling, local_coupling, stimulus) - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(X_next) + if self.state_variable_boundaries is not None: + self.bound_state(X_next) if self.clamped_state_variable_values is not None: self.clamp_state(X_next) return X_next @@ -464,8 +462,8 @@ class SciPySDE(SciPyODEBase): def scheme(self, X, dfun, coupling, local_coupling, stimulus): X_next = self._apply_ode(X, dfun, coupling, local_coupling, stimulus) X_next += self.noise.gfun(X) * self.noise.generate(X.shape) - if self.constraint_state_variable_boundaries is not None: - self.constrain_state(X_next) + if self.state_variable_boundaries is not None: + self.bound_state(X_next) if self.clamped_state_variable_values is not None: self.clamp_state(X_next) return X_next diff --git a/tvb/simulator/models/base.py b/tvb/simulator/models/base.py index abcc7a222..25540d053 100644 --- a/tvb/simulator/models/base.py +++ b/tvb/simulator/models/base.py @@ -60,7 +60,7 @@ class Model(core.Type): _nvar = None number_of_modes = 1 cvar = None - state_variable_constraint = None + state_variable_boundaries = None def _build_observer(self): template = ("def observe(state):\n" @@ -86,21 +86,21 @@ def configure(self): super(Model, self).configure() self.update_derived_parameters() self._build_observer() - # Make sure that if there are any state variable constraints, ... - if isinstance(self.state_variable_constraint, dict): - for sv, constraint in self. state_variable_constraint.items(): + # Make sure that if there are any state variable boundaries, ... + if isinstance(self.state_variable_boundaries, dict): + for sv, bounds in self. state_variable_boundaries.items(): try: - # ...the constraints correspond to model's state variables, + # ...the boundaries correspond to model's state variables, self.state_variables.index(sv) except: # TODO: Add the correct type of error and error message raise # and for every two sided constraint, the left boundary is lower than the right one - if constraint[0] is not None and constraint[1] is not None: - assert constraint[0] <= constraint[1] - elif self.state_variable_constraint is not None: + if bounds[0] is not None and bounds[1] is not None: + assert bounds[0] <= bounds[1] + elif self.state_variable_boundaries is not None: # TODO: Add here a warning or, even, error? - self.state_variable_constraint = None + self.state_variable_boundaries = None @property def nvar(self): diff --git a/tvb/simulator/models/wong_wang.py b/tvb/simulator/models/wong_wang.py index a5c8c8eb4..edd4d5b69 100644 --- a/tvb/simulator/models/wong_wang.py +++ b/tvb/simulator/models/wong_wang.py @@ -37,17 +37,9 @@ @guvectorize([(float64[:],)*11], '(n),(m)' + ',()'*8 + '->(n)', nopython=True) def _numba_dfun(S, c, a, b, d, g, ts, w, j, io, dx): "Gufunc for reduced Wong-Wang model equations." - - if S[0] < 0.0: - S_constraint = 0.0 - elif S[0] > 1.0: - S_constraint = 1.0 - else: - S_constraint = S[0] - - x = w[0]*j[0]*S_constraint + io[0] + j[0]*c[0] + x = w[0]*j[0]*S[0] + io[0] + j[0]*c[0] h = (a[0]*x - b[0]) / (1 - numpy.exp(-d[0]*(a[0]*x - b[0]))) - dx[0] = - (S_constraint / ts[0]) + (1.0 - S_constraint) * h * g[0] + dx[0] = - (S[0] / ts[0]) + (1.0 - S[0]) * h * g[0] class ReducedWongWang(ModelNumbaDfun): @@ -147,8 +139,8 @@ class ReducedWongWang(ModelNumbaDfun): ) # Used for phase-plane axis ranges and to bound random initial() conditions. - state_variable_constraint = basic.Dict( - label="State Variable constraints [lo, hi]", + state_variable_boundaries = basic.Dict( + label="State Variable boundaries [lo, hi]", default={"S": numpy.array([0.0, 1.0])}, doc="""The values for each state-variable should be set to encompass the boundaries of the dynamic range of that state-variable. Set None for one-sided boundaries""", @@ -182,8 +174,7 @@ def _numpy_dfun(self, state_variables, coupling, local_coupling=0.0): """ S = state_variables[0, :] - S[S<0] = 0. - S[S>1] = 1. + c_0 = coupling[0, :] diff --git a/tvb/simulator/models/wong_wang_exc_io_inh_i.py b/tvb/simulator/models/wong_wang_exc_io_inh_i.py index 3c95458f6..7def50e39 100644 --- a/tvb/simulator/models/wong_wang_exc_io_inh_i.py +++ b/tvb/simulator/models/wong_wang_exc_io_inh_i.py @@ -44,30 +44,17 @@ def _numba_dfun(S, c, ae, be, de, ge, te, wp, we, jn, ai, bi, di, gi, ti, wi, ji cc = g[0]*jn[0]*c[0] - if S[0] < 0.0: - S_e = 0.0 - elif S[0] > 1.0: - S_e = 1.0 - else: - S_e = S[0] - - if S[1] < 0.0: - S_i = 0.0 - elif S[1] > 1.0: - S_i = 1.0 - else: - S_i = S[1] - - jnSe = jn[0]*S_e - x = wp[0]*jnSe - ji[0]*S_i + we[0]*io[0] + cc + jnSe = jn[0]*S[0] + + x = wp[0]*jnSe - ji[0]*S[1] + we[0]*io[0] + cc x = ae[0]*x - be[0] h = x / (1 - numpy.exp(-de[0]*x)) - dx[0] = - (S_e / te[0]) + (1.0 - S_e) * h * ge[0] + dx[0] = - (S[0] / te[0]) + (1.0 - S[0]) * h * ge[0] - x = jnSe - S_i + wi[0]*io[0] + l[0]*cc + x = jnSe - S[1] + wi[0]*io[0] + l[0]*cc x = ai[0]*x - bi[0] h = x / (1 - numpy.exp(-di[0]*x)) - dx[1] = - (S_i / ti[0]) + h * gi[0] + dx[1] = - (S[1] / ti[0]) + h * gi[0] class ReducedWongWangExcIOInhI(ModelNumbaDfun): @@ -235,8 +222,8 @@ class ReducedWongWangExcIOInhI(ModelNumbaDfun): ) # Used for phase-plane axis ranges and to bound random initial() conditions. - state_variable_constraint = basic.Dict( - label="State Variable constraints [lo, hi]", + state_variable_boundaries = basic.Dict( + label="State Variable boundaries [lo, hi]", default={"S_e": numpy.array([0.0, 1.0]), "S_i": numpy.array([0.0, 1.0])}, doc="""The values for each state-variable should be set to encompass the boundaries of the dynamic range of that state-variable. Set None for one-sided boundaries""", @@ -274,8 +261,6 @@ def _numpy_dfun(self, state_variables, coupling, local_coupling=0.0): """ S = state_variables[:, :] - S[S < 0] = 0. - S[S > 1] = 1. c_0 = coupling[0, :] diff --git a/tvb/simulator/simulator.py b/tvb/simulator/simulator.py index 709c667bd..2b75f3501 100644 --- a/tvb/simulator/simulator.py +++ b/tvb/simulator/simulator.py @@ -477,8 +477,8 @@ def _configure_history(self, initial_conditions): history[:ic_shape[0], :, :, :] = initial_conditions history = numpy.roll(history, shift, axis=0) self.current_step += ic_shape[0] - 1 - if self.constraint_state_variable_boundaries is not None: - self.integrator.constrain_state(numpy.swapaxes(history, 0, 1)) + if self.integrator.state_variable_boundaries is not None: + self.integrator.bound_state(numpy.swapaxes(history, 0, 1)) LOG.info('Final initial history shape is %r', history.shape) # create initial state from history self.current_state = history[self.current_step % self.horizon].copy() diff --git a/tvb/tests/library/simulator/integrators_test.py b/tvb/tests/library/simulator/integrators_test.py index db118584d..70ddb6c34 100644 --- a/tvb/tests/library/simulator/integrators_test.py +++ b/tvb/tests/library/simulator/integrators_test.py @@ -132,10 +132,10 @@ def test_scipy_dop853(self): def test_scipy_dopri5(self): self._scipy_scheme_tester('Dopri5') - def test_constrain(self): + def test_bound(self): vode = integrators.VODE( - constraint_state_variable_indices=numpy.r_[0, 1, 2, 3], - constraint_state_variable_values=numpy.array([[0.0, 1.0], [None, 1.0], [0.0, None], [None, None]]) + bounded_state_variable_indices=numpy.r_[0, 1, 2, 3], + state_variable_boundaries=numpy.array([[0.0, 1.0], [None, 1.0], [0.0, None], [None, None]]) ) x = numpy.ones((5, 4, 2)) x[:, 0, ] = -x[:, 0, ] @@ -143,7 +143,7 @@ def test_constrain(self): x0 = numpy.array(x) for i in range(10): x = vode.scheme(x, self._dummy_dfun, 0.0, 0.0, 0.0) - for idx, val in zip(vode.constraint_state_variable_indices, vode.constraint_state_variable_indices): + for idx, val in zip(vode.bounded_state_variable_indices, vode.state_variable_boundaries): if idx == 0: assert numpy.all(x[idx] >= val[0]) assert numpy.all(x[idx] <= val[1])