BlueBrain / nmodl

Code Generation Framework For NEURON MODeling Language
https://bluebrain.github.io/nmodl/
Apache License 2.0
53 stars 15 forks source link

sympy returns string "// Not supported in C" when it doesn't know about a variable #156

Open ohm314 opened 5 years ago

ohm314 commented 5 years ago

When sympy does not know about a variable used in an expression and we call sympy.ccode then it will return a multiline string with an error message: example:

// Not supported in C:
// sawtooth
dt*sawtooth(m) + m

Currently in

https://github.com/BlueBrain/nmodl/blob/aeaf72cfcb1e3bb40010f19b6d394b0ce5978b4e/nmodl/ode.py#L374

this is not checked and properly handled.

This is also related to #118

lkeegan commented 5 years ago

@ohm314 this is a good point!

do you have any mod files with examples of this happening?

maybe one solution would be to wrap the ccode call with a function that raises a python exception if the first line is // Not supported in C:? then in SympySolver for derivimplicit if we get a python exception we use the fall-back solver instead?

pramodk commented 5 years ago

do you have any mod files with examples of this happening?

I believe if you run we run SymPy passes without function inlining then we will often see this problem. For example below one:

TITLE T-calcium channel From Migliore CA3

UNITS {
    (mA) = (milliamp)
    (mV) = (millivolt)

    FARADAY = 96520 (coul)
    R = 8.3134 (joule/degC)
    KTOMV = .0853 (mV/degC)
}

NEURON {
    SUFFIX rd_cat
    USEION tca READ etca WRITE itca VALENCE 2
    USEION ca READ cai, cao VALENCE 2
        RANGE gcatbar,cai, itca, etca
}

PARAMETER {
    gcatbar=.003 (mho/cm2)
    cao (mM)
}

STATE {
    m h
}

ASSIGNED {
    v (mV)
    celsius (degC)
    cai (mM)
    itca (mA/cm2)
        gcat (mho/cm2)
    etca (mV)
}

INITIAL {
      m = minf(v)
      h = hinf(v)
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    gcat = gcatbar*m*m*h
    itca = gcat*ghk(v,cai,cao)

}

DERIVATIVE states { : exact when v held constant
    m' = (minf(v) - m)/m_tau(v)
    h' = (hinf(v) - h)/h_tau(v)
}

FUNCTION ghk(v(mV), ci(mM), co(mM)) (mV) {
        LOCAL nu,f

        f = KTF(celsius)/2
        nu = v/f
        ghk=-f*(1. - (ci/co)*exptrap(1,nu))*efun(nu)
}

FUNCTION KTF(celsius (DegC)) (mV) {
        KTF = ((25./293.15)*(celsius + 273.15))
}

FUNCTION efun(z) {
    if (fabs(z) < 1e-4) {
        efun = 1 - z/2
    }else{
        efun = z/(exptrap(2,z) - 1)
    }
}

FUNCTION hinf(v(mV)) {
    LOCAL a,b
    a = 1.e-6*exptrap(3,-v/16.26)
    b = 1/(exptrap(4,(-v+29.79)/10)+1)
    hinf = a/(a+b)
}

FUNCTION minf(v(mV)) {
    LOCAL a,b
    a = 0.2*(-1.0*v+19.26)/(exptrap(5,(-1.0*v+19.26)/10.0)-1.0)
    b = 0.009*exptrap(6,-v/22.03)
    minf = a/(a+b)
}

FUNCTION m_tau(v(mV)) (ms) {
    LOCAL a,b
:   TABLE FROM -150 TO 150 WITH 200
    a = 0.2*(-1.0*v+19.26)/(exptrap(7,(-1.0*v+19.26)/10.0)-1.0)
    b = 0.009*exptrap(8,-v/22.03)
    m_tau = 1/(a+b)
}

FUNCTION h_tau(v(mV)) (ms) {
    LOCAL a,b
    a = 1.e-6*exptrap(9,-v/16.26)
    b = 1/(exptrap(10,(-v+29.79)/10.)+1.)
    h_tau = 1/(a+b)
}

FUNCTION exptrap(loc,x) {
  if (x>=700.0) {
    exptrap = exp(700.0)
  } else {
    exptrap = exp(x)
  }
}

Gives :

 ./bin/nmodl ~/.../mechanisms/tca.mod  passes sympy --analytic
[NMODL] [info] :: Processing /Users/kumbhar/workarena/repos/bbp/reduced_dentate/mechanisms/tca.mod
...
[NMODL] [info] :: Running sympy solve visitor
libc++abi.dylib: terminating with uncaught exception of type std::runtime_error: NMODL Parser Error : syntax error, unexpected / [Location : 1.25]
Abort trap: 6

then in SympySolver for derivimplicit if we get a python exception we use the fall-back solver instead?

:thumbsup: