GssaVoxelPools Class Reference

#include <GssaVoxelPools.h>

Inheritance diagram for GssaVoxelPools:

Collaboration diagram for GssaVoxelPools:

Public Member Functions
void	advance (const ProcInfo *p)
void	advance (const ProcInfo *p, const GssaSystem *g)
double	getReacVelocity (unsigned int r, const double *s) const
	GssaVoxelPools ()
vector< unsigned int >	numFire () const
unsigned int	pickReac ()
void	recalcTime (const GssaSystem *g, double currTime)
bool	refreshAtot (const GssaSystem *g)
void	reinit (const GssaSystem *g)
void	setNumReac (unsigned int n)
void	setStoich (const Stoich *stoichPtr)
void	setVolumeAndDependencies (double vol)
void	updateAllRateTerms (const vector< RateTerm * > &rates, unsigned int numCoreRates)
void	updateDependentMathExpn (const GssaSystem *g, unsigned int rindex, double time)
void	updateDependentRates (const vector< unsigned int > &deps, const Stoich *stoich)
void	updateRateTerms (const vector< RateTerm * > &rates, unsigned int numCoreRates, unsigned int index)
void	updateReacVelocities (const GssaSystem *g, const double *s, vector< double > &v) const
void	xferIn (XferInfo &xf, unsigned int voxelIndex, const GssaSystem *g)
virtual	~GssaVoxelPools ()
Public Member Functions inherited from VoxelPoolsBase
void	addProxyTransferIndex (unsigned int comptIndex, unsigned int transferIndex)
void	addProxyVoxy (unsigned int comptIndex, Id comptId, unsigned int voxel)
	Adds the index of a voxel to which proxy data should go. More...
void	backwardReacVolumeFactor (unsigned int i, double volume)
void	filterCrossRateTerms (const vector< Id > &offSolverReacs, const vector< pair< Id, Id > > &offSolverReacCompts)
void	forwardReacVolumeFactor (unsigned int i, double volume)
double	getDiffConst (unsigned int) const
double	getN (unsigned int) const
double	getNinit (unsigned int) const
double	getVolume () const
	Return the volume of the voxel. More...
double	getXreacScaleProducts (unsigned int i) const
double	getXreacScaleSubstrates (unsigned int i) const
bool	hasXfer (unsigned int comptIndex) const
	True when this voxel has data to be transferred. More...
bool	isVoxelJunctionPresent (Id i1, Id i2) const
void	print () const
	Debugging utility. More...
void	reinit ()
void	resetXreacScale (unsigned int size)
void	resizeArrays (unsigned int totNumPools)
	Allocate # of pools. More...
const double *	S () const
void	scaleVolsBufsRates (double ratio, const Stoich *stoichPtr)
void	setDiffConst (unsigned int, double v)
void	setN (unsigned int i, double v)
void	setNinit (unsigned int, double v)
void	setVolume (double vol)
	Just assigns the volume without any cascading to other values. More...
const double *	Sinit () const
unsigned int	size () const
vector< double > &	Svec ()
	Returns a handle to the mol # vector. More...
double *	varS ()
double *	varSinit ()
	VoxelPoolsBase ()
void	xferIn (const vector< unsigned int > &poolIndex, const vector< double > &values, const vector< double > &lastValues, unsigned int voxelIndex)
void	xferInOnlyProxies (const vector< unsigned int > &poolIndex, const vector< double > &values, unsigned int numProxyPools, unsigned int voxelIndex)
void	xferOut (unsigned int voxelIndex, vector< double > &values, const vector< unsigned int > &poolIndex)
	Assembles data values for sending out for x-compt reacs. More...
virtual	~VoxelPoolsBase ()

Private Attributes
double	atot_
vector< unsigned int >	numFire_
moose::RNG< double >	rng_
	RNG. More...
double	t_
	Time at which next event will occur. More...
vector< double >	v_

Additional Inherited Members
Protected Attributes inherited from VoxelPoolsBase
vector< RateTerm * >	rates_
const Stoich *	stoichPtr_

Detailed Description

Definition at line 15 of file GssaVoxelPools.h.

Constructor & Destructor Documentation

GssaVoxelPools::GssaVoxelPools

(

)

Definition at line 50 of file GssaVoxelPools.cpp.

                               :
    VoxelPoolsBase(), t_( 0.0 ), atot_( 0.0 )
{ ; }

GssaVoxelPools::~GssaVoxelPools ( )

virtual

Definition at line 54 of file GssaVoxelPools.cpp.

References VoxelPoolsBase::rates_.

{
    for ( unsigned int i = 0; i < rates_.size(); ++i )
        delete( rates_[i] );
}

