Non-contiguous subgroups

brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx

@@ -1,31 +1,50 @@
-{# USES_VARIABLES { N, t, rate, _clock_t, _clock_dt, _spikespace,
-                    _num_source_neurons, _source_start, _source_stop } #}
+{# USES_VARIABLES { N, t, rate, _clock_t, _clock_dt, _spikespace} #}
 {% extends 'common_group.pyx' %}
 
 {% block maincode %}
 
     cdef size_t _num_spikes = {{_spikespace}}[_num{{_spikespace}}-1]
-
+    {% if subgroup and not contiguous %}
+    # We use the same data structure as for the eventspace to store the
+    # "filtered" events, i.e. the events that are indexed in the subgroup
+    cdef int[{{source_N}} + 1] _filtered_events
+    cdef size_t _source_index_counter = 0
+    _filtered_events[{{source_N}}] = 0
+    {% endif %}
+    {% if subgroup %}
     # For subgroups, we do not want to record all spikes
+    {% if contiguous %}
     # We assume that spikes are ordered
-    cdef int _start_idx = -1
-    cdef int _end_idx = -1
-    cdef size_t _j
+    _start_idx = _num_spikes
+    _end_idx = _num_spikes
     for _j in range(_num_spikes):
         _idx = {{_spikespace}}[_j]
         if _idx >= _source_start:
             _start_idx = _j
             break
-    if _start_idx == -1:
-        _start_idx = _num_spikes
-    for _j in range(_start_idx, _num_spikes):
+    for _j in range(_num_spikes-1, _start_idx-1, -1):
         _idx = {{_spikespace}}[_j]
-        if _idx >= _source_stop:
-            _end_idx = _j
+        if _idx < _source_stop:
             break
-    if _end_idx == -1:
-        _end_idx =_num_spikes
+        _end_idx = _j
     _num_spikes = _end_idx - _start_idx
+    {% else %}
+    for _j in range(_num_spikes):
+        _idx = {{_spikespace}}[_j]
+        if _idx < {{_source_indices}}[_source_index_counter]:
+            continue
+        while {{_source_indices}}[_source_index_counter] < _idx:
+            _source_index_counter += 1
+        if (_source_index_counter < {{source_N}} and
+                _idx == {{_source_indices}}[_source_index_counter]):
+            _source_index_counter += 1
+            _filtered_events[_filtered_events[{{source_N}}]] = _idx
+            _filtered_events[{{source_N}}] += 1
+        if _source_index_counter == {{source_N}}:
+            break
+    _num_spikes = _filtered_events[{{source_N}}]
+    {% endif %}
+    {% endif %}
 
     # Calculate the new length for the arrays
     cdef size_t _new_len = {{_dynamic_t}}.shape[0] + 1
@@ -36,6 +55,6 @@
 
     # Set the new values
     {{_dynamic_t}}.data[_new_len-1] = {{_clock_t}}
-    {{_dynamic_rate}}.data[_new_len-1] = _num_spikes/{{_clock_dt}}/_num_source_neurons
+    {{_dynamic_rate}}.data[_new_len-1] = _num_spikes/{{_clock_dt}}/{{source_N}}
 
 {% endblock %}

brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx

@@ -1,5 +1,4 @@
-{# USES_VARIABLES { N, _clock_t, count,
-                    _source_start, _source_stop} #}
+{# USES_VARIABLES { N, _clock_t, count} #}
 {% extends 'common_group.pyx' %}
 
 {% block maincode %}
@@ -9,11 +8,20 @@
 
     cdef size_t _num_events = {{_eventspace}}[_num{{_eventspace}}-1]
     cdef size_t _start_idx, _end_idx, _curlen, _newlen, _j
+    {% if subgroup and not contiguous %}
+    # We use the same data structure as for the eventspace to store the
+    # "filtered" events, i.e. the events that are indexed in the subgroup
+    cdef int[{{source_N}} + 1] _filtered_events
+    cdef size_t _source_index_counter = 0
+    _filtered_events[{{source_N}}] = 0
+    {% endif %}
     {% for varname, var in record_variables | dictsort %}
     cdef {{cpp_dtype(var.dtype)}}[:] _{{varname}}_view
     {% endfor %}
     if _num_events > 0:
+        {% if subgroup %}
         # For subgroups, we do not want to record all spikes
+        {% if contiguous %}
         # We assume that spikes are ordered
         _start_idx = _num_events
         _end_idx = _num_events
@@ -28,6 +36,23 @@
                 break
             _end_idx = _j
         _num_events = _end_idx - _start_idx
+        {% else %}
+        for _j in range(_num_events):
+            _idx = {{_eventspace}}[_j]
+            if _idx < {{_source_indices}}[_source_index_counter]:
+                continue
+            while {{_source_indices}}[_source_index_counter] < _idx:
+                _source_index_counter += 1
+            if (_source_index_counter < {{source_N}} and
+                    _idx == {{_source_indices}}[_source_index_counter]):
+                _source_index_counter += 1
+                _filtered_events[_filtered_events[{{source_N}}]] = _idx
+                _filtered_events[{{source_N}}] += 1
+            if _source_index_counter == {{source_N}}:
+                break
+        _num_events = _filtered_events[{{source_N}}]
+        {% endif %}
+        {% endif %}
         if _num_events > 0:
             # scalar code
             _vectorisation_idx = 1
@@ -41,6 +66,8 @@
             _{{varname}}_view = {{get_array_name(var, access_data=False)}}.data
             {% endfor %}
             # Copy the values across
+            {% if subgroup %}
+            {% if contiguous %}
             for _j in range(_start_idx, _end_idx):
                 _idx = {{_eventspace}}[_j]
                 _vectorisation_idx = _idx
@@ -49,4 +76,24 @@
                 _{{varname}}_view [_curlen + _j - _start_idx] = _to_record_{{varname}}
                 {% endfor %}
                 {{count}}[_idx - _source_start] += 1
+            {% else %}
+            for _j in range(_num_events):
+                _idx = _filtered_events[_j]
+                _vectorisation_idx = _idx
+                {{ vector_code|autoindent }}
+                {% for varname in record_variables | sort %}
+                _{{varname}}_view [_curlen + _j] = _to_record_{{varname}}
+                {% endfor %}
+                {{count}}[_to_record_i] += 1
+            {% endif %}
+            {% else %}
+            for _j in range(_num_events):
+                _idx = {{_eventspace}}[_j]
+                _vectorisation_idx = _idx
+                {{ vector_code|autoindent }}
+                {% for varname in record_variables | sort %}
+                _{{varname}}_view [_curlen + _j] = _to_record_{{varname}}
+                {% endfor %}
+                {{count}}[_idx] += 1
+            {% endif %}
 {% endblock %}

brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_

@@ -1,14 +1,19 @@
-{# USES_VARIABLES { rate, t, _spikespace, _num_source_neurons,
-                    _clock_t, _clock_dt, _source_start, _source_stop, N } #}
+{# USES_VARIABLES { rate, t, _spikespace, _clock_t, _clock_dt, N } #}
 {% extends 'common_group.py_' %}
 
 {% block maincode %}
 _spikes = {{_spikespace}}[:{{_spikespace}}[-1]]
+{% if subgroup %}
 # Take subgroups into account
+{% if contiguous %}
 _spikes = _spikes[(_spikes >= _source_start) & (_spikes < _source_stop)]
+{% else %}
+_spikes = _numpy.intersect1d(_spikes, {{_source_indices}}, assume_unique=True)
+{% endif %}
+{% endif %}
 _new_len = {{N}} + 1
 _owner.resize(_new_len)
 {{N}} = _new_len
 {{_dynamic_t}}[-1] = {{_clock_t}}
-{{_dynamic_rate}}[-1] = 1.0 * len(_spikes) / {{_clock_dt}} / _num_source_neurons
+{{_dynamic_rate}}[-1] = 1.0 * len(_spikes) / {{_clock_dt}} / {{source_N}}
 {% endblock %}

brian2/codegen/runtime/numpy_rt/templates/spikemonitor.py_

@@ -1,4 +1,4 @@
-{# USES_VARIABLES {N, count, _clock_t, _source_start, _source_stop, _source_N} #}
+{# USES_VARIABLES {N, count, _clock_t} #}
 {% extends 'common_group.py_' %}
 
 {% block maincode %}
@@ -9,11 +9,15 @@
 _n_events = {{_eventspace}}[-1]
 if _n_events > 0:
     _events = {{_eventspace}}[:_n_events]
+    {% if subgroup %}
     # Take subgroups into account
-    if _source_start != 0 or _source_stop != _source_N:
-        _events = _events[(_events >= _source_start) & (_events < _source_stop)]
-        _n_events = len(_events)
-
+    {% if contiguous %}
+    _events = _events[(_events >= _source_start) & (_events < _source_stop)]
+    {% else %}
+    _events = _numpy.intersect1d(_events, {{_source_indices}}, assume_unique=True)
+    {% endif %}
+    _n_events = len(_events)
+    {% endif %}
     if _n_events > 0:
         _vectorisation_idx = 1
         {{scalar_code|autoindent}}
@@ -28,5 +32,13 @@ if _n_events > 0:
         {% set dynamic_varname = get_array_name(var, access_data=False) %}
         {{dynamic_varname}}[_curlen:_newlen] = _to_record_{{varname}}
         {% endfor %}
+        {% if not subgroup %}
+        {{count}}[_events] += 1
+        {% else %}
+        {% if contiguous %}
         {{count}}[_events - _source_start] += 1
+        {% else %}
+        {{count}}[_to_record_i] += 1
+        {% endif %}
+        {% endif %}
 {% endblock %}

brian2/devices/cpp_standalone/device.py

@@ -17,8 +17,9 @@
 import numpy as np
 
 import brian2
-from brian2.codegen.codeobject import check_compiler_kwds
+from brian2.codegen.codeobject import check_compiler_kwds, create_runner_codeobj
 from brian2.codegen.cpp_prefs import get_compiler_and_args, get_msvc_env
+from brian2.codegen.runtime.numpy_rt import NumpyCodeObject
 from brian2.core.network import Network
 from brian2.devices.device import Device, all_devices, set_device, reset_device
 from brian2.core.functions import Function
@@ -475,7 +476,10 @@ def variableview_set_with_index_array(self, variableview, item,
                                                       indices)
             staticarrayname_value = self.static_array(f"_value_{arrayname}",
                                                       value)
-            self.array_cache[variableview.variable] = None
+            # Put values into the cache
+            cache_variable = self.array_cache[variableview.variable]
+            if cache_variable is not None:
+                cache_variable[indices] = value
             self.main_queue.append(('set_array_by_array', (arrayname,
                                                            staticarrayname_index,
                                                            staticarrayname_value)))
@@ -541,6 +545,36 @@ def variableview_get_with_expression(self, variableview, code,
                                   "variables with string expressions in "
                                   "standalone scripts.")
 
+    def index_wrapper_get_item(self, index_wrapper, item, level):
+        if isinstance(item, str):
+            variables = Variables(None)
+            variables.add_auxiliary_variable('_indices', dtype=np.int32)
+            variables.add_auxiliary_variable('_cond', dtype=np.bool)
+
+            abstract_code = '_cond = ' + item
+            namespace = get_local_namespace(level=level+2)
+            try:
+                codeobj = create_runner_codeobj(index_wrapper.group,
+                                                abstract_code,
+                                                'group_get_indices',
+                                                run_namespace=namespace,
+                                                additional_variables=variables,
+                                                codeobj_class=NumpyCodeObject
+                                                )
+            except NotImplementedError:
+                raise NotImplementedError('Cannot calculate indices with string '
+                                          'expressions in standalone mode if '
+                                          'the expression refers to variable '
+                                          'with values not known before running '
+                                          'the simulation.')
+            indices = codeobj()
+            # Delete code object from device to avoid trying to build it later
+            del self.code_objects[codeobj.name]
+            # Handle subgroups correctly
+            return index_wrapper.indices(indices)
+        else:
+            return index_wrapper.indices(item)
+
     def code_object_class(self, codeobj_class=None, fallback_pref=None):
         """
         Return `CodeObject` class (either `CPPStandaloneCodeObject` class or input)

brian2/devices/cpp_standalone/templates/ratemonitor.cpp

@@ -1,36 +1,66 @@
-{# USES_VARIABLES { N, rate, t, _spikespace, _clock_t, _clock_dt,
-                    _num_source_neurons, _source_start, _source_stop } #}
+{# USES_VARIABLES { N, rate, t, _spikespace, _clock_t, _clock_dt} #}
 {# WRITES_TO_READ_ONLY_VARIABLES { N } #}
 {% extends 'common_group.cpp' %}
 
 {% block maincode %}
     size_t _num_spikes = {{_spikespace}}[_num_spikespace-1];
+    {% if subgroup and not contiguous %}
+    // We use the same data structure as for the eventspace to store the
+    // "filtered" events, i.e. the events that are indexed in the subgroup
+    int32_t _filtered_events[{{source_N}} + 1];
+    _filtered_events[{{source_N}}] = 0;
+    size_t _source_index_counter = 0;
+    {% endif %}
+    {% if subgroup %}
     // For subgroups, we do not want to record all spikes
-    // We assume that spikes are ordered
-    int _start_idx = -1;
-    int _end_idx = -1;
-    for(size_t _j=0; _j<_num_spikes; _j++)
+    size_t _start_idx = _num_spikes;
+    size_t _end_idx = _num_spikes;
+    if (_num_spikes > 0)
     {
-        const size_t _idx = {{_spikespace}}[_j];
-        if (_idx >= _source_start) {
-            _start_idx = _j;
-            break;
+        {% if contiguous %}
+        for(size_t _j=0; _j<_num_spikes; _j++)
+        {
+            const int _idx = {{_spikespace}}[_j];
+            if (_idx >= _source_start) {
+                _start_idx = _j;
+                break;
+            }
         }
-    }
-    if (_start_idx == -1)
-        _start_idx = _num_spikes;
-    for(size_t _j=_start_idx; _j<_num_spikes; _j++)
-    {
-        const size_t _idx = {{_spikespace}}[_j];
-        if (_idx >= _source_stop) {
+        for(size_t _j=_num_spikes-1; _j>=_start_idx; _j--)
+        {
+            const int _idx = {{_spikespace}}[_j];
+            if (_idx < _source_stop) {
+                break;
+            }
             _end_idx = _j;
-            break;
         }
+        _num_spikes = _end_idx - _start_idx;
+        {% else %}
+        for (size_t _j=0; _j<_num_spikes; _j++)
+        {
+            const size_t _idx = {{_spikespace}}[_j];
+            if (_idx < {{_source_indices}}[_source_index_counter])
+                continue;
+            while ({{_source_indices}}[_source_index_counter] < _idx)
+            {
+                _source_index_counter++;
+            }
+            if (_source_index_counter < {{source_N}} &&
+                _idx == {{_source_indices}}[_source_index_counter])
+            {
+                _source_index_counter += 1;
+                _filtered_events[_filtered_events[{{source_N}}]++] = _idx;
+                if (_source_index_counter == {{source_N}})
+                    break;
+            }
+            if (_source_index_counter == {{source_N}})
+                break;
+        }
+        _num_spikes = _filtered_events[{{source_N}}];
+        {% endif %}
     }
-    if (_end_idx == -1)
-        _end_idx =_num_spikes;
-    _num_spikes = _end_idx - _start_idx;
-    {{_dynamic_rate}}.push_back(1.0*_num_spikes/{{_clock_dt}}/_num_source_neurons);
+    {% endif %}
+    {{_dynamic_rate}}.push_back(1.0*_num_spikes/{{_clock_dt}}/{{source_N}});
     {{_dynamic_t}}.push_back({{_clock_t}});
     {{N}}++;
 {% endblock %}