eventbase.h

Go to the documentation of this file.

#ifndef EVENTBASE_H
 
#define EVENTBASE_H
 
#include "algorithm.h"
#include <assert.h>
#include <stdlib.h>
#include <iostream>
#include <cmath>
 
class GslRandomNumberGenerator;
 
// 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(Algorithm *pAlgorithm, 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; }
        double getInternalTimeLeft() const                                                      { return m_Tdiff; }
 
        // 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; }
 
        static bool_t setCheckInverse(bool check);
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;
#ifndef NDEBUG
        static bool s_checkInverse;
#endif // !NDEBUG
};
 
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);
 
#ifndef NDEBUG
        if (s_checkInverse)
        {
                double tmpDT = calculateInternalTimeInterval(pState, t0, dt);
                double diff = std::abs(tmpDT - m_Tdiff);
                if (diff > 1e-8)
                {
                        std::cerr << "ERROR: mismatch 1 between solveForRealTimeInterval and calculateInternalTimeInterval" << std::endl;
                        abort();
                }
        }
#endif // !NDEBUG
 
        m_tEvent = t0 + dt; // automatically marks that we don't need to calculate this again
        assert(m_tEvent >= 0);
 
#ifdef EVENTBASE_ALWAYS_CHECK_NANTIME
        if (m_tEvent != m_tEvent)
                pState->setAbortAlgorithm("NaN detected in internal event time calculation");
#endif // EVENTBASE_ALWAYS_CHECK_NANTIME
 
        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); 
 
#ifndef NDEBUG
        if (s_checkInverse)
        {
                double dt = t1 - m_tLastCalc;
                double tmpDt = solveForRealTimeInterval(pState, dT, m_tLastCalc);
                double diff = std::abs(tmpDt - dt);
                if (diff > 1e-8)
                {
                        std::cerr << "ERROR: mismatch 2 between solveForRealTimeInterval and calculateInternalTimeInterval" << std::endl;
                        abort();
                }
        }
#endif // !NDEBUG
 
        // 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;
 
#ifdef EVENTBASE_ALWAYS_CHECK_NANTIME
        if (m_Tdiff != m_Tdiff)
                pState->setAbortAlgorithm("NaN detected in internal event time calculation");
#endif // EVENTBASE_ALWAYS_CHECK_NANTIME
 
#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