CVode

classes
   ModelDescriptionIssues  event  re_init       statistics     
   active         jacobian       record         stiff          
   atol           maxorder       record_init    store_event    
   atolscale      maxstep        record_remove  use_daspk      
   current_method  minstep       rtol           use_local_dt   
   daspk_interp   netconlist     solve          use_mxb        
   debug_event    order          statename      
   dstates        print_event_queue  states     

SYNTAX

objref cvode
cvode = new CVode()
cvode = new CVode(0 - 1)

DESCRIPTION

Multi order variable time step integration method which may be used in place of the default staggered fixed time step method. The performance benefits can be substantial (factor of more than 10) for problems in which all states vary slowly for long periods of time between fast spikes.

This class interfaces a variant of the CVODE Software Package by Scott D. Cohen and Alan C. Hindmarsh obtained from http://www.netlib.org. The changes to this package (located in $NEURONHOME/src/cvode/cvodemlh.[ch]) have to do with handling abrupt changes in parameters (eg current stimulus pulses, synaptic conductance onset, etc).

When this class is active the finitialize/fadvance calls use the CVode integrator. The argument to the CVode constructor defines whether the instance starts out active (1) or inactive (0). With no argument the instance is inactive. Only one instance of cvode can exist at one time at present. In the default variable step context, the integrator chooses the time step and fadvance returns after one step. Local accuracy is determined by atol and rtol .

The two major energy barriers to using the method are the requirement that hh type models be expressed in a DERIVATIVE block (instead of the explicit exponential integration step commonly implemented in a PROCEDURE) and that my variant of the CVode solver be explicitly notified of the exact time of any discontinuity such as step changes, pulses, and synaptic conductance onset's. These issues are discussed in Channels and Events .

After your mod files are generalized it will probably be convenient to compare the default method with CVode by toggling the UseCVode checkbox in the RunControl .

SEE ALSO

fadvance , finitialize , secondorder , dt .

BUGS

The consequences of solving continuous equations can be sometimes surprising when one is used to discrete fixed time step simulations. For example if one records an action potential (with either fixed or variable time steps) and plays it back into a voltage clamp; the clamp potential is not a discrete function but an exact step function.

Only the SEClamp works with CVode. VClamp cannot be used with this method.

Also .mod authors must take measures to handle step changes in parameters ( Events )

BUGS

Alternative variable step methods such as use_local_dt and use_daspk have been folded into this class and it has become a catchall class for invoking any of the numerical methods. For example, use_mxb is used to switch between the tree structured matrix solver and the general sparse matrix solver. Not all components work together, see current_method for acceptable mixing.


solve

CVode

SYNTAX

cvode.solve()
cvode.solve(tout)

DESCRIPTION

With no argument integrates for one step. All states and assigned variables are consistent at time t. dt is set to the size of the step. With the tout argument, cvode integrates til its step passes tout. Internally cvode returns the interpolated values of the states (at exactly tout) and the CVode class calls the functions necessary to update the assigned variables. Note that cvode.solve(tout) may be called for any value of tout greater than t-dt where dt is the size of its last single step.

For backward compatibility with finitialize/fadvance it is better to use the active method instead of calling solve directly.


statistics

CVode

SYNTAX

cvode.statistics()

DESCRIPTION

Prints information about the number of integration steps, function evaluations, newton iterations, etc.


print_event_queue

CVode

SYNTAX

cvode.print_event_queue()
cvode.print_event_queue(Vector)

DESCRIPTION

With no arg, prints information on the event queue. It should only be called after an finitialize and before changing any aspect of the model structure. Many types of structure changes invalidate pointers used in the event queue.

With a vector argument, the delivery times are copied to the Vector in proper monotonically increasing order.


rtol

CVode

SYNTAX

x = cvode.rtol()
x = cvode.rtol(relative)

DESCRIPTION

Returns the local relative error tolerance. With arg, set the relative tolerance. The default relative tolerance is 0.

The solver attempts to use a step size so that the local error for each state is less than

	rtol*|state| + atol*atolscale_for_state
