NMODL SympySolver - sparse

This notebook describes the implementation of the sparse part of the SympySolverVisitor, which solves the systems of ODEs defined in DERIVATIVE blocks when these ODEs are linear and coupled.

For a higher level overview of the approach to solving ODEs in NMODL, please see the nmodl-odes-overview notebook.

For a more general tutorial on using the NMODL python interface, please see the tutorial notebook.

---

Implementation

The SympySolverVisitor for solver method sparse does the following:

• Construct the implicit Euler equations to convert the sysytem of ODEs to a LINEAR block

• This LINEAR block is then solved as described in nmodl-linear-solver

---

Implementation Tests

The unit tests may be helpful to understand what these functions are doing - SympySolverVisitor tests are located in test/visitor/sympy_solver.cpp, and tests using the sparse solver method have the tag “[sparse]”

Examples

%%capture
! pip install nmodl

import nmodl.dsl as nmodl


def run_sympy_solver(mod_string, cse=False):
    # parse NMDOL file (supplied as a string) into AST
    driver = nmodl.NmodlDriver()
    AST = driver.parse_string(mod_string)
    # run SymtabVisitor to generate Symbol Table
    nmodl.symtab.SymtabVisitor().visit_program(AST)
    # constant folding, inlining & local variable renaming passes
    nmodl.visitor.ConstantFolderVisitor().visit_program(AST)
    nmodl.visitor.InlineVisitor().visit_program(AST)
    nmodl.visitor.LocalVarRenameVisitor().visit_program(AST)
    # run SympySolver visitor
    nmodl.visitor.SympySolverVisitor().visit_program(AST)
    # return contents of new DERIVATIVE block as list of strings
    return nmodl.to_nmodl(
        nmodl.visitor.AstLookupVisitor().lookup(
            AST, nmodl.ast.AstNodeType.DERIVATIVE_BLOCK
        )[0]
    )

Ex. 1

2 coupled linear ODEs

ex1 = """
BREAKPOINT {
    SOLVE states METHOD sparse
}
STATE {
    mc m
}
DERIVATIVE states {
    mc' = -a*mc + b*m
    m' = a*mc - b*m
}
"""
print(run_sympy_solver(ex1, cse=True))

DERIVATIVE states {
    EIGEN_NEWTON_SOLVE[2]{
        LOCAL old_mc, old_m
    }{
        old_mc = mc
        old_m = m
    }{
        nmodl_eigen_x[0] = mc
        nmodl_eigen_x[1] = m
    }{
        nmodl_eigen_f[0] = -nmodl_eigen_x[0]*a*dt-nmodl_eigen_x[0]+nmodl_eigen_x[1]*b*dt+old_mc
        nmodl_eigen_j[0] = -a*dt-1.0
        nmodl_eigen_j[2] = b*dt
        nmodl_eigen_f[1] = nmodl_eigen_x[0]*a*dt-nmodl_eigen_x[1]*b*dt-nmodl_eigen_x[1]+old_m
        nmodl_eigen_j[1] = a*dt
        nmodl_eigen_j[3] = -b*dt-1.0
    }{
        mc = nmodl_eigen_x[0]
        m = nmodl_eigen_x[1]
    }{
    }
}

Ex. 2

5 coupled linear ODEs

ex2 = """
STATE {
    c1 o1 o2 p0 p1
}
BREAKPOINT  {
    SOLVE ihkin METHOD sparse
}
DERIVATIVE ihkin {
    LOCAL alpha, beta, k3p, k4, k1ca, k2
    evaluate_fct(v, cai)
    c1' = (-1*(alpha*c1-beta*o1))
    o1' = (1*(alpha*c1-beta*o1))+(-1*(k3p*o1-k4*o2))
    o2' = (1*(k3p*o1-k4*o2))
    p0' = (-1*(k1ca*p0-k2*p1))
    p1' = (1*(k1ca*p0-k2*p1))
}
"""
print(run_sympy_solver(ex2))

DERIVATIVE ihkin {
    EIGEN_NEWTON_SOLVE[5]{
        LOCAL alpha, beta, k3p, k4, k1ca, k2, old_c1, old_o1, old_o2, old_p0, old_p1
    }{
        evaluate_fct(v, cai)
        old_c1 = c1
        old_o1 = o1
        old_o2 = o2
        old_p0 = p0
        old_p1 = p1
    }{
        nmodl_eigen_x[0] = c1
        nmodl_eigen_x[1] = o1
        nmodl_eigen_x[2] = o2
        nmodl_eigen_x[3] = p0
        nmodl_eigen_x[4] = p1
    }{
        nmodl_eigen_f[0] = -nmodl_eigen_x[0]*alpha*dt-nmodl_eigen_x[0]+nmodl_eigen_x[1]*beta*dt+old_c1
        nmodl_eigen_j[0] = -alpha*dt-1.0
        nmodl_eigen_j[5] = beta*dt
        nmodl_eigen_j[10] = 0
        nmodl_eigen_j[15] = 0
        nmodl_eigen_j[20] = 0
        nmodl_eigen_f[1] = nmodl_eigen_x[0]*alpha*dt-nmodl_eigen_x[1]*beta*dt-nmodl_eigen_x[1]*dt*k3p-nmodl_eigen_x[1]+nmodl_eigen_x[2]*dt*k4+old_o1
        nmodl_eigen_j[1] = alpha*dt
        nmodl_eigen_j[6] = -beta*dt-dt*k3p-1.0
        nmodl_eigen_j[11] = dt*k4
        nmodl_eigen_j[16] = 0
        nmodl_eigen_j[21] = 0
        nmodl_eigen_f[2] = nmodl_eigen_x[1]*dt*k3p-nmodl_eigen_x[2]*dt*k4-nmodl_eigen_x[2]+old_o2
        nmodl_eigen_j[2] = 0
        nmodl_eigen_j[7] = dt*k3p
        nmodl_eigen_j[12] = -dt*k4-1.0
        nmodl_eigen_j[17] = 0
        nmodl_eigen_j[22] = 0
        nmodl_eigen_f[3] = -nmodl_eigen_x[3]*dt*k1ca-nmodl_eigen_x[3]+nmodl_eigen_x[4]*dt*k2+old_p0
        nmodl_eigen_j[3] = 0
        nmodl_eigen_j[8] = 0
        nmodl_eigen_j[13] = 0
        nmodl_eigen_j[18] = -dt*k1ca-1.0
        nmodl_eigen_j[23] = dt*k2
        nmodl_eigen_f[4] = nmodl_eigen_x[3]*dt*k1ca-nmodl_eigen_x[4]*dt*k2-nmodl_eigen_x[4]+old_p1
        nmodl_eigen_j[4] = 0
        nmodl_eigen_j[9] = 0
        nmodl_eigen_j[14] = 0
        nmodl_eigen_j[19] = dt*k1ca
        nmodl_eigen_j[24] = -dt*k2-1.0
    }{
        c1 = nmodl_eigen_x[0]
        o1 = nmodl_eigen_x[1]
        o2 = nmodl_eigen_x[2]
        p0 = nmodl_eigen_x[3]
        p1 = nmodl_eigen_x[4]
    }{
    }
}