nocmodl: diffeq.y:194: difeq_yyerror: Assertion `0' failed.

[davidsterratt]

I'm getting the following (cryptic) error when compiling the following nmodl file:
~/datastore/nrn-5.8/i686/bin/nrnivmodl
"/home/sterratt/datastore/nrn-5.8/i686/bin/nocmodl" psd
Translating psd.mod into psd.c
nocmodl: diffeq.y:194: difeq_yyerror: Assertion `0' failed.
make: *** [psd.lo] Aborted

It seems to be a problem with the DERIVATIVE block (no error when I delete it), but I can't get any further than that. (I've commented out a fair amount to simplify things). Any help would be much appreciated.

David.

psd.mod:

Code: Select all

COMMENT

psd.mod "Postsynaptic density" based on Rubin et al (2005) "Calcium
Time Course as a signal for STDSP".

ENDCOMMENT

NEURON {
	  POINT_PROCESS psd
	  : RANGE  P, V, A, B, D, W
    RANGE  V
    USEION ca READ cai
    GLOBAL pHC, pHN, aHC, aHN, thetav, sigmav, thetad, sigmad, thetab, sigmab, taup, kp, taua, kd, tauv, alphav, taud, alphad, taub, alphab, cp, cd, tauw, p, alphaw, d, betaw
}

UNITS {
    (molar) = (/liter)
    (mM) = (millimolar)
    (uM) = (micromolar)
}

PARAMETER {
    pHC     = 4   (uM)
    pHN     = 4
    aHC     = 0.6 (uM)
    aHN     = 3
    thetav  = 2 (uM)
    sigmav  = -0.05
    thetad  = 2.6 
    sigmad  = -0.01
    thetab  = 0.55
    sigmab  = -0.02
    taup    = 500 (ms)
    kp      = -0.1
    taua    = 5 (ms)
    kd      = -0.002
    tauv    = 10 (ms)
    alphav  = 1.0
    taud    = 250 (ms)
    alphad  = 1.0
    taub    = 40 (ms)
    alphab  = 5.0
    cp      = 5
    cd      = 4
    tauw    = 500 (ms) 
    p       = 0.3
    alphaw  = 0.8
    d       = 0.01
    betaw   = 0.6
}

ASSIGNED {
    cai (mM)
}

: STATE {  P V A B D W }
STATE {  V }

BREAKPOINT {
	  SOLVE state METHOD cnexp
}

DERIVATIVE state {
:    P' = (10     * hill(   cai,pHC,   pHN)    - cp*A*P    )/taup
     V' = (alphav * sigm(cai,thetav,sigmav) - V)/tauv
:    A' = (         hill(   cai,aHC,   aHN)    - A         )/taua
:    B' = (alphab * sigmoid(A,  thetab,sigmab) - B - cd*B*V)/taub
:    D' = (alphad * sigmoid(D,  thetad,sigmad) - D         )/taud
:    W' = (alphaw * sigmoid(P,  p,     kp)     - betaw * sigmoid(D, d, kd) - W)/tauw
}

INITIAL {
:    P = 0
    V = 0 
:    A = 0
:    B = 0 
:    D = 0
:    W = 1
}

FUNCTION hill(x (mM), HC (uM), HN) {
    hill = 1 / ( (x/((0.001)*HC))^HN + 1 )
}

FUNCTION sigm(x, theta, sigma) {
    sigm = 1 / ( 1 + exp((x-theta)/sigma) )
}

[Raj]

Hi David,

This is quite bizar. I stripped it further down but still no hint. The documentation on NMODL however shows many examples in which calling a function is avoided and a procedure setting a variable is called instead. I couldn't find a reason why this was done. The following, however, compiles and the approach in it might work for you:

Code: Select all

NEURON {
	POINT_PROCESS psdproc
}

ASSIGNED{
	myff (1)
}

STATE {
	m
}

BREAKPOINT {
	SOLVE state METHOD cnexp
}

DERIVATIVE state {
	myfunkyfunc()
	m'=myff
}

