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

[John Harrold]

I don't know if the files attached so here are their contents:

<mymod1.c>
/* file mymod1.c */
#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 */
</mymod1.c>

<mymod2.c>
/* file mymod2.c */
#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=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.c>

<testmod.r>
rm(list=ls())

library("deSolve")

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") }

file.copy('mymod1.c', 'mymod.c');
system("R CMD SHLIB mymod.c")
dyn.load("mymod.so")
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(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") }
file.copy('mymod2.c', 'mymod.c');
system("R CMD SHLIB mymod.c")
dyn.load("mymod.so")
dyn.unload(paste("mymod", .Platform$dynlib.ext, sep = ""))
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 Thu, Jan 21, 2016 at 4:38 PM, John Harrold <john.m.harrold at gmail.com>
wrote:

> Howdy,
>
> I accidentally posted this to the R-help list first.
>
> I'm constructing my models in C, and I've run into a problem. If I change
> the length of my parameters (say I increase it from 3 to 4), I get the
> following error:
>
> 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, 4; number in DLL, 3
>
> I suspect any that I'd have the same problem if I changed the number of
> states and that any modifications I make to the structure of the model
> won't take either. It seems that once I've loaded the .so file I'm stuck
> with it. I've tried to unload/load the .so file, but it doesn't seem to
> work.
>
> To illustrate this I'm using one of the compiled code examples. The
> attached script (testmod.r) will reproduce this on my Mac.
>
> I think my question is what is the appropriate way to unload a dynamic
> library such that the so or dll is forced to be reloaded and subsequently
> used?
>
> Thanks
> John
>



-- 
-------------------------------------
John M. Harrold
john.m.harrold _at_gmail
-------------------------------------

	[[alternative HTML version deleted]]