odeproblem.h

Go to the documentation of this file.

// Copyright (C) 2013 - 2026  Metrum Research Group
//
// This file is part of mrgsolve.
//
// mrgsolve is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// mrgsolve is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mrgsolve.  If not, see <http://www.gnu.org/licenses/>.

 
#ifndef ODEPROBLEM_H
#define ODEPROBLEM_H
#include <math.h>
#include <vector>
#include <deque>
#include "RcppInclude.h"
#include "mrgsolv.h"
#include "datarecord.h"
 
// 
// resim functor comes from mrgsolv.h
// so it can get defined in the model
// 
 
class odeproblem;

typedef std::deque<rec_ptr> reclist;

typedef std::vector<reclist> recstack;

typedef void (*init_func)(MRGSOLVE_INIT_SIGNATURE);

typedef void (*table_func)(MRGSOLVE_TABLE_SIGNATURE);

typedef void (*event_func)(MRGSOLVE_EVENT_SIGNATURE);

typedef void (*deriv_func)(MRGSOLVE_ODE_SIGNATURE);

typedef void (*config_func)(MRGSOLVE_CONFIG_SIGNATURE);

//typedef void main_deriv_func(int* neq, double* t,double* y,double* ydot, void* prob);
typedef void (*LSODA_ODE_SYSTEM_TYPE)(double t, double *y, double *dydt, 
              odeproblem* _data);
typedef void main_deriv_func(int* neq, double* t,double* y,double* ydot);
//typedef void (*ODE_FUNC)(double t, double *y, double *dydt, double*);
typedef void (*MRGSOLVE_ODE_FUNC)(int neq, double* t, double* y, double* ydot, std::vector<double>& param);
 
 
#define MRGSOLVE_GET_PRED_CL  (pred[0]) 
#define MRGSOLVE_GET_PRED_VC  (pred[1]) 
#define MRGSOLVE_GET_PRED_KA  (pred[2]) 
#define MRGSOLVE_GET_PRED_Q   (pred[3]) 
#define MRGSOLVE_GET_PRED_VP  (pred[4]) 
#define MRGSOLVE_GET_PRED_Q2  (pred[5]) 
#define MRGSOLVE_GET_PRED_VP2 (pred[6]) 
#define MRGSOLVE_GET_PRED_K10 (pred[0]/pred[1]) 
#define MRGSOLVE_GET_PRED_K12 (pred[3]/pred[1]) 
#define MRGSOLVE_GET_PRED_K21 (pred[3]/pred[4]) 
#define MRGSOLVE_GET_PRED_K13 (pred[5]/pred[1]) 
#define MRGSOLVE_GET_PRED_K31 (pred[5]/pred[6]) 
 
extern "C"{DL_FUNC tofunptr(SEXP a);}
 
main_deriv_func main_derivs;
void neg_istate(int istate);
 
template<typename T,typename type2> void tofunptr(T b, type2 a) {
  b = reinterpret_cast<T>(R_ExternalPtrAddr(a));
}
 
void dosimeta(void*);
void dosimeps(void*);
 

class odeproblem {
 
public:
  odeproblem(Rcpp::List param,Rcpp::NumericVector init, 
             Rcpp::List funs,
             int n_capture_);
 
  ~odeproblem(){};
 
  void do_init_calc(bool answer) {Do_Init_Calc = answer;}
  void advance(double tfrom, double tto, LSODA& solver);
  void call_derivs(double *t, double *y, double *ydot);
  void init(int pos, double value){Init_value[pos] = value;}
  double init(int pos){return Init_value[pos];}
  
  void init_call(const double& time);
  void init_call_record(const double& time);
  void init_derivs(double time);
 
  void y_init(int pos, double value);
  void y_init(Rcpp::NumericVector init);
  void y_add(const unsigned int pos, const double& value);
  
  void table_call();
  void table_init_call();
  void event_call();
  void config_call();
  
  void set_d(rec_ptr this_rec);
 
  void omega(const Rcpp::S4& mod);
  void sigma(const Rcpp::S4& mod);
  void copy_sigma_diagonals();
 
  arma::mat mv_omega(int n);
  arma::mat mv_sigma(int n);
 
  void pass_envir(Rcpp::Environment* x){d.envir=reinterpret_cast<void*>(x);};
  
  bool CFONSTOP(){return d.CFONSTOP;}
  
  const std::vector<double>& param() {return Param;}
  void param(int pos, double value) {Param[pos] = value;}
 
  void check_data_rate(rec_ptr rec, int cmtn);
  void check_modeled_dur(rec_ptr rec); 
  void check_modeled_rate(rec_ptr rec);
 
  void rate(unsigned int pos, double value) {R[pos] = value;}
  double rate(unsigned int pos);
  void rate0(unsigned int pos, double value) {R0[pos] = value;}
  double rate0(unsigned int pos);
 
  int rate_count(unsigned int pos){return infusion_count[pos];}
  void rate_add(unsigned int pos, const double& value);
  void rate_rm(unsigned int pos,  const double& value);
  void rate_bump(const unsigned int pos, const double& value);
  void rate_main(rec_ptr rec, int cmtn);
  void rate_reset();
  
  void dur(unsigned int pos, double value) {D[pos] = value;}
  double dur(unsigned int pos);
  
  void reset_newid(const double id_);
  
  void eta(int pos, double value) {d.ETA[pos] = value;}
  void eps(int pos, double value) {d.EPS[pos] = value;}
  unsigned short int systemoff(){return d.SYSTEMOFF;}
  
  void on(unsigned short int cmt);
  void off(unsigned short int cmt);
  
  int is_on(unsigned int eq_n){return On[eq_n];}
  
  void fbio(unsigned int pos, double value) {F.at(pos) = value;}
  double fbio(unsigned int pos);
  double alag(int cmt);
 
  void time(double time_){d.time = time_;}
  void newind(unsigned int newind_){d.newind = newind_;}
  unsigned int newind(){return d.newind;}
 
  void advan(int x);
  int  advan(){return Advan;}
  void advan2(const double& tfrom, const double& tto);
  void advan4(const double& tfrom, const double& tto);
  void advan11(const double& tfrom, const double& tto);
  
  void set_eta();
  void set_eps();
  int neta(){return Omega.n_rows;}
  int neps(){return Sigma.n_rows;}
  
  void nid(int n) {d.nid=n;}
  void idn(int n) {d.idn=n;}
  
  void nrow(int n) {d.nrow=n;}
  void rown(int n) {d.rown=n;}
  
  void inrow(int n) {d.inrow=n;}
  void irown(int n) {d.irown=n;}
  
  dvec& get_capture() {return Capture;}
  double capture(int i) {return Capture[i];}
  
  void copy_parin(const Rcpp::List& parin, const Rcpp::S4& mod);
  void copy_funs(const Rcpp::List& funs);
  
  bool any_mtime() {return d.mevector.size() > 0;}
  std::vector<mrgsolve::evdata> mtimes(){return d.mevector;}
  void clear_mtime(){d.mevector.clear();}
  void    y(const int pos, const double value){Y[pos] = value;}
  double  y(const int pos){return Y[pos];}
  int  istate(){return Istate;}
  void istate(int value){Istate = value;}
  void lsoda_init(){Istate=1;}
  int npar() {return Npar;}
  int neq() {return Neq;}
  void tol(double atol, double rtol);
  
 
  std::vector<double> Y; 
  std::vector<double> Ydot;  
  std::vector<double> Yout; 
  std::vector<double> Param; 
  std::vector<double> Capture; 
  double Atol; 
  double Rtol; 
  double ssAtol; 
  double ssRtol; //<  relative tolerance when finding steady state
  int Npar; 
  int Neq; 
  int Istate; 
  bool ss_fixed; 
  int ss_n; 
  bool ss_flag; 
  bool system_at_ss; 
  std::vector<int> Ss_cmt; 
 
  bool check_modeled_infusions;
  std::vector<double> R0; 
  std::vector<unsigned int> infusion_count; 
  std::vector<double> R; 
  std::vector<double> D; 
  std::vector<double> Init_value; 
  std::vector<double> Init_dummy; 
  std::vector<double> F; 
  std::vector<double> Alag; 
 
  std::vector<int> On; 
  databox d; 
  
  int Advan;  
  std::vector<double> a;     
  std::vector<double> alpha; 
  
  mrgsolve::resim simeta;
  mrgsolve::resim simeps;
 
  arma::mat Omega;
  arma::mat Sigma;
    
  std::vector<double> pred; 
 
  deriv_func Derivs; 
  init_func Inits; 
  table_func Table; 
  event_func Event; 
  config_func Config; 
  
  bool Do_Init_Calc; 
  int interrupt; 
};
 
 
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);
 
#endif