Object Type:	hsolve

Description:	An object used for implementing the Hines method for the
		implicit solution of dendritic trees, as well as a set of
		routines that maximize speed.  This allows faster, more
		stable numerical integration methods to be used with GENESIS,
                particularly when there are many compartments in a cell.

Author:		U. S. Bhalla and E. De Schutter Caltech 91-93,
                E. De Schutter BBF-UIA 94-99.


------------------------------------------------------------------------------

ELEMENT PARAMETERS

DataStructure:	hsolve_type  [in src/hines/hines_struct.h]

Size:   552 bytes

Fields: path		wildcard list of compartments to be
			integrated; specifies all compartments
			belonging to cell which hsolve element will
			deal with.  The default is
			./##[][TYPE=compartment]

	chanmode	flag that controls channel computations and
			chip-array use.

			0 = default. It is the most general, uses least memory
                        and is slowest (no chip-array).  It takes over the
                        actions of compartments only, computing all other
                        object types as before.  As a consequence, all
                        computed fields of the original elements are updated,
                        and all user-setable fields may be set, just as before
                        the element was taken over by the hsolve element.
                        This means that you can add and delete outgoing
                        messages to compartments or other elements whenever
                        you like and easily change parameters during the
                        course of the simulation, making this the easiest mode
                        to use.  It is also the most compatible mode of
                        operation, and is guaranteed to work with any future
                        new object type.  One significant limitation is that
                        you cannot add or delete AXIAL, RAXIAL or CHANNEL
                        messages, once the hsolve element has been created.
                        In order to calculate the compartment Im field,
                        the compute_Im flag must be set.

			1 = like chanmode 0, but optimizes calculations for
                        tabchannel equations.  This will be faster than
                        chanmode 0 for models containing tabchannels and
                        tab2Dchannels.

			2 = assumes integer exponents (maximum = 6) for
                        tabchannel and tab2Dchannel gate variables.  Chanmodes
                        2-5 are the fastest because of the large 'chip-array'.
                        Under these moses, you can no longer assume that all
                        the fields of the elements that are taken over by
                        hsolve will be updated.  Incoming and outgoing
                        messages to and from the disabled elements will work
                        properly, provided that they were added prior to
                        setting up the hsolve element.  However, the Vm fields
                        of all compartments are automatically updated under
                        chanmode 2, whether or not there were pre-existing
                        messages to output Vm.  Note also that several fields
                        (Gk, Ik, Ek, and Im) are not available for output in
                        chanmodes 2 and 3, so you have to use chanmode 4 or 5
                        with findsolvefield if you want to output these fields.

                        3 = as chanmode 2, but Vm fields of compartments are
                        not automatically updated.  This mode is faster than
                        chanmode 2, at the expense of somewhat greater memory
                        usage.  All disabled elements (including compartments)
                        that had outgoing messages to non-hsolved objects
                        during the SETUP call are updated at a rate determined
                        by outclock.  Thus SAVE and PLOT messages will work
                        properly, except with the fields Gk, Ik, Ek, and Im.

                        4 = as chanmode 3, but more variables (e.g. Gk, Ik,
                        Ek, Im and leak) are stored by hsolve so that they can
                        be accessed by SAVE or PLOT messages (see
                        findsolvefield documentation for more details).  The
                        added field, leak, is available for each compartment
                        to give the current flowing through the membrane
                        resistance Rm.  Im is calculated directly as the sum
                        of the channel currents and the leakage current.

                        5 = as chanmode 4, but normalizes the currents and
                        conductances stored in the hsolve givals array (for
                        display purposes only).

        computeIm       flag to determine whether the compartment Im field
                        is calculated in chanmodes 0 and 1.  If set, Im is
                        approximated by the sum of axial currents and injected
                        current, as done for non-hsolved compartments and
                        symcompartments.

                        0 = the default, does not calculate Im.
                        1 = calculates Im, with a decrease in speed.
                        
        comptmode	flag that selects option for compartment computations

                        0 = uses least memory, and is slightly slower. However
                        the amount of time spent in this part is usually less
                        than 10% of the total, so this should not matter.

			1 = the default, uses a lot of memory and is twice as
                        fast for the actual Hines matrix calculation.  This
                        may be only a 5% difference.  comptmode 1 is
                        required for chanmodes 2-5.

        calcmode        flag affecting operations for chanmodes 2-5.

                        0 = no interpolation will be performed in lookup
                        tables.  This mode is for backward compatibility with
                        older versions of hsolve (prior to GENESIS 2.1) which
                        did not use interpolation.

                        1 = default mode with linear interpolation of values
                        in lookup tables.

        storemode       flag, only valid for chanmodes 4 and 5,  to allow the
                        output of total currents and conductances.  For each
                        type of voltage gated channel in the model, the total
                        currents or conductances are the sum of the
                        corresponding Ik or Gk fields for all compartments
                        where the channel is present.  This assumes that these
                        channels have the same name in each compartment.
                        These are stored in an array called itotal.  When the
                        hsolve element is set up, a message will be output
                        giving a list of channel names and corresponding
                        itotal indices, if the 'silent' command has previously
                        been given with a negative argument. 

                        0 = no sums are stored.
                        1 = total currents are stored.
                        2 = total conductances are stored.

	no_elminfo	= 0/1; flag that controls whether the HPUT and HSET
			actions will work.  Is relevant only for chanmodes
			2-5.  Default is zero (HPUT and HSET work).
			Saves memory if non-zero value is used (HPUT, HSET
			do not work; HRESTORE and HSAVE do).

	outclock	number [0-99] of the clock used for all element
                        updates if chanmode = 3, 4, or 5, will affect SAVE and
                        PLOT messages.  The element update routines are not
                        very efficient, so setting them to a slower clock will
                        speed up things.  Note however that any object that is
                        not computed by the hines solver and that depends on a
                        message from a hsolve-computed object will be affected.

	The other fields displayed with the showobject command are NOT
        to be set by the user.

