Skip to content

[R-sig-dyn-mod] Model development with compiled code

2 messages · Gibbons, Frank, John Harrold

Original
Gibbons, Frank
# 
Hi John,


I found your thread on R-sig-dynamic-models<https://stat.ethz.ch/pipermail/r-sig-dynamic-models/2016q1/000436.html> about the "Confusion over the length of parms" message. Looks like you figured it out, but I'm new at this and can't for the life of me figure it out. Care to share?





Kind regards,



-Frank




PS:
$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"



________________________________

Confidentiality Notice: This message is private and may ...{{dropped:10}}
John Harrold
# 
Howdy Frank,

Well i turns out that I didn't reply because I was kind of ashamed of why I
was making the mistake. I made the correct change in parms[4] but I forget
to make the corresponding change in the initmod function:

    int N=4;

:).

This should work. If you save the text below as a script and run it, it
should work.

<testmod.r>


rm(list=ls())


mymod1 = "
#include <R.h>
static double parms[3];
#define k1 parms[0]
#define k2 parms[1]
#define k3 parms[2]

/* initializer  */
void initmod(void (* odeparms)(int *, double *))
{
    int N=3;
    odeparms(&N, parms);
}

/* Derivatives and 1 output variable */
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip)
{
    if (ip[0] <1) error(\"nout should be at least 1\");
    ydot[0] = -k1*y[0] + k2*y[1]*y[2];
    ydot[2] = k3 * y[1]*y[1];
    ydot[1] = -ydot[0]-ydot[2];

    yout[0] = y[0]+y[1]+y[2];
}

/* The Jacobian matrix */
void jac(int *neq, double *t, double *y, int *ml, int *mu,
           double *pd, int *nrowpd, double *yout, int *ip)
{
  pd[0]               = -k1;
  pd[1]               = k1;
  pd[2]               = 0.0;
  pd[(*nrowpd)]       = k2*y[2];
  pd[(*nrowpd) + 1]   = -k2*y[2] - 2*k3*y[1];
  pd[(*nrowpd) + 2]   = 2*k3*y[1];
  pd[(*nrowpd)*2]     = k2*y[1];
  pd[2*(*nrowpd) + 1] = -k2 * y[1];
  pd[2*(*nrowpd) + 2] = 0.0;
}
/* END file mymod.c */
"

mymod2 = "
#include <R.h>
static double parms[4];
#define k1 parms[0]
#define k2 parms[1]
#define k3 parms[2]
#define k4 parms[3]

/* initializer  */
void initmod(void (* odeparms)(int *, double *))
{
    int N=4;
    odeparms(&N, parms);
}

/* Derivatives and 1 output variable */
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip)
{
    if (ip[0] <1) error(\"nout should be at least 1\");
    ydot[0] = -k1*y[0] + k2*y[1]*y[2];
    ydot[2] = k3 * y[1]*y[1];
    ydot[1] = -ydot[0]-ydot[2];

    yout[0] = y[0]+y[1]+y[2];
}

/* The Jacobian matrix */
void jac(int *neq, double *t, double *y, int *ml, int *mu,
           double *pd, int *nrowpd, double *yout, int *ip)
{
  pd[0]               = -k1;
  pd[1]               = k1;
  pd[2]               = 0.0;
  pd[(*nrowpd)]       = k2*y[2];
  pd[(*nrowpd) + 1]   = -k2*y[2] - 2*k3*y[1];
  pd[(*nrowpd) + 2]   = 2*k3*y[1];
  pd[(*nrowpd)*2]     = k2*y[1];
  pd[2*(*nrowpd) + 1] = -k2 * y[1];
  pd[2*(*nrowpd) + 2] = 0.0;
}
/* END file mymod.c */
"


library("deSolve")

#
# Writing the first file
#

if(file.exists(paste("mymod", .Platform$dynlib.ext, sep = ""))){
  file.remove( paste("mymod", .Platform$dynlib.ext, sep = "")) }
if(file.exists("mymod.c")){
   file.remove("mymod.c") }
if(file.exists("mymod.o")){
   file.remove("mymod.o") }

fileConn<-file("mymod.c")
writeLines(mymod1, fileConn)
close(fileConn)


# Compiling and loading the file
system("R CMD SHLIB mymod.c")
dyn.load(paste("mymod", .Platform$dynlib.ext, sep = ""))

parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7)
Y     <- c(y1 = 1.0, y2 = 0.0, y3 = 0.0)
times <- c(0, 0.4*10^(0:11) )
out <- ode(Y, times, func = "derivs", parms = parms,
           jacfunc = "jac", dllname = "mymod",
           initfunc = "initmod", nout = 1, outnames = "Sum")


if(('mymod' %in% names(getLoadedDLLs()))){
  dyn.unload(getLoadedDLLs()$mymod[["path"]])}

if(file.exists(paste("mymod", .Platform$dynlib.ext, sep = ""))){
  file.remove( paste("mymod", .Platform$dynlib.ext, sep = "")) }
if(file.exists("mymod.o")){
   file.remove("mymod.o") }
if(file.exists("mymod.c")){
   file.remove("mymod.c") }

fileConn<-file("mymod.c")
writeLines(mymod2, fileConn)
close(fileConn)

system("R CMD SHLIB mymod.c")
dyn.load(paste("mymod", .Platform$dynlib.ext, sep = ""))

parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7, k4 = 1.0)
Y     <- c(y1 = 1.0, y2 = 0.0, y3 = 0.0)
times <- c(0, 0.4*10^(0:11) )
out <- ode(Y, times, func = "derivs", parms = parms,
           jacfunc = "jac", dllname = "mymod",
           initfunc = "initmod", nout = 1, outnames = "Sum")

</testmod.r>



On Wed, May 10, 2017 at 10:50 AM Gibbons, Frank <
Frank.Gibbons at astrazeneca.com> wrote: