Add "matexp" solver to nmodl

[ctrl-z-9000-times]

Hello,

I would like to add a new solver the new nmodl translator. The "matexp" solver computes the analytic solution for kinetic models that are linear, which includes all Markov models of ion channels. This solver is perfectly accurate, but runs between 25-50% slower than the "sparse" solver.

It uses Eigen's implementation of the matrix expopnential function, which is in the unsupported section of the library.
https://libeigen.gitlab.io/eigen/docs-nightly/unsupported/group__MatrixFunctions__Module.html#matrixbase_exp

[codecov]

Codecov Report

❌ Patch coverage is 98.43750% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 68.30%. Comparing base (06084a5) to head (e98f2dc).

Files with missing lines	Patch %	Lines
test/nmodl/transpiler/unit/visitor/matexp.cpp	98.41%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #3574      +/-   ##
==========================================
- Coverage   68.30%   68.30%   -0.01%     
==========================================
  Files         689      691       +2     
  Lines      111033   111097      +64     
==========================================
+ Hits        75844    75881      +37     
- Misses      35189    35216      +27     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[ramcdougal]

Interesting contribution, although I don't know the math.

Have you thought about tests?

Generally our tests are in test/

I do have some concerns about using unsupported library functions.

[ctrl-z-9000-times]

The math is deceptively simple:

Given an ODE: x' = A * x
The solution is: x = e^A * x

However A is a matrix and the exponential function e^A is defined over matrices.

---

I was planning on adding tests and documentation, but only after running this by all of you for review.

[nrnhines]

I was planning on adding tests and documentation, but only after running this by all of you for review.

I believe this is a useful addition. I notice the implementation is for the new nmodl and that is fine. nocmodl will be phased out. A possible place for the test is nrn/test/hoctests where you put the mod file there and the python file in tests.
Each file in tests is run separately instead of with pytest. Looks like there is a place in the documentation for it analogous to https://nrn.readthedocs.io/en/latest/nmodl/transpiler/notebooks/nmodl-sympy-solver-cnexp.html#NMODL-SympySolver---cnexp

[ctrl-z-9000-times]

Thank you both for reviewing this!

---

Todo list:

• Investigate second order correctness
• Unit tests
• Documentation

[JCGoran]

What's the advantage of adding this to the codebase? I can see the utility provided that the sparse solver is not adequate (accuracy-wise) in certain scenarios, but otherwise, at least at face value, it doesn't look like we gain much from it.
As you mentioned, it's slower than the sparse one, and adds extra code to maintain, so I'm trying to understand why it would be worth adding it.

[sonarqubecloud]

Quality Gate passed

Issues
12 New issues
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
0.0% Duplication on New Code

See analysis details on SonarQube Cloud

[ctrl-z-9000-times]

Here is my justification: implementing the analytic solution is useful for understanding the accuracy of other solver methods. The sparse solver method is an approximation; its error is proportional to the time step, dt, and the time step affects how fast the simulation runs. With all numeric approximations, there is a trade-off between run-speed and accuracy. The exact solution will help users navigate this tradeoff.

[nrnhines]

" - exact analytic integration if matrix is constant over the dt interval of integration\n"

Note that rates are often functions of voltage. In these cases, the default fixed step method results in an integration accuracy that is first order in dt. With the fixed step method and secondorder=2 integration accuracy is second order in dt. CVode (variable step methods) ignores the METHOD and uses its own implementations for y' = f(y) and evaluation of an approximate jacobian.

[nrnhines]

Why is there a braces around scalar initializer [-Wbraced-scalar-init] warning and how can we get rid of it?

[ctrl-z-9000-times]

I can't reproduce this warning. Maybe adding a comma between the AstNodeType and the closing brace will convince the compiler that this is an initializer list, not a scalar?

[ctrl-z-9000-times]

As we discussed in today's meeting, I will add a test that compares the matexp and cnexp solvers using an HH model.

[JCGoran]

Implementation-wise, I believe this can be done better. With the verbatim blocks everywhere, this approach mixes the implementation (C++ code) with the underlying abstraction (NMOD language), which is not how any of the other visitors in NMODL work.

The usual approach to adding a new feature to NMODL is to check whether the AST needs to be manipulated, and if so, in what way. Then, implement a visitor which does the manipulation of NMODL constructs. After the AST has the correct structure, you can proceed and implement things at the codegen stage, for which no new visitors are necessary (unless you want to generate code for another programming language).

For instance, the current implementation directly inserts a verbatim block with Eigen headers in NMODL; this is completely unnecessary at this stage, and should be done at the codegen step instead. Depending on whether the changes work for coreNEURON or NEURON, one can insert it in codegen_coreneuron_cpp_visitor, codegen_neuron_cpp_visitor, or simply codegen_cpp_visitor if the changes work for both backends.

You can have a look at the following 2 PRs for reference how a new feature can be implemented with the 2-stage approach in NMODL (just an example, plenty of others in that repo):

• BlueBrain/nmodl#1467 (AST manipulation)
• BlueBrain/nmodl#1493 (actual codegen)

Note that opening 2 PRs is not strictly necessary, and they can be combined into one, but it was done this way (at least for this particular feature) to keep the changeset to a manageable size and not overload the reviewer.

[ctrl-z-9000-times]

Hi Everyone,

This PR is ready for another round of review.

Goran: I replaced all of the verbatim with a new AST node and some modifications to the codegen.

Hines: I added a testcase to verify that for ODEs with a single variable the CNEXP solver yields identical results. This test compares two formulations of the same Hodgkin Huxley model:

• HH model with three states and a derivative block, solved by cnexp
• HH model with six states and a kinetic block, solved by matexp

Thanks

[ctrl-z-9000-times]

@JCGoran, I took your advice and reimplemented this PR. I'd really appreciate if you could review this again before you leave?

[JCGoran]

@ctrl-z-9000-times Do you have some docs on how the implementation works conceptually? I see there is a new node type for the codegen, but I'm having a bit of trouble following what's going on due to the size of the PR. Also, shouldn't there be an entry in the SOLVE docs about this new addition?

Some general comments:

• anything not declared static in an implementation file should have a docstring (including the visit_* methods)
• I see regexes, and wonder if they are needed. Usually to compare strings at the NMODL level you can use reindent_text
• at first glance, the mod file tests do not really look like "minimum" working examples
• are the "broken AST" exceptions necessary? If you are doing the replacement correctly, the parser should catch any issues with malformed ASTs
• MatexpVisitor::visit_reaction_statement (maybe others as well) has a lot of variables, maybe it's best to break it down into some internal helper functions

Most of the PR looks good, and I'd be happy to have a final look once the above is answered.

[ctrl-z-9000-times]

@JCGoran thank you for the detailed feedback

The documentation is at:
docs/nmodl/transpiler/notebooks/nmodl-matexp-solver.ipynb

1. Ok, I added static declarations and more documentation strings.
2. Ok, I replaced the regexes with reindent_text.
3. Yes, the mod files are necessary. The two Hodgkin-Huxley models are for verifying the functional output of the solver.
4. Yes, the "broken AST" errors are necessary. Some are null pointer checks and others are in places that should be unreachable.
   EDIT: I replaced the exceptions with standard assert statements.
5. Ok, I broke up the method visit_reaction_statement, splitting out two additional methods transform_decay_statement and transform_reaction_statement

[sonarqubecloud]

Quality Gate passed

Issues
14 New issues
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
17.1% Duplication on New Code

See analysis details on SonarQube Cloud

[sonarqubecloud]

Quality Gate passed

Issues
14 New issues
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
17.1% Duplication on New Code

See analysis details on SonarQube Cloud