Merge 633306d into 1a93613

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

brian2/codegen/runtime/numpy_rt/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);

brian2/codegen/runtime/numpy_rt/templates/statemonitor.py_

@@ -0,0 +1,23 @@
+# USE_SPECIFIERS { _t, _clock_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[-1] = _clock_t
+
+_vectorisation_idx = _indices
+_neuron_idx = _indices[:]
+{% for line in code_lines %}
+{{line}}
+{% endfor %}
+
+{% for _varname in _variable_names %}
+_recorded_{{_varname}}[-1, :] = _to_record_{{_varname}}
+{% endfor %}
+

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 %}

brian2/codegen/runtime/weave_rt/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 %}

brian2/codegen/runtime/weave_rt/templates/statemonitor.cpp

@@ -0,0 +1,66 @@
+{% 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 _neuron_idx = _indices[_idx];
+            const int _vectorisation_idx = _neuron_idx;
+            {% for line in code_lines %}
+            {{line}}
+            {% endfor %}
+
+            // 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}};
+        }
+    }
+    {% endfor %}
+
+{% endmacro %}
+
+{% macro support_code() %}
+{% for line in support_code_lines %}
+{{line}}
+{% endfor %}
+{% endmacro %}

brian2/monitors/__init__.py

@@ -1,2 +1,3 @@
 from spikemonitor import *
 from statemonitor import *
+from ratemonitor import *

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)