Member Function Documentation

void GssaVoxelPools::advance

(

const ProcInfo *

p

)

void GssaVoxelPools::advance	(	const ProcInfo *	p,
		const GssaSystem *	g
	)

Definition at line 173 of file GssaVoxelPools.cpp.

References atot_, ProcInfo::currTime, GssaSystem::dependency, KinSparseMatrix::fireReac(), Stoich::getNumRates(), moose::log(), numFire_, pickReac(), refreshAtot(), rng_, GssaSystem::stoich, VoxelPoolsBase::Svec(), t_, GssaSystem::transposeN, moose::RNG< T >::uniform(), updateDependentRates(), Stoich::updateFuncs(), v_, and VoxelPoolsBase::varS().

{
    double nextt = p->currTime;
    while ( t_ < nextt )
    {
        if ( atot_ <= 0.0 )   // reac system is stuck, will not advance.
        {
            t_ = nextt;
            g->stoich->updateFuncs( varS(), t_ );
            // updateDependentMathExpn( g, 0, t_ );
            return;
        }
        unsigned int rindex = pickReac();
        assert( g->stoich->getNumRates() == v_.size() );
        if ( rindex >= g->stoich->getNumRates() )
        {
            // probably cumulative roundoff error here.
            // Recalculate atot to avoid, and redo.
            if ( !refreshAtot( g ) )   // Stuck state.
            {
                t_ = nextt;
                g->stoich->updateFuncs( varS(), t_ );
                // updateDependentMathExpn( g, 0, t_ );
                return;
            }
            // We had a roundoff error, fixed it, but now need to be sure
            // we only fire a reaction where this is permissible.
            for ( unsigned int i = v_.size(); i > 0; --i )
            {
                if ( fabs( v_[i-1] ) > 0.0 )
                {
                    rindex = i - 1;
                    break;
                }
            }
            assert( rindex < v_.size() );
        }

#if ENABLE_CPP11
        double sign = std::copysign( 1, v_[rindex] );
#else
        double sign = double(v_[rindex] >= 0) - double(0 > v_[rindex] );
#endif

        g->transposeN.fireReac( rindex, Svec(), sign );
        numFire_[rindex]++;

        double r = rng_.uniform();

        while ( r <= 0.0 )
            r = rng_.uniform();

        t_ -= ( 1.0 / atot_ ) * log( r );
        g->stoich->updateFuncs( varS(), t_ );
        // updateDependentMathExpn( g, rindex, t_ );
        updateDependentRates( g->dependency[ rindex ], g->stoich );
    }
}

Here is the call graph for this function:

double GssaVoxelPools::getReacVelocity	(	unsigned int	r,
		const double *	s
	)		const

Definition at line 362 of file GssaVoxelPools.cpp.

References VoxelPoolsBase::rates_.

Referenced by updateDependentRates().

{
    double v = rates_[r]->operator()( s );
    // assert( v >= 0.0 );
    return v;

    // return rates_[r]->operator()( s );
}

Here is the caller graph for this function:

vector< unsigned int > GssaVoxelPools::numFire

(

)

const

Definition at line 298 of file GssaVoxelPools.cpp.

References numFire_.

Referenced by Gsolve::getNumFire().

{
    return numFire_;
}

Here is the caller graph for this function:

unsigned int GssaVoxelPools::pickReac

(

)

Definition at line 104 of file GssaVoxelPools.cpp.

References atot_, rng_, moose::RNG< T >::uniform(), and v_.

Referenced by advance().

