Assert side effects in KINETIC blocks. (#1392)

Effectively the side effects in KINETIC blocks, such as assigning a
STATE to a RANGE variable happen (once more) after the state has been
updated.

This commit adds test for this; and "fixes" the bug by running the
functor method `initialize` again after running Newton.

src/codegen/codegen_cpp_visitor.cpp

@@ -1087,6 +1087,7 @@ void CodegenCppVisitor::visit_eigen_newton_solver_block(const ast::EigenNewtonSo
 
     // assign newton solver results in matrix X to state vars
     print_statement_block(*node.get_update_states_block(), false, false);
+    printer->add_line("newton_functor.initialize(); // TODO mimic calling F again.");
     printer->add_line("newton_functor.finalize();");
 }
 

test/usecases/kinetic/side_effects.mod

@@ -0,0 +1,36 @@
+NEURON {
+    SUFFIX side_effects
+    RANGE x, forward_flux, backward_flux
+    NONSPECIFIC_CURRENT il
+}
+
+ASSIGNED {
+    il
+    x
+    forward_flux
+    backward_flux
+}
+
+STATE {
+    X
+    Y
+}
+
+INITIAL {
+    X = 1
+    Y = 2
+}
+
+BREAKPOINT {
+    SOLVE state METHOD sparse
+
+    il = forward_flux - backward_flux
+}
+
+KINETIC state {
+    ~ X <-> Y (0.4, 0.5)
+    forward_flux = f_flux
+    backward_flux = b_flux
+
+    x = X
+}

test/usecases/kinetic/test_side_effects.py

@@ -0,0 +1,59 @@
+import numpy as np
+
+from neuron import h, gui
+from neuron.units import ms
+
+
+def run_simulation():
+    s = h.Section()
+    s.insert("side_effects")
+
+    mech = s(0.5).side_effects
+    t_hoc = h.Vector().record(h._ref_t)
+    x_hoc = h.Vector().record(mech._ref_x)
+    X_hoc = h.Vector().record(mech._ref_X)
+    Y_hoc = h.Vector().record(mech._ref_Y)
+    forward_flux_hoc = h.Vector().record(mech._ref_forward_flux)
+    backward_flux_hoc = h.Vector().record(mech._ref_backward_flux)
+
+    h.stdinit()
+    h.tstop = 5.0 * ms
+    h.run()
+
+    timeseries = {
+        "t": np.array(t_hoc.as_numpy()),
+        "x": np.array(x_hoc.as_numpy()),
+        "X": np.array(X_hoc.as_numpy()),
+        "Y": np.array(Y_hoc.as_numpy()),
+        "forward_flux": np.array(forward_flux_hoc.as_numpy()),
+        "backward_flux": np.array(backward_flux_hoc.as_numpy()),
+    }
+
+    return timeseries
+
+
+def check_assignment(x, X):
+    # At time t = 0, the side effects aren't applied.
+    np.testing.assert_array_equal(x[1:], X[1:])
+
+
+def check_flux(actual_flux, expected_flux):
+    # At time t = 0, the side effects aren't applied.
+    np.testing.assert_array_almost_equal_nulp(
+        actual_flux[1:], expected_flux[1:], nulp=8
+    )
+
+
+def check_forward_flux(X, actual_flux):
+    check_flux(actual_flux, 0.4 * X)
+
+
+def check_backward_flux(Y, actual_flux):
+    check_flux(actual_flux, 0.5 * Y)
+
+
+if __name__ == "__main__":
+    timeseries = run_simulation()
+    check_assignment(timeseries["x"], timeseries["X"])
+    check_forward_flux(timeseries["X"], timeseries["forward_flux"])
+    check_backward_flux(timeseries["Y"], timeseries["backward_flux"])