HSolvePassive Class Reference

#include <HSolvePassive.h>

Inheritance diagram for HSolvePassive:

Collaboration diagram for HSolvePassive:

Public Member Functions
void	setup (Id seed, double dt)
void	solve ()
Public Member Functions inherited from HinesMatrix
double	getA (unsigned int row, unsigned int col) const
double	getB (unsigned int row) const
unsigned int	getSize () const
double	getVMid (unsigned int row) const
	HinesMatrix ()
void	setup (const vector< TreeNodeStruct > &tree, double dt)

Protected Member Functions
void	backwardSubstitute ()
void	forwardEliminate ()
void	updateMatrix ()

Protected Attributes
vector< CompartmentStruct >	compartment_
vector< Id >	compartmentId_
map< unsigned int, InjectStruct >	inject_
vector< TreeNodeStruct >	tree_
vector< double >	V_
Protected Attributes inherited from HinesMatrix
vector< vdIterator >	backOperand_
double	dt_
vector< double >	HJ_
vector< double >	HJCopy_
vector< double >	HS_
vector< JunctionStruct >	junction_
unsigned int	nCompt_
vector< vdIterator >	operand_
int	stage_
	reached. Used in getA. More...
vector< double >	VMid_
	middle of a time step. More...

Private Member Functions
void	clear ()
double	getV (unsigned int row) const
void	initialize ()
void	storeTree ()
void	walkTree (Id seed)

Additional Inherited Members
Protected Types inherited from HinesMatrix
typedef vector< double >::iterator	vdIterator

Detailed Description

Definition at line 20 of file HSolvePassive.h.

Member Function Documentation

void HSolvePassive::backwardSubstitute ( )

protected

Definition at line 325 of file HSolvePassive.cpp.

References HinesMatrix::backOperand_, HinesMatrix::HS_, HinesMatrix::junction_, HinesMatrix::nCompt_, HinesMatrix::operand_, HinesMatrix::stage_, V_, and HinesMatrix::VMid_.

Referenced by solve(), and HSolveActive::step().

{
    int ic = nCompt_ - 1;
    vector< double >::reverse_iterator ivmid = VMid_.rbegin();
    vector< double >::reverse_iterator iv = V_.rbegin();
    vector< double >::reverse_iterator ihs = HS_.rbegin();
    vector< vdIterator >::reverse_iterator iop = operand_.rbegin();
    vector< vdIterator >::reverse_iterator ibop = backOperand_.rbegin();
    vector< JunctionStruct >::reverse_iterator junction;

    *ivmid = *ihs / *( ihs + 3 );
    *iv = 2 * *ivmid - *iv;
    --ic, ++ivmid, ++iv, ihs += 4;

    int index;
    int rank;
    for ( junction = junction_.rbegin();
            junction != junction_.rend();
            junction++ )
    {
        index = junction->index;
        rank = junction->rank;

        while ( ic > index )
        {
            *ivmid = ( *ihs - *( ihs + 2 ) **( ivmid - 1 ) ) / *( ihs + 3 );
            *iv = 2 * *ivmid - *iv;

            --ic, ++ivmid, ++iv, ihs += 4;
        }

        if ( rank == 1 )
        {
            *ivmid = ( *ihs - **iop ***( iop + 2 ) ) / *( ihs + 3 );

            iop += 3;
        }
        else if ( rank == 2 )
        {
            vdIterator v0 = *( iop );
            vdIterator v1 = *( iop + 2 );
            vdIterator j  = *( iop + 4 );

            *ivmid = ( *ihs
                       - *v0 **( j + 2 )
                       - *v1 **j
                     ) / *( ihs + 3 );

            iop += 5;
        }
        else
        {
            *ivmid = *ihs;
            for ( int i = 0; i < rank; ++i )
            {
                *ivmid -= **ibop ***( ibop + 1 );
                ibop += 2;
            }
            *ivmid /= *( ihs + 3 );

            iop += 3 * rank * ( rank + 1 );
        }

        *iv = 2 * *ivmid - *iv;
        --ic, ++ivmid, ++iv, ihs += 4;
    }

    while ( ic >= 0 )
    {
        *ivmid = ( *ihs - *( ihs + 2 ) **( ivmid - 1 ) ) / *( ihs + 3 );
        *iv = 2 * *ivmid - *iv;

        --ic, ++ivmid, ++iv, ihs += 4;
    }

    stage_ = 2;    // Backward substitution done.
}

Here is the caller graph for this function:

void HSolvePassive::clear ( )

private

Definition at line 35 of file HSolvePassive.cpp.

References compartment_, compartmentId_, HinesMatrix::dt_, inject_, tree_, and V_.

Referenced by setup().

{
    dt_ = 0.0;
    compartment_.clear();
    compartmentId_.clear();
    V_.clear();
    tree_.clear();
    inject_.clear();
}

Here is the caller graph for this function:

void HSolvePassive::forwardEliminate ( )

protected

Definition at line 246 of file HSolvePassive.cpp.

References HinesMatrix::HS_, HinesMatrix::junction_, HinesMatrix::nCompt_, HinesMatrix::operand_, and HinesMatrix::stage_.

Referenced by solve(), and HSolveActive::step().

{
    unsigned int ic = 0;
    vector< double >::iterator ihs = HS_.begin();
    vector< vdIterator >::iterator iop = operand_.begin();
    vector< JunctionStruct >::iterator junction;

    double pivot;
    double division;
    unsigned int index;
    unsigned int rank;
    for ( junction = junction_.begin();
            junction != junction_.end();
            junction++ )
    {
        index = junction->index;
        rank = junction->rank;

        while ( ic < index )
        {
            *( ihs + 4 ) -= *( ihs + 1 ) / *ihs **( ihs + 1 );
            *( ihs + 7 ) -= *( ihs + 1 ) / *ihs **( ihs + 3 );

            ++ic, ihs += 4;
        }

        pivot = *ihs;
        if ( rank == 1 )
        {
            vdIterator j = *iop;
            vdIterator s = *( iop + 1 );

            division    = *( j + 1 ) / pivot;
            *( s )     -= division **j;
            *( s + 3 ) -= division **( ihs + 3 );

            iop += 3;
        }
        else if ( rank == 2 )
        {
            vdIterator j = *iop;
            vdIterator s;

            s           = *( iop + 1 );
            division    = *( j + 1 ) / pivot;
            *( s )     -= division **j;
            *( j + 4 ) -= division **( j + 2 );
            *( s + 3 ) -= division **( ihs + 3 );

            s           = *( iop + 3 );
            division    = *( j + 3 ) / pivot;
            *( j + 5 ) -= division **j;
            *( s )     -= division **( j + 2 );
            *( s + 3 ) -= division **( ihs + 3 );

            iop += 5;
        }
        else
        {
            vector< vdIterator >::iterator
            end = iop + 3 * rank * ( rank + 1 );
            for ( ; iop < end; iop += 3 )
                **iop -= **( iop + 2 ) / pivot ***( iop + 1 );
        }

        ++ic, ihs += 4;
    }

    while ( ic < nCompt_ - 1 )
    {
        *( ihs + 4 ) -= *( ihs + 1 ) / *ihs **( ihs + 1 );
        *( ihs + 7 ) -= *( ihs + 1 ) / *ihs **( ihs + 3 );

        ++ic, ihs += 4;
    }

    stage_ = 1;    // Forward elimination done.
}

Here is the caller graph for this function:

double HSolvePassive::getV ( unsigned int  row ) const

private

Definition at line 407 of file HSolvePassive.cpp.

References V_.

{
    return V_[ row ];
}

void HSolvePassive::initialize ( )

private

Definition at line 100 of file HSolvePassive.cpp.

References CompartmentStruct::CmByDt, compartment_, compartmentId_, HinesMatrix::dt_, CompartmentStruct::EmByRm, Field< A >::get(), inject_, HSolveUtils::leakageChannels(), HinesMatrix::nCompt_, and V_.

Referenced by setup().

{
    nCompt_ = compartmentId_.size();

    double Vm, Cm, Em, Rm, inject;
    double EmLeak, GmLeak;
    double EmGmThev, GmThev;

    vector< Id > leakage;
    vector< Id >::iterator ileakage;

    for ( unsigned int ic = 0; ic < compartmentId_.size(); ++ic )
    {
        Id cc = compartmentId_[ ic ];

        Vm = Field< double >::get( cc, "Vm" );
        Cm = Field< double >::get( cc, "Cm" );
        Em = Field< double >::get( cc, "Em" );
        Rm = Field< double >::get( cc, "Rm" );
        inject = Field< double >::get( cc, "inject" );
        V_.push_back( Vm );

        /*
         *  If there are 'n' leakage channels with reversal potentials:
         *       E[ 1 ] ... E[ n ]
         *  and resistances:
         *       R[ 1 ] ... R[ n ]
         *  respectively, then the membrane circuit can be represented by the
         *  Thevenin equivalent circuit with:
         *
         *  R_thevenin = 1 / sum( 1 / R[ i ] )
         *             = 1 / sum( G[ i ] )
         *  E_thevenin = sum( E[ i ] / R[ i ] ) * R_thevenin
         *             = sum( E[ i ] * G[ i ] ) * R_thevenin
         *
         *  Here we have to add original membrane potential Em as E[ 0 ] and
         *  resistance Rm as R[ 0 ].
         */
        GmThev   = 1.0 / Rm;
        EmGmThev = Em / Rm;

        HSolveUtils::leakageChannels( compartmentId_[ ic ], leakage );
        for ( ileakage = leakage.begin();
                ileakage != leakage.end();
                ileakage++ )
        {
            EmLeak = Field< double >::get( *ileakage, "Ek" );
            GmLeak = Field< double >::get( *ileakage, "Gk" );
            GmThev   += GmLeak;
            EmGmThev += EmLeak * GmLeak;
        }

//        RmThev = 1.0 / GmThev;

        CompartmentStruct compartment;
        compartment.CmByDt = 2.0 * Cm / dt_;
        compartment.EmByRm = EmGmThev;
        compartment_.push_back( compartment );

        if ( inject != 0.0 )
        {
            inject_[ ic ].injectVarying = 0.0;
            inject_[ ic ].injectBasal = inject;
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void HSolvePassive::setup	(	Id	seed,
		double	dt
	)

Definition at line 14 of file HSolvePassive.cpp.

References clear(), HinesMatrix::dt_, initialize(), HinesMatrix::setup(), storeTree(), tree_, and walkTree().

Referenced by HSolveActive::setup().

{
    clear();
    dt_ = dt;
    walkTree( seed );
    initialize();
    storeTree();
    HinesMatrix::setup( tree_, dt_ );
}

Here is the call graph for this function:

Here is the caller graph for this function:

void HSolvePassive::solve

(

)

Definition at line 24 of file HSolvePassive.cpp.

References backwardSubstitute(), forwardEliminate(), and updateMatrix().

{
    updateMatrix();
    forwardEliminate();
    backwardSubstitute();
}

Here is the call graph for this function:

void HSolvePassive::storeTree ( )

private

Definition at line 167 of file HSolvePassive.cpp.

References HSolveUtils::children(), TreeNodeStruct::children, TreeNodeStruct::Cm, compartmentId_, TreeNodeStruct::Em, Field< A >::get(), TreeNodeStruct::initVm, HinesMatrix::nCompt_, TreeNodeStruct::Ra, TreeNodeStruct::Rm, and tree_.

Referenced by setup().

{
    double Ra, Rm, Cm, Em, initVm;

    // Create a map from the MOOSE Id to Hines' index.
    map< Id, unsigned int > hinesIndex;
    for ( unsigned int ic = 0; ic < nCompt_; ++ic )
        hinesIndex[ compartmentId_[ ic ] ] = ic;

    vector< Id > childId;
    vector< Id >::iterator child;

    vector< Id >::iterator ic;

    for ( ic = compartmentId_.begin(); ic != compartmentId_.end(); ++ic )
    {
        childId.clear();

        HSolveUtils::children( *ic, childId );
        Ra = Field< double >::get( *ic, "Ra" );
        Cm = Field< double >::get( *ic, "Cm" );
        Rm = Field< double >::get( *ic, "Rm" );
        Em = Field< double >::get( *ic, "Em" );
        initVm = Field< double >::get( *ic, "initVm" );

        TreeNodeStruct node;
        // Push hines' indices of children
        for ( child = childId.begin(); child != childId.end(); ++child )
            node.children.push_back( hinesIndex[ *child ] );

        node.Ra = Ra;
        node.Rm = Rm;
        node.Cm = Cm;
        node.Em = Em;
        node.initVm = initVm;

        tree_.push_back( node );
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void HSolvePassive::updateMatrix ( )

protected

Definition at line 211 of file HSolvePassive.cpp.

References compartment_, HinesMatrix::HJ_, HinesMatrix::HJCopy_, HinesMatrix::HS_, inject_, InjectStruct::injectBasal, InjectStruct::injectVarying, HinesMatrix::stage_, V_, and value.

Referenced by solve().

{
    /*
     * Copy contents of HJCopy_ into HJ_. Cannot do a vector assign() because
     * iterators to HJ_ get invalidated in MS VC++
     */
    if ( HJ_.size() != 0 )
        memcpy( &HJ_[ 0 ], &HJCopy_[ 0 ], sizeof( double ) * HJ_.size() );

    vector< double >::iterator ihs = HS_.begin();
    vector< double >::iterator iv  = V_.begin();

    vector< CompartmentStruct >::iterator ic;
    for ( ic = compartment_.begin(); ic != compartment_.end(); ++ic )
    {
        *ihs         = *( 2 + ihs );
        *( 3 + ihs ) = *iv * ic->CmByDt + ic->EmByRm;

        ihs += 4, ++iv;
    }

    map< unsigned int, InjectStruct >::iterator inject;
    for ( inject = inject_.begin(); inject != inject_.end(); inject++ )
    {
        unsigned int ic = inject->first;
        InjectStruct& value = inject->second;

        HS_[ 4 * ic + 3 ] += value.injectVarying + value.injectBasal;

        value.injectVarying = 0.0;
    }

    stage_ = 0;    // Update done.
}

Here is the caller graph for this function:

void HSolvePassive::walkTree ( Id  seed )

private

Definition at line 45 of file HSolvePassive.cpp.

References HSolveUtils::adjacent(), and compartmentId_.

Referenced by setup().

{
    //~ // Dirty call to explicitly call the compartments reinitFunc.
    //~ // Should be removed eventually, and replaced with a cleaner way to
    //~ // initialize the model being read.
    //~ HSolveUtils::initialize( seed );

    // Find leaf node
    Id previous;
    vector< Id > adjacent;
    HSolveUtils::adjacent( seed, adjacent );
    if ( adjacent.size() > 1 )
        while ( !adjacent.empty() )
        {
            previous = seed;
            seed = adjacent[ 0 ];

            adjacent.clear();
            HSolveUtils::adjacent( seed, previous, adjacent );
        }

    // Depth-first search
    vector< vector< Id > > cstack;
    Id above;
    Id current;
    cstack.resize( 1 );
    cstack[ 0 ].push_back( seed );
    while ( !cstack.empty() )
    {
        vector< Id >& top = cstack.back();

        if ( top.empty() )
        {
            cstack.pop_back();
            if ( !cstack.empty() )
                cstack.back().pop_back();
        }
        else
        {
            if ( cstack.size() > 1 )
                above = cstack[ cstack.size() - 2 ].back();

            current = top.back();
            compartmentId_.push_back( current );

            cstack.resize( cstack.size() + 1 );
            HSolveUtils::adjacent( current, above, cstack.back() );
        }
    }

    // Compartments get ordered according to their hines' indices once this
    // list is reversed.
    reverse( compartmentId_.begin(), compartmentId_.end() );
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

vector< CompartmentStruct > HSolvePassive::compartment_

protected

Definition at line 36 of file HSolvePassive.h.

Referenced by clear(), HSolve::getIm(), initialize(), and updateMatrix().

vector< Id > HSolvePassive::compartmentId_

protected

Definition at line 37 of file HSolvePassive.h.

Referenced by clear(), initialize(), HSolveActive::manageOutgoingMessages(), HSolve::mapIds(), HSolveActive::readExternalChannels(), HSolveActive::readHHChannels(), HSolveActive::readSynapses(), storeTree(), HSolve::unzombify(), walkTree(), and HSolve::zombify().

map< unsigned int, InjectStruct > HSolvePassive::inject_

protected

inject map. contains the list of compartments that have current injections into them.

Definition at line 45 of file HSolvePassive.h.

Referenced by HSolve::addInject(), clear(), HSolve::getInject(), initialize(), HSolve::setInject(), and updateMatrix().

vector< TreeNodeStruct > HSolvePassive::tree_

protected

Tree info. The tree is used to acquire various values during setup. It contains the user-defined original values of all compartment parameters. Therefore, it is also used during reinit.

Definition at line 41 of file HSolvePassive.h.

Referenced by clear(), HSolve::getCm(), HSolve::getEm(), HSolve::getIm(), HSolve::getInitVm(), HSolve::getRa(), HSolve::getRm(), HSolveActive::reinitCompartments(), HSolve::setCm(), HSolve::setEm(), HSolve::setInitVm(), HSolve::setRa(), HSolve::setRm(), setup(), and storeTree().

vector< double > HSolvePassive::V_

protected

Compartment Vm. V_ is addressed using a compartment index. V_ stores the Vm value of each compartment.

Definition at line 38 of file HSolvePassive.h.

Referenced by backwardSubstitute(), clear(), HSolve::getIk(), HSolve::getIm(), getV(), HSolve::getVm(), initialize(), HSolveActive::readSynapses(), HSolveActive::reinitChannels(), HSolveActive::reinitCompartments(), HSolve::setVm(), and updateMatrix().

---

The documentation for this class was generated from the following files:

• moose-core/hsolve/HSolvePassive.h
• moose-core/hsolve/HSolvePassive.cpp