Merge branch 'master' into fix/aln_input_index

neurolib/models/aln/loadDefaultParams.py

@@ -187,27 +187,6 @@ def loadDefaultParams(Cmat=None, Dmat=None, lookupTableFileName=None, seed=None)
     return params
 
 
-def computeDelayMatrix(lengthMat, signalV, segmentLength=1):
-    """
-    Compute the delay matrix from the fiber length matrix and the signal
-    velocity
-
-        :param lengthMat:       A matrix containing the connection length in
-            segment
-        :param signalV:         Signal velocity in m/s
-        :param segmentLength:   Length of a single segment in mm
-
-        :returns:    A matrix of connexion delay in ms
-    """
-
-    normalizedLenMat = lengthMat * segmentLength
-    if signalV > 0:
-        Dmat = normalizedLenMat / signalV  # Interareal delays in ms
-    else:
-        Dmat = lengthMat * 0.0
-    return Dmat
-
-
 def generateRandomICs(N, seed=None):
     """Generates random Initial Conditions for the interareal network
 

neurolib/models/aln/timeIntegration.py

@@ -1,7 +1,6 @@
 import numpy as np
 import numba
 
-from . import loadDefaultParams as dp
 from ...utils import model_utils as mu
 
 
@@ -53,7 +52,7 @@ def timeIntegration(params):
     if N == 1:
         Dmat = np.ones((N, N)) * params["de"]
     else:
-        Dmat = dp.computeDelayMatrix(
+        Dmat = mu.computeDelayMatrix(
             lengthMat, signalV
         )  # Interareal connection delays, Dmat(i,j) Connnection from jth node to ith (ms)
         Dmat[np.eye(len(Dmat)) == 1] = np.ones(len(Dmat)) * params["de"]
@@ -514,7 +513,7 @@ def timeIntegration_njit_elementwise(
             mufi_rhs = (mui - mufi[no]) / tau_inh
 
             # rate has to be kHz
-            IA_rhs = (a * (Vmean_exc - EA) - IA[no, i - 1] + tauA * b * rates_exc[no, i] * 1e-3) / tauA
+            IA_rhs = (a * (Vmean_exc - EA) - IA[no, i - 1] + tauA * b * rates_exc[no, i-1] * 1e-3) / tauA
 
             # EQ. 4.43
             if distr_delay:

neurolib/models/fhn/loadDefaultParams.py

@@ -79,22 +79,3 @@ def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
     params.y_ext = np.zeros((params.N,))
 
     return params
-
-
-def computeDelayMatrix(lengthMat, signalV, segmentLength=1):
-    """Compute the delay matrix from the fiber length matrix and the signal velocity
-
-    :param lengthMat:       A matrix containing the connection length in segment
-    :param signalV:         Signal velocity in m/s
-    :param segmentLength:   Length of a single segment in mm
-
-    :returns:    A matrix of connexion delay in ms
-    """
-
-    normalizedLenMat = lengthMat * segmentLength
-    # Interareal connection delays, Dmat(i,j) in ms
-    if signalV > 0:
-        Dmat = normalizedLenMat / signalV
-    else:
-        Dmat = lengthMat * 0.0
-    return Dmat

neurolib/models/fhn/timeIntegration.py

@@ -1,13 +1,12 @@
 import numpy as np
 import numba
 
-from . import loadDefaultParams as dp
 from ...utils import model_utils as mu
 
 
 def timeIntegration(params):
     """Sets up the parameters for time integration
-    
+
     :param params: Parameter dictionary of the model
     :type params: dict
     :return: Integrated activity variables of the model
@@ -53,11 +52,10 @@ def timeIntegration(params):
         Dmat = np.zeros((N, N))
     else:
         # Interareal connection delays, Dmat(i,j) Connnection from jth node to ith (ms)
-        Dmat = dp.computeDelayMatrix(lengthMat, signalV)
+        Dmat = mu.computeDelayMatrix(lengthMat, signalV)
         # no self-feedback delay
         Dmat[np.eye(len(Dmat)) == 1] = np.zeros(len(Dmat))
     Dmat_ndt = np.around(Dmat / dt).astype(int)  # delay matrix in multiples of dt
