diff --git a/CheckCXXCompiler.cmake b/CheckCXXCompiler.cmake
index 15a4d730b98679a10bf42e9c5f8df04d6de3fd58..a87acde8d8136fa8d7fc1e09369cce64b3f65cdb 100644
--- a/CheckCXXCompiler.cmake
+++ b/CheckCXXCompiler.cmake
@@ -37,7 +37,7 @@ if(COMPILER_SUPPORTS_CXX11)
elseif(COMPILER_SUPPORTS_CXX0X)
message(STATUS "Your compiler supports c++0x features. Enabling it")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x")
- add_definitions( -DENABLE_CXX11 )
+ add_definitions( -DENABLE_CPP11 )
add_definitions( -DBOOST_NO_CXX11_SCOPED_ENUMS -DBOOST_NO_SCOPED_ENUMS )
if(APPLE)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++" )
diff --git a/basecode/SparseMatrix.h b/basecode/SparseMatrix.h
index 059bb7939345ea6fe88fe0972162cbe6a6e1d7bb..ecbc71eafa78149c8e67e09de69f5a7ceba6a34e 100644
--- a/basecode/SparseMatrix.h
+++ b/basecode/SparseMatrix.h
@@ -582,7 +582,7 @@ public:
void tripletFill( const vector< unsigned int >& row,
const vector< unsigned int >& col,
- const vector< T >& z )
+ const vector< T >& z, bool retainSize = false )
{
unsigned int len = row.size();
if ( len > col.size() ) len = col.size();
@@ -591,16 +591,20 @@ public:
for ( unsigned int i = 0; i < len; ++i )
trip[i]= Triplet< T >(z[i], row[i], col[i] );
sort( trip.begin(), trip.end(), Triplet< T >::cmp );
- unsigned int nr = trip.back().b_ + 1;
- unsigned int nc = 0;
- for ( typename vector< Triplet< T > >::iterator i =
- trip.begin(); i != trip.end(); ++i )
- {
- if ( nc < i->c_ )
- nc = i->c_;
- }
- nc++;
- setSize( nr, nc );
+ unsigned int nr = nrows_;
+ unsigned int nc = ncolumns_;
+ if ( !retainSize ) {
+ nr = trip.back().b_ + 1;
+ nc = 0;
+ for ( typename vector< Triplet< T > >::iterator i =
+ trip.begin(); i != trip.end(); ++i )
+ {
+ if ( nc < i->c_ )
+ nc = i->c_;
+ }
+ nc++;
+ }
+ setSize( nr, nc );
vector< unsigned int > colIndex( nc );
vector< T > entry( nc );
diff --git a/ksolve/Gsolve.cpp b/ksolve/Gsolve.cpp
index a1c576ac18c65fe39c7570050a2f6456fa6f7123..736ca1f84a9499a35a6a43a1a36934aaf501ccc0 100644
--- a/ksolve/Gsolve.cpp
+++ b/ksolve/Gsolve.cpp
@@ -30,195 +30,199 @@
const unsigned int OFFNODE = ~0;
// static function
-SrcFinfo2< Id, vector< double > >* Gsolve::xComptOut() {
- static SrcFinfo2< Id, vector< double > > xComptOut( "xComptOut",
- "Sends 'n' of all molecules participating in cross-compartment "
- "reactions between any juxtaposed voxels between current compt "
- "and another compartment. This includes molecules local to this "
- "compartment, as well as proxy molecules belonging elsewhere. "
- "A(t+1) = (Alocal(t+1) + AremoteProxy(t+1)) - Alocal(t) "
- "A(t+1) = (Aremote(t+1) + Aproxy(t+1)) - Aproxy(t) "
- "Then we update A on the respective solvers with: "
- "Alocal(t+1) = Aproxy(t+1) = A(t+1) "
- "This is equivalent to sending dA over on each timestep. "
- );
- return &xComptOut;
+SrcFinfo2< Id, vector< double > >* Gsolve::xComptOut()
+{
+ static SrcFinfo2< Id, vector< double > > xComptOut( "xComptOut",
+ "Sends 'n' of all molecules participating in cross-compartment "
+ "reactions between any juxtaposed voxels between current compt "
+ "and another compartment. This includes molecules local to this "
+ "compartment, as well as proxy molecules belonging elsewhere. "
+ "A(t+1) = (Alocal(t+1) + AremoteProxy(t+1)) - Alocal(t) "
+ "A(t+1) = (Aremote(t+1) + Aproxy(t+1)) - Aproxy(t) "
+ "Then we update A on the respective solvers with: "
+ "Alocal(t+1) = Aproxy(t+1) = A(t+1) "
+ "This is equivalent to sending dA over on each timestep. "
+ );
+ return &xComptOut;
}
const Cinfo* Gsolve::initCinfo()
{
- ///////////////////////////////////////////////////////
- // Field definitions
- ///////////////////////////////////////////////////////
-
- static ValueFinfo< Gsolve, Id > stoich (
- "stoich",
- "Stoichiometry object for handling this reaction system.",
- &Gsolve::setStoich,
- &Gsolve::getStoich
- );
-
- static ValueFinfo< Gsolve, Id > compartment (
- "compartment",
- "Compartment that contains this reaction system.",
- &Gsolve::setCompartment,
- &Gsolve::getCompartment
- );
-
- static ReadOnlyValueFinfo< Gsolve, unsigned int > numLocalVoxels(
- "numLocalVoxels",
- "Number of voxels in the core reac-diff system, on the "
- "current solver. ",
- &Gsolve::getNumLocalVoxels
- );
- static LookupValueFinfo<
- Gsolve, unsigned int, vector< double > > nVec(
- "nVec",
- "vector of pool counts",
- &Gsolve::setNvec,
- &Gsolve::getNvec
- );
- static ValueFinfo< Gsolve, unsigned int > numAllVoxels(
- "numAllVoxels",
- "Number of voxels in the entire reac-diff system, "
- "including proxy voxels to represent abutting compartments.",
- &Gsolve::setNumAllVoxels,
- &Gsolve::getNumAllVoxels
- );
-
- static ValueFinfo< Gsolve, unsigned int > numPools(
- "numPools",
- "Number of molecular pools in the entire reac-diff system, "
- "including variable, function and buffered.",
- &Gsolve::setNumPools,
- &Gsolve::getNumPools
- );
-
- static ValueFinfo< Gsolve, bool > useRandInit(
- "useRandInit",
- "Flag: True when using probabilistic (random) rounding.\n "
- "Default: True.\n "
- "When initializing the mol# from floating-point Sinit values, "
- "we have two options. One is to look at each Sinit, and round "
- "to the nearest integer. The other is to look at each Sinit, "
- "and probabilistically round up or down depending on the "
- "value. For example, if we had a Sinit value of 1.49, "
- "this would always be rounded to 1.0 if the flag is false, "
- "and would be rounded to 1.0 and 2.0 in the ratio 51:49 if "
- "the flag is true. ",
- &Gsolve::setRandInit,
- &Gsolve::getRandInit
- );
-
- static ValueFinfo< Gsolve, bool > useClockedUpdate(
- "useClockedUpdate",
- "Flag: True to cause all reaction propensities to be updated "
- "on every clock tick.\n"
- "Default: False.\n"
- "This flag should be set when the reaction system "
- "includes a function with a dependency on time or on external "
- "events. It has a significant speed penalty so the flag "
- "should not be set unless there are such functions. " ,
- &Gsolve::setClockedUpdate,
- &Gsolve::getClockedUpdate
- );
- static ReadOnlyLookupValueFinfo<
- Gsolve, unsigned int, vector< unsigned int > > numFire(
- "numFire",
- "Vector of the number of times each reaction has fired."
- "Indexed by the voxel number."
- "Zeroed out at reinit.",
- &Gsolve::getNumFire
- );
-
- ///////////////////////////////////////////////////////
- // DestFinfo definitions
- ///////////////////////////////////////////////////////
-
- static DestFinfo process( "process",
- "Handles process call",
- new ProcOpFunc< Gsolve >( &Gsolve::process ) );
- static DestFinfo reinit( "reinit",
- "Handles reinit call",
- new ProcOpFunc< Gsolve >( &Gsolve::reinit ) );
-
- static DestFinfo voxelVol( "voxelVol",
- "Handles updates to all voxels. Comes from parent "
- "ChemCompt object.",
- new OpFunc1< Gsolve, vector< double > >(
- &Gsolve::updateVoxelVol )
- );
-
- static DestFinfo initProc( "initProc",
- "Handles initProc call from Clock",
- new ProcOpFunc< Gsolve >( &Gsolve::initProc ) );
- static DestFinfo initReinit( "initReinit",
- "Handles initReinit call from Clock",
- new ProcOpFunc< Gsolve >( &Gsolve::initReinit ) );
-
- static DestFinfo xComptIn( "xComptIn",
- "Handles arriving pool 'n' values used in cross-compartment "
- "reactions.",
- new EpFunc2< Gsolve, Id, vector< double > >( &Gsolve::xComptIn )
- );
- ///////////////////////////////////////////////////////
- // Shared definitions
- ///////////////////////////////////////////////////////
- static Finfo* procShared[] = {
- &process, &reinit
- };
- static SharedFinfo proc( "proc",
- "Shared message for process and reinit",
- procShared, sizeof( procShared ) / sizeof( const Finfo* )
- );
-
- static Finfo* initShared[] = {
- &initProc, &initReinit
- };
- static SharedFinfo init( "init",
- "Shared message for initProc and initReinit. This is used"
- " when the system has cross-compartment reactions. ",
- initShared, sizeof( initShared ) / sizeof( const Finfo* )
- );
-
- static Finfo* xComptShared[] = {
- xComptOut(), &xComptIn
- };
- static SharedFinfo xCompt( "xCompt",
- "Shared message for pool exchange for cross-compartment "
- "reactions. Exchanges latest values of all pools that "
- "participate in such reactions.",
- xComptShared, sizeof( xComptShared ) / sizeof( const Finfo* )
- );
- ///////////////////////////////////////////////////////
-
- static Finfo* gsolveFinfos[] =
- {
- &stoich, // Value
- &numLocalVoxels, // ReadOnlyValue
- &nVec, // LookupValue
- &numAllVoxels, // ReadOnlyValue
- &numPools, // Value
- &voxelVol, // DestFinfo
- &proc, // SharedFinfo
- &init, // SharedFinfo
- &xCompt, // SharedFinfo
- // Here we put new fields that were not there in the Ksolve.
- &useRandInit, // Value
- &useClockedUpdate, // Value
- &numFire, // ReadOnlyLookupValue
- };
-
- static Dinfo< Gsolve > dinfo;
- static Cinfo gsolveCinfo(
- "Gsolve",
- Neutral::initCinfo(),
- gsolveFinfos,
- sizeof(gsolveFinfos)/sizeof(Finfo *),
- &dinfo
- );
-
- return &gsolveCinfo;
+ ///////////////////////////////////////////////////////
+ // Field definitions
+ ///////////////////////////////////////////////////////
+
+ static ValueFinfo< Gsolve, Id > stoich (
+ "stoich",
+ "Stoichiometry object for handling this reaction system.",
+ &Gsolve::setStoich,
+ &Gsolve::getStoich
+ );
+
+ static ValueFinfo< Gsolve, Id > compartment (
+ "compartment",
+ "Compartment that contains this reaction system.",
+ &Gsolve::setCompartment,
+ &Gsolve::getCompartment
+ );
+
+ static ReadOnlyValueFinfo< Gsolve, unsigned int > numLocalVoxels(
+ "numLocalVoxels",
+ "Number of voxels in the core reac-diff system, on the "
+ "current solver. ",
+ &Gsolve::getNumLocalVoxels
+ );
+ static LookupValueFinfo<
+ Gsolve, unsigned int, vector< double > > nVec(
+ "nVec",
+ "vector of pool counts",
+ &Gsolve::setNvec,
+ &Gsolve::getNvec
+ );
+ static ValueFinfo< Gsolve, unsigned int > numAllVoxels(
+ "numAllVoxels",
+ "Number of voxels in the entire reac-diff system, "
+ "including proxy voxels to represent abutting compartments.",
+ &Gsolve::setNumAllVoxels,
+ &Gsolve::getNumAllVoxels
+ );
+
+ static ValueFinfo< Gsolve, unsigned int > numPools(
+ "numPools",
+ "Number of molecular pools in the entire reac-diff system, "
+ "including variable, function and buffered.",
+ &Gsolve::setNumPools,
+ &Gsolve::getNumPools
+ );
+
+ static ValueFinfo< Gsolve, bool > useRandInit(
+ "useRandInit",
+ "Flag: True when using probabilistic (random) rounding.\n "
+ "Default: True.\n "
+ "When initializing the mol# from floating-point Sinit values, "
+ "we have two options. One is to look at each Sinit, and round "
+ "to the nearest integer. The other is to look at each Sinit, "
+ "and probabilistically round up or down depending on the "
+ "value. For example, if we had a Sinit value of 1.49, "
+ "this would always be rounded to 1.0 if the flag is false, "
+ "and would be rounded to 1.0 and 2.0 in the ratio 51:49 if "
+ "the flag is true. ",
+ &Gsolve::setRandInit,
+ &Gsolve::getRandInit
+ );
+
+ static ValueFinfo< Gsolve, bool > useClockedUpdate(
+ "useClockedUpdate",
+ "Flag: True to cause all reaction propensities to be updated "
+ "on every clock tick.\n"
+ "Default: False.\n"
+ "This flag should be set when the reaction system "
+ "includes a function with a dependency on time or on external "
+ "events. It has a significant speed penalty so the flag "
+ "should not be set unless there are such functions. " ,
+ &Gsolve::setClockedUpdate,
+ &Gsolve::getClockedUpdate
+ );
+ static ReadOnlyLookupValueFinfo<
+ Gsolve, unsigned int, vector< unsigned int > > numFire(
+ "numFire",
+ "Vector of the number of times each reaction has fired."
+ "Indexed by the voxel number."
+ "Zeroed out at reinit.",
+ &Gsolve::getNumFire
+ );
+
+ ///////////////////////////////////////////////////////
+ // DestFinfo definitions
+ ///////////////////////////////////////////////////////
+
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< Gsolve >( &Gsolve::process ) );
+ static DestFinfo reinit( "reinit",
+ "Handles reinit call",
+ new ProcOpFunc< Gsolve >( &Gsolve::reinit ) );
+
+ static DestFinfo voxelVol( "voxelVol",
+ "Handles updates to all voxels. Comes from parent "
+ "ChemCompt object.",
+ new OpFunc1< Gsolve, vector< double > >(
+ &Gsolve::updateVoxelVol )
+ );
+
+ static DestFinfo initProc( "initProc",
+ "Handles initProc call from Clock",
+ new ProcOpFunc< Gsolve >( &Gsolve::initProc ) );
+ static DestFinfo initReinit( "initReinit",
+ "Handles initReinit call from Clock",
+ new ProcOpFunc< Gsolve >( &Gsolve::initReinit ) );
+
+ static DestFinfo xComptIn( "xComptIn",
+ "Handles arriving pool 'n' values used in cross-compartment "
+ "reactions.",
+ new EpFunc2< Gsolve, Id, vector< double > >( &Gsolve::xComptIn )
+ );
+ ///////////////////////////////////////////////////////
+ // Shared definitions
+ ///////////////////////////////////////////////////////
+ static Finfo* procShared[] =
+ {
+ &process, &reinit
+ };
+ static SharedFinfo proc( "proc",
+ "Shared message for process and reinit",
+ procShared, sizeof( procShared ) / sizeof( const Finfo* )
+ );
+
+ static Finfo* initShared[] =
+ {
+ &initProc, &initReinit
+ };
+ static SharedFinfo init( "init",
+ "Shared message for initProc and initReinit. This is used"
+ " when the system has cross-compartment reactions. ",
+ initShared, sizeof( initShared ) / sizeof( const Finfo* )
+ );
+
+ static Finfo* xComptShared[] =
+ {
+ xComptOut(), &xComptIn
+ };
+ static SharedFinfo xCompt( "xCompt",
+ "Shared message for pool exchange for cross-compartment "
+ "reactions. Exchanges latest values of all pools that "
+ "participate in such reactions.",
+ xComptShared, sizeof( xComptShared ) / sizeof( const Finfo* )
+ );
+ ///////////////////////////////////////////////////////
+
+ static Finfo* gsolveFinfos[] =
+ {
+ &stoich, // Value
+ &numLocalVoxels, // ReadOnlyValue
+ &nVec, // LookupValue
+ &numAllVoxels, // ReadOnlyValue
+ &numPools, // Value
+ &voxelVol, // DestFinfo
+ &proc, // SharedFinfo
+ &init, // SharedFinfo
+ &xCompt, // SharedFinfo
+ // Here we put new fields that were not there in the Ksolve.
+ &useRandInit, // Value
+ &useClockedUpdate, // Value
+ &numFire, // ReadOnlyLookupValue
+ };
+
+ static Dinfo< Gsolve > dinfo;
+ static Cinfo gsolveCinfo(
+ "Gsolve",
+ Neutral::initCinfo(),
+ gsolveFinfos,
+ sizeof(gsolveFinfos)/sizeof(Finfo *),
+ &dinfo
+ );
+
+ return &gsolveCinfo;
}
static const Cinfo* gsolveCinfo = Gsolve::initCinfo();
@@ -228,16 +232,20 @@ static const Cinfo* gsolveCinfo = Gsolve::initCinfo();
//////////////////////////////////////////////////////////////
Gsolve::Gsolve()
- :
- pools_( 1 ),
- startVoxel_( 0 ),
- dsolve_(),
- dsolvePtr_( 0 ),
- useClockedUpdate_( false )
-{;}
+ :
+ pools_( 1 ),
+ startVoxel_( 0 ),
+ dsolve_(),
+ dsolvePtr_( 0 ),
+ useClockedUpdate_( false )
+{
+ ;
+}
Gsolve::~Gsolve()
-{;}
+{
+ ;
+}
//////////////////////////////////////////////////////////////
// Field Access functions
@@ -245,121 +253,131 @@ Gsolve::~Gsolve()
Id Gsolve::getStoich() const
{
- return stoich_;
+ return stoich_;
}
void Gsolve::setCompartment( Id compt )
{
- if ( ( compt.element()->cinfo()->isA( "ChemCompt" ) ) ) {
- compartment_ = compt;
- vector< double > vols =
- Field< vector< double > >::get( compt, "voxelVolume" );
- if ( vols.size() > 0 ) {
- pools_.resize( vols.size() );
- for ( unsigned int i = 0; i < vols.size(); ++i ) {
- pools_[i].setVolume( vols[i] );
- }
- }
- }
+ if ( ( compt.element()->cinfo()->isA( "ChemCompt" ) ) )
+ {
+ compartment_ = compt;
+ vector< double > vols =
+ Field< vector< double > >::get( compt, "voxelVolume" );
+ if ( vols.size() > 0 )
+ {
+ pools_.resize( vols.size() );
+ for ( unsigned int i = 0; i < vols.size(); ++i )
+ {
+ pools_[i].setVolume( vols[i] );
+ }
+ }
+ }
}
Id Gsolve::getCompartment() const
{
- return compartment_;
+ return compartment_;
}
void Gsolve::setStoich( Id stoich )
{
- // This call is done _before_ setting the path on stoich
- assert( stoich.element()->cinfo()->isA( "Stoich" ) );
- stoich_ = stoich;
- stoichPtr_ = reinterpret_cast< Stoich* >( stoich.eref().data() );
- if ( stoichPtr_->getNumAllPools() == 0 ) {
- stoichPtr_ = 0;
- return;
- }
- sys_.stoich = stoichPtr_;
- sys_.isReady = false;
- for ( unsigned int i = 0; i < pools_.size(); ++i )
- pools_[i].setStoich( stoichPtr_ );
+ // This call is done _before_ setting the path on stoich
+ assert( stoich.element()->cinfo()->isA( "Stoich" ) );
+ stoich_ = stoich;
+ stoichPtr_ = reinterpret_cast< Stoich* >( stoich.eref().data() );
+ if ( stoichPtr_->getNumAllPools() == 0 )
+ {
+ stoichPtr_ = 0;
+ return;
+ }
+ sys_.stoich = stoichPtr_;
+ sys_.isReady = false;
+ for ( unsigned int i = 0; i < pools_.size(); ++i )
+ pools_[i].setStoich( stoichPtr_ );
}
unsigned int Gsolve::getNumLocalVoxels() const
{
- return pools_.size();
+ return pools_.size();
}
unsigned int Gsolve::getNumAllVoxels() const
{
- return pools_.size(); // Need to redo.
+ return pools_.size(); // Need to redo.
}
// If we're going to do this, should be done before the zombification.
void Gsolve::setNumAllVoxels( unsigned int numVoxels )
{
- if ( numVoxels == 0 ) {
- return;
- }
- pools_.resize( numVoxels );
- sys_.isReady = false;
+ if ( numVoxels == 0 )
+ {
+ return;
+ }
+ pools_.resize( numVoxels );
+ sys_.isReady = false;
}
vector< double > Gsolve::getNvec( unsigned int voxel) const
{
- static vector< double > dummy;
- if ( voxel < pools_.size() ) {
- return const_cast< GssaVoxelPools* >( &( pools_[ voxel ]) )->Svec();
- }
- return dummy;
+ static vector< double > dummy;
+ if ( voxel < pools_.size() )
+ {
+ return const_cast< GssaVoxelPools* >( &( pools_[ voxel ]) )->Svec();
+ }
+ return dummy;
}
void Gsolve::setNvec( unsigned int voxel, vector< double > nVec )
{
- if ( voxel < pools_.size() ) {
- if ( nVec.size() != pools_[voxel].size() ) {
- cout << "Warning: Gsolve::setNvec: size mismatch ( " <<
- nVec.size() << ", " << pools_[voxel].size() << ")\n";
- return;
- }
- double* s = pools_[voxel].varS();
- for ( unsigned int i = 0; i < nVec.size(); ++i ) {
- s[i] = round( nVec[i] );
- if ( s[i] < 0.0 )
- s[i] = 0.0;
- }
- if ( sys_.isReady )
- pools_[voxel].refreshAtot( &sys_ );
- }
+ if ( voxel < pools_.size() )
+ {
+ if ( nVec.size() != pools_[voxel].size() )
+ {
+ cout << "Warning: Gsolve::setNvec: size mismatch ( " <<
+ nVec.size() << ", " << pools_[voxel].size() << ")\n";
+ return;
+ }
+ double* s = pools_[voxel].varS();
+ for ( unsigned int i = 0; i < nVec.size(); ++i )
+ {
+ s[i] = round( nVec[i] );
+ if ( s[i] < 0.0 )
+ s[i] = 0.0;
+ }
+ if ( sys_.isReady )
+ pools_[voxel].refreshAtot( &sys_ );
+ }
}
vector< unsigned int > Gsolve::getNumFire( unsigned int voxel) const
{
- static vector< unsigned int > dummy;
- if ( voxel < pools_.size() ) {
- return const_cast< GssaVoxelPools* >( &( pools_[ voxel ]) )->numFire();
- }
- return dummy;
+ static vector< unsigned int > dummy;
+ if ( voxel < pools_.size() )
+ {
+ return const_cast< GssaVoxelPools* >( &( pools_[ voxel ]) )->numFire();
+ }
+ return dummy;
}
bool Gsolve::getRandInit() const
{
- return sys_.useRandInit;
+ return sys_.useRandInit;
}
void Gsolve::setRandInit( bool val )
{
- sys_.useRandInit = val;
+ sys_.useRandInit = val;
}
bool Gsolve::getClockedUpdate() const
{
- return useClockedUpdate_;
+ return useClockedUpdate_;
}
void Gsolve::setClockedUpdate( bool val )
{
- useClockedUpdate_ = val;
+ useClockedUpdate_ = val;
}
//////////////////////////////////////////////////////////////
@@ -367,122 +385,143 @@ void Gsolve::setClockedUpdate( bool val )
//////////////////////////////////////////////////////////////
void Gsolve::process( const Eref& e, ProcPtr p )
{
- // cout << stoichPtr_ << " dsolve = " << dsolvePtr_ << endl;
- if ( !stoichPtr_ )
- return;
- // First, handle incoming diffusion values. Note potential for
- // issues with roundoff if diffusion is not integral.
- if ( dsolvePtr_ ) {
- vector< double > dvalues( 4 );
- dvalues[0] = 0;
- dvalues[1] = getNumLocalVoxels();
- dvalues[2] = 0;
- dvalues[3] = stoichPtr_->getNumVarPools();
- dsolvePtr_->getBlock( dvalues );
- // Here we need to convert to integers, just in case. Normally
- // one would use a stochastic (integral) diffusion method with
- // the GSSA, but in mixed models it may be more complicated.
- vector< double >::iterator i = dvalues.begin() + 4;
- for ( ; i != dvalues.end(); ++i ) {
- // cout << *i << " " << round( *i ) << " ";
+ // cout << stoichPtr_ << " dsolve = " << dsolvePtr_ << endl;
+ if ( !stoichPtr_ )
+ return;
+
+ // First, handle incoming diffusion values. Note potential for
+ // issues with roundoff if diffusion is not integral.
+ if ( dsolvePtr_ )
+ {
+ vector< double > dvalues( 4 );
+ dvalues[0] = 0;
+ dvalues[1] = getNumLocalVoxels();
+ dvalues[2] = 0;
+ dvalues[3] = stoichPtr_->getNumVarPools();
+ dsolvePtr_->getBlock( dvalues );
+
+ // Here we need to convert to integers, just in case. Normally
+ // one would use a stochastic (integral) diffusion method with
+ // the GSSA, but in mixed models it may be more complicated.
+ vector< double >::iterator i = dvalues.begin() + 4;
+
+ for ( ; i != dvalues.end(); ++i )
+ {
+ // cout << *i << " " << round( *i ) << " ";
#if SIMPLE_ROUNDING
- *i = round( *i );
+ *i = round( *i );
#else
- double base = floor( *i );
- if ( mtrand() > *i - base )
- *i = base;
- else
- *i = base + 1.0;
+ double base = floor( *i );
+ if ( mtrand() >= (*i - base) )
+ *i = base;
+ else
+ *i = base + 1.0;
#endif
- }
- setBlock( dvalues );
- }
- // Second, take the arrived xCompt reac values and update S with them.
- // Here the roundoff issues are handled by the GssaVoxelPools functions
- for ( unsigned int i = 0; i < xfer_.size(); ++i ) {
- XferInfo& xf = xfer_[i];
- // cout << xfer_.size() << " " << xf.xferVoxel.size() << endl;
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j ) {
- pools_[xf.xferVoxel[j]].xferIn( xf, j, &sys_ );
- }
- }
- // Third, record the current value of pools as the reference for the
- // next cycle.
- for ( unsigned int i = 0; i < xfer_.size(); ++i ) {
- XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j ) {
- pools_[xf.xferVoxel[j]].xferOut( j, xf.lastValues, xf.xferPoolIdx );
- }
- }
-
- // Fourth: Fix the rates if we have had any diffusion or xreacs
- // happening. This is very inefficient at this point, need to fix.
- if ( dsolvePtr_ || xfer_.size() > 0 ) {
- for ( vector< GssaVoxelPools >::iterator
- i = pools_.begin(); i != pools_.end(); ++i ) {
- i->refreshAtot( &sys_ );
- }
- }
- // Fifth, update the mol #s.
- // First we advance the simulation.
- for ( vector< GssaVoxelPools >::iterator
- i = pools_.begin(); i != pools_.end(); ++i ) {
- i->advance( p, &sys_ );
- }
- if ( useClockedUpdate_ ) { // Check if a clocked stim is to be updated
- for ( vector< GssaVoxelPools >::iterator
- i = pools_.begin(); i != pools_.end(); ++i ) {
- i->recalcTime( &sys_, p->currTime );
- }
- }
-
- // Finally, assemble and send the integrated values off for the Dsolve.
- if ( dsolvePtr_ ) {
- vector< double > kvalues( 4 );
- kvalues[0] = 0;
- kvalues[1] = getNumLocalVoxels();
- kvalues[2] = 0;
- kvalues[3] = stoichPtr_->getNumVarPools();
- getBlock( kvalues );
- dsolvePtr_->setBlock( kvalues );
- }
+ }
+ setBlock( dvalues );
+ }
+ // Second, take the arrived xCompt reac values and update S with them.
+ // Here the roundoff issues are handled by the GssaVoxelPools functions
+ for ( unsigned int i = 0; i < xfer_.size(); ++i )
+ {
+ XferInfo& xf = xfer_[i];
+ // cout << xfer_.size() << " " << xf.xferVoxel.size() << endl;
+ for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
+ {
+ pools_[xf.xferVoxel[j]].xferIn( xf, j, &sys_ );
+ }
+ }
+ // Third, record the current value of pools as the reference for the
+ // next cycle.
+ for ( unsigned int i = 0; i < xfer_.size(); ++i )
+ {
+ XferInfo& xf = xfer_[i];
+ for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
+ {
+ pools_[xf.xferVoxel[j]].xferOut( j, xf.lastValues, xf.xferPoolIdx );
+ }
+ }
+
+ // Fourth: Fix the rates if we have had any diffusion or xreacs
+ // happening. This is very inefficient at this point, need to fix.
+ if ( dsolvePtr_ || xfer_.size() > 0 )
+ {
+ for ( vector< GssaVoxelPools >::iterator
+ i = pools_.begin(); i != pools_.end(); ++i )
+ {
+ i->refreshAtot( &sys_ );
+ }
+ }
+ // Fifth, update the mol #s.
+ // First we advance the simulation.
+ for ( vector< GssaVoxelPools >::iterator
+ i = pools_.begin(); i != pools_.end(); ++i )
+ {
+ i->advance( p, &sys_ );
+ }
+ if ( useClockedUpdate_ ) // Check if a clocked stim is to be updated
+ {
+ for ( vector< GssaVoxelPools >::iterator
+ i = pools_.begin(); i != pools_.end(); ++i )
+ {
+ i->recalcTime( &sys_, p->currTime );
+ }
+ }
+
+ // Finally, assemble and send the integrated values off for the Dsolve.
+ if ( dsolvePtr_ )
+ {
+ vector< double > kvalues( 4 );
+ kvalues[0] = 0;
+ kvalues[1] = getNumLocalVoxels();
+ kvalues[2] = 0;
+ kvalues[3] = stoichPtr_->getNumVarPools();
+ getBlock( kvalues );
+ dsolvePtr_->setBlock( kvalues );
+ }
}
void Gsolve::reinit( const Eref& e, ProcPtr p )
{
- if ( !stoichPtr_ )
- return;
- if ( !sys_.isReady )
- rebuildGssaSystem();
- // First reinit concs.
- for ( vector< GssaVoxelPools >::iterator
- i = pools_.begin(); i != pools_.end(); ++i ) {
- i->reinit( &sys_ );
- }
-
- // Second, take the arrived xCompt reac values and update S with them.
- // Here the roundoff issues are handled by the GssaVoxelPools functions
- for ( unsigned int i = 0; i < xfer_.size(); ++i ) {
- const XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j ) {
- pools_[xf.xferVoxel[j]].xferInOnlyProxies(
- xf.xferPoolIdx, xf.values,
- stoichPtr_->getNumProxyPools(), j );
- }
- }
- // Third, record the current value of pools as the reference for the
- // next cycle.
- for ( unsigned int i = 0; i < xfer_.size(); ++i ) {
- XferInfo& xf = xfer_[i];
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j ) {
- pools_[xf.xferVoxel[j]].xferOut( j, xf.lastValues, xf.xferPoolIdx );
- }
- }
- // Fourth, update the atots.
- for ( vector< GssaVoxelPools >::iterator
- i = pools_.begin(); i != pools_.end(); ++i ) {
- i->refreshAtot( &sys_ );
- }
+ if ( !stoichPtr_ )
+ return;
+ if ( !sys_.isReady )
+ rebuildGssaSystem();
+ // First reinit concs.
+ for ( vector< GssaVoxelPools >::iterator
+ i = pools_.begin(); i != pools_.end(); ++i )
+ {
+ i->reinit( &sys_ );
+ }
+
+ // Second, take the arrived xCompt reac values and update S with them.
+ // Here the roundoff issues are handled by the GssaVoxelPools functions
+ for ( unsigned int i = 0; i < xfer_.size(); ++i )
+ {
+ const XferInfo& xf = xfer_[i];
+ for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
+ {
+ pools_[xf.xferVoxel[j]].xferInOnlyProxies(
+ xf.xferPoolIdx, xf.values,
+ stoichPtr_->getNumProxyPools(), j );
+ }
+ }
+ // Third, record the current value of pools as the reference for the
+ // next cycle.
+ for ( unsigned int i = 0; i < xfer_.size(); ++i )
+ {
+ XferInfo& xf = xfer_[i];
+ for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
+ {
+ pools_[xf.xferVoxel[j]].xferOut( j, xf.lastValues, xf.xferPoolIdx );
+ }
+ }
+ // Fourth, update the atots.
+ for ( vector< GssaVoxelPools >::iterator
+ i = pools_.begin(); i != pools_.end(); ++i )
+ {
+ i->refreshAtot( &sys_ );
+ }
}
//////////////////////////////////////////////////////////////
@@ -490,41 +529,46 @@ void Gsolve::reinit( const Eref& e, ProcPtr p )
//////////////////////////////////////////////////////////////
void Gsolve::initProc( const Eref& e, ProcPtr p )
{
- if ( !stoichPtr_ )
- return;
- // vector< vector< double > > values( xfer_.size() );
- for ( unsigned int i = 0; i < xfer_.size(); ++i ) {
- XferInfo& xf = xfer_[i];
- unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
- // values[i].resize( size, 0.0 );
- vector< double > values( size, 0.0 );
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j ) {
- unsigned int vox = xf.xferVoxel[j];
- pools_[vox].xferOut( j, values, xf.xferPoolIdx );
- }
- xComptOut()->sendTo( e, xf.ksolve, e.id(), values );
- }
+ if ( !stoichPtr_ )
+ return;
+ // vector< vector< double > > values( xfer_.size() );
+ for ( unsigned int i = 0; i < xfer_.size(); ++i )
+ {
+ XferInfo& xf = xfer_[i];
+ unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
+ // values[i].resize( size, 0.0 );
+ vector< double > values( size, 0.0 );
+ for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
+ {
+ unsigned int vox = xf.xferVoxel[j];
+ pools_[vox].xferOut( j, values, xf.xferPoolIdx );
+ }
+ xComptOut()->sendTo( e, xf.ksolve, e.id(), values );
+ }
}
void Gsolve::initReinit( const Eref& e, ProcPtr p )
{
- if ( !stoichPtr_ )
- return;
- for ( unsigned int i = 0 ; i < pools_.size(); ++i ) {
- pools_[i].reinit( &sys_ );
- }
- // vector< vector< double > > values( xfer_.size() );
- for ( unsigned int i = 0; i < xfer_.size(); ++i ) {
- XferInfo& xf = xfer_[i];
- unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
- xf.lastValues.assign( size, 0.0 );
- for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j ) {
- unsigned int vox = xf.xferVoxel[j];
- pools_[ vox ].xferOut( j, xf.lastValues, xf.xferPoolIdx );
- // values[i] = xf.lastValues;
- }
- xComptOut()->sendTo( e, xf.ksolve, e.id(), xf.lastValues );
- }
+ if ( !stoichPtr_ )
+ return;
+ for ( unsigned int i = 0 ; i < pools_.size(); ++i )
+ {
+ pools_[i].reinit( &sys_ );
+ }
+ // vector< vector< double > > values( xfer_.size() );
+ for ( unsigned int i = 0; i < xfer_.size(); ++i )
+ {
+ XferInfo& xf = xfer_[i];
+ unsigned int size = xf.xferPoolIdx.size() * xf.xferVoxel.size();
+ xf.lastValues.assign( size, 0.0 );
+ for ( unsigned int j = 0; j < xf.xferVoxel.size(); ++j )
+ {
+ unsigned int vox = xf.xferVoxel[j];
+ pools_[ vox ].xferOut( j, xf.lastValues, xf.xferPoolIdx );
+ // values[i] = xf.lastValues;
+ }
+ xComptOut()->sendTo( e, xf.ksolve, e.id(), xf.lastValues );
+ }
}
//////////////////////////////////////////////////////////////
// Solver setup
@@ -532,26 +576,28 @@ void Gsolve::initReinit( const Eref& e, ProcPtr p )
void Gsolve::rebuildGssaSystem()
{
- stoichPtr_->convertRatesToStochasticForm();
- sys_.transposeN = stoichPtr_->getStoichiometryMatrix();
- sys_.transposeN.transpose();
- sys_.transposeN.truncateRow( stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools() );
- vector< vector< unsigned int > > & dep = sys_.dependency;
- dep.resize( stoichPtr_->getNumRates() );
- for ( unsigned int i = 0; i < stoichPtr_->getNumRates(); ++i ) {
- sys_.transposeN.getGillespieDependence( i, dep[i] );
- }
- fillMmEnzDep();
- fillPoolFuncDep();
- fillIncrementFuncDep();
- makeReacDepsUnique();
- for ( vector< GssaVoxelPools >::iterator
- i = pools_.begin(); i != pools_.end(); ++i ) {
- i->setNumReac( stoichPtr_->getNumRates() );
- i->updateAllRateTerms( stoichPtr_->getRateTerms(),
- stoichPtr_->getNumCoreRates() );
- }
- sys_.isReady = true;
+ stoichPtr_->convertRatesToStochasticForm();
+ sys_.transposeN = stoichPtr_->getStoichiometryMatrix();
+ sys_.transposeN.transpose();
+ sys_.transposeN.truncateRow( stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools() );
+ vector< vector< unsigned int > > & dep = sys_.dependency;
+ dep.resize( stoichPtr_->getNumRates() );
+ for ( unsigned int i = 0; i < stoichPtr_->getNumRates(); ++i )
+ {
+ sys_.transposeN.getGillespieDependence( i, dep[i] );
+ }
+ fillMmEnzDep();
+ fillPoolFuncDep();
+ fillIncrementFuncDep();
+ makeReacDepsUnique();
+ for ( vector< GssaVoxelPools >::iterator
+ i = pools_.begin(); i != pools_.end(); ++i )
+ {
+ i->setNumReac( stoichPtr_->getNumRates() );
+ i->updateAllRateTerms( stoichPtr_->getRateTerms(),
+ stoichPtr_->getNumCoreRates() );
+ }
+ sys_.isReady = true;
}
/**
@@ -562,44 +608,48 @@ void Gsolve::rebuildGssaSystem()
*/
void Gsolve::fillMmEnzDep()
{
- unsigned int numRates = stoichPtr_->getNumRates();
- vector< unsigned int > indices;
-
- // Make a map to look up enzyme RateTerm using
- // the key of the enzyme molecule.
- map< unsigned int, unsigned int > enzMolMap;
- for ( unsigned int i = 0; i < numRates; ++i ) {
- const MMEnzymeBase* mme = dynamic_cast< const MMEnzymeBase* >(
- stoichPtr_->rates( i ) );
- if ( mme ) {
- vector< unsigned int > reactants;
- mme->getReactants( reactants );
- if ( reactants.size() > 1 )
- enzMolMap[ reactants.front() ] = i; // front is enzyme.
- }
- }
-
- // Use the map to fill in deps.
- for ( unsigned int i = 0; i < numRates; ++i ) {
- // Extract the row of all molecules that depend on the reac.
- const int* entry;
- const unsigned int* colIndex;
-
- unsigned int numInRow =
- sys_.transposeN.getRow( i, &entry, &colIndex );
- for( unsigned int j = 0; j < numInRow; ++j ) {
- map< unsigned int, unsigned int >::iterator pos =
- enzMolMap.find( colIndex[j] );
- if ( pos != enzMolMap.end() )
- sys_.dependency[i].push_back( pos->second );
- }
- }
+ unsigned int numRates = stoichPtr_->getNumRates();
+ vector< unsigned int > indices;
+
+ // Make a map to look up enzyme RateTerm using
+ // the key of the enzyme molecule.
+ map< unsigned int, unsigned int > enzMolMap;
+ for ( unsigned int i = 0; i < numRates; ++i )
+ {
+ const MMEnzymeBase* mme = dynamic_cast< const MMEnzymeBase* >(
+ stoichPtr_->rates( i ) );
+ if ( mme )
+ {
+ vector< unsigned int > reactants;
+ mme->getReactants( reactants );
+ if ( reactants.size() > 1 )
+ enzMolMap[ reactants.front() ] = i; // front is enzyme.
+ }
+ }
+
+ // Use the map to fill in deps.
+ for ( unsigned int i = 0; i < numRates; ++i )
+ {
+ // Extract the row of all molecules that depend on the reac.
+ const int* entry;
+ const unsigned int* colIndex;
+
+ unsigned int numInRow =
+ sys_.transposeN.getRow( i, &entry, &colIndex );
+ for( unsigned int j = 0; j < numInRow; ++j )
+ {
+ map< unsigned int, unsigned int >::iterator pos =
+ enzMolMap.find( colIndex[j] );
+ if ( pos != enzMolMap.end() )
+ sys_.dependency[i].push_back( pos->second );
+ }
+ }
}
/**
* Here we fill in the dependencies involving poolFuncs. These are
* the functions that evaluate an expression and assign directly to the
- * # of a target molecule.
+ * # of a target molecule.
* There are two dependencies:
* 1. When a reaction fires, all the Functions that depend on the reactants
* must update their target molecule. This is in sys_.dependentMathExpn[].
@@ -609,55 +659,59 @@ void Gsolve::fillMmEnzDep()
*/
void Gsolve::fillPoolFuncDep()
{
- // create map of funcs that depend on specified molecule.
- vector< vector< unsigned int > > funcMap(
- stoichPtr_->getNumAllPools() );
- unsigned int numFuncs = stoichPtr_->getNumFuncs();
- for ( unsigned int i = 0; i < numFuncs; ++i ) {
- const FuncTerm *f = stoichPtr_->funcs( i );
- vector< unsigned int > molIndex = f->getReactantIndex();
- for ( unsigned int j = 0; j < molIndex.size(); ++j )
- funcMap[ molIndex[j] ].push_back( i );
- }
- // The output of each func is a mol indexed as
- // numVarMols + numBufMols + i
- unsigned int numRates = stoichPtr_->getNumRates();
- sys_.dependentMathExpn.resize( numRates );
- vector< unsigned int > indices;
- for ( unsigned int i = 0; i < numRates; ++i ) {
- vector< unsigned int >& dep = sys_.dependentMathExpn[ i ];
- dep.resize( 0 );
- // Extract the row of all molecules that depend on the reac.
- const int* entry;
- const unsigned int* colIndex;
- unsigned int numInRow =
- sys_.transposeN.getRow( i, &entry, &colIndex );
- for ( unsigned int j = 0; j < numInRow; ++j ) {
- unsigned int molIndex = colIndex[j];
- vector< unsigned int >& funcs = funcMap[ molIndex ];
- dep.insert( dep.end(), funcs.begin(), funcs.end() );
- for ( unsigned int k = 0; k < funcs.size(); ++k ) {
- // unsigned int outputMol = funcs[k] + funcOffset;
- unsigned int outputMol = stoichPtr_->funcs( funcs[k] )->getTarget();
- // Insert reac deps here. Columns are reactions.
- vector< int > e; // Entries: we don't need.
- vector< unsigned int > c; // Column index: the reactions.
- stoichPtr_->getStoichiometryMatrix().
- getRow( outputMol, e, c );
- // Each of the reacs (col entries) depend on this func.
- vector< unsigned int > rdep = sys_.dependency[i];
- rdep.insert( rdep.end(), c.begin(), c.end() );
- }
- }
- }
+ // create map of funcs that depend on specified molecule.
+ vector< vector< unsigned int > > funcMap(
+ stoichPtr_->getNumAllPools() );
+ unsigned int numFuncs = stoichPtr_->getNumFuncs();
+ for ( unsigned int i = 0; i < numFuncs; ++i )
+ {
+ const FuncTerm *f = stoichPtr_->funcs( i );
+ vector< unsigned int > molIndex = f->getReactantIndex();
+ for ( unsigned int j = 0; j < molIndex.size(); ++j )
+ funcMap[ molIndex[j] ].push_back( i );
+ }
+ // The output of each func is a mol indexed as
+ // numVarMols + numBufMols + i
+ unsigned int numRates = stoichPtr_->getNumRates();
+ sys_.dependentMathExpn.resize( numRates );
+ vector< unsigned int > indices;
+ for ( unsigned int i = 0; i < numRates; ++i )
+ {
+ vector< unsigned int >& dep = sys_.dependentMathExpn[ i ];
+ dep.resize( 0 );
+ // Extract the row of all molecules that depend on the reac.
+ const int* entry;
+ const unsigned int* colIndex;
+ unsigned int numInRow =
+ sys_.transposeN.getRow( i, &entry, &colIndex );
+ for ( unsigned int j = 0; j < numInRow; ++j )
+ {
+ unsigned int molIndex = colIndex[j];
+ vector< unsigned int >& funcs = funcMap[ molIndex ];
+ dep.insert( dep.end(), funcs.begin(), funcs.end() );
+ for ( unsigned int k = 0; k < funcs.size(); ++k )
+ {
+ // unsigned int outputMol = funcs[k] + funcOffset;
+ unsigned int outputMol = stoichPtr_->funcs( funcs[k] )->getTarget();
+ // Insert reac deps here. Columns are reactions.
+ vector< int > e; // Entries: we don't need.
+ vector< unsigned int > c; // Column index: the reactions.
+ stoichPtr_->getStoichiometryMatrix().
+ getRow( outputMol, e, c );
+ // Each of the reacs (col entries) depend on this func.
+ vector< unsigned int > rdep = sys_.dependency[i];
+ rdep.insert( rdep.end(), c.begin(), c.end() );
+ }
+ }
+ }
}
/**
* Here we fill in the dependencies involving incrementFuncs. These are
* the functions that evaluate an expression that specifies rate of change
- * of # of a target molecule.
+ * of # of a target molecule.
* There are two dependencies:
- * 1. When a reaction fires, all the incrementFuncs that depend on the
+ * 1. When a reaction fires, all the incrementFuncs that depend on the
* reactants must update their rates. This is added to sys_.dependency[]
* which is the usual handler for reac dependencies. Note that the inputs
* to the incrementFuncs are NOT present in the stoichiometry matrix, so
@@ -668,62 +722,67 @@ void Gsolve::fillPoolFuncDep()
*/
void Gsolve::fillIncrementFuncDep()
{
- // create map of funcs that depend on specified molecule.
- vector< vector< unsigned int > > funcMap(
- stoichPtr_->getNumAllPools() );
- const vector< RateTerm* >& rates = stoichPtr_->getRateTerms();
- vector< FuncRate* > incrementRates;
- vector< unsigned int > incrementRateIndex;
- const vector< RateTerm* >::const_iterator q;
- for ( unsigned int i = 0; i < rates.size(); ++i ) {
- FuncRate *term =
- dynamic_cast< FuncRate* >( rates[i] );
- if (term) {
- incrementRates.push_back( term );
- incrementRateIndex.push_back( i );
- }
- }
-
- for ( unsigned int k = 0; k < incrementRates.size(); ++k ) {
- const vector< unsigned int >& molIndex =
- incrementRates[k]->getFuncArgIndex();
- for ( unsigned int j = 0; j < molIndex.size(); ++j )
- funcMap[ molIndex[j] ].push_back( incrementRateIndex[k] );
- }
-
- unsigned int numRates = stoichPtr_->getNumRates();
- sys_.dependentMathExpn.resize( numRates );
- vector< unsigned int > indices;
- for ( unsigned int i = 0; i < numRates; ++i ) {
- // Algorithm:
- // 1.Go through stoich matrix finding all the poolIndices affected
- // by each Rate Term.
- // 2.Use funcMap to look up FuncRateTerms affected by these indices
- // 3. Add the rateTerm->FuncRateTerm mapping to the dependencies.
-
- const int* entry;
- const unsigned int* colIndex;
- unsigned int numInRow =
- sys_.transposeN.getRow( i, &entry, &colIndex );
- // 1.Go through stoich matrix finding all the poolIndices affected
- // by each Rate Term.
- for ( unsigned int j = 0; j < numInRow; ++j ) {
- unsigned int molIndex = colIndex[j]; // Affected poolIndex
-
- // 2.Use funcMap to look up FuncRateTerms affected by these indices
- vector< unsigned int >& funcs = funcMap[ molIndex ];
- // 3. Add the rateTerm->FuncRateTerm mapping to the dependencies.
- vector< unsigned int >& rdep = sys_.dependency[i];
- rdep.insert( rdep.end(), funcs.begin(), funcs.end() );
- }
- }
+ // create map of funcs that depend on specified molecule.
+ vector< vector< unsigned int > > funcMap(
+ stoichPtr_->getNumAllPools() );
+ const vector< RateTerm* >& rates = stoichPtr_->getRateTerms();
+ vector< FuncRate* > incrementRates;
+ vector< unsigned int > incrementRateIndex;
+ const vector< RateTerm* >::const_iterator q;
+ for ( unsigned int i = 0; i < rates.size(); ++i )
+ {
+ FuncRate *term =
+ dynamic_cast< FuncRate* >( rates[i] );
+ if (term)
+ {
+ incrementRates.push_back( term );
+ incrementRateIndex.push_back( i );
+ }
+ }
+
+ for ( unsigned int k = 0; k < incrementRates.size(); ++k )
+ {
+ const vector< unsigned int >& molIndex =
+ incrementRates[k]->getFuncArgIndex();
+ for ( unsigned int j = 0; j < molIndex.size(); ++j )
+ funcMap[ molIndex[j] ].push_back( incrementRateIndex[k] );
+ }
+
+ unsigned int numRates = stoichPtr_->getNumRates();
+ sys_.dependentMathExpn.resize( numRates );
+ vector< unsigned int > indices;
+ for ( unsigned int i = 0; i < numRates; ++i )
+ {
+ // Algorithm:
+ // 1.Go through stoich matrix finding all the poolIndices affected
+ // by each Rate Term.
+ // 2.Use funcMap to look up FuncRateTerms affected by these indices
+ // 3. Add the rateTerm->FuncRateTerm mapping to the dependencies.
+
+ const int* entry;
+ const unsigned int* colIndex;
+ unsigned int numInRow =
+ sys_.transposeN.getRow( i, &entry, &colIndex );
+ // 1.Go through stoich matrix finding all the poolIndices affected
+ // by each Rate Term.
+ for ( unsigned int j = 0; j < numInRow; ++j )
+ {
+ unsigned int molIndex = colIndex[j]; // Affected poolIndex
+
+ // 2.Use funcMap to look up FuncRateTerms affected by these indices
+ vector< unsigned int >& funcs = funcMap[ molIndex ];
+ // 3. Add the rateTerm->FuncRateTerm mapping to the dependencies.
+ vector< unsigned int >& rdep = sys_.dependency[i];
+ rdep.insert( rdep.end(), funcs.begin(), funcs.end() );
+ }
+ }
}
/*
void Gsolve::fillMathDep()
{
// create map of funcs that depend on specified molecule.
- vector< vector< unsigned int > > funcMap(
+ vector< vector< unsigned int > > funcMap(
stoichPtr_->getNumAllPools() );
unsigned int numFuncs = stoichPtr_->getNumFuncs();
for ( unsigned int i = 0; i < numFuncs; ++i ) {
@@ -732,9 +791,9 @@ void Gsolve::fillMathDep()
for ( unsigned int j = 0; j < molIndex.size(); ++j )
funcMap[ molIndex[j] ].push_back( i );
}
- // The output of each func is a mol indexed as
+ // The output of each func is a mol indexed as
// numVarMols + numBufMols + i
- unsigned int funcOffset =
+ unsigned int funcOffset =
stoichPtr_->getNumVarPools() + stoichPtr_->getNumProxyPools() + stoichPtr_->getNumBufPools();
unsigned int numRates = stoichPtr_->getNumRates();
sys_.dependentMathExpn.resize( numRates );
@@ -745,7 +804,7 @@ void Gsolve::fillMathDep()
// Extract the row of all molecules that depend on the reac.
const int* entry;
const unsigned int* colIndex;
- unsigned int numInRow =
+ unsigned int numInRow =
sys_.transposeN.getRow( i, &entry, &colIndex );
for ( unsigned int j = 0; j < numInRow; ++j ) {
unsigned int molIndex = colIndex[j];
@@ -773,38 +832,39 @@ void Gsolve::fillMathDep()
* Later.
*/
void Gsolve::insertMathDepReacs( unsigned int mathDepIndex,
- unsigned int firedReac )
+ unsigned int firedReac )
{
- /*
- unsigned int molIndex = sumTotals_[ mathDepIndex ].target( S_ );
- vector< unsigned int > reacIndices;
-
- // Extract the row of all reacs that depend on the target molecule
- if ( N_.getRowIndices( molIndex, reacIndices ) > 0 ) {
- vector< unsigned int >& dep = dependency_[ firedReac ];
- dep.insert( dep.end(), reacIndices.begin(), reacIndices.end() );
- }
- */
+ /*
+ unsigned int molIndex = sumTotals_[ mathDepIndex ].target( S_ );
+ vector< unsigned int > reacIndices;
+
+ // Extract the row of all reacs that depend on the target molecule
+ if ( N_.getRowIndices( molIndex, reacIndices ) > 0 ) {
+ vector< unsigned int >& dep = dependency_[ firedReac ];
+ dep.insert( dep.end(), reacIndices.begin(), reacIndices.end() );
+ }
+ */
}
// Clean up dependency lists: Ensure only unique entries.
// Also a reac cannot depend on itself.
void Gsolve::makeReacDepsUnique()
{
- unsigned int numRates = stoichPtr_->getNumRates();
- for ( unsigned int i = 0; i < numRates; ++i ) {
- vector< unsigned int >& dep = sys_.dependency[ i ];
- // Here we want to remove self-entries as well as duplicates.
- sort( dep.begin(), dep.end() );
- vector< unsigned int >::iterator k = dep.begin();
-
- /// STL stuff follows, with the usual weirdness.
- vector< unsigned int >::iterator pos =
- unique( dep.begin(), dep.end() );
- dep.resize( pos - dep.begin() );
- /*
- */
- }
+ unsigned int numRates = stoichPtr_->getNumRates();
+ for ( unsigned int i = 0; i < numRates; ++i )
+ {
+ vector< unsigned int >& dep = sys_.dependency[ i ];
+ // Here we want to remove self-entries as well as duplicates.
+ sort( dep.begin(), dep.end() );
+ vector< unsigned int >::iterator k = dep.begin();
+
+ /// STL stuff follows, with the usual weirdness.
+ vector< unsigned int >::iterator pos =
+ unique( dep.begin(), dep.end() );
+ dep.resize( pos - dep.begin() );
+ /*
+ */
+ }
}
//////////////////////////////////////////////////////////////
@@ -812,31 +872,36 @@ void Gsolve::makeReacDepsUnique()
//////////////////////////////////////////////////////////////
unsigned int Gsolve::getPoolIndex( const Eref& e ) const
{
- return stoichPtr_->convertIdToPoolIndex( e.id() );
+ return stoichPtr_->convertIdToPoolIndex( e.id() );
}
unsigned int Gsolve::getVoxelIndex( const Eref& e ) const
{
- unsigned int ret = e.dataIndex();
- if ( ret < startVoxel_ || ret >= startVoxel_ + pools_.size() )
- return OFFNODE;
- return ret - startVoxel_;
+ unsigned int ret = e.dataIndex();
+ if ( ret < startVoxel_ || ret >= startVoxel_ + pools_.size() )
+ return OFFNODE;
+ return ret - startVoxel_;
}
void Gsolve::setDsolve( Id dsolve )
{
- if ( dsolve == Id () ) {
- dsolvePtr_ = 0;
- dsolve_ = Id();
- } else if ( dsolve.element()->cinfo()->isA( "Dsolve" ) ) {
- dsolve_ = dsolve;
- dsolvePtr_ = reinterpret_cast< ZombiePoolInterface* >(
- dsolve.eref().data() );
- } else {
- cout << "Warning: Gsolve::setDsolve: Object '" << dsolve.path() <<
- "' should be class Dsolve, is: " <<
- dsolve.element()->cinfo()->name() << endl;
- }
+ if ( dsolve == Id () )
+ {
+ dsolvePtr_ = 0;
+ dsolve_ = Id();
+ }
+ else if ( dsolve.element()->cinfo()->isA( "Dsolve" ) )
+ {
+ dsolve_ = dsolve;
+ dsolvePtr_ = reinterpret_cast< ZombiePoolInterface* >(
+ dsolve.eref().data() );
+ }
+ else
+ {
+ cout << "Warning: Gsolve::setDsolve: Object '" << dsolve.path() <<
+ "' should be class Dsolve, is: " <<
+ dsolve.element()->cinfo()->name() << endl;
+ }
}
@@ -846,166 +911,185 @@ void Gsolve::setDsolve( Id dsolve )
void Gsolve::setN( const Eref& e, double v )
{
- unsigned int vox = getVoxelIndex( e );
- if ( vox != OFFNODE ) {
- if ( e.element()->cinfo()->isA( "ZombieBufPool" ) ) {
- // Do NOT round it here, it is folded into rate term.
- pools_[vox].setN( getPoolIndex( e ), v );
- // refresh rates because nInit controls ongoing value of n.
- if ( sys_.isReady )
- pools_[vox].refreshAtot( &sys_ );
- } else {
- pools_[vox].setN( getPoolIndex( e ), round( v ) );
- }
- }
+ unsigned int vox = getVoxelIndex( e );
+ if ( vox != OFFNODE )
+ {
+ if ( e.element()->cinfo()->isA( "ZombieBufPool" ) )
+ {
+ // Do NOT round it here, it is folded into rate term.
+ pools_[vox].setN( getPoolIndex( e ), v );
+ // refresh rates because nInit controls ongoing value of n.
+ if ( sys_.isReady )
+ pools_[vox].refreshAtot( &sys_ );
+ }
+ else
+ {
+ pools_[vox].setN( getPoolIndex( e ), round( v ) );
+ }
+ }
}
double Gsolve::getN( const Eref& e ) const
{
- unsigned int vox = getVoxelIndex( e );
- if ( vox != OFFNODE )
- return pools_[vox].getN( getPoolIndex( e ) );
- return 0.0;
+ unsigned int vox = getVoxelIndex( e );
+ if ( vox != OFFNODE )
+ return pools_[vox].getN( getPoolIndex( e ) );
+ return 0.0;
}
void Gsolve::setNinit( const Eref& e, double v )
{
- unsigned int vox = getVoxelIndex( e );
- if ( vox != OFFNODE ) {
- if ( e.element()->cinfo()->isA( "ZombieBufPool" ) ) {
- // Do NOT round it here, it is folded into rate term.
- pools_[vox].setNinit( getPoolIndex( e ), v );
- // refresh rates because nInit controls ongoing value of n.
- if ( sys_.isReady )
- pools_[vox].refreshAtot( &sys_ );
- } else {
- // I now do the rounding at reinit time. It is better there as
- // it can give a distinct value each cycle. It is also better
- // to keep the full resolution of Ninit for volume scaling.
- // pools_[vox].setNinit( getPoolIndex( e ), round( v ) );
- pools_[vox].setNinit( getPoolIndex( e ), v );
- }
- }
+ unsigned int vox = getVoxelIndex( e );
+ if ( vox != OFFNODE )
+ {
+ if ( e.element()->cinfo()->isA( "ZombieBufPool" ) )
+ {
+ // Do NOT round it here, it is folded into rate term.
+ pools_[vox].setNinit( getPoolIndex( e ), v );
+ // refresh rates because nInit controls ongoing value of n.
+ if ( sys_.isReady )
+ pools_[vox].refreshAtot( &sys_ );
+ }
+ else
+ {
+ // I now do the rounding at reinit time. It is better there as
+ // it can give a distinct value each cycle. It is also better
+ // to keep the full resolution of Ninit for volume scaling.
+ // pools_[vox].setNinit( getPoolIndex( e ), round( v ) );
+ pools_[vox].setNinit( getPoolIndex( e ), v );
+ }
+ }
}
double Gsolve::getNinit( const Eref& e ) const
{
- unsigned int vox = getVoxelIndex( e );
- if ( vox != OFFNODE )
- return pools_[vox].getNinit( getPoolIndex( e ) );
- return 0.0;
+ unsigned int vox = getVoxelIndex( e );
+ if ( vox != OFFNODE )
+ return pools_[vox].getNinit( getPoolIndex( e ) );
+ return 0.0;
}
void Gsolve::setDiffConst( const Eref& e, double v )
{
- ; // Do nothing.
+ ; // Do nothing.
}
double Gsolve::getDiffConst( const Eref& e ) const
{
- return 0;
+ return 0;
}
void Gsolve::setNumPools( unsigned int numPoolSpecies )
{
- sys_.isReady = false;
- unsigned int numVoxels = pools_.size();
- for ( unsigned int i = 0 ; i < numVoxels; ++i ) {
- pools_[i].resizeArrays( numPoolSpecies );
- }
+ sys_.isReady = false;
+ unsigned int numVoxels = pools_.size();
+ for ( unsigned int i = 0 ; i < numVoxels; ++i )
+ {
+ pools_[i].resizeArrays( numPoolSpecies );
+ }
}
unsigned int Gsolve::getNumPools() const
{
- if ( pools_.size() > 0 )
- return pools_[0].size();
- return 0;
+ if ( pools_.size() > 0 )
+ return pools_[0].size();
+ return 0;
}
void Gsolve::getBlock( vector< double >& values ) const
{
- unsigned int startVoxel = values[0];
- unsigned int numVoxels = values[1];
- unsigned int startPool = values[2];
- unsigned int numPools = values[3];
-
- assert( startVoxel >= startVoxel_ );
- assert( numVoxels <= pools_.size() );
- assert( pools_.size() > 0 );
- assert( numPools + startPool <= pools_[0].size() );
- values.resize( 4 + numVoxels * numPools );
-
- for ( unsigned int i = 0; i < numVoxels; ++i ) {
- const double* v = pools_[ startVoxel + i ].S();
- for ( unsigned int j = 0; j < numPools; ++j ) {
- values[ 4 + j * numVoxels + i] = v[ j + startPool ];
- }
- }
+ unsigned int startVoxel = values[0];
+ unsigned int numVoxels = values[1];
+ unsigned int startPool = values[2];
+ unsigned int numPools = values[3];
+
+ assert( startVoxel >= startVoxel_ );
+ assert( numVoxels <= pools_.size() );
+ assert( pools_.size() > 0 );
+ assert( numPools + startPool <= pools_[0].size() );
+ values.resize( 4 + numVoxels * numPools );
+
+ for ( unsigned int i = 0; i < numVoxels; ++i )
+ {
+ const double* v = pools_[ startVoxel + i ].S();
+ for ( unsigned int j = 0; j < numPools; ++j )
+ {
+ values[ 4 + j * numVoxels + i] = v[ j + startPool ];
+ }
+ }
}
void Gsolve::setBlock( const vector< double >& values )
{
- unsigned int startVoxel = values[0];
- unsigned int numVoxels = values[1];
- unsigned int startPool = values[2];
- unsigned int numPools = values[3];
-
- assert( startVoxel >= startVoxel_ );
- assert( numVoxels <= pools_.size() );
- assert( pools_.size() > 0 );
- assert( numPools + startPool <= pools_[0].size() );
-
- for ( unsigned int i = 0; i < numVoxels; ++i ) {
- double* v = pools_[ startVoxel + i ].varS();
- for ( unsigned int j = 0; j < numPools; ++j ) {
- v[ j + startPool ] = values[ 4 + j * numVoxels + i ];
- }
- }
+ unsigned int startVoxel = values[0];
+ unsigned int numVoxels = values[1];
+ unsigned int startPool = values[2];
+ unsigned int numPools = values[3];
+
+ assert( startVoxel >= startVoxel_ );
+ assert( numVoxels <= pools_.size() );
+ assert( pools_.size() > 0 );
+ assert( numPools + startPool <= pools_[0].size() );
+
+ for ( unsigned int i = 0; i < numVoxels; ++i )
+ {
+ double* v = pools_[ startVoxel + i ].varS();
+ for ( unsigned int j = 0; j < numPools; ++j )
+ {
+ v[ j + startPool ] = values[ 4 + j * numVoxels + i ];
+ }
+ }
}
//////////////////////////////////////////////////////////////////////////
void Gsolve::updateVoxelVol( vector< double > vols )
{
- // For now we assume identical numbers of voxels. Also assume
- // identical voxel junctions. But it should not be too hard to
- // update those too.
- if ( vols.size() == pools_.size() ) {
- for ( unsigned int i = 0; i < vols.size(); ++i ) {
- pools_[i].setVolumeAndDependencies( vols[i] );
- }
- stoichPtr_->setupCrossSolverReacVols();
- updateRateTerms( ~0U );
- }
+ // For now we assume identical numbers of voxels. Also assume
+ // identical voxel junctions. But it should not be too hard to
+ // update those too.
+ if ( vols.size() == pools_.size() )
+ {
+ for ( unsigned int i = 0; i < vols.size(); ++i )
+ {
+ pools_[i].setVolumeAndDependencies( vols[i] );
+ }
+ stoichPtr_->setupCrossSolverReacVols();
+ updateRateTerms( ~0U );
+ }
}
void Gsolve::updateRateTerms( unsigned int index )
{
- if ( index == ~0U ) {
- // unsigned int numCrossRates = stoichPtr_->getNumRates() - stoichPtr_->getNumCoreRates();
- for ( unsigned int i = 0 ; i < pools_.size(); ++i ) {
- // pools_[i].resetXreacScale( numCrossRates );
- pools_[i].updateAllRateTerms( stoichPtr_->getRateTerms(),
- stoichPtr_->getNumCoreRates() );
- }
- } else if ( index < stoichPtr_->getNumRates() ) {
- for ( unsigned int i = 0 ; i < pools_.size(); ++i )
- pools_[i].updateRateTerms( stoichPtr_->getRateTerms(),
- stoichPtr_->getNumCoreRates(), index );
- }
+ if ( index == ~0U )
+ {
+ // unsigned int numCrossRates = stoichPtr_->getNumRates() - stoichPtr_->getNumCoreRates();
+ for ( unsigned int i = 0 ; i < pools_.size(); ++i )
+ {
+ // pools_[i].resetXreacScale( numCrossRates );
+ pools_[i].updateAllRateTerms( stoichPtr_->getRateTerms(),
+ stoichPtr_->getNumCoreRates() );
+ }
+ }
+ else if ( index < stoichPtr_->getNumRates() )
+ {
+ for ( unsigned int i = 0 ; i < pools_.size(); ++i )
+ pools_[i].updateRateTerms( stoichPtr_->getRateTerms(),
+ stoichPtr_->getNumCoreRates(), index );
+ }
}
//////////////////////////////////////////////////////////////////////////
VoxelPoolsBase* Gsolve::pools( unsigned int i )
{
- if ( pools_.size() > i )
- return &pools_[i];
- return 0;
+ if ( pools_.size() > i )
+ return &pools_[i];
+ return 0;
}
double Gsolve::volume( unsigned int i ) const
{
- if ( pools_.size() > i )
- return pools_[i].getVolume();
- return 0.0;
+ if ( pools_.size() > i )
+ return pools_[i].getVolume();
+ return 0.0;
}
diff --git a/ksolve/GssaVoxelPools.cpp b/ksolve/GssaVoxelPools.cpp
index b74900aa7bf9183474b51e4144c234bb4eaf183a..810c9481ee376fd7f87732bad02a3ee88c2b99f3 100644
--- a/ksolve/GssaVoxelPools.cpp
+++ b/ksolve/GssaVoxelPools.cpp
@@ -205,10 +205,10 @@ void GssaVoxelPools::advance( const ProcInfo* p, const GssaSystem* g )
numFire_[rindex]++;
double r = rng_.uniform();
- while ( r <= 0.0 )
- {
+
+ while ( r == 0.0 )
r = rng_.uniform();
- }
+
t_ -= ( 1.0 / atot_ ) * log( r );
// g->stoich->updateFuncs( varS(), t_ ); // Handled next line.
updateDependentMathExpn( g, rindex, t_ );
@@ -234,7 +234,7 @@ void GssaVoxelPools::reinit( const GssaSystem* g )
double base = floor( n[i] );
assert( base >= 0.0 );
double frac = n[i] - base;
- if ( rng_.uniform() > frac )
+ if ( rng_.uniform() >= frac )
n[i] = base;
else
n[i] = base + 1.0;
@@ -369,9 +369,9 @@ void GssaVoxelPools::xferIn( XferInfo& xf,
{
double& x = s[*k];
// cout << x << " i = " << *i << *j << " m = " << *m << endl;
- double dx = *i++ - *j++;
+ double dx = *(i++) - *(j++);
double base = floor( dx );
- if ( rng_.uniform() > dx - base )
+ if ( rng_.uniform() >= (dx - base) )
x += base;
else
x += base + 1.0;
@@ -425,7 +425,7 @@ void GssaVoxelPools::xferInOnlyProxies(
if ( *k >= stoichPtr_->getNumVarPools() && *k < proxyEndIndex )
{
double base = floor( *i );
- if ( rng_.uniform() > *i - base )
+ if ( rng_.uniform() >= (*i - base) )
varSinit()[*k] = (varS()[*k] += base );
else
varSinit()[*k] = (varS()[*k] += base + 1.0 );
diff --git a/msg/SparseMsg.cpp b/msg/SparseMsg.cpp
index 12d9da4d8f841fc8a7998e56f02cd912f8b3d50b..46d3fb5e1fb4b0a3da80bb34fc7ecaf0b54c1596 100644
--- a/msg/SparseMsg.cpp
+++ b/msg/SparseMsg.cpp
@@ -41,6 +41,17 @@ const Cinfo* SparseMsg::initCinfo()
"Number of Entries in matrix.",
&SparseMsg::getNumEntries
);
+
+ static ValueFinfo< SparseMsg, vector< unsigned int > > connectionList(
+ "connectionList",
+ "Pairwise specification of connection matrix where each x,y value "
+ "represents a connection from src[x] to dest[y]. "
+ "The (x,y) entries are ordered in a single vector as \n"
+ "(x0, x1,... x_n-1, y0, y1,... y_n-1)\n",
+ &SparseMsg::setEntryPairs,
+ &SparseMsg::getEntryPairs
+ );
+
/// Connection matrix entries to manipulate in Python.
static ReadOnlyValueFinfo< SparseMsg, vector< unsigned int > >
matrixEntry(
@@ -124,6 +135,13 @@ const Cinfo* SparseMsg::initCinfo()
vector< unsigned int > >(
&SparseMsg::tripletFill ) );
+ static DestFinfo tripletFill1( "tripletFill1",
+ "Single contiguous array to fill entire connection matrix using "
+ "triplets of (x,y, fieldindex) ordered as \n"
+ "(x0, x1,... xn-1, y0, y1,... yn-1, fi0, fi1,... fi_n-1)\n",
+ new OpFunc1< SparseMsg, vector< unsigned int > >(
+ &SparseMsg::tripletFill1 ) );
+
////////////////////////////////////////////////////////////////////////
// Assemble it all.
////////////////////////////////////////////////////////////////////////
@@ -132,6 +150,7 @@ const Cinfo* SparseMsg::initCinfo()
&numRows, // readonly value
&numColumns, // readonly value
&numEntries, // readonly value
+ &connectionList, // value
&matrixEntry, // ReadOnlyValue
&columnIndex, // ReadOnlyValue
&rowStart, // ReadOnlyValue
@@ -144,6 +163,7 @@ const Cinfo* SparseMsg::initCinfo()
&transpose, //dest
&pairFill, //dest
&tripletFill, //dest
+ &tripletFill1, //dest
};
static Dinfo< short > dinfo;
@@ -217,6 +237,28 @@ vector< unsigned int > SparseMsg::getRowStart() const
return matrix_.rowStart();
}
+void SparseMsg::setEntryPairs( vector< unsigned int > v )
+{
+ vector< unsigned int > src( v.begin(), v.begin() + v.size()/2 );
+ vector< unsigned int > dest( v.begin() + v.size()/2, v.end() );
+ pairFill( src, dest );
+}
+
+vector< unsigned int > SparseMsg::getEntryPairs() const
+{
+ vector< unsigned int > cols = matrix_.colIndex();
+ vector< unsigned int > y;
+ for ( unsigned int row = 0; row < matrix_.nRows(); ++row ) {
+ unsigned int begin = matrix_.rowStart()[row];
+ unsigned int end = matrix_.rowStart()[row+1];
+ for ( unsigned int j = begin; j < end; ++j )
+ y.push_back( row );
+ }
+ assert( cols.size() == y.size() );
+ y.insert( y.end(), cols.begin(), cols.end() );
+ return y;
+}
+
//////////////////////////////////////////////////////////////////
// DestFields
//////////////////////////////////////////////////////////////////
@@ -256,21 +298,57 @@ void SparseMsg::updateAfterFill()
{
unsigned int startData = e2_->localDataStart();
unsigned int endData = startData + e2_->numLocalData();
- for ( unsigned int i = 0; i < matrix_.nRows(); ++ i ) {
+ SparseMatrix< unsigned int > temp( matrix_);
+ temp.transpose();
+ for ( unsigned int i = 0; i < temp.nRows(); ++ i ) {
const unsigned int* colIndex;
const unsigned int* entry;
- unsigned int num = matrix_.getRow( i, &entry, &colIndex );
+ unsigned int num = temp.getRow( i, &entry, &colIndex );
if ( i >= startData && i < endData ) {
- e2_->resizeField( i - startData, num );
+ // Inefficient. Better to do it in one pass after getting
+ // the max num
+ e2_->resizeField( i - startData, num + 1 );
}
}
e1()->markRewired();
e2()->markRewired();
}
+
void SparseMsg::pairFill( vector< unsigned int > src,
vector< unsigned int> dest )
{
- matrix_.pairFill( src, dest, 0 );
+ // Sanity check
+ vector< unsigned int >::const_iterator i;
+ for ( i = src.begin(); i != src.end(); ++i ) {
+ if (*i >= e1()->numData() ) {
+ cout << "Warning: SparseMsg::pairFill: Src index " << *i <<
+ " exceeds Src array size " << e1()->numData() <<
+ ". Aborting\n";
+ return;
+ }
+ }
+ for ( i = dest.begin(); i != dest.end(); ++i ) {
+ if (*i >= e2()->numData() ) {
+ cout << "Warning: SparseMsg::pairFill: Dest index " << *i <<
+ " exceeds Dest array size " << e2()->numData() <<
+ ". Aborting\n";
+ return;
+ }
+ }
+
+ vector< unsigned int > numAtDest( dest.size(), 0 );
+ vector< unsigned int > fieldIndex( dest.size(), 0 );
+ for ( unsigned int i = 0; i < dest.size(); ++i ) {
+ fieldIndex[i] = numAtDest[ dest[i] ];
+ // Could do on previous line, but clarity
+ ++numAtDest[ dest[i] ];
+ }
+
+ /**
+ * We set retainSize flag to true since the # of src/dest objects
+ * doesn't change. We can explicitly assign it elsewhere if needed.
+ */
+ matrix_.tripletFill( src, dest, fieldIndex, true );
updateAfterFill();
}
@@ -278,10 +356,20 @@ void SparseMsg::tripletFill( vector< unsigned int > src,
vector< unsigned int> destDataIndex,
vector< unsigned int> destFieldIndex )
{
- matrix_.tripletFill( src, destDataIndex, destFieldIndex );
+ // We set retainSize flag to true
+ matrix_.tripletFill( src, destDataIndex, destFieldIndex, true );
updateAfterFill();
}
+void SparseMsg::tripletFill1( vector< unsigned int > v )
+{
+ unsigned int s3 = v.size() / 3;
+ vector< unsigned int > src( v.begin(), v.begin() + s3 );
+ vector< unsigned int > dest( v.begin() + s3, v.begin() + 2 * s3 );
+ vector< unsigned int > fieldIndex( v.begin() + 2 * s3, v.end() );
+ tripletFill( src, dest, fieldIndex );
+}
+
//////////////////////////////////////////////////////////////////
// Here are the actual class functions
//////////////////////////////////////////////////////////////////
@@ -289,7 +377,10 @@ void SparseMsg::tripletFill( vector< unsigned int > src,
SparseMsg::SparseMsg( Element* e1, Element* e2, unsigned int msgIndex )
: Msg( ObjId( managerId_, (msgIndex != 0) ? msgIndex: msg_.size() ),
- e1, e2 )
+ e1, e2 ),
+ numThreads_( 1 ),
+ nrows_( 0 ),
+ p_( 0.0 ), seed_( 0 )
{
unsigned int nrows = 0;
unsigned int ncolumns = 0;
diff --git a/msg/SparseMsg.h b/msg/SparseMsg.h
index ac300cfca851d7e408f4c84e08af860168406490..d6fbf80b04306c1c333a94ff96b1f3649a1fa787 100644
--- a/msg/SparseMsg.h
+++ b/msg/SparseMsg.h
@@ -82,6 +82,9 @@ class SparseMsg: public Msg
long getSeed() const;
void setSeed( long value );
+ vector< unsigned int > getEntryPairs() const;
+ void setEntryPairs( vector< unsigned int > entries );
+
void setEntry( unsigned int row, unsigned int column,
unsigned int value );
@@ -113,6 +116,14 @@ class SparseMsg: public Msg
void tripletFill( vector< unsigned int > src,
vector< unsigned int> dest,
vector< unsigned int > field );
+
+ /**
+ * Fills up the entire message based on triplets of
+ * src,destDataIndex,destFieldIndex, but catenates them all into
+ * a single long vector since PyMoose doesn't know how to handle
+ * multiple vectors.
+ */
+ void tripletFill1( vector< unsigned int > entries );
/**
* Utility function to update all sorts of values after we've
diff --git a/python/moose/__init__.py b/python/moose/__init__.py
index 02f1c5e6349151de0b6a3c96ebf63d730817faef..b5447c6cea7f6509c37bcd5a2f80ed0ba1b65ebf 100644
--- a/python/moose/__init__.py
+++ b/python/moose/__init__.py
@@ -13,6 +13,7 @@ import pydoc
import os
from io import StringIO
from collections import defaultdict
+import textwrap
import moose
from moose._moose import *
@@ -308,6 +309,45 @@ def showmsg(el):
msg.e2.path,
msg.destFieldsOnE2)
+def getfielddoc_(tokens, infoType=''):
+ indent=''
+ classname = tokens[0]
+ # print( classname )
+ while True:
+ try:
+ classelement = _moose.element('/classes/' + classname)
+ # print( classelement )
+ x = ''
+ for finfo in classelement.children:
+ # print (str(finfo)[33:-1])
+ # trial = 'path=/classes[0]/' + classname + '[0]/' + infoType
+ # print (trial)
+ # print (str(finfo)[33:-1] == trial)
+ for fieldelement in finfo:
+ fieldname = fieldelement.fieldName
+ if (str(fieldname).startswith('get')) or (str(fieldname).startswith('set')):
+ continue
+ baseinfo = ''
+ finfotype = fieldelement.name
+ if(infoType == 'destFinfo'):
+ say = 'Destination Message Field'
+ elif(infoType == 'valueFinfo'):
+ say = 'Value Field'
+ elif(infoType == 'srcFinfo'):
+ say = 'Source Message Field'
+ elif(infoType == 'lookupFinfo'):
+ say = 'Lookup Field'
+ elif(infoType == 'sharedFinfo'):
+ say = 'Shared Message Field'
+
+ # finfotype={finfotype}{baseinfo}
+ if (str(finfotype) == infoType):
+ x = x + '{indent}{classname}.{fieldname}: type= {type} finfotype= {say}{baseinfo}\n\t{docs}\n\n'.format(indent=indent, classname=tokens[0], fieldname = fieldname, type=fieldelement.type, say=say, baseinfo=baseinfo, docs=fieldelement.docs)
+
+ return x
+ # classname = classelement.baseClass
+ except ValueError:
+ raise NameError('`%s` has no field called `%s`'% (tokens[0], tokens[1]))
def getfielddoc(tokens, indent=''):
"""Return the documentation for field specified by `tokens`.
@@ -340,6 +380,7 @@ def getfielddoc(tokens, indent=''):
try:
classelement = _moose.element('/classes/' + classname)
for finfo in classelement.children:
+ x = ''
for fieldelement in finfo:
baseinfo = ''
if classname != tokens[0]:
@@ -348,18 +389,11 @@ def getfielddoc(tokens, indent=''):
# The field elements are
# /classes/{ParentClass}[0]/{fieldElementType}[N].
finfotype = fieldelement.name
- return '{indent}{classname}.{fieldname}: type={type}, finfotype={finfotype}{baseinfo}\n\t{docs}\n'.format(
- indent=indent, classname=tokens[0],
- fieldname=fieldname,
- type=fieldelement.type,
- finfotype=finfotype,
- baseinfo=baseinfo,
- docs=fieldelement.docs)
+ x = x + '{indent}{classname}.{fieldname}: type={type}, finfotype={finfotype}{baseinfo}\n\t{docs}\n'.format(indent=indent, classname=tokens[0], fieldname=fieldname, type=fieldelement.type, finfotype=finfotype, baseinfo=baseinfo, docs=fieldelement.docs)
+ return x
classname = classelement.baseClass
except ValueError:
- raise NameError('`%s` has no field called `%s`'
- % (tokens[0], tokens[1]))
-
+ raise NameError('`%s` has no field called `%s`'% (tokens[0], tokens[1]))
def toUnicode(v, encoding='utf8'):
# if isinstance(v, str):
@@ -407,20 +441,98 @@ def getmoosedoc(tokens, inherited=False):
if len(tokens) > 1:
docstring.write(toUnicode(getfielddoc(tokens)))
else:
- docstring.write(toUnicode('%s\n' % (class_element.docs)))
- append_finfodocs(tokens[0], docstring, indent)
+ doc_text = class_element.docs
+ textArray = textwrap.wrap(doc_text, width = 70, replace_whitespace = False, drop_whitespace = False)
+ text = "\n".join(textArray)
+ # print(text)
+ docstring.write(toUnicode('%s\n' % (text)))
+
+ info = '\t Functions and message destinations \n\t------------------------------------\n\t All \'destFinfo\' can be used as destination Field for moose.connect function.\n'
+ docstring.write(toUnicode('%s\n') % (info))
+ docstring.write(toUnicode(getfielddoc_(tokens, 'destFinfo')))
+ # append_finfodocs(tokens[0], docstring, indent)
if inherited:
mro = eval('_moose.%s' % (tokens[0])).mro()
for class_ in mro[1:]:
if class_ == _moose.melement:
break
- docstring.write(toUnicode(
- '\n\n#Inherited from %s#\n' % (class_.__name__)))
- append_finfodocs(class_.__name__, docstring, indent)
+ docstring.write(toUnicode('\n#Inherited from %s#\n' % (class_.__name__)))
+ temp = []
+ temp.append(str(class_.__name__))
+ docstring.write(toUnicode(getfielddoc_(temp, 'destFinfo')))
+ # append_finfodocs(class_.__name__, docstring, indent)
if class_ == _moose.Neutral: # Neutral is the toplevel moose class
break
- return docstring.getvalue()
+ info = '\t Class attributes \n\t -------------------- \n\t All \'valueFinfo\' are attributes of the class, and can be read and written using \n\t standard get and set Functions (unless explicitely mentioned otherwise). \n\t for example, getVm() and setVm() for access to Vm.\n\t All \'sharedFinfo\'...'
+ docstring.write(toUnicode('%s\n') % (info))
+ docstring.write(toUnicode(getfielddoc_(tokens, 'valueFinfo')))
+ # append_finfodocs(tokens[0], docstring, indent)
+ if inherited:
+ mro = eval('_moose.%s' % (tokens[0])).mro()
+ for class_ in mro[1:]:
+ if class_ == _moose.melement:
+ break
+ docstring.write(toUnicode('\n#Inherited from %s#\n' % (class_.__name__)))
+ temp = []
+ temp.append(str(class_.__name__))
+ docstring.write(toUnicode(getfielddoc_(temp, 'valueFinfo')))
+ # append_finfodocs(class_.__name__, docstring, indent)
+ if class_ == _moose.Neutral: # Neutral is the toplevel moose class
+ break
+
+ info = '\t Source Messages \n\t ----------------------- \n\t All \'srcFinfo\' can be used as sourceField for moose.connect function.'
+ docstring.write(toUnicode('%s\n') % (info))
+ docstring.write(toUnicode(getfielddoc_(tokens, 'srcFinfo')))
+ # append_finfodocs(tokens[0], docstring, indent)
+ if inherited:
+ mro = eval('_moose.%s' % (tokens[0])).mro()
+ for class_ in mro[1:]:
+ if class_ == _moose.melement:
+ break
+ docstring.write(toUnicode('\n#Inherited from %s#\n' % (class_.__name__)))
+ temp = []
+ temp.append(str(class_.__name__))
+ docstring.write(toUnicode(getfielddoc_(temp, 'srcFinfo')))
+ # append_finfodocs(class_.__name__, docstring, indent)
+ if class_ == _moose.Neutral: # Neutral is the toplevel moose class
+ break
+
+ info = '\t LookUp Fields \n\t ---------------- \n\t All \'lookupFinfo\' are ..................................\n'
+ docstring.write(toUnicode('%s\n') % (info))
+ docstring.write(toUnicode(getfielddoc_(tokens, 'lookupFinfo')))
+ # append_finfodocs(tokens[0], docstring, indent)
+ if inherited:
+ mro = eval('_moose.%s' % (tokens[0])).mro()
+ for class_ in mro[1:]:
+ if class_ == _moose.melement:
+ break
+ docstring.write(toUnicode('\n#Inherited from %s#\n' % (class_.__name__)))
+ temp = []
+ temp.append(str(class_.__name__))
+ docstring.write(toUnicode(getfielddoc_(temp, 'lookupFinfo')))
+ # append_finfodocs(class_.__name__, docstring, indent)
+ if class_ == _moose.Neutral: # Neutral is the toplevel moose class
+ break
+
+ info = '\t shared fields \n\t ---------------- \n\t All \'sharedFinfo\' are ..................................\n'
+ docstring.write(toUnicode('%s\n') % (info))
+ docstring.write(toUnicode(getfielddoc_(tokens, 'sharedFinfo')))
+ # append_finfodocs(tokens[0], docstring, indent)
+ if inherited:
+ mro = eval('_moose.%s' % (tokens[0])).mro()
+ for class_ in mro[1:]:
+ if class_ == _moose.melement:
+ break
+ docstring.write(toUnicode('\n#Inherited from %s#\n' % (class_.__name__)))
+ temp = []
+ temp.append(str(class_.__name__))
+ docstring.write(toUnicode(getfielddoc_(temp, 'sharedFinfo')))
+ # append_finfodocs(class_.__name__, docstring, indent)
+ if class_ == _moose.Neutral: # Neutral is the toplevel moose class
+ break
+
+ return docstring.getvalue()
def append_finfodocs(classname, docstring, indent):
"""Append list of finfos in class name to docstring"""
@@ -438,7 +550,6 @@ def append_finfodocs(classname, docstring, indent):
except ValueError:
docstring.write(toUnicode('%sNone\n' % (indent)))
-
# the global pager is set from pydoc even if the user asks for paged
# help once. this is to strike a balance between GENESIS user's
# expectation of control returning to command line after printing the
@@ -502,5 +613,4 @@ def doc(arg, inherited=True, paged=True):
if pager:
pager(text)
else:
- print(text)
-
+ print(text, width=80)