Simpact Cyan
Population based event driven simulation using mNRM
eventbase.h
Go to the documentation of this file.
1 #ifndef EVENTBASE_H
2 
3 #define EVENTBASE_H
4 
9 #include "state.h"
10 #include <assert.h>
11 #include <stdlib.h>
12 #include <iostream>
13 #include <cmath>
14 
16 
17 // To perform a check on solveForRealTimeInterval <-> calculateInternalTimeInterval compatibility
18 #ifndef NDEBUG
19  #ifdef SIMPLEMNRM
20  #define EVENTBASE_CHECK_INVERSE
21  #endif
22 #endif // NDEBUG
23 
24 // IMPORTANT: this is only meant for positive times! we use a negative
25 // event time to indicate that a recalculation is necessary
26 
62 class EventBase
63 {
64 public:
65  EventBase();
66  virtual ~EventBase();
67 
68  // None of these public functions should be used directly in a simulation, they
69  // are meant to be used in the implementation of a simulation type
70 
78  virtual void fire(State *pState, double t);
79 
80  // This calls a virtual function so that a derived class can use another distribution for example,
81  // and does not need to limit itself to a poisson process
82  // Normally the state isn't needed, but it may come in handy, especially when
83  // using a different distribution
84  void generateNewInternalTimeDifference(GslRandomNumberGenerator *pRndGen, const State *pState);
85 
86  double getEventTime() const { return m_tEvent; }
87  double getInternalTimeLeft() const { return m_Tdiff; }
88 
89  // always calculate with current state
90  double solveForRealTimeInterval(const State *pState, double t0);
91 
92  // Note that this doesn't need to be called if the propensity hasn't changed. It is
93  // for that reason that we also store m_tLastCalc
94  void subtractInternalTimeInterval(const State *pState, double t1);
95 
96  // May be useful to check for events that can only happen once
97  // (e.g. someone dying)
98  bool isInitialized() const { return !(m_Tdiff < 0); }
99 
100  bool needsEventTimeCalculation() const { return (m_tEvent < 0); }
101  void setNeedEventTimeCalculation() { m_tEvent = -12345; }
102 
105  bool willBeRemoved() const { return m_willBeRemoved; }
106 
109  void setWillBeRemoved(bool f) { m_willBeRemoved = f; }
110 protected:
123  virtual double getNewInternalTimeDifference(GslRandomNumberGenerator *pRndGen, const State *pState);
124 
134  virtual double calculateInternalTimeInterval(const State *pState, double t0, double dt);
135 
146  virtual double solveForRealTimeInterval(const State *pState, double Tdiff, double t0);
147 private:
148  double m_Tdiff;
149  double m_tLastCalc;
150  double m_tEvent; // we'll also use this as a marker to indicate that recalculation is needed
151 
152  bool m_willBeRemoved;
153 };
154 
155 inline double EventBase::solveForRealTimeInterval(const State *pState, double t0)
156 {
157  if (!needsEventTimeCalculation())
158  {
159  double dt = m_tEvent - t0;
160  assert(dt >= 0);
161  return dt;
162  }
163 
164  double dt = solveForRealTimeInterval(pState, m_Tdiff, t0);
165  assert(dt >= 0);
166 
167 #ifdef EVENTBASE_CHECK_INVERSE
168  double tmpDT = calculateInternalTimeInterval(pState, t0, dt);
169  double diff = std::abs(tmpDT - m_Tdiff);
170  if (diff > 1e-8)
171  {
172  std::cerr << "ERROR: mismatch 1 between solveForRealTimeInterval and calculateInternalTimeInterval" << std::endl;
173  abort();
174  }
175 #endif // EVENTBASE_CHECK_INVERSE
176 
177  m_tEvent = t0 + dt; // automatically marks that we don't need to calculate this again
178  assert(m_tEvent >= 0);
179 
180 #ifdef EVENTBASE_ALWAYS_CHECK_NANTIME
181  if (m_tEvent != m_tEvent)
182  {
183  std::cerr << "ERROR: NaN detected in real world event time calculation" << std::endl;
184  abort();
185  }
186 #endif // EVENTBASE_ALWAYS_CHECK_NANTIME
187 
188  m_tLastCalc = t0;
189 
190  return dt;
191 }
192 
193 inline void EventBase::subtractInternalTimeInterval(const State *pState, double t1)
194 {
195  assert(m_Tdiff >= 0); // Could be the case for simultaneous events
196  assert(m_tLastCalc >= 0);
197 
198  double dT = calculateInternalTimeInterval(pState, m_tLastCalc, t1 - m_tLastCalc);
199 
200 #ifdef EVENTBASE_CHECK_INVERSE
201  double dt = t1 - m_tLastCalc;
202  double tmpDt = solveForRealTimeInterval(pState, dT, m_tLastCalc);
203  double diff = std::abs(tmpDt - dt);
204  if (diff > 1e-8)
205  {
206  std::cerr << "ERROR: mismatch 2 between solveForRealTimeInterval and calculateInternalTimeInterval" << std::endl;
207  abort();
208  }
209 #endif // EVENTBASE_CHECK_INVERSE
210 
211  // It's possible that due to numerical inaccuracies dT even becomes
212  // negative. Detect and clip.
213 
214  assert(dT >= 0); // can be zero apparently
215  m_Tdiff -= dT;
216 
217 #ifdef EVENTBASE_ALWAYS_CHECK_NANTIME
218  if (m_Tdiff != m_Tdiff)
219  {
220  std::cerr << "ERROR: NaN detected in internal event time calculation" << std::endl;
221  abort();
222  }
223 #endif // EVENTBASE_ALWAYS_CHECK_NANTIME
224 
225 #ifndef NDEBUG
226  if (!(m_Tdiff > -1e-15))
227  std::cerr << "New m_Tdiff is " << m_Tdiff << std::endl;
228 #endif // NDEBUG
229  assert(m_Tdiff > -1e-15); // Error accumulation does happen!
230  if (m_Tdiff < 0)
231  m_Tdiff = 0;
232 
233  setNeedEventTimeCalculation();
234 }
235 
236 
237 #endif // EVENTBASE_H
This class both describes the simulation state and contains the core algorithm (as shown on the main ...
Definition: state.h:40
void setWillBeRemoved(bool f)
When an event has fired, by default a new internal fire time will be calculated; setting this flag av...
Definition: eventbase.h:109
bool willBeRemoved() const
Check if the event has been marked for deletion, can avoid a call to the random number generator to c...
Definition: eventbase.h:105
virtual double calculateInternalTimeInterval(const State *pState, double t0, double dt)
This function should map the real world time interval onto an internal time interval and return it...
Definition: eventbase.cpp:42
This class allows you to generate random numbers, and uses the GNU Scientific Library for this...
Definition: gslrandomnumbergenerator.h:16
virtual void fire(State *pState, double t)
This function will be called when the event fires, so this should most likely be re-implemented in yo...
Definition: eventbase.cpp:52
This is the base class for events in the mNRM algorithm.
Definition: eventbase.h:62
virtual double getNewInternalTimeDifference(GslRandomNumberGenerator *pRndGen, const State *pState)
This function will be called to generate a new internal time difference.
Definition: eventbase.cpp:28