nlmixr2 3.0
By Matthew Fidler in rxode2 nlmixr2 babelmixr2
September 18, 2024
nlmixr2 3.0 is here
nlmixr2
3.0 has been released and most of the related packages in
the nlmixr2
ecosystem have been updated as well. Since there were a
few non-backward compatible changes and breaking changes, the version
has been bumped from 2 to 3. Most code will run the same, but because
of the breaking change, we changed the major version.
The big changes are:
Non abi binary linkages between every package. This means you will not get messages about needing to re-compile the packages during CRAN updates (and it means much more stable packages and the ability to make more frequent updates). To me, this is as important an achievement to maintainability as removing python’s
sympy
fromnlmixr
.New syntax for lower triangular matrices. I will likely discuss why in a later blog post, but for now I will show the differences:
# The old syntax (which is still supported)
lotri({
a+b ~ c(1,
0.1, 1)
c ~ 1
})
## a b c
## a 1.0 0.1 0
## b 0.1 1.0 0
## c 0.0 0.0 1
# the new syntax, lower triangular matrices can be suplied row by row;
# if restarting it would restart the block matrix
lotri({
a ~ 1
b ~ c(0.1, 1)
c ~ 1
})
## a b c
## a 1.0 0.1 0
## b 0.1 1.0 0
## c 0.0 0.0 1
rxode2
/nlmixr2
will re-order matrices as needed to make them in a banded-matrix format preferred by the estimation tool. This is done when parsing the matrices but can be seen here:
m <- lotri({
a + b + c + d + e + f + g + h + i + j + k + l + m + n + o +
p ~ c(0.4, 0, 0.3, 0, 0, 0, -0.1, 0, 0, 0.2, 0, 0, 0,
0, 0.5, 0, 0, 0, 0, 0, 1.3, 0, 0, 0, 0, 0, -0.6, 0.8,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.9, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.2, 0, 0.3,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.1, 0.2, 0, 0, 0.2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4, 0, 0, 0, 0, 0, -1.1,
0.9, 0, 0, 0, 0, 0, 0, 0, 4.7, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.5, 0, 0.2, 0, 0, 0, 1.9)
})
# note that m is in a non-banded matrix format
m
## a b c d e f g h i j k l m n o p
## a 0.4 0.0 0 -0.1 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.2 0.0 0.0
## b 0.0 0.3 0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## c 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## d -0.1 0.0 0 0.2 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.2 0.0 0.0
## e 0.0 0.0 0 0.0 0.5 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## f 0.0 0.0 0 0.0 0.0 1.3 -0.6 0 0.0 0.0 0.0 0.0 0.0 0.0 -1.1 0.0
## g 0.0 0.0 0 0.0 0.0 -0.6 0.8 0 0.0 0.0 0.0 0.0 0.0 0.0 0.9 0.0
## h 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## i 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## j 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.0 0.9 0.0 -0.2 0.0 0.0 0.0 0.5
## k 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.9 0.0 0.0 0.0 0.0 0.0
## l 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.0 -0.2 0.0 0.3 0.0 0.0 0.0 0.2
## m 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 2.1 0.0 0.0 0.0
## n 0.2 0.0 0 0.2 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.4 0.0 0.0
## o 0.0 0.0 0 0.0 0.0 -1.1 0.9 0 0.0 0.0 0.0 0.0 0.0 0.0 4.7 0.0
## p 0.0 0.0 0 0.0 0.0 0.0 0.0 0 0.0 0.5 0.0 0.2 0.0 0.0 0.0 1.9
# Internally the lotri function rcm (standing for the reverse Cuthill
# McKee (RCM) algorithm algorithm) creates the banded matrix
m <- lotri::rcm(m)
# If you specified this as an ini block, you would then see the
# following code instead:
lotri::lotriAsExpression(m, useIni=TRUE)
## ini({
## p ~ 1.9
## l ~ c(0.2, 0.3)
## j ~ c(0.5, -0.2, 0.9)
## o ~ 4.7
## g ~ c(0.9, 0.8)
## f ~ c(-1.1, -0.6, 1.3)
## n ~ 0.4
## d ~ c(0.2, 0.2)
## a ~ c(0.2, -0.1, 0.4)
## m ~ 2.1
## k ~ 0.9
## i ~ 0.2
## e ~ 0.5
## b ~ 0.3
## h ~ 0
## c ~ 0
## })
# which is both much easier to read and also helps the nlmixr2
# estimation routines.
# If you prefer the old syntax on output you can use
# options(lotri.plusNames=TRUE)
withr::with_options(list(lotri.plusNames=TRUE),{
lotri::lotriAsExpression(m, useIni=TRUE)
})
## ini({
## p + l + j ~ c(1.9, 0.2, 0.3, 0.5, -0.2, 0.9)
## o + g + f ~ c(4.7, 0.9, 0.8, -1.1, -0.6, 1.3)
## n + d + a ~ c(0.4, 0.2, 0.2, 0.2, -0.1, 0.4)
## m ~ 2.1
## k ~ 0.9
## i ~ 0.2
## e ~ 0.5
## b ~ 0.3
## h ~ 0
## c ~ 0
## })
In
rxode2
simulations,iCov
has been refactored to run faster than it did before.In
rxode2
you can now assign a variable to a string (a <- "string"
)Now you can specify per-variable interpolation (
linear(WT)
, ornocb(WT)
etc) inside your model.There is a new type of handling of
mu
-referenced equations, calledmu
3. (We already havemu
2).saem
will now message if it is usingmu
2 ormu
3 referencing while modeling. Standard mu referencing is not called out.There is a new log-likelihood profiling functions exported for use with you
nlmixr2
fits.
I have a technical discussion on why/how things have changed, or jump to the overall changes. The next section discusses why our dependencies have changed.
Why has the number of dependencies changed?
In
2022,
CRAN
requested we split rxode2
, giving us more packages to
maintain. This was mostly because of compilation time, and the
stan
-based, large-compilation time rxode2ll
is still its own
separate package.
In the last release, they requested that we reduce the number of
packages for rxode2
because it took too much of their time to
maintain the split they requested.
This is because the rxode2
, rxode2random
, rxode2et
and
rxode2parse
have binary
inter-dependencies
between all of the packages that were split off.
This means that rxode2
/nlmixr2est
had to depend on exact versions
of the packages rxode2random
, rxode2parse
, and rxode2et
. If
this did not happen, you could crash R. A sentence in their manual
talks about how fragile this is:
NB: this mechanism is fragile, as changes to the interface provided by packA have to be recognised by packB. The consequences of not doing so have included serious corruption to the memory pool of the R session. Either packB has to depend on the exact version of packA or there needs to be a mechanism for packB to test at runtime the version of packA it is linked to matches that it was compiled against.
I can see how this was very painful to maintain on CRAN’s end (and it was difficult to maintain from my end as well).
I must say, that CRAN is a small group of volunteers who manage a
great number of very useful R packages. They make sure
that the packages more or less interact well with each other, which is
actually a non-trivial task. This makes CRAN
just work, when compared
to other repositories where this is not enforced.
So, in short, “CRAN gives, CRAN takes away, blessed be the name of CRAN” (to mis-quote Job).
Note that I have taken other steps to make this easier to maintain on
my side and CRAN’s side too. For those who are interested I will talk
about them in below in the ABI/ABI
section.
Changelog
nlmixr2
- Export likelihood profiling functions
nlmixr2plot
- Changes to work with new
rxode2
nlmixr2extra
New
profile()
method for likelihood profiling (Issue #1)Changes to work with new
rxode2
nlmixr2est
No binary linking to
rxode2
,lbfgsb3c
andn1q1
, which means that updating these will not makenlmixr2est
crash without recompiling.New
mu
3 referencing will take context from the model to see if the algebraic expression can be completed from defined model variables; These variable would have to be unique.
rxode2
Breaking Changes
The model properties was moved from
$params
to$props
so it does not conflict with the low levelrxode2
model$params
Error when specifying
wd
withoutmodName
With Linear and midpoint of a time between two points, how
rxode2
handles missing values has changed. When the missing value is lower than the requested time, it will look backward until it finds the first non-missing value (or if all are missing start looking forward). When the missing value is higher than the requested time, the algorithm will look forward until it finds the first non-missing value (or if all are missing, start looking backward).The order of ODEs is now only determined by the order of
cmt()
andd/dt()
. Compartment properties,tad()
and other compartment related variables no no longer affect compartment sorting. The optionrxode2.syntax.require.ode.first
no longer does anything.The handling of zeros “safely” has changed (see #775)
when
safeZero=TRUE
and the denominator of a division expression is zero, use the Machine’s small number/eps
(you can see this value with.Machine$double.eps
)when
saveLog=TRUE
and the x in thelog(x)
is less than or equal to zero, change this tolog(eps)
when
safePow=TRUE
and the expressionx^y
has a zero forx
and a negative number fory
replacex
witheps
.
Since the protection for divide by zero has changed, the results will also change. This is a more conservative protection mechanism than was applied previously.
Random numbers from
rxode2
are different when usingdop853
,lsoda
orindLin
methods. These now seed the random numbers in the same way asliblsoda
, so the random number provided will be the same with different solving methods.The arguments saved in the
rxSolve
for items likethetaMat
will be the reduced matrices used in solving, not the full matrices (this will likely not break very many items)
Possible breaking changes (though unlikely)
iCov
is no longer merged to the event dataset. This makes solving withiCov
slightly faster (#743)
New features
You can remove covariances for every omega by piping with
%>% ini(diag())
you can be a bit more granular by removing all covariances that have eithereta.ka
oreta.cl
by:%>% ini(diag(eta.ka, eta.cl))
or anything with correlations witheta.cl
with%>% ini(diag(eta.cl))
You can also remove individual covariances by
%>% ini(-cov(a, b))
or%>% ini(-cor(a,b))
.You can specify the type of interpolation applied for added dosing records (or other added records) for columns that are kept with the
keep=
option inrxSolve()
. This new option iskeepInterpolation
and can belocf
for last observation carried forward,nocb
which is the next observation carried backward, as well asNA
which puts aNA
in all imputed data rows. See #756.Note: when interpolation is linear/midpoint for factors/characters it changes to locf with a warning (#759)
Also note, that the default keep interpolation is
na
Now you can specify the interpolation method per covariate in the model:
linear(var1, var2)
says bothvar1
andvar2
would use linear interpolation when they are a time-varying covariate. You could also uselinear(var1)
locf()
declares variables using last observation carried forwardnocb()
declares variables using next observation carried backwardmidpoint()
declares variables using midpoint interpolation
linear()
,locf()
,locb()
,midpoint()
,params()
,cmt()
anddvid()
declarations are now ignored when loading arxode2
model withrxS()
Strings can be assigned to variables in
rxode2
.Strings can now be enclosed with a single quote as well as a double quote. This limitation was only in the rxode2 using string since the R-parser changes single quotes to double quotes. (This has no impact with
rxode2({})
and ui/function form).More robust string encoding for symengine (adapted from
utils::URLencode()
andutils::URLdecode()
)Empty arguments to
rxRename()
give a warning (#688)Promoting from covariates to parameters with model piping (via
ini()
) now allows setting bounds (#692)Added
assertCompartmentName()
,assertCompartmentExists()
,assertCompartmentNew()
,testCompartmentExists()
,assertVariableExists()
testVariableExists()
,assertVariableNew()
,assertVariableName()
, andassertParameterValue()
to verify that a value is a valid nlmixr2 compartment name, nlmixr2 compartment/variable exists in the model, variable name, or parameter value (#726; #733)Added
assertRxUnbounded()
,testRxUnbounded()
,warnRxBounded()
to allownlmixr2
warn about methods that ignore boundaries #760Added functions
tad0()
,tafd0()
,tlast0()
andtfirst0()
that will give0
instead ofNA
when the dose has not been administered yet. This is useful for use in ODEs sinceNA
s will break the solving (so can be used a bit more robustly with models like Weibull absorption).rxode2
is has no more binary link tolotri
, which means that changes in thelotri
package will not requirerxode2
to be recompiled (in most cases) and will not crash the system.rxode2
also has no more binary linkage toPreciseSums
The binary linkage for
dparser
is reduced to C structures only, making changes in dparser less likely to cause segmentation faults inrxode2
if it wasn’t recompiled.A new model property has been added to
$props$cmtProp
and$statePropDf
. Both are data-frames showing which compartment has properties (currentlyini
,f
,alag
,rate
anddur
) in therxode2
ui model. This comes from the lower level model variable$stateProp
which has this information encoded in integers for each state.A new generic method
rxUiDeparse
can be used to deparse meta information into more readable expressions; This currently by default supports lower triangular matrices by lotri, but can be extended to support other types of objects like ’nlmixr2’sfoceiControl()
for instance.
Bug fixes
Fix
ui$props$endpoint
when the ui endpoint is defined in terms of the ode instead of lhs. See #754Fix
ui$props
when the ui is a linear compartment model withoutka
defined.Model extraction
modelExtract()
will now extract model properties. Note that the model property ofalag(cmt)
andlag(cmt)
will give the same value. See #745When assigning reserved variables, the parser will error. See #744
Linear interpolation will now adjust the times as well as the values when
NA
values are observed.Fix when keeping data has
NA
values that it will not crash R; Also fixed some incorrectNA
interpolations. See #756When using
cmt()
sometimes the next statement would be corrupted in the normalized syntax (like for instancelocf
); This bug was fixed (#763)keep
will now error when trying to keep items that are in the rxode2 output data-frame and will be calculated (#764)
Big change
- At the request of CRAN, combine
rxode2parse
,rxode2random
, andrxode2et
into this package; The changes in each of the packages are now placed here:
Fix for merged packages
Fix a bug when simulating nested variables for random simulation (was in
rxode2random
#25)As requested by CRAN remove the
C
codeSET_TYPEOF
which is no longer part of theC
R
API
forrxode2parse
ABI/API
I am writing this technical aside because I couldn’t find a solution
for this problem on the internet or in the writing R extension manual
(my internet searching skills could probably be improved). Also note
that if you are using Rcpp
and exposing the C++
calls, this is
done for you (but not for C
calls).
What is an ABI?
An ABI
stands for application binary interface and is determined by
the compiler and the platform at the time of compile. Each function
has an address in memory. This is why when some down-stream
dependency like rxode2parse
are updated the functions that are
specified change memory addresses. When rxode2
tries to call a
function from the parser, R memory’s is corrupted and R crashes. Of
course you can fix this by recompiling the rxode2
version but
usually waiting for the binary that links was a better
solution.
Addresses of other items other than C functions can also change (like C/C++ structures), so direct access of these items is also something that can cause memory corruption.
This ABI interface is what happens when you register a C
callable in
any package.
The solution – creating an API
An API
is the Application Programming Interface. It has been
popularized by web API
s and can be created for your R package with
packages like plumber
, however it can be any application programming
interface.
Step One – determine and export the C structures and C functions to be used
To create the API
I needed to do two things:
Figure out what functions I would be using (and therefore allow others to use) from
rxode2
.Figure out what structures that I would need to access (like the ODE states, the calculated values of the equations etc). This will also allow other packages to use these functions in their code as well.
The first thing that will need to be done is to create thin accessing
functions for the underlying C
structures. Since we are using R
which is really based on C
, they will need to be C compatible. The
code
below
shows how to access the calculated variables inside rxode2
:
double *getIndLhs(rx_solving_options_ind* ind) {
return ind->lhs;
}
The next step is to create a function that grabs the memory addresses
of these thin C structure access functions as well as any other functions
that you will expose in the API
.
You wrap all of these functions and place them into a R list structure
using R’s function R_MakeExternalPtrFn()
. An example function is
below and is based on the actual rxode2
API pointer
function:
SEXP _rxode2_rxode2Ptr(void) {
int pro = 0; // Counter for the number of PROTECT calls
// Create an external pointer for _lotriLstToMat
SEXP rxode2rxRmvnSEXP = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&_rxode2_rxRmvnSEXP, R_NilValue, R_NilValue)); pro++;
// more code here
// Unprotect all protected objects
UNPROTECT(pro);
// Return the list of external pointers
return ret;
}
After creating the C
function, you register it as you a standard
function R calls from C
. In this case we have:
void R_init_rxode2(DllInfo *info){
R_CallMethodDef callMethods[] = {
{"_rxode2_rxode2Ptr", (DL_FUNC) &_rxode2_rxode2Ptr, 0},
// more registered functions
{NULL, NULL, 0}
};
static const R_CMethodDef cMethods[] = {
{NULL, NULL, 0, NULL}
};
R_registerRoutines(info, cMethods, callMethods, NULL, NULL);
R_useDynamicSymbols(info, FALSE);
}
and of course you need to register the functions to be seen in R, which needs the following line in the NAMESPACE:
useDynLib(rxode2, .registration=TRUE)
Which for rxode2
comes from the roxygen2 line:
#' @useDynLib rxode2, .registration=TRUE
Once all of that is done, you can create an exported R function that lists the function pointers:
#' Get the rxode2 function pointers
#'
#' This function is used to get the function pointers for rxode2. This is
#' used to allow rxode2 to have binary linkage to nlmixr2est.
#'
#' @return a list of function pointers
#' @export
#' @author Matthew L. Fidler
#' @examples
#' .rxode2ptrs()
.rxode2ptrs <- function() {
.Call(`_rxode2_rxode2Ptr`, PACKAGE = "rxode2")
}
Which (of course) works for rxode2
:
rxode2::.rxode2ptrs()[1:3]
## $rxode2rxRmvnSEXP
## <pointer: 0x735bc70c68f0>
##
## $rxode2rxParProgress
## <pointer: 0x735bc72c49c0>
##
## $rxode2getRxSolve_
## <pointer: 0x735bc72c5ba0>
Step 2 create the user interface header
For me, I want the user interface functions to be the same as what I
declared in the API
, used in the second package (like nlmixr2est
).
To do this, I create a header called
rxode2ptr.h
in the inst/include
directory.
In that header, I declare the type definition of each of the functions that will be accessed in the other package as well as an external function pointer to the exported function:
typedef double *(*getIndLhs_t)(rx_solving_options_ind* ind);
extern getIndLhs_t getIndLhs;
Then create a static inline
function (that is it is only defined
when called) that takes the pointers from the rxode2::.rxode2ptrs()
and assigns them the function pointer getIndLhs
:
static inline SEXP iniRxodePtrs0(SEXP p) {
if (_rxode2_rxRmvnSEXP_ == NULL) { // only assign if not already assigned
// more code
getIndLhs = (getIndLhs_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 23));
// more code
}
return R_NilValue;
}
Note that this function is keyed to the number of API elements exposed
and assumes a specific order. Therefore, you should add elements
sequentially and make sure API
needs less elements than exported by
rxode2
(which is coded in as an exception in nlmixr2est
)
Finally we create code that defines the function pointers in a single file and creates an R expression that will register all the C pointer functions.
#define iniRxode2ptr \
_rxode2_rxRmvnSEXP_t _rxode2_rxRmvnSEXP_ = NULL; \
// more code
getIndLhs_t getIndLhs = NULL; \
SEXP iniRxodePtrs(SEXP ptr) { \
return iniRxodePtrs0(ptr); \
} \
This is all that needs to be done for rxode2
to expose the API
Step 3 using the API in nlmixr2est
For most C
(or C++
) files if the header has an accommodation for
C++, simply using:
#include <rxode2ptr.h>
Will register and expose the functions as needed.
In one file, you will also need to declare the function pointers and function pointer with the line:
#include <rxode2ptr.h>
// if in C++ wrap with extern "C" {}
// if you want to change the name of the function you can use #define
// iniRxodePtrs _your_fn_name
iniRxode2ptr // do not include a ; since it will be flagged as a significant warning by CRAN
When registering the Calls inside of nlmixr2est
you will need to add the ``
SEXP iniRxodePtrs(SEXP ptr);
/// some code
void R_init_nlmixr2est(DllInfo *info){
R_CallMethodDef callMethods[] = {
{"iniRxodePtrs", (DL_FUNC) &iniRxodePtrs, 1},
// more registered functions
{NULL, NULL, 0}
};
static const R_CMethodDef cMethods[] = {
{NULL, NULL, 0, NULL}
};
R_registerRoutines(info, cMethods, callMethods, NULL, NULL);
R_useDynamicSymbols(info, FALSE);
}
Last you need to register the API when loading nlmixr2est
:
# This will be saved when compiled
# This way you will know if something has changed in the API
rxode2.api <- names(rxode2::.rxode2ptrs())
.iniRxode2Ptr <- function() {
.ptr <- rxode2::.rxode2ptrs()
.nptr <- names(.ptr)
# Cheks for API changes
if (length(rxode2.api) > length(.nptr)) {
stop("nlmixr2est requires a newer version of rxode2 api, cannot run nlmixr2est\ntry `install.packages(\"rxode2\")` to get a newer version of rxode2", call.=FALSE)
} else {
.nptr <- .nptr[seq_along(rxode2.api)]
if (!identical(rxode2.api, .nptr)) {
.bad <- TRUE
stop("nlmixr2est needs a different version of rxode2 api, cannot run nlmixr2est\ntry `install.packages(\"rxode2\")` to get a newer version of rxode2, or update both packages", call.=FALSE)
}
}
.Call(`iniRxodePtrs`, .ptr,
PACKAGE = "nlmixr2est")
}
.onLoad <- function(libname, pkgname) {
# other code
.iniRxode2Ptr()
}
This same procedure can be used to expose all the API
functions to
R
.
- Posted on:
- September 18, 2024
- Length:
- 19 minute read, 3841 words
- Categories:
- rxode2 nlmixr2 babelmixr2
- See Also: