From 604e63316a762847933c85ce19dcce95a7541ab4 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 21 Feb 2020 18:14:05 +0100
Subject: [PATCH 01/52] WIP: First steps to generalize subgroups to be
 non-contiguous

---
 brian2/groups/neurongroup.py          |  83 +----------
 brian2/groups/subgroup.py             | 200 +++++++++++++++++++++++---
 brian2/spatialneuron/spatialneuron.py |  11 +-
 brian2/tests/test_subgroup.py         |  80 +++++++----
 4 files changed, 234 insertions(+), 140 deletions(-)

diff --git a/brian2/groups/neurongroup.py b/brian2/groups/neurongroup.py
index 5058e1bcc..988c00661 100644
--- a/brian2/groups/neurongroup.py
+++ b/brian2/groups/neurongroup.py
@@ -33,7 +33,7 @@
 from brian2.utils.stringtools import get_identifiers
 
 from .group import Group, CodeRunner, get_dtype
-from .subgroup import Subgroup
+from .subgroup import Subgroup, to_start_stop_or_index
 
 __all__ = ['NeuronGroup']
 
@@ -93,82 +93,6 @@ def check_identifier_pre_post(identifier):
                          'post-synaptic variables in synapses.' % identifier)
 
 
-def to_start_stop(item, N):
-    '''
-    Helper function to transform a single number, a slice or an array of
-    contiguous indices to a start and stop value. This is used to allow for
-    some flexibility in the syntax of specifying subgroups in `.NeuronGroup`
-    and `.SpatialNeuron`.
-
-    Parameters
-    ----------
-    item : slice, int or sequence
-        The slice, index, or sequence of indices to use. Note that a sequence
-        of indices has to be a sorted ascending sequence of subsequent integers.
-    N : int
-        The total number of elements in the group.
-
-    Returns
-    -------
-    start : int
-        The start value of the slice.
-    stop : int
-        The stop value of the slice.
-
-    Examples
-    --------
-    >>> from brian2.groups.neurongroup import to_start_stop
-    >>> to_start_stop(slice(3, 6), 10)
-    (3, 6)
-    >>> to_start_stop(slice(3, None), 10)
-    (3, 10)
-    >>> to_start_stop(5, 10)
-    (5, 6)
-    >>> to_start_stop([3, 4, 5], 10)
-    (3, 6)
-    >>> to_start_stop([3, 5, 7], 10)
-    Traceback (most recent call last):
-        ...
-    IndexError: Subgroups can only be constructed using contiguous indices.
-    
-    '''
-    if isinstance(item, slice):
-        start, stop, step = item.indices(N)
-    elif isinstance(item, numbers.Integral):
-        start = item
-        stop = item + 1
-        step = 1
-    elif (isinstance(item, (Sequence, np.ndarray)) and
-          not isinstance(item, str)):
-        if not (len(item) > 0 and np.all(np.diff(item) == 1)):
-            raise IndexError('Subgroups can only be constructed using '
-                             'contiguous indices.')
-        if not np.issubdtype(np.asarray(item).dtype, np.integer):
-            raise TypeError('Subgroups can only be constructed using integer '
-                            'values.')
-        start = int(item[0])
-        stop = int(item[-1]) + 1
-        step = 1
-    else:
-        raise TypeError('Subgroups can only be constructed using slicing '
-                        'syntax, a single index, or an array of contiguous '
-                        'indices.')
-    if step != 1:
-        raise IndexError('Subgroups have to be contiguous')
-    if start >= stop:
-        raise IndexError('Illegal start/end values for subgroup, %d>=%d' %
-                         (start, stop))
-    if start >= N:
-        raise IndexError('Illegal start value for subgroup, %d>=%d' %
-                         (start, N))
-    if stop > N:
-        raise IndexError('Illegal stop value for subgroup, %d>%d' %
-                         (stop, N))
-    if start < 0:
-        raise IndexError('Indices have to be positive.')
-    return start, stop
-
-
 class StateUpdater(CodeRunner):
     '''
     The `CodeRunner` that updates the state variables of a `NeuronGroup`
@@ -805,9 +729,8 @@ def __setattr__(self, key, value):
                 Group.__setattr__(self, key, value, level=1)
 
     def __getitem__(self, item):
-        start, stop = to_start_stop(item, self._N)
-
-        return Subgroup(self, start, stop)
+        start, stop, indices = to_start_stop_or_index(item, self._N)
+        return Subgroup(self, start, stop, indices)
 
     def _create_variables(self, user_dtype, events):
         '''
diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index 3a9ca6c08..89a1e8909 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -1,3 +1,7 @@
+import numbers
+from collections.abc import Sequence
+
+import numpy as np
 
 from brian2.core.base import weakproxy_with_fallback
 from brian2.core.spikesource import SpikeSource
@@ -9,6 +13,96 @@
 __all__ = ['Subgroup']
 
 
+def to_start_stop_or_index(item, N):
+    '''
+    Helper function to transform a single number, a slice or an array of
+    indices to a start and stop value (if possible), or to an index of positive
+    indices (interpreting negative indices correctly). This is used to allow for
+    some flexibility in the syntax of specifying subgroups in `.NeuronGroup`
+    and `.SpatialNeuron`.
+
+    Parameters
+    ----------
+    item : slice, int or sequence
+        The slice, index, or sequence of indices to use.
+    N : int
+        The total number of elements in the group.
+
+    Returns
+    -------
+    start : int or None
+        The start value of the slice.
+    stop : int or None
+        The stop value of the slice.
+    indices : `np.ndarray` or None
+        The indices.
+
+    Examples
+    --------
+    >>> from brian2.groups.neurongroup import to_start_stop_or_index
+    >>> to_start_stop_or_index(slice(3, 6), 10)
+    (3, 6, None)
+    >>> to_start_stop_or_index(slice(3, None), 10)
+    (3, 10, None)
+    >>> to_start_stop_or_index(5, 10)
+    (5, 6, None)
+    >>> to_start_stop_or_index(slice(None, None, 2), 10)
+    (None, None, array([0, 2, 4, 6, 8]))
+    >>> to_start_stop_or_index([3, 4, 5], 10)
+    (3, 6, None)
+    >>> to_start_stop_or_index([3, 5, 7], 10)
+    (None, None, array([3, 5, 7]))
+    >>> to_start_stop_or_index([-1, -2, -3], 10)
+    (None, None, array([9, 8, 7]))
+    '''
+    start = stop = indices = None
+    if isinstance(item, slice):
+        start, stop, step = item.indices(N)
+        if step != 1:
+            indices = np.arange(start, stop, step)
+            start = stop = None
+    elif isinstance(item, numbers.Integral):
+        start = item
+        stop = item + 1
+        step = 1
+    elif (isinstance(item, (Sequence, np.ndarray)) and
+          not isinstance(item, str)):
+        item = np.asarray(item)
+        if not np.issubdtype(item.dtype, np.integer):
+            raise TypeError('Subgroups can only be constructed using integer '
+                            'values.')
+        if not len(item) > 0:
+            raise IndexError('Cannot create an empty subgroup')
+        sorted_indices = sorted(item)
+        if np.all(np.diff(sorted_indices) == 1) and np.min(sorted_indices) >= 0:
+            start = int(item[0])
+            stop = int(item[-1]) + 1
+        else:
+            if np.min(item) < -N:
+                raise IndexError('Illegal index {} for a group of size '
+                                 '{}'.format(np.min(item), N))
+            elif np.max(item) >= N:
+                raise IndexError('Illegal index {} for a group of size '
+                                 '{}'.format(np.min(item), N))
+            indices = item % N  # interpret negative indices correctly
+    else:
+        raise TypeError("Cannot interpret object of type '{}' as index or "
+                        "slice".format(type(item)))
+
+    if indices is None:
+        if start >= stop:
+            raise IndexError('Illegal start/end values for subgroup, %d>=%d' %
+                             (start, stop))
+        if start >= N:
+            raise IndexError('Illegal start value for subgroup, %d>=%d' %
+                             (start, N))
+        if stop > N:
+            raise IndexError('Illegal stop value for subgroup, %d>%d' %
+                             (stop, N))
+
+    return start, stop, indices
+
+
 class Subgroup(Group, SpikeSource):
     '''
     Subgroup of any `Group`
@@ -17,18 +111,54 @@ class Subgroup(Group, SpikeSource):
     ----------
     source : SpikeSource
         The source object to subgroup.
-    start, stop : int
-        Select only spikes with indices from ``start`` to ``stop-1``.
+    start, stop : int, optional
+        Select only spikes with indices from ``start`` to ``stop-1``. Cannot
+        be specified at the same time as ``indices``.
+    indices : `np.ndarray`, optional
+        The indices of the subgroup. Note that subgroups with non-contiguous
+        indices cannot be used everywhere. Cannot be specified at the same time
+        as ``start`` and ``stop``.
     name : str, optional
         A unique name for the group, or use ``source.name+'_subgroup_0'``, etc.
     '''
-    def __init__(self, source, start, stop, name=None):
+    def __init__(self, source, start=None, stop=None, indices=None, name=None):
+        if start is stop is indices is None:
+            raise TypeError('Need to specify either start and stop or indices.')
+        if start != stop and (start is None or stop is None):
+            raise TypeError('start and stop have to be specified together.')
+        if indices is not None and (start is not None):
+            raise TypeError('Cannot specify both sub_indices and start and '
+                            'stop.')
+        if start is not None:
+            self.contiguous = True
+        else:
+            self.contiguous = False
+            if not len(indices):
+                raise IndexError('Cannot create an empty subgroup.')
+            max_index = np.max(indices)
+            if max_index >= len(source):
+                raise IndexError('Index {} cannot be >= the size of the group '
+                                 '({})'.format(max_index, len(source)))
+            if len(indices) != len(np.unique(indices)):
+                raise IndexError('sub_indices cannot contain repeated values.')
         # First check if the source is itself a Subgroup
         # If so, then make this a Subgroup of the original Group
         if isinstance(source, Subgroup):
-            source = source.source
-            start = start + source.start
-            stop = stop + source.start
+            if source.contiguous:
+                if self.contiguous:
+                    source = source.source
+                    start = start + source.start
+                    stop = stop + source.start
+                else:
+                    if np.max(indices) >= source.stop - source.start:
+                        raise IndexError('Index exceeds range')
+                    indices += source.start
+            else:
+                if self.contiguous:
+                    indices = source.sub_indices[start:stop]
+                    self.contiguous = False
+                else:
+                    indices = source.sub_indices[indices]
             self.source = source
         else:
             self.source = weakproxy_with_fallback(source)
@@ -48,9 +178,13 @@ def __init__(self, source, start, stop, name=None):
                        clock=source._clock,
                        when='thresholds',
                        order=source.order+1, name=name)
-        self._N = stop-start
+        if self.contiguous:
+            self._N = stop-start
+        else:
+            self._N = len(indices)
         self.start = start
         self.stop = stop
+        self.sub_indices = indices
 
         self.events = self.source.events
 
@@ -59,16 +193,25 @@ def __init__(self, source, start, stop, name=None):
         self.variables = Variables(self, default_index='_sub_idx')
 
         # overwrite the meaning of N and i
-        if self.start > 0:
+        if self.contiguous and self.start > 0:
             self.variables.add_constant('_offset', value=self.start)
             self.variables.add_reference('_source_i', source, 'i')
             self.variables.add_subexpression('i',
                                              dtype=source.variables['i'].dtype,
                                              expr='_source_i - _offset',
                                              index='_idx')
-        else:
+        elif self.contiguous:
             # no need to calculate anything if this is a subgroup starting at 0
             self.variables.add_reference('i', source)
+        else:
+            # We need an array to invert the indexing, i.e. an array where you
+            # can use the sub_indices and get back 0, 1, 2, ...
+            inv_idx = np.zeros(np.max(indices) + 1)
+            inv_idx[indices] = np.arange(len(indices))
+            self.variables.add_array('i', size=len(inv_idx),
+                                     dtype=source.variables['i'].dtype,
+                                     values=inv_idx, constant=True,
+                                     read_only=True, unique=True)
 
         self.variables.add_constant('N', value=self._N)
         self.variables.add_constant('_source_N', value=len(source))
@@ -77,13 +220,18 @@ def __init__(self, source, start, stop, name=None):
 
         # Only the variable _sub_idx itself is stored in the subgroup
         # and needs the normal index for this group
-        self.variables.add_arange('_sub_idx', size=self._N, start=self.start,
-                                  index='_idx')
+        if self.contiguous:
+            self.variables.add_arange('_sub_idx', size=self._N, start=self.start,
+                                      index='_idx')
+        else:
+            self.variables.add_array('_sub_idx', size=self._N,
+                                     dtype=np.int32, values=indices,
+                                     index='_idx')
 
         # special indexing for subgroups
         self._indices = Indexing(self, self.variables['_sub_idx'])
 
-        # Deal with special indices
+        # Deal with special sub_indices
         for key, value in self.source.variables.indices.items():
             if value == '0':
                 self.variables.indices[key] = '0'
@@ -91,7 +239,7 @@ def __init__(self, source, start, stop, name=None):
                 continue  # nothing to do, already uses _sub_idx correctly
             else:
                 raise ValueError(('Do not know how to deal with variable %s '
-                                  'using  index %s in a subgroup') % (key,
+                                  'using index %s in a subgroup') % (key,
                                                                       value))
 
         self.namespace = self.source.namespace
@@ -102,20 +250,26 @@ def __init__(self, source, start, stop, name=None):
     spikes = property(lambda self: self.source.spikes)
 
     def __getitem__(self, item):
-        if not isinstance(item, slice):
-            raise TypeError('Subgroups can only be constructed using slicing syntax')
-        start, stop, step = item.indices(self._N)
-        if step != 1:
-            raise IndexError('Subgroups have to be contiguous')
-        if start >= stop:
-            raise IndexError('Illegal start/end values for subgroup, %d>=%d' %
-                             (start, stop))
-        return Subgroup(self.source, self.start + start, self.start + stop)
+        start, stop, indices = to_start_stop_or_index(item, self._N)
+        if self.contiguous:
+            if start is not None:
+                return Subgroup(self.source, self.start + start, self.start + stop)
+            else:
+                return Subgroup(self.source, indices=indices+self.start)
+        else:
+            indices = self.sub_indices[item]
+            return Subgroup(self.source, indices=indices)
 
     def __repr__(self):
-        description = '<{classname} {name} of {source} from {start} to {end}>'
+        if self.contiguous:
+            description = '<{classname} {name} of {source} from {start} to {end}>'
+            str_indices = None
+        else:
+            description = '<{classname} {name} of {source} with indices {indices}>'
+            str_indices = np.array2string(self.sub_indices, threshold=10, separator=', ')
         return description.format(classname=self.__class__.__name__,
                                   name=repr(self.name),
                                   source=repr(self.source.name),
                                   start=self.start,
+                                  indices=str_indices,
                                   end=self.stop)
diff --git a/brian2/spatialneuron/spatialneuron.py b/brian2/spatialneuron/spatialneuron.py
index 563abe4b8..ef702bce2 100644
--- a/brian2/spatialneuron/spatialneuron.py
+++ b/brian2/spatialneuron/spatialneuron.py
@@ -23,7 +23,7 @@
 from brian2.parsing.sympytools import sympy_to_str, str_to_sympy
 from brian2.utils.logger import get_logger
 from brian2.groups.neurongroup import (NeuronGroup, SubexpressionUpdater,
-                                       to_start_stop)
+                                       to_start_stop_or_index)
 from brian2.groups.subgroup import Subgroup
 from brian2.equations.codestrings import Expression
 
@@ -492,13 +492,10 @@ def spatialneuron_segment(neuron, item):
             indices = neuron.morphology.indices[item]
             start, stop = indices[0], indices[-1] + 1
         elif not isinstance(item, slice) and hasattr(item, 'indices'):
-            start, stop = to_start_stop(item.indices[:], neuron._N)
+            start, stop, indices = to_start_stop_or_index(item.indices[:],
+                                                          neuron._N)
         else:
-            start, stop = to_start_stop(item, neuron._N)
-
-        if start >= stop:
-            raise IndexError('Illegal start/end values for subgroup, %d>=%d' %
-                             (start, stop))
+            start, stop, indices = to_start_stop_or_index(item, neuron._N)
 
         return Subgroup(neuron, start, stop)
 
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 627c405c1..302a04e00 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -17,9 +17,12 @@ def test_str_repr():
     '''
     G = NeuronGroup(10, 'v:1')
     SG = G[5:8]
+    SGi = G[[3, 6, 9]]
     # very basic test, only make sure no error is raised
     assert len(str(SG))
     assert len(repr(SG))
+    assert len(str(SGi))
+    assert len(repr(SGi))
 
 
 def test_state_variables():
@@ -85,6 +88,27 @@ def test_state_variables_simple():
     assert_equal(G.c[:], [0, 0, 0, 3, 4, 5, 6, 0, 0, 0])
     assert_equal(G.d[:], [0, 0, 0, 0, 4, 1, 2, 0, 0, 0])
 
+@pytest.mark.standalone_compatible
+def test_state_variables_simple_indexed():
+    G = NeuronGroup(10, '''a : 1
+                           b : 1
+                           c : 1
+                           d : 1
+                           ''')
+    SG = G[[3, 5, 7, 9]]
+    SG.a = 1
+    SG.a['i == 0'] = 2
+    SG.b = 'i'
+    SG.b['i == 3'] = 'i * 2'
+    SG.c = np.arange(3, 7)
+    SG.d[1:2] = 4
+    SG.d[2:4] = [1, 2]
+    run(0*ms)
+    assert_equal(G.a[:], [0, 0, 0, 2, 0, 1, 0, 1, 0, 1])
+    assert_equal(G.b[:], [0, 0, 0, 0, 0, 1, 0, 2, 0, 6])
+    assert_equal(G.c[:], [0, 0, 0, 3, 0, 4, 0, 5, 0, 6])
+    assert_equal(G.d[:], [0, 0, 0, 0, 0, 4, 0, 1, 0, 2])
+
 
 def test_state_variables_string_indices():
     '''
@@ -641,42 +665,27 @@ def test_spike_monitor():
 
 
 @pytest.mark.codegen_independent
-def test_wrong_indexing():
+@pytest.mark.parametrize('item', ['string', slice(10, None), slice(3, 2),
+                                  [9, 10], [10, 11], [2.5, 3.5, 4.5],
+                                  [5, 5, 5], []])
+def test_wrong_indexing(item):
     G = NeuronGroup(10, 'v:1')
-    with pytest.raises(TypeError):
-        G['string']
-
-    with pytest.raises(IndexError):
-        G[10]
-    with pytest.raises(IndexError):
-        G[10:]
-    with pytest.raises(IndexError):
-        G[::2]
-    with pytest.raises(IndexError):
-        G[3:2]
-    with pytest.raises(IndexError):
-        G[[5, 4, 3]]
-    with pytest.raises(IndexError):
-        G[[2, 4, 6]]
-    with pytest.raises(IndexError):
-        G[[-1, 0, 1]]
-    with pytest.raises(IndexError):
-        G[[9, 10, 11]]
-    with pytest.raises(IndexError):
-        G[[9, 10]]
-    with pytest.raises(IndexError):
-        G[[10, 11]]
-    with pytest.raises(TypeError):
-        G[[2.5, 3.5, 4.5]]
+    with pytest.raises((TypeError, IndexError)):
+        G[item]
 
 
 @pytest.mark.codegen_independent
-def test_alternative_indexing():
+@pytest.mark.parametrize('item,expected', [(slice(-3, None), np.array([7, 8, 9])),
+                                           (slice(None, None, 2), np.array([0, 2, 4, 6, 8])),
+                                           (3, np.array([3])),
+                                           ([3, 4, 5], np.array([3, 4, 5])),
+                                           ([3, 5, 7], np.array([3, 5, 7])),
+                                           ([7, 5, 3], np.array([7, 5, 3])),
+                                           ([3, -1], np.array([3, 9]))])
+def test_alternative_indexing(item, expected):
     G = NeuronGroup(10, 'v : integer')
     G.v = 'i'
-    assert_equal(G[-3:].v, np.array([7, 8, 9]))
-    assert_equal(G[3].v, np.array([3]))
-    assert_equal(G[[3, 4, 5]].v, np.array([3, 4, 5]))
+    assert_equal(G[item].v, expected)
 
 
 def test_no_reference_1():
@@ -742,10 +751,21 @@ def test_recursive_subgroup():
     G = NeuronGroup(10, 'v : 1')
     G.v = 'i'
     SG = G[3:8]
+    SGi = G[[3, 5, 7, 9]]
     SG2 = SG[2:4]
+    SGi2 = SGi[2:4]
+    SGii = SGi[[1, 3]]
+    SG2i = SG[[1, 3]]
     assert_equal(SG2.v[:], np.array([5, 6]))
     assert_equal(SG2.v[:], SG.v[2:4])
+    assert_equal(SGi2.v[:], np.array([7, 9]))
+    assert_equal(SGii.v[:], np.array([5, 9]))
+    assert_equal(SG2i.v[:], np.array([4, 6]))
     assert SG2.source.name == G.name
+    assert SGi2.source.name == G.name
+    assert SGii.source.name == G.name
+    assert SG2i.source.name == G.name
+
 
 if __name__ == '__main__':
     test_str_repr()

From d67ca90516bdb423fd19043d5468060b3bfaba8d Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 24 Feb 2020 18:32:25 +0100
Subject: [PATCH 02/52] Update SpikeMonitor for non-contiguous subgroups

---
 .../cython_rt/templates/spikemonitor.pyx      | 52 ++++++++++-
 .../numpy_rt/templates/spikemonitor.py_       | 22 ++++-
 brian2/devices/cpp_standalone/device.py       |  4 +-
 .../cpp_standalone/templates/spikemonitor.cpp | 93 +++++++++++++++----
 brian2/groups/subgroup.py                     |  3 +-
 brian2/monitors/spikemonitor.py               | 27 ++++--
 brian2/tests/test_monitor.py                  | 16 ++++
 brian2/tests/test_subgroup.py                 | 10 +-
 8 files changed, 191 insertions(+), 36 deletions(-)

diff --git a/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx b/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
index b5abb20e0..415787c9a 100644
--- a/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
+++ b/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,24 @@
                 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}}:
+                    break
+            if _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 +67,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 +77,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 %}
diff --git a/brian2/codegen/runtime/numpy_rt/templates/spikemonitor.py_ b/brian2/codegen/runtime/numpy_rt/templates/spikemonitor.py_
index 639328f41..60274dcf1 100644
--- a/brian2/codegen/runtime/numpy_rt/templates/spikemonitor.py_
+++ b/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 %}
diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index d3639ec8b..2dc4885e3 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/brian2/devices/cpp_standalone/device.py
@@ -468,7 +468,9 @@ def variableview_set_with_index_array(self, variableview, item,
                                                       indices)
             staticarrayname_value = self.static_array('_value_'+arrayname,
                                                       value)
-            self.array_cache[variableview.variable] = None
+            # Put values into the cache
+            cache_variable = self.array_cache[variableview.variable]
+            cache_variable[indices] = value
             self.main_queue.append(('set_array_by_array', (arrayname,
                                                            staticarrayname_index,
                                                            staticarrayname_value)))
