KINETIC block with function call result into sympy exception: SympySolverVisitor :: solve_non_lin_system python exception: 'Symbol' object is not callable

[pramodk]

Consider a MOD file glia__dbbs_mod_collection__Na__granule_cell.mod which has:

FUNCTION alfa(v(mV))(/ms){ 
	alfa = Q10*Aalfa*exp(v/Valfa) 
}

FUNCTION beta(v(mV))(/ms){ 
	beta = Q10*Abeta*exp(-v/Vbeta) 
}

FUNCTION teta(v(mV))(/ms){ 
	teta = Q10*Ateta*exp(-v/Vteta) 
}
 

KINETIC kstates {
	: 1 riga
	~ C1 <-> C2 (n1*alfa(v),n4*beta(v))
	~ C2 <-> C3 (n2*alfa(v),n3*beta(v))
	~ C3 <-> C4 (n3*alfa(v),n2*beta(v))
	~ C4 <-> C5 (n4*alfa(v),n1*beta(v))
	~ C5 <-> O  (gamma,delta)
	~  O <-> OB (epsilon,teta(v))
	
	: 2 riga
	~ I1 <-> I2	(n1*alfa(v)*a,n4*beta(v)*b)
	~ I2 <-> I3	(n2*alfa(v)*a,n3*beta(v)*b)
	~ I3 <-> I4	(n3*alfa(v)*a,n2*beta(v)*b)
	~ I4 <-> I5 (n4*alfa(v)*a,n1*beta(v)*b)
	~ I5 <-> I6 (gamma,delta)
	
	: connette 1 riga con 2 riga
	~ C1 <-> I1 (Con,Coff)
	~ C2 <-> I2 (Con*a,Coff*b)
	~ C3 <-> I3 (Con*a^2,Coff*b^2)
	~ C4 <-> I4 (Con*a^3,Coff*b^3)
	~ C5 <-> I5 (Con*a^4,Coff*b^4)
	~  O <-> I6 (Oon,Ooff)
	
	CONSERVE C1+C2+C3+C4+C5+O+OB+I1+I2+I3+I4+I5+I6=1
}

when we run this via NMODl we get:

./bin/nmodl mod2c/glia__dbbs_mod_collection__Na__granule_cell.mod host --c passes --inline
[NMODL] [info] :: Processing ../../../build_nmodl/x86_64/corenrn/mod2c/glia__dbbs_mod_collection__Na__granule_cell.mod
[NMODL] [info] :: Running semantic analysis visitor
[NMODL] [info] :: Running symtab visitor
[NMODL] [info] :: Running CVode to cnexp visitor
[NMODL] [info] :: Running code compatibility checker
[NMODL] [info] :: Running verbatim rename visitor
[NMODL] [info] :: Running KINETIC block visitor
[NMODL] [info] :: Running STEADYSTATE visitor
[NMODL] [info] :: Parsing Units
[NMODL] [info] :: Running nmodl inline visitor
[NMODL] [info] :: Running local variable rename visitor
[NMODL] [info] :: Automatically enable sympy_analytic because it exists solver of type sparse
[NMODL] [info] :: Running sympy solve visitor
[NMODL] [warning] :: SympySolverVisitor :: solve_non_lin_system python exception: 'Symbol' object is not callable
[NMODL] [info] :: Running cnexp visitor
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported
[NMODL] [info] :: Running C backend code generator
[NMODL] [warning] :: RenameVisitor :: Renaming variable v at 121.15 to arg_v
[NMODL] [warning] :: RenameVisitor :: Renaming variable v at 122.23 to arg_v
[NMODL] [warning] :: RenameVisitor :: Renaming variable v at 125.15 to arg_v
[NMODL] [warning] :: RenameVisitor :: Renaming variable v at 126.24 to arg_v
[NMODL] [warning] :: RenameVisitor :: Renaming variable v at 129.15 to arg_v
[NMODL] [warning] :: RenameVisitor :: Renaming variable v at 130.24 to arg_v
[NMODL] [warning] :: RenameVisitor :: Renaming variable alfa at 122.2-5 to ret_alfa
[NMODL] [warning] :: RenameVisitor :: Renaming variable beta at 126.2-5 to ret_beta
[NMODL] [warning] :: RenameVisitor :: Renaming variable teta at 130.2-5 to ret_teta
libc++abi: terminating with uncaught exception of type std::runtime_error: PRIME encountered during code generation, ODEs not solved?

Here some quick comments:

• note the function calls alpha, beta and teta. I think they are not inlined in the above example and probably make sense
• but, if we pass the statements with function calls to sympy, it won't be able to do right thing in the absence of function definition?

The above needs to be verified. Calling sympy / solver expert @cattabiani !

[cattabiani]

ehehe, I need to remember. I will be back to you trying this.

[cattabiani]

probably you can just add the functions like in this example:

nmodl/test/unit/visitor/kinetic_block.cpp

Lines 487 to 493 in af7c772

	KINETIC states {
	c0 = f_flux
	~ x <-> y (a, b)
	c1 = f_flux + b_flux
	~ y <-> z (c, d)
	c2 = f_flux - 2*b_flux
	})";

or use the CONSERVE statement. I forgot what exactly it does. Discard this last bit if it does not make sense.

[pramodk]

