Sampling from ranges with step size ≠ 1 in synapse generator not working correclty (except in numpy)
While working on #1280 I extended our test suite and realized that our current implementation for `sample(..., p=...)` does not give the correct results for ranges with step size ≠ 1 (e.g. `sample(0, N_post, 2, p=0.5)` to only sample from even numbers) for Cython, and neither for ranges going backwards such as `sample(N_post-1, -1, -3, p=0.5)` both for Cython and C++ standalone. Apparently both features are rarely used (otherwise we'd have had complaints about this), but since we provide them, they should of course work. Fixing this is not very difficult (I will do it as part of #1280), but we could also decide that we do not want to support backward ranges at all. On the other hand, checking and raising the error needs more code than supporting it...
diff --git a/brian2/codegen/runtime/cython_rt/templates/common_group.pyx b/brian2/codegen/runtime/cython_rt/templates/common_group.pyx
index fe047cc5..435fd7c2 100644
--- a/brian2/codegen/runtime/cython_rt/templates/common_group.pyx
+++ b/brian2/codegen/runtime/cython_rt/templates/common_group.pyx
@@ -17,6 +17,9 @@ from libc.stdlib cimport abs  # For integers
 from libc.math cimport abs  # For floating point values
 from libc.limits cimport INT_MIN, INT_MAX
 from libcpp cimport bool
+from libcpp.set cimport set
+from cython.operator cimport dereference as _deref, preincrement as _preinc
+cimport cython as _cython
 
 _numpy.import_array()
 cdef extern from "numpy/ndarraytypes.h":
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 e6dd74f0..6e5690e8 100644
--- a/brian2/codegen/runtime/cython_rt/templates/synapses_create_generator.pyx
+++ b/brian2/codegen/runtime/cython_rt/templates/synapses_create_generator.pyx
@@ -38,10 +38,22 @@ cdef void _flush_buffer(buf, dynarr, int buf_len):
     cdef size_t newsize
 
     # The following variables are only used for probabilistic connections
+    {% if iterator_func=='sample' %}
+    cdef int _iter_sign
+    {% if iterator_kwds['sample_size'] == 'fixed' %}
+    cdef bool _selection_algo
+    cdef set[int] _selected_set = set[int]()
+    cdef set[int].iterator _selected_it
+    cdef int _n_selected
+    cdef int _n_dealt_with
+    cdef int _n_total
+    cdef double _U
+    {% else %}
     cdef bool _jump_algo
     cdef double _log1p, _pconst
     cdef size_t _jump
-
+    {% endif %}
+    {% endif %}
     # scalar code
     _vectorisation_idx = 1
     {{scalar_code['setup_iterator']|autoindent}}
@@ -61,25 +73,86 @@ cdef void _flush_buffer(buf, dynarr, int buf_len):
         {% if iterator_func=='range' %}
         for {{iteration_variable}} in range(_iter_low, _iter_high, _iter_step):
         {% elif iterator_func=='sample' %}