diff --git a/brian2/devices/cpp_standalone/templates/spikemonitor.cpp b/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
index ca62a5fe8..8737b609a 100644
--- a/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
+++ b/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
@@ -1,6 +1,5 @@
-{# USES_VARIABLES { N, _clock_t, count,
-                        _source_start, _source_stop} #}
-    {# WRITES_TO_READ_ONLY_VARIABLES { N, count } #}
+{# USES_VARIABLES { N, _clock_t, count} #}
+{# WRITES_TO_READ_ONLY_VARIABLES { N, count } #}
 {% extends 'common_group.cpp' %}
 
 {% block maincode %}
@@ -9,11 +8,20 @@
     {% set _eventspace = get_array_name(eventspace_variable) %}
 
     int32_t _num_events = {{_eventspace}}[_num{{eventspace_variable.name}}-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
+    size_t _start_idx = _num_events;
+    size_t _end_idx = _num_events;
     if (_num_events > 0)
     {
-        size_t _start_idx = _num_events;
-        size_t _end_idx = _num_events;
+        {% if contiguous %}
         for(size_t _j=0; _j<_num_events; _j++)
         {
             const int _idx = {{_eventspace}}[_j];
@@ -31,21 +39,70 @@
             _end_idx = _j;
         }
         _num_events = _end_idx - _start_idx;
-        if (_num_events > 0) {
-            const size_t _vectorisation_idx = 1;
-            {{scalar_code|autoindent}}
-            for(size_t _j=_start_idx; _j<_end_idx; _j++)
+        {% else %}
+        for (size_t _j=0; _j<_num_events; _j++)
+        {
+            const size_t _idx = {{_eventspace}}[_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}})
+                    break;
+            }
+            if (_idx == {{_source_indices}}[_source_index_counter])
             {
-                const size_t _idx = {{_eventspace}}[_j];
-                const size_t _vectorisation_idx = _idx;
-                {{vector_code|autoindent}}
-                {% for varname, var in record_variables | dictsort %}
-                {{get_array_name(var, access_data=False)}}.push_back(_to_record_{{varname}});
-                {% endfor %}
-                {{count}}[_idx-_source_start]++;
+                _source_index_counter += 1;
+                _filtered_events[_filtered_events[{{source_N}}]++] = _idx;
+                if (_source_index_counter == {{source_N}})
+                    break;
             }
-            {{N}} += _num_events;
         }
+        _num_events = _filtered_events[{{source_N}}];
+        {% endif %}
+    }
+    {% endif %}
+    if (_num_events > 0) {
+        const size_t _vectorisation_idx = 1;
+        {{scalar_code|autoindent}}
+        {% if subgroup %}
+        {% if contiguous %}
+        for(size_t _j=_start_idx; _j<_end_idx; _j++)
+        {
+            const size_t _idx = {{_eventspace}}[_j];
+            const size_t _vectorisation_idx = _idx;
+            {{vector_code|autoindent}}
+            {% for varname, var in record_variables | dictsort %}
+            {{get_array_name(var, access_data=False)}}.push_back(_to_record_{{varname}});
+            {% endfor %}
+            {{count}}[_idx-_source_start]++;
+        }
+        {% else %}
+        for(size_t _j=0; _j < _num_events; _j++)
+        {
+            const size_t _idx = _filtered_events[_j];
+            const size_t _vectorisation_idx = _idx;
+            {{vector_code|autoindent}}
+            {% for varname, var in record_variables | dictsort %}
+            {{get_array_name(var, access_data=False)}}.push_back(_to_record_{{varname}});
+            {% endfor %}
+            {{count}}[_to_record_i]++;
+        }
+        {% endif %}
+        {% else %}
+        for (size_t _j=0; _j < _num_events; _j++)
+        {
+            const size_t _idx = {{_eventspace}}[_j];
+            const size_t _vectorisation_idx = _idx;
+            {{ vector_code|autoindent }}
+            {% for varname, var in record_variables | dictsort %}
+            {{get_array_name(var, access_data=False)}}.push_back(_to_record_{{varname}});
+            {% endfor %}
+            {{count}}[_idx]++;
+        }
+        {% endif %}
+        {{N}} += _num_events;
     }
 
 {% endblock %}
diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index 89a1e8909..c956ce64b 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -226,7 +226,8 @@ def __init__(self, source, start=None, stop=None, indices=None, name=None):
         else:
             self.variables.add_array('_sub_idx', size=self._N,
                                      dtype=np.int32, values=indices,
-                                     index='_idx')
+                                     index='_idx', constant=True,
+                                     read_only=True, unique=True)
 
         # special indexing for subgroups
         self._indices = Indexing(self, self.variables['_sub_idx'])
diff --git a/brian2/monitors/spikemonitor.py b/brian2/monitors/spikemonitor.py
index d495a8539..b9047c43f 100644
--- a/brian2/monitors/spikemonitor.py
+++ b/brian2/monitors/spikemonitor.py
@@ -8,6 +8,7 @@
 from brian2.core.spikesource import SpikeSource
 from brian2.units.fundamentalunits import Quantity, DIMENSIONLESS
 from brian2.groups.group import CodeRunner, Group
+from brian2.groups.subgroup import Subgroup
 
 __all__ = ['EventMonitor', 'SpikeMonitor']
 
@@ -135,21 +136,35 @@ def __init__(self, source, event, variables=None, record=True,
                                              dimensions=source_var.dim,
                                              dtype=source_var.dtype,
                                              read_only=True)
+        needed_variables = {eventspace_name} | self.record_variables
+        subgroup = isinstance(source, Subgroup)
+        contiguous = not subgroup or source.contiguous
         self.variables.add_arange('_source_idx', size=len(source))
-        self.variables.add_array('count', size=len(source), dtype=np.int32,
+        self.variables.add_array('count', size=len(source),
+                                 dtype=np.int32,
                                  read_only=True, index='_source_idx')
-        self.variables.add_constant('_source_start', start)
-        self.variables.add_constant('_source_stop', stop)
-        self.variables.add_constant('_source_N', source_N)
         self.variables.add_array('N', size=1, dtype=np.int32, read_only=True,
                                  scalar=True)
+        if subgroup:
+            if contiguous:
+                self.variables.add_constant('_source_start', start)
+                self.variables.add_constant('_source_stop', stop)
+                self.variables.add_constant('_source_N', source_N)
+                needed_variables |= {'_source_start', '_source_stop',
+                                     '_source_N'}
+            else:
+                self.variables.add_reference('_source_indices', source,
+                                             '_sub_idx', index='_source_idx')
+                needed_variables |= {'_source_indices'}
 
         record_variables = {varname: self.variables[varname]
                             for varname in self.record_variables}
         template_kwds = {'eventspace_variable': source.variables[eventspace_name],
                          'record_variables': record_variables,
-                         'record': self.record}
-        needed_variables = {eventspace_name} | self.record_variables
+                         'record': self.record,
+                         'subgroup': subgroup,
+                         'contiguous': contiguous,
+                         'source_N': source.N}
         CodeRunner.__init__(self, group=self, code=code, template='spikemonitor',
                             name=None,  # The name has already been initialized
                             clock=source.clock, when=when,
diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py
index 89a74dfb2..7420b7645 100644
--- a/brian2/tests/test_monitor.py
+++ b/brian2/tests/test_monitor.py
@@ -132,7 +132,12 @@ def test_spike_monitor_subgroups():
     spikes_1 = SpikeMonitor(G[:2])
     spikes_2 = SpikeMonitor(G[2:4])
     spikes_3 = SpikeMonitor(G[4:])
+    spikes_indexed = SpikeMonitor(G[::2])
+    with catch_logs() as l:
+        spikes_indexed_unsorted = SpikeMonitor(G[[4, 0, 2]])  # This should not make a difference
+    assert len(l) == 1 and l[0][1].endswith('.unsorted_subgroup_indices')
     run(defaultclock.dt)
+    # Spikes
     assert_allclose(spikes_all.i, [0, 4, 5])
     assert_allclose(spikes_all.t, [0, 0, 0]*ms)
     assert_allclose(spikes_1.i, [0])
@@ -141,6 +146,17 @@ def test_spike_monitor_subgroups():
     assert len(spikes_2.t) == 0
     assert_allclose(spikes_3.i, [0, 1])  # recorded spike indices are relative
     assert_allclose(spikes_3.t, [0, 0] * ms)
+    assert_allclose(spikes_indexed.i, [0, 2])
+    assert_allclose(spikes_indexed.t, [0, 0]*ms)
+    assert_allclose(spikes_indexed_unsorted.i, [0, 2])
+    assert_allclose(spikes_indexed_unsorted.t, [0, 0]*ms)
+    # Spike count
+    assert_allclose(spikes_all.count, [1, 0, 0, 0, 1, 1])
+    assert_allclose(spikes_1.count, [1, 0])
+    assert_allclose(spikes_2.count, [0, 0])
+    assert_allclose(spikes_3.count, [1, 1])
+    assert_allclose(spikes_indexed.count, [1, 0, 1])
+    assert_allclose(spikes_indexed_unsorted.count, [1, 0, 1])
 
 
 def test_spike_monitor_bug_824():
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 302a04e00..e6bf3a05a 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -639,15 +639,16 @@ def test_run_regularly():
 
 @pytest.mark.standalone_compatible
 def test_spike_monitor():
+    set_device('cpp_standalone', directory='/tmp/testsubgroup')
     G = NeuronGroup(10, 'v:1', threshold='v>1', reset='v=0')
-    G.v[0] = 1.1
-    G.v[2] = 1.1
-    G.v[5] = 1.1
+    G.v[[0, 2, 5]] = 1.1
     SG = G[3:]
     SG2 = G[:3]
+    SGi = G[[3, 5, 7]]
     s_mon = SpikeMonitor(G)
     sub_s_mon = SpikeMonitor(SG)
     sub_s_mon2 = SpikeMonitor(SG2)
+    sub_s_moni = SpikeMonitor(SGi)
     run(defaultclock.dt)
     assert_equal(s_mon.i, np.array([0, 2, 5]))
     assert_equal(s_mon.t_, np.zeros(3))
@@ -655,6 +656,8 @@ def test_spike_monitor():
     assert_equal(sub_s_mon.t_, np.zeros(1))
     assert_equal(sub_s_mon2.i, np.array([0, 2]))
     assert_equal(sub_s_mon2.t_, np.zeros(2))
+    assert_equal(sub_s_moni.i, np.array([1]))
+    assert_equal(sub_s_moni.t, np.zeros(1))
     expected = np.zeros(10, dtype=int)
     expected[[0, 2, 5]] = 1
     assert_equal(s_mon.count, expected)
@@ -662,6 +665,7 @@ def test_spike_monitor():
     expected[[2]] = 1
     assert_equal(sub_s_mon.count, expected)
     assert_equal(sub_s_mon2.count, np.array([1, 0, 1]))
+    assert_equal(sub_s_moni.count, np.array([0, 1, 0]))
 
 
 @pytest.mark.codegen_independent

From eeff84be118485ff9c9d20b8a3121a251a247a11 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Tue, 25 Feb 2020 10:58:26 +0100
Subject: [PATCH 03/52] Sort indices before their use in Subgroup

Warn to notify user if not sorted already
---
 brian2/groups/subgroup.py     | 33 ++++++++++++++++++++++-----------
 brian2/tests/test_subgroup.py |  2 +-
 2 files changed, 23 insertions(+), 12 deletions(-)

diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index c956ce64b..c323d57e3 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -7,11 +7,14 @@
 from brian2.core.spikesource import SpikeSource
 from brian2.core.variables import Variables
 from brian2.units.fundamentalunits import Unit
+from brian2.utils.logger import get_logger
 
 from .group import Group, Indexing
 
 __all__ = ['Subgroup']
 
+logger = get_logger(__name__)
+
 
 def to_start_stop_or_index(item, N):
     '''
@@ -53,7 +56,7 @@ def to_start_stop_or_index(item, N):
     >>> to_start_stop_or_index([3, 5, 7], 10)
     (None, None, array([3, 5, 7]))
     >>> to_start_stop_or_index([-1, -2, -3], 10)
-    (None, None, array([9, 8, 7]))
+    (None, None, array([7, 8, 9]))
     '''
     start = stop = indices = None
     if isinstance(item, slice):
@@ -73,18 +76,25 @@ def to_start_stop_or_index(item, N):
                             'values.')
         if not len(item) > 0:
             raise IndexError('Cannot create an empty subgroup')
-        sorted_indices = sorted(item)
-        if np.all(np.diff(sorted_indices) == 1) and np.min(sorted_indices) >= 0:
+        if np.min(item) < -N:
+            raise IndexError('Illegal index {} for a group of size '
+                             '{}'.format(np.min(item), N))
+        elif np.max(item) >= N:
+            raise IndexError('Illegal index {} for a group of size '
+                             '{}'.format(np.min(item), N))
+        indices = item % N  # interpret negative indices correctly
+
+        if not np.all(item[:-1] <= item[1:]):
+            logger.warn('The indices provided to create the subgroup were '
+                        'not sorted. They will be sorted before use.',
+                        name_suffix='unsorted_subgroup_indices')
+            indices.sort()
+
+        if np.all(np.diff(item) == 1):
             start = int(item[0])
             stop = int(item[-1]) + 1
-        else:
-            if np.min(item) < -N:
-                raise IndexError('Illegal index {} for a group of size '
-                                 '{}'.format(np.min(item), N))
-            elif np.max(item) >= N:
-                raise IndexError('Illegal index {} for a group of size '
-                                 '{}'.format(np.min(item), N))
-            indices = item % N  # interpret negative indices correctly
+            indices = None
+
     else:
         raise TypeError("Cannot interpret object of type '{}' as index or "
                         "slice".format(type(item)))
@@ -141,6 +151,7 @@ def __init__(self, source, start=None, stop=None, indices=None, name=None):
                                  '({})'.format(max_index, len(source)))
             if len(indices) != len(np.unique(indices)):
                 raise IndexError('sub_indices cannot contain repeated values.')
+
         # First check if the source is itself a Subgroup
         # If so, then make this a Subgroup of the original Group
         if isinstance(source, Subgroup):
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index e6bf3a05a..42cf6ff20 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -684,7 +684,7 @@ def test_wrong_indexing(item):
                                            (3, np.array([3])),
                                            ([3, 4, 5], np.array([3, 4, 5])),
                                            ([3, 5, 7], np.array([3, 5, 7])),
-                                           ([7, 5, 3], np.array([7, 5, 3])),
+                                           ([7, 5, 3], np.array([3, 5, 7])),
                                            ([3, -1], np.array([3, 9]))])
 def test_alternative_indexing(item, expected):
     G = NeuronGroup(10, 'v : integer')

From 97659da28737d374fa810bbb307a42bd55e3eec5 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Tue, 25 Feb 2020 16:18:43 +0100
Subject: [PATCH 04/52] Fix logic error in Cython/C++ SpikeMonitor code

---
 .../codegen/runtime/cython_rt/templates/spikemonitor.pyx | 9 ++++-----
 brian2/devices/cpp_standalone/templates/spikemonitor.cpp | 7 ++++---
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx b/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
index 415787c9a..78a74e59f 100644
--- a/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
+++ b/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
@@ -43,14 +43,13 @@
                 continue
             while {{_source_indices}}[_source_index_counter] < _idx:
                 _source_index_counter += 1
-                if _source_index_counter == {{source_N}}:
-                    break
-            if _idx == {{_source_indices}}[_source_index_counter]:
+            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
+            if _source_index_counter == {{source_N}}:
+                break
         _num_events = _filtered_events[{{source_N}}]
         {% endif %}
         {% endif %}
diff --git a/brian2/devices/cpp_standalone/templates/spikemonitor.cpp b/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
index 8737b609a..b94c61cc3 100644
--- a/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
+++ b/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
@@ -48,16 +48,17 @@
             while ({{_source_indices}}[_source_index_counter] < _idx)
             {
                 _source_index_counter++;
-                if (_source_index_counter == {{source_N}})
-                    break;
             }
-            if (_idx == {{_source_indices}}[_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_events = _filtered_events[{{source_N}}];
         {% endif %}

From 4672adf32e948e2497c5efe2b0a63cb86bd5a9db Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Tue, 25 Feb 2020 17:34:09 +0100
Subject: [PATCH 05/52] Make PopulationRateMonitor work with non-contiguous
 subgroups

---
 .../cython_rt/templates/ratemonitor.pyx       | 47 ++++++++----
 .../numpy_rt/templates/ratemonitor.py_        | 11 ++-
 .../cpp_standalone/templates/ratemonitor.cpp  | 74 +++++++++++++------
 brian2/monitors/ratemonitor.py                | 41 ++++++----
 brian2/tests/test_monitor.py                  |  3 +
 5 files changed, 124 insertions(+), 52 deletions(-)

diff --git a/brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx b/brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx
index a96582d04..80b852497 100644
--- a/brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx
+++ b/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 %}
diff --git a/brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_ b/brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_
index 46da9356f..aa4bb1149 100644
--- a/brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_
+++ b/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 %}
\ No newline at end of file
diff --git a/brian2/devices/cpp_standalone/templates/ratemonitor.cpp b/brian2/devices/cpp_standalone/templates/ratemonitor.cpp
index c06362543..f228a4e4e 100644
--- a/brian2/devices/cpp_standalone/templates/ratemonitor.cpp
+++ b/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 %}
diff --git a/brian2/monitors/ratemonitor.py b/brian2/monitors/ratemonitor.py
index 6dd0f8291..4c63beac1 100644
--- a/brian2/monitors/ratemonitor.py
+++ b/brian2/monitors/ratemonitor.py
@@ -7,7 +7,9 @@
 from brian2.core.variables import Variables
 from brian2.units.allunits import second, hertz
 from brian2.units.fundamentalunits import Unit, Quantity, check_units
+from brian2.core.names import Nameable
 from brian2.groups.group import CodeRunner, Group
+from brian2.groups.subgroup import Subgroup
 
 __all__ = ['PopulationRateMonitor']
 
@@ -42,33 +44,46 @@ class PopulationRateMonitor(Group, CodeRunner):
     add_to_magic_network = True
     def __init__(self, source, name='ratemonitor*', codeobj_class=None,
                  dtype=np.float64):
-
+        Nameable.__init__(self, name=name)
         #: The group we are recording from
         self.source = source
 
         self.codeobj_class = codeobj_class
-        CodeRunner.__init__(self, group=self, code='', template='ratemonitor',
-                            clock=source.clock, when='end', order=0, name=name)
-
-        self.add_dependency(source)
 
         self.variables = Variables(self)
+
         # Handle subgroups correctly
-        start = getattr(source, 'start', 0)
-        stop = getattr(source, 'stop', len(source))
-        self.variables.add_constant('_source_start', start)
-        self.variables.add_constant('_source_stop', stop)
-        self.variables.add_reference('_spikespace', source)
+        subgroup = isinstance(source, Subgroup)
+        contiguous = not subgroup or source.contiguous
+        self.variables.add_arange('_source_idx', size=len(source))
+        needed_variables = {}
+        if subgroup:
+            if contiguous:
+                self.variables.add_constant('_source_start', source.start)
+                self.variables.add_constant('_source_stop', source.stop)
+                needed_variables = {'_source_start', '_source_stop'}
+            else:
+                self.variables.add_reference('_source_indices', source,
+                                             '_sub_idx', index='_source_idx')
+                needed_variables = {'_source_indices'}
         self.variables.add_dynamic_array('rate', size=0, dimensions=hertz.dim,
                                          read_only=True, dtype=dtype)
+        CodeRunner.__init__(self, group=self, code='', template='ratemonitor',
+                            clock=source.clock, when='end', order=0, name=name,
+                            needed_variables=needed_variables,
+                            template_kwds={'subgroup': subgroup,
+                                           'contiguous': contiguous,
+                                           'source_N': source.N})
+        self.variables.create_clock_variables(self._clock,
+                                              prefix='_clock_')
         self.variables.add_dynamic_array('t', size=0, dimensions=second.dim,
                                          read_only=True,
                                          dtype=self._clock.variables['t'].dtype)
-        self.variables.add_reference('_num_source_neurons', source, 'N')
         self.variables.add_array('N', dtype=np.int32, size=1,
                                  scalar=True, read_only=True)
-        self.variables.create_clock_variables(self._clock,
-                                              prefix='_clock_')
+        self.variables.add_reference('_spikespace', source)
+
+        self.add_dependency(source)
         self._enable_group_attributes()
 
     def resize(self, new_size):
diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py
index 7420b7645..fd2f0485c 100644
--- a/brian2/tests/test_monitor.py
+++ b/brian2/tests/test_monitor.py
@@ -569,11 +569,14 @@ def test_rate_monitor_subgroups_2():
     rate_1 = PopulationRateMonitor(G[:2])
     rate_2 = PopulationRateMonitor(G[2:4])
     rate_3 = PopulationRateMonitor(G[4:])
+    rate_indexed = PopulationRateMonitor(G[::2])
     run(2*defaultclock.dt)
     assert_allclose(rate_all.rate, 0.5 / defaultclock.dt)
     assert_allclose(rate_1.rate, 0.5 / defaultclock.dt)
     assert_allclose(rate_2.rate, 0*Hz)
     assert_allclose(rate_3.rate, 1 / defaultclock.dt)
+    assert_allclose(rate_indexed.rate, 2/3*(1/defaultclock.dt))  # 2 out of 3
+
 
 if __name__ == '__main__':
     from _pytest.outcomes import Skipped

From 4c677f068c68afa0151266b8c2056048ea386e8a Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 26 Feb 2020 12:00:08 +0100
Subject: [PATCH 06/52] Add subgroups test for StateMonitor

---
 brian2/tests/test_monitor.py | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py
index fd2f0485c..6dbae66e5 100644
--- a/brian2/tests/test_monitor.py
+++ b/brian2/tests/test_monitor.py
@@ -358,6 +358,41 @@ def test_state_monitor():
     assert_allclose(synapse_mon.w[:], np.tile(S.j[:]*nS,
                                               (synapse_mon.w[:].shape[1], 1)).T)
 
+
+@pytest.mark.standalone_compatible
+def test_state_monitor_subgroups():
+    G = NeuronGroup(10, '''v : volt
+                           step_size : volt (constant)''')
+    G.run_regularly('v += step_size')
+    G.step_size = 'i*mV'
+
+    SG1 = G[3:8]
+    SG2 = G[::2]
+
+    state_mon_full = StateMonitor(G, 'v', record=True)
+
+    # monitor subgroups and record from all
+    state_mon1 = StateMonitor(SG1, 'v', record=True)
+    state_mon2 = StateMonitor(SG2, 'v', record=True)
+
+    # monitor subgroup and use (relative) indices
+    state_mon3 = StateMonitor(SG1, 'v', record=[0, 3])
+    state_mon4 = StateMonitor(SG2, 'v', record=[0, 3])
+
+    # monitor full group and use subgroup as indices
+    state_mon5 = StateMonitor(G, 'v', record=SG1)
+    state_mon6 = StateMonitor(G, 'v', record=SG2)
+
+    run(2*defaultclock.dt)
+
+    assert_allclose(state_mon1.v, state_mon_full.v[3:8])
+    assert_allclose(state_mon2.v, state_mon_full.v[::2])
+    assert_allclose(state_mon3.v, state_mon_full.v[3:8][[0, 3]])
+    assert_allclose(state_mon4.v, state_mon_full.v[::2][[0, 3]])
+    assert_allclose(state_mon5.v, state_mon_full.v[3:8])
+    assert_allclose(state_mon6.v, state_mon_full.v[::2])
+
+
 @pytest.mark.standalone_compatible
 @pytest.mark.multiple_runs
 def test_state_monitor_record_single_timestep():

From db2dcd51ca151aa545ed152703a3d79817216d59 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 26 Feb 2020 15:14:39 +0100
Subject: [PATCH 07/52] Fix subgroups of SpatialNeuron

---
 brian2/groups/subgroup.py             | 22 +---------------------
 brian2/spatialneuron/spatialneuron.py | 23 +++++++++++++++++------
 brian2/tests/test_spatialneuron.py    |  4 ++++
 3 files changed, 22 insertions(+), 27 deletions(-)

diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index c323d57e3..01bc41656 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -152,27 +152,7 @@ def __init__(self, source, start=None, stop=None, indices=None, name=None):
             if len(indices) != len(np.unique(indices)):
                 raise IndexError('sub_indices cannot contain repeated values.')
 
-        # First check if the source is itself a Subgroup
-        # If so, then make this a Subgroup of the original Group
-        if isinstance(source, Subgroup):
-            if source.contiguous:
-                if self.contiguous:
-                    source = source.source
-                    start = start + source.start
-                    stop = stop + source.start
-                else:
-                    if np.max(indices) >= source.stop - source.start:
-                        raise IndexError('Index exceeds range')
-                    indices += source.start
-            else:
-                if self.contiguous:
-                    indices = source.sub_indices[start:stop]
-                    self.contiguous = False
-                else:
-                    indices = source.sub_indices[indices]
-            self.source = source
-        else:
-            self.source = weakproxy_with_fallback(source)
+        self.source = weakproxy_with_fallback(source)
 
         # Store a reference to the source's equations (if any)
         self.equations = None
diff --git a/brian2/spatialneuron/spatialneuron.py b/brian2/spatialneuron/spatialneuron.py
index ef702bce2..97232a186 100644
--- a/brian2/spatialneuron/spatialneuron.py
+++ b/brian2/spatialneuron/spatialneuron.py
@@ -461,13 +461,19 @@ def spatialneuron_attribute(neuron, name):
         if name == 'main':  # Main section, without the subtrees
             indices = neuron.morphology.indices[:]
             start, stop = indices[0], indices[-1]
-            return SpatialSubgroup(neuron, start, stop + 1,
-                                   morphology=neuron.morphology)
+            morpho = neuron.morphology
+            if isinstance(neuron, SpatialSubgroup):
+                # For subtrees, make the new Subgroup a child of the original neuron
+                neuron = neuron.source
+            return SpatialSubgroup(neuron, start, stop + 1, morphology=morpho)
         elif (name != 'morphology') and ((name in getattr(neuron.morphology, 'children', [])) or
                                       all([c in 'LR123456789' for c in name])):  # subtree
             morpho = neuron.morphology[name]
             start = morpho.indices[0]
             stop = SpatialNeuron._find_subtree_end(morpho)
+            if isinstance(neuron, SpatialSubgroup):
+                neuron = neuron.source
+
             return SpatialSubgroup(neuron, start, stop + 1, morphology=morpho)
         else:
             return Group.__getattr__(neuron, name)
@@ -489,15 +495,20 @@ def spatialneuron_segment(neuron, item):
                 raise DimensionMismatchError('Start and stop should have units '
                                              'of meter', start, stop)
             # Convert to integers (compartment numbers)
-            indices = neuron.morphology.indices[item]
-            start, stop = indices[0], indices[-1] + 1
+            compartment_indices = neuron.morphology.indices[item]
+            start, stop = compartment_indices[0], compartment_indices[-1] + 1
+            indices = None
         elif not isinstance(item, slice) and hasattr(item, 'indices'):
             start, stop, indices = to_start_stop_or_index(item.indices[:],
                                                           neuron._N)
         else:
             start, stop, indices = to_start_stop_or_index(item, neuron._N)
 
-        return Subgroup(neuron, start, stop)
+        if isinstance(neuron, SpatialSubgroup):
+            # For subtrees, make the new Subgroup a child of the original neuron
+            neuron = neuron.source
+
+        return Subgroup(neuron, start, stop, indices)
 
 
 class SpatialSubgroup(Subgroup):
@@ -519,7 +530,7 @@ class SpatialSubgroup(Subgroup):
 
     def __init__(self, source, start, stop, morphology, name=None):
         self.morphology = morphology
-        Subgroup.__init__(self, source, start, stop, name)
+        Subgroup.__init__(self, source, start, stop, name=name)
 
     def __getattr__(self, name):
         return SpatialNeuron.spatialneuron_attribute(self, name)
diff --git a/brian2/tests/test_spatialneuron.py b/brian2/tests/test_spatialneuron.py
index 56a9b6aea..50bb6eb9a 100644
--- a/brian2/tests/test_spatialneuron.py
+++ b/brian2/tests/test_spatialneuron.py
@@ -659,6 +659,10 @@ def test_spatialneuron_indexing():
     assert len(neuron[0:1].indices[:]) == 1
     assert len(neuron[sec.sec2.indices[:]]) == 16
     assert len(neuron[sec.sec2]) == 16
+    assert len(neuron[:16:2].indices[:]) == 8
+    assert len(neuron[:8][4:].indices[:]) == 4
+    assert len(neuron[:8][::2].indices[:]) == 4
+    assert len(neuron[::2][:8].indices[:]) == 8
 
 
 @pytest.mark.codegen_independent

From aa9237d0e5188b2a97a3b862493f29036e362224 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 26 Feb 2020 16:08:22 +0100
Subject: [PATCH 08/52] Forbid using non-contiguous subgroups in Synapses

---
 brian2/synapses/synapses.py   |  7 +++++++
 brian2/tests/test_synapses.py | 11 +++++++++--
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 5d5ef523d..efb09bb91 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -22,6 +22,7 @@
                                         PARAMETER, INTEGER,
                                         check_subexpressions)
 from brian2.groups.group import Group, CodeRunner, get_dtype
+from brian2.groups.subgroup import Subgroup
 from brian2.groups.neurongroup import (extract_constant_subexpressions,
                                        SubexpressionUpdater,
                                        check_identifier_pre_post)
@@ -1097,6 +1098,9 @@ def _create_variables(self, equations, user_dtype=None):
         if getattr(self.source, 'start', 0) == 0:
             self.variables.add_reference('i', self, '_synaptic_pre')
         else:
+            if isinstance(self.source, Subgroup) and not self.source.contiguous:
+                raise TypeError('Cannot use a non-contiguous subgroup as a '
+                                'source group for Synapses.')
             self.variables.add_reference('_source_i', self.source.source, 'i',
                                          index='_presynaptic_idx')
             self.variables.add_reference('_source_offset', self.source, '_offset')
@@ -1107,6 +1111,9 @@ def _create_variables(self, equations, user_dtype=None):
         if getattr(self.target, 'start', 0) == 0:
             self.variables.add_reference('j', self, '_synaptic_post')
         else:
+            if isinstance(self.target, Subgroup) and not self.target.contiguous:
+                raise TypeError('Cannot use a non-contiguous subgroup as a '
+                                'target group for Synapses.')
             self.variables.add_reference('_target_j', self.target.source, 'i',
                                          index='_postsynaptic_idx')
             self.variables.add_reference('_target_offset', self.target, '_offset')
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index be3da1ad3..5b2bb4a43 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -75,9 +75,16 @@ def test_creation_errors():
     # Check that using pre and on_pre (resp. post/on_post) at the same time
     # raises an error
     with pytest.raises(TypeError):
-        Synapses(G, G, 'w:1', pre='v+=w', on_pre='v+=w', connect=True)
+        Synapses(G, G, 'w:1', pre='v+=w', on_pre='v+=w')
     with pytest.raises(TypeError):
-        Synapses(G, G, 'w:1', post='v+=w', on_post='v+=w', connect=True)
+        Synapses(G, G, 'w:1', post='v+=w', on_post='v+=w')
+    # We do not allow non-contiguous subgroups as source/target groups at the
+    # moment
+    with pytest.raises(TypeError):
+        Synapses(G[::2], G, '')
+    with pytest.raises(TypeError):
+        Synapses(G, G[::2], '')
+
 
 
 @pytest.mark.codegen_independent

From db3ddf4611289b99e0989a7ee5952b82b1084c3b Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 26 Feb 2020 18:26:06 +0100
Subject: [PATCH 09/52] Enable string expressions for subgroups

Needs more testing
---
 brian2/groups/group.py                |  28 +++++--
 brian2/groups/neurongroup.py          |   2 +-
 brian2/groups/subgroup.py             | 111 +++++++++-----------------
 brian2/spatialneuron/spatialneuron.py |   4 +-
 brian2/tests/test_subgroup.py         |  16 +++-
 dev/tools/run_tests.py                |   2 +-
 6 files changed, 78 insertions(+), 85 deletions(-)

diff --git a/brian2/groups/group.py b/brian2/groups/group.py
index b1315328d..f21158846 100644
--- a/brian2/groups/group.py
+++ b/brian2/groups/group.py
@@ -232,6 +232,7 @@ def __call__(self, item=slice(None), index_var=None):
         --------
         SynapticIndexing
         '''
+
         if index_var is None:
             index_var = self.default_index
 
@@ -251,7 +252,7 @@ def __call__(self, item=slice(None), index_var=None):
                     start, stop, step = item.indices(int(self.N.get_value()))
                 else:
                     start, stop, step = item.indices(index_var.size)
-                index_array = np.arange(start, stop, step, dtype=np.int32)
+                index_array = np.arange(start, stop, step)
             else:
                 index_array = np.asarray(item)
                 if index_array.dtype == np.bool:
@@ -261,9 +262,11 @@ def __call__(self, item=slice(None), index_var=None):
                                      'and boolean arrays, not for type '
                                      '%s' % index_array.dtype))
 
+            if index_array.size == 0:
+                return index_array
             if index_var != '_idx':
                 try:
-                    return index_var.get_value()[index_array]
+                    index_array = index_var.get_value()[index_array]
                 except IndexError as ex:
                     # We try to emulate numpy's indexing semantics here:
                     # slices never lead to IndexErrors, instead they return an
@@ -273,7 +276,17 @@ def __call__(self, item=slice(None), index_var=None):
                     else:
                         raise ex
             else:
-                return index_array
+                N = int(self.N.get_value())
+                if np.min(index_array) < -N:
+                    raise IndexError('Illegal index {} for a group of size '
+                                     '{}'.format(np.min(index_array), N))
+                if np.max(index_array) >= N:
+                    raise IndexError('Illegal index {} for a group of size '
+                                     '{}'.format(np.max(index_array), N))
+
+                index_array = index_array % N  # interpret negative indices
+
+            return index_array
 
 
 class IndexWrapper(object):
@@ -288,13 +301,16 @@ def __init__(self, group):
         self.indices = group._indices
 
     def __getitem__(self, item):
+        return self.get_item(item, level=1)
+
+    def get_item(self, 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=1)
+            namespace = get_local_namespace(level=level+1)
             from brian2.devices.device import get_device
             device = get_device()
             codeobj = create_runner_codeobj(self.group,
@@ -304,7 +320,9 @@ def __getitem__(self, item):
                                             additional_variables=variables,
                                             codeobj_class=device.code_object_class(fallback_pref='codegen.string_expression_target')
                                             )
-            return codeobj()
+            indices = codeobj()
+            # Handle subgroups correctly
+            return self.indices(indices)
         else:
             return self.indices(item)
 
diff --git a/brian2/groups/neurongroup.py b/brian2/groups/neurongroup.py
index 988c00661..4024b8690 100644
--- a/brian2/groups/neurongroup.py
+++ b/brian2/groups/neurongroup.py
@@ -729,7 +729,7 @@ def __setattr__(self, key, value):
                 Group.__setattr__(self, key, value, level=1)
 
     def __getitem__(self, item):
-        start, stop, indices = to_start_stop_or_index(item, self._N)
+        start, stop, indices = to_start_stop_or_index(item, self, level=1)
         return Subgroup(self, start, stop, indices)
 
     def _create_variables(self, user_dtype, events):
diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index 01bc41656..c80cd91ec 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -16,7 +16,7 @@
 logger = get_logger(__name__)
 
 
-def to_start_stop_or_index(item, N):
+def to_start_stop_or_index(item, group, level=0):
     '''
     Helper function to transform a single number, a slice or an array of
     indices to a start and stop value (if possible), or to an index of positive
@@ -26,11 +26,11 @@ def to_start_stop_or_index(item, N):
 
     Parameters
     ----------
-    item : slice, int or sequence
-        The slice, index, or sequence of indices to use.
-    N : int
-        The total number of elements in the group.
-
+    item : slice, int, str, or sequence
+        The slice, index, or sequence of indices to use, or a boolean string
+        expression that can be evaluated in the context of the group.
+    group : `Group`
+        The group providing the context for the interpretation.
     Returns
     -------
     start : int or None
@@ -42,73 +42,41 @@ def to_start_stop_or_index(item, N):
 
     Examples
     --------
-    >>> from brian2.groups.neurongroup import to_start_stop_or_index
-    >>> to_start_stop_or_index(slice(3, 6), 10)
+    >>> from brian2.groups.neurongroup import NeuronGroup, to_start_stop_or_index
+    >>> group = NeuronGroup(10, '')
+    >>> to_start_stop_or_index(slice(3, 6), group)
     (3, 6, None)
-    >>> to_start_stop_or_index(slice(3, None), 10)
+    >>> to_start_stop_or_index(slice(3, None), group)
     (3, 10, None)
-    >>> to_start_stop_or_index(5, 10)
+    >>> to_start_stop_or_index(5, group)
     (5, 6, None)
-    >>> to_start_stop_or_index(slice(None, None, 2), 10)
+    >>> to_start_stop_or_index(slice(None, None, 2), group)
     (None, None, array([0, 2, 4, 6, 8]))
-    >>> to_start_stop_or_index([3, 4, 5], 10)
+    >>> to_start_stop_or_index([3, 4, 5], group)
     (3, 6, None)
-    >>> to_start_stop_or_index([3, 5, 7], 10)
+    >>> to_start_stop_or_index([3, 5, 7], group)
     (None, None, array([3, 5, 7]))
-    >>> to_start_stop_or_index([-1, -2, -3], 10)
-    (None, None, array([7, 8, 9]))
+    >>> to_start_stop_or_index([-1, -2, -3], group)
+    (7, 10, None)
     '''
-    start = stop = indices = None
-    if isinstance(item, slice):
-        start, stop, step = item.indices(N)
-        if step != 1:
-            indices = np.arange(start, stop, step)
-            start = stop = None
-    elif isinstance(item, numbers.Integral):
-        start = item
-        stop = item + 1
-        step = 1
-    elif (isinstance(item, (Sequence, np.ndarray)) and
-          not isinstance(item, str)):
-        item = np.asarray(item)
-        if not np.issubdtype(item.dtype, np.integer):
-            raise TypeError('Subgroups can only be constructed using integer '
-                            'values.')
-        if not len(item) > 0:
-            raise IndexError('Cannot create an empty subgroup')
-        if np.min(item) < -N:
-            raise IndexError('Illegal index {} for a group of size '
-                             '{}'.format(np.min(item), N))
-        elif np.max(item) >= N:
-            raise IndexError('Illegal index {} for a group of size '
-                             '{}'.format(np.min(item), N))
-        indices = item % N  # interpret negative indices correctly
-
-        if not np.all(item[:-1] <= item[1:]):
-            logger.warn('The indices provided to create the subgroup were '
-                        'not sorted. They will be sorted before use.',
-                        name_suffix='unsorted_subgroup_indices')
-            indices.sort()
-
-        if np.all(np.diff(item) == 1):
-            start = int(item[0])
-            stop = int(item[-1]) + 1
-            indices = None
-
-    else:
-        raise TypeError("Cannot interpret object of type '{}' as index or "
-                        "slice".format(type(item)))
-
-    if indices is None:
-        if start >= stop:
-            raise IndexError('Illegal start/end values for subgroup, %d>=%d' %
-                             (start, stop))
-        if start >= N:
-            raise IndexError('Illegal start value for subgroup, %d>=%d' %
-                             (start, N))
-        if stop > N:
-            raise IndexError('Illegal stop value for subgroup, %d>%d' %
-                             (stop, N))
+    start = stop = None
+    indices = group.indices.get_item(item, level=level+1)
+    # For convenience, allow subgroups with a single value instead of x:x+1 slice
+    if indices.shape == ():
+        indices = np.array([indices])
+
+    if not np.all(indices[:-1] <= indices[1:]):
+        logger.warn('The indices provided to create the subgroup were '
+                    'not sorted. They will be sorted before use.',
+                    name_suffix='unsorted_subgroup_indices')
+        indices.sort()
+    if not len(indices) > 0:
+        raise IndexError('Cannot create an empty subgroup')
+
+    if np.all(np.diff(indices) == 1):
+        start = int(indices[0])
+        stop = int(indices[-1]) + 1
+        indices = None
 
     return start, stop, indices
 
@@ -242,15 +210,8 @@ def __init__(self, source, start=None, stop=None, indices=None, name=None):
     spikes = property(lambda self: self.source.spikes)
 
     def __getitem__(self, item):
-        start, stop, indices = to_start_stop_or_index(item, self._N)
-        if self.contiguous:
-            if start is not None:
-                return Subgroup(self.source, self.start + start, self.start + stop)
-            else:
-                return Subgroup(self.source, indices=indices+self.start)
-        else:
-            indices = self.sub_indices[item]
-            return Subgroup(self.source, indices=indices)
+        start, stop, indices = to_start_stop_or_index(item, self, level=1)
+        return Subgroup(self.source, start, stop, indices)
 
     def __repr__(self):
         if self.contiguous:
diff --git a/brian2/spatialneuron/spatialneuron.py b/brian2/spatialneuron/spatialneuron.py
index 97232a186..f870dbba5 100644
--- a/brian2/spatialneuron/spatialneuron.py
+++ b/brian2/spatialneuron/spatialneuron.py
@@ -500,9 +500,9 @@ def spatialneuron_segment(neuron, item):
             indices = None
         elif not isinstance(item, slice) and hasattr(item, 'indices'):
             start, stop, indices = to_start_stop_or_index(item.indices[:],
-                                                          neuron._N)
+                                                          neuron)
         else:
-            start, stop, indices = to_start_stop_or_index(item, neuron._N)
+            start, stop, indices = to_start_stop_or_index(item, neuron)
 
         if isinstance(neuron, SpatialSubgroup):
             # For subtrees, make the new Subgroup a child of the original neuron
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 42cf6ff20..760fa48a9 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -153,6 +153,20 @@ def test_state_variables_group_as_index_problematic():
             assert all([entry[1].endswith('ambiguous_string_expression')
                         for entry in l])
 
+@pytest.mark.standalone_compatible
+def test_state_variables_string_group():
+    G = NeuronGroup(10, 'v : 1')
+    G.v = 'i'
+    c = 3
+    SG1 = G['i > 5']
+    SG2 = G['v > c']
+    SG1.v = 'v * 2'
+    run(0*ms)  # for standalone
+    assert_equal(G.v[:], [0, 1, 2, 3, 4, 5, 12, 14, 16, 18])
+    assert_equal(SG1.v[:], [12, 14, 16, 18])
+    assert_equal(SG2.v[:], [4, 5, 12, 14, 16, 18])
+
+
 @pytest.mark.standalone_compatible
 def test_state_monitor():
     G = NeuronGroup(10, 'v : volt')
@@ -669,7 +683,7 @@ def test_spike_monitor():
 
 
 @pytest.mark.codegen_independent
-@pytest.mark.parametrize('item', ['string', slice(10, None), slice(3, 2),
+@pytest.mark.parametrize('item', [slice(10, None), slice(3, 2),
                                   [9, 10], [10, 11], [2.5, 3.5, 4.5],
                                   [5, 5, 5], []])
 def test_wrong_indexing(item):
diff --git a/dev/tools/run_tests.py b/dev/tools/run_tests.py
index 3b5ea6875..7140b52b7 100644
--- a/dev/tools/run_tests.py
+++ b/dev/tools/run_tests.py
@@ -6,5 +6,5 @@
 import brian2
 
 if __name__ == '__main__':
-    if not brian2.test():  # If the test fails, exit with a non-zero error code
+    if not brian2.test(additional_args=['-x']):  # If the test fails, exit with a non-zero error code
         sys.exit(1)

From f0a34d1de361438ee506b6a077f565f87372adcf Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 27 Feb 2020 13:23:35 +0100
Subject: [PATCH 10/52] Make string expressions for subgroups/indexing work in
 standalone

---
 brian2/devices/cpp_standalone/device.py | 36 +++++++++++++++++++++++--
 brian2/groups/group.py                  |  4 ++-
 brian2/tests/test_subgroup.py           |  4 +--
 3 files changed, 38 insertions(+), 6 deletions(-)

diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index 2dc4885e3..f82b84a72 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/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
@@ -470,7 +471,8 @@ def variableview_set_with_index_array(self, variableview, item,
                                                       value)
             # Put values into the cache
             cache_variable = self.array_cache[variableview.variable]
-            cache_variable[indices] = value
+            if cache_variable is not None:
+                cache_variable[indices] = value
             self.main_queue.append(('set_array_by_array', (arrayname,
                                                            staticarrayname_index,
                                                            staticarrayname_value)))
@@ -537,6 +539,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)
diff --git a/brian2/groups/group.py b/brian2/groups/group.py
index f21158846..a717003d1 100644
--- a/brian2/groups/group.py
+++ b/brian2/groups/group.py
@@ -14,7 +14,8 @@
 
 import numpy as np
 
-from brian2.core.base import BrianObject, weakproxy_with_fallback
+from brian2.core.base import (BrianObject, weakproxy_with_fallback,
+                              device_override)
 from brian2.core.names import Nameable
 from brian2.core.preferences import prefs
 from brian2.core.variables import (Variables, Constant, Variable,
@@ -303,6 +304,7 @@ def __init__(self, group):
     def __getitem__(self, item):
         return self.get_item(item, level=1)
 
+    @device_override('index_wrapper_get_item')
     def get_item(self, item, level):
         if isinstance(item, str):
             variables = Variables(None)
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 760fa48a9..5d49d53df 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -158,13 +158,11 @@ def test_state_variables_string_group():
     G = NeuronGroup(10, 'v : 1')
     G.v = 'i'
     c = 3
-    SG1 = G['i > 5']
-    SG2 = G['v > c']
+    SG1 = G['i > 5']  # indexing with constant expressions should even work in standalone
     SG1.v = 'v * 2'
     run(0*ms)  # for standalone
     assert_equal(G.v[:], [0, 1, 2, 3, 4, 5, 12, 14, 16, 18])
     assert_equal(SG1.v[:], [12, 14, 16, 18])
-    assert_equal(SG2.v[:], [4, 5, 12, 14, 16, 18])
 
 
 @pytest.mark.standalone_compatible

From 5acd9175c1754b137f64ccb7886050cfa8a1531f Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 28 Feb 2020 18:10:51 +0100
Subject: [PATCH 11/52] Fix string indexing with references to external
 constants

---
 brian2/groups/group.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/brian2/groups/group.py b/brian2/groups/group.py
index a717003d1..36a9f749a 100644
--- a/brian2/groups/group.py
+++ b/brian2/groups/group.py
@@ -312,7 +312,7 @@ def get_item(self, item, level):
             variables.add_auxiliary_variable('_cond', dtype=np.bool)
 
             abstract_code = '_cond = ' + item
-            namespace = get_local_namespace(level=level+1)
+            namespace = get_local_namespace(level=level+2)  # decorated function
             from brian2.devices.device import get_device
             device = get_device()
             codeobj = create_runner_codeobj(self.group,

From 33f537e9139aebb612d656330061b0f535ad6e80 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 14 Sep 2020 10:51:58 +0200
Subject: [PATCH 12/52] Make doctest compatible with Windows output for 32bit
 ints

---
 brian2/groups/subgroup.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index c80cd91ec..7d17b0043 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -50,12 +50,12 @@ def to_start_stop_or_index(item, group, level=0):
     (3, 10, None)
     >>> to_start_stop_or_index(5, group)
     (5, 6, None)
-    >>> to_start_stop_or_index(slice(None, None, 2), group)
-    (None, None, array([0, 2, 4, 6, 8]))
+    >>> to_start_stop_or_index(slice(None, None, 2), group)  # doctest: +ELLIPSIS
+    (None, None, array([0, 2, 4, 6, 8]...))
     >>> to_start_stop_or_index([3, 4, 5], group)
     (3, 6, None)
-    >>> to_start_stop_or_index([3, 5, 7], group)
-    (None, None, array([3, 5, 7]))
+    >>> to_start_stop_or_index([3, 5, 7], group)  # doctest: +ELLIPSIS
+    (None, None, array([3, 5, 7]...))
     >>> to_start_stop_or_index([-1, -2, -3], group)
     (7, 10, None)
     '''

From b47758395a52547f808f70fe1a5bfca35cd5f503 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 14 Sep 2020 16:27:57 +0200
Subject: [PATCH 13/52] Do not allow unsorted indices in subgroups

Instead of raising a warning
---
 brian2/groups/subgroup.py     | 7 ++-----
 brian2/tests/test_monitor.py  | 8 ++------
 brian2/tests/test_subgroup.py | 3 +--
 3 files changed, 5 insertions(+), 13 deletions(-)

diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index 7d17b0043..214201f00 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -56,7 +56,7 @@ def to_start_stop_or_index(item, group, level=0):
     (3, 6, None)
     >>> to_start_stop_or_index([3, 5, 7], group)  # doctest: +ELLIPSIS
     (None, None, array([3, 5, 7]...))
-    >>> to_start_stop_or_index([-1, -2, -3], group)
+    >>> to_start_stop_or_index([-3, -2, -1], group)
     (7, 10, None)
     '''
     start = stop = None
@@ -66,10 +66,7 @@ def to_start_stop_or_index(item, group, level=0):
         indices = np.array([indices])
 
     if not np.all(indices[:-1] <= indices[1:]):
-        logger.warn('The indices provided to create the subgroup were '
-                    'not sorted. They will be sorted before use.',
-                    name_suffix='unsorted_subgroup_indices')
-        indices.sort()
+        raise TypeError('Subgroups can only be created from ordered indices.')
     if not len(indices) > 0:
         raise IndexError('Cannot create an empty subgroup')
 
diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py
index 6dbae66e5..6b2f4ff6b 100644
--- a/brian2/tests/test_monitor.py
+++ b/brian2/tests/test_monitor.py
@@ -133,9 +133,8 @@ def test_spike_monitor_subgroups():
     spikes_2 = SpikeMonitor(G[2:4])
     spikes_3 = SpikeMonitor(G[4:])
     spikes_indexed = SpikeMonitor(G[::2])
-    with catch_logs() as l:
-        spikes_indexed_unsorted = SpikeMonitor(G[[4, 0, 2]])  # This should not make a difference
-    assert len(l) == 1 and l[0][1].endswith('.unsorted_subgroup_indices')
+    with pytest.raises(TypeError):
+        SpikeMonitor(G[[4, 0, 2]])  # unsorted
     run(defaultclock.dt)
     # Spikes
     assert_allclose(spikes_all.i, [0, 4, 5])
@@ -148,15 +147,12 @@ def test_spike_monitor_subgroups():
     assert_allclose(spikes_3.t, [0, 0] * ms)
     assert_allclose(spikes_indexed.i, [0, 2])
     assert_allclose(spikes_indexed.t, [0, 0]*ms)
-    assert_allclose(spikes_indexed_unsorted.i, [0, 2])
-    assert_allclose(spikes_indexed_unsorted.t, [0, 0]*ms)
     # Spike count
     assert_allclose(spikes_all.count, [1, 0, 0, 0, 1, 1])
     assert_allclose(spikes_1.count, [1, 0])
     assert_allclose(spikes_2.count, [0, 0])
     assert_allclose(spikes_3.count, [1, 1])
     assert_allclose(spikes_indexed.count, [1, 0, 1])
-    assert_allclose(spikes_indexed_unsorted.count, [1, 0, 1])
 
 
 def test_spike_monitor_bug_824():
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 5d49d53df..162afe2fb 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -683,7 +683,7 @@ def test_spike_monitor():
 @pytest.mark.codegen_independent
 @pytest.mark.parametrize('item', [slice(10, None), slice(3, 2),
                                   [9, 10], [10, 11], [2.5, 3.5, 4.5],
-                                  [5, 5, 5], []])
+                                  [5, 5, 5], [], [5, 4, 3]])
 def test_wrong_indexing(item):
     G = NeuronGroup(10, 'v:1')
     with pytest.raises((TypeError, IndexError)):
@@ -696,7 +696,6 @@ def test_wrong_indexing(item):
                                            (3, np.array([3])),
                                            ([3, 4, 5], np.array([3, 4, 5])),
                                            ([3, 5, 7], np.array([3, 5, 7])),
-                                           ([7, 5, 3], np.array([3, 5, 7])),
                                            ([3, -1], np.array([3, 9]))])
 def test_alternative_indexing(item, expected):
     G = NeuronGroup(10, 'v : integer')

From cdb5427b1b5262243d18d81c16d0cfd750730ff0 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 21 Apr 2021 18:36:28 +0200
Subject: [PATCH 14/52] Fix a remaining case of a deprecated numpy scalar

---
 brian2/parsing/bast.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/brian2/parsing/bast.py b/brian2/parsing/bast.py
index a456ba5bc..d792b9f9c 100644
--- a/brian2/parsing/bast.py
+++ b/brian2/parsing/bast.py
@@ -57,7 +57,7 @@ def brian_dtype_from_value(value):
 # so we don't use numpy.issubdtype but instead a precalculated list of all
 # standard types
 
-bool_dtype =numpy.dtype(numpy.bool)
+bool_dtype =numpy.dtype(bool)
 def is_boolean_dtype(obj):
     return numpy.dtype(obj) is bool_dtype
 

From 89430b9814c64de86e30202d310a4f43fe6ea084 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 2 May 2022 11:49:16 +0200
Subject: [PATCH 15/52] minor optimization for non-contiguous subgroups

---
 brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx | 2 ++
 brian2/devices/cpp_standalone/templates/spikemonitor.cpp    | 3 +++
 2 files changed, 5 insertions(+)

diff --git a/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx b/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
index 78a74e59f..1017a671f 100644
--- a/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
+++ b/brian2/codegen/runtime/cython_rt/templates/spikemonitor.pyx
@@ -41,6 +41,8 @@
             _idx = {{_eventspace}}[_j]
             if _idx < {{_source_indices}}[_source_index_counter]:
                 continue
+            if _idx > {{_source_indices}}[{{source_N}}-1]:
+                break
             while {{_source_indices}}[_source_index_counter] < _idx:
                 _source_index_counter += 1
             if (_source_index_counter < {{source_N}} and
diff --git a/brian2/devices/cpp_standalone/templates/spikemonitor.cpp b/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
index b94c61cc3..96ec3e41d 100644
--- a/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
+++ b/brian2/devices/cpp_standalone/templates/spikemonitor.cpp
@@ -40,11 +40,14 @@
         }
         _num_events = _end_idx - _start_idx;
         {% else %}
+        const size_t _max_source_index = {{_source_indices}}[{{source_N}}-1];
         for (size_t _j=0; _j<_num_events; _j++)
         {
             const size_t _idx = {{_eventspace}}[_j];
             if (_idx < {{_source_indices}}[_source_index_counter])
                 continue;
+            if (_idx > _max_source_index)
+                break;
             while ({{_source_indices}}[_source_index_counter] < _idx)
             {
                 _source_index_counter++;

From 269d9e046f1841a1ed98c62c2a0f47aaad5303cd Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 2 May 2022 11:49:38 +0200
Subject: [PATCH 16/52] make PoissonGroup work with non-contiguous subgroups

---
 brian2/input/poissongroup.py | 13 +++----------
 1 file changed, 3 insertions(+), 10 deletions(-)

diff --git a/brian2/input/poissongroup.py b/brian2/input/poissongroup.py
index 7d70bbc4c..afa49f398 100644
--- a/brian2/input/poissongroup.py
+++ b/brian2/input/poissongroup.py
@@ -11,7 +11,7 @@
                                            fail_for_dimension_mismatch)
 from brian2.units.stdunits import Hz
 from brian2.groups.group import Group
-from brian2.groups.subgroup import Subgroup
+from brian2.groups.subgroup import Subgroup, to_start_stop_or_index
 from brian2.groups.neurongroup import Thresholder
 from brian2.utils.stringtools import get_identifiers
 
@@ -97,15 +97,8 @@ def __init__(self, N, rates, dt=None, clock=None, when='thresholds',
             self.rates = rates
 
     def __getitem__(self, item):
-        if not isinstance(item, slice):
-            raise TypeError("Subgroups can only be constructed using slicing syntax")
-        start, stop, step = item.indices(self._N)
-        if step != 1:
-            raise IndexError("Subgroups have to be contiguous")
-        if start >= stop:
-            raise IndexError(f"Illegal start/end values for subgroup, {int(start)}>={int(stop)}")
-
-        return Subgroup(self, start, stop)
+        start, stop, indices = to_start_stop_or_index(item, self, level=1)
+        return Subgroup(self, start, stop, indices)
 
     def before_run(self, run_namespace=None):
         rates_var = self.variables['rates']

From 73cd6d996962dbf3197f411208e81e9e2a0b1e4c Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 2 May 2022 11:50:19 +0200
Subject: [PATCH 17/52] make SpikeMonitor generated code deterministic

---
 brian2/monitors/spikemonitor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/brian2/monitors/spikemonitor.py b/brian2/monitors/spikemonitor.py
index 86871099a..41148b415 100644
--- a/brian2/monitors/spikemonitor.py
+++ b/brian2/monitors/spikemonitor.py
@@ -112,7 +112,7 @@ def __init__(self, source, event, variables=None, record=True,
         # Some dummy code so that code generation takes care of the indexing
         # and subexpressions
         code = [f'_to_record_{v} = _source_{v}'
-                for v in self.record_variables]
+                for v in sorted(self.record_variables)]
         code = '\n'.join(code)
 
         self.codeobj_class = codeobj_class

From d0abc65b14694e28caf01fb6d6bbb23e44046316 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 2 May 2022 14:55:08 +0200
Subject: [PATCH 18/52] Check for illegal subgroup indices in consistent way

---
 brian2/groups/subgroup.py | 31 ++++++++++++++++++-------------
 1 file changed, 18 insertions(+), 13 deletions(-)

diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index e8432d82f..68f13c7db 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -64,11 +64,6 @@ def to_start_stop_or_index(item, group, level=0):
     if indices.shape == ():
         indices = np.array([indices])
 
-    if not np.all(indices[:-1] <= indices[1:]):
-        raise TypeError('Subgroups can only be created from ordered indices.')
-    if not len(indices) > 0:
-        raise IndexError('Cannot create an empty subgroup')
-
     if np.all(np.diff(indices) == 1):
         start = int(indices[0])
         stop = int(indices[-1]) + 1
@@ -97,24 +92,34 @@ class Subgroup(Group, SpikeSource):
     """
     def __init__(self, source, start=None, stop=None, indices=None, name=None):
         if start is stop is indices is None:
-            raise TypeError('Need to specify either start and stop or indices.')
+            raise TypeError("Need to specify either start and stop or indices.")
         if start != stop and (start is None or stop is None):
-            raise TypeError('start and stop have to be specified together.')
+            raise TypeError("start and stop have to be specified together.")
         if indices is not None and (start is not None):
-            raise TypeError('Cannot specify both indices and start and '
-                             'stop.')
+            raise TypeError("Cannot specify both indices and start and "
+                            "stop.")
         if start is not None:
             self.contiguous = True
+            if start < 0:
+                raise IndexError("Start index cannot be negative.")
+            if stop <= start:
+                raise IndexError("Stop index has to be bigger than start.")
+            if stop > len(source):
+                raise IndexError(f"Stop index cannot be > the size of the group "
+                                 f"({stop} > {len(source)}).")
         else:
             self.contiguous = False
             if not len(indices):
                 raise IndexError('Cannot create an empty subgroup.')
+            min_index = np.min(indices)
             max_index = np.max(indices)
+            if min_index < 0:
+                raise IndexError("Indices cannot contain negative values.")
             if max_index >= len(source):
-                raise IndexError('Index {} cannot be >= the size of the group '
-                                '({})'.format(max_index, len(source)))
-            if len(indices) != len(np.unique(indices)):
-                raise IndexError('indices cannot contain repeated values.')
+                raise IndexError(f"Indices cannot be ≥ the size of the group "
+                                 f"({max_index} ≥ {len(source)}).")
+            if not np.all(np.diff(indices) > 0):
+                raise IndexError("indices need to be sorted and cannot contain repeated values.")
 
         # A Subgroup should never be constructed from another Subgroup
         # Instead, use Subgroup(source.source,

From 20e0efbe877d6099ff5cbfce71cf224fdd269d14 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 2 May 2022 15:21:00 +0200
Subject: [PATCH 19/52] revert run_tests script to original version

---
 dev/tools/run_tests.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dev/tools/run_tests.py b/dev/tools/run_tests.py
index 0a74d6a59..61fdfbca8 100644
--- a/dev/tools/run_tests.py
+++ b/dev/tools/run_tests.py
@@ -6,5 +6,5 @@
 import brian2
 
 if __name__ == '__main__':
-    if not brian2.test(additional_args=['-x']):  # If the test fails, exit with a non-zero error code
+    if not brian2.test():  # If the test fails, exit with a non-zero error code
         sys.exit(1)

From 3ae9e1afe0887fe58621946623e841240de7e06b Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 2 May 2022 16:40:36 +0200
Subject: [PATCH 20/52] Add/update tests for subgroups

---
 brian2/tests/test_monitor.py  |  2 +-
 brian2/tests/test_subgroup.py | 60 +++++++++++++++++++++++++++++++++++
 2 files changed, 61 insertions(+), 1 deletion(-)

diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py
index 4eec2e6d0..ff07e0e34 100644
--- a/brian2/tests/test_monitor.py
+++ b/brian2/tests/test_monitor.py
@@ -138,7 +138,7 @@ def test_spike_monitor_subgroups():
     spikes_2 = SpikeMonitor(G[2:4])
     spikes_3 = SpikeMonitor(G[4:])
     spikes_indexed = SpikeMonitor(G[::2])
-    with pytest.raises(TypeError):
+    with pytest.raises(IndexError):
         SpikeMonitor(G[[4, 0, 2]])  # unsorted
     run(defaultclock.dt)
     # Spikes
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index bf0b7e6de..e9088a85f 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -24,6 +24,35 @@ def test_str_repr():
     assert len(str(SGi))
     assert len(repr(SGi))
 
+@pytest.mark.codegen_independent
+def test_creation():
+    G = NeuronGroup(10, '')
+    SG = G[5:8]
+    SGi = G[[3, 6, 9]]
+    SGi2 = G[[3, 4, 5]]
+    assert len(SG) == 3
+    assert SG.contiguous
+    assert len(SGi) == 3
+    assert not SGi.contiguous
+    assert len(SGi2) == 3
+    assert SGi2.contiguous
+
+    with pytest.raises(IndexError):
+        Subgroup(G, start=-1, stop=5)
+    with pytest.raises(IndexError):
+        Subgroup(G, start=1, stop=1)
+    with pytest.raises(IndexError):
+        Subgroup(G, start=1, stop=11)
+    with pytest.raises(IndexError):
+        Subgroup(G, indices=[])
+    with pytest.raises(IndexError):
+        Subgroup(G, indices=[1, 1, 2])
+    with pytest.raises(IndexError):
+        Subgroup(G, indices=[-1, 1, 2])
+    with pytest.raises(IndexError):
+        Subgroup(G, indices=[1, 2, 10])
+    with pytest.raises(IndexError):
+        Subgroup(G, indices=[3, 2, 1])
 
 def test_state_variables():
     """
@@ -95,6 +124,13 @@ def test_state_variables_simple_indexed():
                            c : 1
                            d : 1
                            ''')
+    # Illegal indices:
+    with pytest.raises(IndexError):
+        G[[3, 5, 5, 7, 9]]  # duplicate indices
+    with pytest.raises(IndexError):
+        G[[9, 7, 5, 3]]  # unsorted
+    with pytest.raises(IndexError):
+        G[[8, 10]]  # out of range
     SG = G[[3, 5, 7, 9]]
     SG.a = 1
     SG.a['i == 0'] = 2
@@ -170,16 +206,40 @@ def test_state_monitor():
     G = NeuronGroup(10, 'v : volt')
     G.v = np.arange(10) * volt
     SG = G[5:]
+    SG_indices = G[[2, 4, 6]]
     mon_all = StateMonitor(SG, 'v', record=True)
     mon_0 = StateMonitor(SG, 'v', record=0)
+    mon_all_indices = StateMonitor(SG_indices, 'v', record=True)
+    mon_0_indices = StateMonitor(SG_indices, 'v', record=0)
     run(defaultclock.dt)
 
     assert_allclose(mon_0[0].v, mon_all[0].v)
     assert_allclose(mon_0[0].v, np.array([5]) * volt)
     assert_allclose(mon_all.v.flatten(), np.arange(5, 10) * volt)
 
+    assert_allclose(mon_0_indices[0].v, mon_all_indices[0].v)
+    assert_allclose(mon_0_indices[0].v, np.array([2]) * volt)
+    assert_allclose(mon_all_indices.v.flatten(), np.array([2, 4, 6]) * volt)
+
     with pytest.raises(IndexError):
         mon_all[5]
+    with pytest.raises(IndexError):
+        mon_all_indices[3]
+
+@pytest.mark.standalone_compatible
+def test_rate_monitor():
+    G = NeuronGroup(10, '', threshold='i%2 == 0')
+    SG = G[5:]
+    SG_indices = G[[2, 4, 6]]
+    SG_indices2 = G[[3, 5, 7]]
+    pop_mon = PopulationRateMonitor(SG)
+    pop_mon_indices = PopulationRateMonitor(SG_indices)
+    pop_mon_indices2 = PopulationRateMonitor(SG_indices2)
+    run(defaultclock.dt)
+    r = 1/defaultclock.dt
+    assert_allclose(pop_mon.rate[0], (2*r + 3*0)/5)  # 2 out of 3 neurons are firing
+    assert_allclose(pop_mon_indices.rate[0], r)      # all neurons firing
+    assert_allclose(pop_mon_indices2.rate[0], 0*Hz)  # no neurons firing
 
 
 def test_shared_variable():

From d8707aa00f34a289b549ccdcdd8585edeee6bc34 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 10 Nov 2022 15:06:22 +0100
Subject: [PATCH 21/52] Update documentation for non-contiguous subgroups

[ci skip]
---
 docs_sphinx/user/models.rst    | 16 ++++++----
 docs_sphinx/user/recording.rst | 54 +++++++++-------------------------
 2 files changed, 25 insertions(+), 45 deletions(-)

diff --git a/docs_sphinx/user/models.rst b/docs_sphinx/user/models.rst
index 3d6007e80..a9b8d00bf 100644
--- a/docs_sphinx/user/models.rst
+++ b/docs_sphinx/user/models.rst
@@ -182,11 +182,12 @@ neurons. In general ``G[i:j]`` refers to the neurons with indices from ``i``
 to ``j-1``, as in general in Python.
 
 For convenience, you can also use a single index, i.e. ``G[i]`` is equivalent
-to ``G[i:i+1]``. In some situations, it can be easier to provide a list of
-indices instead of a slice, Brian therefore also allows for this syntax. Note
-that this is restricted to cases that are strictly equivalent with slicing
-syntax, e.g. you can write ``G[[3, 4, 5]]`` instead of ``G[3:6]``, but you
-*cannot* write ``G[[3, 5, 7]]`` or ``G[[5, 4, 3]]``.
+to ``G[i:i+1]``. Brian also allows a simplified form of numpy's
+`integer array indexing <https://numpy.org/doc/stable/user/basics.indexing.html#integer-array-indexing>`_,
+to create "non-contiguous" subgroups (subgroups with "gaps" in them).
+You can for example refer to ``G[[0, 2, 4, 6, 8]]`` or ``G[[-5, -3, -1]]``.
+There are two restrictions to index arrays for subgroups: they cannot contain
+repeated indices, and the indices need to be provided in ascending order.
 
 Subgroups can be used in most places where regular groups are used, e.g. their
 state variables or spiking activity can be recorded using monitors, they can be
@@ -194,6 +195,11 @@ connected via `Synapses`, etc. In such situations, indices (e.g. the indices of
 the neurons to record from in a `StateMonitor`) are relative to the subgroup,
 not to the main group
 
+.. note::
+
+    Non-contiguous subgroups (i.e. subgroups created with the index array syntax)
+    cannot be used as the source/target groups for `Synapses`.
+
 .. admonition:: The following topics are not essential for beginners.
 
     |
diff --git a/docs_sphinx/user/recording.rst b/docs_sphinx/user/recording.rst
index 7190505f6..34b819d6d 100644
--- a/docs_sphinx/user/recording.rst
+++ b/docs_sphinx/user/recording.rst
@@ -214,8 +214,8 @@ monitor::
 
 Note that this technique cannot be applied in :ref:`standalone mode <cpp_standalone>`.
 
-Recording random subsets of neurons
------------------------------------
+Recording subsets of neurons
+----------------------------
 
 In large networks, you might only be interested in the activity of a
 random subset of neurons. While you can specify a ``record`` argument
@@ -227,41 +227,15 @@ by using a :ref:`subgroup <subgroups>`::
     group = NeuronGroup(1000, ...)
     spike_mon = SpikeMonitor(group[:100])  # only record first 100 neurons
 
-It might seem like a restriction that such a subgroup has to be contiguous, but
-the order of neurons in a group does not have any meaning as such; in a randomly
-ordered group of neurons, any contiguous group of neurons can be considered a
-random subset. If some aspects of your model *do* depend on the position of the
-neuron in a group (e.g. a ring model, where neurons are connected based on their
-distance in the ring, or a model where initial values or parameters span a
-range of values in a regular fashion), then this requires an extra step: instead
-of using the order of neurons in the group directly, or depending on the neuron
-index ``i``, create a new, shuffled, index variable as part of the model
-definition and then depend on this index instead::
-
-    group = NeuronGroup(10000, '''....
-                                  index : integer (constant)''')
-    indices = group.i[:]
-    np.random.shuffle(indices)
-    group.index = indices
-    # Then use 'index' in string expressions or use it as an index array
-    # for initial values/parameters defined as numpy arrays
-
-If this solution is not feasible for some reason, there is another approach that
-works for a `SpikeMonitor`/`EventMonitor`. You can add an additional flag to
-each neuron, stating whether it should be recorded or not. Then, you define a
-new :doc:`custom event </advanced/custom_events>` that is identical to the event you are
-interested in, but additionally requires the flag to be set. E.g. to only record
-the spikes of neurons with the ``to_record`` attribute set::
-
-    group = NeuronGroup(..., '''...
-                                to_record : boolean (constant)''',
-                        threshold='...', reset='...',
-                        events={'recorded_spike': '... and to_record'})
-    group.to_record = ...
-    mon_events = EventMonitor(group, 'recorded_spike')
-
-Note that this solution will evaluate the threshold condition for each neuron
-twice, and is therefore slightly less efficient. There's one additional caveat:
-you'll have to manually include ``and not_refractory`` in your ``events``
-definition if your neuron uses refractoriness. This is done automatically
-for the ``threshold`` condition, but not for any user-defined events.
+This also extends to non-contiguous subgroups, e.g. to record every second cell
+you can use::
+
+    group = NeuronGroup(1000, ...)
+    spike_mon = SpikeMonitor(group[::2])  # record every second neuron
+
+Note that this is less efficient than recording from a contiguous subgroup, since
+for every spike in ``group``, Brian needs to check whether its index is part of the
+subgroup (this check is quicker to do when the subgroup is a contiguous range).
+In many models, the order of neurons in a group does not have any meaning as such.
+If your model is randomly ordered, then recording from the first half of neurons
+is equivalent to recording from every second neuron, but more efficient.

From 7e495bbbbb2a8c1055d1b129fac55fe8dde66732 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 10 Nov 2022 15:25:00 +0100
Subject: [PATCH 22/52] Some more basic tests for incorrect subgroup creation

---
 brian2/tests/test_subgroup.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index e9088a85f..9692678cb 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -37,12 +37,20 @@ def test_creation():
     assert len(SGi2) == 3
     assert SGi2.contiguous
 
+    with pytest.raises(ValueError):
+        Subgroup(G)
+    with pytest.raises(ValueError):
+        Subgroup(G, start=3)
+    with pytest.raises(ValueError):
+        Subgroup(G, stop=3)
     with pytest.raises(IndexError):
         Subgroup(G, start=-1, stop=5)
     with pytest.raises(IndexError):
         Subgroup(G, start=1, stop=1)
     with pytest.raises(IndexError):
         Subgroup(G, start=1, stop=11)
+    with pytest.raises(ValueError):
+        Subgroup(G, start=1, stop=11, indices=[2, 4, 6])
     with pytest.raises(IndexError):
         Subgroup(G, indices=[])
     with pytest.raises(IndexError):

From aaacc26332d26c68f4ced6854896c531f1266ae8 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 14 Nov 2022 14:12:25 +0100
Subject: [PATCH 23/52] Fix exception type in subgroup tests
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

ValueError → TypeError
---
 brian2/tests/test_subgroup.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 9692678cb..f1c2db604 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -37,11 +37,11 @@ def test_creation():
     assert len(SGi2) == 3
     assert SGi2.contiguous
 
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         Subgroup(G)
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         Subgroup(G, start=3)
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         Subgroup(G, stop=3)
     with pytest.raises(IndexError):
         Subgroup(G, start=-1, stop=5)
@@ -49,7 +49,7 @@ def test_creation():
         Subgroup(G, start=1, stop=1)
     with pytest.raises(IndexError):
         Subgroup(G, start=1, stop=11)
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         Subgroup(G, start=1, stop=11, indices=[2, 4, 6])
     with pytest.raises(IndexError):
         Subgroup(G, indices=[])

From 22629b8d86676f75b29a3b80302c16ab10fa2a2e Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Mon, 14 Nov 2022 17:31:23 +0100
Subject: [PATCH 24/52] Test+doc for non-contiguous indices in SpatialNeuron

---
 brian2/tests/test_spatialneuron.py      |  5 +++++
 docs_sphinx/user/multicompartmental.rst | 24 +++++++++++++-----------
 2 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/brian2/tests/test_spatialneuron.py b/brian2/tests/test_spatialneuron.py
index 4e080cf09..07073dc82 100644
--- a/brian2/tests/test_spatialneuron.py
+++ b/brian2/tests/test_spatialneuron.py
@@ -661,9 +661,14 @@ def test_spatialneuron_indexing():
     assert len(neuron[sec.sec2.indices[:]]) == 16
     assert len(neuron[sec.sec2]) == 16
     assert len(neuron[:16:2].indices[:]) == 8
+    assert len(neuron['i < 16 and i % 2 == 0'].indices[:]) == 8
     assert len(neuron[:8][4:].indices[:]) == 4
     assert len(neuron[:8][::2].indices[:]) == 4
+    assert len(neuron[:8]['i % 2 == 0'].indices[:]) == 4
     assert len(neuron[::2][:8].indices[:]) == 8
+    assert len(neuron['i % 2 == 0'][:8].indices[:]) == 8
+    assert len(neuron['distance >= 100*um'].indices[:]) == 44
+    assert len(neuron['distance >= 100*um'][:8].indices[:]) == 8
     assert_equal(neuron.sec1.sec11.v, [3, 4, 5, 6]*volt)
     assert_equal(neuron.sec1.sec11[1].v, neuron.sec1.sec11.v[1])
     assert_equal(neuron.sec1.sec11[1:3].v, neuron.sec1.sec11.v[1:3])
diff --git a/docs_sphinx/user/multicompartmental.rst b/docs_sphinx/user/multicompartmental.rst
index 4c3b92d49..5fdcf2671 100644
--- a/docs_sphinx/user/multicompartmental.rst
+++ b/docs_sphinx/user/multicompartmental.rst
@@ -322,15 +322,22 @@ indices of compartments, or with the distance from the root::
     first_compartments = neuron[:3]
     first_compartments = neuron[0*um:30*um]
 
-However, note that this is restricted to contiguous indices which most of the
-time means that all compartments indexed in this way have to be part of the
-same section. Such indices can be acquired directly from the morphology::
+Subgroups can also consist of several compartments that are not directly connected to
+each other and lie in different sections. This can be achieved by using a list of
+indices or a string expression instead of a slice. For example, to create a subgroup
+referring to the compartments with indices 1, 2, 5, and 6, you can use::
 
-    axon = neuron[morpho.axon.indices[:]]
+    subset = neuron[[1, 2, 5, 6]]
 
-or, more concisely::
+The same restrictions as for general :ref:`subgroups` apply: the given indices need to
+be in ascending order without duplicates. Note that you can get indices from the
+`Morphology`, e.g. by asking for ``morpho.axon.indices[:]``.
 
-    axon = neuron[morpho.axon]
+With string indexing, the compartments can be selected by referring to their attributes.
+For example, to select all compartments that are at a distance of more than 100µm from
+the soma, you can use::
+
+    distal = neuron['distance > 100*um']
 
 Synaptic inputs
 ~~~~~~~~~~~~~~~
@@ -409,8 +416,3 @@ Again the location of the threshold can be specified with spatial position::
                            threshold_location=morpho.axon[30*um],
                            refractory='m > 0.4')
 
-Subgroups
-~~~~~~~~~
-
-In the same way that you can refer to a subset of neurons in a `NeuronGroup`,
-you can also refer to a subset of compartments in a `SpatialNeuron`

From 65ec4c3419ea2c430684d9f2242b1201fa1cf18b Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Tue, 22 Nov 2022 11:59:20 +0100
Subject: [PATCH 25/52] Remove deprecated usage of np.bool

---
 brian2/core/variables.py                | 2 +-
 brian2/devices/cpp_standalone/device.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/brian2/core/variables.py b/brian2/core/variables.py
index 8ce33be56..5bdd9d914 100644
--- a/brian2/core/variables.py
+++ b/brian2/core/variables.py
@@ -188,7 +188,7 @@ def __init__(
 
     @property
     def is_boolean(self):
-        return np.issubdtype(self.dtype, np.bool_)
+        return np.issubdtype(self.dtype, bool)
 
     @property
     def is_integer(self):
diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index 0fe765f4f..5986ed610 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/brian2/devices/cpp_standalone/device.py
@@ -588,7 +588,7 @@ 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)
+            variables.add_auxiliary_variable("_cond", dtype=bool)
 
             abstract_code = "_cond = " + item
             namespace = get_local_namespace(level=level + 2)

From 82e3ffc105dfc9f35d912c6b2e1ef2376a14384b Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 23 Mar 2023 15:52:08 +0100
Subject: [PATCH 26/52] Fragile but basic working version of SynapticSubgroup

Needs refactoring and documentation
---
 brian2/synapses/synapses.py | 61 ++++++++++++++++++++++++++-----------
 1 file changed, 44 insertions(+), 17 deletions(-)

diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 6aeb5229a..e3ba88b06 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -529,9 +529,9 @@ def find_synapses(index, synaptic_neuron):
     return synapses
 
 
-class SynapticSubgroup(object):
+class SynapticSubgroup(Group):
     """
-    A simple subgroup of `Synapses` that can be used for indexing.
+    A subgroup of `Synapses` that can be used for indexing and for accessing variables.
 
     Parameters
     ----------
@@ -542,21 +542,39 @@ class SynapticSubgroup(object):
         when new synapses where added after creating this object.
     """
 
-    def __init__(self, synapses, indices):
+    def __init__(self, synapses, indices, name=None):
         self.synapses = weakproxy_with_fallback(synapses)
+        self.source = weakproxy_with_fallback(synapses.source)
+        self.target = weakproxy_with_fallback(synapses.target)
+        self.multisynaptic_index = self.synapses.multisynaptic_index
         self._stored_indices = indices
+        self._N = len(indices)
         self._synaptic_pre = synapses.variables["_synaptic_pre"]
         self._source_N = self._synaptic_pre.size  # total number of synapses
 
-    def _indices(self, index_var="_idx"):
-        if index_var != "_idx":
-            raise AssertionError(f"Did not expect index {index_var} here.")
-        if len(self._synaptic_pre.get_value()) != self._source_N:
-            raise RuntimeError(
-                "Synapses have been added/removed since this "
-                "synaptic subgroup has been created"
-            )
-        return self._stored_indices
+        if name is None:
+            name = f"{self.synapses.name}_subgroup*"
+
+        Group.__init__(
+            self,
+            name=name,
+        )
+        self.variables = Variables(self, default_index="_sub_idx")
+        self.variables.add_references(synapses, list(synapses.variables.keys()))
+        self.variables.add_array(
+            "_sub_idx",
+            size=self._N,
+            dtype=np.int32,
+            values=indices,
+            index="_idx",
+            constant=True,
+            read_only=True,
+            unique=True,
+        )
+
+        self._indices = SynapticIndexing(self, "_sub_idx")
+
+        self._enable_group_attributes()
 
     def __len__(self):
         return len(self._stored_indices)
@@ -569,8 +587,9 @@ def __repr__(self):
 
 
 class SynapticIndexing(object):
-    def __init__(self, synapses):
-        self.synapses = weakref.proxy(synapses)
+    def __init__(self, synapses, default_idx="_idx"):
+        self.synapses = synapses
+        self.default_idx = default_idx
         self.source = weakproxy_with_fallback(self.synapses.source)
         self.target = weakproxy_with_fallback(self.synapses.target)
         self.synaptic_pre = synapses.variables["_synaptic_pre"]
@@ -596,8 +615,13 @@ def __call__(self, index=None, index_var="_idx"):
             if hasattr(index, "_indices"):
                 final_indices = index._indices(index_var=index_var).astype(np.int32)
             elif isinstance(index, slice):
-                start, stop, step = index.indices(len(self.synaptic_pre.get_value()))
-                final_indices = np.arange(start, stop, step, dtype=np.int32)
+                if self.default_idx == "_idx":
+                    start, stop, step = index.indices(
+                        len(self.synaptic_pre.get_value())
+                    )
+                    final_indices = np.arange(start, stop, step, dtype=np.int32)
+                else:
+                    final_indices = self.synapses._stored_indices
             else:
                 final_indices = np.asarray(index)
         elif isinstance(index, tuple):
@@ -658,7 +682,10 @@ def __call__(self, index=None, index_var="_idx"):
         else:
             raise IndexError(f"Unsupported index type {type(index)}")
 
-        if index_var not in ("_idx", "0"):
+        if (
+            index_var not in ("_idx", "0")
+            and not getattr(index_var, "name", None) == "_sub_idx"
+        ):
             return index_var.get_value()[final_indices.astype(np.int32)]
         else:
             return final_indices.astype(np.int32)

From 55094371735ea1b5760fcecb23547fcce2ce866e Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 24 Mar 2023 14:04:15 +0100
Subject: [PATCH 27/52] Remove unnecessary imports

---
 brian2/groups/subgroup.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/brian2/groups/subgroup.py b/brian2/groups/subgroup.py
index 65512be21..9a052afa3 100644
--- a/brian2/groups/subgroup.py
+++ b/brian2/groups/subgroup.py
@@ -1,6 +1,3 @@
-import numbers
-from collections.abc import Sequence
-
 import numpy as np
 
 from brian2.core.base import weakproxy_with_fallback

From fe861b7b9dd13d464b720fcb57232679436577ea Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 24 Mar 2023 17:43:33 +0100
Subject: [PATCH 28/52] Add tests for synaptic subgroups

---
 brian2/synapses/synapses.py   |  1 +
 brian2/tests/test_synapses.py | 56 +++++++++++++++++++++++++++++++++++
 2 files changed, 57 insertions(+)

diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 328c0bad6..1eba3ea71 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -543,6 +543,7 @@ class SynapticSubgroup(Group):
     """
 
     def __init__(self, synapses, indices, name=None):
+        indices = np.atleast_1d(indices)  # Deal with scalar indices
         self.synapses = weakproxy_with_fallback(synapses)
         self.source = weakproxy_with_fallback(synapses.source)
         self.target = weakproxy_with_fallback(synapses.target)
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index 3720705b1..5df61c78b 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -741,6 +741,62 @@ def test_state_variable_indexing():
         S.w.__getitem__(1.5)
 
 
+def test_state_variable_indexing_with_subgroups():
+    G1 = NeuronGroup(5, "v:volt")
+    G1.v = "i*mV"
+    G2 = NeuronGroup(7, "v:volt")
+    G2.v = "10*mV + i*mV"
+    S = Synapses(G1, G2, "w:1", multisynaptic_index="k")
+    S.connect(True, n=2)
+    S.w[:, :, 0] = "5*i + j"
+    S.w[:, :, 1] = "35 + 5*i + j"
+
+    # Slicing
+    assert len(S[:].w) == len(S[:, :].w) == len(S[:, :, :].w) == len(G1) * len(G2) * 2
+    assert len(S[0:, 0:].w) == len(S[0:, 0:, 0:].w) == len(G1) * len(G2) * 2
+    assert len(S[0::2, 0:].w) == 3 * len(G2) * 2
+    assert len(S[0, :].w) == len(S[0, :, :].w) == len(G2) * 2
+    assert len(S[0:2, :].w) == len(S[0:2, :, :].w) == 2 * len(G2) * 2
+    assert len(S[:2, :].w) == len(S[:2, :, :].w) == 2 * len(G2) * 2
+    assert len(S[0:4:2, :].w) == len(S[0:4:2, :, :].w) == 2 * len(G2) * 2
+    assert len(S[:4:2, :].w) == len(S[:4:2, :, :].w) == 2 * len(G2) * 2
+    assert len(S[:, 0].w) == len(S[:, 0, :].w) == len(G1) * 2
+    assert len(S[:, 0:2].w) == len(S[:, 0:2, :].w) == 2 * len(G1) * 2
+    assert len(S[:, :2].w) == len(S[:, :2, :].w) == 2 * len(G1) * 2
+    assert len(S[:, 0:4:2].w) == len(S[:, 0:4:2, :].w) == 2 * len(G1) * 2
+    assert len(S[:, :4:2].w) == len(S[:, :4:2, :].w) == 2 * len(G1) * 2
+    assert len(S[:, :, 0].w) == len(G1) * len(G2)
+    assert len(S[:, :, 0:2].w) == len(G1) * len(G2) * 2
+    assert len(S[:, :, :2].w) == len(G1) * len(G2) * 2
+    assert len(S[:, :, 0:2:2].w) == len(G1) * len(G2)
+    assert len(S[:, :, :2:2].w) == len(G1) * len(G2)
+
+    # 1d indexing is directly indexing synapses!
+    assert len(S[:].w) == len(S[0:].w)
+    assert len(S[[0, 1]].w) == len(S[3:5].w) == 2
+    assert len(S[:].w) == len(S[np.arange(len(G1) * len(G2) * 2)].w)
+    assert S[3].w == S[np.int32(3)].w == S[np.int64(3)].w  # See issue #888
+
+    # Array-indexing (not yet supported for synapse index)
+    assert_equal(S[:, 0:3].w[:], S[:, [0, 1, 2]].w[:])
+    assert_equal(S[:, 0:3].w[:], S[np.arange(len(G1)), [0, 1, 2]].w[:])
+
+    # string-based indexing
+    assert_equal(S[0:3, :].w[:], S["i<3"].w[:])
+    assert_equal(S[:, 0:3].w[:], S["j<3"].w[:])
+    assert_equal(S[:, :, 0].w[:], S["k == 0"].w[:])
+    assert_equal(S[0:3, :].w[:], S["v_pre < 2.5*mV"].w[:])
+    assert_equal(S[:, 0:3].w[:], S["v_post < 12.5*mV"].w[:])
+
+    # invalid indices
+    with pytest.raises(IndexError):
+        S.__getitem__((1, 2, 3, 4))
+    with pytest.raises(IndexError):
+        S.__getitem__(object())
+    with pytest.raises(IndexError):
+        S.__getitem__(1.5)
+
+
 def test_indices():
     G = NeuronGroup(10, "v : 1")
     S = Synapses(G, G, "")

From 73d214f0930f8eb967466dd8afb996ba99270e03 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 24 Mar 2023 18:36:25 +0100
Subject: [PATCH 29/52] Support array indexing for multi-synaptic indices

---
 brian2/synapses/synapses.py   | 20 +++++++++-----------
 brian2/tests/test_synapses.py |  6 +++++-
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 1eba3ea71..08c813946 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -474,7 +474,7 @@ def slice_to_test(x):
         pass
 
     if isinstance(x, slice):
-        if isinstance(x, slice) and x == slice(None):
+        if x == slice(None):
             # No need for testing
             return lambda y: np.repeat(True, len(y))
         start, stop, step = x.start, x.stop, x.step
@@ -668,20 +668,18 @@ def __call__(self, index=None, index_var="_idx"):
                         "'multisynaptic_index' when you create "
                         "the Synapses object."
                     )
+
+                # We want to access the raw arrays here, not go through the Variable
+                synapse_numbers = self.synapse_number.get_value()[matching_synapses]
+
                 if isinstance(k_indices, (numbers.Integral, slice)):
                     test_k = slice_to_test(k_indices)
+                    final_indices = matching_synapses[test_k(synapse_numbers)]
                 else:
-                    raise NotImplementedError(
-                        "Indexing synapses with arrays notimplemented yet"
-                    )
+                    final_indices = matching_synapses[
+                        np.in1d(synapse_numbers, k_indices)
+                    ]
 
-                # We want to access the raw arrays here, not go through the Variable
-                synapse_numbers = self.synapse_number.get_value()[matching_synapses]
-                final_indices = np.intersect1d(
-                    matching_synapses,
-                    np.flatnonzero(test_k(synapse_numbers)),
-                    assume_unique=True,
-                )
         else:
             raise IndexError(f"Unsupported index type {type(index)}")
 
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index 5df61c78b..8ed1bde86 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -721,9 +721,11 @@ def test_state_variable_indexing():
     assert len(S.w[:]) == len(S.w[np.arange(len(G1) * len(G2) * 2)])
     assert S.w[3] == S.w[np.int32(3)] == S.w[np.int64(3)]  # See issue #888
 
-    # Array-indexing (not yet supported for synapse index)
+    # Array-indexing
     assert_equal(S.w[:, 0:3], S.w[:, [0, 1, 2]])
     assert_equal(S.w[:, 0:3], S.w[np.arange(len(G1)), [0, 1, 2]])
+    assert_equal(S.w[:, 0:3], S.w[:, [0, 1, 2], [0, 1]])
+    assert_equal(S.w[:, 0:3, 0], S.w[:, [0, 1, 2], [0]])
 
     # string-based indexing
     assert_equal(S.w[0:3, :], S.w["i<3"])
@@ -780,6 +782,8 @@ def test_state_variable_indexing_with_subgroups():
     # Array-indexing (not yet supported for synapse index)
     assert_equal(S[:, 0:3].w[:], S[:, [0, 1, 2]].w[:])
     assert_equal(S[:, 0:3].w[:], S[np.arange(len(G1)), [0, 1, 2]].w[:])
+    assert_equal(S[:, 0:3].w[:], S[:, [0, 1, 2], [0, 1]].w[:])
+    assert_equal(S[:, 0:3, 0].w[:], S[:, [0, 1, 2], [0]].w[:])
 
     # string-based indexing
     assert_equal(S[0:3, :].w[:], S["i<3"].w[:])

From e51b0ae468306dce4eb7eb49840854efc2736b41 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 29 Mar 2023 18:02:19 +0200
Subject: [PATCH 30/52] Refactor `Group.Indexing` for better readability

---
 brian2/groups/group.py | 117 ++++++++++++++++++++++++-----------------
 1 file changed, 68 insertions(+), 49 deletions(-)

diff --git a/brian2/groups/group.py b/brian2/groups/group.py
index 4b766263a..639b51a8a 100644
--- a/brian2/groups/group.py
+++ b/brian2/groups/group.py
@@ -218,7 +218,8 @@ class Indexing:
     specific indices. Stores strong references to the necessary variables so
     that basic indexing (i.e. slicing, integer arrays/values, ...) works even
     when the respective `VariableOwner` no longer exists. Note that this object
-    does not handle string indexing.
+    does not handle string indexing (handled by `IndexWrapper` and
+    `VariableView.get_with_expression`).
     """
 
     def __init__(self, group, default_index="_idx"):
@@ -256,59 +257,77 @@ def __call__(self, item=slice(None), index_var=None):  # noqa: B008
             raise IndexError(
                 f"Can only interpret 1-d indices, got {len(item)} dimensions."
             )
-        else:
-            if isinstance(item, str) and item == "True":
-                item = slice(None)
-            if isinstance(item, slice):
-                if index_var == "0":
-                    return 0
-                if index_var == "_idx":
-                    start, stop, step = item.indices(int(self.N.get_value()))
-                else:
-                    start, stop, step = item.indices(index_var.size)
-                index_array = np.arange(start, stop, step)
-            else:
-                index_array = np.asarray(item)
-                if index_array.dtype == bool:
-                    index_array = np.nonzero(index_array)[0]
-                elif not np.issubdtype(index_array.dtype, np.signedinteger):
-                    raise TypeError(
-                        "Indexing is only supported for integer "
-                        "and boolean arrays, not for type "
-                        f"{index_array.dtype}"
-                    )
 