------------------------------------------------------------------------------

SIMULATION PARAMETERS

Function:	HinesSolver  [in src/hines/hsolve.c]

Classes:	hsolver

Actions:	RESET	does the standard reset functions, but also for
			chanmodes 2-5, it will update all parameters in
			the chip-array (equivalent to a HRESTORE call) and
			recompute rate factor tables if clocks were changed.

		DELETE	standard actions
		PROCESS	standard actions

		SETUP	setup all internal tables.  Must be called before the
			hsolve element is used.  User accessible fields
			(path, chanmode, etc.) should have been set.  Model
			changes after the SETUP call may have no effect on
			the computations.

		DUPLICATE  does an efficient duplication of an hsolve element.
			Use this in a network simulation for identical
			copies of a neuron. Only the tables containing
			changing parameters (Vm, etc.) will be duplicated.
                        Syntax: call hsolve1 DUPLICATE hsolve2 duplic_path

		HPUT	updates chip-array.  Fields from a single element are
			put into the chip-array (chanmodes 2 or 3).  Do this
			call after a setfield command on the element.
			Syntax: call hsolve HPUT element_path

		HGET	updates hsolve-computed element.  Computed values are
                        put from the chip-array into a single element
                        (chanmodes 2-5).  Do this call before a getfield on
                        the element.  Syntax: call hsolve HGET element_path

		HRESTORE the complete chip-array is updated.  Fields from all
			hsolve-computed elements are put into chip-array
			(chanmodes 2-5).  Do this call after a restore
			command or after multiple setfield commands.
			Syntax: call hsolve HRESTORE


		HSAVE	all hsolve-computed elements are updated.  Computed
			values are put from the chip array into all elements
			(chanmodes 2-5).  Do this call before a save
			command or before multiple getfield commands.
			Syntax: call hsolve HSAVE

Messages:	None.	

Example:

	// Do all the preparatory grunge work
	.
	.
	// create the cell as an hsolve element
	readcell test.p /test -hsolve

        // set the chanmode, comptmode and calcmode (if other than defaults)
        setfield /test chanmode 2

	// set up the arrays and tables for the solver
	call /test/solve SETUP
	// use the Crank-Nicholson method for the hsolve element
	setmethod 11
	// It is essential to call reset (or reschedule) after setting
	// up an hsolver, so that the process list gets updated.
	reset

------------------------------------------------------------------------------

Notes:

The hsolve element completely takes over the calculations for compartment and
symcompartment elements and certain other elements specified as part of the
cell. Each hsolve element should only solve one cell. The timestep for the
integration is determined by the clock assigned to the hsolve element, and the
clocks for these elements which are taken over are ignored.

When using the Hines solver with a neuron, it is best to think of the entire
neuron as a single object since the individual compartments within the cell
are no longer responsible for their own computations.  This loss of object
orientedness is mitigated to some degree by the ability of the solver to
transparently interact with elements utilizing other integration schemes.

