- debugged the implementation of the kuramoto oscillators in `model_t…

…emplates`

- added a simple test for correct behavior of the kuramoto model implementations
- added a kuramoto oscillator version with noise to `model_templates`

model_templates/kuramoto/simple_kuramoto.yml

@@ -1,50 +1,93 @@
 %YAML 1.2
 ---
 
-Op_phase:
+Op_noise:
   base: OperatorTemplate
   equations:
-    - "d/dt * theta = omega + noise_in*sqrt(step_size) + K*net_in"
-    - "theta_wrapped = theta % (2*PI)"
+    - "t_old = t_new"
+    - "t_new = t"
+    - "d/dt * theta = omega + xi*sqrt(t_new-t_old) + K*net_in"
   variables:
     theta:
       default: output
-    theta_wrapped:
-      default: variable
     omega:
       default: 10.0
     K:
       default: 1.0
-    noise_in:
+    t_old:
+      default: 0.0
+    t_new:
+      default: 0.0
+    xi:
       default: input
     net_in:
       default: input
 
+Op_base:
+  base: OperatorTemplate
+  equations:
+    - "d/dt * theta = omega + K*net_in"
+  variables:
+    theta:
+      default: output
+    omega:
+      default: 10.0
+    K:
+      default: 1.0
+    net_in:
+      default: input
+
 Op_coupling:
   base: OperatorTemplate
   equations:
-    - "target_var = sin(source_var-target_var)"
+    - "k = sin(theta_s-theta_t)"
   variables:
-    source_var:
+    theta_s:
+      default: input
+    theta_t:
       default: input
-    target_var:
+    k:
       default: output
 
 Sin_edge:
   base: EdgeTemplate
   operators:
     - Op_coupling
 
-KM_node:
+KM:
   base: NodeTemplate
   operators:
-    - Op_phase
+    - Op_base
+
+KM_noise:
+  base: NodeTemplate
+  operators:
+    - Op_noise
+
+KM_single:
+  base: CircuitTemplate
+  nodes:
+    p1: KM
+
+KM_single_noise:
+  base: CircuitTemplate
+  nodes:
+    p1: KM_noise
+
+KMN:
+  base: CircuitTemplate
+  nodes:
+    p1: KM
+    p2: KM
+  edges:
+    - [{p1: {Op_base/theta: theta_s}, p2: {Op_base/theta: theta_t}}, p2/Op_base/net_in, Sin_edge, {weight: 2.0}]
+    - [{p2: {Op_base/theta: theta_s}, p1: {Op_base/theta: theta_t}}, p1/Op_base/net_in, Sin_edge, {weight: -2.0}]
 
-2KM:
+KMN_noise:
   base: CircuitTemplate
   nodes:
-    O1: KM_node
-    O2: KM_node
+    p1: KM_noise
+    p2: KM_noise
   edges:
-    - [O1/Op_phase/theta, O2/Op_phase/net_in, Sin_edge, {weight: 1.0}]
-    - [O2/Op_phase/theta, O1/Op_phase/net_in, Sin_edge, {weight: 1.0}]
+    - [{p1: {Op_noise/theta: theta_s}, p2: {Op_noise/theta: theta_t}}, p2/Op_noise/net_in, Sin_edge, {weight: 4.0}]
+    - [{p2: {Op_noise/theta: theta_s}, p1: {Op_noise/theta: theta_t}}, p1/Op_noise/net_in, Sin_edge, {weight: -4.0}]

tests/test_implemented_models.py

@@ -165,3 +165,69 @@ def test_3_3_wilson_cowan():
     indices = [np.argmin(np.abs(t - r2.index)) for t in times]
     assert r2.iloc[indices[0], 0] < r2.iloc[indices[1], 0]
     assert r2.iloc[indices[0], 0] - r2.iloc[indices[2], 0] == pytest.approx(0., rel=1e-3, abs=1e-3)
+
+
+def test_3_4_kuramoto():
+    """Tests accurate behavior of kuramoto oscillator model.
+    """
+
+    # assess correct response of single base oscillator
+    ###################################################
+
+    T = 2.0
+    dt = 1e-4
+    dts = 1e-2
+
+    # set up circuit
+    km1 = CircuitIR.from_yaml("model_templates.kuramoto.simple_kuramoto.KM_single")
+    km1 = km1.compile(vectorization=True, backend='numpy', solver='scipy', step_size=dt)
+
+    # perform simulation
+    r1 = km1.run(T, sampling_step_size=dts, outputs={"theta": "p1/Op_base/theta"})
+    km1.clear()
+
+    # test linear oscillator properties
+    omega = 10.0
+    results = np.sin(r1.values[:, 0]*2*np.pi)
+    target = np.sin(omega*2.0*np.pi*r1.index.values)
+    assert results - target == pytest.approx(0., rel=1e-2, abs=1e-2)
+
+    # assess correct response of two coupled oscillators
+    ####################################################
+
+    T = 6.0
+    dt = 1e-4
+    dts = 1e-2
+
+    # set up circuit
+    km2 = CircuitIR.from_yaml("model_templates.kuramoto.simple_kuramoto.KMN")
+    km2 = km2.compile(vectorization=True, backend='numpy', solver='scipy', step_size=dt)
+
+    # perform simulation
+    r2 = km2.run(T, sampling_step_size=dts, outputs={"theta1": "p1/Op_base/theta", "theta2": "p2/Op_base/theta"})
+    km2.clear()
+
+    # test whether oscillator 2 showed a faster phase development than oscillator 1
+    assert r2['theta1'].iloc[-1, 0] < r2['theta2'].iloc[-1, 0]
+
+    # repeat test 2 for two coupled noisy oscillators
+    #################################################
+
+    T = 6.0
+    dt = 1e-4
+    dts = 1e-2
+
+    inp1 = np.random.randn(int(np.round(T / dt))) * 0.5
+    inp2 = np.random.randn(int(np.round(T / dt))) * 0.1
+
+    # set up circuit
+    km3 = CircuitIR.from_yaml("model_templates.kuramoto.simple_kuramoto.KMN_noise")
+    km3 = km3.compile(vectorization=True, backend='numpy', solver='scipy', step_size=dt)
+
+    # perform simulation
+    r3 = km3.run(T, sampling_step_size=dts, outputs={"theta1": "p1/Op_noise/theta", "theta2": "p2/Op_noise/theta"},
+                 inputs={"p1/Op_noise/xi": inp1, "p2/Op_noise/xi": inp2})
+    km3.clear()
+
+    # test whether oscillator 2 showed a faster phase development than oscillator 1
+    assert r3['theta1'].iloc[-1, 0] < r3['theta2'].iloc[-1, 0]