The error test passes if the error in each state, e[i], is such that e[i]/state[i] < rtol OR e[i] < atol*atolscale_for_state (the default atolscale_for_state is 1, see atolscale )


atol

CVode

SYNTAX

x = cvode.atol()
x = cvode.atol(absolute)

DESCRIPTION

Returns the default local absolute error tolerance. With args, set the default absolute tolerance. The default absolute tolerance is 1e-2. A multiplier for specific states may be set with the atolscale function and also may be specified in model descriptions.

The solver attempts to use a step size so that the local error for each state is less than

	rtol*|state| + atol*atolscale_for_state
The error test passes if the error in each state, e[i], is such that e[i]/state[i] < rtol OR e[i] < atol*atolscale_for_state

Therefore states should be scaled (or the absolute tolerance reduced) so that when the value is close to 0, the error is not too large.

(See atolscale for how to set distinct absolute multiplier tolerances for different states.)

Either rtol or atol may be set to 0 but not both. (pure absolute tolerance or pure relative tolerance respectively).


atolscale

CVode

SYNTAX

tol = cvode.atolscale(&var, toleranceMultiplier)
tol = cvode.atolscale(&var)

DESCRIPTION

Specifies the absolute tolerance scale multiplier (default is 1.0) for all STATE's of which the address of var is an instance. Eg. cvode.atolscale(&soma.v(.5), 1e-8) sets the absolute tolerance multiplier for all membrane potentials everywhere. (The syntax for merely specifying a name is admittedly cumbersome but the function is not often needed and it avoids the necessity of explicitly having to parse strings such as "TrigKSyn.G".) The currently specified multiplier for that state name is returned by the function call.

Specification of a particular STATEs absolute tolerance multiplier is only needed if its scale is extremely small or large and is best indicated within the model description file itself using the STATE declaration syntax:n

	state (units) <tolerance>
See nrn/demo/release/cabpump.mod for an example of a model which needs a specific scaling of absolute tolerances (ie, calcium concentration and pump density).


re_init

CVode

SYNTAX

cvode.re_init()

DESCRIPTION

Initializes the integrator. This is done by finitialize when cvode is active .


stiff

CVode

SYNTAX

x = cvode.stiff()
x = cvode.stiff(0-2)

DESCRIPTION

2 is the default. All states computed implicitly.

1 only membrane potential computed implicitly. 0 Adams-Bashforth integration.


active

CVode

SYNTAX

x = cvode.active()
x = cvode.active(0)
x = cvode.active(1)
following two not yet implemented
x = cvode.active(1, dt)
x = cvode.active(tvec)

DESCRIPTION

When CVode is active then finitialize calls re_init and fadvance calls solve .

This function allows one to toggle between the normal integration method and the CVode method with no changes to existing intepreter code. The return value is whether CVode is active.

With only a single 1 arg, the fadvance calls CVode to do a single variable time step.

With the dt arg, fadvance returns at t+dt.

With a Vector tvec argument, CVode is made active and a sequence of calls to fadvance returns at the times given by the elements of tvec. After the last tvec element, fadvance returns after each step.


maxorder

CVode

SYNTAX

x = cvode.maxorder()
x = cvode.maxorder(0 - 12)

DESCRIPTION

Default maximum order for implicit methods is 5. It is usually best to let cvode determine the order. 12 for Adams.


jacobian

CVode

SYNTAX

x = cvode.jacobian()
x = cvode.jacobian(0 - 2)

DESCRIPTION

0 is the default. Linear solvers supplied by NEURON. 1 use dense matrix 2 use diagonal matrix


states

CVode

SYNTAX

objref dest_vector
dest_vector = new Vector()
cvode.states(dest_vector)

DESCRIPTION

Fill the destination Vector with the values of the states. On return dest_vector.size will be the number of states.


dstates

CVode

SYNTAX

cvode.dstates(dest_vector)
Fill the destination Vector with the values of d(state)/dt.

DESCRIPTION