-    params["Dmat_ndt"] = Dmat_ndt
 
     # Additive or diffusive coupling scheme
     coupling = params["coupling"]
@@ -80,8 +78,8 @@ def timeIntegration(params):
     max_global_delay = np.max(Dmat_ndt)
     startind = int(max_global_delay + 1)  # timestep to start integration at
 
-    x_ou = params["x_ou"]
-    y_ou = params["y_ou"]
+    x_ou = params["x_ou"].copy()
+    y_ou = params["y_ou"].copy()
 
     # state variable arrays, have length of t + startind
     # they store initial conditions AND simulated data
@@ -229,13 +227,13 @@ def timeIntegration_njit_elementwise(
                 - ys[no, i - 1]
                 + xs_input_d[no]  # input from other nodes
                 + x_ou[no]  # ou noise
-                + x_ext[no, i-1] # external input
+                + x_ext[no, i - 1]  # external input
             )
             y_rhs = (
                 (xs[no, i - 1] - delta - epsilon * ys[no, i - 1]) / tau
                 + ys_input_d[no]  # input from other nodes
                 + y_ou[no]  # ou noise
-                + y_ext[no, i-1] # external input
+                + y_ext[no, i - 1]  # external input
             )
 
             # Euler integration

neurolib/models/hopf/loadDefaultParams.py

@@ -77,22 +77,3 @@ def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
     params.y_ext = np.zeros((params.N,))
 
     return params
-
-
-def computeDelayMatrix(lengthMat, signalV, segmentLength=1):
-    """Compute the delay matrix from the fiber length matrix and the signal velocity
-
-    :param lengthMat:       A matrix containing the connection length in segment
-    :param signalV:         Signal velocity in m/s
-    :param segmentLength:   Length of a single segment in mm
-
-    :returns:    A matrix of connexion delay in ms
-    """
-
-    normalizedLenMat = lengthMat * segmentLength
-    # Interareal connection delays, Dmat(i,j) in ms
-    if signalV > 0:
-        Dmat = normalizedLenMat / signalV
-    else:
-        Dmat = lengthMat * 0.0
-    return Dmat

neurolib/models/hopf/timeIntegration.py

@@ -1,13 +1,12 @@
 import numpy as np
 import numba
 
