From 0d79ce4f119a1c8674a5c9b8adf864b4e8bbb53a Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 10 Sep 2024 15:17:08 -0500 Subject: [PATCH 1/4] Broken cmt ordering in solve --- tests/testthat/test-cmt-order.R | 124 ++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) diff --git a/tests/testthat/test-cmt-order.R b/tests/testthat/test-cmt-order.R index bcfc75e96..9d7d56d4c 100644 --- a/tests/testthat/test-cmt-order.R +++ b/tests/testthat/test-cmt-order.R @@ -554,4 +554,128 @@ rxTest({ expect_equal(tmp$stateExtra, "cp") }) + test_that("compartment ordering does not affect c-code generation #785", + { + + f1 <- function() { + ini({ + TVCL <- c(0, 0.0308628127403366) + TVV1 <- c(0, 0.132584062223269) + TVQ <- c(0, 0.0107718034646091) + TVV2 <- c(0, 0.0613517707801907) + TVR0 <- c(0, 3.94111989610973) + TVKON <- c(0, 1.31095109796496) + TVKOFF <- c(0, 0.0758072298659602) + TVKDEG <- c(0, 1.07290072658165) + TVKINT <- c(0, 3.49422759360926) + ADD.ERR.PK <- c(0, 0.446231462752061) + ETA.CL ~ 0.0499998244205031 + }) + model({ + CL <- TVCL * exp(ETA.CL) + V1 <- TVV1 + Q <- TVQ + V2 <- TVV2 + R0 <- TVR0 + KON <- TVKON + KOFF <- TVKOFF + KDEG <- TVKDEG + KINT <- TVKINT + KEL <- CL/V1 + K12 <- Q/V1 + K21 <- Q/V2 + KD <- KOFF/KON + KSYN <- R0 * KDEG + A3(0) <- R0 + KEL <- 0 + d/dt(A1) = -(K12 + KEL) * A1 + K21 * A2 - KON * A1 * + A3 + KOFF * A4 * V1 + d/dt(A2) = K12 * A1 - K21 * A2 + d/dt(A3) = KSYN - KDEG * A3 - KON * (A1/V1) * A3 + KOFF * + A4 + d/dt(A4) = KON * (A1/V1) * A3 - (KINT + KOFF) * A4 + CFREE = A1/V1 + CTOT = CFREE + A4 + CP = log(CTOT) + CP ~ add(ADD.ERR.PK) + SPEC.OUT <- SPEC + DGRP.OUT <- DGRP + BW.OUT <- BW + }) + } + + f1 <- f1 %>% zeroRe("omega") + + library(rxode2) + + f2 <- function() { + ini({ + TVCL <- c(0, 0.0308628127403366) + TVV1 <- c(0, 0.132584062223269) + TVQ <- c(0, 0.0107718034646091) + TVV2 <- c(0, 0.0613517707801907) + TVR0 <- c(0, 3.94111989610973) + TVKON <- c(0, 1.31095109796496) + TVKOFF <- c(0, 0.0758072298659602) + TVKDEG <- c(0, 1.07290072658165) + TVKINT <- c(0, 3.49422759360926) + ADD.ERR.PK <- c(0, 0.446231462752061) + ETA.CL ~ 0.0499998244205031 + }) + model({ + cmt(A1) + cmt(A2) + cmt(A3) + cmt(A4) + CL <- TVCL * exp(ETA.CL) + V1 <- TVV1 + Q <- TVQ + V2 <- TVV2 + R0 <- TVR0 + KON <- TVKON + KOFF <- TVKOFF + KDEG <- TVKDEG + KINT <- TVKINT + KEL <- CL/V1 + K12 <- Q/V1 + K21 <- Q/V2 + KD <- KOFF/KON + KSYN <- R0 * KDEG + A3(0) <- R0 + KEL <- 0 + d/dt(A1) = -(K12 + KEL) * A1 + K21 * A2 - KON * A1 * + A3 + KOFF * A4 * V1 + d/dt(A2) = K12 * A1 - K21 * A2 + d/dt(A3) = KSYN - KDEG * A3 - KON * (A1/V1) * A3 + KOFF * + A4 + d/dt(A4) = KON * (A1/V1) * A3 - (KINT + KOFF) * A4 + CFREE = A1/V1 + CTOT = CFREE + A4 + CP = log(CTOT) + CP ~ add(ADD.ERR.PK) + SPEC.OUT <- SPEC + DGRP.OUT <- DGRP + BW.OUT <- BW + }) + } + + f2 <- f2 %>% zeroRe("omega") + + ev <- eventTable() + mybw <- 70 + mymw <- 150000 + mydgrp <- 50 + ev$add.dosing(dose=mydgrp*mybw*1000000/mymw ,nbr.doses= 13, dosing.interval= 7) + ev$add.sampling(seq(0,13*7,length.out=7)) + + ev2 <- ev %>% as.data.frame() %>% + dplyr::mutate(BW=mybw, DGRP=mydgrp, SPEC=4) + + s2 <- rxSolve(f2, ev2, returnType="data.frame") + + s1 <- rxSolve(f1, ev2, returnType="data.frame") + + expect_equal(s1, s2) + }) + }) From cbfc04a4810e312455de5e91ffdff6b48570807f Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 10 Sep 2024 15:22:41 -0500 Subject: [PATCH 2/4] Remove unneded call to rxode2 --- tests/testthat/test-cmt-order.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/testthat/test-cmt-order.R b/tests/testthat/test-cmt-order.R index 9d7d56d4c..bef0e7a29 100644 --- a/tests/testthat/test-cmt-order.R +++ b/tests/testthat/test-cmt-order.R @@ -606,8 +606,6 @@ rxTest({ f1 <- f1 %>% zeroRe("omega") - library(rxode2) - f2 <- function() { ini({ TVCL <- c(0, 0.0308628127403366) From 5dd2b5d13af4ce8a99e401ba7790270f22a7fddb Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 10 Sep 2024 16:34:30 -0500 Subject: [PATCH 3/4] define state order in model --- inst/include/rxode2parse_control.h | 5 +++-- src/codegen.c | 20 +++++++++++++++----- src/genModelVars.c | 9 ++++++--- src/parseCmtProperties.h | 16 ++++++++-------- src/parseDdt.h | 10 +++++----- src/rxData.cpp | 18 ++++++++++++------ tests/testthat/test-cmt-order.R | 4 +++- 7 files changed, 52 insertions(+), 30 deletions(-) diff --git a/inst/include/rxode2parse_control.h b/inst/include/rxode2parse_control.h index 5609c0702..9566d5bbc 100644 --- a/inst/include/rxode2parse_control.h +++ b/inst/include/rxode2parse_control.h @@ -127,8 +127,9 @@ #define RxMv_stateProp 24 #define RxMv_sensProp 25 #define RxMv_normProp 26 -#define RxMv_timeId 27 -#define RxMv_md5 28 +#define RxMv_stateOrd 27 +#define RxMv_timeId 28 +#define RxMv_md5 29 #define RxMvFlag_ncmt 0 #define RxMvFlag_ka 1 #define RxMvFlag_linB 2 diff --git a/src/codegen.c b/src/codegen.c index e19d4df7f..9ffb77515 100644 --- a/src/codegen.c +++ b/src/codegen.c @@ -348,7 +348,7 @@ void codegen(char *model, int show_ode, const char *prefix, const char *libname, // stateExtra sAppendN(&sbOut, " ", 2); doDot(&sbOut, buf); - sAppend(&sbOut, " = __zzStateVar__[%d]*((double)(_ON[%d]));\n", i, i); + sAppend(&sbOut, " = __zzStateVar__[__DDT%d__]*((double)(_ON[__DDT%d__]));\n", i, i); } } else { break; @@ -485,7 +485,7 @@ void codegen(char *model, int show_ode, const char *prefix, const char *libname, for (i = 0; i < tb.de.n; i++) { if (tb.idu[i]) { buf=tb.ss.line[tb.di[i]]; - sAppend(&sbOut, " __zzStateVar__[%d]=((double)(_ON[%d]))*(",i,i); + sAppend(&sbOut, " __zzStateVar__[__DDT%d__]=((double)(_ON[__DDT%d__]))*(",i,i); doDot(&sbOut, buf); sAppendN(&sbOut, ");\n", 3); } @@ -552,7 +552,7 @@ extern SEXP _goodFuns; SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname, SEXP pMd5, SEXP timeId, SEXP mvLast, - SEXP goodFuns){ + SEXP goodFuns) { _goodFuns = PROTECT(goodFuns); _rxode2parse_protected++; if (!sbPm.o || !sbNrm.o){ _rxode2parse_unprotect(); @@ -568,6 +568,7 @@ SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname, } fpIO = fopen(CHAR(STRING_ELT(c_file,0)), "wb"); err_msg((intptr_t) fpIO, "error opening output c file\n", -2); + if (badMd5){ SET_STRING_ELT(VECTOR_ELT(mvLast, RxMv_md5), 0, mkChar("")); } else { @@ -622,11 +623,20 @@ SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname, SET_STRING_ELT(trans, 21, mkChar(buf.s)); // IndF } sPrint(&_mv, "%s", CHAR(STRING_ELT(PROTECT(_rxode2_rxQs(mvLast)), 0))); pro++; - UNPROTECT(pro); sFree(&buf); //SET_STRING_ELT(tran, 0, mkChar()); sFree(&sbOut); sIniTo(&sbOut, (int)((sbPm.sN)*5.3)); + + SEXP stateOrd = PROTECT(VECTOR_ELT(mvLast, RxMv_stateOrd)); pro++; + SEXP stateOrdNames = PROTECT(Rf_getAttrib(stateOrd, R_NamesSymbol)); pro++; + int *stateOrdInt = INTEGER(stateOrd); + sAppend(&sbOut, "// Define translation state order for %d states\n", Rf_length(stateOrd)); + for (int i = 0; i < Rf_length(stateOrd); i++){ + sAppend(&sbOut, "#define __DDT%d__ %d // %s\n", stateOrdInt[i]-1, i, + CHAR(STRING_ELT(stateOrdNames, i))); + } + UNPROTECT(pro); // show_ode = 1 dydt // show_ode = 2 Jacobian // show_ode = 3 Ini statement @@ -670,8 +680,8 @@ SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname, } sAppend(&sbOut, "#define _CENTRAL_ %d\n", tb.statei); } - writeSb(&sbOut, fpIO); } + writeSb(&sbOut, fpIO); gCode(1); // d/dt() gCode(2); // jac gCode(3); // ini() diff --git a/src/genModelVars.c b/src/genModelVars.c index 26520e74e..d565ffcdd 100644 --- a/src/genModelVars.c +++ b/src/genModelVars.c @@ -9,8 +9,8 @@ SEXP generateModelVars(void) { calcNextra(); int pro = 0; - SEXP lst = PROTECT(Rf_allocVector(VECSXP, 27));pro++; - SEXP names = PROTECT(Rf_allocVector(STRSXP, 27));pro++; + SEXP lst = PROTECT(Rf_allocVector(VECSXP, 28));pro++; + SEXP names = PROTECT(Rf_allocVector(STRSXP, 28));pro++; SEXP sNeedSort = PROTECT(Rf_allocVector(INTSXP,1));pro++; int *iNeedSort = INTEGER(sNeedSort); @@ -44,7 +44,6 @@ SEXP generateModelVars(void) { trans_syntax_error_report_fn0(_gbuf.s); return R_NilValue; } - populateStateVectors(state, sens, normState, stateRm, extraState, stateProp, sensProp, normProp, INTEGER(ordF)); SEXP dfdy = PROTECT(Rf_allocVector(STRSXP,tb.ndfdy));pro++; populateDfdy(dfdy); @@ -292,6 +291,10 @@ SEXP generateModelVars(void) { SET_VECTOR_ELT(lst, 26, normProp); SET_STRING_ELT(names, 26, mkChar("normProp")); + Rf_setAttrib(ordF, R_NamesSymbol, state); + SET_VECTOR_ELT(lst, 27, ordF); + SET_STRING_ELT(names, 27, mkChar("stateOrd")); + Rf_setAttrib(tran, R_NamesSymbol, trann); Rf_setAttrib(lst, R_NamesSymbol, names); Rf_setAttrib(model, R_NamesSymbol, modeln); diff --git a/src/parseCmtProperties.h b/src/parseCmtProperties.h index 768e7c71e..8532af61a 100644 --- a/src/parseCmtProperties.h +++ b/src/parseCmtProperties.h @@ -4,8 +4,8 @@ static inline int handleCmtPropertyFbio(nodeInfo ni, char *name, char *v) { if ((tb.dprop[tb.id] & propF) == 0) { tb.dprop[tb.id] += propF; } - sAppend(&sb, "_f[%d] = ", tb.id); - sAppend(&sbDt, "_f[%d] = ", tb.id); + sAppend(&sb, "_f[__DDT%d__] = ", tb.id); + sAppend(&sbDt, "_f[__DDT%d__] = ", tb.id); sAppend(&sbt, "f(%s)=", v); tb.curPropN=tb.id; if (foundF == 0) needSort+=1;// & 1 when F @@ -22,8 +22,8 @@ static inline int handleCmtPropertyAlag(nodeInfo ni, char *name, char *v) { if ((tb.dprop[tb.id] & propAlag) == 0) { tb.dprop[tb.id] += propAlag; } - sAppend(&sb, "_alag[%d] = ", tb.id); - sAppend(&sbDt, "_alag[%d] = ", tb.id); + sAppend(&sb, "_alag[__DDT%d__] = ", tb.id); + sAppend(&sbDt, "_alag[__DDT%d__] = ", tb.id); sAppend(&sbt, "alag(%s)=", v); tb.curPropN=tb.id; if (foundLag == 0) needSort+=2; // & 2 when alag @@ -52,8 +52,8 @@ static inline int handleCmtPropertyDur(nodeInfo ni, char *name, char *v) { tb.dprop[tb.id] += propDur; } sb.o=0;sbDt.o=0; sbt.o=0; - sAppend(&sb, "_dur[%d] = ", tb.id); - sAppend(&sbDt, "_dur[%d] = ", tb.id); + sAppend(&sb, "_dur[__DDT%d__] = ", tb.id); + sAppend(&sbDt, "_dur[__DDT%d__] = ", tb.id); sAppend(&sbt, "dur(%s)=", v); tb.curPropN=tb.id; if (foundDur == 0) needSort+=4;// & 4 when dur @@ -70,8 +70,8 @@ static inline int handleCmtPropertyRate(nodeInfo ni, char *name, char *v) { if ((tb.dprop[tb.id] & propRate) == 0) { tb.dprop[tb.id] += propRate; } - sAppend(&sb, "_rate[%d] = ", tb.id); - sAppend(&sbDt, "_rate[%d] = ", tb.id); + sAppend(&sb, "_rate[__DDT%d__] = ", tb.id); + sAppend(&sbDt, "_rate[__DDT%d__] = ", tb.id); sAppend(&sbt, "rate(%s)=", v); tb.curPropN=tb.id; if (foundRate == 0) needSort+=8;// & 8 when rate diff --git a/src/parseDdt.h b/src/parseDdt.h index b854de47b..5fd6d1adb 100644 --- a/src/parseDdt.h +++ b/src/parseDdt.h @@ -134,11 +134,11 @@ static inline int handleDdtAssign(nodeInfo ni, char *name, int i, D_ParseNode *p /* printf("de[%d]->%s[%d]\n",tb.id,v,tb.ix); */ sb.o =0; sbDt.o =0; if (tb.idu[tb.id] == 0){ - sAppend(&sb, "__DDtStateVar__[%d] = ((double)(_ON[%d]))*(_IR[%d] ", tb.id, tb.id, tb.id); - sAppend(&sbDt, "__DDtStateVar_%d__ = ((double)(_ON[%d]))*(_IR[%d] ", tb.id, tb.id, tb.id); + sAppend(&sb, "__DDtStateVar__[__DDT%d__] = ((double)(_ON[__DDT%d__]))*(_IR[__DDT%d__] ", tb.id, tb.id, tb.id); + sAppend(&sbDt, "__DDtStateVar_%d__ = ((double)(_ON[__DDT%d__]))*(_IR[__DDT%d__] ", tb.id, tb.id, tb.id); } else { - sAppend(&sb, "__DDtStateVar__[%d] = ((double)(_ON[%d]))*(", tb.id, tb.id); - sAppend(&sbDt, "__DDtStateVar_%d__ = ((double)(_ON[%d]))*(", tb.id, tb.id); + sAppend(&sb, "__DDtStateVar__[__DDT%d__] = ((double)(_ON[__DDT%d__]))*(", tb.id, tb.id); + sAppend(&sbDt, "__DDtStateVar_%d__ = ((double)(_ON[__DDT%d__]))*(", tb.id, tb.id); } tb.idu[tb.id]=1; aType(TDDT); @@ -214,7 +214,7 @@ static inline int handleDdtRhs(nodeInfo ni, char *name, D_ParseNode *xpn) { sAppend(&sb, "__DDtStateVar_%d__", tb.id); sAppend(&sbDt, "__DDtStateVar_%d__", tb.id); } else { - sAppend(&sb, "__DDtStateVar__[%d]", tb.id); + sAppend(&sb, "__DDtStateVar__[__DDT%d__]", tb.id); sAppend(&sbDt, "__DDtStateVar_%d__", tb.id); aType(TDDT); } diff --git a/src/rxData.cpp b/src/rxData.cpp index 6e33eaa1f..10d68a989 100644 --- a/src/rxData.cpp +++ b/src/rxData.cpp @@ -738,8 +738,8 @@ List rxModelVars_rxode2(const RObject &obj){ //' //' @noRd List rxModelVars_blank() { - List ret(29); - CharacterVector retN(29); + List ret(30); + CharacterVector retN(30); ret[0] = CharacterVector::create(); // params retN[0] = "params"; ret[1] = CharacterVector::create(); // lhs @@ -830,11 +830,17 @@ List rxModelVars_blank() { ret[26] = normProp; // normProp retN[26] = "normProp"; - ret[27] = IntegerVector::create(0); // timeId - retN[27] = "timeId"; + IntegerVector stateOrd = IntegerVector::create(); + normProp.attr("names") = CharacterVector::create(); + ret[27] = stateOrd; + retN[27] = "stateOrd"; + + ret[28] = IntegerVector::create(0); // timeId + retN[28] = "timeId"; + + ret[29] =CharacterVector::create(_["file_md5"] = "", _["parsed_md5"] = ""); // md5 + retN[29] = "md5"; - ret[28] =CharacterVector::create(_["file_md5"] = "", _["parsed_md5"] = ""); // md5 - retN[28] = "md5"; ret.attr("names") = retN; ret.attr("class") = "rxModelVars"; return ret; diff --git a/tests/testthat/test-cmt-order.R b/tests/testthat/test-cmt-order.R index bef0e7a29..202151830 100644 --- a/tests/testthat/test-cmt-order.R +++ b/tests/testthat/test-cmt-order.R @@ -673,7 +673,9 @@ rxTest({ s1 <- rxSolve(f1, ev2, returnType="data.frame") - expect_equal(s1, s2) + expect_equal(s1[,c("A1", "A2", "A3", "A4")], + s2[,c("A1", "A2", "A3", "A4")]) + }) }) From fec3ab7aa9b0d85253d24f48eb3b1d158ffe15c9 Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 10 Sep 2024 17:45:00 -0500 Subject: [PATCH 4/4] Write only when ddt is defined --- src/codegen.c | 20 +++++++++++++------- tests/testthat/test-convertId.R | 1 + 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/codegen.c b/src/codegen.c index 9ffb77515..e7da220d0 100644 --- a/src/codegen.c +++ b/src/codegen.c @@ -629,12 +629,17 @@ SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname, sIniTo(&sbOut, (int)((sbPm.sN)*5.3)); SEXP stateOrd = PROTECT(VECTOR_ELT(mvLast, RxMv_stateOrd)); pro++; - SEXP stateOrdNames = PROTECT(Rf_getAttrib(stateOrd, R_NamesSymbol)); pro++; - int *stateOrdInt = INTEGER(stateOrd); - sAppend(&sbOut, "// Define translation state order for %d states\n", Rf_length(stateOrd)); - for (int i = 0; i < Rf_length(stateOrd); i++){ - sAppend(&sbOut, "#define __DDT%d__ %d // %s\n", stateOrdInt[i]-1, i, - CHAR(STRING_ELT(stateOrdNames, i))); + int nOrd = Rf_length(stateOrd); + if (nOrd > 0) { + SEXP stateOrdNames = PROTECT(Rf_getAttrib(stateOrd, R_NamesSymbol)); pro++; + int *stateOrdInt = INTEGER(stateOrd); + sAppend(&sbOut, "// Define translation state order for %d states\n", Rf_length(stateOrd)); + for (int i = 0; i < nOrd; i++){ + sAppend(&sbOut, "#define __DDT%d__ %d // %s\n", stateOrdInt[i]-1, i, + CHAR(STRING_ELT(stateOrdNames, i))); + } + writeSb(&sbOut, fpIO); + sbOut.o = 0; } UNPROTECT(pro); // show_ode = 1 dydt @@ -680,8 +685,9 @@ SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname, } sAppend(&sbOut, "#define _CENTRAL_ %d\n", tb.statei); } + writeSb(&sbOut, fpIO); + sbOut.o = 0; } - writeSb(&sbOut, fpIO); gCode(1); // d/dt() gCode(2); // jac gCode(3); // ini() diff --git a/tests/testthat/test-convertId.R b/tests/testthat/test-convertId.R index bd1e07c7e..65179a4ab 100644 --- a/tests/testthat/test-convertId.R +++ b/tests/testthat/test-convertId.R @@ -31,6 +31,7 @@ rxTest({ ) f <- rxSolve(f, dMod) + expect_error(print(f), NA) })