statename

CVode

SYNTAX

strdef dest_string
cvode.statename(i, dest_string)
print dest_string

DESCRIPTION

Return the hoc name of the i'th string in dest_string


netconlist

CVode

SYNTAX

List = cvode.netconlist(precell, postcell, target)
List = cvode.netconlist(precell, postcell, target, list)

DESCRIPTION

Returns a new List (or appends to the list in the 4th argument position and returns a reference to that) of NetCon object references whose precell (or pre), postcell, and target match the pattern specified in the first three arguments. These arguments may each be either an object reference or a string. If an object, then each NetCon appended to the list will match that object exactly. String arguments are regular expressions and the NetCon will match if the name of the object has a substring that is accepted by the regular expression. (Object names are the internal names consisting of the template name followed by an index). An empty string, "", is equivalent to ".*" and matches everything in that field. A template name will match all the objects of that particular class. Note that some of the useful special regular expression characters are ".*+^$<>". The "<>" is used instead of the the standard special characters "[]" to specify a character range and obviates escaping the square bracket characters when attempting to match an array string. ie square brackets are not special and only match themselves.

EXAMPLES

A compact method of iterating over a set of NetCon objects is to create the list iterator
iterator ltr() {local i, cnt
	for i = 0, $o2.count - 1 {
		$o1 = $o2.object(i)
		iterator_statement
	}
}
and then take advantage of the automatic creation and destruction of lists with, for example, to print all the postcells that the given precell connects to:
objref xo
for ltr(xo, cvode.netconlist(precell, "", "")) {
	print xo.postcell
}


record

CVode

SYNTAX

cvode.record(&rangevar, yvec, tvec)
cvode.record(&rangevar, yvec, tvec, 1)

DESCRIPTION

Similar to the Vector record function but also works correctly with the local variable time step method. Limited to recording only range variables of density mechanisms and point processes.

During a run, record the stream of values in the specified range variable into the yvec Vector along with time values into the tvec Vector. Note that each recorded range variable must have a separate tvec which will be different for different cells. On initialization the yvec and tvec Vectors are resized to 1 and the initial value of the range variable and time is stored in the Vectors.

To stop recording into a particular vector, remove all the references either to tvec or yvec or call record_remove .

If the fourth argument is present and equal to 1, the yvec is recorded only at the existing t values in tvec. This option may slow integration since it requires calculation of states at those particular times.


record_init

CVode

SYNTAX

flag = cvode.record_init()
flag = cvode.record_init(anything_but_2)
flag = cvode.record_init(2)

DESCRIPTION

The default value of the flag is 2. In this case a call to finitialize initializes the Vectors being recorded in the list created by record as well as their associated time vectors so that they are all resized to 0 and the time and value is recorded.

With no arg or arg != 2 then initialization of the Vectors happens only on the next finitialize. i.e all results of runs append to the vectors until the next explicit cvode.record_init


record_remove

CVode

SYNTAX

cvode.record_remove(yvec)

DESCRIPTION

Remove yvec (and the corresponding xvec) from the list of recorded vectors. See record .


event

CVode

SYNTAX

cvode.event(t)
cvode.event(t, "statement")

DESCRIPTION

With no argument, an event without a source or target is inserted into the event queue for "delivery" at time t. This has the side effect of causing a return from fadvance() exactly at time t. This is used by the stdrun.hoc file to make sure a simulation stops at tstop or after the approprieate time on pressing "continuerun" or "continuefor"

If the hoc statement argument is present, the statement is executed when the event time arrives. This statement is normally a call to a procedure which may send another cvode.event. Note that since the event queue is cleared upon finitialize() the cvode.event must be sent after that.


minstep

CVode

SYNTAX

hmin = cvode.minstep()
hmin = cvode.minstep(hmin)

DESCRIPTION

Gets (and sets in the arg form) the minimum time step allowed for a CVODE step. Default is 0.0 . An error message is printed if a time step less than the minimum step is used.

BUGS

Not very useful. What we'd really like is a minimum first order implicit step.