+        {% if iterator_kwds['sample_size'] == 'fixed' %}
+        # Note that the following code is written in a slightly convoluted way,
+        # but we have to plug it together with the following code that checks
+        # for the fulfillment of the condition.
+        _n_selected = 0
+        _n_dealt_with = 0
+        with _cython.cdivision(True):
+            if _iter_step > 0:
+                _n_total = (_iter_high - _iter_low - 1) // _iter_step + 1
+            else:
+                _n_total = (_iter_low - _iter_high - 1) // -_iter_step + 1
+            # Value determined by benchmarking, see github PR #1280
+            _selection_algo = 1.0*_iter_size / _n_total > 0.06
+        if _iter_size > _n_total:
+            {% if skip_if_invalid %}
+            _iter_size = _n_total
+            {% else %}
+            raise IndexError('Requested sample size %d is bigger than the '
+                             'population size %d.' % (_iter_size, _n_total))
+            {% endif %}
+        elif _iter_size < 0:
+            {% if skip_if_invalid %}
+            continue
+            {% else %}
+            raise IndexError('Requested sample size %d is negative.' % _iter_size)
+            {% endif %}
+        if _selection_algo:
+            {{iteration_variable}} = _iter_low - _iter_step
+        else:
+            # For the tracking algorithm, we have to first create all values
+            # to make sure they will be iterated in sorted order
+            _selected_set.clear()
+            while _n_selected < _iter_size:
+                _r = <int> (_rand(_vectorisation_idx) * _n_total)
+                while not _selected_set.insert(_r).second:  # .second will be False if duplicate
+                    _r = <int> (_rand(_vectorisation_idx) * _n_total)
+                _n_selected += 1
+            _n_selected = 0
+            _selected_it = _selected_set.begin()
+
+        while _n_selected < _iter_size:
+            if _selection_algo:
+                {{iteration_variable}} += _iter_step
+                # Selection sampling technique
+                # See section 3.4.2 of Donald E. Knuth, AOCP, Vol 2,
+                # Seminumerical Algorithms
+                _n_dealt_with += 1
+                _U = _rand(_vectorisation_idx)
+                if (_n_total - _n_dealt_with) * _U >= _iter_size - _n_selected:
+                    continue
+            else:
+                {{iteration_variable}} = _iter_low + _deref(_selected_it)*_iter_step
+                _preinc(_selected_it)
+            _n_selected += 1
+
+        {% else %}
         if _iter_p==0:
             continue
+        if _iter_step < 0:
+            _iter_sign = -1
+        else:
+            _iter_sign = 1
         _jump_algo = _iter_p<0.25
         if _jump_algo:
             _log1p = log(1-_iter_p)
         else:
             _log1p = 1.0 # will be ignored
         _pconst = 1.0/_log1p
-        {{iteration_variable}} = _iter_low-1
-        while {{iteration_variable}}+1<_iter_high:
-            {{iteration_variable}} += 1
+        {{iteration_variable}} = _iter_low-_iter_step
+        while _iter_sign*({{iteration_variable}} + _iter_step) < _iter_sign*_iter_high:
+            {{iteration_variable}} += _iter_step
             if _jump_algo:
-                _jump = int(floor(log(_rand(_vectorisation_idx))*_pconst))*_iter_step
+                _jump = <int>(log(_rand(_vectorisation_idx))*_pconst)*_iter_step
                 {{iteration_variable}} += _jump
-                if {{iteration_variable}}>=_iter_high:
+                if _iter_sign*{{iteration_variable}} >= _iter_sign*_iter_high:
                     break
             else:
                 if _rand(_vectorisation_idx)>=_iter_p:
                     continue
+        {% endif %}
         {% endif %}
 
             {{vector_code['create_j']|autoindent}}
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 6c7a6bce..b922f655 100644
--- a/brian2/codegen/runtime/numpy_rt/templates/synapses_create_generator.py_
+++ b/brian2/codegen/runtime/numpy_rt/templates/synapses_create_generator.py_
@@ -21,7 +21,7 @@ except:
 
     from numpy.random import binomial as _binomial
     import bisect as _bisect
-    def _sample_without_replacement(low, high, step, p):
+    def _sample_without_replacement(low, high, step, p=None, size=None):
         if p == 0:
             return _numpy.array([], dtype=_numpy.int32)
         if p == 1:
@@ -30,8 +30,30 @@ except:
             n = (high-low-1)//step+1
         else:
             n = (low-high-1)//-step+1
-        if n <= 0:
+        if n <= 0 or (size is not None and size == 0):
             return _numpy.array([], dtype=_numpy.int32)
+        elif size is not None and size < 0:
+            {% if skip_if_invalid %}
+            return _numpy.array([], dtype=_numpy.int32)
+            {% else %}
+            raise IndexError('Requested sample size %d is negative.' % size)
+            {% endif %}
+        elif size is not None and size > n:
+            {% if skip_if_invalid %}
+            size = n
+            {% else %}
+            raise IndexError('Requested sample size %d is bigger than the '
+                             'population size %d.' % (size, n))
+            {% endif %}
+        if n == size:
+            return _numpy.arange(low, high, step)
+        if p is None:
+            if step == 1:
+                _choice = _numpy.random.choice(n,  size=size, replace=False) + low
+            else:
+                _chose_from = _numpy.arange(low, high, step)
+                _choice = _numpy.random.choice(_chose_from, size=size, replace=False)
+            return _numpy.sort(_choice)
         if n <= 1000 or p > 0.25: # based on rough profiling
             samples, = (_np_rand(n)<p).nonzero()
         else:
