Repository: brian-team/brian2 (language: python)

Issue:
Random number generation on C++ standalone is slow
I only just now realized that random number generation on C++ standalone is slow compared to Cython/numpy. E.g., with a trivial example like:
```python
G = NeuronGroup(10000, 'x : 1')
G.run_regularly('x = rand()')
```
the `run_regularly` code object uses ~2.5s on Cython and numpy, and ~7s on C++ standalone. My hunch was that this is because in Cython/numpy, we create a buffer of random numbers at once using `numpy`'s `rand` function, whereas in C++ standalone, we call `rk_double` for each number we need. But it might be more fundamental: numpy switched to a different random number generator a couple of years ago, while in C++ standalone, we are still using the Mersenne-Twister "legacy code"…

Not sure how feasible/easy it is to switch to the same generator in C++ standalone. 

This could also be the time to tackle #19, but I don't see any obvious/simple solution.


Discussion / hints:
Note to self, try TRNG: https://www.numbercrunch.de/trng/
And https://github.com/DEShawResearch/random123
Also see: https://github.com/google/jax/blob/main/docs/jep/263-prng.md

Fix patch:
diff --git a/LICENSE b/LICENSE
index c06e706f5..0534996dd 100644
--- a/LICENSE
+++ b/LICENSE
@@ -659,10 +659,10 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 
--------------------------------------------------------------------------------
+--
 
-The Mersenne-Twister implementation "Randomkit" in brian2.random.randomkit is
-based on code by Jean-Sebastien Roy, provided under the MIT license:
+The `rand()` and `randn()` implementations are based on code by
+Jean-Sebastien Roy, provided under the MIT license:
 
 Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
 
@@ -687,47 +687,6 @@ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 --
 
-The rk_random and rk_seed functions algorithms and the original design of
-the Mersenne Twister RNG:
-
-Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions
-are met:
-
-1. Redistributions of source code must retain the above copyright
-notice, this list of conditions and the following disclaimer.
-
-2. Redistributions in binary form must reproduce the above copyright
-notice, this list of conditions and the following disclaimer in the
-documentation and/or other materials provided with the distribution.
-
-3. The names of its contributors may not be used to endorse or promote
-products derived from this software without specific prior written
-permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
-A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
-CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
-NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
---
-
-Original algorithm for the implementation of rk_interval function from
-Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by
-Magnus Jonsson.
-
---
-
 Constants used in the rk_double implementation by Isaku Wada.
 
 --
diff --git a/brian2/devices/cpp_standalone/codeobject.py b/brian2/devices/cpp_standalone/codeobject.py
index 9778c538d..50fdec9cb 100644
--- a/brian2/devices/cpp_standalone/codeobject.py
+++ b/brian2/devices/cpp_standalone/codeobject.py
@@ -138,15 +138,11 @@ def generate_rand_code(rand_func, owner):
         thread_number = "0"
     else:
         thread_number = "omp_get_thread_num()"
-    if rand_func == "rand":
-        rk_call = "rk_double"
-    elif rand_func == "randn":
-        rk_call = "rk_gauss"
-    else:
+    if rand_func not in ["rand", "randn"]:
         raise AssertionError(rand_func)
     code = """
-           double _%RAND_FUNC%(const int _vectorisation_idx) {
-               return %RK_CALL%(brian::_mersenne_twister_states[%THREAD_NUMBER%]);
+           inline double _%RAND_FUNC%(const int _vectorisation_idx) {
+               return brian::_random_generators[%THREAD_NUMBER%].%RAND_FUNC%();
            }
            """
     code = replace(
@@ -154,7 +150,6 @@ def generate_rand_code(rand_func, owner):
         {
             "%THREAD_NUMBER%": thread_number,
             "%RAND_FUNC%": rand_func,
-            "%RK_CALL%": rk_call,
         },
     )
     return {"support_code": code}
diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py
index a4050a512..0ceed3c15 100644
--- a/brian2/devices/cpp_standalone/device.py
+++ b/brian2/devices/cpp_standalone/device.py
@@ -19,7 +19,6 @@
 
 import numpy as np
 
-import brian2
 from brian2.codegen.codeobject import check_compiler_kwds
 from brian2.codegen.cpp_prefs import get_compiler_and_args, get_msvc_env
 from brian2.codegen.generators.cpp_generator import c_data_type
@@ -228,8 +227,8 @@ def __init__(self):
         self.extra_compile_args = []
         self.define_macros = []
         self.headers = []
-        self.include_dirs = ["brianlib/randomkit"]
-        self.library_dirs = ["brianlib/randomkit"]
+        self.include_dirs = []
+        self.library_dirs = []
         self.runtime_library_dirs = []
         self.run_environment_variables = {}
         if sys.platform.startswith("darwin"):
@@ -873,13 +872,10 @@ def generate_main_source(self, writer):
                     nb_threads = 1
                 main_lines.append(f"for (int _i=0; _i<{nb_threads}; _i++)")
                 if seed is None:  # random
-                    main_lines.append(
-                        "    rk_randomseed(brian::_mersenne_twister_states[_i]);"
-                    )
+                    main_lines.append("    brian::_random_generators[_i].seed();")
                 else:
                     main_lines.append(
-                        f"    rk_seed({seed!r}L + _i,"
-                        " brian::_mersenne_twister_states[_i]);"
+                        f"    brian::_random_generators[_i].seed({seed!r}L + _i);"
                     )
             else:
                 raise NotImplementedError(f"Unknown main queue function type {func}")
@@ -1084,28 +1080,6 @@ def copy_source_files(self, writer, directory):
             os.path.join(directory, "brianlib", "stdint_compat.h"),
         )
 