maxstep

CVode

SYNTAX

hmax = cvode.maxstep()
hmax = cvode.maxstep(hmax)

DESCRIPTION

Gets (and sets in the arg form) the maxium value of the step size allowed for a CVODE step. CVODE will not choose a step size larger than this. The default value is 0 and in this case means infinity.


use_local_dt

CVode

SYNTAX

boolean = cvode.use_local_dt()
boolean = cvode.use_local_dt(boolean)

DESCRIPTION

Gets (and sets) the local variable time step method flag. When CVODE is active this implies a separate CVODE instance for every cell in the simulation. record is the only way at present that variables can be properly obtained when this method is used.

BUGS

Not well integrated with the existing standard run system graphics because cells are generally at different times and an fadvance only changes the variables for the earliest time cell.

use_daspk and use_local_dt cannot both be 1 at present. Toggling one on will toggle the other off.


debug_event

CVode

SYNTAX

cvode.debug_event(1)
cvode.debug_event(2)

DESCRIPTION

Prints information whenever an event is generated or delivered. When the argument is true, information is printed at every integration step as well.


order

CVode

SYNTAX

order = cvode.order()
order = cvode.order(i)

DESCRIPTION

CVODE method order used on the last step. The arg form is for the ith cell instance with the local step method.


use_daspk

CVode

SYNTAX

boolean = cvode.use_daspk()
boolean = cvode.use_daspk(boolean)

DESCRIPTION

Gets (sets for the arg form) the internal flag with regard to whether to use the daspk method when CVode is active If CVode is active and the simulation involves LinearMechanism or extracellular mechanisms then the daspk method is automatic and required.


current_method

CVode

SYNTAX

method = cvode.current_method()

DESCRIPTION

A value that indicates

modeltype + 10*use_sparse13 + 100*methodtype + 1000*localtype

