From faeacdd916271b34f8791e2049bdce1c1ccffd9c Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 19 Sep 2023 10:11:13 +0200 Subject: [PATCH 1/5] Treat some corner cases in modcc (#2216) These tests and fixes cover a few cases related to kinetic blocks: # Non-reaction statement in between reaction statements This leads to the 'normal' statements being hoisted to the top resulting in unexpected semantics. This is now an error. # Reaction statements under conditional This has always been an error; add a test for it. # Reaction/Linear statements in non-KINETIC/LINEAR context This is now an error; did SEGV before. # Derivative statements in KINETIC/LINEAR/PROCEDURE These are now error, they were before, but resulted in a confusing ICE message. --- doc/fileformat/nmodl.rst | 61 ++++++++++++++++++- modcc/expression.cpp | 30 ++++++++- modcc/kineticrewriter.cpp | 1 + .../mod_files/test_alternating_derivative.mod | 21 +++++++ .../test_alternating_differential.mod | 18 ++++++ .../mod_files/test_alternating_reaction.mod | 21 +++++++ .../test_kinetic_under_conditional.mod | 21 +++++++ .../mod_files/test_misplaced_reaction.mod | 15 +++++ test/unit-modcc/test_module.cpp | 50 +++++++++++++++ 9 files changed, 236 insertions(+), 2 deletions(-) create mode 100644 test/unit-modcc/mod_files/test_alternating_derivative.mod create mode 100644 test/unit-modcc/mod_files/test_alternating_differential.mod create mode 100644 test/unit-modcc/mod_files/test_alternating_reaction.mod create mode 100644 test/unit-modcc/mod_files/test_kinetic_under_conditional.mod create mode 100644 test/unit-modcc/mod_files/test_misplaced_reaction.mod diff --git a/doc/fileformat/nmodl.rst b/doc/fileformat/nmodl.rst index 4787c719f9..10b6aa0fcc 100644 --- a/doc/fileformat/nmodl.rst +++ b/doc/fileformat/nmodl.rst @@ -141,7 +141,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: diff --git a/modcc/expression.cpp b/modcc/expression.cpp index 82ec3ac679..add6678f7d 100644 --- a/modcc/expression.cpp +++ b/modcc/expression.cpp @@ -555,6 +555,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 +573,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 diff --git a/modcc/kineticrewriter.cpp b/modcc/kineticrewriter.cpp index b3aef7ab9f..6ad60e76b5 100644 --- a/modcc/kineticrewriter.cpp +++ b/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 dterms; }; diff --git a/test/unit-modcc/mod_files/test_alternating_derivative.mod b/test/unit-modcc/mod_files/test_alternating_derivative.mod new file mode 100644 index 0000000000..ee9ce4a02e --- /dev/null +++ b/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 +} diff --git a/test/unit-modcc/mod_files/test_alternating_differential.mod b/test/unit-modcc/mod_files/test_alternating_differential.mod new file mode 100644 index 0000000000..dd7e86f8ea --- /dev/null +++ b/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) +} diff --git a/test/unit-modcc/mod_files/test_alternating_reaction.mod b/test/unit-modcc/mod_files/test_alternating_reaction.mod new file mode 100644 index 0000000000..77c9f7a674 --- /dev/null +++ b/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) +} diff --git a/test/unit-modcc/mod_files/test_kinetic_under_conditional.mod b/test/unit-modcc/mod_files/test_kinetic_under_conditional.mod new file mode 100644 index 0000000000..50932081a4 --- /dev/null +++ b/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) + } +} diff --git a/test/unit-modcc/mod_files/test_misplaced_reaction.mod b/test/unit-modcc/mod_files/test_misplaced_reaction.mod new file mode 100644 index 0000000000..3ab519525d --- /dev/null +++ b/test/unit-modcc/mod_files/test_misplaced_reaction.mod @@ -0,0 +1,15 @@ +NEURON { SUFFIX test_misplaced_reaction } + +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 + + ~ A <-> B (x, 1/x) +} diff --git a/test/unit-modcc/test_module.cpp b/test/unit-modcc/test_module.cpp index d57e6c4120..6517fc0d63 100644 --- a/test/unit-modcc/test_module.cpp +++ b/test/unit-modcc/test_module.cpp @@ -216,3 +216,53 @@ TEST(Module, net_receive) { EXPECT_FALSE(m.semantic()); } + +TEST(Module, kinetic_footgun_alternating_reac) { + Module m(io::read_all(DATADIR "/mod_files/test_alternating_reaction.mod"), "test_alternating_reaction.mod"); + EXPECT_NE(m.buffer().size(), 0u); + + Parser p(m, false); + EXPECT_TRUE(p.parse()); + + EXPECT_FALSE(m.semantic()); +} + +TEST(Module, kinetic_footgun_alternating_diff) { + Module m(io::read_all(DATADIR "/mod_files/test_alternating_differential.mod"), "test_alternating_differential.mod"); + EXPECT_NE(m.buffer().size(), 0u); + + Parser p(m, false); + EXPECT_TRUE(p.parse()); + + EXPECT_FALSE(m.semantic()); +} + +TEST(Module, derivative_alternating) { + Module m(io::read_all(DATADIR "/mod_files/test_alternating_derivative.mod"), "test_alternating_derivative.mod"); + EXPECT_NE(m.buffer().size(), 0u); + + Parser p(m, false); + EXPECT_TRUE(p.parse()); + + EXPECT_TRUE(m.semantic()); +} + +TEST(Module, misplaced_reaction) { + Module m(io::read_all(DATADIR "/mod_files/test_misplaced_reaction.mod"), "test_misplaced_reaction.mod"); + EXPECT_NE(m.buffer().size(), 0u); + + Parser p(m, false); + EXPECT_TRUE(p.parse()); + + EXPECT_FALSE(m.semantic()); +} + +TEST(Module, reaction_under_conditional) { + Module m(io::read_all(DATADIR "/mod_files/test_kinetic_under_conditional.mod"), "test_kinetic_under_conditional.mod"); + EXPECT_NE(m.buffer().size(), 0u); + + Parser p(m, false); + EXPECT_FALSE(p.parse()); + + EXPECT_FALSE(m.semantic()); +} From b0993917b7a054e4c44f98073a61161ad2d6d2d4 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 19 Sep 2023 10:20:07 +0200 Subject: [PATCH 2/5] Rethrow a more informative error if cell kind is not a `cell_kind` (#2219) Fixes #2138. --- python/recipe.cpp | 3 ++- python/recipe.hpp | 7 ++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/python/recipe.cpp b/python/recipe.cpp index 16f78494f8..2649ba6bdf 100644 --- a/python/recipe.cpp +++ b/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, diff --git a/python/recipe.hpp b/python/recipe.hpp index fe7624138a..8462b88138 100644 --- a/python/recipe.hpp +++ b/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 event_generators(arb::cell_gid_type gid) const override { From 485f65dd2631339904e9d9a5e2ac94e132518984 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Fri, 6 Oct 2023 15:48:43 +0200 Subject: [PATCH 3/5] Matplotlib deprecation. (#2223) Remove deprecated styling form LFP example. --- python/example/probe_lfpykit.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/python/example/probe_lfpykit.py b/python/example/probe_lfpykit.py index d7870f6adf..c3cafc11b8 100644 --- a/python/example/probe_lfpykit.py +++ b/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) From 3998b0176c7256b5e28d0a813b2ff95ef3f13f50 Mon Sep 17 00:00:00 2001 From: Simon Frasch Date: Sat, 7 Oct 2023 09:12:37 +0200 Subject: [PATCH 4/5] CI: Fix spack url parsing (#2225) Github seems to have changed their json tags, which identify the url to download the latest release of spack. --- .github/workflows/test-spack.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test-spack.yml b/.github/workflows/test-spack.yml index 5e2b6e3e96..6a5a963ed1 100644 --- a/.github/workflows/test-spack.yml +++ b/.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: | From 2f4c32598d37f9852978c76952b0a09aeb84385b Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Sat, 7 Oct 2023 13:46:49 +0200 Subject: [PATCH 5/5] Allow taking derivatives of `Xd` (#2221) # Introduction The core problem is that `Xd` is not a `STATE` and only those can appear in derivatives/ODEs. This requires rewriting the idiomatic ``` NEURON { SUFFIX decay USEION x WRITE xd RANGE tau } PARAMETER { tau = 5 } BREAKPOINT { SOLVE dX METHOD cnexp } DERIVATIVE dX { xd' = -tau*xd } ``` as this ``` NEURON { SUFFIX decay USEION x WRITE xd RANGE tau } STATE { F } INITIAL { F = xd } PARAMETER { tau = 5 } BREAKPOINT { SOLVE dX METHOD cnexp xd = F } DERIVATIVE dX { F = xd F' = -tau*F } ``` _However_, `BREAKPOINT` is translated into _two_ C++ calls, not one: 1. compute $F(t)$ from $F'(t)$ 2. update `xd` based on `F` and other things -- like the **solver** -- happen in between. This has proven not only ugly, but also unsound, see #2129 and possibly more. # Fixing the Problem Simply allow `Xd` to be differentiated. As this is not an ionic _current_ nothing bad will happen., right? Right? This could still be awkward when making mechanisms like this ``` BREAKPOINT { SOLVE dX METHOD cnexp xd = xd + 42 } DERIVATIVE dX { xd' = -tau*xd } ``` --- doc/fileformat/nmodl.rst | 4 +--- mechanisms/default/decay.mod | 14 ++++---------- modcc/expression.cpp | 27 +++++++++++++++++---------- modcc/printer/cprinter.cpp | 2 +- modcc/solvers.cpp | 11 +++++------ python/example/diffusion.py | 14 ++++++-------- 6 files changed, 34 insertions(+), 38 deletions(-) diff --git a/doc/fileformat/nmodl.rst b/doc/fileformat/nmodl.rst index 10b6aa0fcc..df5bfedda3 100644 --- a/doc/fileformat/nmodl.rst +++ b/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. diff --git a/mechanisms/default/decay.mod b/mechanisms/default/decay.mod index 3dadbb834b..e7f17a1921 100644 --- a/mechanisms/default/decay.mod +++ b/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 } diff --git a/modcc/expression.cpp b/modcc/expression.cpp index add6678f7d..72ac282981 100644 --- a/modcc/expression.cpp +++ b/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_)); + } /******************************************************************************* diff --git a/modcc/printer/cprinter.cpp b/modcc/printer/cprinter.cpp index 88da752d75..36ec914e6d 100644 --- a/modcc/printer/cprinter.cpp +++ b/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; diff --git a/modcc/solvers.cpp b/modcc/solvers.cpp index 8bd3a992b4..13c312d7e7 100644 --- a/modcc/solvers.cpp +++ b/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 { diff --git a/python/example/diffusion.py b/python/example/diffusion.py index 7ef2d2107c..a4ffae6a29 100644 --- a/python/example/diffusion.py +++ b/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")