Microsimulation API
/home/marcle/src/R/microsimulation/src/calibperson-r.cc
Go to the documentation of this file.
1 
32 #include "microsimulation.h"
33 #include <Rcpp.h>
34 
35 namespace {
36 
37 using namespace ssim;
38 using namespace std;
39 
41 enum stage_t {DiseaseFree,Precursor,PreClinical,Clinical,Death};
42 
44 enum event_t {toPrecursor, toPreClinical, toClinical, toDeath, Count};
45 
47 string stage_names[5] = {"DiseaseFree","Precursor","PreClinical","Clinical","Death"};
48 
50  Rng * rng;
51 
53 class CalibPerson : public cProcess
54 {
55 public:
56  stage_t stage;
57  bool diseasepot;
58  double Lam1,sigm1,p2,lam2,mu3,tau3;
59  double clinTime;
60  int id;
61 
62  // static member(s)
63  static std::map<std::string, std::vector<double> > report;
64 
65  static void resetPopulation ();
66 
67  CalibPerson() {} // default constructor
68 
69  CalibPerson(double *par, int i=0) {
70  Lam1=par[0];
71  sigm1=par[1];
72  p2=par[2];
73  lam2=par[3];
74  mu3=par[4];
75  tau3=par[5];
76  id=i;
77  stage=DiseaseFree;
78  };
79 
80  void init();
81  virtual void handleMessage(const cMessage* msg);
82  virtual Time age() { return now(); }
83 };
84 
85 void CalibPerson::resetPopulation() {
86  report.clear();
87 }
88 
89 // initialise static member(s)
90 std::map<std::string, std::vector<double> > CalibPerson::report;
91 
95 void CalibPerson::init() {
96  if (R::runif(0,1)<p2){
97  diseasepot = true;}
98  else {
99  diseasepot = false;
100  }
101  clinTime=1000;
102  stage=DiseaseFree;
103  double lam1 = exp(R::rnorm(Lam1,sigm1));
104  scheduleAt(R::rexp(lam1), toPrecursor);
105  double x = R::runif(0,1);
106  scheduleAt((65 - 15*log(-log(x))), toDeath); //Gumbel
107  for(int i=10;i<110;i=i+10){
108  scheduleAt(i,Count);
109  }
110 }
111 
115 void CalibPerson::handleMessage(const cMessage* msg) {
116 
117  double ctime[] = {20,40,60,80};
118  int cind;
119 
120 
121  if (msg->kind == toDeath) {
122  stage=Death;
123  clinTime=std::min(clinTime,now());
124 
125  for(unsigned int i=0; i<4 ; i++){
126  if(i < report["TimeAtRisk"].size()){
127  report["TimeAtRisk"][i] += std::min(ctime[i],clinTime);
128  }
129  else {
130  report["TimeAtRisk"].push_back(std::min(ctime[i],clinTime));
131  }
132 
133  if(clinTime < ctime[i]){
134  break;
135  }
136 
137  }
138 
140  }
141 
142  else if (msg->kind == toPrecursor) {
143  stage = Precursor;
144  if (diseasepot){
145  simtime_t dwellTime = now()+ R::rexp(lam2);
146  scheduleAt(dwellTime, toPreClinical);
147  }
148  }
149 
150  else if (msg->kind == toPreClinical) {
151  stage=PreClinical;
152  simtime_t dwellTime = now()+ exp(R::rnorm(mu3,tau3*mu3));
153  scheduleAt(dwellTime, toClinical);
154  }
155 
156  else if (msg->kind == toClinical) {
157  stage=Clinical;
158  clinTime = now();
159  string stagestr = stage_names[stage];
160  }
161 
162  else if (msg->kind == Count){
163  cind = min(9,int(now()/10 - 1));
164  string stagestr = stage_names[stage];
165 
166  if(report.find(stagestr) == report.end()){ //key not found
167  report[stagestr].assign(10,0);
168  }
169  report[stagestr][cind]+=1;
170  }
171 }
172 
173 
174 extern "C" {
175 
176  RcppExport SEXP callCalibrationSimulation(SEXP parms) {
177  Rcpp::List parmsl(parms);
178  int nin = Rcpp::as<int>(parmsl["n"]);
179  std::vector<double> par = Rcpp::as<std::vector<double> >(parmsl["runpar"]);
180 
181  CalibPerson::resetPopulation();
182  rng = new Rng();
183  rng->set();
184 
185  CalibPerson::report.insert(make_pair("TimeAtRisk", std::vector<double>()));
186 
187  CalibPerson person;
188  for (int i = 0; i < nin; i++) {
189  person = CalibPerson(&par[0],0);
190  rng->nextSubstream();
191  Sim::create_process(&person);
193  Sim::clear();
194  }
195 
196  delete rng;
197 
199 
200  }
201 
202 }
203 
204 }
ssim::cMessage
cMessage class for OMNET++ API compatibility. This provides a heavier message class than Sim::Event,...
Definition: microsimulation.h:211
callCalibrationSimulation
SEXP callCalibrationSimulation(SEXP)
ssim::cProcess
cProcess class for OMNET++ API compatibility.
Definition: microsimulation.h:314
Rcpp::wrap
SEXP wrap(const std::vector< std::pair< T1, T2 > > v)
ssim::simtime_t
Time simtime_t
simtime_t typedef for OMNET++ API compatibility
Definition: microsimulation.h:439
stage_t
stage_t
enum of type of disease stage
Definition: person-r-20121231.cc:40
ssim::Sim::create_process
static ProcessId create_process(Process *)
creates a new process
Definition: ssim.cc:108
ssim::Sim::stop_simulation
static void stop_simulation()
stops execution of the simulation
Definition: ssim.cc:251
ssim::now
Time now()
now() function for compatibility with C++SIM
Definition: microsimulation.cc:9
Death
@ Death
Definition: person-r-20121231.cc:41
ssim::Sim::clear
static void clear()
clears out internal data structures
Definition: ssim.cc:115
ssim::Sim::run_simulation
static void run_simulation()
starts execution of the simulation
Definition: ssim.cc:148
ssim::Time
double Time
virtual time type
Definition: ssim.h:75
ssim
name space for the Siena simulator.
Definition: microsimulation.cc:3
illnessDeath::report
EventReport< short, short, double > report
Definition: illness-death.cpp:13
illnessDeath::event_t
event_t
Definition: illness-death.cpp:11
ssim::Rng
Definition: microsimulation.h:548
microsimulation.h
ssim::Rng::set
void set()
Definition: microsimulation.cc:26
ssim::cMessage::kind
short kind
Definition: microsimulation.h:222
ssim::Rng::nextSubstream
void nextSubstream()
Definition: microsimulation.h:560