Fix calculation of stability matrix in Theory class, formatting in...

Fix calculation of stability matrix in Theory class, formatting in MultiAreaModel class, remove deprecated functions from analysis_helpers

multiarea_model/analysis_helpers.py

@@ -327,45 +327,6 @@ Analysis functions
 """
 
 
-def online_hist(fname, tmin, tmax, resolution=1.):
-    """
-    Compute spike histogram from gdf file by reading line after line.
-    Avoids reading in the entire file first, which makes it more suitable
-    for very large spike files.
-
-    Parameters
-    ----------
-    fname : string
-        Name of file to be read.
-    tmin : float
-        Minimal time for the calculation of the histogram.
-    tmax : float
-        Maximal time for the calculation of the histogram.
-    resolution : float, optional
-        Bin width of the histogram. Defaults to 1 ms.
-
-    Returns
-    -------
-    bins : numpy.ndarray
-        Array containing the left edges of the histogram bins
-    valyes : numpy.ndarray
-        Array containing the population rate value in each bin
-    """
-    gdf_file = open(fname, 'r')
-    bins = np.arange(tmin, tmax + resolution, resolution) + resolution / 2.
-    vals = np.zeros_like(bins)
-    current_bin_index = 0
-    for l in gdf_file:
-        data = l.split()
-        if np.logical_and(float(data[1]) > tmin, float(data[1]) < tmax):
-            while float(data[1]) >= bins[current_bin_index + 1]:
-                current_bin_index += 1
-            vals[current_bin_index] += 1
-        elif float(data[1]) > tmax:
-            break
-    return bins[:-1], vals[:-1] / (resolution / 1000.)
-
-
 def pop_rate(data_array, t_min, t_max, num_neur):
     """
     Calculates firing rate of a given array of spikes.

multiarea_model/multiarea_model.py

@@ -139,9 +139,9 @@ class MultiAreaModel:
                 raise NotImplementedError('Stabilization procedure has '
                                           'to be integrated.')
             elif isinstance(self.params['connection_params']['K_stable'], np.ndarray):
-                raise ValueError("Not supported. Please store the "
-                                 "matrix in a file and define the path to the file as "
-                                 "the parameter value.")
+                raise NotImplementedError("Not supported. Please store the "
+                                          "matrix in a file and define the path to the file as "
+                                          "the parameter value.")
             else:  # Assume that the parameter defines a filename containing the matrix
                 K_stable = np.load(self.params['connection_params']['K_stable'])
             ext = {area: {pop: ind[area][pop]['external'] for pop in

multiarea_model/theory.py

@@ -152,7 +152,7 @@ class Theory():
 
         # Loop over all iterations of different initial conditions
         for nsim in range(num_iter):
-            print("Iteration {}".format(nsim))
+            print('Iteration: {}'.format(nsim))
             initial_rates = next(gen)
             for ii in range(dim):
                 nest.SetStatus([neurons[ii]], {'rate': initial_rates[ii]})
@@ -344,7 +344,6 @@ class Theory():
             N_post[:, ii] = N
 
         # Connection probabilities between populations
-        C = 1. - (1.-1./(N_pre * N_post))**(K*N_post)
         mu, sigma = self.mu_sigma(rates)
 
         if np.any(vector_filter is not None):
@@ -370,17 +369,14 @@ class Theory():
         slope_matrix = np.zeros_like(J)
         slope_sigma_matrix = np.zeros_like(J)
         for ii in range(N.size):
-            slope_matrix[:, ii] = slope
-            slope_sigma_matrix[:, ii] = slope_sigma
-        V = C*(1-C)
-        G = (self.NP['tau_m'] * 1e-3)**2 * (slope_matrix*J +
-                                            slope_sigma_matrix*J**2)**2
-        G_N = N_pre * G
-        M = G_N * V
+            slope_mu_matrix[:, ii] = d_nu_d_mu
+            slope_sigma_matrix[:, ii] = d_nu_d_sigma
+        G = self.NP['tau_m'] * 1e-3 * (slope_mu_matrix * K * J +
+                                       slope_sigma_matrix * K * J**2)
         if full_output:
-            return M, slope, slope_sigma, M, C, V, G_N, J, N_pre
+            return G, d_nu_d_mu, d_nu_d_sigma
         else:
-            return M
+            return G
 
     def lambda_max(self, rates, matrix_filter=None,
                    vector_filter=None, full_output=False):