-        # Copy the RandomKit implementation
-        if not os.path.exists(os.path.join(directory, "brianlib", "randomkit")):
-            os.mkdir(os.path.join(directory, "brianlib", "randomkit"))
-        shutil.copy2(
-            os.path.join(
-                os.path.split(inspect.getsourcefile(brian2))[0],
-                "random",
-                "randomkit",
-                "randomkit.c",
-            ),
-            os.path.join(directory, "brianlib", "randomkit", "randomkit.c"),
-        )
-        shutil.copy2(
-            os.path.join(
-                os.path.split(inspect.getsourcefile(brian2))[0],
-                "random",
-                "randomkit",
-                "randomkit.h",
-            ),
-            os.path.join(directory, "brianlib", "randomkit", "randomkit.h"),
-        )
-
     def _insert_func_namespace(self, func, code_object, namespace):
         impl = func.implementations[self.code_object_class()]
         func_namespace = impl.get_namespace(code_object.owner)
@@ -1583,9 +1557,7 @@ def build(
             for codeobj in self.code_objects.values()
             for source_file in codeobj.compiler_kwds.get("sources", [])
         ]
-        additional_source_files += codeobj_source_files + [
-            "brianlib/randomkit/randomkit.c"
-        ]
+        additional_source_files += codeobj_source_files
 
         for d in ["code_objects", "results", "static_arrays"]:
             ensure_directory(os.path.join(directory, d))
@@ -1703,7 +1675,6 @@ def delete(self, code=True, data=True, directory=True, force=False):
             fnames.extend(
                 [
                     os.path.join("brianlib", "spikequeue.h"),
-                    os.path.join("brianlib", "randomkit", "randomkit.h"),
                     os.path.join("brianlib", "stdint_compat.h"),
                 ]
             )
