Merge 45b9352 into 1a6e5fe

.gitignore

@@ -1,3 +1,5 @@
+/dev/ideas/randomkit/CUBA
+/brian2/random/randomkit/randkit.c
 /dev/benchmarks/compare_to_brian1.png
 # Compiled Python files
 *.py[cod]

brian2/codegen/generators/cython_generator.py

@@ -316,20 +316,35 @@ def determine_keywords(self):
 
 
 rand_code = '''
-cdef int _rand_buffer_size = 1024 
-cdef double[:] _rand_buf = _numpy.zeros(_rand_buffer_size, dtype=_numpy.float64)
-cdef int _cur_rand_buf = 0
-cdef double _rand(int _idx):
-    global _cur_rand_buf
-    global _rand_buf
-    if _cur_rand_buf==0:
-        _rand_buf = _numpy.random.rand(_rand_buffer_size)
-    cdef double val = _rand_buf[_cur_rand_buf]
-    _cur_rand_buf = (_cur_rand_buf+1)%_rand_buffer_size
-    return val
+from libc.stdlib cimport malloc
+
+cdef extern from "randomkit.h":
+    ctypedef struct rk_state:
+        unsigned long key[624]
+        int pos
+        int has_gauss
+        double gauss
+    ctypedef enum rk_error:
+        RK_NOERR = 0
+        RK_ENODEV = 1
+        RK_ERR_MAX = 2
+    rk_error rk_randomseed(rk_state *state)
+    double rk_double(rk_state *state)
+
+cdef rk_state* _mtrandstate_rand = NULL
+
+cdef inline double _rand(int _idx):
+    global _mtrandstate_rand
+    cdef rk_error errcode
+    if _mtrandstate_rand==NULL:
+        _mtrandstate_rand = <rk_state *>malloc(sizeof(rk_state))
+        errcode = rk_randomseed(_mtrandstate_rand)
+        if errcode:
+            raise RuntimeError("Cannot initialise random state.")
+    return rk_double(_mtrandstate_rand)
 '''
 
