Skip to content
Snippets Groups Projects
GssaVoxelPools.cpp 14.05 KiB
/**********************************************************************
** This program is part of 'MOOSE', the
** Messaging Object Oriented Simulation Environment.
**           Copyright (C) 2003-2010 Upinder S. Bhalla. and NCBS
** It is made available under the terms of the
** GNU Lesser General Public License version 2.1
** See the file COPYING.LIB for the full notice.
**********************************************************************/

#ifdef ENABLE_CPP11
#include <memory>
#endif

#include "header.h"
#include "RateTerm.h"
#include "FuncTerm.h"
#include "SparseMatrix.h"
#include "KinSparseMatrix.h"
#include "VoxelPoolsBase.h"
#include "../mesh/VoxelJunction.h"
#include "XferInfo.h"
#include "ZombiePoolInterface.h"
#include "Stoich.h"
#include "GssaSystem.h"
#include "GssaVoxelPools.h"
#include "../randnum/RNG.h"
#include "../basecode/global.h"

/**
 * The SAFETY_FACTOR Protects against the total propensity exceeding
 * the cumulative
 * sum of propensities, atot. We do a lot of adding and subtracting of
 * dependency terms from atot. Roundoff error will eventually cause
 * this to drift from the true sum. To guarantee that we never lose
 * the propensity of the last reaction, this safety factor scales the
 * first calculation of atot to be slightly larger. Periodically this
 * will cause the reaction picking step to exceed the last reaction
 * index. This is safe, we just pick another random number.
 * This will happen rather infrequently.
 * That is also a good time to update the cumulative sum.
 * A double should have >15 digits, so cumulative error will be much
 * smaller than this.
 */
const double SAFETY_FACTOR = 1.0 + 1.0e-9;

//////////////////////////////////////////////////////////////
// Class definitions
//////////////////////////////////////////////////////////////

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

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

//////////////////////////////////////////////////////////////
// Solver ops
//////////////////////////////////////////////////////////////

void GssaVoxelPools::updateDependentMathExpn(
    const GssaSystem* g, unsigned int rindex, double time )
{
    // 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 );
}

void GssaVoxelPools::updateDependentRates(
    const vector< unsigned int >& deps, const Stoich* stoich )
{
    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() ) );
    }
}

unsigned int GssaVoxelPools::pickReac()
{
    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();
}

void GssaVoxelPools::setNumReac( unsigned int n )
{
    v_.clear();
    v_.resize( n, 0.0 );
    numFire_.resize( n, 0 );
}

/**
 * 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
 */
bool GssaVoxelPools::refreshAtot( const GssaSystem* g )
{
	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;
}

/**
 * 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
 */
void GssaVoxelPools::recalcTime( const GssaSystem* g, double currTime )
{
    refreshAtot( g );
    assert( t_ > currTime );
    t_ = currTime;
    double r = rng_.uniform( );
    while( r == 0.0 )
        r = rng_.uniform( );
    t_ -= ( 1.0 / atot_ ) * log( r );
}

void GssaVoxelPools::advance( const ProcInfo* p, const GssaSystem* g )
{
    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() );
        }
        double sign = double(v_[rindex] >= 0) - double(0 > v_[rindex] );
        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 );
    }
}

void GssaVoxelPools::reinit( const GssaSystem* g )
{
    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 )
    {
        // round up or down probabilistically depending on fractional
        // num molecules.
        for ( unsigned int i = 0; i < numVarPools; ++i )
        {
            double base = 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;
        }
    }
    else     // Just round to the nearest int.
    {
        for ( unsigned int i = 0; i < numVarPools; ++i )
        {
            n[i] = round( n[i] );
        }
    }
    t_ = 0.0;
    refreshAtot( g );
	numFire_.assign( v_.size(), 0 );
}

vector< unsigned int > GssaVoxelPools::numFire() const
{
	return numFire_;
}

/////////////////////////////////////////////////////////////////////////
// Rate computation functions
/////////////////////////////////////////////////////////////////////////

void GssaVoxelPools::updateAllRateTerms( const vector< RateTerm* >& rates,
        unsigned int numCoreRates )
{
    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 ) );
}

void GssaVoxelPools::updateRateTerms( const vector< RateTerm* >& rates,
                                      unsigned int numCoreRates, unsigned int index 	)
{
    // 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 );
}
/**
 * updateReacVelocities computes the velocity *v* of each reaction.
 * This is a utility function for programs like SteadyState that need
 * to analyze velocity.
 */
void GssaVoxelPools::updateReacVelocities( const GssaSystem* g,
        const double* s, vector< double >& v ) const
{
    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 ) ) );
    }
}

double GssaVoxelPools::getReacVelocity(
    unsigned int r, const double* s ) const
{
    double v = rates_[r]->operator()( s );
    // assert( v >= 0.0 );
    return v;

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

void GssaVoxelPools::setStoich( const Stoich* stoichPtr )
{
    stoichPtr_ = stoichPtr;
}

// Handle volume updates. Inherited virtual func.
void GssaVoxelPools::setVolumeAndDependencies( double vol )
{
    VoxelPoolsBase::setVolumeAndDependencies( vol );
    stoichPtr_->setupCrossSolverReacVols();
    updateAllRateTerms( stoichPtr_->getRateTerms(),
                        stoichPtr_->getNumCoreRates() );
}
//////////////////////////////////////////////////////////////
// Handle cross compartment reactions
//////////////////////////////////////////////////////////////

/*
 * Not sure if we need it. Hold off for now.
static double integralTransfer( double propensity )
{
	double t= floor( propensity );
	if ( rng_.uniform() < propensity - t )
		return t + 1;
	return t;
}
*/

void GssaVoxelPools::xferIn( XferInfo& xf,
                             unsigned int voxelIndex, const GssaSystem* g )
{
    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 );
}

void GssaVoxelPools::xferInOnlyProxies(
    const vector< unsigned int >& poolIndex,
    const vector< double >& values,
    unsigned int numProxyPools,
    unsigned int voxelIndex	)
{
    unsigned int offset = voxelIndex * poolIndex.size();
    vector< double >::const_iterator i = values.begin() + offset;
    unsigned int proxyEndIndex = stoichPtr_->getNumVarPools() +
                                 stoichPtr_->getNumProxyPools();
    for ( vector< unsigned int >::const_iterator
            k = poolIndex.begin(); k != poolIndex.end(); ++k )
    {
        // if ( *k >= size() - numProxyPools )
        if ( *k >= stoichPtr_->getNumVarPools() && *k < proxyEndIndex )
        {
            double base = floor( *i );
            if ( rng_.uniform() >= (*i - base) )
                varSinit()[*k] = (varS()[*k] += base );
            else
                varSinit()[*k] = (varS()[*k] += base + 1.0 );
            // varSinit()[*k] = (varS()[*k] += round( *i ));
        }
        i++;
    }
}