Merge pull request #13 from brian-team/minor_fixes

Minor fixes

brian2modelfitting/fitter.py

@@ -1,4 +1,6 @@
 import abc
+
+from brian2.units.fundamentalunits import DIMENSIONLESS, get_dimensions
 from numpy import ones, array, arange, concatenate, mean, nanmin, reshape
 from brian2 import (NeuronGroup,  defaultclock, get_device, Network,
                     StateMonitor, SpikeMonitor, ms, device, second,
@@ -120,7 +122,7 @@ def __init__(self, dt, model, input, output, input_var, output_var,
         if isinstance(model, str):
             model = Equations(model)
         if input_var not in model.identifiers:
-            raise Exception("%s is not an identifier in the model" % input_var)
+            raise NameError("%s is not an identifier in the model" % input_var)
 
         defaultclock.dt = dt
         self.dt = dt
@@ -144,9 +146,17 @@ def __init__(self, dt, model, input, output, input_var, output_var,
         self.output_var = output_var
         self.model = model
 
+        input_dim = get_dimensions(input)
+        input_dim = '1' if input_dim is DIMENSIONLESS else repr(input_dim)
+        input_eqs = "{} = input_var(t, i % n_traces) : {}".format(input_var,
+                                                                  input_dim)
+        self.model += input_eqs
+
+        input_traces = TimedArray(input.transpose(), dt=dt)
+        self.input_traces = input_traces
+
         # initialization of attributes used later
         self.best_params = None
-        self.input_traces = None
         self.network = None
         self.optimizer = None
         self.metric = None
@@ -427,21 +437,15 @@ def __init__(self, model, input_var, input, output_var, output, dt,
                          param_init)
 
         if output_var not in self.model.names:
-            raise ValueError("%s is not a model variable" % output_var)
+            raise NameError("%s is not a model variable" % output_var)
         if output.shape != input.shape:
-            raise Exception("Input and output must have the same size")
+            raise ValueError("Input and output must have the same size")
 
-        # Replace input variable by TimedArray
         output_traces = TimedArray(output.transpose(), dt=dt)
-        input_traces = TimedArray(input.transpose(), dt=dt)
-        self.model = self.model + Equations(input_var + '= input_var(t, i % n_traces) :\
-                                  ' + "% s" % repr(input.dim))
-
-        self.input_traces = input_traces
 
         # Setup NeuronGroup
         namespace = get_local_namespace(level=level+1)
-        namespace['input_var'] = input_traces
+        namespace['input_var'] = self.input_traces
         namespace['output_var'] = output_traces
         namespace['n_traces'] = self.n_traces
         self.neurons = self.setup_neuron_group(self.n_neurons, namespace)
@@ -493,16 +497,9 @@ def __init__(self, model, input, output, dt, reset, threshold,
                          n_samples, threshold, reset, refractory, method,
                          param_init)
 
-        # Replace input variable by TimedArray
-        input_traces = TimedArray(input.transpose(), dt=dt)
-        self.model = self.model + Equations(input_var + '= input_var(t, i % n_traces) :\
-                                   ' + "% s" % repr(input.dim))
-
-        self.input_traces = input_traces
-
         # Setup NeuronGroup
         namespace = get_local_namespace(level=level+1)
-        namespace['input_var'] = input_traces
+        namespace['input_var'] = self.input_traces
         namespace['n_traces'] = self.n_traces
         self.neurons = self.setup_neuron_group(self.n_neurons, namespace)
 
@@ -555,22 +552,21 @@ def __init__(self, model, input_var, input, output_var, output, dt,
                          param_init)
 
         if output_var not in self.model.names:
-            raise ValueError("%s is not a model variable" % output_var)
+            raise NameError("%s is not a model variable" % output_var)
         if output.shape != input.shape:
-            raise Exception("Input and output must have the same size")
+            raise ValueError("Input and output must have the same size")
 
         # Replace input variable by TimedArray
         output_traces = TimedArray(output.transpose(), dt=dt)
-        input_traces = TimedArray(input.transpose(), dt=dt)
-        self.model = self.model + Equations(input_var + '= input_var(t, i % n_traces) :\
-                                  ' + "% s" % repr(input.dim))
-        self.model = self.model + Equations('total_error : %s' % repr(output.dim**2))
-
-        self.input_traces = input_traces
+        output_dim = get_dimensions(output)
+        squared_output_dim = ('1' if output_dim is DIMENSIONLESS
+                              else repr(output_dim**2))
+        error_eqs = Equations('total_error : {}'.format(squared_output_dim))
+        self.model = self.model + error_eqs
 
         # Setup NeuronGroup
         namespace = get_local_namespace(level=level+1)
-        namespace['input_var'] = input_traces
+        namespace['input_var'] = self.input_traces
         namespace['output_var'] = output_traces
         namespace['n_traces'] = self.n_traces
         self.neurons = self.setup_neuron_group(self.n_neurons, namespace)
@@ -595,7 +591,6 @@ def __init__(self, model, input_var, input, output_var, output, dt,
 
     def calc_errors(self, metric=None):
         """Calculates error in online fashion.To be used inside optim_iter."""
-        errors = self.simulator.network['neurons'].total_error/int((self.duration-self.t_start)/defaultclock.dt)
         errors = self.neurons.total_error/int((self.duration-self.t_start)/defaultclock.dt)
         errors = mean(errors.reshape((self.n_samples, self.n_traces)), axis=1)
         return array(errors)

brian2modelfitting/metric.py

@@ -17,8 +17,8 @@ def firing_rate(spikes):
     return (len(spikes) - 1) / (spikes[-1] - spikes[0])
 
 
-def get_gamma_factor(model, data, delta, time, dt):
-    """
+def get_gamma_factor(model, data, delta, time, dt, rate_correction=True):
+    r"""
     Calculate gamma factor between model and target spike trains,
     with precision delta.
 
@@ -32,11 +32,27 @@ def get_gamma_factor(model, data, delta, time, dt):
         time window
     dt: `~brian2.units.fundamentalunits.Quantity`
         time step
+    time: `~brian2.units.fundamentalunits.Quantity`
+        total time of the simulation
+    rate_correction: bool
+        Whether to include an error term that penalizes differences in firing
+        rate, following `Clopath et al., Neurocomputing (2007)
+        <https://doi.org/10.1016/j.neucom.2006.10.047>`_.
 
     Returns
     -------
-        float
-            The gamma factor
+    float
+        An error based on the Gamma factor. If ``rate_correction`` is used,
+        then the returned error is :math:`2\frac{\lvert r_\mathrm{data} - r_\mathrm{model}\rvert}{r_\mathrm{data}} - \Gamma`
+        (with :math:`r_\mathrm{data}` and :math:`r_\mathrm{model}` being the
+        firing rates in the data/model, and :math:`\Gamma` the coincidence
+        factor). Without ``rate_correction``, the error is
+        :math:`1 - \Gamma`. Note that the coincidence factor :math:`\Gamma`
+        has a maximum value of 1 (when the two spike trains are exactly
+        identical) and a value of 0 if there are only as many coincidences
+        as expected from two homogeneous Poisson processes of the same rate.
+        It can also take negative values if there are fewer coincidences
+        than expected by chance.
     """
     model = array(model)
     data = array(data)
@@ -58,7 +74,7 @@ def get_gamma_factor(model, data, delta, time, dt):
         matched_spikes = (diff <= delta_diff)
         coincidences = sum(matched_spikes)
     elif model_length == 0:
-        return 0
+        coincidences = 0
     else:
         indices = [amin(abs(model - data[i])) <= delta_diff for i in arange(data_length)]
         coincidences = sum(indices)
@@ -68,10 +84,12 @@ def get_gamma_factor(model, data, delta, time, dt):
     norm = .5*(1 - 2 * data_rate * delta)
     gamma = (coincidences - NCoincAvg)/(norm*(model_length + data_length))
 
-    corrected_gamma_factor = 2*abs((data_rate - model_rate)/data_rate - gamma)
-
-    return corrected_gamma_factor
+    if rate_correction:
+        rate_term = 2*abs((data_rate - model_rate)/data_rate)
+    else:
+        rate_term = 1
 
+    return rate_term - gamma
 
 def calc_eFEL(traces, inp_times, feat_list, dt):
     out_traces = []
@@ -431,14 +449,19 @@ class GammaFactor(SpikeMetric):
     Calculate gamma factors between goal and calculated spike trains, with
     precision delta.
 
-    Reference:
-    R. Jolivet et al., 'A benchmark test for a quantitative assessment of
-    simple neuron models',
-    Journal of Neuroscience Methods 169, no. 2 (2008): 417-424.
+    References:
+
+    * `R. Jolivet et al. “A Benchmark Test for a Quantitative Assessment of
+      Simple Neuron Models.” Journal of Neuroscience Methods, 169, no. 2 (2008):
+      417–24. <https://doi.org/10.1016/j.jneumeth.2007.11.006>`_
+    * `C. Clopath et al. “Predicting Neuronal Activity with Simple Models of the
+      Threshold Type: Adaptive Exponential Integrate-and-Fire Model with
+      Two Compartments.” Neurocomputing, 70, no. 10 (2007): 1668–73.
+      <https://doi.org/10.1016/j.neucom.2006.10.047>`_
     """
 
     @check_units(delta=second, time=second, t_start=0*second)
-    def __init__(self, delta, time, t_start=0*second):
+    def __init__(self, delta, time, t_start=0*second, rate_correction=True):
         """
         Initialize the metric with time window delta and time step dt output
 
@@ -448,18 +471,24 @@ def __init__(self, delta, time, t_start=0*second):
             time window
         time: `~brian2.units.fundamentalunits.Quantity`
             total length of experiment
+        rate_correciton: bool
+            Whether to include an error term that penalizes differences in firing
+            rate, following `Clopath et al., Neurocomputing (2007)
+            <https://doi.org/10.1016/j.neucom.2006.10.047>`_.
         """
         super(GammaFactor, self).__init__(t_start=t_start)
         self.delta = delta
         self.time = time
+        self.rate_correction = rate_correction
 
     def get_features(self, traces, output, dt):
         all_gf = []
         for one_sample in traces:
             gf_for_sample = []
             for model_response, target_response in zip(one_sample, output):
                 gf = get_gamma_factor(model_response, target_response,
-                                      self.delta, self.time, dt)
+                                      self.delta, self.time, dt,
+                                      rate_correction=self.rate_correction)
                 gf_for_sample.append(gf)
             all_gf.append(gf_for_sample)
         return array(all_gf)

brian2modelfitting/optimizer.py

@@ -21,8 +21,8 @@ def calc_bounds(parameter_names, **params):
         bounds for each parameter
     """
     for param in parameter_names:
-        if (param not in params):
-            raise Exception("Bounds must be set for parameter %s" % param)
+        if param not in params:
+            raise TypeError("Bounds must be set for parameter %s" % param)
 
     bounds = []
     for name in parameter_names:
@@ -131,9 +131,9 @@ def __init__(self,  method='DE', **kwds):
 
     def initialize(self, parameter_names, popsize, **params):
         for param in params.keys():
-            if (param not in parameter_names):
-                raise Exception("Parameter %s must be defined as a parameter "
-                                "in the model" % param)
+            if param not in parameter_names:
+                raise ValueError("Parameter %s must be defined as a parameter "
+                                 "in the model" % param)
 
         bounds = calc_bounds(parameter_names, **params)
 
@@ -200,9 +200,9 @@ def __init__(self, method='GP', **kwds):
 
     def initialize(self, parameter_names, popsize, **params):
         for param in params.keys():
-            if (param not in parameter_names):
-                raise Exception("Parameter %s must be defined as a parameter "
-                                "in the model" % param)
+            if param not in parameter_names:
+                raise ValueError("Parameter %s must be defined as a parameter "
+                                 "in the model" % param)
 
         bounds = calc_bounds(parameter_names, **params)
 

brian2modelfitting/tests/test_metric.py

@@ -18,16 +18,26 @@ def test_firing_rate():
 def test_get_gamma_factor():
     src = [7, 9, 11] * ms
     src2 = [1, 2, 3] * ms
-    trg = [0, 2, 4, 6, 8] * ms
+    trg = [0, 2, 4, 6, 8, 10] * ms
 
-    gf0 = get_gamma_factor(trg, trg, delta=12*ms, time=12*ms, dt=0.1*ms)
-    gf1 = get_gamma_factor(src2, trg, delta=12*ms, time=12*ms, dt=0.1*ms)
-    gf2 = get_gamma_factor(src, src2, delta=5*ms, time=5*ms, dt=0.1*ms)
+    gf0 = get_gamma_factor(trg, trg, delta=0.5*ms, time=12*ms, dt=0.1*ms)
+    gf1 = get_gamma_factor(src2, trg, delta=0.5*ms, time=12*ms, dt=0.1*ms)
+    gf2 = get_gamma_factor(src, src2, delta=0.5*ms, time=5*ms, dt=0.1*ms)
 
-    assert_equal(gf0, 2.0)
-    assert gf1 > 1.0
-    assert gf2 > 1.0
-    assert gf1 < gf2
+    assert_equal(gf0, -1)
+    assert gf1 > 0  # Since data rate = 2 * model rate
+    assert gf2 > -1
+
+    gf0 = get_gamma_factor(trg, trg, delta=0.5*ms, time=12*ms, dt=0.1*ms,
+                           rate_correction=False)
+    gf1 = get_gamma_factor(src2, trg, delta=0.5*ms, time=12*ms, dt=0.1*ms,
+                           rate_correction=False)
+    gf2 = get_gamma_factor(src, src2, delta=0.5*ms, time=5*ms, dt=0.1*ms,
+                           rate_correction=False)
+
+    assert_equal(gf0, 0)
+    assert gf1 > 0
+    assert gf2 > 0
 
 
 def test_init():
@@ -46,6 +56,7 @@ def test_calc_mse():
                  np.zeros(5))
     assert(np.all(mse.calc(inp, out, 0.1*ms) > 0))
 
+
 def test_calc_mse_t_start():
     mse = MSEMetric(t_start=1*ms)
     out = np.random.rand(2, 20)
@@ -59,20 +70,29 @@ def test_calc_mse_t_start():
     inp[:, :, 10:] = out[None, :, 10:]
     assert_equal(mse.calc(inp, out, 0.1*ms), np.zeros(5))
 
+
 def test_calc_gf():
     assert_raises(TypeError, GammaFactor)
-    assert_raises(DimensionMismatchError, GammaFactor, delta=10)
+    assert_raises(DimensionMismatchError, GammaFactor, delta=10*mV)
     assert_raises(DimensionMismatchError, GammaFactor, time=10)
 
-    inp_gf = np.round(np.sort(np.random.rand(5, 2, 5) * 10), 2)
-    out_gf = np.round(np.sort(np.random.rand(2, 5) * 10), 2)
+    model_spikes = [[np.array([1, 5, 8]), np.array([2, 3, 8, 9])],  # Correct rate
+                    [np.array([1, 5]), np.array([0, 2, 3, 8, 9])]]  # Wrong rate
+    data_spikes = [np.array([0, 5, 9]), np.array([1, 3, 5, 6])]
+
+    gf = GammaFactor(delta=0.5*ms, time=10*ms)
+    errors = gf.calc([data_spikes]*5, data_spikes, 0.1*ms)
+    assert_almost_equal(errors, np.ones(5)*-1)
+    errors = gf.calc(model_spikes, data_spikes, 0.1*ms)
+    assert errors[0] > -1  # correct rate
+    assert errors[1] > errors[0]
+
+    gf = GammaFactor(delta=0.5*ms, time=10*ms, rate_correction=False)
+    errors = gf.calc([data_spikes]*5, data_spikes, 0.1*ms)
+    assert_almost_equal(errors, np.zeros(5))
+    errors = gf.calc(model_spikes, data_spikes, 0.1*ms)
+    assert all(errors > 0)
 
-    gf = GammaFactor(delta=10*ms, time=10*ms)
-    errors = gf.calc(inp_gf, out_gf, 0.1*ms)
-    assert_equal(np.shape(errors), (5,))
-    assert(all(errors > 0))
-    errors = gf.calc([out_gf]*5, out_gf, 0.1*ms)
-    assert_almost_equal(errors, np.ones(5)*2)
 
 def test_get_features_mse():
     mse = MSEMetric()
@@ -100,17 +120,18 @@ def test_get_errors_mse():
 
 
 def test_get_features_gamma():
-    inp_gf = np.round(np.sort(np.random.rand(3, 2, 5) * 10), 2)
-    out_gf = np.round(np.sort(np.random.rand(2, 5) * 10), 2)
+    model_spikes = [[np.array([1, 5, 8]), np.array([2, 3, 8, 9])],  # Correct rate
+                    [np.array([1, 5]), np.array([0, 2, 3, 8, 9])]]  # Wrong rate
+    data_spikes = [np.array([0, 5, 9]), np.array([1, 3, 5, 6])]
 
-    gf = GammaFactor(delta=10*ms, time=10*ms)
-    features = gf.get_features(inp_gf, out_gf, 0.1*ms)
-    assert_equal(np.shape(features), (3, 2))
-    assert(np.all(np.array(features) > 0))
+    gf = GammaFactor(delta=0.5*ms, time=10*ms)
+    features = gf.get_features(model_spikes, data_spikes, 0.1*ms)
+    assert_equal(np.shape(features), (2, 2))
+    assert(np.all(np.array(features) > -1))
 
-    features = gf.get_features([out_gf]*3, out_gf, 0.1*ms)
+    features = gf.get_features([data_spikes]*3, data_spikes, 0.1*ms)
     assert_equal(np.shape(features), (3, 2))
-    assert_almost_equal(features, np.ones((3, 2))*2)
+    assert_almost_equal(features, np.ones((3, 2))*-1)
 
 
 def test_get_errors_gamma():