---

MixedReactionMover.cpp

Go to the documentation of this file.

// MixedReactionMover.cpp: implementation of the CMixedReactionMover class.
//

#include "MixedReactionMover.h"

// Construction/Destruction

CMixedReactionMover::CMixedReactionMover()
{

}

CMixedReactionMover::CMixedReactionMover(int seed, CDifferentialEquationMover *pODEMover, double absTol)
:CReactionMover(seed,1.0)
{
        m_pODEMover = pODEMover;
        m_pRNG = new Rand(seed);
        m_dAbsTol = absTol;
}

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

void CMixedReactionMover::Move(double xInitial, double xFinal, ReactionNetwork *pReactionNetwork)
{
        // downcast to a CMixedReactionNetwork and make sure it is one - you shouldn't
        // pass a non-mixed RN into this mover anyway!
        m_pReactionNetwork = pReactionNetwork;
        m_pMixedReactionNetwork = dynamic_cast<CMixedReactionNetwork *>(pReactionNetwork);
        if(0 == m_pReactionNetwork)
        {
                throw std::runtime_error("Reaction Network is not a mixed network!");
        }
        // bind the MixedNetwork to the TauNetwork
        m_pTauNetwork = new CTauNetwork(m_pMixedReactionNetwork);

        
        int nChem = m_pMixedReactionNetwork->GetNumberOfChemicals();
        int nRates = m_pMixedReactionNetwork->GetNumberOfReactions();
        m_dTime = xInitial;
        double nextTime = m_dTime;
        // may be unsafe if xInitial has a fractional part!
        m_iCount = (int) xInitial;
        double minRatio;
        double trialEndTime = xFinal;

        double *tempChem = new double[nChem];

        // Notify observers at initial time (usually 0)
        Notify();
        nextTime += 1.0;

        while(m_dTime < xFinal)
        {
                // pick a random number
                double R = m_pRNG->uniform();
                // compute total stochastic reaction probability
                double gammaS = m_pMixedReactionNetwork->GetTotalStochasticReactionProbability();
                while(gammaS == 0.0 && m_dTime < nextTime)
                {
                        // XXX debug
                        cout << "GammaS = " << gammaS << endl;
                        cin.get();
                        // compute the putative integration endpoint
                        minRatio = 10000000;
                        for(int i = 0; i < nRates; i++)
                        {
                                double rate = m_pMixedReactionNetwork->GetReaction(i)->GetRate();
                                double tempRatio = m_dAbsTol/rate;
                                minRatio = __min(tempRatio,minRatio);
                        }
                        trialEndTime = __min(m_dTime + m_pODEMover->GetStepSize(), nextTime);
                        trialEndTime = __min(trialEndTime,m_dTime + minRatio);
                        // XXX debug
                        cout << "Trial end time = " << trialEndTime << endl;
                        // now run for this amount of time on the REAL network (all stochastic
                        // reactions have zero rate so there is no danger of the stochastic 
                        // rates changing anything)
                        m_pODEMover->Move(m_dTime,trialEndTime,m_pReactionNetwork);
                        m_pODEMover->ResetStepSize();
                        // now recompute R
                        gammaS = m_pMixedReactionNetwork->GetTotalStochasticReactionProbability();
                        // XXX debug
                        cout << "New GammaS = " << gammaS << endl;
                        double deltaR = 0.5*gammaS*(trialEndTime - m_dTime);
                        // do I need to repeat this until we're below the tolerance?
                        // XXX debug
                        cout << "DeltaR = " << deltaR << ", R*tol = " << R*m_dAbsTol << endl;
                        R = R - deltaR;
                        // XXX debug
                        cout << "New R = " << R << endl;
                        // XXX debug
                        if(R < 0) {cout << "R = " << R << ", Integrator might fail . . ." << endl;}
                        m_dTime = trialEndTime;
                } // end "gammaS == 0.0" while loop
                if(gammaS > 0.0)
                {
                        // call the tau network from 0 to R, get tOfTau at end.  Save current 
                        // chemical concentrations first.
                        for(int i = 0; i < nChem; i++)
                        {
                                tempChem[i] = m_pReactionNetwork->GetChemical(i)->GetAmount();
                        }
                        double deltaT = m_pODEMover->GetStepSize();
                        m_pODEMover->SetStepSize(gammaS*deltaT);
                        m_pODEMover->Move(0,R,m_pTauNetwork);
                        m_pODEMover->ResetStepSize();
                        // get the amount of time that elapsed
                        double tOfTau = m_pTauNetwork->GetTOfTau();
                        cout << "Tau network integrated, time elapsed = " << tOfTau << endl;
                        if(m_dTime + tOfTau < nextTime)
                        {
                                m_dTime += tOfTau;
                                // pick a stochastic move and if OK do it
                                this->DoStochasticMove(R);
                        }
                        else  // inefficient if nextTime << 1/gammaS
                        {
                                // restore old chemical concentrations before move
                                for(int i = 0; i < nChem; i++)
                                {
                                        m_pReactionNetwork->GetChemical(i)->SetAmount(tempChem[i]);
                                }
                                // run on REAL MIXED network (don't want to integrate 
                                // stochastic reactions)
                                m_pODEMover->Move(m_dTime,nextTime,m_pMixedReactionNetwork);
                                m_pODEMover->ResetStepSize();
                                m_dTime = nextTime;
                                // take data
                                // XXX debug
                                cout << "Taking data at time " << m_dTime << endl;
                                m_iCount++;
                                Notify();
                                nextTime += 1.0;
                        }
                }
                else
                {
                        // XXX debug
                        cout << "Taking data at time " << m_dTime << endl;
                        m_iCount++;
                        Notify();
                        nextTime += 1.0;
                }
        
        } // end "while(time < finalTime)" loop

        delete [] tempChem;
        delete m_pTauNetwork;

        return;
}