@@ -1731,7 +1702,6 @@ def delete(self, code=True, data=True, directory=True, force=False):
 
         if directory:
             directories = [
-                os.path.join("brianlib", "randomkit"),
                 "brianlib",
                 "code_objects",
                 "results",
diff --git a/brian2/devices/cpp_standalone/templates/main.cpp b/brian2/devices/cpp_standalone/templates/main.cpp
index 5758fff70..ba57df813 100644
--- a/brian2/devices/cpp_standalone/templates/main.cpp
+++ b/brian2/devices/cpp_standalone/templates/main.cpp
@@ -5,7 +5,6 @@
 {{ openmp_pragma('include') }}
 #include "run.h"
 #include "brianlib/common_math.h"
-#include "randomkit.h"
 
 {% for codeobj in code_objects | sort(attribute='name') %}
 #include "code_objects/{{codeobj.name}}.h"
@@ -36,6 +35,7 @@ void set_from_command_line(const std::vector<std::string> args)
 }
 int main(int argc, char **argv)
 {
+	std::random_device _rd;
 	std::vector<std::string> args(argv + 1, argv + argc);
 	if (args.size() >=2 && args[0] == "--results_dir")
 	{
diff --git a/brian2/devices/cpp_standalone/templates/objects.cpp b/brian2/devices/cpp_standalone/templates/objects.cpp
index 166caab29..fab6bf3a0 100644
--- a/brian2/devices/cpp_standalone/templates/objects.cpp
+++ b/brian2/devices/cpp_standalone/templates/objects.cpp
@@ -20,7 +20,7 @@ set_variable_from_value(name, {{array_name}}, var_size, (char)atoi(s_value.c_str
 #include "brianlib/dynamic_array.h"
 #include "brianlib/stdint_compat.h"
 #include "network.h"
-#include "randomkit.h"
+#include<random>
 #include<vector>
 #include<iostream>
 #include<fstream>
@@ -32,7 +32,11 @@ set_variable_from_value(name, {{array_name}}, var_size, (char)atoi(s_value.c_str
 namespace brian {
 
 std::string results_dir = "results/";  // can be overwritten by --results_dir command line arg
-std::vector< rk_state* > _mersenne_twister_states;
+
+// For multhreading, we need one generator for each thread. We also create a distribution for
+// each thread, even though this is not strictly necessary for the uniform distribution, as
+// the distribution is stateless.
+std::vector< RandomGenerator > _random_generators;
 
 //////////////// networks /////////////////
 {% for net in networks | sort(attribute='name') %}
@@ -76,15 +80,15 @@ template<class T> void set_variable_from_file(std::string varname, T* var_pointe
 
 //////////////// set arrays by name ///////
 void set_variable_by_name(std::string name, std::string s_value) {
-	size_t var_size;
-	size_t data_size;
-	// C-style or Python-style capitalization is allowed for boolean values
+    size_t var_size;
+    size_t data_size;
+    // C-style or Python-style capitalization is allowed for boolean values
     if (s_value == "true" || s_value == "True")
         s_value = "1";
     else if (s_value == "false" || s_value == "False")
         s_value = "0";
-	// non-dynamic arrays
-	{% for var, varname in array_specs | dictsort(by='value') %}
+    // non-dynamic arrays
+    {% for var, varname in array_specs | dictsort(by='value') %}
     {% if not var in dynamic_array_specs and not var.read_only %}
     if (name == "{{var.owner.name}}.{{var.name}}") {
         var_size = {{var.size}};
@@ -167,8 +171,8 @@ const int _num_{{name}} = {{N}};
 // {{S.name}}
 {% for path in S._pathways | sort(attribute='name') %}
 SynapticPathway {{path.name}}(
-		{{dynamic_array_specs[path.synapse_sources]}},
-		{{path.source.start}}, {{path.source.stop}});
+    {{dynamic_array_specs[path.synapse_sources]}},
+    {{path.source.start}}, {{path.source.stop}});
 {% endfor %}
 {% endfor %}
 
@@ -187,107 +191,108 @@ double {{codeobj}}_profiling_info = 0.0;
 
 void _init_arrays()
 {
-	using namespace brian;
+    using namespace brian;
 
     // Arrays initialized to 0
-	{% for var, varname in zero_arrays | sort(attribute='1') %}
-	{% if varname in dynamic_array_specs.values() %}
-	{{varname}}.resize({{var.size}});
-	{% else %}
-	{{varname}} = new {{c_data_type(var.dtype)}}[{{var.size}}];
-	{% endif %}
+    {% for var, varname in zero_arrays | sort(attribute='1') %}
+    {% if varname in dynamic_array_specs.values() %}
+    {{varname}}.resize({{var.size}});
+    {% else %}
+    {{varname}} = new {{c_data_type(var.dtype)}}[{{var.size}}];
+    {% endif %}
     {{ openmp_pragma('parallel-static')}}
-	for(int i=0; i<{{var.size}}; i++) {{varname}}[i] = 0;
+    for(int i=0; i<{{var.size}}; i++) {{varname}}[i] = 0;
 
-	{% endfor %}
+    {% endfor %}
 
-	// Arrays initialized to an "arange"
-	{% for var, varname, start in arange_arrays | sort(attribute='1')%}
-	{% if varname in dynamic_array_specs.values() %}
-	{{varname}}.resize({{var.size}});
-	{% else %}
-	{{varname}} = new {{c_data_type(var.dtype)}}[{{var.size}}];
-	{% endif %}
+    // Arrays initialized to an "arange"
+    {% for var, varname, start in arange_arrays | sort(attribute='1')%}
+    {% if varname in dynamic_array_specs.values() %}
+    {{varname}}.resize({{var.size}});
+    {% else %}
+    {{varname}} = new {{c_data_type(var.dtype)}}[{{var.size}}];
+    {% endif %}
     {{ openmp_pragma('parallel-static')}}
-	for(int i=0; i<{{var.size}}; i++) {{varname}}[i] = {{start}} + i;
-
-	{% endfor %}
-
-	// static arrays
-	{% for (name, dtype_spec, N, filename) in static_array_specs | sort %}
-	{% if name in dynamic_array_specs.values() %}
-	{{name}}.resize({{N}});
-	{% else %}
-	{{name}} = new {{dtype_spec}}[{{N}}];
-	{% endif %}
-	{% endfor %}
-
-	// Random number generator states
-	for (int i=0; i<{{openmp_pragma('get_num_threads')}}; i++)
-	    _mersenne_twister_states.push_back(new rk_state());
+    for(int i=0; i<{{var.size}}; i++) {{varname}}[i] = {{start}} + i;
+
+    {% endfor %}
+
+    // static arrays
+    {% for (name, dtype_spec, N, filename) in static_array_specs | sort %}
+    {% if name in dynamic_array_specs.values() %}
+    {{name}}.resize({{N}});
+    {% else %}
+    {{name}} = new {{dtype_spec}}[{{N}}];
+    {% endif %}
+    {% endfor %}
+
+    // Random number generator states
+    std::random_device rd;
+    for (int i=0; i<{{openmp_pragma('get_num_threads')}}; i++)
+        _random_generators.push_back(RandomGenerator());
 }
 
 void _load_arrays()
 {
-	using namespace brian;
-
-	{% for (name, dtype_spec, N, filename) in static_array_specs | sort %}
-	ifstream f{{name}};
-	f{{name}}.open("static_arrays/{{name}}", ios::in | ios::binary);
-	if(f{{name}}.is_open())
-	{
-	    {% if name in dynamic_array_specs.values() %}
-	    f{{name}}.read(reinterpret_cast<char*>(&{{name}}[0]), {{N}}*sizeof({{dtype_spec}}));
-	    {% else %}
-		f{{name}}.read(reinterpret_cast<char*>({{name}}), {{N}}*sizeof({{dtype_spec}}));
-		{% endif %}
-	} else
-	{
-		std::cout << "Error opening static array {{name}}." << endl;
-	}
-	{% endfor %}
+    using namespace brian;
+
+    {% for (name, dtype_spec, N, filename) in static_array_specs | sort %}
+    ifstream f{{name}};
+    f{{name}}.open("static_arrays/{{name}}", ios::in | ios::binary);
+    if(f{{name}}.is_open())
+    {
+        {% if name in dynamic_array_specs.values() %}
+        f{{name}}.read(reinterpret_cast<char*>(&{{name}}[0]), {{N}}*sizeof({{dtype_spec}}));
+        {% else %}
+        f{{name}}.read(reinterpret_cast<char*>({{name}}), {{N}}*sizeof({{dtype_spec}}));
+        {% endif %}
+    } else
+    {
+        std::cout << "Error opening static array {{name}}." << endl;
+    }
+    {% endfor %}
 }
 
 void _write_arrays()
 {
-	using namespace brian;
-
-	{% for var, varname in array_specs | dictsort(by='value') %}
-	{% if not (var in dynamic_array_specs or var in dynamic_array_2d_specs) %}
-	ofstream outfile_{{varname}};
-	outfile_{{varname}}.open(results_dir + "{{get_array_filename(var)}}", ios::binary | ios::out);
-	if(outfile_{{varname}}.is_open())
-	{
-		outfile_{{varname}}.write(reinterpret_cast<char*>({{varname}}), {{var.size}}*sizeof({{get_array_name(var)}}[0]));
-		outfile_{{varname}}.close();
-	} else
-	{
-		std::cout << "Error writing output file for {{varname}}." << endl;
-	}
-	{% endif %}
-	{% endfor %}
-
-	{% for var, varname in dynamic_array_specs | dictsort(by='value') %}
-	ofstream outfile_{{varname}};
-	outfile_{{varname}}.open(results_dir + "{{get_array_filename(var)}}", ios::binary | ios::out);
-	if(outfile_{{varname}}.is_open())
-	{
+    using namespace brian;
+
+    {% for var, varname in array_specs | dictsort(by='value') %}
+    {% if not (var in dynamic_array_specs or var in dynamic_array_2d_specs) %}
+    ofstream outfile_{{varname}};
+    outfile_{{varname}}.open(results_dir + "{{get_array_filename(var)}}", ios::binary | ios::out);
+    if(outfile_{{varname}}.is_open())
+    {
+        outfile_{{varname}}.write(reinterpret_cast<char*>({{varname}}), {{var.size}}*sizeof({{get_array_name(var)}}[0]));
+        outfile_{{varname}}.close();
+    } else
+    {
+        std::cout << "Error writing output file for {{varname}}." << endl;
+    }
+    {% endif %}
+    {% endfor %}
+
+    {% for var, varname in dynamic_array_specs | dictsort(by='value') %}
+    ofstream outfile_{{varname}};
+    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();
-		}
-	} else
-	{
-		std::cout << "Error writing output file for {{varname}}." << endl;
-	}
-	{% endfor %}
-
-	{% for var, varname in dynamic_array_2d_specs | dictsort(by='value') %}
-	ofstream outfile_{{varname}};
-	outfile_{{varname}}.open(results_dir + "{{get_array_filename(var)}}", ios::binary | ios::out);
-	if(outfile_{{varname}}.is_open())
-	{
+            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;
+    }
+    {% endfor %}
+
+    {% for var, varname in dynamic_array_2d_specs | dictsort(by='value') %}
+    ofstream outfile_{{varname}};
+    outfile_{{varname}}.open(results_dir + "{{get_array_filename(var)}}", ios::binary | ios::out);
+    if(outfile_{{varname}}.is_open())
+    {
         for (int n=0; n<{{varname}}.n; n++)
         {
             if (! {{varname}}(n).empty())
@@ -296,63 +301,63 @@ void _write_arrays()
             }
         }
         outfile_{{varname}}.close();
-	} else
-	{
-		std::cout << "Error writing output file for {{varname}}." << endl;
-	}
-	{% endfor %}
+    } else
+    {
+        std::cout << "Error writing output file for {{varname}}." << endl;
+    }
+    {% endfor %}
     {% if profiled_codeobjects is defined and profiled_codeobjects %}
-	// Write profiling info to disk
-	ofstream outfile_profiling_info;
-	outfile_profiling_info.open(results_dir + "profiling_info.txt", ios::out);
-	if(outfile_profiling_info.is_open())
-	{
-	{% for codeobj in profiled_codeobjects | sort %}
-	outfile_profiling_info << "{{codeobj}}\t" << {{codeobj}}_profiling_info << std::endl;
-	{% endfor %}
-	outfile_profiling_info.close();
-	} else
-	{
-	    std::cout << "Error writing profiling info to file." << std::endl;
-	}
+    // Write profiling info to disk
+    ofstream outfile_profiling_info;
+    outfile_profiling_info.open(results_dir + "profiling_info.txt", ios::out);
+    if(outfile_profiling_info.is_open())
+    {
+    {% for codeobj in profiled_codeobjects | sort %}
+    outfile_profiling_info << "{{codeobj}}\t" << {{codeobj}}_profiling_info << std::endl;
+    {% endfor %}
+    outfile_profiling_info.close();
+    } else
+    {
+        std::cout << "Error writing profiling info to file." << std::endl;
+    }
     {% endif %}
-	// Write last run info to disk
-	ofstream outfile_last_run_info;
-	outfile_last_run_info.open(results_dir + "last_run_info.txt", ios::out);
-	if(outfile_last_run_info.is_open())
-	{
-		outfile_last_run_info << (Network::_last_run_time) << " " << (Network::_last_run_completed_fraction) << std::endl;
-		outfile_last_run_info.close();
-	} else
-	{
-	    std::cout << "Error writing last run info to file." << std::endl;
-	}
+    // Write last run info to disk
+    ofstream outfile_last_run_info;
+    outfile_last_run_info.open(results_dir + "last_run_info.txt", ios::out);
+    if(outfile_last_run_info.is_open())
+    {
+        outfile_last_run_info << (Network::_last_run_time) << " " << (Network::_last_run_completed_fraction) << std::endl;
+        outfile_last_run_info.close();
+    } else
+    {
+        std::cout << "Error writing last run info to file." << std::endl;
+    }
 }
 
 void _dealloc_arrays()
 {
-	using namespace brian;
-
-	{% for var, varname in array_specs | dictsort(by='value') %}
-	{% if varname in dynamic_array_specs.values() %}
-	if({{varname}}!=0)
-	{
-		delete [] {{varname}};
-		{{varname}} = 0;
-	}
-	{% endif %}
-	{% endfor %}
-
-	// static arrays
-	{% for (name, dtype_spec, N, filename) in static_array_specs | sort %}
-	{% if not name in dynamic_array_specs.values() %}
-	if({{name}}!=0)
-	{
-		delete [] {{name}};
-		{{name}} = 0;
-	}
-	{% endif %}
-	{% endfor %}
+    using namespace brian;
+
+    {% for var, varname in array_specs | dictsort(by='value') %}
+    {% if varname in dynamic_array_specs.values() %}
+    if({{varname}}!=0)
+    {
+        delete [] {{varname}};
+        {{varname}} = 0;
+    }
+    {% endif %}
+    {% endfor %}
+
+    // static arrays
+    {% for (name, dtype_spec, N, filename) in static_array_specs | sort %}
+    {% if not name in dynamic_array_specs.values() %}
+    if({{name}}!=0)
+    {
+        delete [] {{name}};
+        {{name}} = 0;
+    }
+    {% endif %}
+    {% endfor %}
 }
 
 {% endmacro %}
@@ -369,15 +374,67 @@ void _dealloc_arrays()
 #include "brianlib/dynamic_array.h"
 #include "brianlib/stdint_compat.h"
 #include "network.h"
-#include "randomkit.h"
+#include<random>
 #include<vector>
 {{ openmp_pragma('include') }}
 
 namespace brian {
 
 extern std::string results_dir;
+
+class RandomGenerator {
+    private:
+        std::mt19937 gen;
+        double stored_gauss;
+        bool has_stored_gauss = false;
+    public:
+        RandomGenerator() {
+            seed();
+        }
+        void seed() {
+            std::random_device rd;
+            gen.seed(rd());
+            has_stored_gauss = false;
+        }
+        void seed(unsigned long seed) {
+            gen.seed(seed);
+            has_stored_gauss = false;
+        }
+        double rand() {
+            /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */
+            const long a = gen() >> 5;
+            const long b = gen() >> 6;
+            return (a * 67108864.0 + b) / 9007199254740992.0;
+        }
+
+        double randn() {
+            if (has_stored_gauss) {
+                const double tmp = stored_gauss;
+                has_stored_gauss = false;
+                return tmp;
+            }
+            else {
+                double f, x1, x2, r2;
+
+                do {
+                    x1 = 2.0*rand() - 1.0;
+                    x2 = 2.0*rand() - 1.0;
+                    r2 = x1*x1 + x2*x2;
+                }
+                while (r2 >= 1.0 || r2 == 0.0);
+
+                /* Box-Muller transform */
+                f = sqrt(-2.0*log(r2)/r2);
+                /* Keep for next call */
+                stored_gauss = f*x1;
+                has_stored_gauss = true;
+                return f*x2;
+            }
+        }
+};
+
 // In OpenMP we need one state per thread
-extern std::vector< rk_state* > _mersenne_twister_states;
+extern std::vector< RandomGenerator > _random_generators;
 
 //////////////// clocks ///////////////////
 {% for clock in clocks | sort(attribute='name') %}
diff --git a/brian2/devices/cpp_standalone/templates/run.cpp b/brian2/devices/cpp_standalone/templates/run.cpp
index 20a429bbf..798493dd0 100644
--- a/brian2/devices/cpp_standalone/templates/run.cpp
+++ b/brian2/devices/cpp_standalone/templates/run.cpp
@@ -2,7 +2,7 @@
 #include<stdlib.h>
 #include "objects.h"
 #include<ctime>
-#include "randomkit.h"
+#include<random>
 
 {% for codeobj in code_objects | sort(attribute='name') %}
 #include "code_objects/{{codeobj.name}}.h"
@@ -22,8 +22,6 @@ void brian_start()
     brian::{{clock.name}}.dt = brian::{{array_specs[clock.variables['dt']]}};
     brian::{{clock.name}}.t = brian::{{array_specs[clock.variables['t']]}};
     {% endfor %}
-    for (int i=0; i<{{openmp_pragma('get_num_threads')}}; i++)
-	    rk_randomseed(brian::_mersenne_twister_states[i]);  // Note that this seed can be potentially replaced in main.cpp
 }
 
 void brian_end()
diff --git a/brian2/random/randomkit/randomkit.c b/brian2/random/randomkit/randomkit.c
deleted file mode 100644
index 169c288dc..000000000
--- a/brian2/random/randomkit/randomkit.c
+++ /dev/null
@@ -1,402 +0,0 @@
-/* Random kit 1.3 */
-
-/*
- * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
- *
- * The rk_random and rk_seed functions algorithms and the original design of
- * the Mersenne Twister RNG:
- *
- *   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
- *   All rights reserved.
- *
- *   Redistribution and use in source and binary forms, with or without
- *   modification, are permitted provided that the following conditions
- *   are met:
- *
- *   1. Redistributions of source code must retain the above copyright
- *   notice, this list of conditions and the following disclaimer.
- *
- *   2. Redistributions in binary form must reproduce the above copyright
- *   notice, this list of conditions and the following disclaimer in the
- *   documentation and/or other materials provided with the distribution.
- *
- *   3. The names of its contributors may not be used to endorse or promote
- *   products derived from this software without specific prior written
- *   permission.
- *
- *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
- *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
- *   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
- *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- *   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- *   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- *   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- *   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- *   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- *   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- *
- * Original algorithm for the implementation of rk_interval function from
- * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by
- * Magnus Jonsson.
- *
- * Constants used in the rk_double implementation by Isaku Wada.
- *
- * Permission is hereby granted, free of charge, to any person obtaining a
- * copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject to
- * the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
- */
-
-/* static char const rcsid[] =
-  "@(#) $Jeannot: randomkit.c,v 1.28 2005/07/21 22:14:09 js Exp $"; */
-#include <stddef.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <errno.h>
-#include <limits.h>
-#include <math.h>
-
-#ifdef _WIN32
-/*
- * Windows
- * XXX: we have to use this ugly defined(__GNUC__) because it is not easy to
- * detect the compiler used in distutils itself
- */
-#if (defined(__GNUC__) && defined(NPY_NEEDS_MINGW_TIME_WORKAROUND))
-
-/*
- * FIXME: ideally, we should set this to the real version of MSVCRT. We need
- * something higher than 0x601 to enable _ftime64 and co
- */
-#define __MSVCRT_VERSION__ 0x0700
-#include <time.h>
-#include <sys/timeb.h>
-
-/*
- * mingw msvcr lib import wrongly export _ftime, which does not exist in the
- * actual msvc runtime for version >= 8; we make it an alias to _ftime64, which
- * is available in those versions of the runtime
- */
-#define _FTIME(x) _ftime64((x))
-#else
-#include <time.h>
-#include <sys/timeb.h>
-#define _FTIME(x) _ftime((x))
-#endif
-
-#ifndef RK_NO_WINCRYPT
-/* Windows crypto */
-#ifndef _WIN32_WINNT
-#define _WIN32_WINNT 0x0400
-#endif
-#include <windows.h>
-#include <wincrypt.h>
-#endif
-
-#else
-/* Unix */
-#include <time.h>
-#include <sys/time.h>
-#include <unistd.h>
-#endif
-
-#include "randomkit.h"
-
-#ifndef RK_DEV_URANDOM
-#define RK_DEV_URANDOM "/dev/urandom"
-#endif
-
-#ifndef RK_DEV_RANDOM
-#define RK_DEV_RANDOM "/dev/random"
-#endif
-
-char *rk_strerror[RK_ERR_MAX] =
-{
-    "no error",
-    "random device unvavailable"
-};
-
-/* static functions */
-static unsigned long rk_hash(unsigned long key);
-
-void
-rk_seed(unsigned long seed, rk_state *state)
-{
-    int pos;
-    seed &= 0xffffffffUL;
-
-    /* Knuth's PRNG as used in the Mersenne Twister reference implementation */
-    for (pos = 0; pos < RK_STATE_LEN; pos++) {
-        state->key[pos] = seed;
-        seed = (1812433253UL * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffUL;
-    }
-    state->pos = RK_STATE_LEN;
-    state->gauss = 0;
-    state->has_gauss = 0;
-    state->has_binomial = 0;
-}
-
-/* Thomas Wang 32 bits integer hash function */
-unsigned long
-rk_hash(unsigned long key)
-{
-    key += ~(key << 15);
-    key ^=  (key >> 10);
-    key +=  (key << 3);
-    key ^=  (key >> 6);
-    key += ~(key << 11);
-    key ^=  (key >> 16);
-    return key;
-}
-
-rk_error
-rk_randomseed(rk_state *state)
-{
-#ifndef _WIN32
-    struct timeval tv;
-#else
-    struct _timeb  tv;
-#endif
-    int i;
-
-    if (rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) {
-        /* ensures non-zero key */
-        state->key[0] |= 0x80000000UL;
-        state->pos = RK_STATE_LEN;
-        state->gauss = 0;
-        state->has_gauss = 0;
-        state->has_binomial = 0;
-
-        for (i = 0; i < 624; i++) {
-            state->key[i] &= 0xffffffffUL;
-        }
-        return RK_NOERR;
-    }
-
-#ifndef _WIN32
-    gettimeofday(&tv, NULL);
-    rk_seed(rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec)
-            ^ rk_hash(clock()), state);
-#else
-    _FTIME(&tv);
-    rk_seed(rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()), state);
-#endif
-
-    return RK_ENODEV;
-}
-
-/* Magic Mersenne Twister constants */
-#define N 624
-#define M 397
-#define MATRIX_A 0x9908b0dfUL
-#define UPPER_MASK 0x80000000UL
-#define LOWER_MASK 0x7fffffffUL
-
-/* Slightly optimised reference implementation of the Mersenne Twister */
-unsigned long
-rk_random(rk_state *state)
-{
-    unsigned long y;
-
-    if (state->pos == RK_STATE_LEN) {
-        int i;
-
-        for (i = 0; i < N - M; i++) {
-            y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK);
-            state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A);
-        }
-        for (; i < N - 1; i++) {
-            y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK);
-            state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A);
-        }
-        y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK);
-        state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
-
-        state->pos = 0;
-    }
-    y = state->key[state->pos++];
-
-    /* Tempering */
-    y ^= (y >> 11);
-    y ^= (y << 7) & 0x9d2c5680UL;
-    y ^= (y << 15) & 0xefc60000UL;
-    y ^= (y >> 18);
-
-    return y;
-}
-
-long
-rk_long(rk_state *state)
-{
-    return rk_ulong(state) >> 1;
-}
-
-unsigned long
-rk_ulong(rk_state *state)
-{
-#if ULONG_MAX <= 0xffffffffUL
-    return rk_random(state);
-#else
-    return (rk_random(state) << 32) | (rk_random(state));
-#endif
-}
-
-unsigned long
-rk_interval(unsigned long max, rk_state *state)
-{
-    unsigned long mask = max, value;
-
-    if (max == 0) {
-        return 0;
-    }
-    /* Smallest bit mask >= max */
-    mask |= mask >> 1;
-    mask |= mask >> 2;
-    mask |= mask >> 4;
-    mask |= mask >> 8;
-    mask |= mask >> 16;
-#if ULONG_MAX > 0xffffffffUL
-    mask |= mask >> 32;
-#endif
-
-    /* Search a random value in [0..mask] <= max */
-#if ULONG_MAX > 0xffffffffUL
-    if (max <= 0xffffffffUL) {
-        while ((value = (rk_random(state) & mask)) > max);
-    }
-    else {
-        while ((value = (rk_ulong(state) & mask)) > max);
-    }
-#else
-    while ((value = (rk_ulong(state) & mask)) > max);
-#endif
-    return value;
-}
-
-double
-rk_double(rk_state *state)
-{
-    /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */
-    long a = rk_random(state) >> 5, b = rk_random(state) >> 6;
-    return (a * 67108864.0 + b) / 9007199254740992.0;
-}
-
-void
-rk_fill(void *buffer, size_t size, rk_state *state)
-{
-    unsigned long r;
-    unsigned char *buf = (unsigned char *)buffer;
-
-    for (; size >= 4; size -= 4) {
-        r = rk_random(state);
-        *(buf++) = r & 0xFF;
-        *(buf++) = (r >> 8) & 0xFF;
-        *(buf++) = (r >> 16) & 0xFF;
-        *(buf++) = (r >> 24) & 0xFF;
-    }
-
-    if (!size) {
-        return;
-    }
-    r = rk_random(state);
-    for (; size; r >>= 8, size --) {
-        *(buf++) = (unsigned char)(r & 0xFF);
-    }
-}
-
-rk_error
-rk_devfill(void *buffer, size_t size, int strong)
-{
-#ifndef _WIN32
-    FILE *rfile;
-    int done;
-
-    if (strong) {
-        rfile = fopen(RK_DEV_RANDOM, "rb");
-    }
-    else {
-        rfile = fopen(RK_DEV_URANDOM, "rb");
-    }
-    if (rfile == NULL) {
-        return RK_ENODEV;
-    }
-    done = fread(buffer, size, 1, rfile);
-    fclose(rfile);
-    if (done) {
-        return RK_NOERR;
-    }
-#else
-
-#ifndef RK_NO_WINCRYPT
-    HCRYPTPROV hCryptProv;
-    BOOL done;
-
-    if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL,
-            CRYPT_VERIFYCONTEXT) || !hCryptProv) {
-        return RK_ENODEV;
-    }
-    done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer);
-    CryptReleaseContext(hCryptProv, 0);
-    if (done) {
-        return RK_NOERR;
-    }
-#endif
-
-#endif
-    return RK_ENODEV;
-}
-
-rk_error
-rk_altfill(void *buffer, size_t size, int strong, rk_state *state)
-{
-    rk_error err;
-
-    err = rk_devfill(buffer, size, strong);
-    if (err) {
-        rk_fill(buffer, size, state);
-    }
-    return err;
-}
-
-double
-rk_gauss(rk_state *state)
-{
-    if (state->has_gauss) {
-        const double tmp = state->gauss;
-        state->gauss = 0;
-        state->has_gauss = 0;
-        return tmp;
-    }
-    else {
-        double f, x1, x2, r2;
-
-        do {
-            x1 = 2.0*rk_double(state) - 1.0;
-            x2 = 2.0*rk_double(state) - 1.0;
-            r2 = x1*x1 + x2*x2;
-        }
-        while (r2 >= 1.0 || r2 == 0.0);
-
-        /* Box-Muller transform */
-        f = sqrt(-2.0*log(r2)/r2);
-        /* Keep for next call */
-        state->gauss = f*x1;
-        state->has_gauss = 1;
-        return f*x2;
-    }
-}
diff --git a/brian2/random/randomkit/randomkit.h b/brian2/random/randomkit/randomkit.h
deleted file mode 100644
index e049488ee..000000000
--- a/brian2/random/randomkit/randomkit.h
+++ /dev/null
@@ -1,189 +0,0 @@
-/* Random kit 1.3 */
-
-/*
- * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
- *
- * Permission is hereby granted, free of charge, to any person obtaining a
- * copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject to
- * the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
- */
-
-/* @(#) $Jeannot: randomkit.h,v 1.24 2005/07/21 22:14:09 js Exp $ */
-
-/*
- * Typical use:
- *
- * {
- *  rk_state state;
- *  unsigned long seed = 1, random_value;
- *
- *  rk_seed(seed, &state); // Initialize the RNG
- *  ...
- *  random_value = rk_random(&state); // Generate random values in [0..RK_MAX]
- * }
- *
- * Instead of rk_seed, you can use rk_randomseed which will get a random seed
- * from /dev/urandom (or the clock, if /dev/urandom is unavailable):
- *
- * {
- *  rk_state state;
- *  unsigned long random_value;
- *
- *  rk_randomseed(&state); // Initialize the RNG with a random seed
- *  ...
- *  random_value = rk_random(&state); // Generate random values in [0..RK_MAX]
- * }
- */
-
-/*
- * Useful macro:
- *  RK_DEV_RANDOM: the device used for random seeding.
- *                 defaults to "/dev/urandom"
- */
-
-#include <stddef.h>
-
-#ifndef _RANDOMKIT_
-#define _RANDOMKIT_
-
-#define RK_STATE_LEN 624
-
-typedef struct rk_state_
-{
-    unsigned long key[RK_STATE_LEN];
-    int pos;
-    int has_gauss; /* !=0: gauss contains a gaussian deviate */
-    double gauss;
-
-    /* The rk_state structure has been extended to store the following
-     * information for the binomial generator. If the input values of n or p
-     * are different than nsave and psave, then the other parameters will be
-     * recomputed. RTK 2005-09-02 */
-
-    int has_binomial; /* !=0: following parameters initialized for
-                              binomial */
-    double psave;
-    long nsave;
-    double r;
-    double q;
-    double fm;
-    long m;
-    double p1;
-    double xm;
-    double xl;
-    double xr;
-    double c;
-    double laml;
-    double lamr;
-    double p2;
-    double p3;
-    double p4;
-
-}
-rk_state;
-
-typedef enum {
-    RK_NOERR = 0, /* no error */
-    RK_ENODEV = 1, /* no RK_DEV_RANDOM device */
-    RK_ERR_MAX = 2
-} rk_error;
-
-/* error strings */
-extern char *rk_strerror[RK_ERR_MAX];
-
-/* Maximum generated random value */
-#define RK_MAX 0xFFFFFFFFUL
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-/*
- * Initialize the RNG state using the given seed.
- */
-extern void rk_seed(unsigned long seed, rk_state *state);
-
-/*
- * Initialize the RNG state using a random seed.
- * Uses /dev/random or, when unavailable, the clock (see randomkit.c).
- * Returns RK_NOERR when no errors occurs.
- * Returns RK_ENODEV when the use of RK_DEV_RANDOM failed (for example because
- * there is no such device). In this case, the RNG was initialized using the
- * clock.
- */
-extern rk_error rk_randomseed(rk_state *state);
-
-/*
- * Returns a random unsigned long between 0 and RK_MAX inclusive
- */
-extern unsigned long rk_random(rk_state *state);
-
-/*
- * Returns a random long between 0 and LONG_MAX inclusive
- */
-extern long rk_long(rk_state *state);
-
-/*
- * Returns a random unsigned long between 0 and ULONG_MAX inclusive
- */
-extern unsigned long rk_ulong(rk_state *state);
-
-/*
- * Returns a random unsigned long between 0 and max inclusive.
- */
-extern unsigned long rk_interval(unsigned long max, rk_state *state);
-
-/*
- * Returns a random double between 0.0 and 1.0, 1.0 excluded.
- */
-extern double rk_double(rk_state *state);
-
-/*
- * fill the buffer with size random bytes
- */
-extern void rk_fill(void *buffer, size_t size, rk_state *state);
-
-/*
- * fill the buffer with randombytes from the random device
- * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is
- * On Unix, if strong is defined, RK_DEV_RANDOM is used. If not, RK_DEV_URANDOM
- * is used instead. This parameter has no effect on Windows.
- * Warning: on most unixes RK_DEV_RANDOM will wait for enough entropy to answer
- * which can take a very long time on quiet systems.
- */
-extern rk_error rk_devfill(void *buffer, size_t size, int strong);
-
-/*
- * fill the buffer using rk_devfill if the random device is available and using
- * rk_fill if is is not
- * parameters have the same meaning as rk_fill and rk_devfill
- * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is
- */
-extern rk_error rk_altfill(void *buffer, size_t size, int strong,
-                            rk_state *state);
-
-/*
- * return a random gaussian deviate with variance unity and zero mean.
- */
-extern double rk_gauss(rk_state *state);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif /* _RANDOMKIT_ */


