Select Git revision
HinesMatrix.cpp 19.58 KiB
/**********************************************************************
** This program is part of 'MOOSE', the
** Messaging Object Oriented Simulation Environment.
** copyright (C) 2003-2007 Upinder S. Bhalla, Niraj Dudani and NCBS
** It is made available under the terms of the
** GNU Lesser General Public License version 2.1
** See the file COPYING.LIB for the full notice.
**********************************************************************/
#include "header.h"
#include "HinesMatrix.h"
#include <sstream>
#include <iomanip>
#include <stdexcept>
HinesMatrix::HinesMatrix()
:
nCompt_( 0 ),
dt_( 0.0 ),
stage_( -1 )
{
;
}
void HinesMatrix::setup( const vector< TreeNodeStruct >& tree, double dt )
{
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();
}
void HinesMatrix::clear()
{
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();}
bool groupCompare(
const vector< unsigned int >& A,
const vector< unsigned int >& B )
{
if ( A.empty() || B.empty() )
return 0;
return A[ 0 ] < B[ 0 ];
}
// Stage 3
void HinesMatrix::makeJunctions()
{
// 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() );
}
// Stage 4
void HinesMatrix::makeMatrix()
{
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() );
}
// Stage 5
void HinesMatrix::makeOperands()
{
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 );
}
}
}
///////////////////////////////////////////////////////////////////////////
// Public interface to matrix
///////////////////////////////////////////////////////////////////////////
unsigned int HinesMatrix::getSize() const
{
return nCompt_;
}
double HinesMatrix::getA( unsigned int row, unsigned int col ) const
{
/*
* 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;
}
}
}
double HinesMatrix::getB( unsigned int row ) const
{
return HS_[ 4 * row + 3 ];
}
double HinesMatrix::getVMid( unsigned int row ) const
{
return VMid_[ row ];
}
///////////////////////////////////////////////////////////////////////////
// Inserting into a stream
///////////////////////////////////////////////////////////////////////////
ostream& operator <<( ostream& s, const HinesMatrix& m )
{
unsigned int size = m.getSize();
s << "\nA:\n";
for ( unsigned int i = 0; i < size; i++ )
{
for ( unsigned int j = 0; j < size; j++ )
s << setw( 12 ) << setprecision( 5 ) << m.getA( i, j );
s << "\n";
}
s << "\n" << "V:\n";
for ( unsigned int i = 0; i < size; i++ )
s << m.getVMid( i ) << "\n";
s << "\n" << "B:\n";
for ( unsigned int i = 0; i < size; i++ )
s << m.getB( i ) << "\n";
return s;
}
///////////////////////////////////////////////////////////////////////////
#ifdef DO_UNIT_TESTS
#include "TestHSolve.h"
void testHinesMatrix()
{
vector< int* > childArray;
vector< unsigned int > childArraySize;
/**
* We test if the Hines' matrix is correctly setup for the following cell:
*
* Soma---> 15 - 14 - 13 - 12
* | |
* | L 11 - 10
* |
* L 16 - 17 - 18 - 19 * |
* L 9 - 8 - 7 - 6 - 5
* | |
* | L 4 - 3
* |
* L 2 - 1 - 0
*
* The numbers are the hines indices of compartments. Compartment X is the
* child of compartment Y if X is one level further away from the soma (#15)
* than Y. So #17 is the parent of #'s 2, 9 and 18.
*/
int childArray_1[ ] =
{
/* c0 */ -1,
/* c1 */ -1, 0,
/* c2 */ -1, 1,
/* c3 */ -1,
/* c4 */ -1, 3,
/* c5 */ -1,
/* c6 */ -1, 5,
/* c7 */ -1, 4, 6,
/* c8 */ -1, 7,
/* c9 */ -1, 8,
/* c10 */ -1,
/* c11 */ -1, 10,
/* c12 */ -1,
/* c13 */ -1, 12,
/* c14 */ -1, 11, 13,
/* c15 */ -1, 14, 16,
/* c16 */ -1, 17,
/* c17 */ -1, 2, 9, 18,
/* c18 */ -1, 19,
/* c19 */ -1,
};
childArray.push_back( childArray_1 );
childArraySize.push_back( sizeof( childArray_1 ) / sizeof( int ) );
/**
* Cell 2:
*
* 3
* |
* Soma---> 2
* / \
* / \
* 1 0
*
*/
int childArray_2[ ] =
{
/* c0 */ -1,
/* c1 */ -1,
/* c2 */ -1, 0, 1, 3,
/* c3 */ -1,
};
childArray.push_back( childArray_2 );
childArraySize.push_back( sizeof( childArray_2 ) / sizeof( int ) );
/**
* Cell 3:
*
* 3
* |
* 2
* / \
* / \ * 1 0 <--- Soma
*
*/
int childArray_3[ ] =
{
/* c0 */ -1, 2,
/* c1 */ -1,
/* c2 */ -1, 1, 3,
/* c3 */ -1,
};
childArray.push_back( childArray_3 );
childArraySize.push_back( sizeof( childArray_3 ) / sizeof( int ) );
/**
* Cell 4:
*
* 3 <--- Soma
* |
* 2
* / \
* / \
* 1 0
*
*/
int childArray_4[ ] =
{
/* c0 */ -1,
/* c1 */ -1,
/* c2 */ -1, 0, 1,
/* c3 */ -1, 2,
};
childArray.push_back( childArray_4 );
childArraySize.push_back( sizeof( childArray_4 ) / sizeof( int ) );
/**
* Cell 5:
*
* 1 <--- Soma
* |
* 2
* / \
* 4 0
* / \
* 3 5
*
*/
int childArray_5[ ] =
{
/* c0 */ -1,
/* c1 */ -1, 2,
/* c2 */ -1, 0, 4,
/* c3 */ -1,
/* c4 */ -1, 3, 5,
/* c5 */ -1,
};
childArray.push_back( childArray_5 );
childArraySize.push_back( sizeof( childArray_5 ) / sizeof( int ) );
/**
* Cell 6:
*
* 3 <--- Soma
* L 4
* L 6 * L 5
* L 2
* L 1
* L 0
*
*/
int childArray_6[ ] =
{
/* c0 */ -1,
/* c1 */ -1,
/* c2 */ -1,
/* c3 */ -1, 4,
/* c4 */ -1, 0, 1, 2, 5, 6,
/* c5 */ -1,
/* c6 */ -1,
};
childArray.push_back( childArray_6 );
childArraySize.push_back( sizeof( childArray_6 ) / sizeof( int ) );
/**
* Cell 7: Single compartment
*/
int childArray_7[ ] =
{
/* c0 */ -1,
};
childArray.push_back( childArray_7 );
childArraySize.push_back( sizeof( childArray_7 ) / sizeof( int ) );
/**
* Cell 8: 3 compartments; soma is in the middle.
*/
int childArray_8[ ] =
{
/* c0 */ -1,
/* c1 */ -1, 0, 2,
/* c2 */ -1,
};
childArray.push_back( childArray_8 );
childArraySize.push_back( sizeof( childArray_8 ) / sizeof( int ) );
/**
* Cell 9: 3 compartments; first compartment is soma.
*/
int childArray_9[ ] =
{
/* c0 */ -1, 1,
/* c1 */ -1, 2,
/* c2 */ -1,
};
childArray.push_back( childArray_9 );
childArraySize.push_back( sizeof( childArray_9 ) / sizeof( int ) );
////////////////////////////////////////////////////////////////////////////
// Run tests
////////////////////////////////////////////////////////////////////////////
HinesMatrix H;
vector< TreeNodeStruct > tree;
double dt = 1.0;
/*
* This is the full reference matrix which will be compared to its sparse * implementation.
*/
vector< vector< double > > matrix;
double epsilon = 1e-17;
unsigned int i;
unsigned int j;
unsigned int nCompt;
int* array;
unsigned int arraySize;
for ( unsigned int cell = 0; cell < childArray.size(); cell++ )
{
array = childArray[ cell ];
arraySize = childArraySize[ cell ];
nCompt = count( array, array + arraySize, -1 );
// Prepare cell
tree.clear();
tree.resize( nCompt );
for ( i = 0; i < nCompt; ++i )
{
tree[ i ].Ra = 15.0 + 3.0 * i;
tree[ i ].Rm = 45.0 + 15.0 * i;
tree[ i ].Cm = 500.0 + 200.0 * i * i;
}
int count = -1;
for ( unsigned int a = 0; a < arraySize; a++ )
if ( array[ a ] == -1 )
count++;
else
tree[ count ].children.push_back( array[ a ] );
// Prepare local matrix
makeFullMatrix( tree, dt, matrix );
// Prepare sparse matrix
H.setup( tree, dt );
// Compare matrices
for ( i = 0; i < nCompt; ++i )
for ( j = 0; j < nCompt; ++j )
{
ostringstream error;
error << "Testing Hines' Matrix: Cell# "
<< cell + 1 << ", entry (" << i << ", " << j << ")";
ASSERT(
fabs( matrix[ i ][ j ] - H.getA( i, j ) ) < epsilon,
error.str()
);
}
}
cout << "." << flush;
}
#endif // DO_UNIT_TESTS