PROCEDURE myfunkyfunc() {
	myff=1
}

and the minimal variant showing your problem:

Code: Select all

NEURON {
POINT_PROCESS psdminimal
}

STATE {
m
}

BREAKPOINT {
SOLVE state METHOD cnexp
}

DERIVATIVE state {
m'=myfunkyfunc()
}

FUNCTION myfunkyfunc() {
myfunkyfunc=1
}

It would be nice to get some understanding of this behaviour, the comments in the diffeq.* files about the 0.0 treatment which should be changed when version 2 of a certain something is reached might be an important hint for those who can read yacc.

Good luck!

[ted]

Excellent question, David, and thoughtful dissection of the intermediate files, Raj.
What happened here is a result of the fact that "cnexp is limited to linear equations
that are very easy to parse" (to quote an old email from Michael Hines). After all, it
was intended for use with HH-style mechanisms, which involve very simple ODEs.

The workaround is basically a trick: use LOCAL variables to keep the DEs simple.
Changing the DERIVATIVE block to this

Code: Select all

DERIVATIVE state {
  LOCAL xV
:    P' = (10     * hill(   cai,pHC,   pHN)    - cp*A*P    )/taup

:     V' = (alphav * sigm(cai,thetav,sigmav) - V)/tauv
     xV = (alphav * sigm(cai,thetav,sigmav) - V)/tauv
     V' = xV

:    A' = (         hill(   cai,aHC,   aHN)    - A         )/taua
:    B' = (alphab * sigmoid(A,  thetab,sigmab) - B - cd*B*V)/taub
:    D' = (alphad * sigmoid(D,  thetad,sigmad) - D         )/taud
:    W' = (alphaw * sigmoid(P,  p,     kp)     - betaw * sigmoid(D, d, kd) - W)/tauw
}

lets the code compile without complaint.

One more caveat about integration methods--instead of cnexp, you must use
SOLVE states METHOD derivimplicit
if any of the following applies:
--a DE is nonlinear
--a DE's RHS depends on any states
--two or more DEs are coupled

--Ted

[davidsterratt]

Excellent question, David, and thoughtful dissection of the intermediate files, Raj.
What happened here is a result of the fact that "cnexp is limited to linear equations
that are very easy to parse" (to quote an old email from Michael Hines). After all, it
was intended for use with HH-style mechanisms, which involve very simple ODEs.

Thanks very much for the speedy response Raj and Ted.

The workaround is basically a trick: use LOCAL variables to keep the DEs simple.
Changing the DERIVATIVE block to this

Code: Select all

DERIVATIVE state {
  LOCAL xV
:    P' = (10     * hill(   cai,pHC,   pHN)    - cp*A*P    )/taup

:     V' = (alphav * sigm(cai,thetav,sigmav) - V)/tauv
     xV = (alphav * sigm(cai,thetav,sigmav) - V)/tauv
     V' = xV

:    A' = (         hill(   cai,aHC,   aHN)    - A         )/taua
:    B' = (alphab * sigmoid(A,  thetab,sigmab) - B - cd*B*V)/taub
:    D' = (alphad * sigmoid(D,  thetad,sigmad) - D         )/taud
:    W' = (alphaw * sigmoid(P,  p,     kp)     - betaw * sigmoid(D, d, kd) - W)/tauw
}

lets the code compile without complaint.

Yes this seems to work fine. I suppose that this trick is numerically equivalent to the one that Raj mentioned? I ask as there are quite a few mod files that set the RHSs of the DEs in a PROCEDURE block.

One more caveat about integration methods--instead of cnexp, you must use
SOLVE states METHOD derivimplicit
if any of the following applies:
--a DE is nonlinear
--a DE's RHS depends on any states
--two or more DEs are coupled

I tried just changing the integration method to derivimplicit (i.e. not using local variables or procedures as above). However, this didn't affect the parsing, so it looks like the strict parsing of the RHS doesn't depend on which integration method is specified. If derivimplicit really can cope with anything on the RHS, would it be possible to alter nrnivmodl so that it only does the strict parsing if cnexp (or other restrictive integration methods) are specified, and if the parse fails, gives an error message suggesting the use of derivimplicit?

