[Rcpp-devel] Using Rcpp to solve ODEs: LSODE (deSolve) says "Confusion over the length of parms"

[Gibbons, Frank]

Hi,

Just joined the list. I've searched online for help, and could find little mention of the problem I'm having, and zero solutions. I'm getting a strange error message when using Rcpp to combine a C-implementation of some simple linear ODEs with deSolve in package "popED". (It's a package for optimal design of experiments. Because it relies on repeated solution of the ODEs at each iterative step, a pure-R implementation is too slow for me on non-trivial problems.)

I've attached the R code (it implements a pure-R/deSolve solution, then a deSolve/Rcpp implementation). The problem is one of pharmacokinetics: what's the optimal time (or dose) to use in order to best estimate parameters that describe drug absorption, distribution and elimination from the body, using the parameters KA, V, and CL respectively. The design can be constrained in a number of ways. The script does two things for each model implementation: first it evaluates an initial design, second it attempts to optimize that design by parameter search, which requires repeated evaluation of the (ODE-based) model. Strange thing is: it all works fine for the initial model evaluation, so the code works (and Rcpp is 3x faster, expect bigger gains in the optimization step). But then it falls over in the model optimization step:

> tic(); output <- poped_optim(poped.db.compiled, opt_xt =TRUE, parallel=F, method = c("LS")); toc()
Error in lsoda(y, times, func, parms, ...) :
  Confusion over the length of parms
In addition: Warning message:
In lsoda(y, times, func, parms, ...) :
  Number of parameters passed to solver, 3; number in DLL, 5
Called from: lsoda(y, times, func, parms, ...)

How can it run to evaluate the model, then fail when it comes to optimization? Since the error is coming from LSODA (the solver that is wrapped by deSolve), I'm guessing the problem lies in interface between deSolve and the code I've got. I'm no stranger to language extension (years ago, I used SWIG to extend some Python libraries with C code), but I'm struggling to understand where to look for the source of this problem. In particular, it's telling me the DLL has 5 parameters, but I just don't see where that could be coming from: I explicitly pass in three parameters (the KA, V, and CL from above). Does Rcpp inject extras that are somehow being exposed inappropriately?

One last thing: you may look at this model and ask why I'm not using the obvious closed-form solution for the concentration profile - that's the fastest way to the answer. The answer is that this isn't really the problem I'm interested in, and those don't have an easy closed-form solution.

Would really appreciate any pointers you can give. Many thanks and kind regards,

-Frank







P.S.
> R.Version()
$platform
[1] "x86_64-w64-mingw32"

$arch
[1] "x86_64"

$os
[1] "mingw32"

$system
[1] "x86_64, mingw32"

$status
[1] ""

$major
[1] "3"

$minor
[1] "3.2"

$year
[1] "2016"

$month
[1] "10"

$day
[1] "31"

$`svn rev`
[1] "71607"

$language
[1] "R"

$version.string
[1] "R version 3.3.2 (2016-10-31)"

$nickname
[1] "Sincere Pumpkin Patch"



>




Frank Gibbons, PhD
Principal Scientist
_____________________________________________________________________________________________
DMPK Modeling & Simulation | IMED Oncology | AstraZeneca
35 Gatehouse Lane, Waltham MA  02451, USA
+1-781-839-4032   frank.gibbons at astrazeneca.com<mailto:frank.gibbons at astrazeneca.com>


________________________________

Confidentiality Notice: This message is private and may contain confidential and proprietary information. If you have received this message in error, please notify us and remove it from your system and note that you must not copy, distribute or take any action in reliance on it. Any unauthorized use or disclosure of the contents of this message is not permitted and may be unlawful.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20170120/e43c17cb/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: myExampleODE.compiled.R
Type: application/octet-stream
Size: 4343 bytes
Desc: myExampleODE.compiled.R
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20170120/e43c17cb/attachment-0001.obj>