-            if index_array.size == 0:
-                return index_array
-            if index_var != "_idx":
-                try:
-                    index_array = index_var.get_value()[index_array]
-                except IndexError as ex:
-                    # We try to emulate numpy's indexing semantics here:
-                    # slices never lead to IndexErrors, instead they return an
-                    # empty array if they don't match anything
-                    if isinstance(item, slice):
-                        return np.array([], dtype=np.int32)
-                    else:
-                        raise ex
-            else:
-                N = int(self.N.get_value())
-                if np.min(index_array) < -N:
-                    raise IndexError(
-                        "Illegal index {} for a group of size {}".format(
-                            np.min(index_array), N
-                        )
+        index_array = self._to_index_array(item, index_var)
+
+        if index_array.size == 0:
+            return index_array
+
+        if index_var not in ("_idx", "0"):
+            try:
+                index_array = index_var.get_value()[index_array]
+            except IndexError as ex:
+                # We try to emulate numpy's indexing semantics here:
+                # slices never lead to IndexErrors, instead they return an
+                # empty array if they don't match anything
+                if isinstance(item, slice):
+                    return np.array([], dtype=np.int32)
+                else:
+                    raise ex
+        else:
+            N = int(self.N.get_value())
+            if np.min(index_array) < -N:
+                raise IndexError(
+                    "Illegal index {} for a group of size {}".format(
+                        np.min(index_array), N
                     )
-                if np.max(index_array) >= N:
-                    raise IndexError(
-                        "Illegal index {} for a group of size {}".format(
-                            np.max(index_array), N
-                        )
+                )
+            if np.max(index_array) >= N:
+                raise IndexError(
+                    "Illegal index {} for a group of size {}".format(
+                        np.max(index_array), N
                     )
+                )
 
-                index_array = index_array % N  # interpret negative indices
+            index_array = index_array % N  # interpret negative indices
 
-            return index_array
+        return index_array
+
+    def _to_index_array(self, item, index_var):
+        """
+        Convert slices, integer/boolean arrays to an integer array of indices.
+
+        Parameters
+        ----------
+        item: slice, array, int
+            The indices to translate.
+        index_var : `ArrayVariable`, str
+            The index variable.
+        Returns
+        -------
+        indices : `numpy.ndarray`
+            The flat indices corresponding to the indices given in `item`.
+        """
+        if isinstance(item, slice):
+            if index_var == "0":
+                index_array = np.array(0)
+            else:
+                index_size = (
+                    int(self.N.get_value()) if index_var == "_idx" else index_var.size
+                )
+                start, stop, step = item.indices(index_size)
+                index_array = np.arange(start, stop, step)
+        else:  # array, sequence, or single value
+            index_array = np.asarray(item)
+            if index_array.dtype == bool:
+                index_array = np.flatnonzero(index_array)
+            elif not np.issubdtype(index_array.dtype, np.signedinteger):
+                raise TypeError(
+                    "Indexing is only supported for integer "
+                    "and boolean arrays, not for type "
+                    f"{index_array.dtype}"
+                )
+        return index_array
 
 
 class IndexWrapper:

From 6ceba162243b77489364a6fc1778848816a8f222 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 29 Mar 2023 18:03:07 +0200
Subject: [PATCH 31/52] Add new tests for synaptic subgroups

---
 brian2/tests/test_synapses.py | 34 ++++++++++++++++++++++++++++++++--
 1 file changed, 32 insertions(+), 2 deletions(-)

diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index 8ed1bde86..f6f8779e5 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -3563,17 +3563,47 @@ def test_synaptic_subgroups():
     from_3 = syn[3, :]
     assert len(from_3) == 3
     assert all(syn.i[from_3] == 3)
+    assert all(from_3.i == 3)
     assert_array_equal(syn.j[from_3], np.arange(3))
+    assert_array_equal(from_3.j, np.arange(3))
 
     to_2 = syn[:, 2]
     assert len(to_2) == 5
     assert all(syn.j[to_2] == 2)
-    assert_array_equal(syn.i[to_2], np.arange(5))
+    assert all(to_2.j == 2)
 
     mixed = syn[1:3, :2]
     assert len(mixed) == 4
     connections = {(i, j) for i, j in zip(syn.i[mixed], syn.j[mixed])}
-    assert connections == {(1, 0), (1, 1), (2, 0), (2, 1)}
+    expected = {(1, 0), (1, 1), (2, 0), (2, 1)}
+    # assert connections == expected
+    connections = {(i, j) for i, j in zip(mixed.i[:], mixed.j[:])}
+    assert connections == expected
+
+    string_based = syn["i > 1 and j % 2 == 0"]
+    assert len(string_based) == 6
+    connections = {(i, j) for i, j in zip(syn.i[string_based], syn.j[string_based])}
+    expected = {(2, 0), (2, 2), (3, 0), (3, 2), (4, 0), (4, 2)}
+    assert connections == expected
+    connections = {(i, j) for i, j in zip(string_based.i[:], string_based.j[:])}
+    assert connections == expected
+
+    array_based_2d = syn[[2, 3, 4], [0, 2]]  # same as above
+    assert len(array_based_2d) == 6
+    connections = {(i, j) for i, j in zip(syn.i[array_based_2d], syn.j[array_based_2d])}
+    expected = {(2, 0), (2, 2), (3, 0), (3, 2), (4, 0), (4, 2)}
+    assert connections == expected
+    connections = {(i, j) for i, j in zip(array_based_2d.i[:], array_based_2d.j[:])}
+    assert connections == expected
+
+    indices = np.flatnonzero((syn.i > 1) & (syn.j % 2 == 0))
+    array_based_1d = syn[indices]  # same as above
+    assert len(array_based_1d) == 6
+    connections = {(i, j) for i, j in zip(syn.i[array_based_1d], syn.j[array_based_1d])}
+    expected = {(2, 0), (2, 2), (3, 0), (3, 2), (4, 0), (4, 2)}
+    assert connections == expected
+    connections = {(i, j) for i, j in zip(array_based_1d.i[:], array_based_1d.j[:])}
+    assert connections == expected
 
 
 @pytest.mark.codegen_independent

From dbbb74163c4fa2044bc7a27832e453cd0217701e Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 30 Mar 2023 18:23:52 +0200
Subject: [PATCH 32/52] Small refactoring for synaptic indexing

Also fixing a few issues related for standalone
---
 brian2/devices/cpp_standalone/device.py |   9 +-
 brian2/synapses/synapses.py             | 156 +++++++++++-------------
 brian2/tests/test_synapses.py           |   8 +-
 3 files changed, 81 insertions(+), 92 deletions(-)

diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index 25c54f13a..ad554c27f 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/brian2/devices/cpp_standalone/device.py
@@ -462,8 +462,8 @@ def resize(self, var, new_size):
         self.main_queue.append(("resize_array", (array_name, new_size)))
 
     def variableview_set_with_index_array(self, variableview, item, value, check_units):
-        if isinstance(item, slice) and item == slice(None):
-            item = "True"
+        if item == "True":
+            item = slice(None)
         value = Quantity(value)
 
         if (
@@ -478,7 +478,7 @@ def variableview_set_with_index_array(self, variableview, item, value, check_uni
                 ("set_by_single_value", (array_name, item, value_str))
             )
         # Simple case where we don't have to do any indexing
-        elif item == "True" and variableview.index_var in ("_idx", "0"):
+        elif item == slice(None) and variableview.index_var in ("_idx", "0"):
             self.fill_with_array(variableview.variable, value)
         else:
             # We have to calculate indices. This will not work for synaptic
@@ -717,6 +717,9 @@ def code_object(
                         cache.get(var, np.empty(0, dtype=int)), value
                     )
                     do_not_invalidate.add(var)
+                # Also update the size attribute
+                variables["_synaptic_pre"].size += variables["sources"].size
+                variables["_synaptic_post"].size += variables["targets"].size
 
         codeobj = super().code_object(
             owner,
diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 08c813946..c1fd29a47 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -25,7 +25,7 @@
     Equations,
     check_subexpressions,
 )
-from brian2.groups.group import CodeRunner, Group, get_dtype
+from brian2.groups.group import CodeRunner, Group, Indexing, get_dtype
 from brian2.groups.neurongroup import (
     SubexpressionUpdater,
     check_identifier_pre_post,
@@ -561,6 +561,7 @@ def __init__(self, synapses, indices, name=None):
             name=name,
         )
         self.variables = Variables(self, default_index="_sub_idx")
+        self.variables.add_constant("N", self._N)
         self.variables.add_references(synapses, list(synapses.variables.keys()))
         self.variables.add_array(
             "_sub_idx",
@@ -573,7 +574,7 @@ def __init__(self, synapses, indices, name=None):
             unique=True,
         )
 
-        self._indices = SynapticIndexing(self, "_sub_idx")
+        self._indices = SynapticIndexing(self, self.variables["_sub_idx"])
 
         self._enable_group_attributes()
 
@@ -587,12 +588,11 @@ def __repr__(self):
         )
 
 
-class SynapticIndexing:
+class SynapticIndexing(Indexing):
     def __init__(self, synapses, default_idx="_idx"):
-        self.synapses = synapses
-        self.default_idx = default_idx
-        self.source = weakproxy_with_fallback(self.synapses.source)
-        self.target = weakproxy_with_fallback(self.synapses.target)
+        super().__init__(synapses, default_idx)
+        self.source = weakproxy_with_fallback(synapses.source)
+        self.target = weakproxy_with_fallback(synapses.target)
         self.synaptic_pre = synapses.variables["_synaptic_pre"]
         self.synaptic_post = synapses.variables["_synaptic_post"]
         if synapses.multisynaptic_index is not None:
@@ -600,96 +600,82 @@ def __init__(self, synapses, default_idx="_idx"):
         else:
             self.synapse_number = None
 
-    def __call__(self, index=None, index_var="_idx"):
+    def __call__(self, item=slice(None), index_var=None):  # noqa: B008
         """
         Returns synaptic indices for `index`, which can be a tuple of indices
         (including arrays and slices), a single index or a string.
 
         """
-        if index is None or (isinstance(index, str) and index == "True"):
-            index = slice(None)
-
-        if not isinstance(index, (tuple, str)) and (
-            isinstance(index, (numbers.Integral, np.ndarray, slice, Sequence))
-            or hasattr(index, "_indices")
-        ):
-            if hasattr(index, "_indices"):
-                final_indices = index._indices(index_var=index_var).astype(np.int32)
-            elif isinstance(index, slice):
-                if self.default_idx == "_idx":
-                    start, stop, step = index.indices(
-                        len(self.synaptic_pre.get_value())
-                    )
-                    final_indices = np.arange(start, stop, step, dtype=np.int32)
-                else:
-                    final_indices = self.synapses._stored_indices
-            else:
-                final_indices = np.asarray(index)
-        elif isinstance(index, tuple):
-            if len(index) == 2:  # two indices (pre- and postsynaptic cell)
-                index = (index[0], index[1], slice(None))
-            elif len(index) > 3:
-                raise IndexError(f"Need 1, 2 or 3 indices, got {len(index)}.")
-
-            i_indices, j_indices, k_indices = index
-            # Convert to absolute indices (e.g. for subgroups)
-            # Allow the indexing to fail, we'll later return an empty array in
-            # that case
-            try:
-                if hasattr(
-                    i_indices, "_indices"
-                ):  # will return absolute indices already
-                    i_indices = i_indices._indices()
-                else:
-                    i_indices = self.source._indices(i_indices)
-                pre_synapses = find_synapses(i_indices, self.synaptic_pre.get_value())
-            except IndexError:
-                pre_synapses = np.array([], dtype=np.int32)
-            try:
-                if hasattr(j_indices, "_indices"):
-                    j_indices = j_indices._indices()
-                else:
-                    j_indices = self.target._indices(j_indices)
-                post_synapses = find_synapses(j_indices, self.synaptic_post.get_value())
-            except IndexError:
-                post_synapses = np.array([], dtype=np.int32)
+        if index_var is None:
+            index_var = self.default_index
+
+        if not isinstance(item, tuple):
+            # 1d indexing = synaptic indices
+            if hasattr(item, "_indices"):
+                item = item._indices()
+            final_indices = self._to_index_array(item, index_var)
+        else:
+            # 2d or 3d indexing = pre-/post-synaptic indices and (optionally) synapse number
+            if len(item) == 2:  # two indices (pre- and postsynaptic cell)
+                item = (item[0], item[1], slice(None))
+            elif len(item) > 3:
+                raise IndexError(f"Need 1, 2 or 3 indices, got {len(item)}.")
+
+            i_indices, j_indices, k_indices = item
+
+            if k_indices != slice(None) and self.synapse_number is None:
+                raise IndexError(
+                    "To index by the third dimension you need "
+                    "to switch on the calculation of the "
+                    "'multisynaptic_index' when you create "
+                    "the Synapses object."
+                )
 
-            matching_synapses = np.intersect1d(
-                pre_synapses, post_synapses, assume_unique=True
+            final_indices = self._from_3d_indices_to_index_array(
+                i_indices, j_indices, k_indices
             )
 
-            if isinstance(k_indices, slice) and k_indices == slice(None):
-                final_indices = matching_synapses
-            else:
-                if self.synapse_number is None:
-                    raise IndexError(
-                        "To index by the third dimension you need "
-                        "to switch on the calculation of the "
-                        "'multisynaptic_index' when you create "
-                        "the Synapses object."
-                    )
-
-                # We want to access the raw arrays here, not go through the Variable
-                synapse_numbers = self.synapse_number.get_value()[matching_synapses]
+        if index_var not in ("_idx", "0"):
+            return index_var.get_value()[final_indices.astype(np.int32)]
+        else:
+            return final_indices.astype(np.int32)
 
-                if isinstance(k_indices, (numbers.Integral, slice)):
-                    test_k = slice_to_test(k_indices)
-                    final_indices = matching_synapses[test_k(synapse_numbers)]
-                else:
-                    final_indices = matching_synapses[
-                        np.in1d(synapse_numbers, k_indices)
-                    ]
+    def _from_3d_indices_to_index_array(self, i_indices, j_indices, k_indices):
+        # Convert to absolute indices (e.g. for subgroups)
+        # Allow the indexing to fail, we'll later return an empty array in
+        # that case
+        try:
+            if hasattr(i_indices, "_indices"):  # will return absolute indices already
+                i_indices = i_indices._indices()
+            else:
+                i_indices = self.source._indices(i_indices)
+            pre_synapses = find_synapses(i_indices, self.synaptic_pre.get_value())
+        except IndexError:
+            pre_synapses = np.array([], dtype=np.int32)
+        try:
+            if hasattr(j_indices, "_indices"):
+                j_indices = j_indices._indices()
+            else:
+                j_indices = self.target._indices(j_indices)
+            post_synapses = find_synapses(j_indices, self.synaptic_post.get_value())
+        except IndexError:
+            post_synapses = np.array([], dtype=np.int32)
+        matching_synapses = np.intersect1d(
+            pre_synapses, post_synapses, assume_unique=True
+        )
 
+        if k_indices == slice(None):
+            final_indices = matching_synapses
         else:
-            raise IndexError(f"Unsupported index type {type(index)}")
+            # We want to access the raw arrays here, not go through the Variable
+            synapse_numbers = self.synapse_number.get_value()[matching_synapses]
 
-        if (
-            index_var not in ("_idx", "0")
-            and not getattr(index_var, "name", None) == "_sub_idx"
-        ):
-            return index_var.get_value()[final_indices.astype(np.int32)]
-        else:
-            return final_indices.astype(np.int32)
+            if isinstance(k_indices, (numbers.Integral, slice)):
+                test_k = slice_to_test(k_indices)
+                final_indices = matching_synapses[test_k(synapse_numbers)]
+            else:
+                final_indices = matching_synapses[np.in1d(synapse_numbers, k_indices)]
+        return final_indices
 
 
 class Synapses(Group):
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index f6f8779e5..3c7f6475b 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -737,9 +737,9 @@ def test_state_variable_indexing():
     # invalid indices
     with pytest.raises(IndexError):
         S.w.__getitem__((1, 2, 3, 4))
-    with pytest.raises(IndexError):
+    with pytest.raises(TypeError):
         S.w.__getitem__(object())
-    with pytest.raises(IndexError):
+    with pytest.raises(TypeError):
         S.w.__getitem__(1.5)
 
 
@@ -795,9 +795,9 @@ def test_state_variable_indexing_with_subgroups():
     # invalid indices
     with pytest.raises(IndexError):
         S.__getitem__((1, 2, 3, 4))
-    with pytest.raises(IndexError):
+    with pytest.raises(TypeError):
         S.__getitem__(object())
-    with pytest.raises(IndexError):
+    with pytest.raises(TypeError):
         S.__getitem__(1.5)
 
 

From 2ff3cb116043a567720877718660f612d7013755 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 31 Mar 2023 20:40:21 +0200
Subject: [PATCH 33/52] Fix a test

---
 brian2/tests/test_synapses.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index 3c7f6475b..ce33cbed5 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -1250,8 +1250,8 @@ def test_transmission_scalar_delay_different_clocks():
             run(2 * ms)
             assert len(l) == 1, "expected a warning, got %d" % len(l)
             assert l[0][1].endswith("synapses_dt_mismatch")
-
-    run(0 * ms)
+    else:
+        run(2 * ms)
     assert_allclose(mon[0].v[mon.t < 0.5 * ms], 0)
     assert_allclose(mon[0].v[mon.t >= 0.5 * ms], 1)
     assert_allclose(mon[1].v[mon.t < 1.5 * ms], 0)

From fe0add72821ba584999add6f9b82d9b64b79876c Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 31 Mar 2023 20:41:57 +0200
Subject: [PATCH 34/52] Write sizes of dynamic arrays to disk after run in
 standalone mode

---
 brian2/devices/cpp_standalone/device.py       | 50 +++++++++++--------
 .../cpp_standalone/templates/objects.cpp      | 23 +++++++--
 2 files changed, 47 insertions(+), 26 deletions(-)

diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index ad554c27f..ddb5dc6da 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/brian2/devices/cpp_standalone/device.py
@@ -528,7 +528,14 @@ def get_value(self, var, access_data=True):
         # however (values that have been set with explicit values and not
         # changed in code objects)
         if self.array_cache.get(var, None) is not None:
-            return self.array_cache[var]
+            values = self.array_cache[var]
+            if isinstance(var, DynamicArrayVariable):
+                # Make sure that size information is up-to-date as well
+                if var.ndim == 2:
+                    var.size = values.shape
+                else:
+                    var.size = values.size
+            return values
         else:
             # After the network has been run, we can retrieve the values from
             # disk
@@ -537,27 +544,16 @@ def get_value(self, var, access_data=True):
                 fname = os.path.join(self.project_dir, self.get_array_filename(var))
                 with open(fname, "rb") as f:
                     data = np.fromfile(f, dtype=dtype)
-                # This is a bit of an heuristic, but our 2d dynamic arrays are
-                # only expanding in one dimension, we assume here that the
-                # other dimension has size 0 at the beginning
-                if isinstance(var.size, tuple) and len(var.size) == 2:
-                    if var.size[0] * var.size[1] == len(data):
-                        size = var.size
-                    elif var.size[0] == 0:
-                        size = (len(data) // var.size[1], var.size[1])
-                    elif var.size[1] == 0:
-                        size = (var.size[0], len(data) // var.size[0])
-                    else:
-                        raise IndexError(
-                            "Do not now how to deal with 2d "
-                            f"array of size {var.size!s}, the array on "
-                            f"disk has length {len(data)}."
-                        )
+                if (
+                    isinstance(var.size, tuple) and len(var.size) == 2
+                ):  # 2D dynamic array
+                    values = data.reshape(var.size)
+                else:
+                    values = data
+                # assert (np.atleast_1d(values.shape) == np.atleast_1d(var.size)).all(), f"{values.shape} ≠ {var.size} ({var.name})"
+                self.array_cache[var] = values
+                return values
 
-                    var.size = size
-                    return data.reshape(var.size)
-                var.size = len(data)
-                return data
             raise NotImplementedError(
                 "Cannot retrieve the values of state "
                 "variables in standalone code before the "
@@ -1301,6 +1297,18 @@ def run(self, directory, with_output, run_args):
             except ReferenceError:
                 pass
 
+        # Make sure the size information is up-to-date for dynamic variables
+        # Note that the actual data will only be loaded on demand
+        with in_directory(directory):
+            for var in self.dynamic_arrays:
+                fname = self.get_array_filename(var) + "_size"
+                new_size = int(np.loadtxt(fname, delimiter=" ", dtype=np.int32))
+                var.size = new_size
+            for var in self.dynamic_arrays_2d:
+                fname = self.get_array_filename(var) + "_size"
+                new_size = tuple(np.loadtxt(fname, delimiter=" ", dtype=np.int32))
+                var.size = new_size
+
     def build(
         self,
         directory="output",
diff --git a/brian2/devices/cpp_standalone/templates/objects.cpp b/brian2/devices/cpp_standalone/templates/objects.cpp
index a6fc2d5a0..801d023ba 100644
--- a/brian2/devices/cpp_standalone/templates/objects.cpp
+++ b/brian2/devices/cpp_standalone/templates/objects.cpp
@@ -157,15 +157,20 @@ void _write_arrays()
 	outfile_{{varname}}.open("{{get_array_filename(var) | replace('\\', '\\\\')}}", ios::binary | ios::out);
 	if(outfile_{{varname}}.is_open())
 	{
-        if (! {{varname}}.empty() )
-        {
-			outfile_{{varname}}.write(reinterpret_cast<char*>(&{{varname}}[0]), {{varname}}.size()*sizeof({{varname}}[0]));
-		    outfile_{{varname}}.close();
-		}
+		outfile_{{varname}}.write(reinterpret_cast<char*>(&{{varname}}[0]), {{varname}}.size()*sizeof({{varname}}[0]));
+		outfile_{{varname}}.close();
 	} else
 	{
 		std::cout << "Error writing output file for {{varname}}." << endl;
 	}
+	outfile_{{varname}}.open("{{get_array_filename(var) | replace('\\', '\\\\')}}_size", ios::out);
+	if (outfile_{{varname}}.is_open()) {
+	    outfile_{{varname}} << {{varname}}.size();
+	    outfile_{{varname}}.close();
+	} else
+	{
+		std::cout << "Error writing size file for {{varname}}." << endl;
+	}
 	{% endfor %}
 
 	{% for var, varname in dynamic_array_2d_specs | dictsort(by='value') %}
@@ -185,6 +190,14 @@ void _write_arrays()
 	{
 		std::cout << "Error writing output file for {{varname}}." << endl;
 	}
+	outfile_{{varname}}.open("{{get_array_filename(var) | replace('\\', '\\\\')}}_size", ios::out);
+	if (outfile_{{varname}}.is_open()) {
+	    outfile_{{varname}} << {{varname}}.n << " " << {{varname}}.m;
+	    outfile_{{varname}}.close();
+	} else
+	{
+		std::cout << "Error writing size file for {{varname}}." << endl;
+	}
 	{% endfor %}
     {% if profiled_codeobjects is defined and profiled_codeobjects %}
 	// Write profiling info to disk

From eaf1efa8e2a75bb1b8665ecdce85b5b941dba8de Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 7 Apr 2023 17:16:16 +0200
Subject: [PATCH 35/52] Simplify RateMonitor template for non-contiguous
 subgroup

No need to store the individual spikes, we only care about the total number
---
 .../cython_rt/templates/ratemonitor.pyx       | 24 +++++++------------
 .../cpp_standalone/templates/ratemonitor.cpp  | 21 +++++++---------
 2 files changed, 17 insertions(+), 28 deletions(-)

diff --git a/brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx b/brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx
index 80b852497..7b0cfab8d 100644
--- a/brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx
+++ b/brian2/codegen/runtime/cython_rt/templates/ratemonitor.pyx
@@ -2,18 +2,12 @@
 {% 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 %}
+    cdef int32_t _filtered_spikes = 0
+    cdef size_t _source_index_counter = 0
+    {% if contiguous %}{# contiguous subgroup #}
     # We assume that spikes are ordered
     _start_idx = _num_spikes
     _end_idx = _num_spikes
@@ -27,8 +21,8 @@
         if _idx < _source_stop:
             break
         _end_idx = _j
-    _num_spikes = _end_idx - _start_idx
-    {% else %}
+    _filtered_spikes = _end_idx - _start_idx
+    {% else %}{# non-contiguous subgroup #}
     for _j in range(_num_spikes):
         _idx = {{_spikespace}}[_j]
         if _idx < {{_source_indices}}[_source_index_counter]:
@@ -38,12 +32,12 @@
         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
+            _filtered_spikes += 1
+
         if _source_index_counter == {{source_N}}:
             break
-    _num_spikes = _filtered_events[{{source_N}}]
     {% endif %}
+    _num_spikes = _filtered_spikes
     {% endif %}
     
     # Calculate the new length for the arrays
diff --git a/brian2/devices/cpp_standalone/templates/ratemonitor.cpp b/brian2/devices/cpp_standalone/templates/ratemonitor.cpp
index f228a4e4e..9e22eaf12 100644
--- a/brian2/devices/cpp_standalone/templates/ratemonitor.cpp
+++ b/brian2/devices/cpp_standalone/templates/ratemonitor.cpp
@@ -4,20 +4,15 @@
 
 {% 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
+    int32_t _filtered_spikes = 0;
+    size_t _source_index_counter = 0;
     size_t _start_idx = _num_spikes;
     size_t _end_idx = _num_spikes;
     if (_num_spikes > 0)
     {
-        {% if contiguous %}
+        {% if contiguous %} {# contiguous subgroup #}
+        // We filter the spikes, making use of the fact that they are sorted
         for(size_t _j=0; _j<_num_spikes; _j++)
         {
             const int _idx = {{_spikespace}}[_j];
@@ -34,8 +29,8 @@
             }
             _end_idx = _j;
         }
-        _num_spikes = _end_idx - _start_idx;
-        {% else %}
+        _filtered_spikes = _end_idx - _start_idx;
+        {% else %} {# non-contiguous subgroup #}
         for (size_t _j=0; _j<_num_spikes; _j++)
         {
             const size_t _idx = {{_spikespace}}[_j];
@@ -49,15 +44,15 @@
                 _idx == {{_source_indices}}[_source_index_counter])
             {
                 _source_index_counter += 1;
-                _filtered_events[_filtered_events[{{source_N}}]++] = _idx;
+                _filtered_spikes += 1;
                 if (_source_index_counter == {{source_N}})
                     break;
             }
             if (_source_index_counter == {{source_N}})
                 break;
         }
-        _num_spikes = _filtered_events[{{source_N}}];
         {% endif %}
+        _num_spikes = _filtered_spikes;
     }
     {% endif %}
     {{_dynamic_rate}}.push_back(1.0*_num_spikes/{{_clock_dt}}/{{source_N}});

From 3759bf73034419370396dccca8718ecbe9cc8033 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Fri, 7 Apr 2023 17:24:10 +0200
Subject: [PATCH 36/52] Minor reformatting in test

---
 brian2/tests/test_monitor.py | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py
index 51eaca837..b23dbb953 100644
--- a/brian2/tests/test_monitor.py
+++ b/brian2/tests/test_monitor.py
@@ -17,7 +17,7 @@ def test_spike_monitor():
     G = NeuronGroup(
         3,
         """dv/dt = rate : 1
-                          rate: Hz""",
+           rate: Hz""",
         threshold="v>1",
         reset="v=0",
     )
@@ -90,8 +90,8 @@ def test_spike_monitor_variables():
     G = NeuronGroup(
         3,
         """dv/dt = rate : 1
-                          rate : Hz
-                          prev_spikes : integer""",
+           rate : Hz
+           prev_spikes : integer""",
         threshold="v>1",
         reset="v=0; prev_spikes += 1",
     )
@@ -123,8 +123,8 @@ def test_spike_monitor_get_states():
     G = NeuronGroup(
         3,
         """dv/dt = rate : 1
-                          rate : Hz
-                          prev_spikes : integer""",
+           rate : Hz
+           prev_spikes : integer""",
         threshold="v>1",
         reset="v=0; prev_spikes += 1",
     )
@@ -193,7 +193,7 @@ def test_event_monitor():
     G = NeuronGroup(
         3,
         """dv/dt = rate : 1
-                          rate: Hz""",
+           rate: Hz""",
         events={"my_event": "v>1"},
     )
     G.run_on_event("my_event", "v=0")
@@ -321,8 +321,8 @@ def test_state_monitor():
     G = NeuronGroup(
         2,
         """dv/dt = -v / (10*ms) : 1
-                          f = clip(v, 0.1, 0.9) : 1
-                          rate: Hz""",
+           f = clip(v, 0.1, 0.9) : 1
+           rate: Hz""",
         threshold="v>1",
         reset="v=0",
         refractory=2 * ms,
@@ -398,7 +398,7 @@ def test_state_monitor_subgroups():
     G = NeuronGroup(
         10,
         """v : volt
-                           step_size : volt (constant)""",
+           step_size : volt (constant)""",
     )
     G.run_regularly("v += step_size")
     G.step_size = "i*mV"
@@ -497,8 +497,8 @@ def test_state_monitor_get_states():
     G = NeuronGroup(
         2,
         """dv/dt = -v / (10*ms) : 1
-                          f = clip(v, 0.1, 0.9) : 1
-                          rate: Hz""",
+           f = clip(v, 0.1, 0.9) : 1
+           rate: Hz""",
         threshold="v>1",
         reset="v=0",
         refractory=2 * ms,
@@ -688,7 +688,7 @@ def test_rate_monitor_subgroups():
     G = NeuronGroup(
         4,
         """dv/dt = rate : 1
-                          rate : Hz""",
+           rate : Hz""",
         threshold="v>0.999",
         reset="v=0",
     )

From 2c4787b2a5a8c3b61aa368633b34d1b4d6738f6a Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Wed, 19 Apr 2023 17:13:49 +0200
Subject: [PATCH 37/52] Fail for out-of-range 1d array indices in Synapses

---
 brian2/synapses/synapses.py   | 23 ++++++++++++++++++++++-
 brian2/tests/test_synapses.py | 16 +++++++++++-----
 2 files changed, 33 insertions(+), 6 deletions(-)

diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index c1fd29a47..26ea99abf 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -614,6 +614,18 @@ def __call__(self, item=slice(None), index_var=None):  # noqa: B008
             if hasattr(item, "_indices"):
                 item = item._indices()
             final_indices = self._to_index_array(item, index_var)
+            if len(final_indices):
+                if np.min(final_indices) < 0:
+                    raise IndexError("Negative indices are not allowed.")
+                try:
+                    N = self.N.get_value()
+                    if np.max(final_indices) >= N:
+                        raise IndexError(
+                            f"Index {np.max(final_indices)} is out of bounds for "
+                            f"{N} synapses."
+                        )
+                except NotImplementedError:
+                    logger.warn("Cannot check synaptic indices for correctness")
         else:
             # 2d or 3d indexing = pre-/post-synaptic indices and (optionally) synapse number
             if len(item) == 2:  # two indices (pre- and postsynaptic cell)
@@ -636,7 +648,16 @@ def __call__(self, item=slice(None), index_var=None):  # noqa: B008
             )
 
         if index_var not in ("_idx", "0"):
-            return index_var.get_value()[final_indices.astype(np.int32)]
+            try:
+                return index_var.get_value()[final_indices.astype(np.int32)]
+            except IndexError as ex:
+                # We try to emulate numpy's indexing semantics here:
+                # slices never lead to IndexErrors, instead they return an
+                # empty array if they don't match anything
+                if isinstance(item, slice):
+                    return np.array([], dtype=np.int32)
+                else:
+                    raise ex
         else:
             return final_indices.astype(np.int32)
 
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index ce33cbed5..295880993 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -773,13 +773,17 @@ def test_state_variable_indexing_with_subgroups():
     assert len(S[:, :, 0:2:2].w) == len(G1) * len(G2)
     assert len(S[:, :, :2:2].w) == len(G1) * len(G2)
 
-    # 1d indexing is directly indexing synapses!
-    assert len(S[:].w) == len(S[0:].w)
-    assert len(S[[0, 1]].w) == len(S[3:5].w) == 2
-    assert len(S[:].w) == len(S[np.arange(len(G1) * len(G2) * 2)].w)
+    # 1d indexing is directly indexing synapses and should be equivalent to numpy syntax!
+    assert len(S[:].w) == len(S[0:].w) == len(S.w[:])
+    assert len(S[[0, 1]].w) == len(S[3:5].w) == len(S.w[:][[0, 1]]) == 2
+    assert (
+        len(S[:].w)
+        == len(S[np.arange(len(G1) * len(G2) * 2)].w)
+        == len(S.w[:][np.arange(len(G1) * len(G2) * 2)])
+    )
     assert S[3].w == S[np.int32(3)].w == S[np.int64(3)].w  # See issue #888
 
-    # Array-indexing (not yet supported for synapse index)
+    # Array-indexing
     assert_equal(S[:, 0:3].w[:], S[:, [0, 1, 2]].w[:])
     assert_equal(S[:, 0:3].w[:], S[np.arange(len(G1)), [0, 1, 2]].w[:])
     assert_equal(S[:, 0:3].w[:], S[:, [0, 1, 2], [0, 1]].w[:])
@@ -799,6 +803,8 @@ def test_state_variable_indexing_with_subgroups():
         S.__getitem__(object())
     with pytest.raises(TypeError):
         S.__getitem__(1.5)
+    with pytest.raises(IndexError):
+        S.__getitem__(np.arange(len(G1) * len(G2) * 2 + 1))
 
 
 def test_indices():

From 53ae0d87ec7e8f9f161c575fbfcedf5511fc3ed7 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 20 Apr 2023 10:51:52 +0200
Subject: [PATCH 38/52] Handle scalar synaptic indices correctly

---
 brian2/synapses/synapses.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 26ea99abf..5af19e8b4 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -614,7 +614,7 @@ def __call__(self, item=slice(None), index_var=None):  # noqa: B008
             if hasattr(item, "_indices"):
                 item = item._indices()
             final_indices = self._to_index_array(item, index_var)
-            if len(final_indices):
+            if final_indices.size > 0:
                 if np.min(final_indices) < 0:
                     raise IndexError("Negative indices are not allowed.")
                 try:

From 6d27f5364895cddfc85b891245a530f4ce6a7e98 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 20 Apr 2023 13:43:05 +0200
Subject: [PATCH 39/52] Raise errors for boolean indices of incorrect shape

---
 brian2/groups/group.py           | 5 +++++
 brian2/tests/test_neurongroup.py | 9 ++++++++-
 2 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/brian2/groups/group.py b/brian2/groups/group.py
index 639b51a8a..54009def2 100644
--- a/brian2/groups/group.py
+++ b/brian2/groups/group.py
@@ -320,6 +320,11 @@ def _to_index_array(self, item, index_var):
         else:  # array, sequence, or single value
             index_array = np.asarray(item)
             if index_array.dtype == bool:
+                if not index_array.shape == (self.N.get_value(),):
+                    raise IndexError(
+                        "Boolean index did not match shape of indexed array;shape is"
+                        f" {(self.N.get_value(), )}, but index is {index_array.shape}"
+                    )
                 index_array = np.flatnonzero(index_array)
             elif not np.issubdtype(index_array.dtype, np.signedinteger):
                 raise TypeError(
diff --git a/brian2/tests/test_neurongroup.py b/brian2/tests/test_neurongroup.py
index b307957a9..a38ce93c3 100644
--- a/brian2/tests/test_neurongroup.py
+++ b/brian2/tests/test_neurongroup.py
@@ -1175,19 +1175,26 @@ def test_state_variable_access():
     assert_allclose(np.asarray(G.v[:]), np.arange(10))
     assert have_same_dimensions(G.v[:], volt)
     assert_allclose(np.asarray(G.v[:]), G.v_[:])
-    # Accessing single elements, slices and arrays
+    # Accessing single elements, slices and integer and boolean arrays
     assert G.v[5] == 5 * volt
     assert G.v_[5] == 5
     assert_allclose(G.v[:5], np.arange(5) * volt)
     assert_allclose(G.v_[:5], np.arange(5))
     assert_allclose(G.v[[0, 5]], [0, 5] * volt)
     assert_allclose(G.v_[[0, 5]], np.array([0, 5]))
+    assert_allclose(G.v[[True, False] * 5], np.arange(10)[::2] * volt)
 
     # Illegal indexing
     with pytest.raises(IndexError):
         G.v[0, 0]
     with pytest.raises(IndexError):
         G.v_[0, 0]
+    with pytest.raises(IndexError):
+        G.v[[0, 10]]  # out of range array indices
+    with pytest.raises(IndexError):
+        G.v[[True, False]]  # too few boolean indices
+    with pytest.raises(IndexError):
+        G.v[[True] * 11]  # too many boolean indices
     with pytest.raises(TypeError):
         G.v[object()]
     with pytest.raises(TypeError):

From adc17ec90d782e58f07f314ceea3791a813ce5df Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 20 Apr 2023 16:02:57 +0200
Subject: [PATCH 40/52] Handle IndexErrors correctly for synapses as well

---
 brian2/synapses/synapses.py   | 36 +++++++++++++++++------------------
 brian2/tests/test_subgroup.py |  2 --
 brian2/tests/test_synapses.py | 18 ++++++++++++++++++
 3 files changed, 35 insertions(+), 21 deletions(-)

diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 5af19e8b4..5f0a89a48 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -663,24 +663,18 @@ def __call__(self, item=slice(None), index_var=None):  # noqa: B008
 
     def _from_3d_indices_to_index_array(self, i_indices, j_indices, k_indices):
         # Convert to absolute indices (e.g. for subgroups)
-        # Allow the indexing to fail, we'll later return an empty array in
-        # that case
-        try:
-            if hasattr(i_indices, "_indices"):  # will return absolute indices already
-                i_indices = i_indices._indices()
-            else:
-                i_indices = self.source._indices(i_indices)
-            pre_synapses = find_synapses(i_indices, self.synaptic_pre.get_value())
-        except IndexError:
-            pre_synapses = np.array([], dtype=np.int32)
-        try:
-            if hasattr(j_indices, "_indices"):
-                j_indices = j_indices._indices()
-            else:
-                j_indices = self.target._indices(j_indices)
-            post_synapses = find_synapses(j_indices, self.synaptic_post.get_value())
-        except IndexError:
-            post_synapses = np.array([], dtype=np.int32)
+        if hasattr(i_indices, "_indices"):  # will return absolute indices already
+            i_indices = i_indices._indices()
+        else:
+            i_indices = self.source._indices(i_indices)
+        pre_synapses = find_synapses(i_indices, self.synaptic_pre.get_value())
+
+        if hasattr(j_indices, "_indices"):
+            j_indices = j_indices._indices()
+        else:
+            j_indices = self.target._indices(j_indices)
+        post_synapses = find_synapses(j_indices, self.synaptic_post.get_value())
+
         matching_synapses = np.intersect1d(
             pre_synapses, post_synapses, assume_unique=True
         )
@@ -690,11 +684,15 @@ def _from_3d_indices_to_index_array(self, i_indices, j_indices, k_indices):
         else:
             # We want to access the raw arrays here, not go through the Variable
             synapse_numbers = self.synapse_number.get_value()[matching_synapses]
-
             if isinstance(k_indices, (numbers.Integral, slice)):
                 test_k = slice_to_test(k_indices)
                 final_indices = matching_synapses[test_k(synapse_numbers)]
             else:
+                k_indices = np.asarray(k_indices)
+                if k_indices.dtype == bool:
+                    raise NotImplementedError(
+                        "Boolean indexing not supported for synapse number"
+                    )
                 final_indices = matching_synapses[np.in1d(synapse_numbers, k_indices)]
         return final_indices
 
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 8036a676c..7e0389dc8 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -516,7 +516,6 @@ def test_synapse_access():
     S.w = "2*j"
     assert all(S.w[:, 1] == 2)
 
-    assert len(S.w[:, 10]) == 0
     assert len(S.w["j==10"]) == 0
 
     # Test referencing pre- and postsynaptic variables
@@ -531,7 +530,6 @@ def test_synapse_access():
     assert len(S) == len(S.w[SG1, SG2])
     assert_equal(S.w[SG1, 1], S.w[:, 1])
     assert_equal(S.w[1, SG2], S.w[1, :])
-    assert len(S.w[SG1, 10]) == 0
 
 
 def test_synapses_access_subgroups():
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index 295880993..718a48173 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -727,6 +727,18 @@ def test_state_variable_indexing():
     assert_equal(S.w[:, 0:3], S.w[:, [0, 1, 2], [0, 1]])
     assert_equal(S.w[:, 0:3, 0], S.w[:, [0, 1, 2], [0]])
 
+    # Array indexing with boolean arrays
+    assert_equal(
+        S.w[:, 0:3], S.w[:, np.array([True, True, True, False, False, False, False])]
+    )
+    assert_equal(
+        S.w[:, 0:3],
+        S.w[
+            [True, True, True, True, True],
+            [True, True, True, False, False, False, False],
+        ],
+    )
+
     # string-based indexing
     assert_equal(S.w[0:3, :], S.w["i<3"])
     assert_equal(S.w[:, 0:3], S.w["j<3"])
@@ -737,6 +749,12 @@ def test_state_variable_indexing():
     # invalid indices
     with pytest.raises(IndexError):
         S.w.__getitem__((1, 2, 3, 4))
+    with pytest.raises(IndexError):
+        S.w[[0, 5], :]  # out-of-range array index
+    with pytest.raises(IndexError):
+        S.w[[True, False], :]  # too few boolean indices
+    with pytest.raises(IndexError):
+        S.w[[True, False, True, False, True, False], :]  # too many boolean indices
     with pytest.raises(TypeError):
         S.w.__getitem__(object())
     with pytest.raises(TypeError):

From fed868e8bb732249e4819973245c81b4270c2114 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 20 Apr 2023 22:07:19 +0200
Subject: [PATCH 41/52] Install coveralls with pip everywhere

---
 .github/workflows/testsuite.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/testsuite.yml b/.github/workflows/testsuite.yml
index 61000adfb..da5718c40 100644
--- a/.github/workflows/testsuite.yml
+++ b/.github/workflows/testsuite.yml
@@ -91,7 +91,7 @@ jobs:
         if: ${{ startsWith(matrix.os, 'ubuntu-') && matrix.python-version == needs.get_python_versions.outputs.max-python }}
         run: |
           cp  $GITHUB_WORKSPACE/../.coverage .
-          conda install --quiet --yes coveralls
+          pip3 install --upgrade coveralls
           coveralls --service=github --rcfile=$GITHUB_WORKSPACE/.coveragerc
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
@@ -107,7 +107,7 @@ jobs:
     - name: Finished
       run: |
         pip3 install --upgrade coveralls
-        coveralls --finish
+        coveralls --service=github --finish
       env:
         GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 

From a7f6cbe193bb289e615585d02570fd03289fe4e3 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@inserm.fr>
Date: Thu, 15 Jun 2023 21:52:56 +0200
Subject: [PATCH 42/52] Fix incorrect merge

---
 brian2/groups/group.py | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/brian2/groups/group.py b/brian2/groups/group.py
index d664e0e42..bd36e655c 100644
--- a/brian2/groups/group.py
+++ b/brian2/groups/group.py
@@ -257,18 +257,25 @@ def __call__(self, item=slice(None), index_var=None):  # noqa: B008
             raise IndexError(
                 f"Can only interpret 1-d indices, got {len(item)} dimensions."
             )
-        else:
-            if isinstance(item, str) and item == "True":
-                item = slice(None)
-            if isinstance(item, slice):
-                if index_var == "0":
-                    return 0
-                if index_var == "_idx":
-                    start, stop, step = item.indices(self.N.item())
+
+        index_array = self._to_index_array(item, index_var)
+
+        if index_array.size == 0:
+            return index_array
+
+        if index_var not in ("_idx", "0"):
+            try:
+                index_array = index_var.get_value()[index_array]
+            except IndexError as ex:
+                # We try to emulate numpy's indexing semantics here:
+                # slices never lead to IndexErrors, instead they return an
+                # empty array if they don't match anything
+                if isinstance(item, slice):
+                    return np.array([], dtype=np.int32)
                 else:
                     raise ex
         else:
-            N = int(self.N.get_value())
+            N = self.N.item()
             if np.min(index_array) < -N:
                 raise IndexError(
                     "Illegal index {} for a group of size {}".format(

From a7cd7f3e87bfbc563b3389e7742d49309b93635f Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Thu, 26 Oct 2023 18:43:16 +0200
Subject: [PATCH 43/52] Do not force float64 in TimedArray

---
 brian2/input/timedarray.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/brian2/input/timedarray.py b/brian2/input/timedarray.py
index 447b92703..ee1e4addc 100644
--- a/brian2/input/timedarray.py
+++ b/brian2/input/timedarray.py
@@ -228,12 +228,21 @@ class TimedArray(Function, Nameable, CacheKey):
 
     @check_units(dt=second)
     def __init__(self, values, dt, name=None):
+        from brian2.core.preferences import prefs
+
         if name is None:
             name = "_timedarray*"
         Nameable.__init__(self, name)
         dimensions = get_dimensions(values)
         self.dim = dimensions
-        values = np.asarray(values, dtype=np.float64)
+        values = np.asarray(values)  # infer dtype
+        if values.dtype == np.object:
+            raise TypeError("TimedArray does not support arrays with dtype 'object'")
+        elif (
+            values.dtype == np.float64 and prefs.core.default_float_dtype != np.float64
+        ):
+            # Reduce the precision of the values array to the default scalar type
+            values = values.astype(prefs.core.default_float_dtype)
         self.values = values
         dt = float(dt)
         self.dt = dt

From 437e3fe0f478e330caca7bc536252bd38c655e37 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Thu, 26 Oct 2023 18:44:29 +0200
Subject: [PATCH 44/52] Invalidate caches on each run

---
 brian2/devices/cpp_standalone/device.py | 20 ++++++++++++++++----
 brian2/tests/test_cpp_standalone.py     |  5 +++--
 2 files changed, 19 insertions(+), 6 deletions(-)

diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index cb28a5078..17c25e068 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/brian2/devices/cpp_standalone/device.py
@@ -246,6 +246,8 @@ def __init__(self):
             self.libraries += ["advapi32"]
         self.extra_link_args = []
         self.writer = None
+        # List of variables to invalidate
+        self.invalidate_on_run = []
 
     def reinit(self):
         # Remember the build_on_run setting and its options -- important during
@@ -785,7 +787,7 @@ def code_object(
                 and var not in do_not_invalidate
                 and (not var.read_only or var in written_readonly_vars)
             ):
-                self.array_cache[var] = None
+                self.invalidate_on_run.append(var)
 
         return codeobj
 
@@ -1315,7 +1317,6 @@ def run(
                 list_rep.append(f"{name}={string_value}")
 
             run_args = list_rep
-
         # Invalidate array cache for all variables set on the command line
         for arg in run_args:
             s = arg.split("=")
@@ -1326,6 +1327,9 @@ def run(
                         and var.owner.name + "." + var.name == s[0]
                     ):
                         self.array_cache[var] = None
+        # Invalidate array cache for all variables written by code objects
+        for var in self.invalidate_on_run:
+            self.array_cache[var] = None
         run_args = ["--results_dir", self.results_dir] + run_args
         # Invalidate the cached end time of the clock and network, to deal with stopped simulations
         for clock in self.clocks:
@@ -1416,7 +1420,11 @@ def _prepare_variableview_run_arg(self, key, value):
         value_ar = np.asarray(value, dtype=key.dtype)
         if value_ar.ndim == 0 or value_ar.size == 1:
             # single value, give value directly on command line
-            string_value = repr(value_ar.item())
+            value = value_ar.item()
+            if value is True or value is None:  # translate True/False to 1/0
+                string_value = "1" if value else "0"
+            else:
+                string_value = repr(value)
             value_name = None
         else:
             if value_ar.ndim != 1 or (
@@ -1439,7 +1447,11 @@ def _prepare_timed_array_run_arg(self, key, value):
         value_ar = np.asarray(value, dtype=key.values.dtype)
         if value_ar.ndim == 0 or value_ar.size == 1:
             # single value, give value directly on command line
-            string_value = repr(value_ar.item())
+            value = value_ar.item()
+            if value is True or value is None:  # translate True/False to 1/0
+                string_value = "1" if value else "0"
+            else:
+                string_value = repr(value)
             value_name = None
         elif value_ar.shape == key.values.shape:
             value_name = f"init_{key.name}_values_{md5(value_ar.data).hexdigest()}.dat"
diff --git a/brian2/tests/test_cpp_standalone.py b/brian2/tests/test_cpp_standalone.py
index 0d4db56b8..f593e4e51 100644
--- a/brian2/tests/test_cpp_standalone.py
+++ b/brian2/tests/test_cpp_standalone.py
@@ -844,6 +844,7 @@ def test_change_parameter_without_recompile_dict_syntax():
             G.n: ar,
             G.b: np.arange(10) % 2 != 0,
             stim: ar2,
+            on_off: [False, True, False],
         }
     )
     assert array_equal(G.x, ar)
@@ -854,9 +855,9 @@ def test_change_parameter_without_recompile_dict_syntax():
         mon.s.T / nA,
         np.array(
             [
-                [0, 2, 4, 6, 8, 10, 12, 14, 16, 18],
                 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # on_off(t) == False
-                [40, 42, 44, 46, 48, 50, 52, 54, 56, 58],
+                [20, 22, 24, 26, 28, 30, 32, 34, 36, 38],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # on_off(t) == False
             ]
         ),
     )

From b476451f3b2869c5cf014d98493b0cec49a9b237 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Mon, 30 Oct 2023 12:09:14 +0100
Subject: [PATCH 45/52] Fix cache invalidation

---
 brian2/devices/cpp_standalone/codeobject.py |  4 ++++
 brian2/devices/cpp_standalone/device.py     | 16 +++++++++++-----
 2 files changed, 15 insertions(+), 5 deletions(-)

diff --git a/brian2/devices/cpp_standalone/codeobject.py b/brian2/devices/cpp_standalone/codeobject.py
index 9778c538d..8c83372cc 100644
--- a/brian2/devices/cpp_standalone/codeobject.py
+++ b/brian2/devices/cpp_standalone/codeobject.py
@@ -108,8 +108,12 @@ def __init__(self, *args, **kwds):
         super().__init__(*args, **kwds)
         #: Store whether this code object defines before/after blocks
         self.before_after_blocks = []
+        #: Store variables that are updated by this code object and therefore invalidate the cache
+        self.invalidate_cache_variables = set()
 
     def __call__(self, **kwds):
+        for var in self.invalidate_cache_variables:
+            get_device().array_cache[var] = None
         return self.run()
 
     def compile_block(self, block):
diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index 17c25e068..4465a23e7 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/brian2/devices/cpp_standalone/device.py
@@ -246,8 +246,6 @@ def __init__(self):
             self.libraries += ["advapi32"]
         self.extra_link_args = []
         self.writer = None
-        # List of variables to invalidate
-        self.invalidate_on_run = []
 
     def reinit(self):
         # Remember the build_on_run setting and its options -- important during
@@ -787,7 +785,7 @@ def code_object(
                 and var not in do_not_invalidate
                 and (not var.read_only or var in written_readonly_vars)
             ):
-                self.invalidate_on_run.append(var)
+                codeobj.invalidate_cache_variables.add(var)
 
         return codeobj
 
@@ -1327,9 +1325,12 @@ def run(
                         and var.owner.name + "." + var.name == s[0]
                     ):
                         self.array_cache[var] = None
+
         # Invalidate array cache for all variables written by code objects
-        for var in self.invalidate_on_run:
-            self.array_cache[var] = None
+        for codeobj in self.code_objects.values():
+            for var in codeobj.invalidate_cache_variables:
+                self.array_cache[var] = None
+
         run_args = ["--results_dir", self.results_dir] + run_args
         # Invalidate the cached end time of the clock and network, to deal with stopped simulations
         for clock in self.clocks:
@@ -2006,6 +2007,11 @@ def network_run(
 
         net.after_run()
 
+        # Invalidate array cache for all variables written by code objects
+        for _, codeobj in code_objects:
+            for var in codeobj.invalidate_cache_variables:
+                self.array_cache[var] = None
+
         # Manually set the cache for the clocks, simulation scripts might
         # want to access the time (which has been set in code and is therefore
         # not accessible by the normal means until the code has been built and

From 285dffb7aa0720b1998faaee47e9cc3b361e9be8 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Mon, 30 Oct 2023 12:11:10 +0100
Subject: [PATCH 46/52] Correctly the return type of TimedArray functions in
 generated code

---
 brian2/input/timedarray.py | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/brian2/input/timedarray.py b/brian2/input/timedarray.py
index ee1e4addc..4bdb27e16 100644
--- a/brian2/input/timedarray.py
+++ b/brian2/input/timedarray.py
@@ -5,6 +5,7 @@
 
 import numpy as np
 
+from brian2.codegen.generators import c_data_type
 from brian2.core.clocks import defaultclock
 from brian2.core.functions import Function
 from brian2.core.names import Nameable
@@ -45,7 +46,7 @@ def cpp_impl(owner):
         K = _find_K(owner.clock.dt_, dt)
         code = (
             """
-        static inline double %NAME%(const double t)
+        static inline %TYPE% %NAME%(const double t)
         {
             const double epsilon = %DT% / %K%;
             int i = (int)((t/epsilon + 0.5)/%K%);
@@ -61,6 +62,7 @@ def cpp_impl(owner):
             .replace("%DT%", f"{dt:.18f}")
             .replace("%K%", str(K))
             .replace("%NUM_VALUES%", str(len(values)))
+            .replace("%TYPE%", c_data_type(values.dtype))
         )
 
         return code
@@ -72,7 +74,7 @@ def _generate_cpp_code_2d(values, dt, name):
     def cpp_impl(owner):
         K = _find_K(owner.clock.dt_, dt)
         support_code = """
-        static inline double %NAME%(const double t, const int i)
+        static inline %TYPE% %NAME%(const double t, const int i)
         {
             const double epsilon = %DT% / %K%;
             if (i < 0 || i >= %COLS%)
@@ -93,6 +95,7 @@ def cpp_impl(owner):
                 "%K%": str(K),
                 "%COLS%": str(values.shape[1]),
                 "%ROWS%": str(values.shape[0]),
+                "%TYPE%": c_data_type(values.dtype),
             },
         )
         return code
@@ -105,7 +108,7 @@ def cython_impl(owner):
         K = _find_K(owner.clock.dt_, dt)
         code = (
             """
-        cdef double %NAME%(const double t):
+        cdef %TYPE% %NAME%(const double t):
             global _namespace%NAME%_values
             cdef double epsilon = %DT% / %K%
             cdef int i = (int)((t/epsilon + 0.5)/%K%)
@@ -120,6 +123,7 @@ def cython_impl(owner):
             .replace("%DT%", f"{dt:.18f}")
             .replace("%K%", str(K))
             .replace("%NUM_VALUES%", str(len(values)))
+            .replace("%TYPE%", c_data_type(values.dtype))
         )
 
         return code
@@ -131,7 +135,7 @@ def _generate_cython_code_2d(values, dt, name):
     def cython_impl(owner):
         K = _find_K(owner.clock.dt_, dt)
         code = """
-        cdef double %NAME%(const double t, const int i):
+        cdef %TYPE% %NAME%(const double t, const int i):
             global _namespace%NAME%_values
             cdef double epsilon = %DT% / %K%
             if i < 0 or i >= %COLS%:
@@ -151,6 +155,7 @@ def cython_impl(owner):
                 "%K%": str(K),
                 "%COLS%": str(values.shape[1]),
                 "%ROWS%": str(values.shape[0]),
+                "%TYPE%": c_data_type(values.dtype),
             },
         )
         return code
@@ -236,7 +241,7 @@ def __init__(self, values, dt, name=None):
         dimensions = get_dimensions(values)
         self.dim = dimensions
         values = np.asarray(values)  # infer dtype
-        if values.dtype == np.object:
+        if values.dtype == object:
             raise TypeError("TimedArray does not support arrays with dtype 'object'")
         elif (
             values.dtype == np.float64 and prefs.core.default_float_dtype != np.float64
@@ -347,7 +352,9 @@ def unitless_timed_array_func(t, i):
         self.implementations.add_dynamic_implementation(
             "numpy", create_numpy_implementation
         )
-        values_flat = self.values.astype(np.double, order="C", copy=False).ravel()
+        values_flat = self.values.astype(
+            self.values.dtype, order="C", copy=False
+        ).ravel()
         namespace = lambda owner: {f"{self.name}_values": values_flat}
 
         for target, (_, func_2d) in TimedArray.implementations.items():

From 4cb865b5e3bc66c29ace16a3e4282b2a8e8f96b3 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Thu, 2 Nov 2023 18:33:12 +0100
Subject: [PATCH 47/52] WIP: Allow non-contiguous subgroups for synapses

Only C++ standalone supported until now
---
 .../templates/synapses_create_array.cpp       |   8 ++
 .../templates/synapses_create_generator.cpp   |  16 +++
 brian2/synapses/synapses.py                   | 110 +++++++++++++-----
 brian2/tests/test_subgroup.py                 |  48 ++++++++
 brian2/tests/test_synapses.py                 |   6 -
 5 files changed, 150 insertions(+), 38 deletions(-)

diff --git a/brian2/devices/cpp_standalone/templates/synapses_create_array.cpp b/brian2/devices/cpp_standalone/templates/synapses_create_array.cpp
index 2d652bc70..b13881ab5 100644
--- a/brian2/devices/cpp_standalone/templates/synapses_create_array.cpp
+++ b/brian2/devices/cpp_standalone/templates/synapses_create_array.cpp
@@ -16,8 +16,16 @@ const size_t _new_num_synapses = _old_num_synapses + _numsources;
 constants or scalar arrays#}
 const size_t _N_pre = {{constant_or_scalar('N_pre', variables['N_pre'])}};
 const size_t _N_post = {{constant_or_scalar('N_post', variables['N_post'])}};
+{% if "_target_sub_idx" in variables %}
+{{_dynamic_N_incoming}}.resize({{get_array_name(variables['_target_sub_idx'])}}[_num_target_sub_idx - 1] + 1);
+{% else %}
 {{_dynamic_N_incoming}}.resize(_N_post + _target_offset);
+{% endif %}
+{% if "_source_sub_idx" in variables %}
+{{_dynamic_N_outgoing}}.resize({{get_array_name(variables['_source_sub_idx'])}}[_num_source_sub_idx - 1] + 1);
+{% else %}
 {{_dynamic_N_outgoing}}.resize(_N_pre + _source_offset);
+{% endif %}
 
 for (size_t _idx=0; _idx<_numsources; _idx++) {
     {# After this code has been executed, the arrays _real_sources and
diff --git a/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp b/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp
index 04cfe1c17..bd2fab5a8 100644
--- a/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp
+++ b/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp
@@ -18,8 +18,16 @@
     constants or scalar arrays#}
     const size_t _N_pre = {{constant_or_scalar('N_pre', variables['N_pre'])}};
     const size_t _N_post = {{constant_or_scalar('N_post', variables['N_post'])}};
+    {% if "_target_sub_idx" in variables %}
+    {{_dynamic_N_incoming}}.resize({{get_array_name(variables['_target_sub_idx'])}}[_num_target_sub_idx - 1] + 1);
+    {% else %}
     {{_dynamic_N_incoming}}.resize(_N_post + _target_offset);
+    {% endif %}
+    {% if "_source_sub_idx" in variables %}
+    {{_dynamic_N_outgoing}}.resize({{get_array_name(variables['_source_sub_idx'])}}[_num_source_sub_idx - 1] + 1);
+    {% else %}
     {{_dynamic_N_outgoing}}.resize(_N_pre + _source_offset);
+    {% endif %}
     size_t _raw_pre_idx, _raw_post_idx;
     {# For a connect call j='k+i for k in range(0, N_post, 2) if k+i < N_post'
     "j" is called the "result index" (and "_post_idx" the "result index array", etc.)
@@ -35,7 +43,11 @@
     for(size_t _{{outer_index}}=0; _{{outer_index}}<_{{outer_index_size}}; _{{outer_index}}++)
     {
         bool __cond, _cond;
+        {% if outer_contiguous %}
         _raw{{outer_index_array}} = _{{outer_index}} + {{outer_index_offset}};
+        {% else %}
+        _raw{{outer_index_array}} = {{get_array_name(variables[outer_sub_idx])}}[_{{outer_index}}];
+        {% endif %}
         {% if not result_index_condition %}
         {
             {{vector_code['create_cond']|autoindent}}
@@ -181,7 +193,11 @@
             }
             _{{result_index}} = __{{result_index}}; // make the previously locally scoped var available
             {{outer_index_array}} = _{{outer_index_array}};
+            {% if result_contiguous %}
             _raw{{result_index_array}} = _{{result_index}} + {{result_index_offset}};
+            {% else %}
+            _raw{{result_index_array}} = {{get_array_name(variables[result_sub_idx])}}[_{{result_index}}];
+            {% endif %}
             {% if result_index_condition %}
             {
                 {% if result_index_used %}
diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index c051633c1..88b211edc 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -1363,45 +1363,55 @@ def _create_variables(self, equations, user_dtype=None):
         self.variables.add_reference("_presynaptic_idx", self, "_synaptic_pre")
         self.variables.add_reference("_postsynaptic_idx", self, "_synaptic_post")
 
-        # Except for subgroups (which potentially add an offset), the "i" and
-        # "j" variables are simply equivalent to `_synaptic_pre` and
-        # `_synaptic_post`
-        if getattr(self.source, "start", 0) == 0:
-            self.variables.add_reference("i", self, "_synaptic_pre")
-        else:
-            if isinstance(self.source, Subgroup) and not self.source.contiguous:
-                raise TypeError(
-                    "Cannot use a non-contiguous subgroup as a "
-                    "source group for Synapses."
-                )
+        # Except for subgroups, the "i" and "j" variables are simply equivalent to
+        # `_synaptic_pre` and `_synaptic_post`
+        if (
+            isinstance(self.source, Subgroup)
+            and not getattr(self.source, "start", -1) == 0
+        ):
             self.variables.add_reference(
                 "_source_i", self.source.source, "i", index="_presynaptic_idx"
             )
-            self.variables.add_reference("_source_offset", self.source, "_offset")
-            self.variables.add_subexpression(
-                "i",
-                dtype=self.source.source.variables["i"].dtype,
-                expr="_source_i - _source_offset",
-                index="_presynaptic_idx",
-            )
-        if getattr(self.target, "start", 0) == 0:
-            self.variables.add_reference("j", self, "_synaptic_post")
-        else:
-            if isinstance(self.target, Subgroup) and not self.target.contiguous:
-                raise TypeError(
-                    "Cannot use a non-contiguous subgroup as a "
-                    "target group for Synapses."
+            if getattr(self.source, "contiguous", True):
+                # Contiguous subgroup simply shift the indices
+                self.variables.add_subexpression(
+                    "i",
+                    dtype=self.source.source.variables["i"].dtype,
+                    expr="_source_i - _source_offset",
+                    index="_presynaptic_idx",
+                )
+            else:
+                # Non-contiguous subgroups need a full translation
+                self.variables.add_reference(
+                    "i", self.source, "i", index="_presynaptic_idx"
                 )
+        else:
+            # No subgroup or subgroup starting at 0
+            self.variables.add_reference("i", self, "_synaptic_pre")
+
+        if (
+            isinstance(self.target, Subgroup)
+            and not getattr(self.target, "start", -1) == 0
+        ):
             self.variables.add_reference(
                 "_target_j", self.target.source, "i", index="_postsynaptic_idx"
             )
-            self.variables.add_reference("_target_offset", self.target, "_offset")
-            self.variables.add_subexpression(
-                "j",
-                dtype=self.target.source.variables["i"].dtype,
-                expr="_target_j - _target_offset",
-                index="_postsynaptic_idx",
-            )
+            if getattr(self.target, "contiguous", True):
+                # Contiguous subgroup simply shift the indices
+                self.variables.add_subexpression(
+                    "j",
+                    dtype=self.target.source.variables["i"].dtype,
+                    expr="_target_j - _target_offset",
+                    index="_postsynaptic_idx",
+                )
+            else:
+                # Non-contiguous subgroups need a full translation
+                self.variables.add_reference(
+                    "j", self.target, "i", index="_postsynaptic_idx"
+                )
+        else:
+            # No subgroup or subgroup starting at 0
+            self.variables.add_reference("j", self, "_synaptic_post")
 
         # Add the standard variables
         self.variables.add_array(
@@ -1934,11 +1944,21 @@ def _add_synapses_from_arrays(self, sources, targets, n, p, namespace=None):
         if "_offset" in self.source.variables:
             variables.add_reference("_source_offset", self.source, "_offset")
             abstract_code += "_real_sources = sources + _source_offset\n"
+        elif not getattr(self.source, "contiguous", True):
+            variables.add_reference(
+                "_source_sub_idx", self.source, "_sub_idx", index="sources"
+            )
+            abstract_code += "_real_sources = _source_sub_idx\n"
         else:
             abstract_code += "_real_sources = sources\n"
         if "_offset" in self.target.variables:
             variables.add_reference("_target_offset", self.target, "_offset")
             abstract_code += "_real_targets = targets + _target_offset\n"
+        elif not getattr(self.target, "contiguous", True):
+            variables.add_reference(
+                "_target_sub_idx", self.target, "_sub_idx", index="targets"
+            )
+            abstract_code += "_real_targets = _target_sub_idx\n"
         else:
             abstract_code += "_real_targets = targets"
         logger.debug(
@@ -2022,11 +2042,25 @@ def _add_synapses_generator(
         outer_index_size = "N_pre" if over_presynaptic else "N_post"
         outer_index_array = "_pre_idx" if over_presynaptic else "_post_idx"
         outer_index_offset = "_source_offset" if over_presynaptic else "_target_offset"
+        outer_sub_idx = "_source_sub_idx" if over_presynaptic else "_target_sub_idx"
+        outer_contiguous = (
+            getattr(self.source, "contiguous", True)
+            if over_presynaptic
+            else getattr(self.target, "contiguous", True)
+        )
+
         result_index = "j" if over_presynaptic else "i"
         result_index_size = "N_post" if over_presynaptic else "N_pre"
         target_idx = "_postsynaptic_idx" if over_presynaptic else "_presynaptic_idx"
         result_index_array = "_post_idx" if over_presynaptic else "_pre_idx"
         result_index_offset = "_target_offset" if over_presynaptic else "_source_offset"
+        result_sub_idx = "_target_sub_idx" if over_presynaptic else "_source_sub_idx"
+        result_contiguous = (
+            getattr(self.target, "contiguous", True)
+            if over_presynaptic
+            else getattr(self.source, "contiguous", True)
+        )
+
         result_index_name = "postsynaptic" if over_presynaptic else "presynaptic"
         template_kwds.update(
             {
@@ -2034,11 +2068,15 @@ def _add_synapses_generator(
                 "outer_index_size": outer_index_size,
                 "outer_index_array": outer_index_array,
                 "outer_index_offset": outer_index_offset,
+                "outer_sub_idx": outer_sub_idx,
+                "outer_contiguous": outer_contiguous,
                 "result_index": result_index,
                 "result_index_size": result_index_size,
                 "result_index_name": result_index_name,
                 "result_index_array": result_index_array,
                 "result_index_offset": result_index_offset,
+                "result_sub_idx": result_sub_idx,
+                "result_contiguous": result_contiguous,
             }
         )
         abstract_code = {
@@ -2126,11 +2164,19 @@ def _add_synapses_generator(
         else:
             variables.add_constant("_source_offset", value=0)
 
+        if not getattr(self.source, "contiguous", True):
+            variables.add_reference("_source_sub_idx", self.source, "_sub_idx")
+            needed_variables.append("_source_sub_idx")
+
         if "_offset" in self.target.variables:
             variables.add_reference("_target_offset", self.target, "_offset")
         else:
             variables.add_constant("_target_offset", value=0)
 
+        if not getattr(self.target, "contiguous", True):
+            variables.add_reference("_target_sub_idx", self.target, "_sub_idx")
+            needed_variables.append("_target_sub_idx")
+
         variables.add_auxiliary_variable("_raw_pre_idx", dtype=np.int32)
         variables.add_auxiliary_variable("_raw_post_idx", dtype=np.int32)
 
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 7e0389dc8..69e125474 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -378,6 +378,54 @@ def test_synapse_creation_generator():
     assert all(S5.v_pre[:] < 25)
 
 
+@pytest.mark.standalone_compatible
+def test_synapse_creation_generator_non_contiguous():
+    G1 = NeuronGroup(10, "v:1")
+    G2 = NeuronGroup(20, "v:1")
+    G1.v = "i"
+    G2.v = "10 + i"
+    SG1 = G1[[0, 2, 4, 6, 8]]
+    SG2 = G2[1::2]
+    S = Synapses(SG1, SG2, "w:1")
+    S.connect(j="i*2 + k for k in range(2)")  # diverging connections
+
+    # connect based on pre-/postsynaptic state variables
+    S2 = Synapses(SG1, SG2, "w:1")
+    S2.connect(j="k for k in range(N_post) if v_pre > 2")
+
+    S3 = Synapses(SG1, SG2, "w:1")
+    S3.connect(j="k for k in range(N_post) if v_post < 25")
+
+    S4 = Synapses(SG2, SG1, "w:1")
+    S4.connect(j="k for k in range(N_post) if v_post > 2")
+
+    S5 = Synapses(SG2, SG1, "w:1")
+    S5.connect(j="k for k in range(N_post) if v_pre < 25")
+
+    run(0 * ms)  # for standalone
+
+    # Internally, the "real" neuron indices should be used
+    assert_equal(S._synaptic_pre[:], np.arange(0, 10, 2).repeat(2))
+    assert_equal(S._synaptic_post[:], np.arange(1, 20, 2))
+    # For the user, the subgroup-relative indices should be presented
+    assert_equal(S.i[:], np.arange(5).repeat(2))
+    assert_equal(S.j[:], np.arange(10))
+
+    # N_incoming and N_outgoing should also be correct
+    assert all(S.N_outgoing[:] == 2)
+    assert all(S.N_incoming[:] == 1)
+
+    assert len(S2) == 3 * len(SG2), str(len(S2))
+    assert all(S2.v_pre[:] > 2)
+    assert len(S3) == 7 * len(SG1), f"{len(S3)} != {7 * len(SG1)} "
+    assert all(S3.v_post[:] < 25)
+
+    assert len(S4) == 3 * len(SG2), str(len(S4))
+    assert all(S4.v_post[:] > 2)
+    assert len(S5) == 7 * len(SG1), f"{len(S5)} != {7 * len(SG1)} "
+    assert all(S5.v_pre[:] < 25)
+
+
 @pytest.mark.standalone_compatible
 def test_synapse_creation_generator_multiple_synapses():
     G1 = NeuronGroup(10, "v:1")
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index 718a48173..fbc5b3a73 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -88,12 +88,6 @@ def test_creation_errors():
         Synapses(G, G, "w:1", pre="v+=w", on_pre="v+=w")
     with pytest.raises(TypeError):
         Synapses(G, G, "w:1", post="v+=w", on_post="v+=w")
-    # We do not allow non-contiguous subgroups as source/target groups at the
-    # moment
-    with pytest.raises(TypeError):
-        Synapses(G[::2], G, "")
-    with pytest.raises(TypeError):
-        Synapses(G, G[::2], "")
 
 
 @pytest.mark.codegen_independent

From bb7ddb98204e0e8aeaf8f11b3a0d8909b5e3e4c9 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Tue, 14 Nov 2023 17:15:52 +0100
Subject: [PATCH 48/52] Allow non-contiguous subgroups for synapses

Implementation for Cython and numpy
---
 .../templates/synapses_create_generator.pyx   | 10 +++++--
 .../templates/synapses_create_generator.py_   | 14 +++++++++-
 brian2/synapses/synapses.py                   | 28 ++++++++++---------
 3 files changed, 36 insertions(+), 16 deletions(-)

diff --git a/brian2/codegen/runtime/cython_rt/templates/synapses_create_generator.pyx b/brian2/codegen/runtime/cython_rt/templates/synapses_create_generator.pyx
index eb9200616..ce313e090 100644
--- a/brian2/codegen/runtime/cython_rt/templates/synapses_create_generator.pyx
+++ b/brian2/codegen/runtime/cython_rt/templates/synapses_create_generator.pyx
@@ -68,8 +68,11 @@ cdef void _flush_buffer(buf, dynarr, int buf_len):
     {{scalar_code['update']|autoindent}}
 
     for _{{outer_index}} in range({{outer_index_size}}):
+        {% if outer_contiguous %}
         _raw{{outer_index_array}} = _{{outer_index}} + {{outer_index_offset}}
-
+        {% else %}
+        _raw{{outer_index_array}} = {{get_array_name(variables[outer_sub_idx])}}[_{{outer_index}}]
+        {% endif %}
         {% if not result_index_condition %}
         {{vector_code['create_cond']|autoindent}}
         if not _cond:
@@ -162,8 +165,11 @@ cdef void _flush_buffer(buf, dynarr, int buf_len):
         {% endif %}
 
             {{vector_code['generator_expr']|autoindent}}
+            {% if result_contiguous %}
             _raw{{result_index_array}} = _{{result_index}} + {{result_index_offset}}
-
+            {% else %}
+            _raw{{result_index_array}} = {{get_array_name(variables[result_sub_idx])}}[_{{result_index}}]
+            {% endif %}
             {% if result_index_condition %}
             {% if result_index_used %}
             {# The condition could index outside of array range #}
diff --git a/brian2/codegen/runtime/numpy_rt/templates/synapses_create_generator.py_ b/brian2/codegen/runtime/numpy_rt/templates/synapses_create_generator.py_
index 0e201765c..86d1a1cd6 100644
--- a/brian2/codegen/runtime/numpy_rt/templates/synapses_create_generator.py_
+++ b/brian2/codegen/runtime/numpy_rt/templates/synapses_create_generator.py_
@@ -106,7 +106,11 @@ _vectorisation_idx = 1
 "k" is called the inner variable #}
 
 for _{{outer_index}} in range({{outer_index_size}}):
+    {% if outer_contiguous %}
     _raw{{outer_index_array}} = _{{outer_index}} + {{outer_index_offset}}
+    {% else %}
+    _raw{{outer_index_array}} = {{get_array_name(variables[outer_sub_idx])}}[_{{outer_index}}]
+    {% endif %}
     {% if not result_index_condition %}
     {{vector_code['create_cond']|autoindent}}
     if not _cond:
@@ -126,7 +130,11 @@ for _{{outer_index}} in range({{outer_index_size}}):
     _vectorisation_idx = {{inner_variable}}
     {{vector_code['generator_expr']|autoindent}}
     _vectorisation_idx = _{{result_index}}
-    _raw{{result_index_array}} = _{{result_index}} + {{result_index_offset}}
+    {% if result_contiguous %}
+    _raw{{result_index_array}} = _{{result_index}} + {{result_index_offset}};
+    {% else %}
+    _raw{{result_index_array}} = {{get_array_name(variables[result_sub_idx])}}[_{{result_index}}]
+    {% endif %}
     {% if result_index_condition %}
     {% if result_index_used %}
     {# The condition could index outside of array range #}
@@ -184,7 +192,11 @@ for _{{outer_index}} in range({{outer_index_size}}):
     {% endif %}
 
     _vectorisation_idx = _{{result_index}}
+    {% if result_contiguous %}
     _raw{{result_index_array}} = _{{result_index}} + {{result_index_offset}}
+    {% else %}
+    _raw{{result_index_array}} = {{get_array_name(variables[result_sub_idx])}}[_{{result_index}}]
+    {% endif %}
     {{vector_code['update']|autoindent}}
 
     if not _numpy.isscalar(_n):
diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 88b211edc..9e57c1aaf 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -1331,10 +1331,14 @@ def _create_variables(self, equations, user_dtype=None):
             self.variables.add_reference("_target_offset", self.target, "_offset")
         else:
             self.variables.add_constant("_target_offset", value=0)
+        if "_sub_idx" in self.source.variables:
+            self.variables.add_reference("_source_sub_idx", self.source, "_sub_idx")
         if "_offset" in self.source.variables:
             self.variables.add_reference("_source_offset", self.source, "_offset")
         else:
             self.variables.add_constant("_source_offset", value=0)
+        if "_sub_idx" in self.target.variables:
+            self.variables.add_reference("_target_sub_idx", self.target, "_sub_idx")
         # To cope with connections to/from other synapses, N_incoming/N_outgoing
         # will be resized when synapses are created
         self.variables.add_dynamic_array(
@@ -1823,13 +1827,19 @@ def _resize(self, number):
         self.variables["N"].set_value(number)
 
     def _update_synapse_numbers(self, old_num_synapses):
-        source_offset = self.variables["_source_offset"].get_value()
-        target_offset = self.variables["_target_offset"].get_value()
-        # This resizing is only necessary if we are connecting to/from synapses
-        post_with_offset = self.variables["N_post"].item() + target_offset
-        pre_with_offset = self.variables["N_pre"].item() + source_offset
+        if "_source_sub_idx" in self.variables:
+            pre_with_offset = self.variables["_source_sub_idx"].get_value()[-1] + 1
+        else:
+            source_offset = self.variables["_source_offset"].get_value()
+            pre_with_offset = self.variables["N_pre"].item() + source_offset
+        if "_target_sub_idx" in self.variables:
+            post_with_offset = self.variables["_target_sub_idx"].get_value()[-1] + 1
+        else:
+            target_offset = self.variables["_target_offset"].get_value()
+            post_with_offset = self.variables["N_post"].item() + target_offset
         self.variables["N_incoming"].resize(post_with_offset)
         self.variables["N_outgoing"].resize(pre_with_offset)
+
         N_outgoing = self.variables["N_outgoing"].get_value()
         N_incoming = self.variables["N_incoming"].get_value()
         synaptic_pre = self.variables["_synaptic_pre"].get_value()
@@ -1945,9 +1955,6 @@ def _add_synapses_from_arrays(self, sources, targets, n, p, namespace=None):
             variables.add_reference("_source_offset", self.source, "_offset")
             abstract_code += "_real_sources = sources + _source_offset\n"
         elif not getattr(self.source, "contiguous", True):
-            variables.add_reference(
-                "_source_sub_idx", self.source, "_sub_idx", index="sources"
-            )
             abstract_code += "_real_sources = _source_sub_idx\n"
         else:
             abstract_code += "_real_sources = sources\n"
@@ -1955,9 +1962,6 @@ def _add_synapses_from_arrays(self, sources, targets, n, p, namespace=None):
             variables.add_reference("_target_offset", self.target, "_offset")
             abstract_code += "_real_targets = targets + _target_offset\n"
         elif not getattr(self.target, "contiguous", True):
-            variables.add_reference(
-                "_target_sub_idx", self.target, "_sub_idx", index="targets"
-            )
             abstract_code += "_real_targets = _target_sub_idx\n"
         else:
             abstract_code += "_real_targets = targets"
@@ -2165,7 +2169,6 @@ def _add_synapses_generator(
             variables.add_constant("_source_offset", value=0)
 
         if not getattr(self.source, "contiguous", True):
-            variables.add_reference("_source_sub_idx", self.source, "_sub_idx")
             needed_variables.append("_source_sub_idx")
 
         if "_offset" in self.target.variables:
@@ -2174,7 +2177,6 @@ def _add_synapses_generator(
             variables.add_constant("_target_offset", value=0)
 
         if not getattr(self.target, "contiguous", True):
-            variables.add_reference("_target_sub_idx", self.target, "_sub_idx")
             needed_variables.append("_target_sub_idx")
 
         variables.add_auxiliary_variable("_raw_pre_idx", dtype=np.int32)

From ce5bb3b7f5d21419c57fd4d44f38ebeafb53decd Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Wed, 15 Nov 2023 00:15:30 +0100
Subject: [PATCH 49/52] Non-contiguous subgroups for synapses: summed variables

---
 .../cython_rt/templates/summed_variable.pyx   |  4 +++
 .../numpy_rt/templates/summed_variable.py_    |  4 +++
 brian2/core/network.py                        | 34 ++++++++++++-------
 .../templates/summed_variable.cpp             |  5 ++-
 brian2/synapses/synapses.py                   | 11 ++++--
 brian2/tests/test_subgroup.py                 | 22 ++++++++----
 6 files changed, 58 insertions(+), 22 deletions(-)

diff --git a/brian2/codegen/runtime/cython_rt/templates/summed_variable.pyx b/brian2/codegen/runtime/cython_rt/templates/summed_variable.pyx
index 4d2b06827..0911b14a1 100644
--- a/brian2/codegen/runtime/cython_rt/templates/summed_variable.pyx
+++ b/brian2/codegen/runtime/cython_rt/templates/summed_variable.pyx
@@ -8,7 +8,11 @@
 
     # Set all the target variable values to zero
     for _target_idx in range({{_target_size_name}}):
+        {% if _target_contiguous %}
         {{_target_var_array}}[_target_idx + {{_target_start}}] = 0
+        {% else %}
+        {{_target_var_array}}[{{_target_indices}}[_target_idx]] = 0
+        {% endif %}
 
     # scalar code
     _vectorisation_idx = 1
diff --git a/brian2/codegen/runtime/numpy_rt/templates/summed_variable.py_ b/brian2/codegen/runtime/numpy_rt/templates/summed_variable.py_
index b925f66f3..734497a39 100644
--- a/brian2/codegen/runtime/numpy_rt/templates/summed_variable.py_
+++ b/brian2/codegen/runtime/numpy_rt/templates/summed_variable.py_
@@ -17,6 +17,9 @@ _vectorisation_idx = LazyArange(N)
 # We write to the array, using the name provided as a keyword argument to the
 # template
 # Note that for subgroups, we do not want to overwrite the full array
+{% if not _target_contiguous %}
+{{_target_var_array}}[{{_target_indices}}] = _numpy.broadcast_to(_synaptic_var, (N, ))
+{% else %}
 {% if _target_start > 0 %}
 _indices = {{_index_array}} - {{_target_start}}
 {% else %}
@@ -32,4 +35,5 @@ _length = _target_stop - {{_target_start}}
 {{_target_var_array}}[{{_target_start}}:_target_stop] = _numpy.bincount(_indices,
                                  minlength=_length,
                                  weights=_numpy.broadcast_to(_synaptic_var, (N, )))
+{% endif %}
 {% endblock %}
diff --git a/brian2/core/network.py b/brian2/core/network.py
index e54073fc6..c2af08578 100644
--- a/brian2/core/network.py
+++ b/brian2/core/network.py
@@ -15,6 +15,8 @@
 from collections import Counter, defaultdict, namedtuple
 from collections.abc import Mapping, Sequence
 
+import numpy as np
+
 from brian2.core.base import BrianObject, BrianObjectException
 from brian2.core.clocks import Clock, defaultclock
 from brian2.core.names import Nameable
@@ -312,18 +314,26 @@ def _check_multiple_summed_updaters(objects):
                         "the target group instead."
                     )
                     raise NotImplementedError(msg)
-                elif (
-                    obj.target.start < other_target.stop
-                    and other_target.start < obj.target.stop
-                ):
-                    # Overlapping subgroups
-                    msg = (
-                        "Multiple 'summed variables' target the "
-                        f"variable '{obj.target_var.name}' in overlapping "
-                        f"groups '{other_target.name}' and '{obj.target.name}'. "
-                        "Use separate variables in the target groups instead."
-                    )
-                    raise NotImplementedError(msg)
+                else:
+                    if getattr(obj.target, "contiguous", True):
+                        target_indices = np.arange(obj.target.start, obj.target.stop)
+                    else:
+                        target_indices = obj.target.indices[:]
+                    if getattr(other_target, "contiguous", True):
+                        other_indices = np.arange(other_target.start, other_target.stop)
+                    else:
+                        other_indices = other_target.indices[:]
+                    if np.intersect1d(
+                        target_indices, other_indices, assume_unique=True
+                    ).size:
+                        # Overlapping subgroups
+                        msg = (
+                            "Multiple 'summed variables' target the "
+                            f"variable '{obj.target_var.name}' in overlapping "
+                            f"groups '{other_target.name}' and '{obj.target.name}'. "
+                            "Use separate variables in the target groups instead."
+                        )
+                        raise NotImplementedError(msg)
             summed_targets[obj.target_var] = obj.target
 
 
diff --git a/brian2/devices/cpp_standalone/templates/summed_variable.cpp b/brian2/devices/cpp_standalone/templates/summed_variable.cpp
index 9967d0fda..1a9f4d273 100644
--- a/brian2/devices/cpp_standalone/templates/summed_variable.cpp
+++ b/brian2/devices/cpp_standalone/templates/summed_variable.cpp
@@ -7,12 +7,15 @@
     //// MAIN CODE ////////////
     {# This enables summed variables for connections to a synapse #}
     const int _target_size = {{constant_or_scalar(_target_size_name, variables[_target_size_name])}};
-
     // Set all the target variable values to zero
     {{ openmp_pragma('parallel-static') }}
     for (int _target_idx=0; _target_idx<_target_size; _target_idx++)
     {
+        {% if _target_contiguous %}
         {{_target_var_array}}[_target_idx + {{_target_start}}] = 0;
+        {% else %}
+        {{_target_var_array}}[{{_target_indices}}[_target_idx]] = 0;
+        {% endif %}
     }
 
     // scalar code
diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 9e57c1aaf..73ac587c1 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -122,7 +122,7 @@ class SummedVariableUpdater(CodeRunner):
     def __init__(
         self, expression, target_varname, synapses, target, target_size_name, index_var
     ):
-        # Handling sumped variables using the standard mechanisms is not
+        # Handling summed variables using the standard mechanisms is not
         # possible, we therefore also directly give the names of the arrays
         # to the template.
 
@@ -139,14 +139,21 @@ def __init__(
             "_index_var": synapses.variables[index_var],
             "_target_start": getattr(target, "start", 0),
             "_target_stop": getattr(target, "stop", -1),
+            "_target_contiguous": True,
         }
+        needed_variables = [target_varname, target_size_name, index_var]
+        self.variables = Variables(synapses)
+        if not getattr(target, "contiguous", True):
+            self.variables.add_reference("_target_indices", target, "_sub_idx")
+            needed_variables.append("_target_indices")
+            template_kwds["_target_contiguous"] = False
 
         CodeRunner.__init__(
             self,
             group=synapses,
             template="summed_variable",
             code=code,
-            needed_variables=[target_varname, target_size_name, index_var],
+            needed_variables=needed_variables,
             # We want to update the summed variable before
             # the target group gets updated
             clock=target.clock,
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 69e125474..6ca472912 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -644,20 +644,28 @@ def test_synapses_access_subgroups_problematic():
 @pytest.mark.standalone_compatible
 def test_subgroup_summed_variable():
     # Check in particular that only neurons targeted are reset to 0 (see github issue #925)
-    source = NeuronGroup(1, "")
-    target = NeuronGroup(5, "Iin : 1")
+    source = NeuronGroup(1, "x : 1")
+    target = NeuronGroup(
+        7,
+        """Iin : 1
+                                          x : 1""",
+    )
+    source.x = 5
     target.Iin = 10
+    target.x = "i"
     target1 = target[1:2]
-    target2 = target[3:]
-
-    syn1 = Synapses(source, target1, "Iin_post = 5 : 1  (summed)")
+    target2 = target[3:5]
+    target3 = target[[0, 6]]
+    syn1 = Synapses(source, target1, "Iin_post = x_pre + x_post : 1  (summed)")
     syn1.connect(True)
-    syn2 = Synapses(source, target2, "Iin_post = 1 : 1  (summed)")
+    syn2 = Synapses(source, target2, "Iin_post = x_pre + x_post : 1  (summed)")
     syn2.connect(True)
+    syn3 = Synapses(source, target3, "Iin_post = x_pre + x_post : 1  (summed)")
+    syn3.connect(True)
 
     run(2 * defaultclock.dt)
 
-    assert_array_equal(target.Iin, [10, 5, 10, 1, 1])
+    assert_array_equal(target.Iin, [5, 6, 10, 8, 9, 10, 11])
 
 
 def test_subexpression_references():

From 16e2ffc44118089201a7d52cf0b42ac10530b027 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Thu, 16 Nov 2023 12:02:02 +0100
Subject: [PATCH 50/52] Test summed variable error for overlapping subgroups

---
 brian2/tests/test_subgroup.py | 69 +++++++++++++++++++++++++++--------
 1 file changed, 54 insertions(+), 15 deletions(-)

diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 6ca472912..e67683559 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -648,7 +648,7 @@ def test_subgroup_summed_variable():
     target = NeuronGroup(
         7,
         """Iin : 1
-                                          x : 1""",
+                  x : 1""",
     )
     source.x = 5
     target.Iin = 10
@@ -668,6 +668,45 @@ def test_subgroup_summed_variable():
     assert_array_equal(target.Iin, [5, 6, 10, 8, 9, 10, 11])
 
 
+@pytest.mark.codegen_independent
+def test_subgroup_summed_variable_overlap():
+    # Check that overlapping subgroups raise an error
+    source = NeuronGroup(1, "")
+    target = NeuronGroup(10, "Iin : 1")
+    target1 = target[1:3]
+    target2 = target[2:5]
+    target3 = target[[1, 6]]
+    target4 = target[[4, 6]]
+
+    syn1 = Synapses(source, target1, "Iin_post = 1 : 1  (summed)")
+    syn1.connect(True)
+
+    syn2 = Synapses(source, target2, "Iin_post = 2 : 1  (summed)")
+    syn2.connect(True)
+
+    syn3 = Synapses(source, target3, "Iin_post = 3 : 1  (summed)")
+    syn3.connect(True)
+
+    syn4 = Synapses(source, target4, "Iin_post = 4 : 1  (summed)")
+    syn4.connect(True)
+
+    net1 = Network(source, target, syn1, syn2)  # overlap between contiguous subgroups
+    with pytest.raises(NotImplementedError):
+        net1.run(0 * ms)
+
+    net2 = Network(
+        source, target, syn1, syn3
+    )  # overlap between contiguous and non-contiguous subgroups
+    with pytest.raises(NotImplementedError):
+        net2.run(0 * ms)
+
+    net3 = Network(
+        source, target, syn3, syn4
+    )  # overlap between non-contiguous subgroups
+    with pytest.raises(NotImplementedError):
+        net3.run(0 * ms)
+
+
 def test_subexpression_references():
     """
     Assure that subexpressions in targeted groups are handled correctly.
@@ -675,7 +714,7 @@ def test_subexpression_references():
     G = NeuronGroup(
         10,
         """v : 1
-                           v2 = 2*v : 1""",
+                  v2 = 2*v : 1""",
     )
     G.v = np.arange(10)
     SG1 = G[:5]
@@ -685,8 +724,8 @@ def test_subexpression_references():
         SG1,
         SG2,
         """w : 1
-                          u = v2_post + 1 : 1
-                          x = v2_pre + 1 : 1""",
+           u = v2_post + 1 : 1
+           x = v2_pre + 1 : 1""",
     )
     S1.connect("i==(5-1-j)")
     assert_equal(S1.i[:], np.arange(5))
@@ -698,8 +737,8 @@ def test_subexpression_references():
         G,
         SG2,
         """w : 1
-                             u = v2_post + 1 : 1
-                             x = v2_pre + 1 : 1""",
+           u = v2_post + 1 : 1
+           x = v2_pre + 1 : 1""",
     )
     S2.connect("i==(5-1-j)")
     assert_equal(S2.i[:], np.arange(5))
@@ -711,8 +750,8 @@ def test_subexpression_references():
         SG1,
         G,
         """w : 1
-                             u = v2_post + 1 : 1
-                             x = v2_pre + 1 : 1""",
+           u = v2_post + 1 : 1
+           x = v2_pre + 1 : 1""",
     )
     S3.connect("i==(10-1-j)")
     assert_equal(S3.i[:], np.arange(5))
@@ -729,7 +768,7 @@ def test_subexpression_no_references():
     G = NeuronGroup(
         10,
         """v : 1
-                           v2 = 2*v : 1""",
+                  v2 = 2*v : 1""",
     )
     G.v = np.arange(10)
 
@@ -739,8 +778,8 @@ def test_subexpression_no_references():
         G[:5],
         G[5:],
         """w : 1
-                          u = v2_post + 1 : 1
-                          x = v2_pre + 1 : 1""",
+           u = v2_post + 1 : 1
+           x = v2_pre + 1 : 1""",
     )
     S1.connect("i==(5-1-j)")
     assert_equal(S1.i[:], np.arange(5))
@@ -752,8 +791,8 @@ def test_subexpression_no_references():
         G,
         G[5:],
         """w : 1
-                             u = v2_post + 1 : 1
-                             x = v2_pre + 1 : 1""",
+           u = v2_post + 1 : 1
+           x = v2_pre + 1 : 1""",
     )
     S2.connect("i==(5-1-j)")
     assert_equal(S2.i[:], np.arange(5))
@@ -765,8 +804,8 @@ def test_subexpression_no_references():
         G[:5],
         G,
         """w : 1
-                             u = v2_post + 1 : 1
-                             x = v2_pre + 1 : 1""",
+           u = v2_post + 1 : 1
+           x = v2_pre + 1 : 1""",
     )
     S3.connect("i==(10-1-j)")
     assert_equal(S3.i[:], np.arange(5))

From f1ea1b52d7b5435edbf0d4d70d1a7015915b41ac Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Fri, 13 Sep 2024 13:28:20 +0200
Subject: [PATCH 51/52] minor fixes

---
 brian2/codegen/runtime/numpy_rt/numpy_rt.py | 2 +-
 brian2/tests/test_subgroup.py               | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/brian2/codegen/runtime/numpy_rt/numpy_rt.py b/brian2/codegen/runtime/numpy_rt/numpy_rt.py
index d56e8b855..e94211ef8 100644
--- a/brian2/codegen/runtime/numpy_rt/numpy_rt.py
+++ b/brian2/codegen/runtime/numpy_rt/numpy_rt.py
@@ -127,7 +127,7 @@ def __iter__(self):
             return iter(self.indices)
 
     # Allow conversion to a proper array with np.array(...)
-    def __array__(self, dtype=None, copy=None):
+    def __array__(self, dtype=np.int32, copy=None):
         if copy is False:
             raise ValueError("LazyArange does not support copy=False")
         if self.indices is None:
diff --git a/brian2/tests/test_subgroup.py b/brian2/tests/test_subgroup.py
index 045dab041..a141fa916 100644
--- a/brian2/tests/test_subgroup.py
+++ b/brian2/tests/test_subgroup.py
@@ -951,7 +951,6 @@ def test_run_regularly():
 
 @pytest.mark.standalone_compatible
 def test_spike_monitor():
-    set_device("cpp_standalone", directory="/tmp/testsubgroup")
     G = NeuronGroup(10, "v:1", threshold="v>1", reset="v=0")
     G.v[[0, 2, 5]] = 1.1
     SG = G[3:]

From c21ec60b4839526ad628102040f42c04f472add3 Mon Sep 17 00:00:00 2001
From: Marcel Stimberg <marcel.stimberg@sorbonne-universite.fr>
Date: Fri, 13 Sep 2024 15:59:26 +0200
Subject: [PATCH 52/52] Re-add write sizes of dynamic arrays to disk

Changes from commit fe0add72821 got lost in merge
---
 .../cpp_standalone/templates/objects.cpp      | 25 +++++++++++++++----
 1 file changed, 20 insertions(+), 5 deletions(-)

diff --git a/brian2/devices/cpp_standalone/templates/objects.cpp b/brian2/devices/cpp_standalone/templates/objects.cpp
index fab6bf3a0..411a29500 100644
--- a/brian2/devices/cpp_standalone/templates/objects.cpp
+++ b/brian2/devices/cpp_standalone/templates/objects.cpp
@@ -277,15 +277,21 @@ void _write_arrays()
     outfile_{{varname}}.open(results_dir + "{{get_array_filename(var)}}", ios::binary | ios::out);
     if(outfile_{{varname}}.is_open())
     {
-        if (! {{varname}}.empty() )
-        {
-            outfile_{{varname}}.write(reinterpret_cast<char*>(&{{varname}}[0]), {{varname}}.size()*sizeof({{varname}}[0]));
-            outfile_{{varname}}.close();
-        }
+        outfile_{{varname}}.write(reinterpret_cast<char*>(&{{varname}}[0]), {{varname}}.size()*sizeof({{varname}}[0]));
+        outfile_{{varname}}.close();
     } else
     {
         std::cout << "Error writing output file for {{varname}}." << endl;
     }
+    outfile_{{varname}}.open("{{get_array_filename(var) | replace('\\', '\\\\')}}_size", ios::out);
+        if (outfile_{{varname}}.is_open())
+        {
+            outfile_{{varname}} << {{varname}}.size();
+            outfile_{{varname}}.close();
+        } else
+        {
+                std::cout << "Error writing size file for {{varname}}." << endl;
+        }
     {% endfor %}
 
     {% for var, varname in dynamic_array_2d_specs | dictsort(by='value') %}
@@ -305,6 +311,15 @@ void _write_arrays()
     {
         std::cout << "Error writing output file for {{varname}}." << endl;
     }
+    outfile_{{varname}}.open("{{get_array_filename(var) | replace('\\', '\\\\')}}_size", ios::out);
+        if (outfile_{{varname}}.is_open()) {
+            outfile_{{varname}} << {{varname}}.n << " " << {{varname}}.m;
+            outfile_{{varname}}.close();
+        } else
+        {
+                std::cout << "Error writing size file for {{varname}}." << endl;
+        }
+
     {% endfor %}
     {% if profiled_codeobjects is defined and profiled_codeobjects %}
     // Write profiling info to disk