@@ -89,8 +111,12 @@ for _i in range(N_pre):
     {% if iterator_func=='range' %}
     {{iteration_variable}} = _numpy.arange(_iter_low, _iter_high, _iter_step)
     {% elif iterator_func=='sample' %}
+    {% if iterator_kwds['sample_size'] == 'fixed' %}
+    {{iteration_variable}} = _sample_without_replacement(_iter_low, _iter_high, _iter_step, None, size=_iter_size)
+    {% else %}
     {{iteration_variable}} = _sample_without_replacement(_iter_low, _iter_high, _iter_step, _iter_p)
     {% endif %}
+    {% endif %}
 
     _vectorisation_idx = {{iteration_variable}}
     {{vector_code['create_j']|autoindent}}
diff --git a/brian2/core/base.py b/brian2/core/base.py
index 35d7e863..9ee249c4 100644
--- a/brian2/core/base.py
+++ b/brian2/core/base.py
@@ -1,6 +1,7 @@
 '''
 All Brian objects should derive from `BrianObject`.
 '''
+import functools
 import weakref
 import traceback
 import os
@@ -277,8 +278,8 @@ def device_override_decorated_function(*args, **kwds):
             else:
                 return func(*args, **kwds)
 
-        device_override_decorated_function.__doc__ = func.__doc__
         device_override_decorated_function.original_function = func
+        functools.update_wrapper(device_override_decorated_function, func)
 
         return device_override_decorated_function
 
diff --git a/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp b/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp
index 9c335fcf..047208a4 100644
--- a/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp
+++ b/brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp
@@ -53,6 +53,7 @@
         long _uiter_high;
         long _uiter_step;
         {% if iterator_func=='sample' %}
