Merge remote-tracking branch 'origin/master' into qa/clean-up-load-ba…

…lance

.github/workflows/test-spack.yml

@@ -38,8 +38,8 @@ jobs:
       - name: Get Spack (latest_release)
         if: ${{ matrix.spack-version == 'latest_release' }}
         run: |
-          wget "$(curl -sH "Accept: application/vnd.github.v3+json" https://api.github.com/repos/spack/spack/releases/latest | grep browser_download_url |  cut -d '"' -f 4)"
-          tar xfz spack*.tar.gz
+          wget -O latest_spack.tar.gz "$(curl -sH "Accept: application/vnd.github.v3+json" https://api.github.com/repos/spack/spack/releases/latest | grep tarball_url |  cut -d '"' -f 4)"
+          tar xfz latest_spack.tar.gz
           mv spack*/ spack
       - name: Prep
         run: |

doc/fileformat/nmodl.rst

@@ -84,9 +84,7 @@ Ions
   not made ``STATE`` they may be written to, but their values will be reset to
   their initial values every time step.
 * The diffusive concentration ``Xd`` does not share this semantics. It will not
-  be reset, even if not in ``STATE``, and may freely be written. This comes at the
-  cost of awkward treatment of ODEs for ``Xd``, see the included ``decay.mod`` for
-  an example.
+  be reset, even if not in ``STATE``, and may freely be written.
 * ``Xd`` is present on all cables iff its associated diffusivity is set to a
   non-zero value.
 
@@ -141,7 +139,66 @@ Unsupported features
 * ``LOCAL`` variables outside blocks are not supported.
 * free standing blocks are not supported.
 * ``INDEPENDENT`` variables are not supported.
-* arrays and pointers are not supported by Arbor.
+* loops, arrays, and pointers are not supported by Arbor.
+
+Alternating normal and reaction statements in KINETIC
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This is not so much a feature as a bug that Arbor's ``modcc`` does not reproduce
+from NEURON. Given a file like:
+
+.. code-block::
+
+   NEURON { SUFFIX test_kinetic_alternating_reaction }
+
+   STATE { A B C D }
+
+   BREAKPOINT {
+     SOLVE foobar METHOD sparse
+   }
+
+   KINETIC foobar {
+     LOCAL x, y
+
+     x = 23*v
+     y = 42*v
+
+     ~ A <-> B (x, y)
+
+     x = sin(y)
+     y = cos(x)
+
+     ~ C <-> D (x, y)
+   }
+
+one might expect that the reaction between ``C`` and ``D`` occurs with a rate
+proportional to ``sin(x)`` and ``cos(y)``. However, this file is equivalent to
+
+.. code-block::
+
+   NEURON { SUFFIX test_kinetic_alternating_reaction }
+
+   STATE { A B C D }
+
+   BREAKPOINT {
+     SOLVE foobar METHOD sparse
+   }
+
+   KINETIC foobar {
+     LOCAL x, y
+
+     x = 23*v
+     y = 42*v
+     x = sin(y)
+     y = cos(x)
+
+     ~ A <-> B (x, y)
+     ~ C <-> D (x, y)
+   }
+
+which is almost never what the author would expect. Thus, this construction constitutes
+an error in ``modcc``; if want this particular behaviour, please state this explicit by
+writing code similar to the second example.
 
 .. _arbornmodl:
 

mechanisms/default/decay.mod

@@ -1,21 +1,15 @@
 NEURON {
     SUFFIX decay
     USEION x WRITE xd
-    RANGE F, tau
+    RANGE tau
 }
 
 PARAMETER { tau = 5 }
 
-INITIAL { F = xd }
-
-STATE { F }
-
 BREAKPOINT {
-   SOLVE dF METHOD cnexp
-   xd = F
+   SOLVE dX METHOD cnexp
 }
 
