HinesMatrix Class Reference

#include <HinesMatrix.h>

Inheritance diagram for HinesMatrix:

Collaboration diagram for HinesMatrix:

Public Member Functions
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 Types
typedef vector< double >::iterator	vdIterator

Protected Attributes
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 ()
void	makeJunctions ()
void	makeMatrix ()
void	makeOperands ()
	elimination easier. More...

Private Attributes
vector< vector< unsigned int > >	coupled_
vector< double >	Ga_
map< unsigned int, unsigned int >	groupNumber_
map< unsigned int, vdIterator >	operandBase_
const vector< TreeNodeStruct > *	tree_
	setup. More...

Detailed Description

Definition at line 55 of file HinesMatrix.h.

Member Typedef Documentation

typedef vector< double >::iterator HinesMatrix::vdIterator

protected

Definition at line 68 of file HinesMatrix.h.

Constructor & Destructor Documentation

HinesMatrix::HinesMatrix

(

)

Definition at line 16 of file HinesMatrix.cpp.

    :
    nCompt_( 0 ),
    dt_( 0.0 ),
    stage_( -1 )
{
    ;
}

Member Function Documentation

void HinesMatrix::clear ( )

private

Definition at line 53 of file HinesMatrix.cpp.

References backOperand_, coupled_, dt_, Ga_, groupNumber_, HJ_, HJCopy_, HS_, junction_, nCompt_, operand_, operandBase_, stage_, tree_, and VMid_.

Referenced by setup().

{
    nCompt_ = 0;
    dt_ = 0.0;
    junction_.clear();
    HS_.clear();
    HJ_.clear();
    HJCopy_.clear();
    VMid_.clear();
    operand_.clear();
    backOperand_.clear();
    stage_ = 0;

    tree_ = 0;
    Ga_.clear();
    coupled_.clear();
    operandBase_.clear();
    groupNumber_.clear();
}

Here is the caller graph for this function:

double HinesMatrix::getA	(	unsigned int	row,
		unsigned int	col
	)		const

Definition at line 375 of file HinesMatrix.cpp.

References coupled_, group(), groupNumber_, HJ_, HS_, nCompt_, and stage_.

Referenced by operator<<().

