diff --git a/R/linCmt.R b/R/linCmt.R index 60eb8c68..d73cf2dc 100644 --- a/R/linCmt.R +++ b/R/linCmt.R @@ -88,8 +88,7 @@ linCmt <- function(data, ..., lagdepot=0, lagcentral=0, addlKeepsCov=FALSE, addlDropSs=TRUE, ssAtDoseTime=TRUE, scale=1, gradient=FALSE, - keep=NULL, - covsInterpolation = c("locf", "linear", "nocb", "midpoint")) { + keep=NULL) { checkmate::assertNumeric(fdepot, len=1, lower=0, finite=TRUE, any.missing=FALSE) checkmate::assertNumeric(fcentral, len=1, lower=0, finite=TRUE, any.missing=FALSE) checkmate::assertNumeric(ratedepot, len=1, lower=0, finite=TRUE, any.missing=FALSE) @@ -103,8 +102,8 @@ linCmt <- function(data, ..., checkmate::assertLogical(addlDropSs, len=1, any.missing=FALSE) checkmate::assertLogical(ssAtDoseTime, len=1, any.missing=FALSE) checkmate::assertLogical(gradient, len=1, any.missing=FALSE) - covsInterpolation <- match.arg(covsInterpolation) # First take care of the possible inputs to lag time + .linNamesData <- names(data) .w <- grepl(.rxLinInfoKaReg, .linNamesData, ignore.case=TRUE) if (length(.w) > 0L) { @@ -189,6 +188,7 @@ linCmt <- function(data, ..., "dur(central)=durcentral\n", "alag(depot)=lagdepot\n", "alag(central)=lagcentral\n"), linear=TRUE) + } data$rxRowNum <- seq_along(data[,1]) .data <- as.data.frame(etTransParse(data, .mv, dropUnits=TRUE, allTimeVar = TRUE, @@ -210,33 +210,14 @@ linCmt <- function(data, ..., return(rep(0.0, .l)) } } - .i <- seq_along(.dati[,.t]) - .v <- .dati[,.t] - .w <- which(!is.na(.v)) - .yleft <- .v[.w[1]] - .yright <- .v[.w[length(.w)]] - if (covsInterpolation == "locf") { - fun <- approxfun(.i[.w], .v[.w], method="constant", - yleft = .yleft, yright = .yright) - } else if (covsInterpolation == "nocb") { - fun <- approxfun(.i[.w], .v[.w], method="constant", - yleft = .yleft, yright = .yright, f = 1) - } else if (covsInterpolation == "midpoint") { - fun <- approxfun(.i[.w], .v[.w], method="constant", - yleft = .yleft, yright = .yright, f = 0.5) - } else if (covsInterpolation == "linear") { - fun <- approxfun(.i[.w], .v[.w], yleft = .yleft, yright = .yright) - } else { - stop("unknown covsInterpolation", call.=FALSE) - } - fun(.i) + .dati[,.t] }), names(.trans$str))) - .ret <- .Call(`_rxode2parse_compC`, - list(.dati[, c("TIME", "EVID", "AMT", "II")], .extra, .trans$trans), .mv) - .ret$ID <- i - .ret + .pars <- + list(.dati[, c("TIME", "EVID", "AMT", "II")], + .extra, + .trans$trans) }) - dplyr::as_tibble(do.call(`rbind`, .ret)) + .ret } #' Calculate the lambdas and coefficients of the two compartment model #' diff --git a/inst/include/rxode2parseGetTime.h b/inst/include/rxode2parseGetTime.h index 9671f968..2c3af735 100644 --- a/inst/include/rxode2parseGetTime.h +++ b/inst/include/rxode2parseGetTime.h @@ -29,21 +29,14 @@ extern t_calc_mtime calc_mtime; } \ } + static inline double getLag(rx_solving_options_ind *ind, int id, int cmt, double time) { rx_solving_options *op = &op_global; returnBadTime(time); if (ind->wh0 == EVID0_SS0 || ind->wh0 == EVID0_SS20) { return time; } - double ret = NA_REAL; - if (ind->cmt < 0) { - } else if (ind->cmt < op->neq) { - ret = LAG(id, cmt, time); - } else if (ind->cmt == op->neq) { - ret = ind->linCmtConstants[linCmt_tlag]; - } else if (ind->cmt == op->neq+1) { - ret = ind->linCmtConstants[linCmt_tlag2]; - } + double ret = LAG(id, cmt, time); if (ISNA(ret)) { op->badSolve=1; op->naTime = 1; @@ -54,15 +47,7 @@ static inline double getLag(rx_solving_options_ind *ind, int id, int cmt, double static inline double getRate(rx_solving_options_ind *ind, int id, int cmt, double dose, double t){ rx_solving_options *op = &op_global; returnBadTime(t); - double ret = NA_REAL; - if (ind->cmt < 0) { - } else if (ind->cmt < op->neq) { - ret = RATE(id, cmt, dose, t); - } else if (ind->cmt == op->neq) { - ret = ind->linCmtConstants[linCmt_rate1]; - } else if (ind->cmt == op->neq+1) { - ret = ind->linCmtConstants[linCmt_rate2]; - } + double ret = RATE(id, cmt, dose, t); if (ISNA(ret)){ op->badSolve=1; op->naTime = 1; @@ -74,15 +59,7 @@ static inline double getDur(rx_solving_options_ind *ind, int id, int cmt, double rx_solving_options *op = &op_global; returnBadTime(t); if (ISNA(t)) return t; - double ret = NA_REAL; - if (ind->cmt < 0) { - } else if (ind->cmt < op->neq) { - ret = DUR(id, cmt, dose, t); - } else if (ind->cmt == op->neq) { - ret = ind->linCmtConstants[linCmt_dur1]; - } else if (ind->cmt == op->neq+1) { - ret = ind->linCmtConstants[linCmt_dur2]; - } + double ret = DUR(id, cmt, dose, t); if (ISNA(ret)){ op->badSolve=1; op->naTime = 1; diff --git a/inst/include/rxode2parseHandleEvid.h b/inst/include/rxode2parseHandleEvid.h index 660c662d..0d5d7064 100644 --- a/inst/include/rxode2parseHandleEvid.h +++ b/inst/include/rxode2parseHandleEvid.h @@ -32,30 +32,6 @@ extern "C" { #define rxErrModeledFss2n3 4194304 #define rxErrRate02 8388608 -#define linCmt_tlag 0 -#define linCmt_F 1 -#define linCmt_rate1 2 -#define linCmt_dur1 3 -#define linCmt_tlag2 4 -#define linCmt_F2 5 -#define linCmt_rate2 6 -#define linCmt_dur2 7 - - static inline void updateLinCmtVars(rx_solving_options_ind *ind, - double *d_tlag, double *d_F, double *d_rate1, double *d_dur1, - // Oral parameters - double *d_tlag2, double *d_F2, double *d_rate2, double *d_dur2) { - double *l=ind->linCmtConstants; - l[linCmt_tlag] = *d_tlag; - l[linCmt_tlag] = *d_tlag2; - l[linCmt_F] = *d_F; - l[linCmt_F2] = *d_F2; - l[linCmt_rate1] = *d_rate1; - l[linCmt_dur1] = *d_dur1; - l[linCmt_rate2] = *d_rate2; - l[linCmt_dur2] = *d_dur2; - } - #if defined(__cplusplus) #define FLOOR(x) std::floor(x) extern "C" { @@ -344,15 +320,7 @@ extern t_F AMT; static inline double getAmt(rx_solving_options_ind *ind, int id, int cmt, double dose, double t, double *y) { - double ret = NA_REAL; - rx_solving_options *op = &op_global; - if (ind->cmt < op->neq) { - ret = AMT(id, cmt, dose, t, y); - } else if (ind->cmt == op->neq) { - ret = ind->linCmtConstants[linCmt_F]*dose; - } else if (ind->cmt == op->neq+1) { - ret = ind->linCmtConstants[linCmt_F2]*dose; - } + double ret = AMT(id, cmt, dose, t, y); if (ISNA(ret)){ rx_solving_options *op = &op_global; op->badSolve=1; diff --git a/inst/include/rxode2parseStruct.h b/inst/include/rxode2parseStruct.h index b6a6ca99..ca98af0a 100644 --- a/inst/include/rxode2parseStruct.h +++ b/inst/include/rxode2parseStruct.h @@ -290,7 +290,6 @@ typedef struct { //double *extraDoseIi; // ii doses unsupported bool lastIsSs2; double *timeThread; - double *linCmtConstants; } rx_solving_options_ind; typedef struct { diff --git a/man/linCmt.Rd b/man/linCmt.Rd index d5da13f2..ac5cdfe1 100644 --- a/man/linCmt.Rd +++ b/man/linCmt.Rd @@ -20,8 +20,7 @@ linCmt( ssAtDoseTime = TRUE, scale = 1, gradient = FALSE, - keep = NULL, - covsInterpolation = c("locf", "linear", "nocb", "midpoint") + keep = NULL ) } \arguments{ @@ -103,12 +102,13 @@ having to compile any \code{rxode2} models \examples{ d <- nlmixr2data::nmtest -names(d) <- sub("lagt", "lagcentral", - sub("bioav", "fdepot", - sub("rat2", "ratecentral", - sub("dur2", "durcentral", names(d))))) +names(d) <- names(d) \%>\% + sub("lagt", "lagcentral",.) \%>\% + sub("bioav", "fdepot", .) \%>\% + sub("rat2", "ratecentral", .) \%>\% + sub("dur2", "durcentral", .) -ret <- linCmt(d, cl=1.1, v=20, ka=1.5) +linCmt(d, cl=1.1, v=20, ka=1.5) } \author{ diff --git a/src/comp.c b/src/comp.c index 4a491dd0..b3c88a39 100644 --- a/src/comp.c +++ b/src/comp.c @@ -253,8 +253,6 @@ double linCmtCompA(rx_solve *rx, unsigned int id, double _t, int linCmt, double d_ka, double d_tlag2, double d_F2, double d_rate2, double d_dur2) { lin_context_c_t lin; rx_solving_options_ind *ind = &(rx->subjects[id]); - updateLinCmtVars(ind, &d_tlag, &d_F, &d_rate1, &d_dur1, - &d_tlag2, &d_F2, &d_rate2, &d_dur2); rx_solving_options *op = rx->op; double t = _t - ind->curShift; unsigned int ncmt=0; @@ -341,7 +339,7 @@ double linCmtCompA(rx_solve *rx, unsigned int id, double _t, int linCmt, } if (!op->badSolve){ ind->idx = i; - if (getEvid(ind, ind->ix[i]) == 3) { + if (getEvid(ind, ind->ix[i]) == 3){ ind->curShift -= rx->maxShift; for (int j = op->neq + op->extraCmt; j--;) { ind->InfusionRate[j] = 0; @@ -366,23 +364,17 @@ double linCmtCompA(rx_solve *rx, unsigned int id, double _t, int linCmt, return(yp[oral0]/lin.v); } -extern double getTimeParse(int idx, rx_solving_options_ind *ind) { - return getTime__(idx, ind, 0); -} SEXP _rxode2parse_compC(SEXP in, SEXP mv) { - rx_solve *rx=&rx_global; - rx->op = &op_global; + rx_solve *rx=(&rx_global); rx_solving_options *op = rx->op; rx_solving_options_ind *oldInd = rx->subjects; - t_getTime curGetTime = _rxode2parse_getTime; - _rxode2parse_getTime = getTimeParse; iniSolvingRx(rx); iniSolvingOptions(op); int pro = 0; SEXP dat = PROTECT(VECTOR_ELT(in, 0)); pro++; SEXP par = PROTECT(VECTOR_ELT(in, 1)); pro++; - int trans = INTEGER(VECTOR_ELT(in, 2))[0]; + int trans = INTEGER(VECTOR_ELT(in, 2)); double rate[2]; rate[0] = rate[1] = 0.0; int cnt = 0; @@ -392,10 +384,8 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { indR.dadt_counter = &cnt; indR.jac_counter = &cnt; indR.InfusionRate = rate; - int BadDose[2]; - BadDose[0] = BadDose[1] = 0; - indR.BadDose = BadDose; - indR.nBadDose = 0; + //indR.BadDose + // indR.nBadDose indR.HMAX = 0.0; // Determined by diff indR.curDose = NA_REAL; indR.dosenum = 0; @@ -433,8 +423,8 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { // ..$ TIME: num [1:135] 0 0 1 2 3 4 5 6 7 8 ... - indR.all_times = REAL(VECTOR_ELT(dat, 0)); - indR.n_all_times = Rf_length(VECTOR_ELT(dat, 0)); + indR.all_times = REAL(VECTOR_ELT(dat, 1)); + indR.n_all_times = Rf_length(VECTOR_ELT(dat, 1)); // ..$ EVID: int [1:135] 101 0 0 0 0 0 0 0 0 0 ... indR.evid = INTEGER(VECTOR_ELT(dat, 1)); // $EVID; rx->nall = Rf_length(VECTOR_ELT(dat, 1)); @@ -482,12 +472,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { op->kind = 0; op->extraCmt= INTEGER(VECTOR_ELT(mv, RxMv_extraCmt))[0]; op->nDisplayProgress = 0; - double ssAtol[2]; - ssAtol[0] = ssAtol[1] = 1.0e-8; - double ssRtol[2]; - ssRtol[0] = ssRtol[1] = 1.0e-6; - op->ssAtol=ssAtol; - op->ssRtol=ssRtol; SEXP flagsS = VECTOR_ELT(mv, RxMv_flags); int *flags = INTEGER(flagsS); rx->linKa = flags[RxMvFlag_ka]; @@ -561,7 +545,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { int *idose = malloc(2*indR.n_all_times*sizeof(int)); if (idose == NULL) { - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } /* int *ix; */ @@ -580,7 +563,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { int *tmpI = (int*)malloc(EVID_EXTRA_SIZE* sizeof(int)); if (tmpI == NULL) { free(idose); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.ignoredDoses = tmpI; @@ -588,7 +570,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { if (tmpI == NULL) { free(idose); free(indR.ignoredDoses); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.ignoredDosesN = tmpI; @@ -608,7 +589,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { free(idose); free(indR.ignoredDoses); free(indR.ignoredDosesN); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.pendingDoses = tmpI; @@ -618,7 +598,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { free(indR.ignoredDoses); free(indR.ignoredDosesN); free(indR.pendingDoses); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.extraDoseTimeIdx = tmpI; @@ -629,7 +608,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { free(indR.ignoredDosesN); free(indR.pendingDoses); free(indR.extraDoseTimeIdx); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.extraDoseEvid = tmpI; @@ -640,11 +618,9 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { free(indR.ignoredDosesN); free(indR.pendingDoses); free(indR.extraDoseTimeIdx); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.extraDoseTime = tmpD; - tmpD = (double*)malloc(EVID_EXTRA_SIZE* sizeof(double)); if (tmpD == NULL) { free(idose); free(indR.ignoredDoses); @@ -652,7 +628,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { free(indR.pendingDoses); free(indR.extraDoseTimeIdx); free(indR.extraDoseTime); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.extraDoseDose = tmpD; @@ -660,7 +635,7 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { indR.extraDoseNewXout = NA_REAL; indR.idxExtra = 0; indR.extraSorted = 0; - tmpD = (double*)calloc(indR.n_all_times*(rx->linNcmt+rx->linKa+1), sizeof(double)); + tmpD = (double*)calloc(indR.n_all_times*(rx->linNcmt+rx->linKa), sizeof(double)); if (tmpD == NULL) { free(idose); free(indR.ignoredDoses); @@ -669,16 +644,9 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { free(indR.extraDoseTimeIdx); free(indR.extraDoseTime); free(indR.extraDoseDose); - _rxode2parse_getTime = curGetTime; Rf_errorcall(R_NilValue, _("ran out of memory")); } indR.solve = tmpD; - indR.timeThread = tmpD + indR.n_all_times*(rx->linNcmt+rx->linKa); - double linCmtConstants[8]; - linCmtConstants[0]=linCmtConstants[1]=linCmtConstants[2]= - linCmtConstants[3]=linCmtConstants[4]=linCmtConstants[5]= - linCmtConstants[6]=linCmtConstants[7]=0.0; - indR.linCmtConstants = linCmtConstants; indR.ixds = indR.idx=0; rx_solving_options_ind* ind = &indR; @@ -732,87 +700,6 @@ SEXP _rxode2parse_compC(SEXP in, SEXP mv) { free(indR.extraDoseDose); free(indR.solve); rx->subjects = oldInd; - - SEXP dfNames = PROTECT(Rf_allocVector(STRSXP, 19)); pro++; - SEXP dfVals = PROTECT(Rf_allocVector(VECSXP, 19)); pro++; - SEXP rnVals = PROTECT(Rf_allocVector(INTSXP, 2)); pro++; - int *rnI = INTEGER(rnVals); - - rnI[0] = NA_INTEGER; - rnI[1] = -indR.n_all_times; - - //"TIME" - SET_STRING_ELT(dfNames, 0, Rf_mkChar("TIME")); - SET_VECTOR_ELT(dfVals, 0, VECTOR_ELT(dat, 0)); - - // "EVID" - SET_STRING_ELT(dfNames, 1, Rf_mkChar("EVID")); - SET_VECTOR_ELT(dfVals, 1, VECTOR_ELT(dat, 1)); - - // "AMT" - SET_STRING_ELT(dfNames, 2, Rf_mkChar("AMT")); - SET_VECTOR_ELT(dfVals, 2, VECTOR_ELT(dat, 2)); - - // "II" - SET_STRING_ELT(dfNames, 3, Rf_mkChar("II")); - SET_VECTOR_ELT(dfVals, 3, VECTOR_ELT(dat, 3)); - - // "Cc" - SET_STRING_ELT(dfNames, 3, Rf_mkChar("Cc")); - SET_VECTOR_ELT(dfVals, 3, CcSxp); - - SET_STRING_ELT(dfNames, 4, Rf_mkChar("p1")); - SET_VECTOR_ELT(dfVals, 4, VECTOR_ELT(par, 0)); - - SET_STRING_ELT(dfNames, 5, Rf_mkChar("v1")); - SET_VECTOR_ELT(dfVals, 5, VECTOR_ELT(par, 1)); - - SET_STRING_ELT(dfNames, 6, Rf_mkChar("p2")); - SET_VECTOR_ELT(dfVals, 6, VECTOR_ELT(par, 2)); - - SET_STRING_ELT(dfNames, 7, Rf_mkChar("p3")); - SET_VECTOR_ELT(dfVals, 7, VECTOR_ELT(par, 3)); - - SET_STRING_ELT(dfNames, 8, Rf_mkChar("p4")); - SET_VECTOR_ELT(dfVals, 8, VECTOR_ELT(par, 4)); - - SET_STRING_ELT(dfNames, 9, Rf_mkChar("p5")); - SET_VECTOR_ELT(dfVals, 9, VECTOR_ELT(par, 5)); - - SET_STRING_ELT(dfNames, 10, Rf_mkChar("lagdepot")); - SET_VECTOR_ELT(dfVals, 10, VECTOR_ELT(par, 6)); - - SET_STRING_ELT(dfNames, 11, Rf_mkChar("fdepot")); - SET_VECTOR_ELT(dfVals, 11, VECTOR_ELT(par, 7)); - - SET_STRING_ELT(dfNames, 12, Rf_mkChar("ratedepot")); - SET_VECTOR_ELT(dfVals, 12, VECTOR_ELT(par, 8)); - - SET_STRING_ELT(dfNames, 13, Rf_mkChar("durdepot")); - SET_VECTOR_ELT(dfVals, 13, VECTOR_ELT(par, 9)); - - SET_STRING_ELT(dfNames, 14, Rf_mkChar("ka")); - SET_VECTOR_ELT(dfVals, 14, VECTOR_ELT(par, 10)); - - SET_STRING_ELT(dfNames, 15, Rf_mkChar("lagcentral")); - SET_VECTOR_ELT(dfVals, 15, VECTOR_ELT(par, 11)); - - SET_STRING_ELT(dfNames, 16, Rf_mkChar("fcentral")); - SET_VECTOR_ELT(dfVals, 16, VECTOR_ELT(par, 12)); - - SET_STRING_ELT(dfNames, 17, Rf_mkChar("ratecentral")); - SET_VECTOR_ELT(dfVals, 17, VECTOR_ELT(par, 13)); - - SET_STRING_ELT(dfNames, 18, Rf_mkChar("durcentral")); - SET_VECTOR_ELT(dfVals, 18, VECTOR_ELT(par, 14)); - - // R_NameSymbol - SEXP cls = PROTECT(Rf_allocVector(STRSXP, 1)); pro++; - SET_STRING_ELT(cls, 0, Rf_mkChar("data.frame")); - Rf_setAttrib(dfVals, R_NamesSymbol, dfNames); - Rf_setAttrib(dfVals, R_RowNamesSymbol, rnVals); - Rf_setAttrib(dfVals, R_ClassSymbol, cls); - _rxode2parse_getTime = curGetTime; UNPROTECT(pro); - return dfVals; + return R_NilValue; } diff --git a/src/init.c b/src/init.c index 1f557e62..34dece67 100644 --- a/src/init.c +++ b/src/init.c @@ -99,8 +99,6 @@ SEXP _rxode2parse_chin(SEXP x, SEXP table); SEXP _rxode2parse_getForder(void); int _rxode2parse_useForder(void); -SEXP _rxode2parse_compC(SEXP in, SEXP mv); - int get_sexp_uniqueL( SEXP s ); SEXP _rxode2parse_funPtrs(void) { @@ -139,7 +137,6 @@ SEXP _rxode2parse_solComp3(SEXP sK10, SEXP sK12, SEXP sK21, SEXP sK13, SEXP sK31 SEXP _rxode2parse_solComp2(SEXP sK10, SEXP sK12, SEXP sK21); void R_init_rxode2parse(DllInfo *info){ R_CallMethodDef callMethods[] = { - {"_rxode2parse_compC", (DL_FUNC) &_rxode2parse_compC, 2}, {"_rxode2parse_solComp3",(DL_FUNC) &_rxode2parse_solComp3, 5}, {"_rxode2parse_solComp2",(DL_FUNC) &_rxode2parse_solComp2, 3}, {"_rxode2parse_linCmtR1", (DL_FUNC) &_rxode2parse_linCmtR1, 2},