{
    double r = rng_.uniform( ) * atot_;
    double sum = 0.0;

    // This is an inefficient way to do it. Can easily get to
    // log time or thereabouts by doing one or two levels of
    // subsidiary tables. Too many levels causes slow-down because
    // of overhead in managing the tree.
    // Slepoy, Thompson and Plimpton 2008
    // report a linear time version.
    for ( vector< double >::const_iterator
            i = v_.begin(); i != v_.end(); ++i )
    {
        if ( r < ( sum += fabs( *i ) ) )
        {
            // double vel = fabs( getReacVelocity( i - v_.begin(), S() ) );
            // assert( vel >= 0 );
            return static_cast< unsigned int >( i - v_.begin() );
        }
    }
    return v_.size();
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GssaVoxelPools::recalcTime	(	const GssaSystem *	g,
		double	currTime
	)

Recalculates the time for the next event. Used when we have a clocked update of rates due to timed functions. In such cases the propensities may change invisibly, so we need to update time estimates

Definition at line 162 of file GssaVoxelPools.cpp.

References atot_, moose::log(), refreshAtot(), rng_, t_, and moose::RNG< T >::uniform().

{
    refreshAtot( g );
    assert( t_ > currTime );
    t_ = currTime;
    double r = rng_.uniform( );
    while( r == 0.0 )
        r = rng_.uniform( );
    t_ -= ( 1.0 / atot_ ) * log( r );
}

Here is the call graph for this function:

bool GssaVoxelPools::refreshAtot

(

const GssaSystem *

g

)

Cleans out all reac rates and recalculates atot. Needed whenever a mol conc changes, or if there is a roundoff error. Returns true if OK, returns false if it is in a stuck state and atot<=0

Definition at line 140 of file GssaVoxelPools.cpp.

References atot_, VoxelPoolsBase::S(), SAFETY_FACTOR, GssaSystem::stoich, t_, Stoich::updateFuncs(), updateReacVelocities(), v_, and VoxelPoolsBase::varS().

Referenced by advance(), recalcTime(), reinit(), and xferIn().

{
    g->stoich->updateFuncs( varS(), t_ );
    updateReacVelocities( g, S(), v_ );
    atot_ = 0;
    for ( vector< double >::const_iterator
            i = v_.begin(); i != v_.end(); ++i )
        atot_ += fabs(*i);
    atot_ *= SAFETY_FACTOR;
    // Check if the system is in a stuck state. If so, terminate.
    if ( atot_ <= 0.0 )
    {
        return false;
    }
    return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GssaVoxelPools::reinit

(

const GssaSystem *

g

)

Builds the gssa system as needed.

Definition at line 232 of file GssaVoxelPools.cpp.

References moose::__rng_seed__, moose::error, Stoich::getNumVarPools(), numFire_, refreshAtot(), VoxelPoolsBase::reinit(), rng_, moose::RNG< T >::setSeed(), GssaSystem::stoich, t_, moose::RNG< T >::uniform(), Stoich::updateFuncs(), GssaSystem::useRandInit, v_, and VoxelPoolsBase::varS().

{
    rng_.setSeed( moose::__rng_seed__ );
    VoxelPoolsBase::reinit(); // Assigns S = Sinit;
    unsigned int numVarPools = g->stoich->getNumVarPools();
    g->stoich->updateFuncs( varS(), 0 );

    double* n = varS();

    if( g->useRandInit )
    {
        vector<double> error(numVarPools, 0.0);
        map<double, vector<Eref>> groupByVal;

        for ( unsigned int i = 0; i < numVarPools; ++i )
        {
            error[i] = n[i];
            double base = std::floor( n[i] );
            assert( base >= 0.0 );
            double frac = n[i] - base;
            if ( rng_.uniform() >= frac )
                n[i] = base;
            else
                n[i] = base + 1.0;

            error[i] -= n[i];

            //if( true )
            //{
            //    //NOTE: Thats how I get the name of the pool at this index.
            //    Eref e = g->stoich->getPoolByIndex( i ).eref();
            //    groupByVal[n[i]].push_back( e );

            //    // Guess the fix the error.
            //}


        }

        double extra = std::accumulate( error.begin(), error.end(), 0.0 );
        if( std::abs(extra) > 0.1 )
        {
            cout << "WARNING: Extra " << extra
                << " molecules in system after converting fractional to integer e.g. 1.1 = 1 (~90% of times) or 2 (~10% of times)." << endl;
        }
    }
    else     // Just round to the nearest int.
    {
        for ( unsigned int i = 0; i < numVarPools; ++i )
        {
#if ENABLE_CPP11
            // Just like rint but does not raise exception.
            // See http://en.cppreference.com/w/cpp/numeric/math/nearbyint for
            // details.
            n[i] = std::nearbyint(n[i]);
#else
            n[i] = round( n[i] );
#endif
        }
    }

    t_ = 0.0;
    refreshAtot( g );
    numFire_.assign( v_.size(), 0 );
}

Here is the call graph for this function:

void GssaVoxelPools::setNumReac

(

unsigned int

n

)

Definition at line 128 of file GssaVoxelPools.cpp.

References numFire_, and v_.

{
    v_.clear();
    v_.resize( n, 0.0 );
    numFire_.resize( n, 0 );
}

void GssaVoxelPools::setStoich

(

const Stoich *

stoichPtr

)

Definition at line 372 of file GssaVoxelPools.cpp.

References VoxelPoolsBase::stoichPtr_.

{
    stoichPtr_ = stoichPtr;
}

void GssaVoxelPools::setVolumeAndDependencies ( double  vol )

virtual

Assign the volume, and handle the cascading effects by scaling all the dependent values of nInit and rates if applicable.

Reimplemented from VoxelPoolsBase.

Definition at line 378 of file GssaVoxelPools.cpp.

References Stoich::getNumCoreRates(), Stoich::getRateTerms(), VoxelPoolsBase::setVolumeAndDependencies(), VoxelPoolsBase::stoichPtr_, and updateAllRateTerms().

{
    VoxelPoolsBase::setVolumeAndDependencies( vol );
    updateAllRateTerms( stoichPtr_->getRateTerms(),
                        stoichPtr_->getNumCoreRates() );
}

Here is the call graph for this function:

void GssaVoxelPools::updateAllRateTerms ( const vector< RateTerm * > &  rates, unsigned int  numCoreRates  )

virtual

Reassign entire rate vector.

Implements VoxelPoolsBase.

Definition at line 307 of file GssaVoxelPools.cpp.

References VoxelPoolsBase::getVolume(), VoxelPoolsBase::getXreacScaleProducts(), VoxelPoolsBase::getXreacScaleSubstrates(), and VoxelPoolsBase::rates_.

Referenced by setVolumeAndDependencies().

{
    for ( unsigned int i = 0; i < rates_.size(); ++i )
        delete( rates_[i] );
    rates_.resize( rates.size() );

    for ( unsigned int i = 0; i < numCoreRates; ++i )
        rates_[i] = rates[i]->copyWithVolScaling( getVolume(), 1, 1 );
    for ( unsigned int i = numCoreRates; i < rates.size(); ++i )
        rates_[i] = rates[i]->copyWithVolScaling( getVolume(),
                    getXreacScaleSubstrates(i - numCoreRates),
                    getXreacScaleProducts(i - numCoreRates ) );
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GssaVoxelPools::updateDependentMathExpn	(	const GssaSystem *	g,
		unsigned int	rindex,
		double	time
	)

Definition at line 64 of file GssaVoxelPools.cpp.

References GssaSystem::stoich, Stoich::updateFuncs(), and VoxelPoolsBase::varS().

{
    // The issue is that if the expression depends on t, we really need
    // to update it every timestep. But then a cascading set of reacs
    // should also be updated.
    // The lower commented block has all funcs updated every time step,
    // but this too
    // doesn't update the cascading reacs. So this is now handled by the
    // useClockedUpdate flag, and we use the lower block here instead.
    /*
    const vector< unsigned int >& deps = g->dependentMathExpn[ rindex ];
    for( vector< unsigned int >::const_iterator
            i = deps.begin(); i != deps.end(); ++i ) {
            g->stoich->funcs( *i )->evalPool( varS(), time );
    }
    */
    /*
    unsigned int numFuncs = g->stoich->getNumFuncs();
    for( unsigned int i = 0; i < numFuncs; ++i )
    {
        g->stoich->funcs( i )->evalPool( varS(), time );
    }
    */
    // This function is equivalent to the loop above.
    g->stoich->updateFuncs( varS(), time );
}

Here is the call graph for this function:

void GssaVoxelPools::updateDependentRates	(	const vector< unsigned int > &	deps,
		const Stoich *	stoich
	)

Definition at line 92 of file GssaVoxelPools.cpp.

References atot_, getReacVelocity(), VoxelPoolsBase::S(), and v_.

Referenced by advance().

{
    for ( vector< unsigned int >::const_iterator
            i = deps.begin(); i != deps.end(); ++i )
    {
        atot_ -= fabs( v_[ *i ] );
        // atot_ += ( v[ *i ] = ( *rates_[ *i ] )( S() );
        atot_ += fabs( v_[ *i ] = getReacVelocity( *i, S() ) );
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GssaVoxelPools::updateRateTerms ( const vector< RateTerm * > &  rates, unsigned int  numCoreRates, unsigned int  index  )

virtual

Update specified index on rate terms. The entire rates vector is passed in for the source values, the index specifies which entry on the local rate vector is to be updated.

Implements VoxelPoolsBase.

Definition at line 322 of file GssaVoxelPools.cpp.

References VoxelPoolsBase::getVolume(), VoxelPoolsBase::getXreacScaleProducts(), VoxelPoolsBase::getXreacScaleSubstrates(), and VoxelPoolsBase::rates_.

{
    // During setup or expansion of the reac system, this might be called
    // before the local rates_ term is assigned. If so, ignore.
    if ( index >= rates_.size() )
        return;
    delete( rates_[index] );
    if ( index >= numCoreRates )
        rates_[index] = rates[index]->copyWithVolScaling(
                            getVolume(),
                            getXreacScaleSubstrates(index - numCoreRates),
                            getXreacScaleProducts(index - numCoreRates ) );
    else
        rates_[index] = rates[index]->copyWithVolScaling(
                            getVolume(), 1.0, 1.0 );
}

Here is the call graph for this function:

void GssaVoxelPools::updateReacVelocities	(	const GssaSystem *	g,
		const double *	s,
		vector< double > &	v
	)		const

updateReacVelocities computes the velocity v of each reaction. This is a utility function for programs like SteadyState that need to analyze velocity.

Definition at line 344 of file GssaVoxelPools.cpp.

References Stoich::getStoichiometryMatrix(), SparseMatrix< T >::nColumns(), VoxelPoolsBase::rates_, and GssaSystem::stoich.

Referenced by refreshAtot().

{
    const KinSparseMatrix& N = g->stoich->getStoichiometryMatrix();
    assert( N.nColumns() == rates_.size() );

    vector< RateTerm* >::const_iterator i;
    v.clear();
    v.resize( rates_.size(), 0.0 );
    vector< double >::iterator j = v.begin();

    for ( i = rates_.begin(); i != rates_.end(); i++)
    {
        *j++ = (**i)( s );
        assert( !std::isnan( *( j-1 ) ) );
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GssaVoxelPools::xferIn	(	XferInfo &	xf,
		unsigned int	voxelIndex,
		const GssaSystem *	g
	)

Digests incoming data values for cross-compt reactions. Sums the changes in the values onto the specified pools.

Definition at line 400 of file GssaVoxelPools.cpp.

References XferInfo::lastValues, refreshAtot(), rng_, XferInfo::subzero, moose::RNG< T >::uniform(), XferInfo::values, VoxelPoolsBase::varS(), and XferInfo::xferPoolIdx.

{
    unsigned int offset = voxelIndex * xf.xferPoolIdx.size();
    vector< double >::const_iterator i = xf.values.begin() + offset;
    vector< double >::const_iterator j = xf.lastValues.begin() + offset;
    vector< double >::iterator m = xf.subzero.begin() + offset;
    double* s = varS();
    bool hasChanged = false;
    for ( vector< unsigned int >::const_iterator
            k = xf.xferPoolIdx.begin(); k != xf.xferPoolIdx.end(); ++k )
    {
        double& x = s[*k];
        // cout << x << "   i = " << *i << *j << "  m = " << *m << endl;
        double dx = *i++ - *j++;
        double base = floor( dx );
        if ( rng_.uniform() >= (dx - base) )
            x += base;
        else
            x += base + 1.0;

        // x += round( *i++ - *j );
        if ( x < *m )
        {
            *m -= x;
            x = 0;
        }
        else
        {
            x -= *m;
            *m = 0;
        }
        /*
        double y = fabs( x - *j );
        // hasChanged |= ( fabs( x - *j ) < 0.1 ); // they are all integers.
        hasChanged |= ( y > 0.1 ); // they are all integers.
        */
        // j++;
        m++;
    }
    // If S has changed, then I should update rates of all affected reacs.
    // Go through stoich matrix to find affected reacs for each mol
    // Store in list, since some may be hit more than once in this func.
    // When done, go through and update all affected ones.
    /*
    if ( hasChanged ) {
        refreshAtot( g );
    }
    */
    // Does this fix the problem of negative concs?
    refreshAtot( g );
}

Here is the call graph for this function:

Member Data Documentation

double GssaVoxelPools::atot_

private

Total propensity of all the reactions in the system

Definition at line 82 of file GssaVoxelPools.h.

Referenced by advance(), pickReac(), recalcTime(), refreshAtot(), and updateDependentRates().

vector< unsigned int > GssaVoxelPools::numFire_

private

Definition at line 92 of file GssaVoxelPools.h.

Referenced by advance(), numFire(), reinit(), and setNumReac().

moose::RNG<double> GssaVoxelPools::rng_

private

RNG.

Definition at line 97 of file GssaVoxelPools.h.

Referenced by advance(), pickReac(), recalcTime(), reinit(), and xferIn().

double GssaVoxelPools::t_

private

Time at which next event will occur.

Definition at line 77 of file GssaVoxelPools.h.

Referenced by advance(), recalcTime(), refreshAtot(), and reinit().

vector< double > GssaVoxelPools::v_

private

State vector of reaction velocities. Only a subset are recalculated on each step.

Definition at line 88 of file GssaVoxelPools.h.

Referenced by advance(), pickReac(), refreshAtot(), reinit(), setNumReac(), and updateDependentRates().

---

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

• moose-core/ksolve/GssaVoxelPools.h
• moose-core/ksolve/GssaVoxelPools.cpp