define state order in model

nlmixr2 · Sep 10, 2024 · 5dd2b5d · 5dd2b5d
1 parent cbfc04a
commit 5dd2b5d
Show file tree

Hide file tree

Showing 7 changed files with 52 additions and 30 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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")])
+
   })
 
 })