-from . import loadDefaultParams as dp
 from ...utils import model_utils as mu
 
 
 def timeIntegration(params):
     """Sets up the parameters for time integration
-    
+
     :param params: Parameter dictionary of the model
     :type params: dict
     :return: Integrated activity variables of the model
@@ -59,10 +58,9 @@ def timeIntegration(params):
         Dmat = np.zeros((N, N))
     else:
         # Interareal connection delays, Dmat(i,j) Connnection from jth node to ith (ms)
-        Dmat = dp.computeDelayMatrix(lengthMat, signalV)
+        Dmat = mu.computeDelayMatrix(lengthMat, signalV)
         Dmat[np.eye(len(Dmat)) == 1] = np.zeros(len(Dmat))
     Dmat_ndt = np.around(Dmat / dt).astype(int)  # delay matrix in multiples of dt
-    params["Dmat_ndt"] = Dmat_ndt
     # ------------------------------------------------------------------------
 
     # Initialization
@@ -74,8 +72,8 @@ def timeIntegration(params):
     max_global_delay = np.max(Dmat_ndt)
     startind = int(max_global_delay + 1)  # timestep to start integration at
 
-    x_ou = params["x_ou"]
-    y_ou = params["y_ou"]
+    x_ou = params["x_ou"].copy()
+    y_ou = params["y_ou"].copy()
 
     # state variable arrays, have length of t + startind
     # they store initial conditions AND simulated data
@@ -208,14 +206,14 @@ def timeIntegration_njit_elementwise(
                 - w * ys[no, i - 1]
                 + xs_input_d[no]  # input from other nodes
                 + x_ou[no]  # ou noise
-                + x_ext[no, i-1]  # external input
+                + x_ext[no, i - 1]  # external input
             )
             y_rhs = (
                 (a - xs[no, i - 1] ** 2 - ys[no, i - 1] ** 2) * ys[no, i - 1]
                 + w * xs[no, i - 1]
                 + ys_input_d[no]  # input from other nodes
                 + y_ou[no]  # ou noise
-                + y_ext[no, i-1]  # external input
+                + y_ext[no, i - 1]  # external input
             )
 
             # Euler integration

neurolib/models/kuramoto/__init__.py

@@ -0,0 +1 @@
+from .model import KuramotoModel

neurolib/models/kuramoto/loadDefaultParams.py

@@ -0,0 +1,62 @@
+import numpy as np
+
+from neurolib.utils.collections import dotdict
+
+
+def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
+    """Load default parameters for the Kuramoto Model model
+
+    :param Cmat: Structural connectivity matrix (adjacency matrix) of coupling strengths, will be normalized to 1. If not given, then a single node simulation will be assumed, defaults to None
+    :type Cmat: numpy.ndarray, optional
+    :param Dmat: Fiber length matrix, will be used for computing the delay matrix together with the signal transmission speed parameter `signalV`, defaults to None
+    :type Dmat: numpy.ndarray, optional
+    :param seed: Seed for the random number generator, defaults to None
+    :type seed: int, optional
+
+    :return: A dictionary with the default parameters of the model
+    :rtype: dict
+    """
+    params = dotdict({})
+
+    ### runtime parameters
+
+    params.dt = 0.1 
+    params.duration = 2000  
+
+    np.random.seed(seed)  
+    params.seed = seed
+
+    # model parameters
+    params.N = 1
+    params.k = 2
+
+    # connectivity
+    if Cmat is None:
+        params.N = 1
+        params.Cmat = np.zeros((1, 1))
+        params.lengthMat = np.zeros((1, 1))
+    else:
+        params.Cmat = Cmat.copy()  # coupling matrix
+        np.fill_diagonal(params.Cmat, 0)  # no self connections
+        params.N = len(params.Cmat)  # override number of nodes
+        params.lengthMat = Dmat
+
+    params.omega = np.ones((params.N,)) * np.pi
+
+    params.signalV = 20.0
+
+    # Ornstein-Uhlenbeck process
+    params.tau_ou = 5.0  # ms Timescale of the Ornstein-Uhlenbeck noise process
+    params.sigma_ou = 0.0  # 1/ms/sqrt(ms) noise intensity
+
+    # init values
+    params.theta_init = np.random.uniform(low=0, high=2*np.pi, size=(params.N, 1))
+
+    # Ornstein-Uhlenbeck process
+    params.theta_ou = np.zeros((params.N,))
+
+    # external input
+    params.theta_ext = np.zeros((params.N,))
+
+    return params
+

neurolib/models/kuramoto/model.py

@@ -0,0 +1,35 @@
+from . import loadDefaultParams as dp
+from . import timeIntegration as ti
+
+from neurolib.models.model import Model
+
+
+class KuramotoModel(Model):
+    """
+    Kuramoto Model
+
+    Based on:
+    Kuramoto, Yoshiki (1975). H. Araki (ed.). Lecture Notes in Physics, International Symposium on Mathematical Problems in Theoretical Physics.
+    """
+
+    name = "kuramoto"
+    description = "Kuramoto Model"
+
+    init_vars = ['theta_init', 'theta_ou']
+    state_vars = ['theta', 'theta_ou'] 
+    output_vars = ['theta']
+    default_output = 'theta'
+    input_vars = ['theta_ext']
+    default_input = 'theta_ext'
+
+    def __init__(self, params=None, Cmat=None, Dmat=None, seed=None):
+        self.Cmat = Cmat
+        self.Dmat = Dmat
+        self.seed = seed
+
+        integration = ti.timeIntegration
+
+        if params is None:
+            params = dp.loadDefaultParams(Cmat=self.Cmat, Dmat=self.Dmat, seed=self.seed)
+
+        super().__init__(params=params, integration=integration)