Remove integrate_siegert_python function introduced in 4a2c3dd6

multiarea_model/theory.py

@@ -64,7 +64,7 @@ class Theory():
     def __hash__(self):
         return hash(self.label)
 
-    def integrate_siegert_nest(self, full_output=True):
+    def integrate_siegert(self, full_output=True):
         """
         Integrate siegert formula to obtain stationary rates. See Eq. (3)
         and following in Schuecker, Schmidt et al. (2017).
@@ -187,99 +187,6 @@ class Theory():
         else:
             return self.network.structure_vec, rates
 
-    def integrate_siegert_python(self, full_output=True, parallel=True):
-        """
-        Integrate siegert formula to obtain stationary rates. See Eq. (3)
-        and following in Schuecker, Schmidt et al. (2017).
-
-        Use Runge-Kutta in Python.
-        """
-        dt = self.params['dt']
-        T = self.params['T']
-        K = copy(self.network.K_matrix)
-        tau = self.NP['tau_m'] * 1e-3
-        dim = np.shape(K)[0]
-
-        # Set initial rates of neurons:
-        # If defined as None, set them 0
-        if self.params['initial_rates'] is None:
-            num_iter = 1
-            gen = (np.zeros(dim) for ii in range(num_iter))
-        # iterate over different initial conditions drawn from a random distribution
-        elif self.params['initial_rates'] == 'random_uniform':
-                gen = self.initial_rates(self.params['initial_rates_iter'],
-                                         dim,
-                                         mode=self.params['initial_rates'],
-                                         rate_max=1000.)
-                num_iter = self.params['initial_rates_iter']
-        # initial rates are explicitly defined in self.params
-        elif isinstance(self.params['initial_rates'], np.ndarray):
-                num_iter = 1
-                gen = (self.params['initial_rates'] for ii in range(num_iter))
-
-        rates = []
-        # Loop over all iterations of different initial conditions
-        for nsim in range(num_iter):
-            print("Iteration {}".format(nsim))
-            initial_rates = next(gen)
-
-            self.t = np.arange(0, T, dt)
-            y = np.zeros((len(self.t), dim), dtype=float)
-            y[0, :] = initial_rates
-
-            # Integration loop
-            for i, tl in enumerate(self.t[:-1]):
-                print(tl)
-                delta_y, new_mu, new_sigma = self.Phi(
-                    y[i, :], parallel=parallel, return_mu_sigma=True)
-                k1 = delta_y
-                k2 = self.Phi(
-                    y[i, :] + dt * 1e-3 / tau / 2 * k1, parallel=parallel, return_mu_sigma=False)
-                k3 = self.Phi(
-                    y[i, :] + dt * 1e-3 / tau / 2 * k2, parallel=parallel, return_mu_sigma=False)
-                k4 = self.Phi(
-                    y[i, :] + dt * 1e-3 / tau * k3, parallel=parallel, return_mu_sigma=False)
-                y[i + 1, :] = y[i, :].copy() + dt * 1e-3 / 6. * (k1 + 2 * k2 + 2 * k3 + k4) / tau
-
-            if full_output:
-                rates.append(y.T)
-            else:
-                # Keep only initial and final rates
-                rates.append(y.T[:, [0, -1]])
-
-        if num_iter == 1:
-            return self.network.structure_vec, rates[0]
-        else:
-            return self.network.structure_vec, rates
-
-    def nu0_fb(self, arg):
-        mu = arg[0]
-        sigma = arg[1]
-        return nu0_fb(mu, sigma,
-                      self.NP['tau_m'] * 1e-3,
-                      self.NP['tau_syn'] * 1e-3,
-                      self.NP['t_ref'] * 1e-3,
-                      self.NP['theta'],
-                      self.NP['V_reset'])
-
-    def Phi(self, rates, parallel=False,
-            return_mu_sigma=False, return_leak=True):
-        mu, sigma = self.mu_sigma(rates, external=True)
-        if parallel:
-            pool = multiprocessing.Pool(processes=4)
-            new_rates = np.array(
-                pool.map(self.nu0_fb, zip(mu, sigma)))
-            pool.close()
-            pool.join()
-        else:
-            new_rates = [self.nu0_fb((m, s)) for m, s in zip(mu, sigma)]
-        if return_leak:
-            new_rates -= rates
-        if return_mu_sigma:
-            return (new_rates), mu, sigma
-        else:
-            return (new_rates)
-
     def replace_cc_input(self):
         """
         Helper function to replace cortico-cortical input by different variants.

