From f336c7154ea11f3093db982f26edd30fcaaf83ee Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Fri, 12 Jul 2013 14:40:39 +0200 Subject: [PATCH 01/15] Switch SpikeMonitor implementation to code generation --- brian2/codegen/languages/base.py | 4 ++ .../languages/cpp/templates/monitor.cpp | 22 ------- .../languages/cpp/templates/spikemonitor.cpp | 31 +++++++++ .../languages/python/templates/monitor.py_ | 6 -- .../python/templates/spikemonitor.py_ | 15 +++++ brian2/monitors/spikemonitor.py | 64 +++++++++++++------ brian2/tests/test_monitor.py | 36 +++++++++++ 7 files changed, 130 insertions(+), 48 deletions(-) delete mode 100644 brian2/codegen/languages/cpp/templates/monitor.cpp create mode 100644 brian2/codegen/languages/cpp/templates/spikemonitor.cpp delete mode 100644 brian2/codegen/languages/python/templates/monitor.py_ create mode 100644 brian2/codegen/languages/python/templates/spikemonitor.py_ create mode 100644 brian2/tests/test_monitor.py diff --git a/brian2/codegen/languages/base.py b/brian2/codegen/languages/base.py index 34447a04a..ec6e74742 100644 --- a/brian2/codegen/languages/base.py +++ b/brian2/codegen/languages/base.py @@ -148,6 +148,10 @@ def template_reset(self): def template_threshold(self): return self.templater.threshold + @property + def template_spikemonitor(self): + return self.templater.spikemonitor + @property def template_synapses(self): return self.templater.synapses diff --git a/brian2/codegen/languages/cpp/templates/monitor.cpp b/brian2/codegen/languages/cpp/templates/monitor.cpp deleted file mode 100644 index 3eee1d47c..000000000 --- a/brian2/codegen/languages/cpp/templates/monitor.cpp +++ /dev/null @@ -1,22 +0,0 @@ -{% macro main() %} - // Get the current length and new length of t and i arrays - int _curlen = t_arr.attr("shape")[0]; - int _newlen = _curlen+_num_spikes; - // Resize the arrays - py::tuple _newlen_tuple(1); - _newlen_tuple[0] = _newlen; - t_arr.mcall("resize", _newlen_tuple); - i_arr.mcall("resize", _newlen_tuple); - // Get the potentially newly created underlying data arrays - double *_t_arr_data = (double*)(((PyArrayObject*)(PyObject*)t_arr.attr("data"))->data); - int *_i_arr_data = (int*)(((PyArrayObject*)(PyObject*)i_arr.attr("data"))->data); - // Copy the values across - for(int _i=0; _i<_num_spikes; _i++) - { - _t_arr_data[_curlen+_i] = t; - _i_arr_data[_curlen+_i] = _spikes[_i]; - } -{% endmacro %} - -{% macro support_code() %} -{% endmacro %} diff --git a/brian2/codegen/languages/cpp/templates/spikemonitor.cpp b/brian2/codegen/languages/cpp/templates/spikemonitor.cpp new file mode 100644 index 000000000..f54c40478 --- /dev/null +++ b/brian2/codegen/languages/cpp/templates/spikemonitor.cpp @@ -0,0 +1,31 @@ +{% macro main() %} + + // USE_SPECIFIERS { _t, _i, t, _spikes, _count } + + if (_num_spikes > 0) + { + // Get the current length and new length of t and i arrays + const int _curlen = _t.attr("shape")[0]; + const int _newlen = _curlen + _num_spikes; + // Resize the arrays + py::tuple _newlen_tuple(1); + _newlen_tuple[0] = _newlen; + _t.mcall("resize", _newlen_tuple); + _i.mcall("resize", _newlen_tuple); + // Get the potentially newly created underlying data arrays + double *_t_data = (double*)(((PyArrayObject*)(PyObject*)_t.attr("data"))->data); + // TODO: How to get the correct datatype automatically here? + npy_int64 *_i_data = (npy_int64*)(((PyArrayObject*)(PyObject*)_i.attr("data"))->data); + // Copy the values across + for(int _idx=0; _idx<_num_spikes; _idx++) + { + const int _neuron_idx = _spikes[_idx]; + _t_data[_curlen + _idx] = t; + _i_data[_curlen + _idx] = _neuron_idx; + _count[_neuron_idx]++; + } + } +{% endmacro %} + +{% macro support_code() %} +{% endmacro %} diff --git a/brian2/codegen/languages/python/templates/monitor.py_ b/brian2/codegen/languages/python/templates/monitor.py_ deleted file mode 100644 index 97d1ff745..000000000 --- a/brian2/codegen/languages/python/templates/monitor.py_ +++ /dev/null @@ -1,6 +0,0 @@ -_curlen = len(t_arr) -_newlen = len(t_arr)+len(_spikes) -t_arr.resize(_newlen) -i_arr.resize(_newlen) -t_arr[_curlen:_newlen] = t -i_arr[_curlen:_newlen] = _spikes diff --git a/brian2/codegen/languages/python/templates/spikemonitor.py_ b/brian2/codegen/languages/python/templates/spikemonitor.py_ new file mode 100644 index 000000000..35d125755 --- /dev/null +++ b/brian2/codegen/languages/python/templates/spikemonitor.py_ @@ -0,0 +1,15 @@ +# { USE_SPECIFIERS _i, _t, _spikes, _count, _num_source_neurons, t } + +_n_spikes = len(_spikes) +if _n_spikes > 0: + import numpy as np + + _curlen = len(_t) + _newlen = _curlen + _n_spikes + _t.resize(_newlen) + _i.resize(_newlen) + _t[_curlen:_newlen] = t + _i[_curlen:_newlen] = _spikes + + # This is slow but correctly handles multiple spikes per neuron + _count += np.bincount(_spikes, minlength=_num_source_neurons); \ No newline at end of file diff --git a/brian2/monitors/spikemonitor.py b/brian2/monitors/spikemonitor.py index 1f799c007..59d4c8909 100644 --- a/brian2/monitors/spikemonitor.py +++ b/brian2/monitors/spikemonitor.py @@ -1,15 +1,19 @@ import weakref -from numpy import zeros +import numpy as np from brian2.core.base import BrianObject from brian2.core.preferences import brian_prefs from brian2.core.scheduler import Scheduler +from brian2.core.specifiers import ReadOnlyValue, AttributeValue, ArrayVariable +from brian2.codegen.languages.python import PythonLanguage from brian2.memory.dynamicarray import DynamicArray1D from brian2.units.allunits import second +from brian2.units.fundamentalunits import Unit __all__ = ['SpikeMonitor'] + class SpikeMonitor(BrianObject): ''' Record spikes from a `NeuronGroup` or other spike source @@ -28,10 +32,15 @@ class SpikeMonitor(BrianObject): A unique name for the object, otherwise will use ``source.name+'_spikemonitor_0'``, etc. ''' - def __init__(self, source, record=True, when=None, name='spikemonitor*'): + def __init__(self, source, record=True, when=None, name='spikemonitor*', + language=None): self.source = weakref.proxy(source) self.record = bool(record) + if language is None: + language = PythonLanguage() + self.language = language + # run by default on source clock at the end scheduler = Scheduler(when) if not scheduler.defined_clock: @@ -42,7 +51,27 @@ def __init__(self, source, record=True, when=None, name='spikemonitor*'): # create data structures self.reinit() - + + self.specifiers = {'t': AttributeValue('t', second, np.float64, + self.clock, 't'), + '_spikes': AttributeValue('_spikes', Unit(1), + self.source.spikes.dtype, + self.source, 'spikes'), + # The template needs to have access to the + # DynamicArray here, having access to the underlying + # array is not enough since we want to do the resize + # in the template + '_i': ReadOnlyValue('_i', Unit(1), self._i.dtype, + self._i), + '_t': ReadOnlyValue('_t', Unit(1), self._t.dtype, + self._t), + '_count': ArrayVariable('_count', Unit(1), + self.count.dtype, + self.count, ''), + '_num_source_neurons': ReadOnlyValue('_num_source_neurons', + Unit(1), np.int, + len(self.source))} + def reinit(self): ''' Clears all recorded spikes @@ -52,24 +81,19 @@ def reinit(self): dtype=brian_prefs['core.default_scalar_dtype']) #: Array of the number of times each source neuron has spiked - self.count = zeros(len(self.source), dtype=int) - + self.count = np.zeros(len(self.source), dtype=int) + + def pre_run(self, namespace): + template = self.language.template_spikemonitor + self.codeobj = self.language.create_codeobj(self.name, + '', # No model-specific code + {}, # no namespace + self.specifiers, + template, + indices={}) + def update(self): - spikes = self.source.spikes - nspikes = len(spikes) - if nspikes: - if self.record: - # update i, t arrays - i = self._i - t = self._t - oldsize = len(i) - newsize = oldsize+nspikes - i.resize(newsize) - t.resize(newsize) - i.data[oldsize:] = spikes - t.data[oldsize:] = self.clock.t_ - # update count - self.count[spikes] += 1 + self.codeobj() @property def i(self): diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py new file mode 100644 index 000000000..14feab65d --- /dev/null +++ b/brian2/tests/test_monitor.py @@ -0,0 +1,36 @@ +import numpy as np +from numpy.testing.utils import assert_allclose, assert_equal + +from brian2 import * + +# We can only test C++ if weave is availabe +try: + import scipy.weave + languages = [PythonLanguage(), CPPLanguage()] +except ImportError: + # Can't test C++ + languages = [PythonLanguage()] + + +def test_spike_monitor(): + for language in languages: + defaultclock.t = 0*second + G = NeuronGroup(2, '''dv/dt = rate : 1 + rate: Hz''', threshold='v>1', reset='v=0', + language=language) + # We don't use 100 and 1000Hz, because then the membrane potential would + # be exactly at 1 after 10 resp. 100 timesteps. Due to floating point + # issues this will not be exact, + G.rate = [101, 1001] * Hz + + mon = SpikeMonitor(G, language=language) + net = Network(G, mon) + net.run(10*ms) + + assert_allclose(mon.t[mon.i == 0], [9.9]*ms) + assert_allclose(mon.t[mon.i == 1], np.arange(10)*ms + 0.9*ms) + assert_equal(mon.count, np.array([1, 10])) + +if __name__ == '__main__': + test_spike_monitor() + From 551d9ed614b9b0947b1027703bff2122438fc6c6 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Mon, 22 Jul 2013 17:21:32 +0200 Subject: [PATCH 02/15] Start implementing StateMonitor with code generation (not working yet) --- brian2/codegen/languages/base.py | 4 + .../python/templates/statemonitor.py_ | 23 ++++ brian2/monitors/statemonitor.py | 119 +++++++++--------- 3 files changed, 88 insertions(+), 58 deletions(-) create mode 100644 brian2/codegen/languages/python/templates/statemonitor.py_ diff --git a/brian2/codegen/languages/base.py b/brian2/codegen/languages/base.py index ec6e74742..e271e94cb 100644 --- a/brian2/codegen/languages/base.py +++ b/brian2/codegen/languages/base.py @@ -152,6 +152,10 @@ def template_threshold(self): def template_spikemonitor(self): return self.templater.spikemonitor + @property + def template_statemonitor(self): + return self.templater.statemonitor + @property def template_synapses(self): return self.templater.synapses diff --git a/brian2/codegen/languages/python/templates/statemonitor.py_ b/brian2/codegen/languages/python/templates/statemonitor.py_ new file mode 100644 index 000000000..5cff06134 --- /dev/null +++ b/brian2/codegen/languages/python/templates/statemonitor.py_ @@ -0,0 +1,23 @@ +# USE_SPECIFIERS { _t, _values, t, _indices } + +# Resize dynamic arrays +new_len = len(_t) + 1 +_t.resize(new_len) + +{% for _varname in _variable_names %} +_recorded_{{_varname}}.resize((_new_len, _num_indices)) +{% endfor %} + +# Store values +_t[new_len] = t + +_vectorisation_idx = _indices +_record_idx = _indices +{% for line in code_lines %} +{{line}} +{% endfor %} + +{% for _varname in _variable_names %} +_recorded_{{_varname}}[-1, :] = _to_record_{{_varname}} +{% endfor %} + diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index 5ede29a63..c84f18edc 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -1,22 +1,28 @@ import weakref -from numpy import array, arange +import numpy as np -from brian2.core.specifiers import Value +from brian2.core.specifiers import (Value, DynamicArrayVariable, + ReadOnlyValue, ArrayVariable, + Index) from brian2.core.base import BrianObject from brian2.core.scheduler import Scheduler +from brian2.codegen.languages.python import PythonLanguage from brian2.groups.group import Group +from brian2.units.fundamentalunits import Unit from brian2.units.allunits import second +from brian2.memory.dynamicarray import DynamicArray, DynamicArray1D __all__ = ['StateMonitor'] + class MonitorVariable(Value): def __init__(self, name, unit, dtype, monitor): Value.__init__(self, name, unit, dtype) self.monitor = weakref.proxy(monitor) def get_value(self): - return array(self.monitor._values[self.name]) + return self.monitor._values[self.name] class StateMonitor(BrianObject, Group): @@ -62,19 +68,16 @@ class StateMonitor(BrianObject, Group): run(100*ms) plot(M.t, M.V) show() - - Notes - ----- - TODO: multiple features, below: - - * Cacheing extracted values (t, V, etc.) - * Improve efficiency by using dynamic arrays instead of lists? ''' def __init__(self, source, variables, record=None, when=None, - name='statemonitor*'): + name='statemonitor*', language=None): self.source = weakref.proxy(source) + if language is None: + language = PythonLanguage() + self.language = language + # run by default on source clock at the end scheduler = Scheduler(when) if not scheduler.defined_clock: @@ -91,12 +94,14 @@ def __init__(self, source, variables, record=None, when=None, self.variables = variables # record should always be an array of ints + self.record_all = False if record is None or record is False: - record = array([], dtype=int) + record = np.array([], dtype=int) elif record is True: - record = arange(len(source)) + self.record_all = True + record = np.arange(len(source)) else: - record = array(record, dtype=int) + record = np.array(record, dtype=int) #: The array of recorded indices self.indices = record @@ -106,61 +111,59 @@ def __init__(self, source, variables, record=None, when=None, # initialise Group access self.specifiers = {} - for variable in variables: + for idx, variable in enumerate(variables): spec = source.specifiers[variable] + self.specifiers['_source_' + variable] = ArrayVariable('_source_'+variable, + spec.unit, + spec.dtype, + spec.array, + '_record_idx', + group=spec.group, + constant=spec.constant, + scalar=spec.scalar, + is_bool=spec.is_bool) self.specifiers[variable] = MonitorVariable(variable, spec.unit, spec.dtype, - self) + self) + self.specifiers['_recorded_'+variable] = ReadOnlyValue('_recorded_'+variable, Unit(1), + self._values[variable].dtype, + self._values[variable]) + self.specifiers['t'] = DynamicArrayVariable('t', second, dtype=np.float64, + array=self._t, index='', + group=self) + + self.specifiers['_t'] = ReadOnlyValue('_t', Unit(1), self._t.dtype, + self._t), + Group.__init__(self) - + + def reinit(self): - self._values = dict((var, []) for var in self.variables) - self._t = [] + self._values = dict((v, DynamicArray((0, len(self.variables)))) + for v in self.variables) + self._t = DynamicArray1D(0) + def pre_run(self, namespace): + template = self.language.template_statemonitor + # Some dummy code so that code generation takes care of the indexing + # and subexpressions + code = ['_to_record_%s = _source_%s' % (v, v) + for v in self.variables] + code = '\n'.join(code) + self.codeobj = self.language.create_codeobj(self.name, + code, + {}, # no namespace + self.specifiers, + template, + indices={'_record_idx': Index('_record_idx', self.record_all)}, + template_kwds={'_variable_names': self.variables}) + def update(self): - for var in self.variables: - self._values[var].append(self.source.state_(var)[self.indices]) - self._t.append(self.clock.t_) - - @property - def t(self): - ''' - Array of record times. - ''' - return array(self._t)*second - - @property - def t_(self): - ''' - Array of record times (without units). - ''' - return array(self._t) + self.codeobj() def __repr__(self): description = '<{classname}, recording {variables} from {source}>' return description.format(classname=self.__class__.__name__, variables=repr(self.variables), source=self.source.name) - - -if __name__=='__main__': - from pylab import * - from brian2 import * - from brian2.codegen.languages import * - import time - - N = 100 - tau = 10*ms - eqs = ''' - dV/dt = (2*volt-V)/tau : volt - ''' - threshold = 'V>1' - reset = 'V = 0' - G = NeuronGroup(N, eqs, threshold=threshold, reset=reset) - G.V = rand(N)*volt - M = StateMonitor(G, True, record=range(5)) - run(100*ms) - print M.V.shape - plot(M.t, M.V) - show() From 7ee916bc6753fdb69ccc6872173e0bf32ca1148e Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 10:40:43 +0200 Subject: [PATCH 03/15] Python codegeneration working for StateMonitor --- .../python/templates/statemonitor.py_ | 8 +++---- brian2/monitors/statemonitor.py | 22 +++++++++++++------ 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/brian2/codegen/languages/python/templates/statemonitor.py_ b/brian2/codegen/languages/python/templates/statemonitor.py_ index 5cff06134..ab287297a 100644 --- a/brian2/codegen/languages/python/templates/statemonitor.py_ +++ b/brian2/codegen/languages/python/templates/statemonitor.py_ @@ -1,15 +1,15 @@ -# USE_SPECIFIERS { _t, _values, t, _indices } +# USE_SPECIFIERS { _t, _clock_t, _indices } # Resize dynamic arrays -new_len = len(_t) + 1 -_t.resize(new_len) +_new_len = len(_t) + 1 +_t.resize(_new_len) {% for _varname in _variable_names %} _recorded_{{_varname}}.resize((_new_len, _num_indices)) {% endfor %} # Store values -_t[new_len] = t +_t[-1] = _clock_t _vectorisation_idx = _indices _record_idx = _indices diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index c84f18edc..458ff1bc1 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -4,7 +4,7 @@ from brian2.core.specifiers import (Value, DynamicArrayVariable, ReadOnlyValue, ArrayVariable, - Index) + AttributeValue, Index) from brian2.core.base import BrianObject from brian2.core.scheduler import Scheduler from brian2.codegen.languages.python import PythonLanguage @@ -25,6 +25,14 @@ def get_value(self): return self.monitor._values[self.name] +class MonitorTime(Value): + def __init__(self, monitor): + Value.__init__(self, 't', second, np.float64) + self.monitor = weakref.proxy(monitor) + + def get_value(self): + return self.monitor._t[:] + class StateMonitor(BrianObject, Group): ''' Record values of state variables during a run @@ -129,15 +137,15 @@ def __init__(self, source, variables, record=None, when=None, self.specifiers['_recorded_'+variable] = ReadOnlyValue('_recorded_'+variable, Unit(1), self._values[variable].dtype, self._values[variable]) - self.specifiers['t'] = DynamicArrayVariable('t', second, dtype=np.float64, - array=self._t, index='', - group=self) - self.specifiers['_t'] = ReadOnlyValue('_t', Unit(1), self._t.dtype, - self._t), + self._t) - Group.__init__(self) + self.specifiers['_clock_t'] = AttributeValue('t', second, np.float64, + self.clock, 't_') + self.specifiers['t'] = MonitorTime(self) + + Group.__init__(self) def reinit(self): self._values = dict((v, DynamicArray((0, len(self.variables)))) From b4a61f6d932dc3d1596db12270b1a23bd16776a2 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 15:07:19 +0200 Subject: [PATCH 04/15] C++ codegeneration working for StateMonitor --- .../languages/cpp/templates/statemonitor.cpp | 62 +++++++++++++++++++ brian2/monitors/statemonitor.py | 2 +- 2 files changed, 63 insertions(+), 1 deletion(-) create mode 100644 brian2/codegen/languages/cpp/templates/statemonitor.cpp diff --git a/brian2/codegen/languages/cpp/templates/statemonitor.cpp b/brian2/codegen/languages/cpp/templates/statemonitor.cpp new file mode 100644 index 000000000..380edbf7b --- /dev/null +++ b/brian2/codegen/languages/cpp/templates/statemonitor.cpp @@ -0,0 +1,62 @@ +{% macro main() %} + + // USE_SPECIFIERS { _t, _clock_t, _indices } + + ////// SUPPORT CODE /// + {% for line in support_code_lines %} + //{{line}} + {% endfor %} + + ////// HANDLE DENORMALS /// + {% for line in denormals_code_lines %} + {{line}} + {% endfor %} + + ////// HASH DEFINES /////// + {% for line in hashdefine_lines %} + {{line}} + {% endfor %} + + ///// POINTERS //////////// + {% for line in pointers_lines %} + {{line}} + {% endfor %} + + // Get the current length and new length of t and value arrays + const int _curlen = _t.attr("shape")[0]; + const int _new_len = _curlen + 1; + // Resize the arrays + PyObject_CallMethod(_t, "resize", "i", _new_len); + {% for _varname in _variable_names %} + + PyObject_CallMethod(_recorded_{{_varname}}, "resize", "((ii))", + _new_len, _num_indices); + {% endfor %} + + // Get the potentially newly created underlying data arrays and copy the + // data + double *_t_data = (double*)(((PyArrayObject*)(PyObject*)_t.attr("data"))->data); + _t_data[_new_len - 1] = _clock_t; + + {% for _varname in _variable_names %} + PyArrayObject *_record_data = (((PyArrayObject*)(PyObject*)_recorded_{{_varname}}.attr("data"))); + const npy_intp* _record_strides = _record_data->strides; + for (int _idx=0; _idx < _num_indices; _idx++) + { + const int _record_idx = _indices[_idx]; + const int _vectorisation_idx = _record_idx; + {% for line in code_lines %} + {{line}} + {% endfor %} + double *recorded_entry = ((double*)(_record_data->data + (_new_len - 1)*_record_strides[0] + _idx*_record_strides[1])); + *recorded_entry = _to_record_{{_varname}}; + } + {% endfor %} + +{% endmacro %} + +{% macro support_code() %} +{% for line in support_code_lines %} +{{line}} +{% endfor %} +{% endmacro %} diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index 458ff1bc1..1e0814563 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -121,7 +121,7 @@ def __init__(self, source, variables, record=None, when=None, self.specifiers = {} for idx, variable in enumerate(variables): spec = source.specifiers[variable] - self.specifiers['_source_' + variable] = ArrayVariable('_source_'+variable, + self.specifiers['_source_' + variable] = ArrayVariable(variable, spec.unit, spec.dtype, spec.array, From 3626cedc7c60825ff5503d999c8b1d15d339be57 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 15:20:12 +0200 Subject: [PATCH 05/15] Fix recording several state variables with C++ code --- brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp b/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp index 380edbf7b..94e0d27db 100644 --- a/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp +++ b/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp @@ -39,6 +39,7 @@ _t_data[_new_len - 1] = _clock_t; {% for _varname in _variable_names %} + { PyArrayObject *_record_data = (((PyArrayObject*)(PyObject*)_recorded_{{_varname}}.attr("data"))); const npy_intp* _record_strides = _record_data->strides; for (int _idx=0; _idx < _num_indices; _idx++) @@ -51,6 +52,7 @@ double *recorded_entry = ((double*)(_record_data->data + (_new_len - 1)*_record_strides[0] + _idx*_record_strides[1])); *recorded_entry = _to_record_{{_varname}}; } + } {% endfor %} {% endmacro %} From baacef4711a4272a08f01b64e73bf07a4d693957 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 16:03:26 +0200 Subject: [PATCH 06/15] Use numpy_resize for DynamicArrays --- brian2/monitors/statemonitor.py | 10 +++++++--- brian2/synapses/synapses.py | 18 ++++++++++++------ 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index f9c82f766..c68f9350c 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -8,6 +8,7 @@ AttributeValue, Index) from brian2.core.base import BrianObject from brian2.core.scheduler import Scheduler +from brian2.core.preferences import brian_prefs from brian2.groups.group import Group from brian2.units.fundamentalunits import Unit from brian2.units.allunits import second @@ -31,7 +32,7 @@ def __init__(self, monitor): self.monitor = weakref.proxy(monitor) def get_value(self): - return self.monitor._t[:] + return self.monitor._t.data.copy() class StateMonitor(BrianObject, Group): ''' @@ -147,9 +148,12 @@ def __init__(self, source, variables, record=None, when=None, Group.__init__(self) def reinit(self): - self._values = dict((v, DynamicArray((0, len(self.variables)))) + self._values = dict((v, DynamicArray((0, len(self.variables)), + use_numpy_resize=True, + dtype=self.source.specifiers[v].dtype)) for v in self.variables) - self._t = DynamicArray1D(0) + self._t = DynamicArray1D(0, use_numpy_resize=True, + dtype=brian_prefs['core.default_scalar_dtype']) def pre_run(self, namespace): # Some dummy code so that code generation takes care of the indexing diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py index 3a773bc41..4c7c032a2 100644 --- a/brian2/synapses/synapses.py +++ b/brian2/synapses/synapses.py @@ -136,7 +136,9 @@ def __init__(self, synapses, code, prepost, objname=None): indices = {'_neuron_idx': Index('_neuron_idx', False), '_postsynaptic_idx': Index('_postsynaptic_idx', False), '_presynaptic_idx': Index('_presynaptic_idx', False)} - self._delays = DynamicArray1D(len(synapses.indices), dtype=np.float64) + self._delays = DynamicArray1D(len(synapses.indices), + use_numpy_resize=True, + dtype=np.float64) self.queue = SpikeQueue() self.spiking_synapses = [] self.specifiers = {'_spiking_synapses': AttributeValue('_spiking_synapses', @@ -311,11 +313,15 @@ def __init__(self, synapses): self.target_len = len(synapses.target) self.synapses = weakref.proxy(synapses) dtype = smallest_inttype(MAX_SYNAPSES) - self.synaptic_pre = DynamicArray1D(0, dtype=dtype) - self.synaptic_post = DynamicArray1D(0, dtype=dtype) - self.pre_synaptic = [DynamicArray1D(0, dtype=dtype) + self.synaptic_pre = DynamicArray1D(0, use_numpy_resize=True, + dtype=dtype) + self.synaptic_post = DynamicArray1D(0, use_numpy_resize=True, + dtype=dtype) + self.pre_synaptic = [DynamicArray1D(0, use_numpy_resize=True, + dtype=dtype) for _ in xrange(self.source_len)] - self.post_synaptic = [DynamicArray1D(0, dtype=dtype) + self.post_synaptic = [DynamicArray1D(0, use_numpy_resize=True, + dtype=dtype) for _ in xrange(self.target_len)] self.i = IndexView(self, self.synaptic_pre) self.j = IndexView(self, self.synaptic_post) @@ -920,7 +926,7 @@ def _allocate_memory(self, dtype=None): curdtype = dtype if curdtype is None: curdtype = brian_prefs['core.default_scalar_dtype'] - arrays[name] = DynamicArray1D(0) + arrays[name] = DynamicArray1D(0, use_numpy_resize=True) logger.debug("NeuronGroup memory allocated successfully.") return arrays From 48cc7754788d96340dbb913e6f9dd121dcf13710 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 17:13:15 +0200 Subject: [PATCH 07/15] Make StateMonitor work again for subexpressions --- brian2/monitors/statemonitor.py | 38 ++++++++++++++++++++------------- brian2/tests/test_monitor.py | 14 +++++++----- 2 files changed, 32 insertions(+), 20 deletions(-) diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index c68f9350c..78efbd76b 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -2,8 +2,7 @@ import numpy as np -from brian2.codegen.codeobject import create_codeobject -from brian2.core.specifiers import (Value, +from brian2.core.specifiers import (Value, Subexpression, ReadOnlyValue, ArrayVariable, AttributeValue, Index) from brian2.core.base import BrianObject @@ -13,6 +12,7 @@ from brian2.units.fundamentalunits import Unit from brian2.units.allunits import second from brian2.memory.dynamicarray import DynamicArray, DynamicArray1D +from brian2.groups.group import create_runner_codeobj __all__ = ['StateMonitor'] @@ -119,17 +119,22 @@ def __init__(self, source, variables, record=None, when=None, # initialise Group access self.specifiers = {} - for idx, variable in enumerate(variables): + for variable in variables: spec = source.specifiers[variable] - self.specifiers['_source_' + variable] = ArrayVariable(variable, - spec.unit, - spec.dtype, - spec.array, - '_record_idx', - group=spec.group, - constant=spec.constant, - scalar=spec.scalar, - is_bool=spec.is_bool) + if isinstance(spec, ArrayVariable): + self.specifiers['_source_' + variable] = ArrayVariable(variable, + spec.unit, + spec.dtype, + spec.array, + '_record_idx', + group=spec.group, + constant=spec.constant, + scalar=spec.scalar, + is_bool=spec.is_bool) + elif isinstance(spec, Subexpression): + self.specifiers['_source_' + variable] = weakref.proxy(spec) + else: + raise TypeError('Variable %s cannot be recorded.' % variable) self.specifiers[variable] = MonitorVariable(variable, spec.unit, spec.dtype, @@ -160,11 +165,14 @@ def pre_run(self, namespace): # and subexpressions code = ['_to_record_%s = _source_%s' % (v, v) for v in self.variables] + code += ['_recorded_%s = _recorded_%s' % (v, v) + for v in self.variables] code = '\n'.join(code) - self.codeobj = create_codeobject(self.name, + self.codeobj = create_runner_codeobj(self.source, code, - {}, # no namespace - self.specifiers, + name=self.name, + additional_specifiers=self.specifiers, + additional_namespace=namespace, template_name='statemonitor', indices={'_record_idx': Index('_record_idx', self.record_all)}, template_kwds={'_variable_names': self.variables}, diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py index 14feab65d..a32c923ec 100644 --- a/brian2/tests/test_monitor.py +++ b/brian2/tests/test_monitor.py @@ -6,24 +6,25 @@ # We can only test C++ if weave is availabe try: import scipy.weave - languages = [PythonLanguage(), CPPLanguage()] + languages = ['numpy', 'weave'] except ImportError: # Can't test C++ - languages = [PythonLanguage()] + languages = ['numpy'] def test_spike_monitor(): + language_before = brian_prefs.codegen.target for language in languages: + brian_prefs.codegen.target = language defaultclock.t = 0*second G = NeuronGroup(2, '''dv/dt = rate : 1 - rate: Hz''', threshold='v>1', reset='v=0', - language=language) + rate: Hz''', threshold='v>1', reset='v=0') # We don't use 100 and 1000Hz, because then the membrane potential would # be exactly at 1 after 10 resp. 100 timesteps. Due to floating point # issues this will not be exact, G.rate = [101, 1001] * Hz - mon = SpikeMonitor(G, language=language) + mon = SpikeMonitor(G) net = Network(G, mon) net.run(10*ms) @@ -31,6 +32,9 @@ def test_spike_monitor(): assert_allclose(mon.t[mon.i == 1], np.arange(10)*ms + 0.9*ms) assert_equal(mon.count, np.array([1, 10])) + brian_prefs.codegen.target = language_before + + if __name__ == '__main__': test_spike_monitor() From 1f2642caa42fb0a00cf13c6d8e7b0d01ac1c9ce7 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 17:36:42 +0200 Subject: [PATCH 08/15] Switch off refchecks for synaptic indices -- this way we can use the faster resize. No code should ever store a reference to the underlying data attribute. --- brian2/synapses/synapses.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py index 4c7c032a2..08a716433 100644 --- a/brian2/synapses/synapses.py +++ b/brian2/synapses/synapses.py @@ -314,14 +314,14 @@ def __init__(self, synapses): self.synapses = weakref.proxy(synapses) dtype = smallest_inttype(MAX_SYNAPSES) self.synaptic_pre = DynamicArray1D(0, use_numpy_resize=True, - dtype=dtype) + dtype=dtype, refcheck=False) self.synaptic_post = DynamicArray1D(0, use_numpy_resize=True, - dtype=dtype) + dtype=dtype, refcheck=False) self.pre_synaptic = [DynamicArray1D(0, use_numpy_resize=True, - dtype=dtype) + dtype=dtype, refcheck=False) for _ in xrange(self.source_len)] self.post_synaptic = [DynamicArray1D(0, use_numpy_resize=True, - dtype=dtype) + dtype=dtype, refcheck=False) for _ in xrange(self.target_len)] self.i = IndexView(self, self.synaptic_pre) self.j = IndexView(self, self.synaptic_post) @@ -926,7 +926,9 @@ def _allocate_memory(self, dtype=None): curdtype = dtype if curdtype is None: curdtype = brian_prefs['core.default_scalar_dtype'] - arrays[name] = DynamicArray1D(0, use_numpy_resize=True) + arrays[name] = DynamicArray1D(0, use_numpy_resize=True, + dtype=curdtype, + refcheck=False) logger.debug("NeuronGroup memory allocated successfully.") return arrays From ac84555a9fff91daf43bdcf12c02c46a003114d4 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 18:12:05 +0200 Subject: [PATCH 09/15] Use numpy's concatenate instead of hstack, it's apparently faster --- brian2/synapses/synapses.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py index 08a716433..b76b8ea84 100644 --- a/brian2/synapses/synapses.py +++ b/brian2/synapses/synapses.py @@ -194,8 +194,8 @@ def pre_update(self): # Push new spikes into the queue spikes = self.source.spikes if len(spikes): - indices = np.hstack((self.synapse_indices[spike] - for spike in spikes)).astype(np.int32) + indices = np.concatenate([self.synapse_indices[spike] + for spike in spikes]).astype(np.int32) if len(indices): if len(self._delays) > 1: delays = np.round(self._delays[indices] / self.dt).astype(int) From 0af38b0d789ae655054bdf9e8e46f1777e781e81 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Tue, 23 Jul 2013 19:40:15 +0200 Subject: [PATCH 10/15] Fix bugs and add tests for StateMonitor, no longer derive from Groups (this doesn't play nicely with subexpressions) --- .../numpy_rt/templates/statemonitor.py_ | 2 +- .../weave_rt/templates/statemonitor.cpp | 8 +- brian2/monitors/statemonitor.py | 90 ++++++++----------- brian2/tests/test_monitor.py | 45 ++++++++++ 4 files changed, 90 insertions(+), 55 deletions(-) diff --git a/brian2/codegen/runtime/numpy_rt/templates/statemonitor.py_ b/brian2/codegen/runtime/numpy_rt/templates/statemonitor.py_ index ab287297a..640f432f1 100644 --- a/brian2/codegen/runtime/numpy_rt/templates/statemonitor.py_ +++ b/brian2/codegen/runtime/numpy_rt/templates/statemonitor.py_ @@ -12,7 +12,7 @@ _recorded_{{_varname}}.resize((_new_len, _num_indices)) _t[-1] = _clock_t _vectorisation_idx = _indices -_record_idx = _indices +_neuron_idx = _indices[:] {% for line in code_lines %} {{line}} {% endfor %} diff --git a/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp b/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp index 94e0d27db..0dbeb8d04 100644 --- a/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp +++ b/brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp @@ -44,12 +44,14 @@ const npy_intp* _record_strides = _record_data->strides; for (int _idx=0; _idx < _num_indices; _idx++) { - const int _record_idx = _indices[_idx]; - const int _vectorisation_idx = _record_idx; + const int _neuron_idx = _indices[_idx]; + const int _vectorisation_idx = _neuron_idx; {% for line in code_lines %} {{line}} {% endfor %} - double *recorded_entry = ((double*)(_record_data->data + (_new_len - 1)*_record_strides[0] + _idx*_record_strides[1])); + + // FIXME: This will not work for variables with other data types + double *recorded_entry = (double*)(_record_data->data + (_new_len - 1)*_record_strides[0] + _idx*_record_strides[1]); *recorded_entry = _to_record_{{_varname}}; } } diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index 78efbd76b..be8da2d62 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -8,33 +8,16 @@ from brian2.core.base import BrianObject from brian2.core.scheduler import Scheduler from brian2.core.preferences import brian_prefs -from brian2.groups.group import Group -from brian2.units.fundamentalunits import Unit +from brian2.units.fundamentalunits import Unit, Quantity from brian2.units.allunits import second from brian2.memory.dynamicarray import DynamicArray, DynamicArray1D from brian2.groups.group import create_runner_codeobj +from units.fundamentalunits import have_same_dimensions __all__ = ['StateMonitor'] -class MonitorVariable(Value): - def __init__(self, name, unit, dtype, monitor): - Value.__init__(self, name, unit, dtype) - self.monitor = weakref.proxy(monitor) - - def get_value(self): - return self.monitor._values[self.name] - - -class MonitorTime(Value): - def __init__(self, monitor): - Value.__init__(self, 't', second, np.float64) - self.monitor = weakref.proxy(monitor) - - def get_value(self): - return self.monitor._t.data.copy() - -class StateMonitor(BrianObject, Group): +class StateMonitor(BrianObject): ''' Record values of state variables during a run @@ -100,16 +83,19 @@ def __init__(self, source, variables, record=None, when=None, elif isinstance(variables, str): variables = [variables] self.variables = variables - + + if len(self.variables) == 0: + raise ValueError('No variables to record') + # record should always be an array of ints self.record_all = False if record is None or record is False: - record = np.array([], dtype=int) + record = np.array([], dtype=np.int32) elif record is True: self.record_all = True - record = np.arange(len(source)) + record = np.arange(len(source), dtype=np.int32) else: - record = np.array(record, dtype=int) + record = np.array(record, dtype=np.int32) #: The array of recorded indices self.indices = record @@ -117,43 +103,27 @@ def __init__(self, source, variables, record=None, when=None, # create data structures self.reinit() - # initialise Group access + # Setup specifiers self.specifiers = {} for variable in variables: spec = source.specifiers[variable] - if isinstance(spec, ArrayVariable): - self.specifiers['_source_' + variable] = ArrayVariable(variable, - spec.unit, - spec.dtype, - spec.array, - '_record_idx', - group=spec.group, - constant=spec.constant, - scalar=spec.scalar, - is_bool=spec.is_bool) - elif isinstance(spec, Subexpression): - self.specifiers['_source_' + variable] = weakref.proxy(spec) - else: - raise TypeError('Variable %s cannot be recorded.' % variable) - self.specifiers[variable] = MonitorVariable(variable, - spec.unit, - spec.dtype, - self) + self.specifiers[variable] = weakref.proxy(spec) + self.specifiers['_recorded_'+variable] = ReadOnlyValue('_recorded_'+variable, Unit(1), self._values[variable].dtype, self._values[variable]) + self.specifiers['_t'] = ReadOnlyValue('_t', Unit(1), self._t.dtype, self._t) - self.specifiers['_clock_t'] = AttributeValue('t', second, np.float64, self.clock, 't_') - - self.specifiers['t'] = MonitorTime(self) - - Group.__init__(self) + self.specifiers['_indices'] = ArrayVariable('_indices', Unit(1), + np.int32, self.indices, + index='', group=None, + constant=True) def reinit(self): - self._values = dict((v, DynamicArray((0, len(self.variables)), + self._values = dict((v, DynamicArray((0, len(self.indices)), use_numpy_resize=True, dtype=self.source.specifiers[v].dtype)) for v in self.variables) @@ -163,7 +133,7 @@ def reinit(self): def pre_run(self, namespace): # Some dummy code so that code generation takes care of the indexing # and subexpressions - code = ['_to_record_%s = _source_%s' % (v, v) + code = ['_to_record_%s = %s' % (v, v) for v in self.variables] code += ['_recorded_%s = _recorded_%s' % (v, v) for v in self.variables] @@ -174,13 +144,31 @@ def pre_run(self, namespace): additional_specifiers=self.specifiers, additional_namespace=namespace, template_name='statemonitor', - indices={'_record_idx': Index('_record_idx', self.record_all)}, + indices={'_neuron_idx': Index('_neuron_idx', self.record_all)}, template_kwds={'_variable_names': self.variables}, codeobj_class=self.codeobj_class) def update(self): self.codeobj() + def __getattr__(self, item): + # TODO: Decide about the interface + if item == 't': + return Quantity(self._t.data.copy(), dim=second.dim) + elif item == 't_': + return self._t.data.copy() + elif item in self.variables: + unit = self.specifiers[item].unit + if have_same_dimensions(unit, 1): + return self._values[item].data.copy() + else: + return Quantity(self._values[item].data.copy(), + dim=unit.dim) + elif item.endswith('_') and item[:-1] in self.variables: + return self._values[item[:-1]].data.copy() + else: + getattr(super(StateMonitor, self), item) + def __repr__(self): description = '<{classname}, recording {variables} from {source}>' return description.format(classname=self.__class__.__name__, diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py index a32c923ec..073c0a4af 100644 --- a/brian2/tests/test_monitor.py +++ b/brian2/tests/test_monitor.py @@ -35,6 +35,51 @@ def test_spike_monitor(): brian_prefs.codegen.target = language_before +def test_state_monitor(): + # Record all neurons + language_before = brian_prefs.codegen.target + for language in languages: + brian_prefs.codegen.target = language + defaultclock.t = 0*second + # Check that all kinds of variables can be recorded + G = NeuronGroup(2, '''dv/dt = -v / (10*ms) : 1 + f = clip(v, 0.1, 0.9) : 1 + rate: Hz''', threshold='v>1', reset='v=0') + G.rate = [100, 1000] * Hz + G.v = 1 + + # Use a single StateMonitor + v_mon = StateMonitor(G, 'v', record=True) + v_mon1 = StateMonitor(G, 'v', record=[1]) + + # Use a StateMonitor for specified variables + multi_mon = StateMonitor(G, ['v', 'f', 'rate'], record=True) + multi_mon1 = StateMonitor(G, ['v', 'f', 'rate'], record=[1]) + + net = Network(G, v_mon, v_mon1, + multi_mon, multi_mon1) + net.run(10*ms) + + # Check v recording + assert_allclose(v_mon.v, + np.exp(np.tile(-v_mon.t - defaultclock.dt, (2, 1)).T / (10*ms))) + assert_allclose(v_mon.v_, + np.exp(np.tile(-v_mon.t_ - defaultclock.dt_, (2, 1)).T / float(10*ms))) + assert_equal(v_mon.v, multi_mon.v) + assert_equal(v_mon.v_, multi_mon.v_) + assert_equal(v_mon.v[:, 1:2], v_mon1.v) + assert_equal(multi_mon.v[:, 1:2], multi_mon1.v) + + # Other variables + assert_equal(multi_mon.rate_, np.tile(np.atleast_2d(G.rate_), + (multi_mon.rate.shape[0], 1))) + assert_equal(multi_mon.rate[:, 1:2], multi_mon1.rate) + assert_allclose(np.clip(multi_mon.v, 0.1, 0.9), multi_mon.f) + assert_allclose(np.clip(multi_mon1.v, 0.1, 0.9), multi_mon1.f) + + brian_prefs.codegen.target = language_before + if __name__ == '__main__': test_spike_monitor() + test_state_monitor() From 2ad5792be02890fe6032cf3426b9f5a55312ceca Mon Sep 17 00:00:00 2001 From: mstimberg Date: Tue, 23 Jul 2013 21:30:51 +0200 Subject: [PATCH 11/15] Fix an incorrect import that sneaked in --- brian2/monitors/statemonitor.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index be8da2d62..4d1a4c469 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -2,17 +2,15 @@ import numpy as np -from brian2.core.specifiers import (Value, Subexpression, - ReadOnlyValue, ArrayVariable, +from brian2.core.specifiers import (ReadOnlyValue, ArrayVariable, AttributeValue, Index) from brian2.core.base import BrianObject from brian2.core.scheduler import Scheduler from brian2.core.preferences import brian_prefs -from brian2.units.fundamentalunits import Unit, Quantity +from brian2.units.fundamentalunits import Unit, Quantity, have_same_dimensions from brian2.units.allunits import second from brian2.memory.dynamicarray import DynamicArray, DynamicArray1D from brian2.groups.group import create_runner_codeobj -from units.fundamentalunits import have_same_dimensions __all__ = ['StateMonitor'] From b0f3183a04c1465ec9205df49f9e0c86651eea3f Mon Sep 17 00:00:00 2001 From: mstimberg Date: Tue, 23 Jul 2013 22:58:28 +0200 Subject: [PATCH 12/15] Make the __getattr__ mechanism in StateMonitor Python3-compatible (we should document/implement a standard mechanism for this, we are using it in several places) --- brian2/monitors/statemonitor.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index 4d1a4c469..e6cb2525b 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -73,8 +73,9 @@ def __init__(self, source, variables, record=None, when=None, scheduler.clock = source.clock if not scheduler.defined_when: scheduler.when = 'end' + BrianObject.__init__(self, when=scheduler, name=name) - + # variables should always be a list of strings if variables is True: variables = source.equations.names @@ -120,6 +121,8 @@ def __init__(self, source, variables, record=None, when=None, index='', group=None, constant=True) + self._group_attribute_access_active = True + def reinit(self): self._values = dict((v, DynamicArray((0, len(self.indices)), use_numpy_resize=True, @@ -150,6 +153,17 @@ def update(self): self.codeobj() def __getattr__(self, item): + # We do this because __setattr__ and __getattr__ are not active until + # _group_attribute_access_active attribute is set, and if it is set, + # then __getattr__ will not be called. Therefore, if getattr is called + # with this name, it is because it hasn't been set yet and so this + # method should raise an AttributeError to agree that it hasn't been + # called yet. + if item == '_group_attribute_access_active': + raise AttributeError + if not hasattr(self, '_group_attribute_access_active'): + raise AttributeError + # TODO: Decide about the interface if item == 't': return Quantity(self._t.data.copy(), dim=second.dim) @@ -165,7 +179,7 @@ def __getattr__(self, item): elif item.endswith('_') and item[:-1] in self.variables: return self._values[item[:-1]].data.copy() else: - getattr(super(StateMonitor, self), item) + raise AttributeError('Unknown attribute %s' % item) def __repr__(self): description = '<{classname}, recording {variables} from {source}>' From 6ca3cf092b3c89cc26ff5d1f36c9de398e120ac6 Mon Sep 17 00:00:00 2001 From: mstimberg Date: Tue, 23 Jul 2013 23:08:48 +0200 Subject: [PATCH 13/15] More tests for the monitors --- brian2/monitors/statemonitor.py | 3 --- brian2/tests/test_monitor.py | 23 ++++++++++++++++++++++- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/brian2/monitors/statemonitor.py b/brian2/monitors/statemonitor.py index e6cb2525b..c9b876fa5 100644 --- a/brian2/monitors/statemonitor.py +++ b/brian2/monitors/statemonitor.py @@ -83,9 +83,6 @@ def __init__(self, source, variables, record=None, when=None, variables = [variables] self.variables = variables - if len(self.variables) == 0: - raise ValueError('No variables to record') - # record should always be an array of ints self.record_all = False if record is None or record is False: diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py index 073c0a4af..0cf5effe6 100644 --- a/brian2/tests/test_monitor.py +++ b/brian2/tests/test_monitor.py @@ -30,8 +30,17 @@ def test_spike_monitor(): assert_allclose(mon.t[mon.i == 0], [9.9]*ms) assert_allclose(mon.t[mon.i == 1], np.arange(10)*ms + 0.9*ms) + assert_allclose(mon.t_[mon.i == 0], np.array([9.9*float(ms)])) + assert_allclose(mon.t_[mon.i == 1], (np.arange(10) + 0.9)*float(ms)) assert_equal(mon.count, np.array([1, 10])) + i, t = mon.it + i_, t_ = mon.it_ + assert_equal(i, mon.i) + assert_equal(i, i_) + assert_equal(t, mon.t) + assert_equal(t_, mon.t_) + brian_prefs.codegen.target = language_before @@ -48,6 +57,10 @@ def test_state_monitor(): G.rate = [100, 1000] * Hz G.v = 1 + # A bit peculiar, but in principle one should be allowed to record + # nothing except for the time + nothing_mon = StateMonitor(G, [], record=True) + # Use a single StateMonitor v_mon = StateMonitor(G, 'v', record=True) v_mon1 = StateMonitor(G, 'v', record=[1]) @@ -56,10 +69,18 @@ def test_state_monitor(): multi_mon = StateMonitor(G, ['v', 'f', 'rate'], record=True) multi_mon1 = StateMonitor(G, ['v', 'f', 'rate'], record=[1]) - net = Network(G, v_mon, v_mon1, + net = Network(G, nothing_mon, + v_mon, v_mon1, multi_mon, multi_mon1) net.run(10*ms) + # Check time recordings + assert_equal(nothing_mon.t, v_mon.t) + assert_equal(nothing_mon.t_, np.asarray(nothing_mon.t)) + assert_equal(nothing_mon.t_, v_mon.t_) + assert_allclose(nothing_mon.t, + np.arange(len(nothing_mon.t)) * defaultclock.dt) + # Check v recording assert_allclose(v_mon.v, np.exp(np.tile(-v_mon.t - defaultclock.dt, (2, 1)).T / (10*ms))) From b05aeb8c5e0d16c5208baefcca4d43d7c9512f4c Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Wed, 24 Jul 2013 11:58:12 +0200 Subject: [PATCH 14/15] Add a simple PopulationRateMonitor (no binning, no smoothing) --- .../numpy_rt/templates/ratemonitor.py_ | 7 + .../weave_rt/templates/ratemonitor.cpp | 22 ++++ brian2/monitors/__init__.py | 1 + brian2/monitors/ratemonitor.py | 123 ++++++++++++++++++ brian2/tests/test_monitor.py | 34 ++++- 5 files changed, 183 insertions(+), 4 deletions(-) create mode 100644 brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_ create mode 100644 brian2/codegen/runtime/weave_rt/templates/ratemonitor.cpp create mode 100644 brian2/monitors/ratemonitor.py diff --git a/brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_ b/brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_ new file mode 100644 index 000000000..cc966b95e --- /dev/null +++ b/brian2/codegen/runtime/numpy_rt/templates/ratemonitor.py_ @@ -0,0 +1,7 @@ +# { USE_SPECIFIERS _rate, _t, _spikes, _num_source_neurons, t, dt } + +_new_len = len(_t) + 1 +_t.resize(_new_len) +_rate.resize(_new_len) +_t[-1] = t +_rate[-1] = 1.0 * len(_spikes) / dt / _num_source_neurons diff --git a/brian2/codegen/runtime/weave_rt/templates/ratemonitor.cpp b/brian2/codegen/runtime/weave_rt/templates/ratemonitor.cpp new file mode 100644 index 000000000..f2dd28133 --- /dev/null +++ b/brian2/codegen/runtime/weave_rt/templates/ratemonitor.cpp @@ -0,0 +1,22 @@ +{% macro main() %} + + // USE_SPECIFIERS { _t, _rate, t, dt, _spikes } + + // Calculate the new length for the arrays + const npy_int _new_len = (npy_int)(_t.attr("shape")[0]) + 1; + + // Resize the arrays + PyObject_CallMethod(_t, "resize", "i", _new_len); + PyObject_CallMethod(_rate, "resize", "i", _new_len); + // Get the potentially newly created underlying data arrays + double *_t_data = (double*)(((PyArrayObject*)(PyObject*)_t.attr("data"))->data); + double *_rate_data = (double*)(((PyArrayObject*)(PyObject*)_rate.attr("data"))->data); + + //Set the new values + _t_data[_new_len - 1] = t; + _rate_data[_new_len - 1] = 1.0 * _num_spikes / (double)dt / _num_source_neurons; + +{% endmacro %} + +{% macro support_code() %} +{% endmacro %} diff --git a/brian2/monitors/__init__.py b/brian2/monitors/__init__.py index 9b2211944..bffb6a0d4 100644 --- a/brian2/monitors/__init__.py +++ b/brian2/monitors/__init__.py @@ -1,2 +1,3 @@ from spikemonitor import * from statemonitor import * +from ratemonitor import * \ No newline at end of file diff --git a/brian2/monitors/ratemonitor.py b/brian2/monitors/ratemonitor.py new file mode 100644 index 000000000..57e11ac36 --- /dev/null +++ b/brian2/monitors/ratemonitor.py @@ -0,0 +1,123 @@ +import weakref + +import numpy as np + +from brian2.codegen.codeobject import create_codeobject +from brian2.core.base import BrianObject +from brian2.core.preferences import brian_prefs +from brian2.core.scheduler import Scheduler +from brian2.core.specifiers import ReadOnlyValue, AttributeValue, ArrayVariable +from brian2.memory.dynamicarray import DynamicArray1D +from brian2.units.allunits import second, hertz +from brian2.units.fundamentalunits import Unit, Quantity + +__all__ = ['PopulationRateMonitor'] + + +class PopulationRateMonitor(BrianObject): + ''' + Record instantaneous firing rates, averaged across neurons from a + `NeuronGroup` or other spike source. + + Parameters + ---------- + source : (`NeuronGroup`, `SpikeSource`) + The source of spikes to record. + when : `Scheduler`, optional + When to record the spikes, by default uses the clock of the source + and records spikes in the slot 'end'. + name : str, optional + A unique name for the object, otherwise will use + ``source.name+'_ratemonitor_0'``, etc. + codeobj_class : class, optional + The `CodeObject` class to run code with. + ''' + def __init__(self, source, when=None, name='ratemonitor*', + codeobj_class=None): + self.source = weakref.proxy(source) + + # run by default on source clock at the end + scheduler = Scheduler(when) + if not scheduler.defined_clock: + scheduler.clock = source.clock + if not scheduler.defined_when: + scheduler.when = 'end' + + self.codeobj_class = codeobj_class + BrianObject.__init__(self, when=scheduler, name=name) + + # create data structures + self.reinit() + + self.specifiers = {'t': AttributeValue('t', second, np.float64, + self.clock, 't'), + 'dt': AttributeValue('dt', second, np.float64, + self.clock, 'dt'), + '_spikes': AttributeValue('_spikes', Unit(1), + self.source.spikes.dtype, + self.source, 'spikes'), + # The template needs to have access to the + # DynamicArray here, having access to the underlying + # array is not enough since we want to do the resize + # in the template + '_rate': ReadOnlyValue('_rates', Unit(1), self._rate.dtype, + self._rate), + '_t': ReadOnlyValue('_t', Unit(1), self._t.dtype, + self._t), + '_num_source_neurons': ReadOnlyValue('_num_source_neurons', + Unit(1), np.int, + len(self.source))} + + def reinit(self): + ''' + Clears all recorded rates + ''' + self._rate = DynamicArray1D(0, use_numpy_resize=True, + dtype=brian_prefs['core.default_scalar_dtype']) + self._t = DynamicArray1D(0, use_numpy_resize=True, + dtype=getattr(self.clock.t, 'dtype', + np.dtype(type(self.clock.t)))) + + def pre_run(self, namespace): + self.codeobj = create_codeobject(self.name, + '', # No model-specific code + {}, # no namespace + self.specifiers, + template_name='ratemonitor', + indices={}) + + def update(self): + self.codeobj() + + @property + def rate(self): + ''' + Array of recorded rates (in units of Hz). + ''' + return Quantity(self._rate.data.copy(), dim=hertz.dim) + + @property + def rate_(self): + ''' + Array of recorded rates (unitless). + ''' + return self._rate.data.copy() + + @property + def t(self): + ''' + Array of recorded time points (in units of second). + ''' + return Quantity(self._t.data.copy(), dim=second.dim) + + @property + def t_(self): + ''' + Array of recorded time points (unitless). + ''' + return self._t.data.copy() + + def __repr__(self): + description = '<{classname}, recording {source}>' + return description.format(classname=self.__class__.__name__, + source=self.source.name) diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py index 0cf5effe6..e33df4b5e 100644 --- a/brian2/tests/test_monitor.py +++ b/brian2/tests/test_monitor.py @@ -45,7 +45,6 @@ def test_spike_monitor(): def test_state_monitor(): - # Record all neurons language_before = brian_prefs.codegen.target for language in languages: brian_prefs.codegen.target = language @@ -100,7 +99,34 @@ def test_state_monitor(): brian_prefs.codegen.target = language_before -if __name__ == '__main__': - test_spike_monitor() - test_state_monitor() +def test_rate_monitor(): + language_before = brian_prefs.codegen.target + for language in languages: + brian_prefs.codegen.target = language + G = NeuronGroup(5, 'v : 1', threshold='v>1') # no reset + G.v = 1.1 # All neurons spike every time step + rate_mon = PopulationRateMonitor(G) + net = Network(G, rate_mon) + net.run(10*defaultclock.dt) + + assert_allclose(rate_mon.t, np.arange(10) * defaultclock.dt) + assert_allclose(rate_mon.t_, np.arange(10) * float(defaultclock.dt)) + assert_allclose(rate_mon.rate, np.ones(10) / defaultclock.dt) + assert_allclose(rate_mon.rate_, np.asarray(np.ones(10) / defaultclock.dt)) + + defaultclock.t = 0*ms + G = NeuronGroup(10, 'v : 1', threshold='v>1') # no reset + G.v[:5] = 1.1 # Half of the neurons fire every time step + rate_mon = PopulationRateMonitor(G) + net = Network(G, rate_mon) + net.run(10*defaultclock.dt) + assert_allclose(rate_mon.rate, 0.5 * np.ones(10) / defaultclock.dt) + assert_allclose(rate_mon.rate_, 0.5 *np.asarray(np.ones(10) / defaultclock.dt)) + + brian_prefs.codegen.target = language_before + +if __name__ == '__main__': + #est_spike_monitor() + #test_state_monitor() + test_rate_monitor() From 633306d110087307740bcfe81ae8486cfc613cc2 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Wed, 24 Jul 2013 12:01:12 +0200 Subject: [PATCH 15/15] Don't comment out tests in the main block of test_monitor --- brian2/tests/test_monitor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/brian2/tests/test_monitor.py b/brian2/tests/test_monitor.py index e33df4b5e..9f07d3253 100644 --- a/brian2/tests/test_monitor.py +++ b/brian2/tests/test_monitor.py @@ -127,6 +127,6 @@ def test_rate_monitor(): brian_prefs.codegen.target = language_before if __name__ == '__main__': - #est_spike_monitor() - #test_state_monitor() + test_spike_monitor() + test_state_monitor() test_rate_monitor()