Skip to content

Commit

Permalink
Merge pull request #785 from nlmixr2/785-new-order-affects-c-code
Browse files Browse the repository at this point in the history 
New ordering of states can interfere with C code generation
  • Loading branch information
@mattfidler
mattfidler authored Sep 11, 2024
2 parents 6105d95 + fec3ab7 commit 7dbe392
Show file tree
Hide file tree
Showing 8 changed files with 179 additions and 28 deletions.
5 changes: 3 additions & 2 deletions inst/include/rxode2parse_control.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
24 changes: 20 additions & 4 deletions src/codegen.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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();
Expand All @@ -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 {
Expand Down Expand Up @@ -622,11 +623,25 @@ 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++;
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
// show_ode = 2 Jacobian
// show_ode = 3 Ini statement
Expand Down Expand Up @@ -671,6 +686,7 @@ SEXP _rxode2_codegen(SEXP c_file, SEXP prefix, SEXP libname,
sAppend(&sbOut, "#define _CENTRAL_ %d\n", tb.statei);
}
writeSb(&sbOut, fpIO);
sbOut.o = 0;
}
gCode(1); // d/dt()
gCode(2); // jac
Expand Down
9 changes: 6 additions & 3 deletions src/genModelVars.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
16 changes: 8 additions & 8 deletions src/parseCmtProperties.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
10 changes: 5 additions & 5 deletions src/parseDdt.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
}
Expand Down
18 changes: 12 additions & 6 deletions src/rxData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
124 changes: 124 additions & 0 deletions tests/testthat/test-cmt-order.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")

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[,c("A1", "A2", "A3", "A4")],
s2[,c("A1", "A2", "A3", "A4")])

})

})
1 change: 1 addition & 0 deletions tests/testthat/test-convertId.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ rxTest({
)

f <- rxSolve(f, dMod)

expect_error(print(f), NA)

})
Expand Down

0 comments on commit 7dbe392

Please sign in to comment.