Merge pull request #601 from nlmixr2/dTheta

D theta
nlmixr2 · Oct 27, 2023 · f95fb27 · f95fb27
2 parents 4e761bb + d98e7be
commit f95fb27
Show file tree

Hide file tree

Showing 12 changed files with 51 additions and 30 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -75,11 +75,12 @@ rxExpandSens2_ <- function(state, s1, s2) {
 #' @param state is the state to expand
 #' @param neta is the number of etas
 #' @param pred type of prediction
+#' @param isTheta logical, is the expansion actually for thetas instead of etas
 #' @keywords internal
 #' @return String of symengine expressions to evaluate to calculate df/deta
 #' @export
-rxExpandFEta_ <- function(state, neta, pred) {
-    .Call(`_rxode2_rxExpandFEta_`, state, neta, pred)
+rxExpandFEta_ <- function(state, neta, pred, isTheta = FALSE) {
+    .Call(`_rxode2_rxExpandFEta_`, state, neta, pred, isTheta)
 }
 
 #' Rep R0 for foce

diff --git a/R/rxode2.R b/R/rxode2.R
@@ -96,7 +96,7 @@ NA_LOGICAL <- NA # nolint
 #' @includeRmd man/rmdhunks/rxode2-create-models.Rmd
 #'
 #' @includeRmd man/rmdhunks/rxode2-syntax-hunk.Rmd
-#' 
+#'
 #' @return An object (environment) of class `rxode2` (see Chambers and Temple Lang (2001))
 #'      consisting of the following list of strings and functions:
 #'
@@ -1275,6 +1275,7 @@ rxCompile <- function(model, dir, prefix, force = FALSE, modName = NULL,
       .cc <- gsub("\n", "", .cc)
       .cflags <- rawToChar(sys::exec_internal(file.path(R.home("bin"), "R"), c("CMD", "config", "CFLAGS"))$stdout)
       .cflags <- gsub("\n", "", .cflags)
+      .cflags <- paste0(.cflags, " -O", getOption("rxode2.compile.O", "2"))
       .shlibCflags <- rawToChar(sys::exec_internal(file.path(R.home("bin"), "R"), c("CMD", "config", "SHLIB_CFLAGS"))$stdout)
       .shlibCflags <- gsub("\n", "", .shlibCflags)
       .cpicflags <- rawToChar(sys::exec_internal(file.path(R.home("bin"), "R"), c("CMD", "config", "CPICFLAGS"))$stdout)

diff --git a/R/symengine.R b/R/symengine.R
@@ -67,7 +67,11 @@ regIfOrElse <- rex::rex(or(regIf, regElse))
   "rxGt" = c("(", ">", ")"),
   "rxLt" = c("(", "<", ")"),
   "rxAnd" = c("(", "&&", ")"),
-  "rxOr" = c("(", "||", ")")
+  "rxOr" = c("(", "||", ")"),
+  "R_pow"=c("(", ")^(", ")"),
+  "R_pow_di"=c("(", ")^(", ")"),
+  "Rx_pow"=c("(", ")^(", ")"),
+  "Rx_pow_di"=c("(", ")^(", ")")
 )
 
 ## atan2

diff --git a/inst/include/rxode2_RcppExports.h b/inst/include/rxode2_RcppExports.h
@@ -89,17 +89,17 @@ namespace rxode2 {
         return Rcpp::as<List >(rcpp_result_gen);
     }
 
-    inline List rxExpandFEta_(CharacterVector state, int neta, int pred) {
-        typedef SEXP(*Ptr_rxExpandFEta_)(SEXP,SEXP,SEXP);
+    inline List rxExpandFEta_(CharacterVector state, int neta, int pred, bool isTheta = false) {
+        typedef SEXP(*Ptr_rxExpandFEta_)(SEXP,SEXP,SEXP,SEXP);
         static Ptr_rxExpandFEta_ p_rxExpandFEta_ = NULL;
         if (p_rxExpandFEta_ == NULL) {
-            validateSignature("List(*rxExpandFEta_)(CharacterVector,int,int)");
+            validateSignature("List(*rxExpandFEta_)(CharacterVector,int,int,bool)");
             p_rxExpandFEta_ = (Ptr_rxExpandFEta_)R_GetCCallable("rxode2", "_rxode2_rxExpandFEta_");
         }
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_rxExpandFEta_(Shield<SEXP>(Rcpp::wrap(state)), Shield<SEXP>(Rcpp::wrap(neta)), Shield<SEXP>(Rcpp::wrap(pred)));
+            rcpp_result_gen = p_rxExpandFEta_(Shield<SEXP>(Rcpp::wrap(state)), Shield<SEXP>(Rcpp::wrap(neta)), Shield<SEXP>(Rcpp::wrap(pred)), Shield<SEXP>(Rcpp::wrap(isTheta)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();

diff --git a/man/reexports.Rd b/man/reexports.Rd
diff --git a/man/rxExpandFEta_.Rd b/man/rxExpandFEta_.Rd
diff --git a/man/rxode2.Rd b/man/rxode2.Rd
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -144,22 +144,23 @@ RcppExport SEXP _rxode2_rxExpandSens2_(SEXP stateSEXP, SEXP s1SEXP, SEXP s2SEXP)
     return rcpp_result_gen;
 }
 // rxExpandFEta_
-List rxExpandFEta_(CharacterVector state, int neta, int pred);
-static SEXP _rxode2_rxExpandFEta__try(SEXP stateSEXP, SEXP netaSEXP, SEXP predSEXP) {
+List rxExpandFEta_(CharacterVector state, int neta, int pred, bool isTheta);
+static SEXP _rxode2_rxExpandFEta__try(SEXP stateSEXP, SEXP netaSEXP, SEXP predSEXP, SEXP isThetaSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::traits::input_parameter< CharacterVector >::type state(stateSEXP);
     Rcpp::traits::input_parameter< int >::type neta(netaSEXP);
     Rcpp::traits::input_parameter< int >::type pred(predSEXP);
-    rcpp_result_gen = Rcpp::wrap(rxExpandFEta_(state, neta, pred));
+    Rcpp::traits::input_parameter< bool >::type isTheta(isThetaSEXP);
+    rcpp_result_gen = Rcpp::wrap(rxExpandFEta_(state, neta, pred, isTheta));
     return rcpp_result_gen;
 END_RCPP_RETURN_ERROR
 }
-RcppExport SEXP _rxode2_rxExpandFEta_(SEXP stateSEXP, SEXP netaSEXP, SEXP predSEXP) {
+RcppExport SEXP _rxode2_rxExpandFEta_(SEXP stateSEXP, SEXP netaSEXP, SEXP predSEXP, SEXP isThetaSEXP) {
     SEXP rcpp_result_gen;
     {
         Rcpp::RNGScope rcpp_rngScope_gen;
-        rcpp_result_gen = PROTECT(_rxode2_rxExpandFEta__try(stateSEXP, netaSEXP, predSEXP));
+        rcpp_result_gen = PROTECT(_rxode2_rxExpandFEta__try(stateSEXP, netaSEXP, predSEXP, isThetaSEXP));
     }
     Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
     if (rcpp_isInterrupt_gen) {
@@ -1948,7 +1949,7 @@ static int _rxode2_RcppExport_validate(const char* sig) {
         signatures.insert("List(*rxExpandGrid_)(RObject&,RObject&,RObject&)");
         signatures.insert("List(*rxExpandSens_)(CharacterVector,CharacterVector)");
         signatures.insert("List(*rxExpandSens2_)(CharacterVector,CharacterVector,CharacterVector)");
-        signatures.insert("List(*rxExpandFEta_)(CharacterVector,int,int)");
+        signatures.insert("List(*rxExpandFEta_)(CharacterVector,int,int,bool)");
         signatures.insert("std::string(*rxRepR0_)(int)");
         signatures.insert("List(*rxExpandNesting)(const RObject&,List&,bool)");
         signatures.insert("bool(*rxIs)(const RObject&,std::string)");

diff --git a/src/expandGrid.cpp b/src/expandGrid.cpp
@@ -248,49 +248,52 @@ List rxExpandSens2_(CharacterVector state, CharacterVector s1, CharacterVector s
 //' @param state is the state to expand
 //' @param neta is the number of etas
 //' @param pred type of prediction
+//' @param isTheta logical, is the expansion actually for thetas instead of etas
 //' @keywords internal
 //' @return String of symengine expressions to evaluate to calculate df/deta
 //' @export
 //[[Rcpp::export]]
-List rxExpandFEta_(CharacterVector state, int neta, int pred){
+List rxExpandFEta_(CharacterVector state, int neta, int pred, bool isTheta=false){
   CharacterVector fe(neta);
   CharacterVector calcS(neta);
   int nstate = state.size();
+	std::string th = "";
+	if (isTheta) th = "TH";
   for (int i = 0; i < neta; i++){
     std::string etaN = std::to_string(i+1);
     std::string feta;
     std::string calc;
 
     switch(pred){
     case 2:
-      feta = "rx__sens_rx_pred__BY_ETA_" +
+      feta = "rx__sens_rx_pred__BY_" + th + "ETA_" +
         etaN + "___";
-      calc = "assign(\"" + feta + "\",with(.s,-D(rx_pred_, ETA_" +
+      calc = "assign(\"" + feta + "\",with(.s,-D(rx_pred_, " + th + "ETA_" +
         etaN + "_)";
       for (int j = nstate; j--;){
-        calc += "-rx__sens_" + as<std::string>(state[j]) + "_BY_ETA_" + etaN +
+        calc += "-rx__sens_" + as<std::string>(state[j]) + "_BY_" + th + "ETA_" + etaN +
           "___*D(rx_pred_,\""+ symengineRes(as<std::string>(state[j])) + "\")";
       }
       calc += "), envir=.s)";
       break;
     case 1:
-      feta = "rx__sens_rx_pred__BY_ETA_" +
+      feta = "rx__sens_rx_pred__BY_" + th + "ETA_" +
         etaN + "___";
-      calc = "assign(\"" + feta + "\",with(.s,D(rx_pred_, ETA_" +
+      calc = "assign(\"" + feta + "\",with(.s,D(rx_pred_, " + th + "ETA_" +
         etaN + "_)";
       for (int j = nstate; j--;){
-        calc += "+rx__sens_" + as<std::string>(state[j]) + "_BY_ETA_" + etaN +
+        calc += "+rx__sens_" + as<std::string>(state[j]) + "_BY_" + th + "ETA_" + etaN +
           "___*D(rx_pred_,"+ symengineRes(as<std::string>(state[j])) + ")";
       }
       calc += "), envir=.s)";
       break;
     case 0:
-      feta = "rx__sens_rx_r__BY_ETA_" +
+      feta = "rx__sens_rx_r__BY_" + th + "ETA_" +
         etaN + "___";
-      calc = "assign(\"" + feta + "\",with(.s,D(rx_r_,ETA_" +
+      calc = "assign(\"" + feta + "\",with(.s,D(rx_r_," + th +"ETA_" +
         etaN + "_)";
       for (int j = nstate; j--;){
-        calc += "+rx__sens_" + as<std::string>(state[j]) + "_BY_ETA_" + etaN +
+        calc += "+rx__sens_" + as<std::string>(state[j]) + "_BY_" + th + "ETA_" + etaN +
           "___*D(rx_r_,"+ symengineRes(as<std::string>(state[j])) + ")";
       }
       calc += "), envir=.s)";

diff --git a/src/init.c b/src/init.c
@@ -168,7 +168,7 @@ SEXP _rxode2_rxAllowUnload(SEXP);
 SEXP _rxode2_rxExpandGrid_(SEXP, SEXP, SEXP);
 SEXP _rxode2_rxExpandSens_(SEXP, SEXP);
 SEXP _rxode2_rxExpandSens2_(SEXP, SEXP, SEXP);
-SEXP _rxode2_rxExpandFEta_(SEXP, SEXP, SEXP);
+SEXP _rxode2_rxExpandFEta_(SEXP, SEXP, SEXP, SEXP);
 SEXP _rxode2_rxRepR0_(SEXP);
 SEXP _rxode2_rLKJ1(SEXP, SEXP, SEXP);
 SEXP _rxode2_rLKJcv1(SEXP, SEXP);
@@ -424,7 +424,7 @@ void R_init_rxode2(DllInfo *info){
     {"_rxode2_rxExpandGrid_", (DL_FUNC) &_rxode2_rxExpandGrid_, 3},
     {"_rxode2_rxExpandSens_", (DL_FUNC) &_rxode2_rxExpandSens_, 2},
     {"_rxode2_rxExpandSens2_",(DL_FUNC) &_rxode2_rxExpandSens2_, 3},
-    {"_rxode2_rxExpandFEta_", (DL_FUNC) &_rxode2_rxExpandFEta_, 3},
+    {"_rxode2_rxExpandFEta_", (DL_FUNC) &_rxode2_rxExpandFEta_, 4},
     {"_rxode2_rxRepR0_", (DL_FUNC) &_rxode2_rxRepR0_, 1},
     {"_rxode2_rxOptRep_", (DL_FUNC) &_rxode2_rxOptRep_, 1},
     {"_rxode2_rxSetSilentErr", (DL_FUNC) &_rxode2_rxSetSilentErr, 1},

diff --git a/src/rxData.cpp b/src/rxData.cpp
@@ -2486,7 +2486,7 @@ LogicalVector rxSolveFree(){
   _globals.gomega = NULL;
   if (_globals.gsigma != NULL) free(_globals.gsigma);
   _globals.gsigma = NULL;
-  
+
   _globals.zeroTheta = false;
   _globals.zeroOmega = false;
   _globals.zeroSigma = false;

diff --git a/tests/testthat/test-dsl.R b/tests/testthat/test-dsl.R
@@ -553,4 +553,12 @@ rxTest({
     expect_equal(rxToSE("tad"), "(t-tlast())")
     expect_equal(rxFromSE("tlast(NaN)"), "tlast()")
   })
+
+  test_that("parsing errors", {
+
+    test <- "E0=THETA[1];\nEm=0.5;\nE50=THETA[2];\ng=2;\nv=E0+Em*t^g/(E50^g+t^g);\nrx_yj_~152;\nrx_lambda_~1;\nrx_low_~0;\nrx_hi_~1;\nrx_r_~0;\nrx_pred_~DV*v-log(1+exp(v));\n"
+
+    expect_error(rxS(test), NA)
+
+  })
 })