phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
popul_proc.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 // POPUL_PROC class
3 
4 #ifndef _POPUL_PROC_H_
5 #define _POPUL_PROC_H_
6 
7 #include <string>
8 #include <cstring>
9 
10 #include "internal.h"
11 
13 
19 template <class STATE, class PARAMETERS, size_t NEVENT, size_t NDEME = 1>
20 class popul_proc_t {
21 
22 protected:
23  // TYPE DEFINITIONS
24  typedef STATE state_t;
25  typedef PARAMETERS parameters_t;
26  static const size_t nevent = NEVENT;
27  static const size_t ndeme = NDEME;
28 
29 protected:
30  // MEMBER DATA
31  slate_t next; // time of next event
32  size_t event; // mark of next event
33  slate_t current; // current time
34  state_t state; // current state
35  parameters_t params; // model parameters
36 
37 private:
38 
39  void clean (void) {}; // memory cleanup
40 
41 public:
42 
43  // SERIALIZATION
45  size_t bytesize (void) const {
46  return 2*sizeof(slate_t) + sizeof(size_t)
47  + sizeof(state_t) + sizeof(parameters_t);
48  };
50  friend raw_t* operator>> (const popul_proc_t &X, raw_t *o) {
51  slate_t A[2]; A[0] = X.current; A[1] = X.next;
52  memcpy(o,A,sizeof(A)); o += sizeof(A);
53  memcpy(o,&X.event,sizeof(size_t)); o += sizeof(size_t);
54  memcpy(o,&X.state,sizeof(state_t)); o += sizeof(state_t);
55  memcpy(o,&X.params,sizeof(parameters_t)); o += sizeof(parameters_t);
56  return o;
57  };
59  friend raw_t* operator>> (raw_t *o, popul_proc_t &X) {
60  X.clean();
61  slate_t A[2];
62  memcpy(A,o,sizeof(A)); o += sizeof(A);
63  X.current = A[0]; X.next = A[1];
64  memcpy(&X.event,o,sizeof(size_t)); o += sizeof(size_t);
65  memcpy(&X.state,o,sizeof(state_t)); o += sizeof(state_t);
66  memcpy(&X.params,o,sizeof(parameters_t)); o += sizeof(parameters_t);
67  return o;
68  };
69 
70 public:
71  // CONSTRUCTORS, ETC.
74  popul_proc_t (double t0 = 0) {
75  clean();
76  next = current = slate_t(t0);
77  event = 0;
78  };
81  o >> *this;
82  };
84  popul_proc_t (const popul_proc_t & X) {
85  raw_t *o = new raw_t[X.bytesize()];
86  X >> o;
87  o >> *this;
88  delete[] o;
89  };
92  clean();
93  raw_t *o = new raw_t[X.bytesize()];
94  X >> o;
95  o >> *this;
96  delete[] o;
97  return *this;
98  };
100  popul_proc_t (popul_proc_t &&) = delete;
104  ~popul_proc_t (void) {
105  clean();
106  };
107 
108 public:
109 
110  // INFORMATION EXTRACTORS
112  slate_t time (void) const {
113  return current;
114  };
115 
116 public:
117 
118  virtual void valid (void) const {};
119 
120 public:
121 
123  void update_params (double*, int);
125  void update_IVPs (double*, int);
127  double event_rates (double *rate, int n) const;
129  virtual void rinit (void) = 0;
131  virtual void jump (int e) = 0;
133  std::string yaml (std::string tab) const;
134 
135 public:
136 
138  void update_clocks (void) {
139  double rate[nevent];
140  double total_rate = event_rates(rate,nevent);
141  if (total_rate > 0) {
142  next = current+rexp(1/total_rate);
143  } else {
144  next = R_PosInf;
145  }
146  double u = runif(0,total_rate);
147  event = 0;
148  while (u > rate[event] && event < nevent) {
149  if (rate[event] < 0)
150  err("in '%s' (%s line %d): invalid negative rate[%zd]=%lg", // #nocov
151  __func__,__FILE__,__LINE__,event,rate[event]); // #nocov
152  u -= rate[event];
153  event++;
154  }
155  assert(event < nevent);
156  };
157 
160  int play (double tfin) {
161  int count = 0;
162  if (current > tfin)
163  err("cannot simulate backward! (current t=%lg, requested t=%lg)",current,tfin);
164  while (next < tfin) {
165  current = next;
166  jump(event);
167  update_clocks();
168  count++;
169  R_CheckUserInterrupt();
170  }
171  if (next > tfin) current = tfin; // relies on Markov property
172  return count;
173  };
174 
175 };
176 
177 #define PARAM_SET(X) if (!ISNA(p[m])) params.X = p[m]; \
178  m++;
179 #define RATE_CALC(X) total += rate[m++] = (X);
180 #define YAML_PARAM(X) (t + #X + ": " + std::to_string(params.X) + "\n")
181 #define YAML_STATE(X) (t + #X + ": " + std::to_string(state.X) + "\n")
182 
183 #endif
Population process class.
Definition: popul_proc.h:20
virtual void valid(void) const
Definition: popul_proc.h:118
popul_proc_t(const popul_proc_t &X)
copy constructor
Definition: popul_proc.h:84
void update_IVPs(double *, int)
set initial-value parameters
Definition: lbdp.cc:48
static const size_t nevent
Definition: popul_proc.h:26
void update_params(double *, int)
set parameters
Definition: lbdp.cc:39
popul_proc_t(double t0=0)
Definition: popul_proc.h:74
slate_t current
Definition: popul_proc.h:33
popul_proc_t(raw_t *o)
constructor from serialized binary form
Definition: popul_proc.h:80
void update_clocks(void)
updates clock and next event
Definition: popul_proc.h:138
STATE state_t
Definition: popul_proc.h:24
slate_t time(void) const
get current time.
Definition: popul_proc.h:112
virtual void jump(int e)=0
makes a jump
void clean(void)
Definition: popul_proc.h:39
std::string yaml(std::string tab) const
machine/human readable info
Definition: lbdp.cc:26
parameters_t params
Definition: popul_proc.h:35
static const size_t ndeme
Definition: popul_proc.h:27
double event_rates(double *rate, int n) const
compute event rates
Definition: lbdp.cc:55
state_t state
Definition: popul_proc.h:34
~popul_proc_t(void)
destructor
Definition: popul_proc.h:104
popul_proc_t(popul_proc_t &&)=delete
move constructor
virtual void rinit(void)=0
initialize the state
size_t event
Definition: popul_proc.h:32
size_t bytesize(void) const
size of serialized binary form
Definition: popul_proc.h:45
PARAMETERS parameters_t
Definition: popul_proc.h:25
int play(double tfin)
Definition: popul_proc.h:160
friend raw_t * operator>>(const popul_proc_t &X, raw_t *o)
binary serialization
Definition: popul_proc.h:50
slate_t next
Definition: popul_proc.h:31
popul_proc_t & operator=(const popul_proc_t &X)
copy assignment operator
Definition: popul_proc.h:91
Rbyte raw_t
Definition: internal.h:43
#define err(...)
Definition: internal.h:18
double slate_t
Definition: internal.h:44
#define n
Definition: lbdp_pomp.c:8