STATEs which are arrays

NMODL and the Channel Builder.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: STATEs which are arrays

Post by hines »

You changed

Code: Select all

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) ;
to

Code: Select all

 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) ;
and forgot to change the last token in the line from mh to mh[0]
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: STATEs which are arrays

Post by patoorio »

Doh!
Now it works. Thanks a lot!!
patoorio
Posts: 81
Joined: Wed Jan 30, 2008 12:46 pm

Re: STATEs which are arrays

Post by patoorio »

Well, it compiled but didn't work properly. There were more errors in the .c file that permitted compilation but made a completely different mechanism.

To illustrate it better, I experimented with a simpler model; 4 states in a line. In the 'non-array' version, the DERIVATIVE block is

Code: Select all

DERIVATIVE states {
	rates(v)
	m0' = -3*am*m0 + bm*m1
	m1' = -(2*am+bm)*m1 + 3*am*m0 + 2*bm*m2
	m2' = -(am+2*bm)*m2 + 2*am*m1 + 3*bm*m3
	m3' = -3*bm*m3 + am*m2
}
while the DERIVATIVE block in the 'array' version is

Code: Select all

DERIVATIVE states {
	rates(v)
	m'[0] = -3*am*m[0] + bm*m[1]
	m'[1] = -(2*am+bm)*m[1] + 3*am*m[0] + 2*bm*m[2]
	m'[2] = -(am+2*bm)*m[2] + 2*am*m[1] + 3*bm*m[3]
	m'[3] = -3*bm*m[3] + am*m[2]
}
Quite similar, aren't they? Well, the _ode_matsol1 and states functions generated for the 'non-array' mechanism are

Code: Select all

 static int _ode_matsol1 (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) {
 rates ( _threadargscomma_ v ) ;
 Dm0 = Dm0  / (1. - dt*( (- 3.0 * am)*(1.0) )) ;
 Dm1 = Dm1  / (1. - dt*( (- ( 2.0 * am + bm ))*(1.0) )) ;
 Dm2 = Dm2  / (1. - dt*( (- ( am + 2.0 * bm ))*(1.0) )) ;
 Dm3 = Dm3  / (1. - dt*( (- 3.0 * bm)*(1.0) )) ;
}
 /*END CVODE*/
 static int states (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) { {
   rates ( _threadargscomma_ v ) ;
    m0 = m0 + (1. - exp(dt*((- 3.0 * am)*(1.0))))*(- ( (bm)*(m1) ) / ( ((- 3.0)*(am))*(1.0) ) - m0) ;
    m1 = m1 + (1. - exp(dt*((- ( 2.0 * am + bm ))*(1.0))))*(- ( ((3.0)*(am))*(m0) + ((2.0)*(bm))*(m2) ) / ( (- ( (2.0)*(am) + bm ))*(1.0) ) - m1) ;
    m2 = m2 + (1. - exp(dt*((- ( am + 2.0 * bm ))*(1.0))))*(- ( ((2.0)*(am))*(m1) + ((3.0)*(bm))*(m3) ) / ( (- ( am + (2.0)*(bm) ))*(1.0) ) - m2) ;
    m3 = m3 + (1. - exp(dt*((- 3.0 * bm)*(1.0))))*(- ( (am)*(m2) ) / ( ((- 3.0)*(bm))*(1.0) ) - m3) ;
   }
  return 0;
}
and for the 'array' mechanism I get

Code: Select all

 static int _ode_matsol1 (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) {
 rates ( _threadargscomma_ v ) ;
 Dm [ 0 ] = Dm  / (1. - dt*( (- 3.0 * am)*(1.0) + (bm)*(1.0) )) ;
 Dm [ 1 ] = Dm  / (1. - dt*( (- ( 2.0 * am + bm ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0) )) ;
 Dm [ 2 ] = Dm  / (1. - dt*( (- ( am + 2.0 * bm ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0) )) ;
 Dm [ 3 ] = Dm  / (1. - dt*( (- 3.0 * bm)*(1.0) + (am)*(1.0) )) ;
}
 /*END CVODE*/
 static int states (double* _p, Datum* _ppvar, Datum* _thread, _NrnThread* _nt) { {
   rates ( _threadargscomma_ v ) ;
   Dm [ 0  m = m + (1. - exp(dt*((- 3.0 * am)*(1.0) + (bm)*(1.0))))*(- ( 0.0 ) / ( ((- 3.0)*(am))*(1.0) + (bm)*(1.0) ) - m) ;
   Dm [ 1  m = m + (1. - exp(dt*((- ( 2.0 * am + bm ))*(1.0) + (3.0 * am)*(1.0) + (2.0 * bm)*(1.0))))*(- ( 0.0 ) / ( (- ( (2.0)*(am) + bm ))*(1.0) + ((3.0)*(am))*(1.0) + ((2.0)*(bm))*(1.0) ) - m) ;
   Dm [ 2  m = m + (1. - exp(dt*((- ( am + 2.0 * bm ))*(1.0) + (2.0 * am)*(1.0) + (3.0 * bm)*(1.0))))*(- ( 0.0 ) / ( (- ( am + (2.0)*(bm) ))*(1.0) + ((2.0)*(am))*(1.0) + ((3.0)*(bm))*(1.0) ) - m) ;
   Dm [ 3  m = m + (1. - exp(dt*((- 3.0 * bm)*(1.0) + (am)*(1.0))))*(- ( 0.0 ) / ( ((- 3.0)*(bm))*(1.0) + (am)*(1.0) ) - m) ;
   }
  return 0;
}
Besides the already discussed errors, there are several extra terms in the right side of the 'array' functions that are not present in the 'non-array' ones. So it actually takes a little more time to fix things (specially when dealing with an array of 8 states with a noise term).
Now why did I start bothering with all this? Because once I got the impression (I don't remember how) that when writing the things in the array form the mechanism would run faster. And it did in its compiled-but-incorrect form, but once I really made both mechanisms to be the same there was no speed-up. :(

One last question: Are all those parentheses necessary? And the multiplications by (1.0)? Would I see some speed improvement if I remove some of them?

Thanks a lot and best regards,

Patricio
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: STATEs which are arrays

Post by hines »

arrays are not faster because all variables are defines into a p array. (See the beginning of the translated c file).
The extra parentheses are not necessary and the multiply by (1.0) is not necessary. The hope is that the optimizer will
simplify to an optimum style. Whether the -O2 optimization in fact is doing anything useful is an experimental question.
Post Reply