-DERIVATIVE dF {
-   F = xd
-   F' = -tau*F
+DERIVATIVE dX {
+   xd' = -tau*xd
 }

modcc/expression.cpp

@@ -129,7 +129,6 @@ void IdentifierExpression::semantic(scope_ptr scp) {
             return;
         }
     }
-
     // save the symbol
     symbol_ = s;
 }
@@ -156,18 +155,26 @@ expression_ptr DerivativeExpression::clone() const {
 }
 
 void DerivativeExpression::semantic(scope_ptr scp) {
+    // Check for semantic errors first
     error_ = false;
-
     IdentifierExpression::semantic(scp);
-    if (has_error()) {
-        return;
-    }
-    auto v = symbol_->is_variable();
-    if (!v || !v->is_state()) {
-        error( pprintf("the variable '%' must be a state variable to be differentiated",
-                        yellow(spelling_), location_));
-        return;
+    if (has_error()) return;
+
+    // STATE is ok to take derivatives of.
+    if (auto var = symbol_->is_variable(); var && var->is_state()) return;
+
+    // Diffusive concentrations may also be differentiated
+    if (auto local = symbol_->is_local_variable()) {
+        if (auto ext = local->external_variable();
+            ext && ext->data_source() == sourceKind::ion_diffusive) {
+            return;
+        }
     }
+
+    // Anything else raises an error
+    error(pprintf("The variable '%' must be a state variable or diffusive concentration to be differentiated.",
+                  yellow(spelling_), location_));
+
 }
 
 /*******************************************************************************
@@ -555,6 +562,8 @@ void ProcedureExpression::semantic(scope_ptr scp) {
     error_ = false;
     scope_ = scp;
 
+    auto is_kinetic = [](const auto& it) { return (it->is_reaction() || it->is_conserve()); };
+
     // assert that the symbol is already visible in the global_symbols
     if(scope_->find_global(name()) == nullptr) {
         throw compiler_exception(
@@ -571,8 +580,34 @@ void ProcedureExpression::semantic(scope_ptr scp) {
 
     // this loop could be used to then check the types of statements in the body
     for(auto& e : *(body_->is_block())) {
-        if(e->is_initial_block())
+        if (e->is_initial_block()) {
             error("INITIAL block not allowed inside "+::to_string(kind_)+" definition");
+        }
+        if (kind_ != procedureKind::kinetic && is_kinetic(e)) {
+            error("reaction statement not allowed inside " + ::to_string(kind_)+" definition");
+        }
+        if (kind_ != procedureKind::linear && e->is_linear()) {
+            error("linear statement not allowed inside "+::to_string(kind_)+" definition");
+        }
+        if (kind_ != procedureKind::linear && e->is_derivative()) {
+            error("derivative statement not allowed inside "+::to_string(kind_)+" definition");
+        }
+    }
+
+    // We start a new loop here for preserving our sanity
+    if (kind_ == procedureKind::kinetic) {
+        auto it = body_->is_block()->begin();
+        auto end = body_->is_block()->end();
+        // skip all 'normal' statements
+        for (; it != end && !is_kinetic(*it); ++it) {}
+        // skip all 'reaction' statements
+        for (; it != end && is_kinetic(*it); ++it) {}
+        // We have trailing 'normal' statements
+        if (it != end) {
+            error("Found alternating reaction (A <-> B ...) and normal statements. "
+                  "This is allowed by NMODL, but likely not what you want; see: "
+                  "https://docs.arbor-sim.org/en/latest/fileformat/nmodl.html#unsupported-features");
+        }
     }
 
     // perform semantic analysis for each expression in the body

modcc/kineticrewriter.cpp

@@ -26,6 +26,7 @@ class KineticRewriter : public BlockRewriterBase {
     virtual void finalize() override;
 
 private:
+
     // Acccumulated terms for derivative expressions, keyed by id name.
     std::map<std::string, expression_ptr> dterms;
 };

modcc/printer/cprinter.cpp

@@ -308,7 +308,7 @@ ARB_LIBMODCC_API std::string emit_cpp_source(const Module& module_, const printe
     out << popindent << "}\n\n";
 
     out << "static void advance_state(arb_mechanism_ppack* pp) {\n" << indent;
-    emit_body(state_api);
+    emit_body(state_api, true);
     out << popindent << "}\n\n";
 
     out << "static void compute_currents(arb_mechanism_ppack* pp) {\n" << indent;

modcc/solvers.cpp

@@ -68,14 +68,13 @@ void CnexpSolverVisitor::visit(AssignmentExpression *e) {
         return;
     }
     else if (r.is_homogeneous) {
-        // s' = a*s becomes s = s*exp(a*dt); use a_ as a local variable
-        // for the coefficient.
+        // s' = a*s becomes s = s*exp(a*dt); use a_ as a local variable for the
+        // coefficient.
         auto local_a_term = make_unique_local_assign(scope, coef, "a_");
+        auto s_update = pprintf("% = %*exp_pade_11(%*dt)",
+                                s, s, local_a_term.id->is_identifier()->spelling());
         statements_.push_back(std::move(local_a_term.local_decl));
         statements_.push_back(std::move(local_a_term.assignment));
-        auto a_ = local_a_term.id->is_identifier()->spelling();
-
-        std::string s_update = pprintf("% = %*exp_pade_11(%*dt)", s, s, a_);
         statements_.push_back(Parser{s_update}.parse_line_expression());
         return;
     }
@@ -1054,7 +1053,7 @@ class UnusedVisitor : public Visitor {
     }
 
     virtual void visit(IdentifierExpression* e) override {
-        if (lhs_sym && lhs_sym->is_local_variable()) {
+        if (lhs_sym && lhs_sym->is_local_variable() && !lhs_sym->is_local_variable()->external_variable()) {
             deps.insert({lhs_sym->name(), e->name()});
         }
         else {

python/example/diffusion.py

@@ -11,7 +11,6 @@ def __init__(self, cell, probes):
         self.the_cell = cell
         self.the_probes = probes
         self.the_props = A.neuron_cable_properties()
-        self.the_props.catalogue = A.default_catalogue()
 
     def num_cells(self):
         return 1
@@ -63,15 +62,14 @@ def event_generators(self, gid):
         # Plot
         for lbl, ix in zip(m, range(1, d.shape[1])):
             ax.plot(d[:, 0], d[:, ix], label=lbl)
+        W = 8
         # Table
         print("Sodium concentration (NaD/mM)")
-        print("|-" + "-+-".join("-" * 20 for _ in range(d.shape[1])) + "-|")
-        print(
-            "| Time (ms)            | " + " | ".join(f"{str(l):<20}" for l in m) + " |"
-        )
-        print("|-" + "-+-".join("-" * 20 for _ in range(d.shape[1])) + "-|")
+        print("|-" + "-+-".join("-" * W for _ in range(d.shape[1])) + "-|")
+        print(f"| {'Time/ms':<{W}} | " + " | ".join(f"{l.prox:<{W}}" for l in m) + " |")
+        print("|-" + "-+-".join("-" * W for _ in range(d.shape[1])) + "-|")
         for ix in range(d.shape[0]):
-            print("| " + " | ".join(f"{v:>20.3f}" for v in d[ix, :]) + " |")
-        print("|-" + "-+-".join("-" * 20 for _ in range(d.shape[1])) + "-|")
+            print("| " + " | ".join(f"{v:>{W}.3f}" for v in d[ix, :]) + " |")
+        print("|-" + "-+-".join("-" * W for _ in range(d.shape[1])) + "-|")
 ax.legend()
 fg.savefig("results.pdf")

python/example/probe_lfpykit.py

@@ -392,10 +392,6 @@ def colorbar(
 # show predictions at the last time point of simulation
 time_index = -1
 
-# use seaborn style
-plt.style.use("seaborn")
-plt.rcParams["image.cmap"] = "viridis"
-
 # create figure and axis
 fig, ax = plt.subplots(1, 1, figsize=(16, 6), dpi=100)
 

python/recipe.cpp

@@ -188,7 +188,8 @@ void register_recipe(pybind11::module& m) {
         .def("cell_description", &py_recipe::cell_description, pybind11::return_value_policy::copy,
             "gid"_a,
             "High level description of the cell with global identifier gid.")
-        .def("cell_kind", &py_recipe::cell_kind,
+        .def("cell_kind",
+             &py_recipe::cell_kind,
             "gid"_a,
             "The kind of cell with global identifier gid.")
         .def("event_generators", &py_recipe::event_generators,

python/recipe.hpp

@@ -64,7 +64,12 @@ class py_recipe_trampoline: public py_recipe {
     }
 
     arb::cell_kind cell_kind(arb::cell_gid_type gid) const override {
-        PYBIND11_OVERRIDE_PURE(arb::cell_kind, py_recipe, cell_kind, gid);
+        try {
+            PYBIND11_OVERRIDE_PURE(arb::cell_kind, py_recipe, cell_kind, gid);
+        }
+        catch (const pybind11::cast_error& e) {
+            throw pybind11::type_error{"Couldn't convert return value of recipe.cell_kind(gid) to a cell_kind. Please check your recipe."};
+        }
     }
 
     std::vector<pybind11::object> event_generators(arb::cell_gid_type gid) const override {

test/unit-modcc/mod_files/test_alternating_derivative.mod

@@ -0,0 +1,21 @@
+NEURON { SUFFIX test_derivative_alternating }
+
+STATE { A B }
+
+BREAKPOINT { : NOTE we need a newline here :/
+  SOLVE foobar METHOD sparse
+}
+
+: NOTE this is OK to do, while alternating normal and reaction statements is not!
+
+DERIVATIVE foobar {
+  LOCAL x
+
+  x = 23*v
+
+  A' = B*x
+
+  x = sin(x)
+
+  B' = x*A
+}

test/unit-modcc/mod_files/test_alternating_differential.mod

@@ -0,0 +1,18 @@
+NEURON { SUFFIX test_kinetic_alternating_differential }
+
+STATE { A B C D }
+
+BREAKPOINT { : NOTE we need a newline here :/
+  SOLVE foobar METHOD sparse
+}
+
+KINETIC foobar {
+  LOCAL x, y
+
+  x = 23*v
+  y = 42*v
+
+  ~ A <-> B (x, y)
+  C' = 0.1
+  ~ C <-> D (x, y)
+}

test/unit-modcc/mod_files/test_alternating_reaction.mod

@@ -0,0 +1,21 @@
+NEURON { SUFFIX test_kinetic_alternating_reaction }
+
+STATE { A B C D }
+
+BREAKPOINT { : NOTE we need a newline here :/
+  SOLVE foobar METHOD sparse
+}
+
+KINETIC foobar {
+  LOCAL x, y
+
+  x = 23*v
+  y = 42*v
+
+  ~ A <-> B (x, y)
+
+  x = sin(y)
+  y = cos(x)
+
+  ~ C <-> D (x, y)
+}

test/unit-modcc/mod_files/test_kinetic_under_conditional.mod

@@ -0,0 +1,21 @@
+NEURON { SUFFIX test_kinetic_under_conditional }
+
+STATE { A B }
+
+BREAKPOINT { : NOTE we need a newline here :/
+  SOLVE foobar METHOD sparse
+}
+
+KINETIC foobar {
+  LOCAL x, y
+
+  x = 23*v
+  y = 42*v
+
+  if (v<0) {
+    ~ A <-> B (x, y)
+  }
+  else {
+    ~ A <-> B (y, x)
+  }
+}