After discussion with @cattabiani, I will put a simple example showing what we want to do here:

 NEURON {
     SUFFIX glia
     RANGE alfa, beta
     RANGE Aalfa, Valfa, Abeta, Vbeta
 }

 PARAMETER {
     Aalfa = 353.91 ( /ms)
     Valfa = 13.99 ( /mV)
     Abeta = 1.272  ( /ms)
     Vbeta = 13.99 ( /mV)
     n1 = 5.422
     n4 = 0.738
 }

 STATE {
     C1
     C2
 }

 BREAKPOINT {
     SOLVE kstates METHOD sparse
 }

 FUNCTION alfa(v(mV)) {
     alfa = Q10*Aalfa*exp(v/Valfa)
 }

 FUNCTION beta(v(mV)) {
     beta = Q10*Abeta*exp(-v/Vbeta)
 }

 KINETIC kstates {
     ~ C1 <-> C2 (n1*alfa(v),n4*beta(v))
 }

If we run above example via nmodl today then we get:

./bin/nmodl test.mod
[NMODL] [info] :: Processing test.mod
[NMODL] [info] :: Running symtab visitor
[NMODL] [info] :: Running semantic analysis visitor
[NMODL] [info] :: Running CVode to cnexp visitor
[NMODL] [info] :: Running code compatibility checker
[NMODL] [info] :: Running verbatim rename visitor
[NMODL] [info] :: Running KINETIC block visitor
[NMODL] [info] :: Running STEADYSTATE visitor
[NMODL] [info] :: Parsing Units
[NMODL] [info] :: Running local variable rename visitor
[NMODL] [info] :: Automatically enable sympy_analytic because it exists solver of type sparse
[NMODL] [info] :: Running sympy solve visitor
[NMODL] [warning] :: Python bindings are not available
[NMODL] [warning] :: SympySolverVisitor :: solve_non_lin_system python exception: 'Symbol' object is not callable
....

Note that the problematic part is that the python exception with 'Symbol' object is not callable. When we pass a line containing function call beta(v) to sympy:

 ~ C1 <-> C2 (n1*alfa(v),n4*beta(v))

then sympy is not aware of function definition and hence error.

So what we want to do in this case is avoid function calls by transforming KINETIC block like below:

 KINETIC kstates {
     LOCAL x, y
     x = alfa(v)
     y = beta(v)
     ~ C1 <-> C2 (n1*x,n4*y)
 }

i.e. we avoid function calls from all reaction_statement.

With this, I was wondering why we can't just use regular Inlining pass here i.e. using passes --inline should have already worked here!! When I looked into main.cpp then I realised that we run InlineVisitor after KineticBlockVisitor. I don't know / remember if there is any reason for this. I think we can simply run inlining pass bit before.

So I tried:

diff --git a/src/main.cpp b/src/main.cpp
index 1e2a4387..9a388eb9 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -434,6 +434,12 @@ int main(int argc, const char* argv[]) {
             SymtabVisitor(update_symtab).visit_program(*ast);
         }

+        if (nmodl_inline) {
+            logger->info("Running nmodl inline visitor");
+            InlineVisitor().visit_program(*ast);
+            ast_to_nmodl(*ast, filepath("inline"));
+        }
+
         /// note that we can not symtab visitor in update mode as we
         /// replace kinetic block with derivative block of same name
         /// in global scope
@@ -470,12 +476,6 @@ int main(int argc, const char* argv[]) {
         /// that old symbols (e.g. prime variables) are not lost
         update_symtab = true;

-        if (nmodl_inline) {
-            logger->info("Running nmodl inline visitor");
-            InlineVisitor().visit_program(*ast);
-            ast_to_nmodl(*ast, filepath("inline"));
-        }

With this, I was hoping ./bin/nmodl test.mod passes --inline to inline the alfa and beta function calls but it doesn't! 🤔
Using passes --show-symtab I realise that alpha, beta are also declared as RANGE variables! Our current logic in InlineVisitor::visit_function_call tries to access function using symbol table. As RANGE variable is declared fist, the node we are getting here will be RANGE alfa than FUNCTION alfa. If I comment out RANGE alfa, beta then everything seems to work as expected.

So action items could be:

• I am not sure if we should throw an error if someone is declaring RANGE variable as FUNCTION/PROCEDURE as well. But warning make sense!
• In inlining visitor, rather than relying on symbol table to get corresponding node AST node for function/procedure, we can do:
  * use has_any_property() to see if the symbol has function or procedure type
  * if the above is true then use collect_nodes() or something similar to find AST node of type PROCEDURE/FUNCTION of the given name.

@alkino : I have assigned this to you in case you want to tackle this. I or @iomaganaris can also help to clarify details.

[pramodk]

@cattabiani : I think we were a bit wrong i.e. we forgot that that in sympy we already pass the list of user-defined functions here

nmodl/nmodl/ode.py

Line 199 in ea1acbd

function_calls: set of function calls used in the ODE

and unit tests are already here

nmodl/test/unit/visitor/sympy_solver.cpp

Line 389 in ea1acbd

GIVEN("Derivative block with calling external functions passes sympy") {

[alkino]

So here is a lot of problems:

• Name of a range variable can conflict with the name of a function: Store all nodes for a specific symbol #934
• When not inlining usage of locally defined functions can lead to crash: Apply user_functions to vecFcode too #935
• Some locally defined functions have protected names Apply user_functions to vecFcode too #935

[alkino]

There is three tickets on this. See my last comment