-randn_code = rand_code.replace('rand', 'randn').replace('randnom', 'random')
+randn_code = rand_code.replace('rand', 'randn').replace('randnom', 'random').replace('rk_double', 'rk_gauss')
 
 DEFAULT_FUNCTIONS['rand'].implementations.add_implementation(CythonCodeGenerator,
                                                              code=rand_code,

brian2/codegen/runtime/cython_rt/extension_manager.py

@@ -30,6 +30,7 @@
 except ImportError:
     Cython = None
 
+import brian2
 from brian2.utils.logger import std_silent, get_logger
 from brian2.utils.stringtools import deindent
 from brian2.core.preferences import prefs
@@ -174,22 +175,30 @@ def _load_module(self, module_path, include_dirs, library_dirs,
                 import numpy
                 c_include_dirs.append(numpy.get_include())
 
-            # TODO: We should probably have a special folder just for header
-            # files that are shared between different codegen targets
             import brian2.synapses as synapses
-            synapses_dir = os.path.dirname(synapses.__file__)
+            synapses_dir = os.path.abspath(os.path.dirname(synapses.__file__))
             c_include_dirs.append(synapses_dir)
 
+            brian2dir, _ = os.path.split(brian2.__file__)
+            rkdir = os.path.abspath(os.path.join(brian2dir, 'utils',
+                                                 'random', 'randomkit'))
+            randomkitc = os.path.join(rkdir, 'randomkit.c')
+            c_include_dirs.append(rkdir)
+
+
             pyx_file = os.path.join(lib_dir, module_name + '.pyx')
             # ignore Python 3 unicode stuff for the moment
             #pyx_file = py3compat.cast_bytes_py2(pyx_file, encoding=sys.getfilesystemencoding())
             #with io.open(pyx_file, 'w') as f:#, encoding='utf-8') as f:
             #    f.write(code)
             open(pyx_file, 'w').write(code)
 
+            if sys.platform=='win32':
+                libraries = libraries+['advapi32'] # needed for randomkit
+
             extension = Extension(
                 name=module_name,
-                sources=[pyx_file],
+                sources=[pyx_file, randomkitc],
                 include_dirs=c_include_dirs,
                 library_dirs=library_dirs,
                 runtime_library_dirs=runtime_library_dirs,

brian2/codegen/runtime/weave_rt/weave_rt.py

@@ -18,6 +18,7 @@
         # No weave for Python 3
         weave = None
 
+import brian2
 from brian2.core.variables import (DynamicArrayVariable, ArrayVariable,
                                    AuxiliaryVariable, Subexpression)
 from brian2.core.preferences import prefs
@@ -101,12 +102,23 @@ def __init__(self, owner, code, variables, variable_indices,
         # TODO: We should probably have a special folder just for header
         # files that are shared between different codegen targets
         import brian2.synapses as synapses
-        synapses_dir = os.path.dirname(synapses.__file__)
+        synapses_dir = os.path.abspath(os.path.dirname(synapses.__file__))
+        brian2dir, _ = os.path.split(brian2.__file__)
+        rkdir = os.path.abspath(os.path.join(brian2dir, 'utils',
+                                             'random', 'randomkit'))
+        randomkitc = os.path.join(rkdir, 'randomkit.c')
+
         self.include_dirs.append(synapses_dir)
+        self.include_dirs.append(rkdir)
         self.library_dirs = list(prefs['codegen.cpp.library_dirs'])
         self.runtime_library_dirs = list(prefs['codegen.cpp.runtime_library_dirs'])
         self.libraries = list(prefs['codegen.cpp.libraries'])
-        self.headers = ['<algorithm>', '<limits>', '"stdint_compat.h"'] + prefs['codegen.cpp.headers']
+        if sys.platform == 'win32':
+            self.libraries.append('advapi32')  # needed for randomkit
+        self.headers = (['<algorithm>', '<limits>',
+                         '"stdint_compat.h"', '"randomkit.h"'] +
+                        prefs['codegen.cpp.headers'])
+        self.sources = [randomkitc]
         self.annotated_code = self.code.main+'''
 /*
 The following code is just compiler options for the call to weave.inline.
@@ -249,6 +261,7 @@ def run(self):
                 extra_link_args=self.extra_link_args,
                 include_dirs=self.include_dirs,
                 library_dirs=self.library_dirs,
+                sources=self.sources,
                 verbose=0)
             with std_silent():
                 ret_val = weave.inline(*self._inline_args, **self._inline_kwds)
@@ -262,85 +275,46 @@ def run(self):
     codegen_targets.add(WeaveCodeObject)
 
 
-# Use a special implementation for the randn function that makes use of numpy's
-# randn
+# Use RandomKit for random number generation (same algorithm as numpy)
+# Note that we create a new random state for each codeobject, but this
+# is seeded with a few hundred bytes from OS-level entropy generators
+# and so the random numbers produced will be independent. If these
+# RNGs aren't able to produce enough data to seed the RNG it will raise
+# an error.
 randn_code = {'support_code': '''
-        #define BUFFER_SIZE 1024
-        // A randn() function that returns a single random number. Internally
-        // it asks numpy's randn function for BUFFER_SIZE
-        // random numbers at a time and then returns one number from this
-        // buffer.
-        // It needs a reference to the numpy_randn object (the original numpy
-        // function), because this is otherwise only available in
-        // compiled_function (where is is automatically handled by weave).
-        //
-        double _randn(const int _vectorisation_idx) {
-            // the _vectorisation_idx argument is unused for now, it could in
-            // principle be used to get reproducible random numbers when using
-            // OpenMP etc.
-            static PyArrayObject *randn_buffer = NULL;
-            static double *buf_pointer = NULL;
-            static npy_int curbuffer = 0;
-            if(curbuffer==0)
-            {
-                if(randn_buffer) Py_DECREF(randn_buffer);
-                py::tuple args(1);
-                args[0] = BUFFER_SIZE;
-                randn_buffer = (PyArrayObject *)PyArray_FromAny(_namespace_numpy_randn.call(args),
-                                                                NULL, 1, 1, 0, NULL);
-                buf_pointer = (double*)PyArray_GETPTR1(randn_buffer, 0);
+        rk_state *_mtrandstate_randn = NULL;
+        inline double _randn(const int _vectorisation_idx) {
+            if(_mtrandstate_randn==NULL) {
+                _mtrandstate_randn = new rk_state;
+                rk_error errcode = rk_randomseed(_mtrandstate_randn);
+                if(errcode)
+                {
+                    PyErr_SetString(PyExc_RuntimeError, "Cannot initialise random state");
+                    throw 1;
+                }
             }
-            double number = buf_pointer[curbuffer];
-            curbuffer = curbuffer+1;
-            if (curbuffer == BUFFER_SIZE)
-                // This seems to be safer then using (curbuffer + 1) % BUFFER_SIZE, we might run into
-                // an integer overflow for big networks, otherwise.
-                curbuffer = 0;
-            return number;
+            return rk_gauss(_mtrandstate_randn);
         }
         '''}
 DEFAULT_FUNCTIONS['randn'].implementations.add_implementation(WeaveCodeObject,
                                                               code=randn_code,
-                                                              name='_randn',
-                                                              namespace={'_numpy_randn': numpy.random.randn})
+                                                              name='_randn')
 
-# Also use numpy for rand
 rand_code = {'support_code': '''
-        #define BUFFER_SIZE 1024
-        // A rand() function that returns a single random number. Internally
-        // it asks numpy's rand function for BUFFER_SIZE
-        // random numbers at a time and then returns one number from this
-        // buffer.
-        // It needs a reference to the numpy_rand object (the original numpy
-        // function), because this is otherwise only available in
-        // compiled_function (where is is automatically handled by weave).
-        //
-        double _rand(const int _vectorisation_idx) {
-            // the _vectorisation_idx argument is unused for now, it could in
-            // principle be used to get reproducible random numbers when using
-            // OpenMP etc.
-            static PyArrayObject *rand_buffer = NULL;
-            static double *buf_pointer = NULL;
-            static npy_int curbuffer = 0;
-            if(curbuffer==0)
-            {
-                if(rand_buffer) Py_DECREF(rand_buffer);
-                py::tuple args(1);
-                args[0] = BUFFER_SIZE;
-                rand_buffer = (PyArrayObject *)PyArray_FromAny(_namespace_numpy_rand.call(args),
-                                                               NULL, 1, 1, 0, NULL);
-                buf_pointer = (double*)PyArray_GETPTR1(rand_buffer, 0);
+        rk_state *_mtrandstate_rand = NULL;
+        inline double _rand(const int _vectorisation_idx) {
+            if(_mtrandstate_rand==NULL) {
+                _mtrandstate_rand = new rk_state;
+                rk_error errcode = rk_randomseed(_mtrandstate_rand);
+                if(errcode)
+                {
+                    PyErr_SetString(PyExc_RuntimeError, "Cannot initialise random state");
+                    throw 1;
+                }
             }
-            double number = buf_pointer[curbuffer];
-            curbuffer = curbuffer+1;
-            if (curbuffer == BUFFER_SIZE)
-                // This seems to be safer then using (curbuffer + 1) % BUFFER_SIZE, we might run into
-                // an integer overflow for big networks, otherwise.
-                curbuffer = 0;
-            return number;
+            return rk_double(_mtrandstate_rand);
         }
         '''}
 DEFAULT_FUNCTIONS['rand'].implementations.add_implementation(WeaveCodeObject,
                                                              code=rand_code,
-                                                             namespace={'_numpy_rand': numpy.random.rand},
                                                              name='_rand')

brian2/devices/cpp_standalone/codeobject.py

@@ -8,6 +8,9 @@
                                                      c_data_type)
 from brian2.devices.device import get_device
 from brian2.core.preferences import prefs
+from brian2.core.functions import DEFAULT_FUNCTIONS, Function
+
+import numpy
 
 __all__ = ['CPPStandaloneCodeObject']
 
@@ -93,7 +96,8 @@ class CPPStandaloneCodeObject(CodeObject):
                           env_globals={'c_data_type': c_data_type,
                                        'openmp_pragma': openmp_pragma,
                                        'constant_or_scalar': constant_or_scalar,
-                                       'prefs': prefs})
+                                       'prefs': prefs,
+                                       'zip': zip})
     generator_class = CPPCodeGenerator
 
     def __call__(self, **kwds):
@@ -103,3 +107,21 @@ def run(self):
         get_device().main_queue.append(('run_code_object', (self,)))
 
 codegen_targets.add(CPPStandaloneCodeObject)
+
+rand_code = {'support_code': '''
+        inline double _rand(const int _vectorisation_idx) {
+            return rk_double(&brian::_mersenne_twister_state);
+        }
+        '''}
+DEFAULT_FUNCTIONS['rand'].implementations.add_implementation(CPPStandaloneCodeObject,
+                                                             code=rand_code,
+                                                             name='_rand')
+
+randn_code = {'support_code': '''
+        inline double _randn(const int _vectorisation_idx) {
+            return rk_gauss(&brian::_mersenne_twister_state);
+        }
+        '''}
+DEFAULT_FUNCTIONS['randn'].implementations.add_implementation(CPPStandaloneCodeObject,
+                                                              code=randn_code,
+                                                              name='_randn')

brian2/devices/cpp_standalone/device.py

@@ -4,6 +4,7 @@
 import os
 import shutil
 import subprocess
+import sys
 import inspect
 import platform
 from collections import defaultdict, Counter
@@ -14,6 +15,8 @@
 import numpy as np
 from cpuinfo import cpuinfo
 
+import brian2
+
 from brian2.codegen.cpp_prefs import get_compiler_and_args
 from brian2.core.network import Network
 from brian2.devices.device import Device, all_devices, set_device, reset_device
@@ -77,7 +80,7 @@ def __init__(self, project_dir):
 
     def write(self, filename, contents):
         logger.diagnostic('Writing file %s:\n%s' % (filename, contents))
-        if filename.lower().endswith('.cpp'):
+        if filename.lower().endswith('.cpp') or filename.lower().endswith('.c'):
             self.source_files.append(filename)
         elif filename.lower().endswith('.h'):
             self.header_files.append(filename)
@@ -668,9 +671,10 @@ def generate_makefile(self, writer, compiler, compiler_flags, linker_flags, nb_t
             else:
                 openmp_flag = ''
             # Generate the visual studio makefile
-            source_bases = [fname.replace('.cpp', '').replace('/', '\\') for fname in writer.source_files]
+            source_bases = [fname.replace('.cpp', '').replace('.c', '').replace('/', '\\') for fname in writer.source_files]
             win_makefile_tmp = CPPStandaloneCodeObject.templater.win_makefile(
                 None, None,
+                source_files=writer.source_files,
                 source_bases=source_bases,
                 compiler_flags=compiler_flags,
                 linker_flags=linker_flags,
@@ -707,7 +711,17 @@ def copy_source_files(self, writer, directory):
                      os.path.join(directory, 'brianlib', 'spikequeue.h'))
         shutil.copy2(os.path.join(os.path.split(inspect.getsourcefile(Synapses))[0], 'stdint_compat.h'),
                      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],
+                                  'utils', 'random', 'randomkit', 'randomkit.c'),
+                     os.path.join(directory, 'brianlib', 'randomkit', 'randomkit.c'))
+        shutil.copy2(os.path.join(os.path.split(inspect.getsourcefile(brian2))[0],
+                                  'utils', 'random', 'randomkit', 'randomkit.h'),
+                     os.path.join(directory, 'brianlib', 'randomkit', 'randomkit.h'))
+
     def write_static_arrays(self, directory):
         # # Find np arrays in the namespaces and convert them into static
         # # arrays. Hopefully they are correctly used in the code: For example,
@@ -891,14 +905,20 @@ def build(self, directory='output',
         compiler, extra_compile_args = get_compiler_and_args()
         compiler_obj = ccompiler.new_compiler(compiler=compiler)
         compiler_flags = (ccompiler.gen_preprocess_options(prefs['codegen.cpp.define_macros'],
-                                                           prefs['codegen.cpp.include_dirs']) +
+                                                           prefs['codegen.cpp.include_dirs']+['brianlib/randomkit']) +
                           extra_compile_args)
+        if sys.platform=='win32':
+            wincrypt = ['advapi32']
+        else:
+            wincrypt = []
         linker_flags = (ccompiler.gen_lib_options(compiler_obj,
-                                                  library_dirs=prefs['codegen.cpp.library_dirs'],
+                                                  library_dirs=prefs['codegen.cpp.library_dirs']+['brianlib/randomkit'],
                                                   runtime_library_dirs=prefs['codegen.cpp.runtime_library_dirs'],
-                                                  libraries=prefs['codegen.cpp.libraries']) +
+                                                  libraries=prefs['codegen.cpp.libraries']+wincrypt) +
                         prefs['codegen.cpp.extra_link_args'])
 
+        additional_source_files.append('brianlib/randomkit/randomkit.c')
+
         for d in ['code_objects', 'results', 'static_arrays']:
             ensure_directory(os.path.join(directory, d))
 

brian2/devices/cpp_standalone/templates/main.cpp

@@ -5,6 +5,7 @@
 {{ openmp_pragma('include') }}
 #include "run.h"
 #include "brianlib/common_math.h"
+#include "randomkit.h"
 
 {% for codeobj in code_objects %}
 #include "code_objects/{{codeobj.name}}.h"

brian2/devices/cpp_standalone/templates/makefile

@@ -2,6 +2,7 @@ PROGRAM = main
 SRCS = {{source_files}}
 H_SRCS = {{header_files}}
 OBJS = ${SRCS:.cpp=.o}
+OBJS := ${OBJS:.c=.o}
 CC = @g++
 DEBUG = -g
 OPTIMISATIONS = {{ compiler_flags }}