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

[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:

> 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:
>
> > 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"
>
>
>
> ________________________________
>
> Confidentiality Notice: This message is private and may ...{{dropped:10}}
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>

	[[alternative HTML version deleted]]