Merge pull request #116 from brian-team/refractoriness_always_in_thre…

…shold.py

Refractoriness always in threshold.py

brian2/groups/neurongroup.py

@@ -2,22 +2,20 @@
 This model defines the `NeuronGroup`, the core of most simulations.
 '''
 import numpy as np
-from numpy import array
 import sympy
 
 from brian2.equations.equations import (Equations, DIFFERENTIAL_EQUATION,
                                         STATIC_EQUATION, PARAMETER)
 from brian2.equations.refractory import add_refractoriness
 from brian2.stateupdaters.base import StateUpdateMethod
-from brian2.memory import allocate_array
+from brian2.devices.device import get_device
 from brian2.core.preferences import brian_prefs
 from brian2.core.base import BrianObject
 from brian2.core.namespace import create_namespace
-from brian2.core.variables import (Variable, AttributeVariable, ArrayVariable,
-                                    StochasticVariable, Subexpression)
+from brian2.core.variables import (ArrayVariable, StochasticVariable,
+                                   Subexpression)
 from brian2.core.spikesource import SpikeSource
 from brian2.core.scheduler import Scheduler
-from brian2.devices.device import get_device
 from brian2.parsing.expressions import (parse_expression_unit,
                                         is_boolean_expression)
 from brian2.utils.logger import get_logger
@@ -32,6 +30,37 @@
 logger = get_logger(__name__)
 
 
+def get_refractory_code(group):
+    ref = group._refractory
+    if ref is None:
+        # No refractoriness
+        abstract_code = ''
+    elif isinstance(ref, Quantity):
+        abstract_code = 'not_refractory = 1*((t - lastspike) > %f)\n' % ref
+    else:
+        namespace = group.namespace
+        unit = parse_expression_unit(str(ref), namespace, group.variables)
+        if have_same_dimensions(unit, second):
+            abstract_code = 'not_refractory = 1*((t - lastspike) > %s)\n' % ref
+        elif have_same_dimensions(unit, Unit(1)):
+            if not is_boolean_expression(str(ref), namespace,
+                                         group.variables):
+                raise TypeError(('Refractory expression is dimensionless '
+                                 'but not a boolean value. It needs to '
+                                 'either evaluate to a timespan or to a '
+                                 'boolean value.'))
+            # boolean condition
+            # we have to be a bit careful here, we can't just use the given
+            # condition as it is, because we only want to *leave*
+            # refractoriness, based on the condition
+            abstract_code = 'not_refractory = 1*(not_refractory or not (%s))\n' % ref
+        else:
+            raise TypeError(('Refractory expression has to evaluate to a '
+                             'timespan or a boolean value, expression'
+                             '"%s" has units %s instead') % (ref, unit))
+    return abstract_code
+
+
 class StateUpdater(GroupCodeRunner):
     '''
     The `GroupCodeRunner` that updates the state variables of a `NeuronGroup`
@@ -62,33 +91,7 @@ def __init__(self, group, method):
     def update_abstract_code(self):
 
         # Update the not_refractory variable for the refractory period mechanism
-        ref = self.group._refractory
-        if ref is None:
-            # No refractoriness
-            self.abstract_code = ''
-        elif isinstance(ref, Quantity):
-            self.abstract_code = 'not_refractory = (t - lastspike) > %f\n' % ref
-        else:
-            namespace = self.group.namespace
-            unit = parse_expression_unit(str(ref), namespace, self.group.variables)
-            if have_same_dimensions(unit, second):
-                self.abstract_code = 'not_refractory = (t - lastspike) > %s\n' % ref
-            elif have_same_dimensions(unit, Unit(1)):
-                if not is_boolean_expression(str(ref), namespace,
-                                             self.group.variables):
-                    raise TypeError(('Refractory expression is dimensionless '
-                                     'but not a boolean value. It needs to '
-                                     'either evaluate to a timespan or to a '
-                                     'boolean value.'))
-                # boolean condition
-                # we have to be a bit careful here, we can't just use the given
-                # condition as it is, because we only want to *leave*
-                # refractoriness, based on the condition
-                self.abstract_code = 'not_refractory = not_refractory or not (%s)\n' % ref
-            else:
-                raise TypeError(('Refractory expression has to evaluate to a '
-                                 'timespan or a boolean value, expression'
-                                 '"%s" has units %s instead') % (ref, unit))
+        self.abstract_code = get_refractory_code(self.group)
 
         self.abstract_code += self.method(self.group.equations,
                                           self.group.variables)
