Replace randomkit by standard library implementation

Closes #1483

LICENSE

@@ -661,77 +661,6 @@ THE SOFTWARE.
 
 -------------------------------------------------------------------------------
 
-The Mersenne-Twister implementation "Randomkit" in brian2.random.randomkit is
-based on code by Jean-Sebastien Roy, provided under the MIT license:
-
-Copyright (c) 2003-2005, Jean-Sebastien Roy ([email protected])
-
-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.
-
---
-
-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.
-
---
-
 The floor division and modulo implementations for integer numbers in C++ standalone
 are based on the implementation in Cython (https://github.com/cython/cython), published
 under the Apache License 2.0:

brian2/devices/cpp_standalone/codeobject.py

@@ -139,14 +139,14 @@ def generate_rand_code(rand_func, owner):
     else:
         thread_number = "omp_get_thread_num()"
     if rand_func == "rand":
-        rk_call = "rk_double"
+        rk_call = "_uniform_random"
     elif rand_func == "randn":
-        rk_call = "rk_gauss"
+        rk_call = "_normal_random"
     else:
         raise AssertionError(rand_func)
     code = """
            double _%RAND_FUNC%(const int _vectorisation_idx) {
-               return %RK_CALL%(brian::_mersenne_twister_states[%THREAD_NUMBER%]);
+               return brian::%RK_CALL%(brian::_mersenne_twister_generators[%THREAD_NUMBER%]);
            }
            """
     code = replace(

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"):
@@ -874,12 +873,11 @@ def generate_main_source(self, writer):
                 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]);"
+                        "    brian::_mersenne_twister_generators[_i] = std::mt19937(_rd());"
                     )
                 else:
                     main_lines.append(
-                        f"    rk_seed({seed!r}L + _i,"
-                        " brian::_mersenne_twister_states[_i]);"
+                        f"brian::_mersenne_twister_generators[_i].seed({seed!r}L + _i);"
                     )
             else:
                 raise NotImplementedError(f"Unknown main queue function type {func}")
@@ -1084,28 +1082,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 +1559,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 +1677,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 +1704,6 @@ def delete(self, code=True, data=True, directory=True, force=False):
 
         if directory:
             directories = [
-                os.path.join("brianlib", "randomkit"),
                 "brianlib",
                 "code_objects",
                 "results",

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")
 	{

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,9 @@ 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;
+std::vector< std::mt19937 > _mersenne_twister_generators;
+std::uniform_real_distribution<double> _uniform_random;
+std::normal_distribution<double> _normal_random;
 
 //////////////// networks /////////////////
 {% for net in networks | sort(attribute='name') %}
@@ -223,8 +225,9 @@ void _init_arrays()
 	{% endfor %}
 
 	// Random number generator states
+	std::random_device rd;
 	for (int i=0; i<{{openmp_pragma('get_num_threads')}}; i++)
-	    _mersenne_twister_states.push_back(new rk_state());
+	    _mersenne_twister_generators.push_back(std::mt19937(rd()));
 }
 
 void _load_arrays()
@@ -369,15 +372,17 @@ 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;
 // In OpenMP we need one state per thread
-extern std::vector< rk_state* > _mersenne_twister_states;
+extern std::vector< std::mt19937 > _mersenne_twister_generators;
+extern std::uniform_real_distribution<double> _uniform_random;
+extern std::normal_distribution<double> _normal_random;
 
 //////////////// clocks ///////////////////
 {% for clock in clocks | sort(attribute='name') %}

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()