Merge pull request #817 from nlmixr2/815-parsing

Fix parsing for tad() and podo()

NEWS.md

@@ -1,5 +1,10 @@
 # rxode2 (development version)
 
+- Rework the `tad()` and related functions so they use the same
+  interface as compartments (this way they do not depend on the order
+  of compartments); See #815.  Also allow dummy variables to ignore
+  state requirements (for parsing mu-referenced equations).
+
 - Add `getRxNpars` to api.  This allows the development version of
   `babelmixr2` to better check what model is loaded and unload/reload
   as necessary.

src/codegen.c

@@ -676,8 +676,8 @@ SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname,
     if (tb.hasKa == 1) {
       sAppend(&sbOut, "#define _DEPOT_ %d\n", tb.statei);
       sAppend(&sbOut, "#define _CENTRAL_ %d\n", tb.statei+1);
-    } else if (tb.hasCentral == 1) {
-      if (tb.hasDepot){
+    } else if (tb.hasCentralCmt == 1) {
+      if (tb.hasDepotCmt){
         fclose(fpIO);
         _rxode2parse_unprotect();
         err_trans("linCmt() does not have 'depot' compartment without a 'ka'");

src/genModelVars.h

@@ -105,13 +105,13 @@ static inline void calcNparamsNlhsNslhs(void) {
 
 static inline void calcNextra(void) {
   int offCmt=0,nExtra = 0;
-  char *buf, buf2[200];
+  char *buf=NULL, buf2[200];
   for (int i = 0; i < tb.statei; i++){
     if (offCmt == 0 && tb.idu[i] == 0){
+      buf=tb.ss.line[tb.di[i]];
       offCmt = 1;
       nExtra++;
-      buf=tb.ss.line[tb.di[i]];
-    } else if (offCmt == 1 && tb.idu[i] == 1){
+    } else if (offCmt == 1 && tb.idu[i] == 1) {
       // There is an compartment that doesn't have a derivative
       if (tb.linCmt == 0){
         char *v = rc_dup_str(buf, 0);
@@ -200,28 +200,60 @@ static inline SEXP calcIniVals(void) {
 SEXP orderForderS1(SEXP ordIn);
 
 static inline int sortStateVectorsErrHandle(int prop, int pass, int i) {
-  if (prop == 0 || pass == 1) {
+  if (prop == 0 || pass == 1 || tb.dummyLhs == 1) {
     return 1;
   }
+  char *buf = NULL;
+  buf = tb.ss.line[tb.di[i]];
   if ((prop & prop0) != 0) {
-    sAppend(&sbt, "'%s(0)', ", tb.ss.line[tb.di[i]]);
+    sAppend(&sbt, "'%s(0)', ", buf);
   }
   if ((prop & propF) != 0) {
-    sAppend(&sbt, "'f(%s)', ", tb.ss.line[tb.di[i]]);
+    sAppend(&sbt, "'f(%s)', ", buf);
   }
   if ((prop & propAlag) != 0) {
-    sAppend(&sbt, "'alag(%s)', ", tb.ss.line[tb.di[i]]);
+    sAppend(&sbt, "'alag(%s)', ", buf);
   }
   if ((prop & propRate) != 0) {
-    sAppend(&sbt, "'rate(%s)', ", tb.ss.line[tb.di[i]]);
+    sAppend(&sbt, "'rate(%s)', ", buf);
   }
   if ((prop & propDur) != 0) {
-    sAppend(&sbt, "'dur(%s)', ", tb.ss.line[tb.di[i]]);
+    sAppend(&sbt, "'dur(%s)', ", buf);
+  }
+  if ((prop & propTad) != 0) {
+    sAppend(&sbt, "'tad(%s)', ", buf);
+  }
+  if ((prop & propTad0) != 0) {
+    sAppend(&sbt, "'tad0(%s)', ", buf);
+  }
+  if ((prop & propTafd) != 0) {
+    sAppend(&sbt, "'tafd(%s)', ", buf);
+  }
+  if ((prop & propTafd0) != 0) {
+    sAppend(&sbt, "'tafd0(%s)', ", buf);
+  }
+  if ((prop & propTlast) != 0) {
+    sAppend(&sbt, "'tlast(%s)', ", buf);
+  }
+  if ((prop & propTlast0) != 0) {
+    sAppend(&sbt, "'tlast0(%s)', ", buf);
+  }
+  if ((prop & propTfirst) != 0) {
+    sAppend(&sbt, "'tfirst(%s)', ", buf);
+  }
+  if ((prop & propTfirst0) != 0) {
+    sAppend(&sbt, "'tfirst0(%s)', ", buf);
+  }
+  if ((prop & propPodo) != 0) {
+    sAppend(&sbt, "'podo(%s)', ", buf);
+  }
+  if ((prop & propDose) != 0) {
+    sAppend(&sbt, "'dose(%s)', ", buf);
   }
   // Take off trailing "',
   sbt.o -= 2;
   sbt.s[sbt.o] = 0;
-  sAppend(&sbt, " present, but d/dt(%s) not defined\n", tb.ss.line[tb.di[i]]);
+  sAppend(&sbt, " present, but d/dt(%s) not defined\n", buf);
   return 0;
 }
 
@@ -236,7 +268,7 @@ static inline SEXP sortStateVectors(SEXP ordS) {
     int cur = tb.didx[i];
     int prop = tb.dprop[i];
     int pass = 0;
-    if (tb.linCmt){
+    if (tb.linCmt) {
       if (tb.hasDepotCmt == 1 && !strcmp("depot", tb.ss.line[tb.di[i]])){
         pass = 1;
       } else if ((tb.hasCentralCmt == 1 || tb.hasDepotCmt == 1)  &&

src/parseCmtProperties.h

@@ -183,6 +183,9 @@ static inline int handleRemainingAssignmentsCalcPropComplexAssign(nodeInfo ni, c
       // New assignment
       tb.ixL = tb.ix;
       tb.lh[tb.ix] = isLHS;
+      if (!strcmp(v, "rxdummyLhs")) {
+        tb.dummyLhs = 1;
+      }
     } else if (tb.ix < 0){
       if (!strcmp("rxlin___", v)) {
         tb.ixL=-1;
@@ -216,6 +219,9 @@ static inline int handleRemainingAssignmentsCalcPropComplexAssign(nodeInfo ni, c
         tb.lh[tb.ix] = isLHSparam;
       } else {
         tb.lh[tb.ix] = isLHS;
+        if (!strcmp(v, "rxdummyLhs")) {
+          tb.dummyLhs = 1;
+        }
       }
       tb.ixL=-1;
     }
@@ -299,6 +305,9 @@ static inline int handleRemainingAssignmentsCalcPropIni(nodeInfo ni, char *name,
         /* Rprintf("Duplicate %s; %d %d\n", v, tb.lh[tb.ix], tb.ini0[tb.ix]); */
         if (tb.lh[tb.ix] != isLHS){
           tb.lh[tb.ix] = isLHS;
+          if (!strcmp(v, "rxdummyLhs")) {
+            tb.dummyLhs = 1;
+          }
           if (nodeHas(ini0) && tb.ini0[tb.ix] == 1){
             sPrint(&_gbuf,"cannot have conditional initial conditions for '%s'",v);
             updateSyntaxCol();
@@ -331,7 +340,6 @@ static inline int handleRemainingAssignmentsCalcProps(nodeInfo ni, char *name, i
   return 0;
 }
 
-
 static inline int finalizeLineParam(nodeInfo ni, char *name) {
   if (nodeHas(param_statement)) {
     sbDt.o = 0; sbt.o = 0;

src/parseDdt.h

@@ -40,7 +40,7 @@ static inline int new_de(const char *s, int fromWhere) {
 
 static inline int isCmtLhsStatement(nodeInfo ni, char *name, char *v) {
   int hasLhs = 0;
-  if (nodeHas(cmt_statement)){
+  if (nodeHas(cmt_statement)) {
     new_or_ith(v);
     if (tb.lh[tb.ix] || tb.ini[tb.ix]){
       hasLhs=1;
@@ -73,7 +73,10 @@ static inline int add_deCmtProp(nodeInfo ni, char *name, char *v, int hasLhs, in
 }
 
 static inline int add_deState(nodeInfo ni, char *name, char *v, int hasLhs, int fromWhere) {
-  new_or_ith(v);
+  if (new_or_ith(v)) {
+    addSymbolStr(v);
+    new_or_ith(v);
+  }
   if (((tb.ini[tb.ix] == 1 && tb.ini0[tb.ix] == 0) ||
        (tb.lh[tb.ix] == isLHS || tb.lh[tb.ix] == isLHSparam))){
     updateSyntaxCol();

src/parseFunsDosing.h

@@ -57,23 +57,51 @@ static inline int handleFunctionTadSingleStateCcode(transFunctions *tf,char *v2)
 	sAppend(&sb, "_%s1(", tf->v);
 	sAppend(&sbDt, "_%s1(", tf->v);
 	if (new_de(v2, 0)){
-		if (!strcmp("depot", v2)){
-			tb.hasDepot = 1;
-			aAppendN("_DEPOT_)", 8);
-		} else if (!strcmp("central", v2)){
-			tb.hasCentral = 1;
+		if (tb.linCmt && !strcmp("depot", v2)) {
+      tb.hasDepotCmt = 1;
+      aAppendN("_DEPOT_)", 8);
+			return 1;
+		} else if (tb.linCmt && !strcmp("central", v2)) {
+			tb.hasCentralCmt = 1;
 			aAppendN("_CENTRAL_)", 10);
+			return 1;
 		} else {
-			tb.statei++;
-			sAppend(&sb, "%d)", tb.de.n);
-			sAppend(&sbDt, "%d)", tb.de.n);
+      // cannot be lhs statements in tad style assignments
+      // also cannot be from anywhere
+      // temporarily turn off that this is a function
+      // This is not a function
+      int fn = tb.fn;
+      tb.fn = 0;
+      add_de(tf->ni, tf->name, v2, 0, 0);
+      // turn back on that this is a function
+      tb.fn = fn;
 		}
 	} else {
 		new_or_ith(v2);
-		sAppend(&sb, "%d)", tb.id);
-		sAppend(&sbDt, "%d)", tb.id);
 	}
-	// tad(cmt)
+  sAppend(&sb, "__DDT%d__)", tb.id);
+  sAppend(&sbDt, "__DDT%d__)", tb.id);
+	if (tf->isTad && (tb.dprop[tb.id] & propTad) == 0) {
+    tb.dprop[tb.id] += propTad;
+  } else if (tf->isTad0 && (tb.dprop[tb.id] & propTad0) == 0) {
+    tb.dprop[tb.id] += propTad0;
+  } else if (tf->isTafd && (tb.dprop[tb.id] & propTafd) == 0) {
+		tb.dprop[tb.id] += propTafd;
+	} else if (tf->isTafd0 && (tb.dprop[tb.id] & propTafd0) == 0) {
+		tb.dprop[tb.id] += propTafd0;
+	} else if (tf->isTlast && (tb.dprop[tb.id] & propTlast) == 0) {
+		tb.dprop[tb.id] += propTlast;
+	} else if (tf->isTlast0 && (tb.dprop[tb.id] & propTlast0) == 0) {
+		tb.dprop[tb.id] += propTlast0;
+	} else if (tf->isTfirst && (tb.dprop[tb.id] & propTfirst) == 0) {
+		tb.dprop[tb.id] += propTfirst;
+	} else if (tf->isTfirst0 && (tb.dprop[tb.id] & propTfirst0) == 0) {
+		tb.dprop[tb.id] += propTfirst0;
+	} else if (tf->isDose && (tb.dprop[tb.id] & propDose) == 0) {
+		tb.dprop[tb.id] += propDose;
+	} else if (tf->isPodo && (tb.dprop[tb.id] & propPodo) == 0) {
+		tb.dprop[tb.id] += propPodo;
+	}
 	return 1;
 }
 

src/parseVars.h

@@ -115,6 +115,7 @@ static inline int skipReservedVariables(const char *s) {
 /* new symbol? if no, find it's ith */
 static inline int new_or_ith(const char *s) {
   int i;
+  tb.ix=-2;
   if (tb.fn) {tb.ix=-2; return 0;}
   if (!strcmp("lhs", s)){tb.ix=-1; return 0;}
   if (assertForbiddenVariables(s) == 0) return 0;

src/tran.c

@@ -405,8 +405,6 @@ void reset(void) {
   tb.maxtheta   = 0;
   tb.hasCmt     = 0;
   tb.maxeta     = 0;
-  tb.hasDepot   = 0;
-  tb.hasCentral = 0;
   tb.hasDepotCmt = 0;
   tb.hasCentralCmt = 0;
   tb.hasKa      = 0;
@@ -424,6 +422,7 @@ void reset(void) {
   tb.linExtra   = false;
   tb.nwhile     = 0;
   tb.lvlStr     = 0;
+  tb.dummyLhs   = 0;
   tb.nInd       = 0;
   tb.simflg     = 0;
   tb.nLlik      = 0;

src/tran.h

@@ -84,8 +84,6 @@ lhs symbols?
   int maxtheta;
   int hasCmt;
   int maxeta;
-  int hasDepot;
-  int hasCentral;
   int hasDepotCmt;
   int hasCentralCmt;
   int hasKa;
@@ -109,6 +107,7 @@ lhs symbols?
   int lastDdt;
   int nLlik;
   int lvlStr;
+  int dummyLhs;
 } symtab;
 
 extern symtab tb;
@@ -369,10 +368,25 @@ void _rxode2parse_unprotect(void);
 
 char *getLine (char *src, int line, int *lloc);
 
-#define prop0    1
-#define propF    2
-#define propAlag 4
-#define propRate 8
-#define propDur  16
+#define prop0     1
+#define propF     2
+#define propAlag  4
+#define propRate  8
+#define propDur   16
+// tad
+#define propTad   32
+#define propTad0  64
+// tafd
+#define propTafd  128
+#define propTafd0 256
+// tlast
+#define propTlast 512
+#define propTlast0 1024
+// tfirst
+#define propTfirst 2048
+#define propTfirst0 4096
+// podo and dose
+#define propPodo 8192
+#define propDose 16384
 
 #endif // __TRAN_H__

tests/testthat/test-tad.R

@@ -519,4 +519,60 @@ rxTest({
     expect_false(isTRUE(all.equal(x$tad, x$tade)))
 
   })
+
+  test_that("tad parsing", {
+
+    mod2 <- function() {
+      ini({
+        ## Table 3 from Savic 2007
+        cl  <- 17.2 # (L/hr)
+        vc  <- 45.1 # L
+        ka  <- 0.38 # 1/hr
+        mtt <- 1.37 # hr
+        f2 <-0.5    # Fraction of 1st Order portion
+        n   <- 20.1
+      })
+      model({
+        k           <- cl/vc
+        bio <- 1-f2
+        ktr = (n+1)/mtt
+        ## note that lgammafn is the same as lgamma in R.
+        d/dt(depot1) = exp(log(bio*podo(depot))+
+                             log(ktr)+n*log(ktr*tad(depot))-
+                             ktr*tad(depot)-lgammafn(n+1))-ka*depot1
+        d/dt(depot2) <- -ka*depot2
+        f(depot2) <-f2
+        d/dt(cen)   <- ka*depot1 + ka*depot2-k*cen
+      })
+    }
+
+    expect_error(mod2())
+
+    mod2 <- function() {
+      ini({
+        ## Table 3 from Savic 2007
+        cl  <- 17.2 # (L/hr)
+        vc  <- 45.1 # L
+        ka  <- 0.38 # 1/hr
+        mtt <- 1.37 # hr
+        f2 <-0.5    # Fraction of 1st Order portion
+        n   <- 20.1
+      })
+      model({
+        k           <- cl/vc
+        bio <- 1-f2
+        ktr = (n+1)/mtt
+        ## note that lgammafn is the same as lgamma in R.
+        d/dt(depot1) = exp(log(bio*podo(depot1))+
+                             log(ktr)+n*log(ktr*tad(depot1))-
+                             ktr*tad(depot1)-lgammafn(n+1))-ka*depot1
+        d/dt(depot2) <- -ka*depot2
+        f(depot2) <-f2
+        d/dt(cen)   <- ka*depot1 + ka*depot2-k*cen
+      })
+    }
+
+    mod2 <- mod2()
+
+  })
 })