@@ -118,9 +121,8 @@ def __init__(self, group):
         self.update_abstract_code()
         check_code_units(self.abstract_code, self.group, ignore_keyerrors=True)
 
-
     def update_abstract_code(self):
-        self.abstract_code = '_cond = ' + self.group.threshold
+        self.abstract_code = '_cond = (%s) and not_refractory' % self.group.threshold
 
 
 class Resetter(GroupCodeRunner):

brian2/groups/poissongroup.py

@@ -3,15 +3,15 @@
 from brian2.core.base import BrianObject
 from brian2.core.namespace import create_namespace
 from brian2.core.spikesource import SpikeSource
-from brian2.core.scheduler import Scheduler
 from brian2.core.variables import ArrayVariable
 from brian2.devices.device import get_device
+from brian2.equations import Equations
 from brian2.units.fundamentalunits import check_units, Unit
 from brian2.units.allunits import second
 from brian2.units.stdunits import Hz
 
 from .group import Group
-from .neurongroup import Thresholder
+from .neurongroup import Thresholder, StateUpdater
 
 __all__ = ['PoissonGroup']
 
@@ -62,7 +62,8 @@ def __init__(self, N, rates, clock=None, name='poissongroup*',
         # for more complex use cases.
 
         #: The array storing the refractoriness information (not used, currently)
-        self._not_refractory = get_device().array(self, '_not_refractory', N, 1, dtype=np.bool)
+        self._not_refractory = get_device().array(self, '_not_refractory', N, 1,
+                                                  dtype=np.bool)
         self._lastspike = get_device().array(self, '_lastspike', N, 1)
 
         self.variables = Group._create_variables(self)
@@ -88,6 +89,12 @@ def __init__(self, N, rates, clock=None, name='poissongroup*',
         self.thresholder = Thresholder(self)
         self.contained_objects.append(self.thresholder)
 
+        # This is quite inefficient, we need a state updater to reset
+        # not_refractory after every time step
+        self.equations = Equations([])
+        self._refractory = False
+        self.state_updater = StateUpdater(self, method='independent')
+        self.contained_objects.append(self.state_updater)
         Group.__init__(self)
 
     @property

brian2/tests/test_refractory.py

@@ -3,6 +3,7 @@
 
 from brian2 import *
 from brian2.equations.refractory import add_refractoriness
+
 
 # We can only test C++ if weave is availabe
 try:
@@ -72,14 +73,14 @@ def test_refractoriness_threshold():
             dv/dt = 200*Hz : 1
             ref : second
             ref_no_unit : 1
-            ''', threshold='not_refractory and (v > 1)',
+            ''', threshold='v > 1',
                             reset='v=0', refractory=ref_time,
                             codeobj_class=codeobj_class)
             G.ref = 10*ms
             G.ref_no_unit = 10
             # The neuron should spike after 5ms but then not spike for the next
-            # 10ms. The state variable should continue to integrate so there
-            # should be a spike after 15ms
+            # 10ms. The state variable should continue to integrate so there should
+            # be a spike after 15ms
             spike_mon = SpikeMonitor(G)
             net = Network(G, spike_mon)
             net.run(16*ms)

docs_sphinx/user/refractoriness.rst

@@ -87,12 +87,14 @@ stop being updated during refractoriness can be marked with the
 In the above model, the ``v`` variable is clamped at 0 for 2ms after a spike but
 the adaptation variable ``w`` continues to update during this time.
 
-The second behaviour (ignore threshold crossings) can be implemented by
-including the special variable ``not_refractory`` in the threshold condition.
+The second behaviour (ignore threshold crossings) is already implemented
+automatically -- the threshold condition can only evaluate to ``True`` if
+the ``not_refractory`` variable is ``True``. This is achieved by automatically
+appending ``and not_refractory`` to the threshold condition.
 In the following model, the variable ``v`` continues to update during the
 refractory period but it does not elicit a spike if it crosses the threshold::
 
     G = NeuronGroup(N, 'dv/dt = -v / tau : 1',
-                    threshold='not_refractory and (v > 1)', reset='v=0',
+                    threshold='v > 1', reset='v=0',
                     refractory=2*ms)