There are 2 basic modes for hsolve operation: without chip-array (chanmode 0
or 1), or with chip-array (chanmodes 2-5).  Without chip-array is the most
compatible mode, but is the slowest.  With the chip-array, hsolve is much
faster because the original elements are no longer used, instead all
simulation parameters are stored in a huge array (this improves memory access
times).  Unfortunately you can no longer expect that it will update the fields
in computed objects (like Vm in a compartment or Ca in Ca_concen) so that
graphic or file output might not work.  Vice versa, if you change a field
(like inject in compartment) it might not affect the simulation.  However,
there are methods available to get values in and out of the chip-array, which
involve special settings and/or the use of the HPUT, HGET, HSAVE and HRESTORE
actions listed above.  This is ilustrated in another example below.  The
documentation for findsolvefield describes another method, introduced in
GENESIS 2.1, that also allows access to these values.

The example above illustrates the process of setting up the hsolve element.
First, one should create the cell as an hsolve element at the root of the cell
element tree.  This may be done in one step by using the "-hsolve" option with
readcell.  Note that, starting with GENESIS version 2.2, it is necessary to
provide the full path to the cell when using readcell to directly create an
hsolve.  (e.g.  you can't use "readcell test.p test -hsolve", even if "/" is
the current working element.)

If readcell is not used, then the hsolved cell should be created with
statements like

    create hsolve /test
    create compartment /test/soma
    create tabchannel /test/soma/Na_channel
    ...

It is required that only compartments (or symcompartments) be children of the
cell, and that channel, concentration, etc. elements should be children or
grandchildren of the compartment to which they are attached.

Next, set any fields of the hsolve which are needed to specify non-default
values of chanmode, comptmode, calcmode, storemode, etc.  Then, one needs to
call the SETUP action in order to tell the hsolve element to create all the
solution arrays and tables.  Finally, chose either method 10 (backwards Euler)
or 11 (Crank-Nicholson) as the method to be used.  All elements that lie
outside this tree will continue to be treated by the previous explicit method
(typically exponential Euler).  The Scripts/examples/hines directory contains
an example that which demonstrates both chanmodes 0 and 2.

Although it is not presently necessary to create the hsolve element as the
root of the element tree, future versions of hsolve may require this.  The
older method of creating the hsolve is to create the cell as a neutral
element, create the hsolve element as a child element, and then set the hsolve
path field to indicate which elements will be taken over.  Typically, a
wildcard path is used to refer to all compartment elements in the simulation.
This includes all hsolvable sublelements of these compartments.  For example,

        // readcell will place the compartments below the neutral '/test'
	readcell test.p /test

	// create the hines solver element below the cell '/test'
	create hsolve /test/solve

	// Specify the path for the solver
	setfield /test/solve path /test/##[][TYPE=compartment]


RESTRICTIONS

It is important to be aware of some of the restrictions imposed by the use
of hsolve:

Only the objects compartment, symcompartment, tabchannel, tab2Dchannel,
tabcurrent, spikegen, Ca_concen, concpool, nernst, Mg_block, ghk, taupump,
mmpump, hillpump, difshell, fixbuffer, difbuffer, dif2buffer, fura2, synchan,
and the oldconn library channels channelC2 and channelC3 are handled by
hsolve.  If your simulation uses these listed objects, the use of hsolve will
increase speed significantly.  Note that as of GENESIS version 2.2,
symcompartment objects are handled by hsolve.

An important restriction introduced in GENESIS version 2.1 is that, for
chanmodes 2-5, the element tree of your cell (or other element tree to be
taken over by hsolve) must not contain any non-hsolvable elements other than
neutral elements.  For existing simulations which violate this restriction, it
will be easiest to switch to chanmode 0 or 1.

It is not guaranteed that messages from unlisted objects to listed objects
will work, though solve should alert you about any such incompatibilities.
Most incoming messages will work if you use chanmodes 2-5, however you will
receive a warning if the source of the message is not hsolved.  This is to let
you know that the source element is being solved by a less accurate explicit
method.  Messages from listed to unlisted objects will only work for certain
chanmode settings.

hsolve does not keep track of changes made to the model after the hsolve SETUP
command (particularly the following commands: create, copy, delete, addmsg,
deletemsg).  In chanmode 0 or 1 this only affects the listed objects and
messages between them; in chanmodes 2-5 this affects the listed objects and
ALL messages going in or out from them.  The hsolver should be used in these
modes only for finished models (you can use it for parameter searches); do not
use it if you are still constructing and testing a model.

hsolve may miscalculate if you change clocks used by tabchannels,
tab2Dchannels, synchans, or channelC2/C3 objects without doing a reset (only
in chanmodes 2-5).