David.

[hines]

Sorry. The intention is to handle that list of equations and automatically
switch from the requested cnexp method to the derivimplicit method.
My error was that I programmed the parser for single argument
functions and so it fails on multiargument functions. If you remove
the extra arguments from each function (which is nonsense but allows
parsing) you would get the result

Code: Select all

nocmodl temp1
Translating temp1.mod into temp1.c
Could not translate using cnexp method; using derivimplicit

The cnexp method can only be used with equations of the form
y' = f(y) in whch f(y) is linear in y and there are no other states involved
in f(y). So it handles

Code: Select all

A' = (         hill(   cai)    - A         )/taua

just fine and automatically switches to derivimplicit when it sees

Code: Select all

 D' = (alphad * sigmoid(D) - D         )/taud

Unfortunately it fails to realize that

Code: Select all

P' = (10     * hill(   cai)    - cp*A*P    )/taup

involves another state so fails to switch. So if only the A' and P'
equations are present the switch from cnexp to derivimplicit fails to
occur and the simulation would not be as numerically stable as it
would be with derivimplicit. I will fix both the multiargument
function problem and this multiple state in a linear equation problem
after I return from Brazil on Sunday.

To perhaps confuse matters a bit, I was surprised that an explicit change to derivimplicit did not fix the problem because there is no reason in that case to use the cnexp specific parser. Then I realized that the parser is also used to determine the diagonal jacobian elements for the cvode variable step method. So the parser is being used for ALL derivitive blocks and cannot be avoided. This is for efficiency reasons
in which if the equation is nonlinear, y' = f(y) one need to compute
f(y+dy) - f(y), but if the equation is linear in y, e.g y'=x*y*z
it is still linear in y even though it involves other states and
dy'/dy = x*z
which is much faster to compute than (x*(y+.001)*z - x*y*z)/.001

Anyway, DO NOT replace multiargument functions involving the state
mentioned on the left hand side with local variables. That will ruin the numerical stability properties of the translated equations. i.e

Code: Select all

D' = (alphad * sigmoid(D,  thetad,sigmad) - D         )/taud

cannot be written as

Code: Select all

xxx = alphad*sigmoid(D, thetad, sigmad)
D' = (xxx - D)/taud

because xxx involves D and the diagonal jacobian element needs that
contribution.

[davidsterratt]

Thanks in advance for the fix and thanks for the explanation Mike -- I now understand a bit more about the numerical solution methods in NEURON. Just to check I've understood it right, when the fix has been applied, the diagaonal element of the Jacobian corresponding to
<code>
V' = (alphav * sigm(cai,thetav,sigmav) - V)/tauv
</code>
will be computed analitically as -1/tauv? If this is the case, I suppose the derivatives of V wrt to the other states could be very efficiently computed too, since they are 0. Does NEURON do this?

David.

[hines]

The multiargument bug and a few related problems (no arguments or one argument which is 0) is fixed in the 5.8.117 versions at http://www.neuron.yale.edu/ftp/neuron/versions/alpha/
See the 498 change set log message at
http://www.neuron.yale.edu/viewcvs-bin/ ... iew=branch
Note that this does not fix the "not switching from cnexp to derivimplicit when other states besides the diagonal appears" problem. I think I will let that alone for now because I am not entirely sure it is a bug.
You are correct in supposing that the diagonal jacobian element is -1/tauv and that is what the parser decides. You are also correct in imagining that all the off diagonal jacobian elements could be computed in the same way; but, although the parser can do that also, it is being used in the cvode case to only do the diagonal elements. If you want to use the complete analytic jacobian then it would be best to use the kinetic scheme representation, which is usually the precursor rep for sets of diffeqs. I should mention that when a set of diffeqs uses the derivimplicit method, the full (not sparse) jacobian matrix is numerically computed. This isn't too bad for small sets of equations but again there are good efficiency reasons to use either cnexp or kinetic schemes.