diff --git a/NEWS.md b/NEWS.md index eac8ccd00..3febe2b7e 100644 --- a/NEWS.md +++ b/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. diff --git a/src/codegen.c b/src/codegen.c index e7da220d0..d59e3b275 100644 --- a/src/codegen.c +++ b/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'"); diff --git a/src/genModelVars.h b/src/genModelVars.h index 2395c2634..60efbbc4f 100644 --- a/src/genModelVars.h +++ b/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) && diff --git a/src/parseCmtProperties.h b/src/parseCmtProperties.h index 8532af61a..7bb3b4549 100644 --- a/src/parseCmtProperties.h +++ b/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; diff --git a/src/parseDdt.h b/src/parseDdt.h index 5fd6d1adb..4b5ee62ce 100644 --- a/src/parseDdt.h +++ b/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(); diff --git a/src/parseFunsDosing.h b/src/parseFunsDosing.h index 118e14885..2ad5cfbe8 100644 --- a/src/parseFunsDosing.h +++ b/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; } diff --git a/src/parseVars.h b/src/parseVars.h index 5689be1b5..97f82cf97 100644 --- a/src/parseVars.h +++ b/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; diff --git a/src/tran.c b/src/tran.c index abeef91e6..641d33115 100644 --- a/src/tran.c +++ b/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; diff --git a/src/tran.h b/src/tran.h index 7cb1434bc..38819f336 100644 --- a/src/tran.h +++ b/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__ diff --git a/tests/testthat/test-tad.R b/tests/testthat/test-tad.R index bbb7ebfb9..be8ccd3e0 100644 --- a/tests/testthat/test-tad.R +++ b/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() + + }) })