The Hines solver utilizes a considerable amount of memory.  Memory use
increases as the result of comptmode + chanmode - no_elminfo flags.  In
particular, going from chanmode 0 or 1 to chanmode 2, 3, or 4 causes a big
jump in memory usage.  Approximately 100 bytes are required per compartment if
the speed-optimized version of the hsolve element is used. Most of this
storage is for tables describing the sequence of calculations. These tables
can, however, be shared between cells which are identical in topology and
differ only in parameter values.

FURTHER EXAMPLES

* Copying cells with the hsolve element

When a cell using the hsolve element is copied, the pointers in the
duplicated hsolve element are unchanged. In other words, they refer back to
the original cell. One could simply change the path to refer to the current
cell and call SETUP again, but that would unnecessarily duplicate a lot of
tables.  The efficient option would be to use the original tables where they
are identical, and construct new ones where they refer to the current cell.
This is done using the DUPLICATE sction, which takes the name of the new
hsolve element and a wildcard path duplic_path as arguments.  duplic_path
points to all compartments that should be taken over by the new solver.
For example,

	// Copy the original cell to /test2
	copy /test /test2

        // create the hines solver element below the cell '/test'
        create hsolve /test/solve

        // Specify the path for the solver and the chanmode
        setfield /test/solve path /test/##[][TYPE=compartment] chanmode 4

        call /test/solve SETUP

        // Duplicate the hsolver
        call /test/solve DUPLICATE /test2/solve /test2/##[][TYPE=compartment]

        reset

There are several commands, including createmap and cellsheet, that copy
cells. The DUPLICATE command should be issued for the hsolvers on each copy.
When using chanmodes 3 or higher, findsolvefield must be used with messages in
order to access fields of duplicated hsolves, as shown in the documentation
for findsolvefield.

* Deleting hsolve elements

When an hsolve element is set up it removes the relevant channels and
compartments from the list of elements to be processed, by setting a bit
(0x100) on the flag field of those elements. These elements are re-enabled
and the bit set back to 0 when the hsolve element is deleted.

WARNING: in the present implementation, the hsolve element does not check
whether it has had duplicates made using the DUPLICATE command. This makes it
possible to delete the tables for all of the copies by deleting any one of
them, resulting in segmentation violations and similar amusements. In other
words, do not delete hsolves on duplicated cells.

* Accessing data fields in chanmodes 2-5.

A script similar to this example might be used for running a long simulation
as a background job, using no graphics.  The HPUT and HSAVE actions are used
in chanmode 2 or 3 to allow the simulation results and the final state of
the simulation to be saved to disk.  A disk_out element is used to output the
Vm of a soma compartment in a compressed binary format.  However, it is
equally possible to output Vm from every compartment, or from an array of
cells.  Later, a script using a disk_in element can display the results using
messages from the disk_in to an xgraph or xview.  The save command is used
here to save field values of all elements except nernst elements and the
hsolve element.  This will allow the simulation to be resumed from its final
state by using the restore command.

For more details, see Scripts/examples/XODUS/fileview and the documentation
for disk_out, disk_in, save, and restore.

	/* load scripts and global variables */
	include defaults
	include other_stuff

	pushe /library
	make_my_stuff
	pope

	/* make the model */
	readcell my_cell {my_cellpath} -hsolve

	/* set the clocks */
	setclock 0 1e-5	  /* integration */
	setclock 1 1e-4   /* output */

	/* create the output or graphics elements */

	/* note that if you want to output Vm only, but from multiple
	** compartments, you should use chanmode == 2; if you want to
	**  output other fields or Vm from only one compartment, 
	** use chanmode == 3 */

	create disk_out /output/disk
	useclock /output/disk 1
	addmsg {my_cellpath}/soma /output/disk Vm
	addmsg ...

	/* create any other elements */
	....

	/* setup the hines solver */
	setfield {my_cellpath} \
		comptmode   1 \
		chanmode    3 \
		outclock    1
	call {my_cellpath} SETUP
	setmethod 11		  // Crank-Nicholson integration method

	/* initialize output */
	setfield /output/plot_out filename {filename} initialize 1
	setfield /output/plot_out filename {filename} append 1 leave_open 1

	reset

	step 0.10 -t
	/* do a current injection */
	setfield {my_cellpath}/soma inject 1.0e-9         /* in Amps */
	/* update chip array */
	call {my_cellpath} HPUT {my_cellpath}/soma 
	step 2.00 -t

	/* save the integration values to disk */
	call {my_cellpath} HSAVE
	save {my_cellpath}/##[][TYPE!=nernst][TYPE!=hsolve] {savename} 

	quit

