Switch CTRNN integration from forward Euler to exponential Euler (ETD1)

The exponential Euler method integrates the linear decay term -u/tau
exactly, making it unconditionally stable regardless of dt/tau ratio.
Forward Euler required dt < 2*tau for stability. The nonlinear forcing
term is still held constant per step (same assumption as before).

Update formula: u_new = exp(-dt/tau) * u_old + (1 - exp(-dt/tau)) * z
Test reference values regenerated for the new integrator.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: CodeReclaimers

neat/ctrnn/__init__.py

@@ -1,5 +1,6 @@
 """Handles the continuous-time recurrent neural network implementation."""
 
+import math
 
 from neat.graphs import required_for_output
 
@@ -78,7 +79,8 @@ def advance(self, inputs, advance_time, time_step=None):
                 node_inputs = [ivalues[i] * w for i, w in ne.links]
                 s = ne.aggregation(node_inputs)
                 z = ne.activation(ne.bias + ne.response * s)
-                ovalues[node_key] += dt / ne.time_constant * (-ovalues[node_key] + z)
+                decay = math.exp(-dt / ne.time_constant)
+                ovalues[node_key] = decay * ovalues[node_key] + (1.0 - decay) * z
 
             self.time_seconds += dt
 

tests/test_ctrnn.py

@@ -48,12 +48,15 @@ def test_basic_two_neuron_dynamics():
 
     # Check specific reference values at selected timesteps to guard against
     # regressions in the CTRNN integration behavior.
+    # These values use exponential Euler (ETD1) integration, which integrates
+    # the linear decay term exactly: u(t+h) = decay*u(t) + (1-decay)*z
+    # where decay = exp(-h/tau).
     reference = {
         0: (0.0, 0.0),
-        100: (0.3746852284, 0.8115273872),
-        500: (0.5634208426, 0.1952743492),
-        1149: (0.4092483263, 0.8024978195),
-        1249: (0.7531862678, 0.3112381247),
+        100: (0.3087387628, 0.7906422519),
+        500: (0.7965631789, 0.4505240611),
+        1149: (0.4952237403, 0.1974657233),
+        1249: (0.2539079366, 0.7016150055),
     }
     tol = 1e-6
     for idx, (exp0, exp1) in reference.items():