{
    /*
     * If forward elimination is done, or backward substitution is done, and
     * if (row, col) is in the lower triangle, then return 0.
     */
    if ( ( stage_ == 1 || stage_ == 2 ) && row > col )
        return 0.0;

    if ( row >= nCompt_ || col >= nCompt_ )
        return 0.0;

    if ( row == col )
        return HS_[ 4 * row ];

    unsigned int smaller = row < col ? row : col;
    unsigned int bigger = row > col ? row : col;

    if ( groupNumber_.find( smaller ) == groupNumber_.end() )
    {
        if ( bigger - smaller == 1 )
            return HS_[ 4 * smaller + 1 ];
        else
            return 0.0;
    }
    else
    {
        // We could use: groupNumber = groupNumber_[ smaller ], but this is a
        // const function
        unsigned int groupNumber = groupNumber_.find( smaller )->second;
        const vector< unsigned int >& group = coupled_[ groupNumber ];
        unsigned int location, size;
        unsigned int smallRank, bigRank;

        if ( find( group.begin(), group.end(), bigger ) != group.end() )
        {
            location = 0;
            for ( int i = 0; i < static_cast< int >( groupNumber ); ++i )
            {
                size = coupled_[ i ].size();
                location += size * ( size - 1 );
            }

            size = group.size();
            smallRank = group.end() - find( group.begin(), group.end(), smaller ) - 1;
            bigRank = group.end() - find( group.begin(), group.end(), bigger ) - 1;
            location += size * ( size - 1 ) - smallRank * ( smallRank + 1 );
            location += 2 * ( smallRank - bigRank - 1 );

            if ( row == smaller )
                return HJ_[ location ];
            else
                return HJ_[ location + 1 ];
        }
        else
        {
            return 0.0;
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

double HinesMatrix::getB

(

unsigned int

row

)

const

Definition at line 436 of file HinesMatrix.cpp.

References HS_.

Referenced by operator<<().

{
    return HS_[ 4 * row + 3 ];
}

Here is the caller graph for this function:

unsigned int HinesMatrix::getSize

(

)

const

Definition at line 370 of file HinesMatrix.cpp.

References nCompt_.

Referenced by operator<<().

{
    return nCompt_;
}

Here is the caller graph for this function:

double HinesMatrix::getVMid

(

unsigned int

row

)

const

Definition at line 441 of file HinesMatrix.cpp.

References VMid_.

Referenced by operator<<().

{
    return VMid_[ row ];
}

Here is the caller graph for this function:

void HinesMatrix::makeJunctions ( )

private

This function creates junction structs to be stored in junction_. It does so by first looking through all compartments and finding those which have more than zero or one child. The zero-one child compts are left alone. coupled_ is populated with the rest of the compartments. Each element of coupled_ is sorted (they are all unsigned ints) and coupled_ itself is sorted by the first element in each element (groupCompare does this comparison in HinesMatrix.cpp). Note: the children themselves are unsigned ints that store the Hines index of the corresponding child compartment. So essentially, at each branch, a JunctionStruct is created for each child, which contains the hines index of that child and its rank, which is group-size() - childIndex - 1.

Definition at line 84 of file HinesMatrix.cpp.

References coupled_, group(), groupCompare(), groupNumber_, junction_, and nCompt_.

Referenced by setup().

{
    // 3.1
    for ( unsigned int i = 0; i < nCompt_; ++i )
    {
        const vector< unsigned int >& c = ( *tree_ )[ i ].children;

        if ( c.size() == 0 )
            continue;

        if ( c.size() == 1 )
        {
            int diff = ( int )( c[ 0 ] ) - i;

            if ( diff == 1 || diff == -1 )
                continue;
        }

        // "coupled" contains a list of all children..
        coupled_.push_back( c );
        // ..and the parent compartment itself.
        coupled_.back().push_back( i );
    }

    // 3.2
    vector< vector< unsigned int > >::iterator group;
    for ( group = coupled_.begin(); group != coupled_.end(); ++group )
        sort( group->begin(), group->end() );

    sort( coupled_.begin(), coupled_.end(), groupCompare );

    // 3.3
    unsigned int index;
    unsigned int rank;
    for ( group = coupled_.begin(); group != coupled_.end(); ++group )
        // Loop uptil penultimate compartment in group
        for ( unsigned int c = 0; c < group->size() - 1; ++c )
        {
            index = ( *group )[ c ];
            rank = group->size() - c - 1;
            junction_.push_back( JunctionStruct( index, rank ) );

            groupNumber_[ index ] = group - coupled_.begin();
        }

    sort( junction_.begin(), junction_.end() );
}

Here is the call graph for this function:

Here is the caller graph for this function:

void HinesMatrix::makeMatrix ( )

private

Populates HS_ and HJ_. All of the electrical circuit analysis goes into this one single function (and updateMatrix, of course).

Definition at line 133 of file HinesMatrix.cpp.

References coupled_, dt_, Ga_, group(), HJ_, HJCopy_, HS_, junction_, nCompt_, operandBase_, and tree_.

Referenced by setup().

{
    const vector< TreeNodeStruct >& node = *tree_;

    // Setting up HS
    HS_.resize( 4 * nCompt_, 0.0 );
    for ( unsigned int i = 0; i < nCompt_; ++i )
        HS_[ 4 * i + 2 ] =
            node[ i ].Cm / ( dt_ / 2.0 ) +
            1.0 / node[ i ].Rm;

    double gi, gj, gij;
    vector< JunctionStruct >::iterator junction = junction_.begin();
    for ( unsigned int i = 0; i < nCompt_ - 1; ++i )
    {
        if ( !junction_.empty() &&
                junction < junction_.end() &&
                i == junction->index )
        {
            ++junction;
            continue;
        }

        gi = Ga_[ i ];
        gj = Ga_[ i + 1 ];
        gij = gi * gj / ( gi + gj );

        HS_[ 4 * i + 1 ] = -gij;
        HS_[ 4 * i + 2 ] += gij;
        HS_[ 4 * i + 6 ] += gij;
    }

    vector< vector< unsigned int > >::iterator group;
    vector< unsigned int >::iterator i;
    for ( group = coupled_.begin(); group != coupled_.end(); ++group )
    {
        double gsum = 0.0;

        for ( i = group->begin(); i != group->end(); ++i )
            gsum += Ga_[ *i ];

        for ( i = group->begin(); i != group->end(); ++i )
        {
            gi = Ga_[ *i ];

            HS_[ 4 * *i + 2 ] += gi * ( 1.0 - gi / gsum );
        }
    }

    // Setting up HJ
    vector< unsigned int >::iterator j;
    unsigned int size = 0;
    unsigned int rank;
    for ( group = coupled_.begin(); group != coupled_.end(); ++group )
    {
        rank = group->size() - 1;
        size += rank * ( rank + 1 );
    }

    HJ_.reserve( size );

    for ( group = coupled_.begin(); group != coupled_.end(); ++group )
    {
        double gsum = 0.0;

        for ( i = group->begin(); i != group->end(); ++i )
            gsum += Ga_[ *i ];

        for ( i = group->begin(); i != group->end() - 1; ++i )
        {
            int base = HJ_.size();

            for ( j = i + 1; j != group->end(); ++j )
            {
                gij = Ga_[ *i ] * Ga_[ *j ] / gsum;

                HJ_.push_back( -gij );
                HJ_.push_back( -gij );
            }

            //~ operandBase_[ *i ] = &HJ_[ base ];
            operandBase_[ *i ] = HJ_.begin() + base;
        }
    }

    // Copy diagonal elements into their final locations
    for ( unsigned int i = 0; i < nCompt_; ++i )
        HS_[ 4 * i ] = HS_[ 4 * i + 2 ];
    // Create copy of HJ
    HJCopy_.assign( HJ_.begin(), HJ_.end() );
}

Here is the call graph for this function:

Here is the caller graph for this function:

void HinesMatrix::makeOperands ( )

private

elimination easier.

Makes operands in order to make forward

Definition at line 226 of file HinesMatrix.cpp.

References backOperand_, coupled_, group(), groupNumber_, HS_, junction_, nCompt_, operand_, operandBase_, and VMid_.

Referenced by setup().

{
    unsigned int index;
    unsigned int rank;
    unsigned int farIndex;
    vdIterator base;
    vector< JunctionStruct >::iterator junction;

    // Allocate space in VMid. Needed, since we will store pointers to its
    // elements below.
    VMid_.resize( nCompt_ );

    // Operands for forward-elimination
    for ( junction = junction_.begin(); junction != junction_.end(); ++junction )
    {
        index = junction->index;
        rank = junction->rank;
        base = operandBase_[ index ];

        // This is the list of compartments connected at a junction.
        const vector< unsigned int >& group =
            coupled_[ groupNumber_[ index ] ];

        if ( rank == 1 )
        {
            operand_.push_back( base );

            // Select last member.
            farIndex = group[ group.size() - 1 ];
            operand_.push_back( HS_.begin() + 4 * farIndex );
            operand_.push_back( VMid_.begin() + farIndex );
        }
        else if ( rank == 2 )
        {
            operand_.push_back( base );

            // Select 2nd last member.
            farIndex = group[ group.size() - 2 ];
            operand_.push_back( HS_.begin() + 4 * farIndex );
            operand_.push_back( VMid_.begin() + farIndex );

            // Select last member.
            farIndex = group[ group.size() - 1 ];
            operand_.push_back( HS_.begin() + 4 * farIndex );
            operand_.push_back( VMid_.begin() + farIndex );
        }
        else
        {
            // Operations on diagonal elements and elements from B (as in Ax = B).
            int start = group.size() - rank;
            for ( unsigned int j = 0; j < rank; ++j )
            {
                farIndex = group[ start + j ];

                // Diagonal elements
                operand_.push_back( HS_.begin() + 4 * farIndex );
                operand_.push_back( base + 2 * j );
                operand_.push_back( base + 2 * j + 1 );

                // Elements from B
                operand_.push_back( HS_.begin() + 4 * farIndex + 3 );
                operand_.push_back( HS_.begin() + 4 * index + 3 );
                operand_.push_back( base + 2 * j + 1 );
            }

            // Operations on off-diagonal elements.
            vdIterator left;
            vdIterator above;
            vdIterator target;

            // Upper triangle elements
            left = base + 1;
            target = base + 2 * rank;
            for ( unsigned int i = 1; i < rank; ++i )
            {
                above = base + 2 * i;
                for ( unsigned int j = 0; j < rank - i; ++j )
                {
                    operand_.push_back( target );
                    operand_.push_back( above );
                    operand_.push_back( left );

                    above += 2;
                    target += 2;
                }
                left += 2;
            }

            // Lower triangle elements
            target = base + 2 * rank + 1;
            above = base;
            for ( unsigned int i = 1; i < rank; ++i )
            {
                left = base + 2 * i + 1;
                for ( unsigned int j = 0; j < rank - i; ++j )
                {
                    operand_.push_back( target );
                    operand_.push_back( above );
                    operand_.push_back( left );

                    /*
                     * This check required because the MS VC++ compiler is
                     * paranoid about iterators going out of bounds, even if
                     * they are never used after that.
                     */
                    if ( i == rank - 1 && j == rank - i - 1 )
                        continue;

                    target += 2;
                    left += 2;
                }
                above += 2;
            }
        }
    }

    // Operands for backward substitution
    for ( junction = junction_.begin(); junction != junction_.end(); ++junction )
    {
        if ( junction->rank < 3 )
            continue;

        index = junction->index;
        rank = junction->rank;
        base = operandBase_[ index ];

        // This is the list of compartments connected at a junction.
        const vector< unsigned int >& group =
            coupled_[ groupNumber_[ index ] ];

        unsigned int start = group.size() - rank;
        for ( unsigned int j = 0; j < rank; ++j )
        {
            farIndex = group[ start + j ];

            backOperand_.push_back( base + 2 * j );
            backOperand_.push_back( VMid_.begin() + farIndex );
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void HinesMatrix::setup	(	const vector< TreeNodeStruct > &	tree,
		double	dt
	)

Definition at line 25 of file HinesMatrix.cpp.

References clear(), dt_, Ga_, makeJunctions(), makeMatrix(), makeOperands(), nCompt_, and tree_.

Referenced by HSolvePassive::setup().

{
    clear();

    nCompt_ = tree.size();

#if  SANITY_CHECK
    stringstream ss;
    if(nCompt_ <= 0)
    {
        ss << "Horror, horror! Trying to create a matrix with size " << nCompt_
           << endl;
        dump(ss.str(), "ERROR");
        throw range_error("Expected greater than 0.");
    }
#endif     /* -----  not STRICT_CHECK  ----- */

    dt_ = dt;
    tree_ = &tree;

    for ( unsigned int i = 0; i < nCompt_; i++ )
        Ga_.push_back( 2.0 / tree[ i ].Ra );

    makeJunctions();
    makeMatrix();
    makeOperands();
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

vector< vdIterator > HinesMatrix::backOperand_

protected

Definition at line 86 of file HinesMatrix.h.

Referenced by HSolvePassive::backwardSubstitute(), clear(), and makeOperands().

vector< vector< unsigned int > > HinesMatrix::coupled_

private

Contains a list of all children of a given compt. Also contains the parent itself. i.e., for each compartment that has more than one child, coupled_ stores a vector containing the children of the compartment and the compartment itself. coupled_ is therefore a vector of such vectors.

Definition at line 115 of file HinesMatrix.h.

Referenced by clear(), getA(), makeJunctions(), makeMatrix(), and makeOperands().

double HinesMatrix::dt_

protected

Definition at line 71 of file HinesMatrix.h.

Referenced by HSolvePassive::clear(), clear(), HSolvePassive::initialize(), makeMatrix(), HSolveActive::readCalcium(), HSolvePassive::setup(), and setup().

vector< double > HinesMatrix::Ga_

private

Definition at line 114 of file HinesMatrix.h.

Referenced by clear(), makeMatrix(), and setup().

map< unsigned int, unsigned int > HinesMatrix::groupNumber_

private

Tells you the index of a compartment's group within coupled_, given the compartment's Hines index.

Definition at line 124 of file HinesMatrix.h.

Referenced by clear(), getA(), makeJunctions(), and makeOperands().

vector< double > HinesMatrix::HJ_

protected

Hines, junctions. Flattened array containing the off-diagonal elements of the Hines matrix

Definition at line 79 of file HinesMatrix.h.

Referenced by clear(), getA(), makeMatrix(), and HSolvePassive::updateMatrix().

vector< double > HinesMatrix::HJCopy_

protected

Definition at line 82 of file HinesMatrix.h.

Referenced by clear(), makeMatrix(), and HSolvePassive::updateMatrix().

vector< double > HinesMatrix::HS_

protected

Hines, series. Flattened array containing the tridiagonal of the approximately tridiagonal Hines matrix, stacked against the column vector "b" that appears on the RHS of the equation that we're trying to solve: Ax=b.

Definition at line 74 of file HinesMatrix.h.

Referenced by HSolvePassive::backwardSubstitute(), clear(), HSolvePassive::forwardEliminate(), getA(), getB(), makeMatrix(), makeOperands(), and HSolvePassive::updateMatrix().

vector< JunctionStruct > HinesMatrix::junction_

protected

Definition at line 73 of file HinesMatrix.h.

Referenced by HSolvePassive::backwardSubstitute(), clear(), HSolvePassive::forwardEliminate(), makeJunctions(), makeMatrix(), and makeOperands().

unsigned int HinesMatrix::nCompt_

protected

Definition at line 70 of file HinesMatrix.h.

Referenced by HSolve::addInject(), HSolvePassive::backwardSubstitute(), clear(), HSolvePassive::forwardEliminate(), getA(), HSolve::getIm(), HSolve::getInject(), getSize(), HSolvePassive::initialize(), makeJunctions(), makeMatrix(), makeOperands(), HSolveActive::readCalcium(), HSolveActive::readHHChannels(), HSolveActive::readSynapses(), HSolveActive::reinitCompartments(), HSolve::setInject(), setup(), and HSolvePassive::storeTree().

vector< vdIterator > HinesMatrix::operand_

protected

Definition at line 85 of file HinesMatrix.h.

Referenced by HSolvePassive::backwardSubstitute(), clear(), HSolvePassive::forwardEliminate(), and makeOperands().

map< unsigned int, vdIterator > HinesMatrix::operandBase_

private

Contains iterators into HJ_ demarcating where a child's neighbours begin. Used for iterating through HJ_ along with junction_.

Definition at line 121 of file HinesMatrix.h.

Referenced by clear(), makeMatrix(), and makeOperands().

int HinesMatrix::stage_

protected

reached. Used in getA.

Which stage the simulation has

Definition at line 87 of file HinesMatrix.h.

Referenced by HSolvePassive::backwardSubstitute(), clear(), HSolvePassive::forwardEliminate(), getA(), and HSolvePassive::updateMatrix().

const vector< TreeNodeStruct >* HinesMatrix::tree_

private

setup.

Stores compt info for

Definition at line 112 of file HinesMatrix.h.

Referenced by clear(), makeMatrix(), and setup().

vector< double > HinesMatrix::VMid_

protected

middle of a time step.

Compartment voltage at the

Definition at line 83 of file HinesMatrix.h.

Referenced by HSolvePassive::backwardSubstitute(), clear(), getVMid(), and makeOperands().

---

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

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