---

ContinuousTimeMonteCarloMover.cpp

Go to the documentation of this file.

// ContinuousTimeMonteCarloMover.cpp: implementation of the ContinuousTimeMonteCarloMover class.
//

#include "ContinuousTimeMonteCarloMover.h"

// Construction/Destruction

CContinuousTimeMonteCarloMover::CContinuousTimeMonteCarloMover()
{

}

CContinuousTimeMonteCarloMover::CContinuousTimeMonteCarloMover(int seed, double frequency)
:CReactionMover(seed, frequency)
{
        m_pRNG = new Rand();
        m_pRNG->seed(seed);
}

CContinuousTimeMonteCarloMover::~CContinuousTimeMonteCarloMover()
{
        delete m_pRNG;
}

void CContinuousTimeMonteCarloMover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        if(this->MoveTimeIsZero(xFinal-xInitial)) {return;}

        // set the member pointer to the local pointer that is passed in
        m_pReactionNetwork = pReactionNetwork;
        m_dTime = xInitial;
        m_iCount = (int) xInitial;
        double nextTime = m_dTime;
        
        double rand1, rand2;
        double dcount = 0.0, scount = 0.0;
        double timeIncrement = 0.0;

        m_iCount++;
        Notify();
        nextTime += 1.0;

        while(m_dTime < xFinal)
        {
                // Choose a move and get the time it will take
                rand1 = m_pRNG->uniform();
                rand2 = m_pRNG->uniform();
                timeIncrement = this->ChooseMove(rand1,rand2);
                m_dTime += timeIncrement;

                // if the next move carries over the current data taking time,
                // take the data first then do the move
                if( m_dTime >= nextTime)
                {
                        // notify                               
                        m_iCount++;
                        Notify();
                        nextTime += 1.0;
                }
                        
                // perform the currently selected move
                this->DoMove();

        }
}

// ChooseMove() should be rewritten to make more sense!

double CContinuousTimeMonteCarloMover::ChooseMove(double rand1, double rand2)
{
        double totalReactionProbability = 0.0;
        int nReactions = m_pReactionNetwork->GetNumberOfReactions();

        double *reactionChannelRate = new double[nReactions];
        double *reactionChannelRateSum = new double[nReactions+1];
        
        // reactionChannelRateSum gives sum up to but not including i
        reactionChannelRateSum[0] = 0.0;
        
        // get a copy of the reaction rates
        std::vector<double> *reactionRates = m_pReactionNetwork->GetReactionRates();

        // compute reaction channel rate sums
        for (int i = 0; i < nReactions; i++)
        {
                reactionChannelRate[i] = reactionRates->at(i);
                reactionChannelRateSum[i+1] = reactionChannelRateSum[i] + reactionChannelRate[i];
                totalReactionProbability += reactionChannelRate[i];
        }

        // choose reaction channel and time to next reaction
        double reactionChannelProbability = rand1*totalReactionProbability;
        // avoid division by zero error
        if(rand2 == 0.0 ) rand2 = 1.0;
        double timeToNextReaction=(1.0/totalReactionProbability)*log(1.0/rand2);
                
        int channelChosen = 0;
        while ((channelChosen < (nReactions - 1))
                                && (reactionChannelProbability > reactionChannelRateSum[channelChosen+1])) 
        {
                channelChosen++;
        }
                
        // set the correct reaction channel
        m_dChannelChosen = channelChosen;

        // free memory

        delete [] reactionChannelRate;
        delete [] reactionChannelRateSum;
        
        return timeToNextReaction;
}

void CContinuousTimeMonteCarloMover::DoMove()
{
        m_pReactionNetwork->GetReaction(m_dChannelChosen)->DoReactionOnce();
}

---