where modeltype has the value: 0 if there are no sections or LinearMechanisms (i.e. empty model) 2 if the extracellular mechanism or LinearMechanism is present. (in this case the fully implicit fixed step or daspk methods are required and cvode cannot be used. 1 otherwise

use_sparse13 is 0 if the tree structured matrix solver is used and 1 if the general sparse matrix solver is used. The latter is required for daspk and not allowed for cvode. The fixed step methods can use either. The latter takes about twice as much time as the former.

methodtype = secondorder if CVode is not active. It equals 3 if CVODE is being used and 4 is DASPK is used.

localtype = 1 if the local step method is used. This implies methodtype==3


use_mxb

CVode

SYNTAX

boolean = cvode.use_mxb()
boolean = cvode.use_mxb(boolean)

DESCRIPTION

Switch between the tree structured matrix solver (0) and the general sparse matrix solver (1). Either is acceptable for fixed step methods. For CVODE only the tree structured solver is allowed. For DASPK only the general sparse solver is allowed.


daspk_interp

CVode

SYNTAX

cvode.daspk_interp(tout)

DESCRIPTION

Moves the simulation to tout which must be in the interval of the last daspk step. Do not use this.


store_event

CVode

SYNTAX

cvode.store_event(vec)

DESCRIPTION

Accumulates all the sent events as adjacent pairs in the vector. The pairs are the time at which the event was sent and the time it is to be delivered. The user should do a vec.resize(0) before starting a run. Cvode will stop storing with cvode.store_event(). This is primarily for gathering data to design more efficient priority queues. It may be eliminated when the tq-exper branch is merged back to the main branch. Notice that there is no info about event type or where the event is coming from or going to.


ModelDescriptionIssues

CVode
   Channels       Concentrations  Events        
The following aspects of model descriptions (.mod files) are relevant to their use with CVode.

KINETIC block - No changes required.

DERIVATIVE block - No changes required. The Jacobian is approximated as a diagonal matrix. If the states are linear in state' = f(state) the diagonal elements are calculated analytically, otherwise the diagonal elements are calculated using the numerical derivative (f(s+.01) - f(s))/.001 .

LINEAR, NONLINEAR blocks - No changes required. However, at this time they can only be SOLVED from a PROCEDURE or FUNCTION, not from the BREAKPOINT block. The nrn/src/nrnoc/vclmp.mod file gives an example of correct usage in which the function icur is called from the BREAKPOINT block and in turn SOLVE's a LINEAR block. If desired, it will be a simple matter to allow these blocks to be solved from the BREAKPOINT block.

SOLVE PROCEDURE within a BREAKPOINT block - Changes probably required. Such a procedure is called once after each return from CVode.solve().


Channels

ModelDescriptionIssues
The SOLVE PROCEDURE form was often used to implement the exponential integration method for HH like states and was very efficient in the context of the Crank-Nicholson like staggered time step approach historically used by NEURON. Furthermore the exponential integration often used tables of rates which were calculated under the assumption of a fixed time step, dt. Although it can still be used under some circumstances, the usage to integrate states should be considered obsolete and converted to a DERIVATIVE form. To do this, 1) replace the PROCEDURE block with a DERIVATIVE block, eg. DERIVATIVE states { m' = (minf - m)/mtau ... } 2) replace the SOLVE statement in the BREAKPOINT block with SOLVE states METHOD cnexp 3) if using tables, store mtau instead of (1 - exp(-dt/mtau)) The nmodl translator will emit c code for both the staggered time step and high order variable time step methods. The only downside is slightly less efficiency with the staggered time step method since the exp(-dt...) is calculated instead of looked up in tables.

In summary, no model should anymore depend on dt.


Concentrations

ModelDescriptionIssues


Events

ModelDescriptionIssues

How does one handle events? This is really the only serious difficulty in writing models that work properly in the context of a variable time step method. All models which involve discontinuous functions of time, eg steps, pulses, synaptic onset, require special provision to notify the integrator that an event has occurred within this time step, ie between t-dt and t. If this is not done, the time step may be so large that it completely misses a pulse or synaptic event. And if it does see the effect of the event, there is a huge inefficiency involved in the variable step method's search for the location of the event and the concomitant tremendous reduction in size of dt.

So, if you change any variable discontinuously in the model at some time tevent, call call

        at_time(tevent)
The user may check the return value of this function to decide if something needs changing. Examples of the two styles of usage are:

1) Just notify and do the logic separately.

	at_time(del)
	at_time(del + dur)
	if (t >= del && t <= dur) {
		istim = on_value
	}else{
		istim = 0
	}

2) Use the at_time return value to do the logic.

INITIAL {
	istim = 0
}
...
	if (at_time(del)) {
		istim = on_value
	}
	if (at_time(del + dur)) {
		istim = 0
	}
Notice the requirement of initialization or else if the previous run was stopped before del + dur the value of istim would be on_value at the beginning of the next run.

What happens internally when at_time(tevent) is called?

The interesting case (t-dt < tevent <= t) --- First, at_time returns 0. Then CVode changes its step size to (tevent - (t-dt) - epsilon) and redoes the step starting at t-dt. Note that this should be safely prior to the event (so at_time still returns 0), but if not then the above process will repeat until a step size is found for which there is no event. CVode then re-initializes it's internal state and restarts from a new initial condition at tevent+epsilon. Now when at_time is called, it returns 1. Note that in its single step mode, CVode.solve() will return at t = tevent-epsilon, the subsequent call will start the initial condition at t = tevent + epsilon and return after a normal step (usually quite small).

The case (tevent <= t - dt) --- at_time returns 0.

The case (tevent > t) --- at_time returns 0.

Note that an action potential model with axonal delay delivering a "message" to a synaptic model may or may not think it worthwhile to call at_time at the time of threshold (I would just do my own interpolation to set t_threshold) but will certainly call at_time(t_threshold + delay) (and possibly not allow t_threshold to change again until it returns a 1);

I am sorry that the variable time step method requires that the model author take careful account of events but I see no way to have them automatically taken care of.


neuron/neuron/classes/cvode.hel : 21918 Mar 31