void CMixedReactionMover::DoStochasticMove(double dRand)
{
        int i;
        double totalRxnProb = 0.0;
        int nReactions = m_pMixedReactionNetwork->GetNumberOfStochasticReactions();
        
        double *rxnChannelRate = new double[nReactions];
        double *rxnChannelRateSum = new double[nReactions+1];

        // rxnChannelRateSum gives sum up to but not including i
        rxnChannelRateSum[0] = 0.0;
        
        // get a copy of the stochastic reaction rates 
        std::vector<double> *pRxnRates = m_pMixedReactionNetwork->GetStochasticReactionRates();

        // compute reaction channel rate sums
        for (i = 0; i < nReactions; i++)
        {
                rxnChannelRate[i] = pRxnRates->at(i);
                rxnChannelRateSum[i+1] = rxnChannelRateSum[i] + rxnChannelRate[i];
                totalRxnProb += rxnChannelRate[i];
        }

        // choose reaction channel
        double rxnChannelProb = dRand*totalRxnProb;
        
        int chosen = 0;
        while ((chosen < (nReactions - 1)) && (rxnChannelProb > rxnChannelRateSum[chosen+1])) 
        {
                chosen++;
        }
        // if the reaction would change the chemical by an amount less than that present
        // do not do the reaction; otherwise it's OK
        std::vector<Chemical *> *pRxnChem = m_pMixedReactionNetwork->GetStochasticReaction(chosen)->GetChemicals();
        for(i = 0; i < pRxnChem->size(); i++)
        {
                int whichChem = pRxnChem->at(i)->GetChemicalNumber();
                double currAmount = pRxnChem->at(i)->GetAmount();
                double deltaAmount = m_pMixedReactionNetwork->GetStochasticReaction(chosen)->GetNumberChangedByReaction(whichChem);
                if(currAmount + deltaAmount < 0.0)
                {
                        // fail the reaction and return
                        delete [] rxnChannelRate;
                        delete [] rxnChannelRateSum;
                        return;
                }
        }
        m_pMixedReactionNetwork->GetStochasticReaction(chosen)->DoReactionOnce();

        delete [] rxnChannelRate;
        delete [] rxnChannelRateSum;
}

---