doxygen/odeproblem_8h_source.html

// Copyright (C) 2013 - 2019  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 "RcppInclude.h"

#include "mrgsolv.h"

#include "datarecord.h"

//#include "LSODA.h"


//

// resim functor comes from mrgsolv.h

// so it can get defined in the model

//


class odeproblem;


typedef std::vector<rec_ptr> reclist;


typedef std::vector<reclist> recstack;


typedef void (*init_func)(MRGSOLVE_INIT_SIGNATURE);


typedef void (*table_func)(MRGSOLVE_TABLE_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_K10 (pred[0]/pred[1])

#define MRGSOLVE_GET_PRED_K12 (pred[3]/pred[1])

#define MRGSOLVE_GET_PRED_K21 (pred[3]/pred[4])


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::NumericVector 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 config_call();


  void set_d(rec_ptr this_rec);


  void omega(Rcpp::NumericMatrix& x);

  void sigma(Rcpp::NumericMatrix& x);


  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 rate(unsigned int pos, double value) {R[pos] = value;}

  double rate(unsigned int pos) {return R[pos];}

  void rate0(unsigned int pos, double value) {R0[pos] = value;}

  double rate0(unsigned int pos){return R0[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);

  void rate_reset();


  void dur(unsigned int pos, double value) {D[pos] = value;}

  double dur(unsigned int pos){return D[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 neta(const int n);

  void neps(const int n);


  void nid(int n) {d.nid = n;}

  void nrow(int n) {d.nrow = n;}

  void idn(int n) {d.idn = n;}

  void rown(int n) {d.rown=n;}


  dvec& get_capture() {return Capture;}

  double capture(int i) {return Capture[i];}


  void copy_parin(const Rcpp::List& parin);

  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;

  std::vector<int> Ss_cmt;


  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;

  config_func Config;


  bool Do_Init_Calc;

};


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);


Rcpp::List TOUCH_FUNS(const Rcpp::NumericVector& lparam,

                      const Rcpp::NumericVector& linit,

                      const Rcpp::CharacterVector& capture,

                      const Rcpp::List& funs);


#endif