eventbase.h

Go to the documentation of this file.

#ifndef EVENTBASE_H

#define EVENTBASE_H

#include "state.h"
#include <assert.h>
#include <stdlib.h>
#include <iostream>
#include <cmath>

class GslRandomNumberGenerator;

// To perform a check on solveForRealTimeInterval <-> calculateInternalTimeInterval compatibility
#ifndef NDEBUG
        #define EVENTBASE_CHECK_INVERSE
#endif // NDEBUG

// IMPORTANT: this is only meant for positive times! we use a negative
// event time to indicate that a recalculation is necessary

class EventBase
{
public:
        EventBase();
        virtual ~EventBase();

        // None of these public functions should be used directly in a simulation, they
        // are meant to be used in the implementation of a simulation type 
        
        virtual void fire(State *pState, double t);

        // This calls a virtual function so that a derived class can use another distribution for example,
        // and does not need to limit itself to a poisson process
        // Normally the state isn't needed, but it may come in handy, especially when
        // using a different distribution
        void generateNewInternalTimeDifference(GslRandomNumberGenerator *pRndGen, const State *pState);

        double getEventTime() const                                                             { return m_tEvent; }

        // always calculate with current state
        double solveForRealTimeInterval(const State *pState, double t0);

        // Note that this doesn't need to be called if the propensity hasn't changed. It is
        // for that reason that we also store m_tLastCalc
        void subtractInternalTimeInterval(const State *pState, double t1);

        // May be useful to check for events that can only happen once
        // (e.g. someone dying)
        bool isInitialized() const                                                              { return !(m_Tdiff < 0); }

        bool needsEventTimeCalculation() const                                                  { return (m_tEvent < 0); }
        void setNeedEventTimeCalculation()                                                      { m_tEvent = -12345; }

        bool willBeRemoved() const                                                              { return m_willBeRemoved; }
        void setWillBeRemoved(bool f)                                                           { m_willBeRemoved = f; }
protected:
        virtual double getNewInternalTimeDifference(GslRandomNumberGenerator *pRndGen, const State *pState);

        virtual double calculateInternalTimeInterval(const State *pState, double t0, double dt);

        virtual double solveForRealTimeInterval(const State *pState, double Tdiff, double t0);
private:
        double m_Tdiff;
        double m_tLastCalc;
        double m_tEvent; // we'll also use this as a marker to indicate that recalculation is needed

        bool m_willBeRemoved;
};

inline double EventBase::solveForRealTimeInterval(const State *pState, double t0)
{
        if (!needsEventTimeCalculation())
        {
                double dt = m_tEvent - t0;
                assert(dt >= 0);
                return dt;
        }

        double dt = solveForRealTimeInterval(pState, m_Tdiff, t0);
        assert(dt >= 0);

#ifdef EVENTBASE_CHECK_INVERSE
        double tmpDT = calculateInternalTimeInterval(pState, t0, dt);
        double diff = std::abs(tmpDT - m_Tdiff);
        if (diff > 1e-10)
        {
                std::cerr << "ERROR: mismatch 1 between solveForRealTimeInterval and calculateInternalTimeInterval" << std::endl;
                exit(-1);
        }
#endif // EVENTBASE_CHECK_INVERSE

        m_tEvent = t0 + dt; // automatically marks that we don't need to calculate this again
        assert(m_tEvent >= 0);

        m_tLastCalc = t0;

        return dt;
}

inline void EventBase::subtractInternalTimeInterval(const State *pState, double t1)
{ 
        assert(m_Tdiff >= 0); // Could be the case for simultaneous events
        assert(m_tLastCalc >= 0);

        double dT = calculateInternalTimeInterval(pState, m_tLastCalc, t1 - m_tLastCalc); 

#ifdef EVENTBASE_CHECK_INVERSE
        double dt = t1 - m_tLastCalc;
        double tmpDt = solveForRealTimeInterval(pState, dT, m_tLastCalc);
        double diff = std::abs(tmpDt - dt);
        if (diff > 1e-10)
        {
                std::cerr << "ERROR: mismatch 2 between solveForRealTimeInterval and calculateInternalTimeInterval" << std::endl;
                exit(-1);
        }
#endif // EVENTBASE_CHECK_INVERSE

        // It's possible that due to numerical inaccuracies dT even becomes
        // negative. Detect and clip. 

        assert(dT >= 0); // can be zero apparently
        m_Tdiff -= dT;

#ifndef NDEBUG
        if (!(m_Tdiff > -1e-15))
                std::cerr << "New m_Tdiff is " << m_Tdiff << std::endl;
#endif // NDEBUG
        assert(m_Tdiff > -1e-15); // Error accumulation does happen!
        if (m_Tdiff < 0)
                m_Tdiff = 0;

        setNeedEventTimeCalculation();
}


#endif // EVENTBASE_H