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

@@ -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.
 

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_));
+
 }
 
 /*******************************************************************************

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")