Page 1 of 2
STATEs which are arrays
Posted: Sat Apr 05, 2008 4:02 pm
by sec6
I'm attempting to write a mechanism in NMODL, using a STATE which is an array.
Can someone tell me what is wrong with the following code?
Code: Select all
NEURON {
POINT_PROCESS arrTest
RANGE outVar
}
STATE {
outVar[1]
}
ASSIGNED {
vec[1]
}
INITIAL {
outVar[0] = 0
vec[0] = 0
}
BREAKPOINT {
SOLVE states METHOD cnexp
}
DERIVATIVE states {
outVar[0]' = 0
:outVar' = vec
:outVar' = 0
:outVar'[0] = 0
}
When I try to compile the toy program above, I get the following error message:
Code: Select all
syntax error:
Illegal block at line 24 in file arrTest.mod
outVar[0]' = 0
^
make: *** [arrTest.c] Error 1
Three alternative versions of what seems to be the culprit line (commented out, in the DERIVATIVE block in the program above) also won't compile.
For the first and second of these alternatives, the error message is:
Code: Select all
variable needs an index:
Illegal block at line 25 in file arrTest.mod
outVar' = vec
^
make: *** [arrTest.c] Error 1
and
Code: Select all
variable needs an index:
Illegal block at line 26 in file arrTest.mod
outVar' = 0
^
make: *** [arrTest.c] Error 1
For the third alternative, I get a cascade of error messages suggesting that that one isn't even close.
Posted: Sat Apr 05, 2008 6:40 pm
by hines
The syntax is
However I see that the cvode translation is broken for arrays so you cannot use it til I fix the nocmodl translator.
Posted: Sat Apr 05, 2008 7:17 pm
by sec6
However I see that the cvode translation is broken for arrays so you cannot use it til I fix the nocmodl translator.
Drat. Any idea when the nocmodl translater'll be fixed? I need to decide whether to write the program in some less elegant way, vs. waiting for the nocmodl translater to be fixed.
Also, what I really want to do is something like:
Code: Select all
DERIVATIVE {
FOR cnt = 0 TO N {
outVar'[cnt] = fooBar(cnt)
}
}
Except I can't actually put a FOR in a DERIVATIVE block because a DERIVATIVE block is a set of equations to be solved, and not a sequence of instructions to be executed right?
Is there a right way to accomplish what I'm trying to do?
Posted: Sat Apr 05, 2008 8:56 pm
by hines
I see that the general handling of arrays
in DERIVATIVE blocks for cvode, cnexp, and derivimplicit is a lot more difficult than I envisioned and the fix could be a long time coming. In your case, if you don't
mind doing some manual repair to the generated c file it is possible to get it
working. First , avoid the loop index
as an arg to foobar. i.e
Code: Select all
DERIVATIVE states {
LOCAL x
FROM i = 0 TO 2 {
x = i
outVar'[i] = -foobar(x)
}
}
FUNCTION foobar(a) {
foobar = a
}
will emit the c code:
Code: Select all
/*CVODE*/
static int _ode_spec1 () {_reset=0;
{
double _lx ;
{int _li ;for ( _li = 0 ; _li <= 2 ; _li ++ ) {
_lx = ((double) _li ) ;
DoutVar [ _li ] = - foobar ( _lx ) ;
} }
}
return _reset;
}
static int _ode_matsol1() {
double _lx ;
{int _li ;for ( _li = 0 ; _li <= 2 ; _li ++ ) {
_lx = ((double) _li ) ;
DoutVar [ _li ] = DoutVar / (1. - dt*( 0.0 )) ;
}
/*END CVODE*/
static int states () {_reset=0;
{
double _lx ;
{int _li ;for ( _li = 0 ; _li <= 2 ; _li ++ ) {
_lx = ((double) _li ) ;
DoutVar [ _li outVar = outVar - dt*(- ( - foobar ( _lx ) ) ) ;
} }
}
return 0;
}
double foobar ( _la )
double _la ;
{
double _lfoobar;
_lfoobar = _la ;
return _lfoobar;
}
_ode_spec1 is ok as is.
_ode_matsol1 has two problems. First
the second occurrence of DoutVar is missing the [_li] suffix. Second, before the closing } of the function needs to
have three, total. i.e replace the closing
} with }}}. Actually, in the case you gave, the whole body is useless since the computation is DoutVar[_li] = DoutVar_li];
The main line in the states () function up to the -dt token is completely botched and the line should be
Code: Select all
outVar[_li] = outVar[_li] - dt*(.....
If you remove the
from your mod file statements you can see the proper form of the translated c code. The form will be different if one of the args to foobar is the dependent variable of the equation.
By the way, the c code is in the subdirectory (e.g. i686) created by nrnivmodl. You can change the c code
and rerun nrnivmodl and the translator will not overwrite it (unless the mod file
gets changed and has a date later than the c file.)
Posted: Sat Apr 05, 2008 9:15 pm
by hines
I should mention that another way to calculate y'
= f(i) is as a kinetic scheme
(for which state arrays do work). ie.
Code: Select all
PARAMETER {dummy = -1}
STATE {outVar[3]}
BREAKPOINT {
SOLVE states METHOD sparse
}
KINETIC states {
FROM i = 0 TO 2 {
~ outVar[i] <-> dummy (0, foobar(i))
}
}
again, if foobar(i) is a function of outVar(i), then divide the foobar body by outVar(i) and write
~ outVar
<-> dummy (foobar(i), 0)
But note that if outVar really does
appear in the denominator (as in a michaelis menton equation the whole thing is likely to be numerically unstable and one would need a better kinetic scheme with an extra state to simulate enzyme kinetics.
Posted: Sun Apr 06, 2008 8:38 am
by sec6
Hmm. Both suggestions would be a bit difficult, but it's the right kind of difficult, i.e. worth the effort and teaches me something useful.
The kinetic scheme is tempting ('cause I can keep my hands off the C code) but would require rethinking the model at a conceptual level. If I do that, can I
1) Have rate constants which change over time, with discontinuities (the discontinuities are "kosher" i.e. occur in the NET_RECEIVE block, using net_send, and so on).
2) Have a) one-way reactions and b) infinite sources & sinks, e.g. A<=>B, where k_ab >0, k_ba =0, and the concentration of A never changes no matter how long the reaction runs (but the concentration of B does change)?
However, your comment about numerical instability frightens me: I understand very little of such issues, and try to stay far far away from any such problems.
I think I'll take a crack at hand-patching the C code. It may be a little while before I make enough progress to post anything substantive to this thread.
Posted: Sun Apr 06, 2008 8:58 am
by hines
In kinetic schemes, rate constants can
change with time (with discontinuities
in both states and parameters that are
changed within NET_RECEIVE).
Any reactant that is a PARAMETER or ASSIGNED is not a part of the list of dependent variables (states) and is constant so they play the role of infinite sources and sinks. Don't worry about numerical instabilitites if foobar(i) is not a function of states.
Posted: Thu Apr 24, 2008 2:08 pm
by sec6
sec6 wrote:Hmm. Both suggestions would be a bit difficult, but it's the right kind of difficult, i.e. worth the effort and teaches me something useful.
The kinetic scheme is tempting ('cause I can keep my hands off the C code) but would require rethinking the model at a conceptual level. If I do that, can I
1) Have rate constants which change over time, with discontinuities (the discontinuities are "kosher" i.e. occur in the NET_RECEIVE block, using net_send, and so on).
2) Have a) one-way reactions and b) infinite sources & sinks, e.g. A<=>B, where k_ab >0, k_ba =0, and the concentration of A never changes no matter how long the reaction runs (but the concentration of B does change)?
However, your comment about numerical instability frightens me: I understand very little of such issues, and try to stay far far away from any such problems.
I think I'll take a crack at hand-patching the C code. It may be a little while before I make enough progress to post anything substantive to this thread.
Hand-patching the C code proved to be a workable (albeit tedious) solution.
Re: STATEs which are arrays
Posted: Tue Nov 09, 2010 10:53 am
by patoorio
Hi,
I'm trying to the same as sec6. When compiling, after mknrndll stops there is a "mymodfile.c" file left but it is empty. Where am I supposed to edit and fix the c code?
PS: I'm working in Windows. Please don't tell me that fixing the c code is only possible under Linux.
Regards.
Re: STATEs which are arrays
Posted: Wed Nov 10, 2010 7:41 am
by hines
Is it empty if you start an rxvt terminal window, cd to the proper folder, and type
nocmodl mymodfile
In that folder it should translate the mymodfile.mod to mymodfile.c
If it is empty it means that the failure occurrd not as a mistranslation but due to an error in the mod file. If it is not clear what the error is due to the message generated please send the mod file to me at michael dot hines at yale dot edu.
Re: STATEs which are arrays
Posted: Wed Nov 10, 2010 8:13 am
by patoorio
Thanks!
There was indeed an error in the mod file. I fixed it, and now both mknrndll and nocmodl create a non-empty .c file.
I see where and how to fix the c file it but I'm having another problem: unlike nrnivmodl, mknrndll does overwrites the .c file. It doesn't seem to care about file dates, it starts over and the error is generated again.
What commands should I issue in the rxvt terminal window?
Regards.
Re: STATEs which are arrays
Posted: Wed Nov 10, 2010 10:03 am
by hines
I've been looking at the gnu make reference about implicit rules but I have not yet been able to
make the most obvious strategy work, ie. modify c:/nrn72/lib/mknrndll.mak use the implicit rules
%.c : %.mod
nocmodl $*
%.o : %.c
gcc $(CFLAGS) -c $*.c
No matter what I do, it seems to continue to use a built-in implicit rule that uses m2c.
Anyway, if you can work around that problem, then you can also modify
c:/nrn72/lib/mknrndl2.sh
to eliminate the fragment:
for i in $prefixes
do
rm -f $i.c
done
Sorry I cannot be more definitive. I remember working on this some years ago and
finally settling on creating the .o file directly from the .mod file with a rule
that removes the intermediate .c file. I have no idea why cygwin make is causing this problem.
On linux the old style
.mod.c:
"$(bindir)/nocmodl" $*
.c.lo:
...
works well, but I can't get that to work in the cygwin version either.
Re: STATEs which are arrays
Posted: Wed Nov 10, 2010 12:15 pm
by patoorio
Thanks a lot, Mike.
I'll give it a try. I still can use a Linux machine, though ;)
Re: STATEs which are arrays
Posted: Wed Nov 10, 2010 2:19 pm
by hines
A bit more experimenting and reading the gnu make manual allowed resolution of the implicit rule problem. i.e. see
http://www.neuron.yale.edu/hg/neuron/nr ... 4ad065edb0
so that mknrndll no longer has to delete the *.c files derived from the *.mod files.
Since mknrndll.mak.in is used by configure to create mknrndll.mak I created a new setup.exe with this change.
http://www.neuron.yale.edu/ftp/neuron/v ... -setup.exe
Re: STATEs which are arrays
Posted: Wed Nov 10, 2010 3:05 pm
by patoorio
Almost, almost there!!
The DERIVATIVE block is:
Code: Select all
DERIVATIVE states {
mh'[0] = (-3*am-ah)*mh[0] + bm*mh[1] + bh*mh[4]
mh'[1] = (-2*am-bm-ah)*mh[1] + 3*am*mh[0] + 2*bm*mh[2] + bh*mh[5]
mh'[2] = (-am-2*bm-ah)*mh[2] + 2*am*mh[1] + 3*bm*mh[3] + bh*mh[6]
mh'[3] = (-3*bm-ah)*mh[3] + am*mh[2] + bh*mh[7]
mh'[4] = (-3*am-bh)*mh[4] + bm*mh[5] + ah*mh[0]
mh'[5] = (-2*am-bm-bh)*mh[5] + 3*am*mh[4] + 2*bm*mh[6] + ah*mh[1]
mh'[6] = (-am-2*bm-bh)*mh[6] + 2*am*mh[5] + 3*bm*mh[7] + ah*mh[2]
mh'[7] = (-3*bm-bh)*mh[7] + am*mh[6] + ah*mh[3]
}
The relevant lines in the .c file after nrnoc (or mknrndll or nrnivmodl) are
Code: Select all
static int _ode_matsol1 (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) {
Dmh [ 0 ] = Dmh / (1. - dt*( (( - 3.0 * am - ah ))*(1.0) + (bm)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 1 ] = Dmh / (1. - dt*( (( - 2.0 * am - bm - ah ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 2 ] = Dmh / (1. - dt*( (( - am - 2.0 * bm - ah ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 3 ] = Dmh / (1. - dt*( (( - 3.0 * bm - ah ))*(1.0) + (am)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 4 ] = Dmh / (1. - dt*( (( - 3.0 * am - bh ))*(1.0) + (bm)*(1.0) + (ah)*(1.0) )) ;
Dmh [ 5 ] = Dmh / (1. - dt*( (( - 2.0 * am - bm - bh ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (ah)*(1.0) )) ;
Dmh [ 6 ] = Dmh / (1. - dt*( (( - am - 2.0 * bm - bh ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (ah)*(1.0) )) ;
Dmh [ 7 ] = Dmh / (1. - dt*( (( - 3.0 * bm - bh ))*(1.0) + (am)*(1.0) + (ah)*(1.0) )) ;
}
/*END CVODE*/
static int states (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) { {
Dmh [ 0 mh = mh + (1. - exp(dt*((( - 3.0 * am - ah ))*(1.0) + (bm)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(am) - ah ))*(1.0) + (bm)*(1.0) + (bh)*(1.0) ) - mh) ;
Dmh [ 1 mh = mh + (1. - exp(dt*((( - 2.0 * am - bm - ah ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( (- 2.0)*(am) - bm - ah ))*(1.0) + ((3.0)*(am))*(1.0) + ((2.0)*(bm))*(1.0) + (bh)*(1.0) ) - mh) ;
Dmh [ 2 mh = mh + (1. - exp(dt*((( - am - 2.0 * bm - ah ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( - am - (2.0)*(bm) - ah ))*(1.0) + ((2.0)*(am))*(1.0) + ((3.0)*(bm))*(1.0) + (bh)*(1.0) ) - mh) ;
Dmh [ 3 mh = mh + (1. - exp(dt*((( - 3.0 * bm - ah ))*(1.0) + (am)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(bm) - ah ))*(1.0) + (am)*(1.0) + (bh)*(1.0) ) - mh) ;
Dmh [ 4 mh = mh + (1. - exp(dt*((( - 3.0 * am - bh ))*(1.0) + (bm)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(am) - bh ))*(1.0) + (bm)*(1.0) + (ah)*(1.0) ) - mh) ;
Dmh [ 5 mh = mh + (1. - exp(dt*((( - 2.0 * am - bm - bh ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( (- 2.0)*(am) - bm - bh ))*(1.0) + ((3.0)*(am))*(1.0) + ((2.0)*(bm))*(1.0) + (ah)*(1.0) ) - mh) ;
Dmh [ 6 mh = mh + (1. - exp(dt*((( - am - 2.0 * bm - bh ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( - am - (2.0)*(bm) - bh ))*(1.0) + ((2.0)*(am))*(1.0) + ((3.0)*(bm))*(1.0) + (ah)*(1.0) ) - mh) ;
Dmh [ 7 mh = mh + (1. - exp(dt*((( - 3.0 * bm - bh ))*(1.0) + (am)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(bm) - bh ))*(1.0) + (am)*(1.0) + (ah)*(1.0) ) - mh) ;
}
return 0;
}
Which then I change to
Code: Select all
static int _ode_matsol1 (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) {
Dmh [ 0 ] = Dmh [ 0 ] / (1. - dt*( (( - 3.0 * am - ah ))*(1.0) + (bm)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 1 ] = Dmh [ 1 ] / (1. - dt*( (( - 2.0 * am - bm - ah ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 2 ] = Dmh [ 2 ] / (1. - dt*( (( - am - 2.0 * bm - ah ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 3 ] = Dmh [ 3 ] / (1. - dt*( (( - 3.0 * bm - ah ))*(1.0) + (am)*(1.0) + (bh)*(1.0) )) ;
Dmh [ 4 ] = Dmh [ 4 ] / (1. - dt*( (( - 3.0 * am - bh ))*(1.0) + (bm)*(1.0) + (ah)*(1.0) )) ;
Dmh [ 5 ] = Dmh [ 5 ] / (1. - dt*( (( - 2.0 * am - bm - bh ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (ah)*(1.0) )) ;
Dmh [ 6 ] = Dmh [ 6 ] / (1. - dt*( (( - am - 2.0 * bm - bh ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (ah)*(1.0) )) ;
Dmh [ 7 ] = Dmh [ 7 ] / (1. - dt*( (( - 3.0 * bm - bh ))*(1.0) + (am)*(1.0) + (ah)*(1.0) )) ;
}
/*END CVODE*/
static int states (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) { {
mh [ 0 ] = mh [ 0 ] + (1. - exp(dt*((( - 3.0 * am - ah ))*(1.0) + (bm)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(am) - ah ))*(1.0) + (bm)*(1.0) + (bh)*(1.0) ) - mh) ;
mh [ 1 ] = mh [ 1 ] + (1. - exp(dt*((( - 2.0 * am - bm - ah ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( (- 2.0)*(am) - bm - ah ))*(1.0) + ((3.0)*(am))*(1.0) + ((2.0)*(bm))*(1.0) + (bh)*(1.0) ) - mh) ;
mh [ 2 ] = mh [ 2 ] + (1. - exp(dt*((( - am - 2.0 * bm - ah ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( - am - (2.0)*(bm) - ah ))*(1.0) + ((2.0)*(am))*(1.0) + ((3.0)*(bm))*(1.0) + (bh)*(1.0) ) - mh) ;
mh [ 3 ] = mh [ 3 ] + (1. - exp(dt*((( - 3.0 * bm - ah ))*(1.0) + (am)*(1.0) + (bh)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(bm) - ah ))*(1.0) + (am)*(1.0) + (bh)*(1.0) ) - mh) ;
mh [ 4 ] = mh [ 4 ] + (1. - exp(dt*((( - 3.0 * am - bh ))*(1.0) + (bm)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(am) - bh ))*(1.0) + (bm)*(1.0) + (ah)*(1.0) ) - mh) ;
mh [ 5 ] = mh [ 5 ] + (1. - exp(dt*((( - 2.0 * am - bm - bh ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( (- 2.0)*(am) - bm - bh ))*(1.0) + ((3.0)*(am))*(1.0) + ((2.0)*(bm))*(1.0) + (ah)*(1.0) ) - mh) ;
mh [ 6 ] = mh [ 6 ] + (1. - exp(dt*((( - am - 2.0 * bm - bh ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( - am - (2.0)*(bm) - bh ))*(1.0) + ((2.0)*(am))*(1.0) + ((3.0)*(bm))*(1.0) + (ah)*(1.0) ) - mh) ;
mh [ 7 ] = mh [ 7 ] + (1. - exp(dt*((( - 3.0 * bm - bh ))*(1.0) + (am)*(1.0) + (ah)*(1.0))))*(- ( 0.0 ) / ( (( (- 3.0)*(bm) - bh ))*(1.0) + (am)*(1.0) + (ah)*(1.0) ) - mh) ;
}
return 0;
}
and I get the following error:
Code: Select all
"/usr/local/nrn/share/nrn/libtool" --mode=compile gcc -DHAVE_CONFIG_H -I. -I.. -I"/usr/local/nrn/include/nrn" -I"/usr/local/nrn/x86_64/lib" -g -O2 -c -o hh58Fa.lo `test -f 'hh58Fa.c' || echo '/'`hh58Fa.c
gcc -DHAVE_CONFIG_H -I. -I.. -I/usr/local/nrn/include/nrn -I/usr/local/nrn/x86_64/lib -g -O2 -c hh58Fa.c -fPIC -DPIC -o .libs/hh58Fa.o
hh58Fa.c: In function 'states':
hh58Fa.c:246: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c:247: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c:248: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c:249: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c:250: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c:251: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c:252: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c:253: error: invalid operands to binary - (have 'double' and 'double *')
hh58Fa.c: In function '_thread_cleanup':
hh58Fa.c:327: warning: incompatible implicit declaration of built-in function 'free'
make: *** [hh58Fa.lo] Error 1
I've done this both in Linux and Windows, getting the same result. Did I do something wrong fixing the c file? Apparently, ode_matsol1 is OK, the problem seems to be in states.
Thanks a lot