+        long _uiter_size;
         double _uiter_p;
         {% endif %}
         {
@@ -61,13 +62,87 @@
             _uiter_high = _iter_high;
             _uiter_step = _iter_step;
             {% if iterator_func=='sample' %}
+            {% if iterator_kwds['sample_size'] == 'fixed' %}
+            _uiter_size = _iter_size;
+            {% else %}
             _uiter_p = _iter_p;
             {% endif %}
+            {% endif %}
         }
         {% if iterator_func=='range' %}
         for(long {{iteration_variable}}=_uiter_low; {{iteration_variable}}<_uiter_high; {{iteration_variable}}+=_uiter_step)
         {
         {% elif iterator_func=='sample' %}
+        const int _iter_sign = _uiter_step > 0 ? 1 : -1;
+        {% if iterator_kwds['sample_size'] == 'fixed' %}
+        std::set<int> _selected_set = std::set<int>();
+        std::set<int>::iterator _selected_it;
+        int _n_selected = 0;
+        int _n_dealt_with = 0;
+        int _n_total;
+        if (_uiter_step > 0)
+            _n_total = (_uiter_high - _uiter_low - 1) / _uiter_step + 1;
+        else
+            _n_total = (_uiter_low - _uiter_high - 1) / -_uiter_step + 1;
+        // Value determined by benchmarking, see github PR #1280
+        const bool _selection_algo = 1.0*_uiter_size / _n_total > 0.06;
+        if (_uiter_size > _n_total)
+        {
+            {% if skip_if_invalid %}
+            _uiter_size = _n_total;
+            {% else %}
+            cout << "Error: Requested sample size " << _uiter_size << " is bigger than the " <<
+                    "population size " << _n_total << "." << endl;
+            exit(1);
+            {% endif %}
+        } else if (_uiter_size < 0)
+        {
+            {% if skip_if_invalid %}
+            continue;
+            {% else %}
+            cout << "Error: Requested sample size " << _uiter_size << " is negative." << endl;
+            exit(1);
+            {% endif %}
+        } else if (_uiter_size == 0)
+            continue;
+        long {{iteration_variable}};
+
+        if (_selection_algo)
+        {
+            {{iteration_variable}} = _uiter_low - _uiter_step;
+        } else
+        {
+            // For the tracking algorithm, we have to first create all values
+            // to make sure they will be iterated in sorted order
+            _selected_set.clear();
+            while (_n_selected < _uiter_size)
+            {
+                int _r = (int)(_rand(_vectorisation_idx) * _n_total);
+                while (! _selected_set.insert(_r).second)
+                    _r = (int)(_rand(_vectorisation_idx) * _n_total);
+                _n_selected++;
+            }
+            _n_selected = 0;
+            _selected_it = _selected_set.begin();
+        }
+        while (_n_selected < _uiter_size)
+        {
+            if (_selection_algo)
+            {
+                // Selection sampling technique
+                // See section 3.4.2 of Donald E. Knuth, AOCP, Vol 2, Seminumerical Algorithms
+                {{iteration_variable}} += _uiter_step;
+                _n_dealt_with++;
+                const double _U = _rand(_vectorisation_idx);
+                if ((_n_total - _n_dealt_with) * _U >= _uiter_size - _n_selected)
+                    continue;
+            } else
+            {
+                {{iteration_variable}} = _uiter_low + (*_selected_it)*_uiter_step;
+                _selected_it++;
+            }
+            _n_selected++;
+        {% else %}
         if(_uiter_p==0) continue;
         const bool _jump_algo = _uiter_p<0.25;
         double _log1p;
@@ -76,17 +151,18 @@
         else
             _log1p = 1.0; // will be ignored
         const double _pconst = 1.0/log(1-_uiter_p);
-        for(long {{iteration_variable}}=_uiter_low; {{iteration_variable}}<_uiter_high; {{iteration_variable}}++)
+        for(long {{iteration_variable}}=_uiter_low; _iter_sign*{{iteration_variable}}<_iter_sign*_uiter_high; {{iteration_variable}} += _uiter_step)
         {
             if(_jump_algo) {
                 const double _r = _rand(_vectorisation_idx);
                 if(_r==0.0) break;
                 const int _jump = floor(log(_r)*_pconst)*_uiter_step;
                 {{iteration_variable}} += _jump;
-                if({{iteration_variable}}>=_uiter_high) continue;
+                if (_iter_sign*{{iteration_variable}} >= _iter_sign * _uiter_high) continue;
             } else {
-                if(_rand(_vectorisation_idx)>=_uiter_p) continue;
+                if (_rand(_vectorisation_idx)>=_uiter_p) continue;
             }
+        {% endif %}
         {% endif %}
             long __j, _j, _pre_idx, __pre_idx;
             {
diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 466224ee..95d25227 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -1347,6 +1347,7 @@ def connect(self, condition=None, i=None, j=None, p=1., n=1,
         >>> S.connect(j='k for k in range(i+1)') # Connect neuron i to all j with 0<=j<=i
         >>> S.connect(j='i+(-1)**k for k in range(2) if i>0 and i<N_pre-1') # connect neuron i to its neighbours if it has both neighbours
         >>> S.connect(j='k for k in sample(N_post, p=i*1.0/(N_pre-1))') # neuron i connects to j with probability i/(N-1)
+        >>> S.connect(j='k for k in sample(N_post, size=i//2)') # Each neuron connects to i//2 other neurons (chosen randomly)
         '''
         # check types
         if condition is not None and not isinstance(condition, (bool,
@@ -1696,10 +1697,6 @@ def _add_synapses_generator(self, j, n, skip_if_invalid=False, namespace=None, l
         template_kwds.update(parsed)
         template_kwds['skip_if_invalid'] = skip_if_invalid
 
-        if (parsed['iterator_func'] == 'sample' and
-                    parsed['iterator_kwds']['sample_size']=='fixed'):
-            raise NotImplementedError("Fixed sample size not implemented yet.")
-
         abstract_code = {'setup_iterator': '',
                          'create_j': '',
                          'create_cond': '',
@@ -1799,7 +1796,6 @@ def _add_synapses_generator(self, j, n, skip_if_invalid=False, namespace=None, l
                       "using generator '%s'") % (self.source.name,
                                                  self.target.name,
                                                  parsed['original_expression']))
-
         codeobj = create_runner_codeobj(self,
                                         abstract_code,
                                         'synapses_create_generator',
diff --git a/brian2/tests/test_synapses.py b/brian2/tests/test_synapses.py
index c36cd5ce..4443d06e 100644
--- a/brian2/tests/test_synapses.py
+++ b/brian2/tests/test_synapses.py
@@ -2438,12 +2438,274 @@ def test_synapse_generator_random():
 
     assert len(S1) == 0
     _compare(S2, np.ones((len(G), len(G2))))
-    assert 0 <= len(S2) <= len(G) * len(G2)
+    assert 0 <= len(S3) <= len(G) * len(G2)
+    assert len(S4) == 7
+    assert_equal(S4.i, np.ones(7)*2)
+    assert_equal(S4.j, np.arange(7))
+
+
+@pytest.mark.standalone_compatible
+def test_synapse_generator_random_positive_steps():
+    # Test generator with sampling from stepped ranges (e.g. all even numbers)
+    G = NeuronGroup(4, '''v: 1
+                          x : integer''', threshold='False')
+    G.x = 'i'
+    G2 = NeuronGroup(7, '''v: 1
+                           y : 1''', threshold='False')
+    G2.y = '1.0*i/N'
+
+    S1 = Synapses(G, G2, 'w:1', 'v+=w')
+    S1.connect(j='k for k in sample(2, N_post, 2, p=0)')
+
+    S2 = Synapses(G, G2, 'w:1', 'v+=w')
+    S2.connect(j='k for k in sample(2, N_post, 2, p=1)')
+
+    # Just make sure using values between 0 and 1 work in principle (note that
+    # 0.25 is the cutoff between the general method and the "jump method", so
+    # we test a value above and below
+    S3 = Synapses(G, G2, 'w:1', 'v+=w')
+    S3.connect(j='k for k in sample(2, N_post, 2, p=0.2)')
+
+    S3b = Synapses(G, G2, 'w:1', 'v+=w')
+    S3b.connect(j='k for k in sample(2, N_post, 2, p=0.3)')
+
+    # Use pre-/post-synaptic variables for "stochastic" connections that are
+    # actually deterministic
+    S4 = Synapses(G, G2, 'w:1', on_pre='v+=w')
+    S4.connect(j='k for k in sample(2, N_post, 2, p=int(x_pre==2)*1.0)')
+
+    with catch_logs() as _:  # Ignore warnings about empty synapses
+        run(0*ms)  # for standalone
+
+    assert len(S1) == 0
+    S2_comp = np.zeros((len(G), len(G2)))
+    S2_comp[:, 2::2] = 1
+    _compare(S2, S2_comp)
+    assert 0 <= len(S3) <= len(G) * 3
+    assert all(S3.j[:] % 2 == 0)
+    assert all(S3.j >= 2)
+    assert 0 <= len(S3b) <= len(G) * 3
+    assert all(S3b.j[:] % 2 == 0)
+    assert all(S3b.j >= 2)
+    assert len(S4) == 3
+    assert_equal(S4.i, np.ones(3)*2)
+    assert_equal(S4.j, np.arange(2, 7, 2))
+
+
+@pytest.mark.standalone_compatible
+def test_synapse_generator_random_negative_steps():
+    # Test generator with sampling from stepped ranges (e.g. all even numbers)
+    # going backwards
+    G = NeuronGroup(4, '''v: 1
+                          x : integer''', threshold='False')
+    G.x = 'i'
+    G2 = NeuronGroup(7, '''v: 1
+                           y : 1''', threshold='False')
+    G2.y = '1.0*i/N'
+
+    S1 = Synapses(G, G2, 'w:1', 'v+=w')
+    S1.connect(j='k for k in sample(N_post-1, 0, -2, p=0)')
+
+    S2 = Synapses(G, G2, 'w:1', 'v+=w')
+    S2.connect(j='k for k in sample(N_post-1, 0, -2, p=1)')
+
+    # Just make sure using values between 0 and 1 work in principle (note that
+    # 0.25 is the cutoff between the general method and the "jump method", so
+    # we test a value above and below
+    S3 = Synapses(G, G2, 'w:1', 'v+=w')
+    S3.connect(j='k for k in sample(N_post-1, 0, -2, p=0.2)')
+
+    S3b = Synapses(G, G2, 'w:1', 'v+=w')
+    S3b.connect(j='k for k in sample(N_post-1, 0, -2, p=0.3)')
+
+    # Use pre-/post-synaptic variables for "stochastic" connections that are
+    # actually deterministic
+    S4 = Synapses(G, G2, 'w:1', on_pre='v+=w')
+    S4.connect(j='k for k in sample(N_post-1, 0, -2, p=int(x_pre==2)*1.0)')
+
+    with catch_logs() as _:  # Ignore warnings about empty synapses
+        run(0*ms)  # for standalone
+
+    assert len(S1) == 0
+    S2_comp = np.zeros((len(G), len(G2)))
+    S2_comp[:, 2::2] = 1
+    _compare(S2, S2_comp)
+    assert 0 <= len(S3) <= len(G) * 3
+    assert all(S3.j[:] % 2 == 0)
+    assert all(S3.j >= 2)
+    assert 0 <= len(S3b) <= len(G) * 3
+    assert all(S3b.j[:] % 2 == 0)
+    assert all(S3b.j >= 2)
+    assert len(S4) == 3
+    assert_array_equal(S4.i, np.ones(3)*2)
+    assert_array_equal(S4.j, [6, 4, 2])
+
+
+@pytest.mark.standalone_compatible
+def test_synapse_generator_fixed_random():
+    # Random samples with fixed size
+    G = NeuronGroup(4, '''v: 1
+                          x : integer''', threshold='False')
+    G.x = 'i'
+    G2 = NeuronGroup(7, '''v: 1
+                           y : 1''', threshold='False')
+    G2.y = '1.0*i/N'
+
+    S1 = Synapses(G, G2, 'w:1', 'v+=w')
+    S1.connect(j='k for k in sample(N_post, size=0)')
+
+    S2 = Synapses(G, G2, 'w:1', 'v+=w')
+    S2.connect(j='k for k in sample(N_post, size=N_post)')
+
+    S3 = Synapses(G, G2, 'w:1', 'v+=w')
+    S3.connect(j='k for k in sample(N_post, size=3)')
+
+    # Use pre-/post-synaptic variables for "stochastic" connections that are
+    # actually deterministic
+    S4 = Synapses(G, G2, 'w:1', on_pre='v+=w')
+    S4.connect(j='k for k in sample(N_post, size=int(x_pre==2)*N_post)')
+
+    with catch_logs() as _:  # Ignore warnings about empty synapses
+        run(0*ms)  # for standalone
+
+    assert len(S1) == 0
+    _compare(S2, np.ones((len(G), len(G2))))
+    # Each neuron should have 3 outgoing connections
+    assert_array_equal(S3.N_outgoing_pre, np.ones(4)*3)
+    # Synapses should be sorted and unique
+    for source_idx in range(4):
+        assert len(set(S3.j[source_idx, :])) == 3
+        assert all(S3.j[source_idx, :] == sorted(S3.j[source_idx, :]))
     assert len(S4) == 7
     assert_equal(S4.i, np.ones(7)*2)
     assert_equal(S4.j, np.arange(7))
 
 
+@pytest.mark.standalone_compatible
+def test_synapse_generator_fixed_random_positive_steps():
+    # Test generator with fixed-size sampling from stepped ranges (e.g. all
+    # even numbers)
+    G = NeuronGroup(4, '''v: 1
+                          x : integer''', threshold='False')
+    G.x = 'i'
+    G2 = NeuronGroup(7, '''v: 1
+                           y : 1''', threshold='False')
+    G2.y = '1.0*i/N'
+
+    S1 = Synapses(G, G2, 'w:1', 'v+=w')
+    S1.connect(j='k for k in sample(2, N_post, 2, size=0)')
+
+    S2 = Synapses(G, G2, 'w:1', 'v+=w')
+    S2.connect(j='k for k in sample(2, N_post, 2, size=3)')
+
+    # Just make sure using values between 0 and 1 work in principle
+    S3 = Synapses(G, G2, 'w:1', 'v+=w')
+    S3.connect(j='k for k in sample(2, N_post, 2, size=2)')
+
+    # Use pre-/post-synaptic variables for "stochastic" connections that are
+    # actually deterministic
+    S4 = Synapses(G, G2, 'w:1', on_pre='v+=w')
+    S4.connect(j='k for k in sample(2, N_post, 2, size=int(x_pre==2)*3)')
+
+    with catch_logs() as _:  # Ignore warnings about empty synapses
+        run(0*ms)  # for standalone
+
+    assert len(S1) == 0
+    S2_comp = np.zeros((len(G), len(G2)))
+    S2_comp[:, 2::2] = 1
+    _compare(S2, S2_comp)
+    assert len(S3) == len(G) * 2
+    assert all(S3.N_outgoing_pre == 2)
+    assert all(S3.j[:] % 2 == 0)
+    assert all(S3.j >= 2)
+    assert all([len(S3.j[x, :]) == len(set(S3.j[x, :]))
+                for x in range(len(G))])
+    assert len(S4) == 3
+    assert_equal(S4.i, np.ones(3)*2)
+    assert_equal(S4.j, np.arange(2, 7, 2))
+
+
+@pytest.mark.standalone_compatible
+def test_synapse_generator_fixed_random_negative_steps():
+    # Test generator with fixed-size sampling from stepped ranges (e.g. all
+    # even numbers) going backwards
+    G = NeuronGroup(4, '''v: 1
+                          x : integer''', threshold='False')
+    G.x = 'i'
+    G2 = NeuronGroup(7, '''v: 1
+                           y : 1''', threshold='False')
+    G2.y = '1.0*i/N'
+
+    S1 = Synapses(G, G2, 'w:1', 'v+=w')
+    S1.connect(j='k for k in sample(N_post-1, 0, -2, size=0)')
+
+    S2 = Synapses(G, G2, 'w:1', 'v+=w')
+    S2.connect(j='k for k in sample(N_post-1, 0, -2, size=3)')
+
+    # Just make sure using intermediate values between 0 and 1 work in principle
+    S3 = Synapses(G, G2, 'w:1', 'v+=w')
+    S3.connect(j='k for k in sample(N_post-1, 0, -2, size=2)')
+
+    # Use pre-/post-synaptic variables for "stochastic" connections that are
+    # actually deterministic
+    S4 = Synapses(G, G2, 'w:1', on_pre='v+=w')
+    S4.connect(j='k for k in sample(N_post-1, 0, -2, size=int(x_pre==2)*3)')
+
+    with catch_logs() as _:  # Ignore warnings about empty synapses
+        run(0*ms)  # for standalone
+
+    assert len(S1) == 0
+    S2_comp = np.zeros((len(G), len(G2)))
+    S2_comp[:, 2::2] = 1
+    _compare(S2, S2_comp)
+    assert len(S3) == len(G) * 2
+    assert all(S3.N_outgoing_pre == 2)
+    assert all(S3.j[:] % 2 == 0)
+    assert all(S3.j >= 2)
+    assert all([len(S3.j[x, :]) == len(set(S3.j[x, :]))
+                for x in range(len(G))])
+    assert len(S4) == 3
+    assert_equal(S4.i, np.ones(3) * 2)
+    assert_equal(S4.j, np.arange(6, 0, -2))
+
+@pytest.mark.standalone_compatible
+def test_synapse_generator_fixed_random_error1():
+    set_device('cpp_standalone')
+    G = NeuronGroup(5, '')
+    G2 = NeuronGroup(7, '')
+    S = Synapses(G, G2)
+    with pytest.raises((BrianObjectException, IndexError, RuntimeError)):
+        # Won't work for i=4
+        S.connect(j='k for k in sample(N_post, size=i+4)')
+        run(0*ms)  # for standalone
+
+
+@pytest.mark.standalone_compatible
+def test_synapse_generator_fixed_random_error2():
+    G = NeuronGroup(5, '')
+    G2 = NeuronGroup(7, '')
+    S = Synapses(G, G2)
+    with pytest.raises((BrianObjectException, IndexError, RuntimeError)):
+        # Won't work for i=4
+        S.connect(j='k for k in sample(N_post, size=3-i)')
+        run(0*ms)  # for standalone
+
+
+@pytest.mark.standalone_compatible
+def test_synapse_generator_fixed_random_skip_if_invalid():
+    G = NeuronGroup(5, '')
+    G2 = NeuronGroup(7, '')
+    S1 = Synapses(G, G2)
+    S2 = Synapses(G, G2)
+    # > N_post for i=4
+    S1.connect(j='k for k in sample(N_post, size=i+4)', skip_if_invalid=True)
+    # < 0 for i=4
+    S2.connect(j='k for k in sample(N_post, size=3-i)', skip_if_invalid=True)
+    run(0*ms)  # for standalone
+    assert_array_equal(S1.N_outgoing_pre, [4, 5, 6, 7, 7])
+    assert_array_equal(S2.N_outgoing_pre, [3, 2, 1, 0, 0])
+
+
 @pytest.mark.standalone_compatible
 def test_synapse_generator_random_with_condition():
     G = NeuronGroup(4, 'v: 1', threshold='False')
diff --git a/docs_sphinx/user/synapses.rst b/docs_sphinx/user/synapses.rst
index 064bd24b..070f3e90 100644
--- a/docs_sphinx/user/synapses.rst
+++ b/docs_sphinx/user/synapses.rst
@@ -377,13 +377,18 @@ created. Finally, ``RANGE`` can be either:
 1. a Python ``range``, e.g. ``range(N)`` is the integers from
    0 to N-1, ``range(A, B)`` is the integers from A to B-1,
    ``range(low, high, step)`` is the integers from ``low`` to
-   ``high-1`` with steps of size ``step``, or
-2. it can be a random sample ``sample(N, p=0.1)`` gives a
+   ``high-1`` with steps of size ``step``;
+2. a random sample ``sample(N, p=0.1)`` gives a
    random sample of integers from 0 to N-1 with 10% probability
    of each integer appearing in the sample. This can have extra
    arguments like range, e.g. ``sample(low, high, step, p=0.1)``
    will give each integer in ``range(low, high, step)`` with
-   probability 10%.
+   probability 10%;
+3. a random sample ``sample(N, size=10)`` with a fixed size,
+   in this example 10 values chosen (without replacement) from
+   the integers from 0 to N-1. As for the random sample based on
+   a probability, the ``sample`` expression can take additional
+   arguments to sample from a restricted range.
 
 If you try to create an invalid synapse (i.e. connecting
 neurons that are outside the correct range) then you will get
@@ -398,6 +403,11 @@ that are outside the valid range::
 
     S.connect(j='i+(-1)**k for k in range(2)', skip_if_invalid=True)
 
+You can also use this argument to deal with random samples of
+incorrect size, i.e. a negative size or a size bigger than the
+total population size. With ``skip_if_invalid=True``, no error will
+be raised and a size of 0 or the population size will be used.
+
 Summed variables
 ----------------
 In many cases, the postsynaptic neuron has a variable that represents a sum of variables over all
@@ -567,7 +577,9 @@ interpreted in the following way:
     |        If uniform random number between 0 and 1 < p(i, j):
     |            Create n(i, j) synapses for (i, j)
 
-With the generator syntax ``j='EXPR for VAR in RANGE if COND'``, the interpretation is:
+With the generator syntax ``j='EXPR for VAR in RANGE if COND'`` (where the
+``RANGE`` can be a full range or a random sample as described above), the interpretation
+is:
 
     | For every i:
     |     for every VAR in RANGE:
['brian2/codegen/runtime/numpy_rt/templates/synapses_create_generator.py_', 'brian2/tests/test_synapses.py', 'brian2/codegen/runtime/cython_rt/templates/synapses_create_generator.pyx', 'brian2/devices/cpp_standalone/templates/synapses_create_generator.cpp', 'brian2/synapses/synapses.py', 'docs_sphinx/user/synapses.rst', 'brian2/codegen/runtime/cython_rt/templates/common_group.pyx', 'brian2/core/base.py']