---------------------------------------------------------------------------
DETAILED DESCRIPTION

The Hines library provides elements and functions for the efficient implicit
solution of the systems of differential equations arising in single-neuron
models. The sparse matrix arising from the the branched structure of neurons
is ordered by the method described by Michael Hines, which permits it to be
solved in order N operations using Gaussian elimination without pivoting.
The nonlinear equations resulting from the Hodgkin-Huxley description of ion
channels are treated as conditionally linear, and also solved in an
efficient, second-order manner.  Compatibility with other integration
schemes is maintained so that mixed integration schemes are feasible.

The element sets up a data structure whereby the only change needed for
other identical cells is to reassign pointers to the compartment and channel
elements.  This avoids having to reallocate and rederive the solving
scheme.  Solutions are done using gaussian forward and backward elimination
without pivoting. The numbering and evaluation sequence ensures that no new
off-diagonal terms are formed which might mess up the sparse matrix
representation.  The sparse matrix is represented by a single 'values' array
indexed using the ri and cip arrays.  The actual solution system is stored
as a linear sequence of operations on the 'values' array, in an attempt to
speed things up a bit.

The implicit method leads to numerically stable solutions for all time
steps.  Numerical accuracy, however, imposes limits on the integration time
step which must be empirically evaluated for each simulation.  Typical
values are 100 microsecond time steps for purely passive membranes and 20
microsecond time steps for membranes with active channels.  This confers a
considerable speed up in detailed neuronal simulations where the presence of
electrotonically small compartments (such as dendritic spines) leads to
systems of equations which are numerically very stiff. In order to obtain
similar accuracy using explicit integration methods, time steps in the
nanosecond range may be required.

In addition to the substantial speed advantage gained from the use of an
implicit solution, the coding of the Hines solver element has been done with
a view to computational efficiency.  This results in a speed up of around 50%
compared to the same model with the same time step, using the conventional
explicit Euler integration scheme.

Code changes needed for the Hines library

This section is probably of interest only to individuals who have developed
their own code libraries and who are deeply entangled in GENESIS coding
issues. For most users the changes are transparent, once the .simrc and
schedule.g file have been updated to reflect the addition of a new library.

The incorporation of a completely new integration scheme to the simulator can
be expected to have numerous repercussions. Fortunately, the changes needed for
the Hines library turned out to be relatively painless. No changes were needed
for the code of any of the existing elements. There is a slight rearrangement
of the order of fields in the various channel objects, and a new globally
#define'd type called CHAN_TYPE has been created so that all these channel
objects have the following common fields at identical memory offsets: Ik, Gk
and Ek (See src/sim/struct_defs.h). All channels whose structures obey this
convention can be used in cells utilizing the hsolve element.  A very small
number of channel objects (such as channelC) do not share these fields and
therefore cannot be part of cell models which use the hines solver.

A new flag (0x100) has been employed for the Hines solver, which causes the
removal of compartment and tabchannel elements from the action list (the list
of elements whose action functions are to be called according to the clocks).
This flag is similar to the disable option except that it leaves the children
of these elements intact.

In the interests of speed (which is, after all, what the Hines method is all
about) I have 'unrolled' the sparse matrix solution into a single giant
function array. This is done at setup time, so the actual solution does not
involve any conditionals and minimises array lookups. This is moderately
expensive in terms of memory. The size of the function array allocated is

     S = FA * sizeof (int) ;  FA = 10 + 1.5 * M^2 / N

where S is the size of the array, FA is the number of functions allocated, M
is the number of non-zero coefficients in the solution matrix and N is the
number of compartments. M depends strongly on the branching pattern of the
cell. FA is an approximation to, and is greater than F, which is the actual
number of locations needed in the function table.  For example, in a mitral
cell model, which has limited branching:

     N = 286 ; M = 856 ; FA = 3853 ; S = 15412 ; F = 3710

In a granule cell model with numerous dendritic spines:

     N = 944 ; M = 2830 ; FA = 12736 ; S = 50944 ; F = 12264

A reasonable approximation is that each compartment requires about 54 bytes in
the function table, which is not too bad.

References:	Hines, M. (1984) Efficient computation of branched nerve
		equations.  Int. J. Bio-Med. Comp.  15: 69-76

		Mascagni, M.V. (1989) Methods in Neuronal Modeling Ed: Koch
		and Segev. Chapter 13. 439-484

See also:	setmethod, findsolvefield
