Syntax:
cvode = h.CVode()
- 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.
Although for historical reasons, this class is called CVode at the hoc level, in fact it is an interface to a family of methods which are implemented on top of the CVODES and IDA integrators of the SUNDIALS package,
SUite of Nonlinear and DIfferential/ALgebraic equation Solvers Release 2.0.1, January 2005 Peter Brown, Aaron Collier, Keith Grant, Alan Hindmarsh, Steve Lee, Radu Serban, Dan Shumaker, Carol Woodward Center for Applied Scientific Computing, LLNL(see :meth:`CVode.use_local_dt` and :meth:`CVode.use_daspk`)
When this class is :meth:`CVode.active` the finitialize/fadvance calls use the CVode integrator. In the default variable step context, the integrator chooses the time step and fadvance returns after one step. Local accuracy is determined by :meth:`CVode.atol` and :meth:`CVode.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 the 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 :ref:`Channels <ModelDescriptionIssues_Channels>` and :ref:`Events <ModelDescriptionIssues_Events>`.
After your mod files are generalized it will probably be convenient to compare the default method with CVode by toggling the Use variable dt checkbox in the :ref:`VariableStepControl` panel :menuselection:`NEURON Main Menu --> Tools --> VariableStepControl`.
.. seealso::
:func:`fadvance`, :func:`finitialize`, :data:`secondorder`, :data:`dt`
Warning
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 (:ref:`Events <ModelDescriptionIssues_Events>`)
Warning
Alternative variable step methods such as :meth:`CVode.use_local_dt` and :meth:`CVode.use_daspk` have been folded into this class and it has become a catchall class for invoking any of the numerical methods. For example, :meth:`CVode.use_mxb` is used to switch between the tree structured matrix solver and the general sparse matrix solver. Not all components work together, see :meth:`CVode.current_method` for acceptable mixing.
Note
A from neuron import gui or h.load_file('stdrun.hoc') will create an instance called
h.cvode. Although this class is not strictly speaking a singleton, there is only one
integrator and it may be controlled and queried by any instance.
.. method:: CVode.solve
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 :func:`finitialize`/:func:`fadvance`
it is better to use the :meth:`CVode.active` method instead of calling
solve directly.
.. method:: CVode.statistics
Syntax:
``cvode.statistics()``
Description:
Prints information about the number of integration steps, function evaluations,
newton iterations, etc.
.. seealso::
:meth:`CVode.spike_stat`
.. method:: CVode.spike_stat
Syntax:
``cvode.spike_stat(vector)``
Description:
Similar to :meth:`CVode.statistics` but returns statistics information in the
passed :class:`Vector` argument. The vector will be resized to length
11 and the elements are:
.. code-block::
none
0 total number of equations (0 unless cvode has been active).
1 number of NetCon objects.
2 total number of events delivered.
3 number of NetCon events delivered.
4 number of PreSyn events put onto queue.
5 number of SelfEvents delivered.
6 number of SelfEvents put onto queue (net_send from mod files).
7 number of SelfEvents moved (net_move from mod files).
8 number of items inserted into event queue.
9 number of items moved to a new time in the event queue.
10 number of items removed from event queue.
.. note::
``vector`` must be an instance of :class:`Vector`
.. method:: CVode.print_event_queue
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 :class:`Vector` in
proper monotonically increasing order.
.. method:: CVode.event_queue_info
Syntax:
``cvode.event_queue_info(2, tvec, list)``
``cvode.event_queue_info(3, tvec, flagvec, list)``
Description:
Returns NetCon (2) or SelfEvent (3) information currently on the event queue.
If the type is 2, NetCon information currently on the event queue
is returned: delivery times are returned in tvec and the corresponding
NetCon objects are returned in the :class:`List` arg. If the type is 3,
SelfEvent information is returned: delivery times are returned in tvec,
the flags are returned in flagvec, and the SelfEvent targets
(ArtificialCells are PointProcesses) returned in the List arg.
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.
The delivery times are copied to the Vector in
proper monotonically increasing order.
.. note::
``list`` must be an instance of :class:`List`; you cannot use a Python list ``[]``.
.. method:: CVode.rtol
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
.. math::
(\mathrm{rtol}) |\mathrm{state}| + (\mathrm{atol})(\mathrm{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 :meth:`atolscale` )
.. method:: CVode.atol
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 :meth:`CVode.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
.. math::
(\mathrm{rtol}) |\mathrm{state}| + (\mathrm{atol})(\mathrm{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 :meth:`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).
.. method:: CVode.atolscale
Syntax:
**only works when called from HOC**
``tol = cvode.atolscale(ptr_var, toleranceMultiplier)``
``tol = cvode.atolscale(ptr_var)``
**works for both HOC and Python**
``tol = cvode.atolscale("basename" [, toleranceMultiplier])``
Description:
Specifies the absolute tolerance scale multiplier (default is 1.0)
for all STATE's of which the address
of var is an instance.
**Only the last form is currently supported in Python**; the first two forms
work from HOC but not Python.
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:
.. code-block::
none
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).
The "basename" form is simpler than the pointer form and was added to
simplify the implementation of the AtolTool. The pointer form required
the state to actually exist at the specified location. Base names are
``v``, ``vext``, state_suffix such as ``m_hh``, and PointProcessName.state such
as ``ExpSyn.g``.
.. method:: CVode.re_init
Syntax:
``cvode.re_init()``
Description:
Initializes the integrator. This is done by :func:`finitialize` when cvode
is :meth:`~CVode.active`.
.. method:: CVode.stiff
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.
.. method:: CVode.active
Syntax:
``x = cvode.active()``
``x = cvode.active(False or True)``
``x = cvode.active(0 or 1)``
**following two not yet implemented**
``x = cvode.active(True, dt)``
``x = cvode.active(tvec)``
Description:
When CVode is active then :func:`finitialize`
calls :meth:`CVode.re_init` and :func:`fadvance` calls :meth:`CVode.solve`.
This function allows one to toggle between the normal integration
method and the CVode method with no changes to existing interpreter
code. The return value is True is CVode is active; otherwise it is
False.
With only a single True (or 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.
.. method:: CVode.maxorder
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.
.. method:: CVode.jacobian
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
.. method:: CVode.states
Syntax:
.. code-block::
python
states_copy = h.Vector()
cvode.states(states_copy)
Description:
Fill the destination ``states_copy`` :class:`Vector` with the values of the states.
On return ``states_copy.size()`` will be the number of states.
.. method:: CVode.dstates
Syntax:
``cvode.dstates(dest_vector)``
Description:
Fill the destination :class:`Vector` with the values of d(state)/dt.
.. method:: CVode.f
Syntax:
``cvode.f(t, yvec, ypvec)``
Description:
returns f(yvec, t) in the :class:`Vector` ypvec. f is the existing model.
Size of yvec must be equal to the number of states ( ie vector size
returned by :meth:`CVode.states`). ypvec will be resized to the proper size.
Note that the order of the states in the vector is indicated by the
names returned by :meth:`CVode.statename`
.. warning::
Works only for global variable time step method.
Works only with single thread.
.. method:: CVode.yscatter
Syntax:
``cvode.yscatter(yvec)``
Description:
Fills the state variables with the values specified in the :class:`Vector` yvec.
Size of yvec must be equal to the number of states ( ie vector size
returned by :meth:`CVode.states`). Note that active CVode requires a subsequent
:meth:`CVode.re_init` if one wishes to integrate from the yvec state point.
.. warning::
Works only for global variable time step method.
Works only with single thread.
.. note::
``yvec`` must be a NEURON :class:`Vector` object. To scatter from an arbitrary Python iterable
``data`` (at the cost of an extra copy), use, e.g.
.. code-block::
python
h.CVode().yscatter(h.Vector(data))
.. method:: CVode.ygather
Syntax:
``cvode.ygather(yvec)``
Description:
Fills the :class:`Vector` yvec with the state variables (will be resized to the number of
states). This is analogous to :meth:`CVode.states` after a :meth:`CVode.re_init`.
.. warning::
Works only for global variable time step method.
Works only with single thread.
.. method:: CVode.fixed_step
Syntax:
``cvode.fixed_step()``
Description:
Uses the fixed step method to advance the simulation by :data:`dt` .
The initial condition is whatever state values are present (eg subsequent
to a previous integration step or :meth:`CVode.yscatter` or :meth:`CVode.f` or explicitly
user modified state values). The model state values are those after the
fixed step integration (but are NOT the same as the current state defined
by CVode and returned by :meth:`CVode.states` (that would be the case only after
a subsequent :meth:`CVode.re_init`)) To get the new current states in CVode
vector order, use :meth:`CVode.ygather`.
Valid under all circumstances. This is basically an :func:`fadvance` using
the fixed step method and avoids the overhead of
.. code-block::
python
h.CVode().active(False)
h.fadvance()
h.CVode().active(True)
in order to allow the use of the CVode functions assigning state and
evaluating states and dstates/dt; use via:
.. code-block::
python
h.CVode().fixed_step()
.. warning::
:meth:`CVode.dstates` are invalid and should be determined by a call to
:meth:`CVode.f` using the current state from :meth:`CVode.ygather` .
.. method:: CVode.error_weights
Syntax:
``cvode.error_weights(dest_vector)``
Description:
Fill the destination :class:`Vector` with the values of the weights used
to compute the norm of the local error in cvodes and ida.
.. method:: CVode.acor
Syntax:
``cvode.acor(dest_vector)``
Description:
Fill the destination :class:`Vector` with the values of the local errors
on the last step.
.. method:: CVode.statename
Syntax:
``cvode.statename(i, dest_string)``
``cvode.statename(i, dest_string, style)``
Description:
Return the HOC name of the i'th string in ``dest_string``, a NEURON string reference.
The default style, 0, is to attempt to specify the name in terms of
object references such as cell[3].syn[2].g. Style 1 specifies the name
in terms of the object id, eg. ExpSyn[25].g or Cell[25].soma.v(.5).
Style 2 returns the basename, e.g. v, or ExpSyn.g .
Example:
.. code-block::
python
from neuron import h
h.load_file('stdrun.hoc') # defines h.cvode
result = h.ref('')
soma = h.Section(name='soma')
h.cvode_active(True)
h.cvode.statename(0, result)
print(result[0])
The above code displays: ``soma.v(0.5)``
.. warning::
``dest_string`` must be a NEURON string reference (e.g. ``dest_string = h.ref('')``)
not a Python string, as those are immutable.
.. method:: CVode.netconlist
Syntax:
``List = cvode.netconlist(precell, postcell, target)``
``List = cvode.netconlist(precell, postcell, target, list)``
Description:
Returns a new :class:`List` (or appends to the list in the 4th argument
position and returns a reference to that) of :class:`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.
Example:
To print all the postcells that the given ``precell`` connects to:
.. code-block::
python
for nc in h.CVode().netconlist(precell, '', ''):
print(nc.postcell())
.. method:: CVode.record
Syntax:
``cvode.record(_ref_rangevar, yvec, tvec)``
``cvode.record(_ref_rangevar, yvec, tvec, 1)``
Description:
Similar to the Vector :meth:`~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 :class:`Vector` along with time values into the tvec :class:`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 :func:`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.
.. method:: CVode.record_remove
Syntax:
``cvode.record_remove(yvec)``
Description:
Remove yvec (and the corresponding xvec)
from the list of recorded :class:`Vector`s. See :meth:`record`.
.. method:: CVode.event
Syntax:
``cvode.event(t)``
``cvode.event(t, function)``
``cvode.event(t, function, pointprocess, re_init)``
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 :func:`fadvance` (or :meth:`CVode.solve` or :meth:`ParallelContext.psolve` or :func:`batch_run`
exactly at time t. This is used by the stdrun.hoc file
to make sure a simulation stops at tstop or after the appropriate
time on pressing "continuerun" or "continuefor". When :meth:`CVode.use_local_dt`
is active, all cells are interpolated to the event time.
If the hoc statement argument is present, the statement is executed (in
the object context of the call to cvode.event) 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 :func:`finitialize` the cvode.event must be sent after that.
Multiple threads and/or the local variable time step method, sometimes require
a bit of extra thought about the purpose of the statement. Should it be executed
only in the context of a single thread, should it be executed only in the
context of a single cell, and should only the integrator associated with that
cell be initialized due to a state change caused by the statement?
When the third arg is absent, then before the statement is executed, all cells
of all threads are interpolated to time t, all threads
join at time t, and the statement is executed by the main thread. A call to
:meth:`CVode.re_init` is allowed. If the third arg (a POINT_PROCESS object) is
present, then, the integrator of the cell (if lvardt) containing the POINT_PROCESS
is interpolated to time t, and the statement is executed by the thread
containing the POINT_PROCESS. Meanwhile, the other threads keep executing.
The statement should only access states and parameters associated with the
cell containing the POINT_PROCESS. If any states or parameters are changed,
then the fourth arg should be set to 1 to cause a re-initialization of only
the integrator managing the cell (:meth:`CVode.re_init` is nonsense in this context).
Example:
.. code-block::
python
from neuron import h, gui
def hi():
print('hello from hi, h.t = %g' % h.t)
h.finitialize(-65)
h.CVode().event(1.3, hi)
h.continuerun(2)
.. method:: CVode.minstep
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.
.. warning::
Not very useful. What we'd really like is a minimum first order implicit step.
.. method:: CVode.maxstep
Syntax:
``hmax = cvode.maxstep()``
``hmax = cvode.maxstep(hmax)``
Description:
Gets (and sets in the arg form) the maximum 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.
.. method:: CVode.use_local_dt
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 :meth:`~CVode.active`, this implies a separate CVODE
instance for every cell in the simulation. :meth:`CVode.record` is the only way
at present that variables can be properly obtained when this method is used.
.. warning::
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.
:meth:`CVode.use_daspk` and use_local_dt cannot both be 1 at present. Toggling one
on will toggle the other off.
.. method:: CVode.debug_event
Syntax:
``cvode.debug_event(1)``
``cvode.debug_event(2)``
Description:
Prints information whenever an event is generated or delivered. When the
argument is 2, information is printed at every integration step as well.
.. method:: CVode.use_long_double
Syntax:
``boolean = cvode.use_long_double()``
``booelan = cvode.use_long_double(boolean)``
Description:
When true, vector methods involving sums over the elements are accumulated
in a long double variable. This is useful in debugging when the
global variable time step method gives different results for different
:meth:`ParallelContext.nthread` or numbers of processes. It may be the case that the difference is
due to differences in round-off error due to the non-associativity of
computer addition. I.e when threads are used each thread adds up its own
group of numbers and then the group results are added together. When
a long double is used as the accumulator for addition, the round off error
is much more likely to be the same regardless of the order of addition. Note that
this DOES NOT make the simulation more accurate --- just more likely to be identical for
different numbers of threads or processes (if the difference without it was due to
round off errors during summation).
.. method:: CVode.order
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.
.. method:: CVode.use_daspk
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 IDA method when CVode is :meth:`~CVode.active`. If CVode is active
and the simulation involves :func:`LinearMechanism` or :func:`extracellular` mechanisms
then the IDA method is automatic and required.
Daspk refers to the Differential Algebraic Solver with the Preconditioned
Krylov method. The SUNDIALS package now calls this the IDA (Integrator
for Differential-Algebraic problems) integrator but it is really the same
thing.
.. method:: CVode.condition_order
Syntax:
``order = cvode.condition_order()``
``order = cvode.condition_order(1or2)``
Description:
When condition_order is 1 then :func:`NetCon` threshold detection takes place at a time
step boundary. This is the default. When condition_order is 2 then
NetCon threshold detection times are linearly interpolated within the
integration step interval for which the threshold occurred. Second order
threshold is limited to variable step methods and is ignored for the
fixed step methods. Note that second order threshold detection time may change
due to synaptic events within the interval or even be abandoned.
It is useful for cells with approach threshold very slowly or with large
time steps.
.. method:: CVode.dae_init_dteps
Syntax:
``eps = cvode.dae_init_dteps()``
``eps = cvode.dae_init_dteps(eps)``
``eps = cvode.dae_init_dteps(eps, style)``
Description:
The size of the "infinitesimal" fixed fully implicit step used for
initialization of the DAE solver, see :func:`use_daspk` , in order to
meet the the initial condition requirement of f(y',y,t)=0. The default
is 1e-9 ms.
The default heuristic for meeting the initial condition requirement based
on the pre-initialization value of all the states and an initialization time
of t0 is:
t = t0 Vector.play continuous.
Two dteps voltage solve steps. (does not change t, or membrane mechanism
states but changes v,vext).
The initial value of y is the present value of the
states.
t = t0 + dteps Vector.play continuous
One dteps step without changing y but it does determine dy/dt of the
v, vext portion of states.
t = t0 determine the dy/dt of the membrane mechanism states.
(note: membrane mechanism states are all derivative or kinetic
scheme states)
.. warning::
A number of things can go wrong with the heuristics used to provide
the integrator with a consistent initial condition. When this happens
the default behavior is to stop. However one can modify the error
handling and/or choose a second
initialization heuristic that might work by setting the style method.
The working values of style are 0,1,2, 8,9,10. the latter style group
(010 bit set) chooses the alternative heuristic. This alternative
is very similar to the default except the third dteps step that determines
y' also is allowed to change y. This may be more reliable when the user
is not using Vector.play continuous.
If the 1 or 2 bit is
set, a warning is printed instead of an error and the sim continues.
If the 2 bit is set, then for the next 1e-6 ms, the integrator solves the
equation f(y', y, t)*(1 - exp(-1e-7(t - t0)) where t0 is the initialization
time. I call this parasitic since it is supposed to be
analogous to every voltage having a small capacitance to ground.
It has not been determined if the parasitic
heuristic has a reliable mathematical basis and the user should investigate
the state change patterns in the neighborhood of the initialization time.
.. method:: CVode.simgraph_remove
Syntax:
``cvode.simgraph_remove()``
Description:
Removes all items from the list of Graph lines recorded during
a local variable step simulation. Graph lines would have been added to this
list with :ref:`gui_graph`.
.. method:: CVode.state_magnitudes
Syntax:
``cvode.state_magnitudes(integer)``
``cvode.state_magnitudes(Vector, integer)``
``maxstate = cvode.state_magnitudes("basename", _ref_maxacor)``
Description:
``cvode.state_magnitudes(1)`` activates the calculation of the
running maximum magnitudes of states and acor. 0 turns it off.
``cvode.state_magnitudes(2)`` creates an internal
list of the maximum of the maximum states and acors
according to the state basename currently in the model. Statenames not
in use have a maximum magnitude state and acor value of -1e9.
``maxstate = cvode.state_magnitudes("basename", _ref_maxacor)``
returns the maxstate and maxacor for the state type, e.g. "v" or
"ExpSyn.g", or "m_hh". Note: state type names can be determined from
MechanismType and MechanismStandard
``cvode.state_magnitudes(Vector, 0)`` returns all the maximum magnitudes for
each state in the Vector. This is analogous to cvode.states(Vector).
``cvode.state_magnitudes(Vector, 1)`` returns the maximum magnitudes for
each acor in the Vector.
.. method:: CVode.current_method
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 = :data:`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
.. method:: CVode.use_mxb
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.
.. method:: CVode.use_fast_imem
Syntax:
``boolean = cvode.use_fast_imem()``
``boolean = cvode.use_fast_imem(boolean)``
Description:
When true, compute i_membrane\_ for all segments during a simulation.
This is closely related to i_membrane which is computed when the
extracellular mechanism is inserted. However, i_membrane\_ (note
the trailing '\_'), has dimensions of nA instead of mA/cm2 (ie. total
membrane current out of the segment), is available
at 0 area nodes (locations 0 and 1 of every section), does not require
that extracellular be inserted (and so is much faster), and works
during parallel simulations with variable step methods. (ie. does not
require IDA which is currently not available in parallel).
i_membrane\_ exists as a range variable only when ``use_fast_imem`` has
been called with an argument of 1. Conversely, i_membrane\_ is
not computed when ``use_fast_imem`` is not called or with an
argument of 0.
i_membrane\_ include capacity current and all transmembrane
ionic currents but not stimulus currents. POINT_PROCESS synaptic
currents are considered ionic currents and so are included
in i_membrane\_. From charge conservation
a fundamental property is that the sum of all i_membrane\_ is
identical to the sum of all ELECTRODE_CURRENT (Current cannot
flow axially out of a cell since the root and leaves of each
cell tree have sealed end boundary conditions.)
The following tests this conservation law, assuming that the only
ELECTRODE_CURRENTs are IClamp. Note the idiom that visits all segments
of a model but only once each segment to sum up i_membrane\_
.. code:: python
from neuron import h
h.CVode().use_fast_imem(1)
def assert_whole_model_charge_conservation():
# sum over all membrane current
total_imem = 0.0
for sec in h.allsec():
for seg in sec.allseg(): # also the 0 area nodes at 0 and 1
if seg.x == sec.orientation() and sec.parentseg() is not None:
continue # skip segment shared with parent
total_imem += seg.i_membrane_
# sum over all ELECTRODE_CURRENT (if only using IClamp)
total_iclamp_cur = sum(ic.i for ic in h.List('IClamp'))
print("total_imem=%g total_iclamp_cur=%g" % (total_imem, total_iclamp_cur))
assert(abs(total_imem - total_iclamp_cur) < 1e-12)
In the above fragment ``sec.parentseg()`` is needed to count
the root and use of ``sec.trueparentseg()`` would count all sections
that connect to the root section at 0 because all those sections have
a trueparentseg of None.
Also, although an extremely rare edge case, ``sec.orientation()``
is needed to match which segment is closest to root.
.. method:: CVode.store_events
Syntax:
``cvode.store_events(vec)``
Description:
Accumulates all the sent events as adjacent pairs in the :class:`Vector` vec.
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.
.. method:: CVode.queue_mode
Syntax:
``mode = cvode.queue_mode(boolean use_fixed_step_bin_queue, boolean use_self_queue)``
Description:
Normally, there is one event queue for all pending events. However, for the
fixed step method one can obtain marginally better queue performance through
the use of a bin queue for NetCon events. This utilizes a queue with
bins of size dt which has a very fast insertion time and every time step
all the events in a bin are delivered to their targets. Note that the
numerics of the simulation will differ compared to the default splay
tree queue (which stores double precision delivery times) if
NetCon.delay values are not integer multiples of dt. Also, even with
the fixed step method and and delays as integer multiples of dt, results
can differ at the double precision round off level due to the different order
that same time events can be received by the NET_RECEIVE block.
The optional "use_self_queue" (default 0) argument can only be used if the
the simulation is run with :meth:`~ParallelContext.psolve` method
of the :class:`ParallelContext` and must be selected prior to a call of
:meth:`ParallelContext.set_maxstep` since this special technique requires a
computation of the global minimum :meth:`NetCon.delay` (not just the
minimum interprocessor NetCon delay) and that delay must be
greater than 0. The technique avoids the use of the normal splay tree queue
for self events for ARTIFICIAL_CELLs (events initiated by the net_send call
and which may be manipulated by the net_move call in the NET_RECEIVE block).
It may thus be considerably faster. However, every minimum NetCon delay interval,
all the ARTIFICIAL_CELLS must be iterated to see if there are any outstanding
net_send events that need to be handled. Thus it is likely to have a beneficial
performance impact only for large numbers of ARTIFICIAL_CELLs which receive
many external input events per reasonable minimum delay interval. This method
has not receive much testing and the results should be compared with the
default queuing method.
Returns ``2*use_self_queue + use_fixed_step_bin_queue``.
.. seealso::
:meth:`ParallelContext.spike_compress`
.. method:: CVode.structure_change_count
Syntax:
``intcnt = cvode.structure_change_count()``
Description:
Returns the integer internal value of structure_change_cnt.
Structure_change_cnt is internally incremented whenever the
low level computable structures of the model have been setup
due to a change in number of segments, sections, topology, etc,
and some internal function requires that the computable structures
are consistent with the user level description of the model such as
finitialize, fadvance, define_shape, and many others.
.. method:: CVode.diam_change_count
Syntax:
``cnt = cvode.diam_change_count()``
Description:
Returns the integer internal value of diam_change_cnt.
Diam_change_cnt is internally incremented whenever some internal
function checks the diam_changed flag and calls the internal
recalc_diam() function.
.. method:: CVode.extra_scatter_gather
Syntax:
``cvode.extra_scatter_gather(direction, pycallable)``
Description:
If the direction is 0, the pycallable is
called immediately AFTER cvode has scattered its state variables.
If the direction is 1, the pycallable is called immediately BEFORE
cvode gathers the values of the state variables.
For the fixed step method, the direction 0 pycallable is called
after voltages have been updated and immediately before the
nonvint part (before DERIVATIVE, KINETIC, etc. blocks). It is also
called during cvode.re_init() when cvode is inactive.
.. warning::
Works only for fixed and global variable time step methods.
Works only with single thread.
Example of setting and removing, with arguments:
.. code::
from neuron import h
def hello1(cort_secs):
print('hello1')
cort_secs.append('corticalcell')
def hello2(arg):
print('hello2', arg)
cort_secs = []
recording_callback = (hello1, cort_secs)
# declaring a function to run with every fadvance
h.CVode().extra_scatter_gather(0, recording_callback)
h.finitialize(-65)
h.fadvance()
h.fadvance()
# removing the previous function
h.CVode().extra_scatter_gather_remove(recording_callback)
print('---')
# declaring a new function to run with each fadvance
recording_callback = (hello2, cort_secs)
h.CVode().extra_scatter_gather(0, recording_callback)
h.finitialize(-65)
h.fadvance()
h.fadvance()
.. method:: CVode.extra_scatter_gather_remove
Syntax:
``cvode.extra_scatter_gather_remove(pycallable)``
Description:
Removes the pycallable from list of callbacks used when cvode
scatters its state variables or gathers its dstate variable.
.. method:: CVode.cache_efficient
Syntax:
``mode = cvode.cache_efficient(True or False)``
Description:
When set, G*v = R matrix and vectors are reallocated in tree order so that
all the elements of each type are contiguous in memory. Pointers to these
elements used by the GUI, Vector, Pointer, etc. are updated.
Much of the implementation was contributed by Hubert Eichner <eichnerh@in.tum.de>
:meth:`ParallelContext.multisplit` automatically sets h.CVode().cache_efficient(True)
0 or 1 can be used instead of ``False`` or ``True``, respectively.
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().
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-Nicolson 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,
replace the PROCEDURE block with a DERIVATIVE block, eg.
DERIVATIVE states { m' = (minf - m)/mtau ... }replace the SOLVE statement in the BREAKPOINT block with
SOLVE states METHOD cnexpif using tables, store mtau instead of (1 -\exp(-dt/m_{tau})) 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 :data:`dt`.
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:
Just notify and do the logic separately.
at_time(del) at_time(del + dur) if t >= del and t <= del + dur: istim = on_value else: istim = 0Use 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.