mrgsolve
odeproblem.h
Go to the documentation of this file.
1 // Copyright (C) 2013 - 2019 Metrum Research Group
2 //
3 // This file is part of mrgsolve.
4 //
5 // mrgsolve is free software: you can redistribute it and/or modify it
6 // under the terms of the GNU General Public License as published by
7 // the Free Software Foundation, either version 2 of the License, or
8 // (at your option) any later version.
9 //
10 // mrgsolve is distributed in the hope that it will be useful, but
11 // WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with mrgsolve. If not, see <http://www.gnu.org/licenses/>.
17 
23 #ifndef ODEPROBLEM_H
24 #define ODEPROBLEM_H
25 #include <math.h>
26 #include <vector>
27 #include "RcppInclude.h"
28 #include "mrgsolv.h"
29 #include "datarecord.h"
30 //#include "LSODA.h"
31 
32 
33 //
34 // resim functor comes from mrgsolv.h
35 // so it can get defined in the model
36 //
37 
38 class odeproblem;
39 
41 typedef std::vector<rec_ptr> reclist;
42 
44 typedef std::vector<reclist> recstack;
45 
48 
51 
54 
57 
59 //typedef void main_deriv_func(int* neq, double* t,double* y,double* ydot, void* prob);
60 typedef void (*LSODA_ODE_SYSTEM_TYPE)(double t, double *y, double *dydt,
61  odeproblem* _data);
62 typedef void main_deriv_func(int* neq, double* t,double* y,double* ydot);
63 //typedef void (*ODE_FUNC)(double t, double *y, double *dydt, double*);
64 typedef void (*MRGSOLVE_ODE_FUNC)(int neq, double* t, double* y, double* ydot, std::vector<double>& param);
65 
66 
67 #define MRGSOLVE_GET_PRED_CL (pred[0])
68 #define MRGSOLVE_GET_PRED_VC (pred[1])
69 #define MRGSOLVE_GET_PRED_KA (pred[2])
70 #define MRGSOLVE_GET_PRED_Q (pred[3])
71 #define MRGSOLVE_GET_PRED_VP (pred[4])
72 #define MRGSOLVE_GET_PRED_K10 (pred[0]/pred[1])
73 #define MRGSOLVE_GET_PRED_K12 (pred[3]/pred[1])
74 #define MRGSOLVE_GET_PRED_K21 (pred[3]/pred[4])
75 
76 extern "C"{DL_FUNC tofunptr(SEXP a);}
77 
79 void neg_istate(int istate);
80 
81 template<typename T,typename type2> void tofunptr(T b, type2 a) {
82  b = reinterpret_cast<T>(R_ExternalPtrAddr(a));
83 }
84 
85 void dosimeta(void*);
86 void dosimeps(void*);
87 
88 
94 class odeproblem {
95 
96 public:
97  odeproblem(Rcpp::NumericVector param,Rcpp::NumericVector init,
98  Rcpp::List funs,
99  int n_capture_);
100 
102 
103  void do_init_calc(bool answer) {Do_Init_Calc = answer;}
104  void advance(double tfrom, double tto, LSODA& solver);
105  void call_derivs(double *t, double *y, double *ydot);
106  void init(int pos, double value){Init_value[pos] = value;}
107  double init(int pos){return Init_value[pos];}
108 
109  void init_call(const double& time);
110  void init_call_record(const double& time);
111  void init_derivs(double time);
112 
113  void y_init(int pos, double value);
114  void y_init(Rcpp::NumericVector init);
115  void y_add(const unsigned int pos, const double& value);
116 
117  void table_call();
118  void table_init_call();
119  void config_call();
120 
121  void set_d(rec_ptr this_rec);
122 
123  void omega(Rcpp::NumericMatrix& x);
124  void sigma(Rcpp::NumericMatrix& x);
125 
126  arma::mat mv_omega(int n);
127  arma::mat mv_sigma(int n);
128 
129  void pass_envir(Rcpp::Environment* x){d.envir=reinterpret_cast<void*>(x);};
130 
131  bool CFONSTOP(){return d.CFONSTOP;}
132 
133  const std::vector<double>& param() {return Param;}
134  void param(int pos, double value) {Param[pos] = value;}
135 
136  void rate(unsigned int pos, double value) {R[pos] = value;}
137  double rate(unsigned int pos) {return R[pos];}
138  void rate0(unsigned int pos, double value) {R0[pos] = value;}
139  double rate0(unsigned int pos){return R0[pos];}
140 
141  int rate_count(unsigned int pos){return infusion_count[pos];}
142  void rate_add(unsigned int pos, const double& value);
143  void rate_rm(unsigned int pos, const double& value);
144  void rate_bump(const unsigned int pos, const double& value);
145  void rate_main(rec_ptr rec);
146  void rate_reset();
147 
148  void dur(unsigned int pos, double value) {D[pos] = value;}
149  double dur(unsigned int pos){return D[pos];}
150 
151  void reset_newid(const double id_);
152 
153  void eta(int pos, double value) {d.ETA[pos] = value;}
154  void eps(int pos, double value) {d.EPS[pos] = value;}
155  unsigned short int systemoff(){return d.SYSTEMOFF;}
156 
157  void on(unsigned short int cmt);
158  void off(unsigned short int cmt);
159 
160  int is_on(unsigned int eq_n){return On[eq_n];}
161 
162  void fbio(unsigned int pos, double value) {F.at(pos) = value;}
163  double fbio(unsigned int pos);
164  double alag(int cmt);
165 
166  void time(double time_){d.time = time_;}
167  void newind(unsigned int newind_){d.newind = newind_;}
168  unsigned int newind(){return d.newind;}
169 
170  void advan(int x);
171  int advan(){return Advan;}
172  void advan2(const double& tfrom, const double& tto);
173  void advan4(const double& tfrom, const double& tto);
174 
175  void neta(const int n);
176  void neps(const int n);
177 
178  void nid(int n) {d.nid = n;}
179  void nrow(int n) {d.nrow = n;}
180  void idn(int n) {d.idn = n;}
181  void rown(int n) {d.rown=n;}
182 
183  dvec& get_capture() {return Capture;}
184  double capture(int i) {return Capture[i];}
185 
187  void copy_parin(const Rcpp::List& parin);
188  void copy_funs(const Rcpp::List& funs);
189 
190  bool any_mtime() {return d.mevector.size() > 0;}
191  std::vector<mrgsolve::evdata> mtimes(){return d.mevector;}
192  void clear_mtime(){d.mevector.clear();}
193  void y(const int pos, const double value){Y[pos] = value;}
194  double y(const int pos){return Y[pos];}
195  int istate(){return Istate;}
196  void istate(int value){Istate = value;}
197  void lsoda_init(){Istate=1;}
199  int npar() {return Npar;}
201  int neq() {return Neq;}
203  void tol(double atol, double rtol);
204 
205 
206  std::vector<double> Y;
207  std::vector<double> Ydot;
208  std::vector<double> Yout;
209  std::vector<double> Param;
210  std::vector<double> Capture;
211  double Atol;
212  double Rtol;
213  double ssAtol;
214  double ssRtol; //< relative tolerance when finding steady state
215  int Npar;
216  int Neq;
217  int Istate;
218  bool ss_fixed;
219  int ss_n;
220  bool ss_flag;
221  std::vector<int> Ss_cmt;
222 
223  std::vector<double> R0;
224  std::vector<unsigned int> infusion_count;
225  std::vector<double> R;
226  std::vector<double> D;
227  std::vector<double> Init_value;
228  std::vector<double> Init_dummy;
229  std::vector<double> F;
230  std::vector<double> Alag;
231 
232  std::vector<int> On;
234 
235  int Advan;
236  std::vector<double> a;
237  std::vector<double> alpha;
238 
241 
242  arma::mat Omega;
243  arma::mat Sigma;
244 
245  std::vector<double> pred;
246 
251 
253 };
254 
255 
256 double PolyExp(const double& x,
257  const double& dose,
258  const double& rate,
259  const double& xinf,
260  const double& tau,
261  const bool ss,
262  const std::vector<double>& a,
263  const std::vector<double>& alpha,
264  const int n);
265 
266 Rcpp::List TOUCH_FUNS(const Rcpp::NumericVector& lparam,
267  const Rcpp::NumericVector& linit,
268  const Rcpp::CharacterVector& capture,
269  const Rcpp::List& funs);
270 
271 #endif
RcppInclude.h
odeproblem::rate_count
int rate_count(unsigned int pos)
Definition: odeproblem.h:141
databox::rown
int rown
current output row number
Definition: mrgsolv.h:101
odeproblem::rown
void rown(int n)
sets the currenw data set row number
Definition: odeproblem.h:181
odeproblem::eta
void eta(int pos, double value)
Definition: odeproblem.h:153
odeproblem::systemoff
unsigned short int systemoff()
Definition: odeproblem.h:155
odeproblem::Rtol
double Rtol
relative tolerance used by ODE solver
Definition: odeproblem.h:212
odeproblem::advance
void advance(double tfrom, double tto, LSODA &solver)
Definition: odeproblem.cpp:325
databox::SYSTEMOFF
unsigned short int SYSTEMOFF
flag to stop advancing system for current ID
Definition: mrgsolv.h:94
odeproblem::Istate
int Istate
istate value
Definition: odeproblem.h:217
capture
double capture
Definition: housemodel-mread-header.h:21
odeproblem::R0
std::vector< double > R0
acutal current infusion rate
Definition: odeproblem.h:223
odeproblem::is_on
int is_on(unsigned int eq_n)
Definition: odeproblem.h:160
odeproblem::Sigma
arma::mat Sigma
variance/covariance matrix for within-subject variability
Definition: odeproblem.h:243
odeproblem::alpha
std::vector< double > alpha
used for advan 1/2/3/4 calculation
Definition: odeproblem.h:237
odeproblem::d
databox d
various data passed to model functions
Definition: odeproblem.h:233
main_derivs
main_deriv_func main_derivs
Definition: odeproblem.h:78
recstack
std::vector< reclist > recstack
vector of reclist vectors comprising a data set
Definition: odeproblem.h:44
odeproblem::Config
config_func Config
$PREAMBLE function
Definition: odeproblem.h:250
MRGSOLVE_CONFIG_SIGNATURE
#define MRGSOLVE_CONFIG_SIGNATURE
signature for $PREAMBLE
Definition: mrgsolv.h:129
odeproblem::pred
std::vector< double > pred
brings clearances, volumes, and & for advan 1/2/3/4
Definition: odeproblem.h:245
odeproblem::reset_newid
void reset_newid(const double id_)
Reset odeproblem object for new individual.
Definition: odeproblem.cpp:272
odeproblem::istate
void istate(int value)
Definition: odeproblem.h:196
odeproblem::Neq
int Neq
number of equations
Definition: odeproblem.h:216
odeproblem::param
const std::vector< double > & param()
Definition: odeproblem.h:133
odeproblem::mv_sigma
arma::mat mv_sigma(int n)
Definition: odeproblem.cpp:743
odeproblem::pass_envir
void pass_envir(Rcpp::Environment *x)
Definition: odeproblem.h:129
odeproblem::init_call
void init_call(const double &time)
Definition: odeproblem.cpp:197
odeproblem::advan2
void advan2(const double &tfrom, const double &tto)
Definition: odeproblem.cpp:352
odeproblem::ssRtol
double ssRtol
Definition: odeproblem.h:214
rec_ptr
std::shared_ptr< datarecord > rec_ptr
Definition: datarecord.h:29
odeproblem::Ss_cmt
std::vector< int > Ss_cmt
vector of compartments to consider for SS
Definition: odeproblem.h:221
odeproblem::newind
void newind(unsigned int newind_)
Definition: odeproblem.h:167
odeproblem::do_init_calc
void do_init_calc(bool answer)
Definition: odeproblem.h:103
LSODA_ODE_SYSTEM_TYPE
void(* LSODA_ODE_SYSTEM_TYPE)(double t, double *y, double *dydt, odeproblem *_data)
function to hand off to DLSODA
Definition: odeproblem.h:60
main_deriv_func
void main_deriv_func(int *neq, double *t, double *y, double *ydot)
Definition: odeproblem.h:62
MRGSOLVE_TABLE_SIGNATURE
#define MRGSOLVE_TABLE_SIGNATURE
signature for $TABLE
Definition: mrgsolv.h:121
odeproblem::D
std::vector< double > D
receive user input for infusion duration
Definition: odeproblem.h:226
odeproblem::Yout
std::vector< double > Yout
used to hold Y values during solving
Definition: odeproblem.h:208
odeproblem::simeps
mrgsolve::resim simeps
functor for resimulating epsilons
Definition: odeproblem.h:240
dvec
std::vector< double > dvec
vector of doubles
Definition: mrgsolv.h:114
odeproblem::Table
table_func Table
$TABLE function
Definition: odeproblem.h:249
odeproblem::ssAtol
double ssAtol
absolute tolerance when finding steady state
Definition: odeproblem.h:213
odeproblem::Npar
int Npar
number of parameters
Definition: odeproblem.h:215
odeproblem::off
void off(unsigned short int cmt)
Definition: odeproblem.cpp:317
init_func
void(* init_func)(MRGSOLVE_INIT_SIGNATURE)
$MAIN function
Definition: odeproblem.h:47
odeproblem::mv_omega
arma::mat mv_omega(int n)
Definition: odeproblem.cpp:739
databox::time
double time
current simulation time
Definition: mrgsolv.h:92
odeproblem::table_init_call
void table_init_call()
odeproblem::neps
void neps(const int n)
set number of EPSs in the model
Definition: odeproblem.cpp:126
odeproblem::infusion_count
std::vector< unsigned int > infusion_count
number of active infusions
Definition: odeproblem.h:224
odeproblem::nid
void nid(int n)
sets the number of IDs
Definition: odeproblem.h:178
databox::newind
unsigned int newind
new individual flag
Definition: mrgsolv.h:91
odeproblem::rate_add
void rate_add(unsigned int pos, const double &value)
Definition: odeproblem.cpp:292
odeproblem::R
std::vector< double > R
receive user input for infusion rate
Definition: odeproblem.h:225
odeproblem::Atol
double Atol
absolute tolerance used by ODE solver
Definition: odeproblem.h:211
odeproblem::Do_Init_Calc
bool Do_Init_Calc
Flag regulating whether or not initials are taken from $MAIN
Definition: odeproblem.h:252
mrgsolv.h
odeproblem::npar
int npar()
returns the number of parameters
Definition: odeproblem.h:199
databox::idn
int idn
current ID number
Definition: mrgsolv.h:99
odeproblem::rate0
double rate0(unsigned int pos)
Definition: odeproblem.h:139
odeproblem::Advan
int Advan
simulation mode: 1/2/3/4 (PK models) or 13 (odes)
Definition: odeproblem.h:235
odeproblem::param
void param(int pos, double value)
Definition: odeproblem.h:134
config_func
void(* config_func)(MRGSOLVE_CONFIG_SIGNATURE)
$PREAMBLE function
Definition: odeproblem.h:56
odeproblem::Param
std::vector< double > Param
parameter vector
Definition: odeproblem.h:209
odeproblem::rate0
void rate0(unsigned int pos, double value)
Definition: odeproblem.h:138
odeproblem::rate
void rate(unsigned int pos, double value)
Definition: odeproblem.h:136
table_func
void(* table_func)(MRGSOLVE_TABLE_SIGNATURE)
$TABLE function
Definition: odeproblem.h:50
databox::EPS
std::vector< double > EPS
vector of EPS values
Definition: mrgsolv.h:90
odeproblem::y
double y(const int pos)
Definition: odeproblem.h:194
odeproblem::fbio
void fbio(unsigned int pos, double value)
Definition: odeproblem.h:162
odeproblem::copy_funs
void copy_funs(const Rcpp::List &funs)
Definition: odeproblem.cpp:663
odeproblem::nrow
void nrow(int n)
sets the number of data set rows
Definition: odeproblem.h:179
odeproblem::istate
int istate()
Definition: odeproblem.h:195
odeproblem::mtimes
std::vector< mrgsolve::evdata > mtimes()
Definition: odeproblem.h:191
odeproblem::tol
void tol(double atol, double rtol)
sets the absolute and relative tolerances
Definition: odeproblem.cpp:130
odeproblem::y_add
void y_add(const unsigned int pos, const double &value)
add value to compartment pos
Definition: odeproblem.cpp:158
odeproblem::On
std::vector< int > On
compartment on/off indicator
Definition: odeproblem.h:232
odeproblem::lsoda_init
void lsoda_init()
Definition: odeproblem.h:197
odeproblem::capture
double capture(int i)
Definition: odeproblem.h:184
odeproblem::neq
int neq()
returns the number of state variables
Definition: odeproblem.h:201
odeproblem::y_init
void y_init(int pos, double value)
Definition: odeproblem.cpp:143
databox::CFONSTOP
bool CFONSTOP
carry forward on stop indicator
Definition: mrgsolv.h:102
datarecord.h
PolyExp
double PolyExp(const double &x, const double &dose, const double &rate, const double &xinf, const double &tau, const bool ss, const std::vector< double > &a, const std::vector< double > &alpha, const int n)
Definition: odeproblem.cpp:539
databox::ETA
std::vector< double > ETA
vector of ETA values
Definition: mrgsolv.h:89
reclist
std::vector< rec_ptr > reclist
vector of datarecord objects for one ID
Definition: odeproblem.h:38
databox
member functions mevent and tad come in via housemodel; see inst/base/databox.cpp
Definition: mrgsolv.h:87
odeproblem::Init_dummy
std::vector< double > Init_dummy
initial conditions for user input
Definition: odeproblem.h:228
odeproblem::Ydot
std::vector< double > Ydot
dxdt values
Definition: odeproblem.h:207
odeproblem::advan4
void advan4(const double &tfrom, const double &tto)
Definition: odeproblem.cpp:425
odeproblem::Alag
std::vector< double > Alag
dosing lag time
Definition: odeproblem.h:230
odeproblem::ss_n
int ss_n
Max number of doses during SS advance before warning is issued.
Definition: odeproblem.h:219
databox::mevector
std::vector< mrgsolve::evdata > mevector
a collection of model events to pass back
Definition: mrgsolv.h:107
odeproblem::get_capture
dvec & get_capture()
Definition: odeproblem.h:183
mrgsolve::resim
Resim functor.
Definition: mrgsolv.h:38
odeproblem::set_d
void set_d(rec_ptr this_rec)
Definition: odeproblem.cpp:185
odeproblem::ss_fixed
bool ss_fixed
If true, then no warning is issued if SS not reached in ss_n doses.
Definition: odeproblem.h:218
MRGSOLVE_INIT_SIGNATURE
#define MRGSOLVE_INIT_SIGNATURE
signature for $MAIN
Definition: mrgsolv.h:117
dosimeps
void dosimeps(void *)
Definition: odeproblem.cpp:41
odeproblem::Omega
arma::mat Omega
variance/covariance matrix for between-subject variability
Definition: odeproblem.h:242
odeproblem::on
void on(unsigned short int cmt)
Definition: odeproblem.cpp:313
odeproblem::odeproblem
odeproblem(Rcpp::NumericVector param, Rcpp::NumericVector init, Rcpp::List funs, int n_capture_)
Definition: odeproblem.cpp:49
odeproblem::copy_parin
void copy_parin(const Rcpp::List &parin)
copies items passed in through parin into the odeproblem object
Definition: odeproblem.cpp:647
odeproblem::config_call
void config_call()
Call $PREAMBLE function.
Definition: odeproblem.cpp:232
TOUCH_FUNS
Rcpp::List TOUCH_FUNS(const Rcpp::NumericVector &lparam, const Rcpp::NumericVector &linit, const Rcpp::CharacterVector &capture, const Rcpp::List &funs)
odeproblem::Init_value
std::vector< double > Init_value
initial conditions
Definition: odeproblem.h:227
odeproblem::init_call_record
void init_call_record(const double &time)
Definition: odeproblem.cpp:221
odeproblem::call_derivs
void call_derivs(double *t, double *y, double *ydot)
Definition: odeproblem.cpp:174
dosimeta
void dosimeta(void *)
Definition: odeproblem.cpp:33
odeproblem::Inits
init_func Inits
$MAIN function
Definition: odeproblem.h:248
odeproblem::rate_main
void rate_main(rec_ptr rec)
Definition: odeproblem.cpp:244
odeproblem::neta
void neta(const int n)
set number of ETAs in the model
Definition: odeproblem.cpp:121
odeproblem::rate_bump
void rate_bump(const unsigned int pos, const double &value)
Definition: odeproblem.cpp:297
odeproblem::sigma
void sigma(Rcpp::NumericMatrix &x)
Definition: odeproblem.cpp:735
odeproblem::dur
double dur(unsigned int pos)
Definition: odeproblem.h:149
odeproblem::omega
void omega(Rcpp::NumericMatrix &x)
Definition: odeproblem.cpp:731
databox::nrow
int nrow
number of rows in output data set
Definition: mrgsolv.h:100
LSODA
Definition: LSODA.h:42
odeproblem::y
void y(const int pos, const double value)
Definition: odeproblem.h:193
odeproblem::rate_rm
void rate_rm(unsigned int pos, const double &value)
Definition: odeproblem.cpp:301
odeproblem::CFONSTOP
bool CFONSTOP()
Definition: odeproblem.h:131
tofunptr
DL_FUNC tofunptr(SEXP a)
odeproblem::Derivs
deriv_func Derivs
$ODE function
Definition: odeproblem.h:247
neg_istate
void neg_istate(int istate)
odeproblem::any_mtime
bool any_mtime()
Definition: odeproblem.h:190
odeproblem::alag
double alag(int cmt)
Definition: odeproblem.cpp:114
deriv_func
void(* deriv_func)(MRGSOLVE_ODE_SIGNATURE)
$ODE function
Definition: odeproblem.h:53
odeproblem::~odeproblem
~odeproblem()
Definition: odeproblem.h:101
odeproblem::simeta
mrgsolve::resim simeta
functor for resimulating etas
Definition: odeproblem.h:239
odeproblem::rate_reset
void rate_reset()
Reset all infusion rates.
Definition: odeproblem.cpp:237
MRGSOLVE_ODE_FUNC
void(* MRGSOLVE_ODE_FUNC)(int neq, double *t, double *y, double *ydot, std::vector< double > &param)
Definition: odeproblem.h:64
odeproblem::rate
double rate(unsigned int pos)
Definition: odeproblem.h:137
odeproblem::F
std::vector< double > F
bioavability
Definition: odeproblem.h:229
odeproblem::dur
void dur(unsigned int pos, double value)
Definition: odeproblem.h:148
odeproblem::clear_mtime
void clear_mtime()
Definition: odeproblem.h:192
odeproblem::eps
void eps(int pos, double value)
Definition: odeproblem.h:154
odeproblem::a
std::vector< double > a
used for advan 1/2/3/4 calculations
Definition: odeproblem.h:236
odeproblem::newind
unsigned int newind()
Definition: odeproblem.h:168
odeproblem::idn
void idn(int n)
sets the current ID number
Definition: odeproblem.h:180
MRGSOLVE_ODE_SIGNATURE
#define MRGSOLVE_ODE_SIGNATURE
signature for $ODE
Definition: mrgsolv.h:125
databox::envir
void * envir
model environment
Definition: mrgsolv.h:103
odeproblem
Definition: odeproblem.h:94
odeproblem::ss_flag
bool ss_flag
flag indicating when the system is advancing to SS
Definition: odeproblem.h:220
odeproblem::table_call
void table_call()
Call $TABLE function.
Definition: odeproblem.cpp:227
odeproblem::init
double init(int pos)
Definition: odeproblem.h:107
databox::nid
int nid
number of IDs in the data set
Definition: mrgsolv.h:98
odeproblem::time
void time(double time_)
Definition: odeproblem.h:166
odeproblem::Capture
std::vector< double > Capture
captured data items
Definition: odeproblem.h:210
odeproblem::advan
int advan()
Definition: odeproblem.h:171
odeproblem::init
void init(int pos, double value)
Definition: odeproblem.h:106
odeproblem::Y
std::vector< double > Y
compartment amounts
Definition: odeproblem.h:206
odeproblem::init_derivs
void init_derivs(double time)
Definition: odeproblem.cpp:181