Test patch:
diff --git a/brian2/tests/test_neurongroup.py b/brian2/tests/test_neurongroup.py
index dc96a8526..28b094f11 100644
--- a/brian2/tests/test_neurongroup.py
+++ b/brian2/tests/test_neurongroup.py
@@ -11,7 +11,7 @@
 from brian2.core.network import Network
 from brian2.core.preferences import prefs
 from brian2.core.variables import linked_var
-from brian2.devices.device import device, seed
+from brian2.devices.device import device, get_device, seed
 from brian2.equations.equations import Equations
 from brian2.groups.group import get_dtype
 from brian2.groups.neurongroup import NeuronGroup
@@ -1986,6 +1986,60 @@ def test_random_values_fixed_seed():
     assert_allclose(G.v1[:], G.v2[:])
 
 
+_random_values = {
+    ("RuntimeDevice", "numpy", None): (
+        [0.1636023, 0.76229608, 0.74945305, 0.82121212, 0.82669968],
+        [-0.7758696, 0.13295831, 0.87360834, -1.21879122, 0.62980314],
+    ),
+    ("RuntimeDevice", "cython", None): (
+        [0.1636023, 0.76229608, 0.74945305, 0.82121212, 0.82669968],
+        [-0.7758696, 0.13295831, 0.87360834, -1.21879122, 0.62980314],
+    ),
+    ("CPPStandaloneDevice", "cython", 1): (
+        [0.1636023, 0.76229608, 0.74945305, 0.82121212, 0.82669968],
+        [-0.7758696, 0.13295831, 0.87360834, -1.21879122, 0.62980314],
+    ),
+    ("CPPStandaloneDevice", None, 4): (
+        [0.1636023, 0.76229608, 0.27318909, 0.44124824, 0.69454226],
+        [0.36643979, -1.53883951, 0.07274151, 1.34278769, 0.63249739],
+    ),
+}
+
+
+def _config_tuple():
+    config = [
+        get_device().__class__.__name__,
+        prefs.codegen.target,
+        prefs.devices.cpp_standalone.openmp_threads,
+    ]
+    if config[0] == "RuntimeDevice":
+        config[2] = None
+    else:
+        config[1] = None
+    return tuple(config)
+
+
+@pytest.mark.standalone_compatible
+def test_random_values_fixed_seed_numbers():
+    # Verify a subset of random numbers, to make sure these numbers stay the same across updates
+    G = NeuronGroup(
+        100,
+        """
+        v1 : 1
+        v2 : 1
+        """,
+    )
+    seed(9876)
+    G.v1 = "rand()"
+    G.v2 = "randn()"
+    run(0 * ms)  # for standalone
+    expected_values = _random_values.get(_config_tuple(), None)
+    if expected_values is None:
+        pytest.skip("Random values not known for this configuration")
+    assert_allclose(G.v1[::20], expected_values[0])
+    assert_allclose(G.v2[::20], expected_values[1])
+
+
 @pytest.mark.standalone_compatible
 @pytest.mark.multiple_runs
 def test_random_values_fixed_and_random():