tests/test_meanfield.py

@@ -10,4 +10,4 @@ def test_meanfield():
     network_params = {}
     theory_params = {}
     M0 = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
-    p, r0 = M0.theory.integrate_siegert_nest()
+    p, r0 = M0.theory.integrate_siegert()

tests/test_network_scaling.py

@@ -19,7 +19,7 @@ K0 = M0.K_matrix
 W0 = M0.W_matrix
 N0 = M0.N_vec
 syn0 = M0.syn_matrix
-p, r0 = M0.theory.integrate_siegert_nest()
+p, r0 = M0.theory.integrate_siegert()
 
 d = vector_to_dict(r0[:, -1],
                    M0.area_list,
@@ -39,7 +39,7 @@ K = M.K_matrix
 W = M.W_matrix
 N = M.N_vec
 syn = M.syn_matrix
-p, r = M.theory.integrate_siegert_nest()
+p, r = M.theory.integrate_siegert()
 
 assert(np.allclose(K, network_params['K_scaling'] * K0))
 assert(np.allclose(N, network_params['N_scaling'] * N0))
@@ -58,8 +58,24 @@ sigma0 = np.sqrt(1e-3 * tau_m * np.dot(M0.K_matrix * M0.J_matrix**2, r0_extend))
 sigma = np.sqrt(1e-3 * tau_m * np.dot(M.K_matrix * M.J_matrix**2, r0_extend))
 
 
-p, r0_py = M0.theory.integrate_siegert_python(parallel=False)
-p, r_py = M.theory.integrate_siegert_python(parallel=False)
+r0_alex_loaded = np.load('test_network_scaling_r0.npy')
+r_alex_loaded = np.load('test_network_scaling_r.npy')
+
+network_params = {}
+theory_params = {'initial_rates': r0_alex_loaded[:, -1],
+                 'T': 50.}
+M0_alex = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
+p, r0_alex = M0_alex.theory.integrate_siegert()
+
+network_params = {}
+theory_params = {'initial_rates': r_alex_loaded[:, -1],
+                 'T': 50.}
+M0_alex = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
+p, r_alex = M0_alex.theory.integrate_siegert()
+
+
+# p, r0_py = M0.theory.integrate_siegert_python(parallel=False)
+# p, r_py = M.theory.integrate_siegert_python(parallel=False)
 
 assert(np.allclose(mu, mu0))
 assert(np.allclose(sigma, sigma0))

tests/test_replace_cc.py

@@ -16,7 +16,7 @@ def test_het_poisson_stat_mf():
     network_params = {}
     theory_params = {}
     M0 = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
-    p, r0 = M0.theory.integrate_siegert_nest()
+    p, r0 = M0.theory.integrate_siegert()
 
     rates = vector_to_dict(r0[:, -1], M0.area_list, M0.structure)
     with open('mf_rates.json', 'w') as f:
@@ -26,7 +26,7 @@ def test_het_poisson_stat_mf():
                                             'replace_cc_input_source': 'mf_rates.json'}}
     theory_params = {}
     M = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
-    p, r = M.theory.integrate_siegert_nest()
+    p, r = M.theory.integrate_siegert()
 
     assert(np.allclose(r0[:, -1], r[:, -1]))
 
@@ -35,7 +35,7 @@ def test_hom_poisson_stat_mf():
     network_params = {'connection_params': {'replace_cc': 'hom_poisson_stat'}}
     theory_params = {}
     M = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
-    p, r = M.theory.integrate_siegert_nest()
+    p, r = M.theory.integrate_siegert()
 
     mu, sigma = M.theory.replace_cc_input()
     # Test for V1