[R-sig-dyn-mod] running C code with multiple - plain text

[Andras Farkas]

Dear All,

I have two questions about an error message I received when runnin the following deSolve routine in C compiled code. 

#C Code 
/* file modelFME3.c */ 

#include <R.h> 

static double parms[9]; 
static double forc[5]; 
#define K parms[0] 
#define V parms[1] 
#define K12 parms[2] 
#define K21 parms[3] 
#define K13 parms[4] 
#define K31 parms[5] 
#define VP parms[6] 
#define DOSE parms[7] 
#define TAU parms[8] 
#define importr2 forc[0] 
#define importr3 forc[1] 
#define importr4 forc[2] 
#define importr5 forc[3] 
#define importr6 forc[4] 

/* initializer  */ 

void initmod(void (* odeparms)(int *, double *)) 

{ 
int N=9; 

odeparms(&N, parms); 

} 

void forcc(void (* odeforcs)(int *, double *)) 
{ 
int N=5; 
odeforcs(&N, forc); 
} 


/* Derivatives and 6 output variable */ 

void derivs (int *neq, double *t, double *y, double *ydot, 

double *yout, int *ip) 

{ 

double In; 
if ( fmod(*t,TAU) < 0.25 ){ 
In =  DOSE/0.25; 
} 
else { 
In  =  0; 
} 
if (ip[0] <1) error("nout should be at least 1"); 

ydot[0] = -(K12+K+K13*importr2+K13*importr3)*y[0]+(K21*y[1])+ 
(K31*y[2]*importr2)+(K31*y[3]*importr3)+(K31*y[4]*importr4) 
+(K31*y[5]*importr5)+(K31*y[6]*importr6); 

ydot[1] = K12*y[0]-K21*y[1]; 
ydot[2] = In*importr2 + K13*y[0]*importr2 - K31*y[2]*importr2; 
ydot[3] = K13*y[0]*importr3-K31*y[3]*importr3; 
ydot[4] = K13*y[0]*importr4-K31*y[4]*importr4; 
ydot[5] = K13*y[0]*importr5-K31*y[5]*importr5; 
ydot[6] = In*importr6 + K13*y[0]*importr6-K31*y[6]*importr6; 

yout[0] = y[0]/V; 
yout[1] = y[2]/VP; 
yout[2] = y[3]/VP; 
yout[3] = y[4]/VP; 
yout[4] = y[5]/VP; 
yout[5] = y[6]/VP; 

} 


/* END file modelFME3.c */ 


#R code
DOSE <-45 
TAU <-6 
pars <-c(K=0.01,V=12,K12=5,K21=5,K13=5, 
K31=5,VP=2.5,DOSE=DOSE,TAU=TAU) 

t    <- seq(0, 31,by=0.01) 

rtimes<-c(0.00 , 6.00, 12.00, 18.00, 20.25, 20.25, 
24.00, 24.15, 24.50, 25.50, 27.58, 27.60, 29.50 ,29.67) 
r2<-c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
r3<-c(0 ,1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
r4<-c(0 ,0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) 
r5<-c(0 ,0, 0,1,1,1, 0, 0, 0, 0, 0, 0, 0, 0) 
r6<-c(0 ,0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1) 


times  <- t 
signal <- as.data.frame(list(times = times, 
importr2 = rep(0, length(times)))) 
border1 <-ifelse(rtimes[sum(r2)]==max(rtimes), 
rtimes[sum(r2)],rtimes[sum(r2)+1]) 
signal$importr2[signal$times <=border1] <- 1 
ftime  <- seq(0, 31, 0.01) 
sigimp <- approxfun(signal$times, signal$importr2, rule = 2) 
Sigimp <- approx(signal$times, signal$importr2, 
xout=ftime ,rule = 2)$y 
forcingsr2 <- cbind(ftime, Sigimp) 

signal <- as.data.frame(list(times = times, 
importr3 = rep(0, length(times)))) 
border2 <-ifelse(rtimes[sum(r2,r3)]==max(rtimes), 
rtimes[sum(r2,r3)],rtimes[sum(r2,r3)+1]) 
signal$importr3[signal$times >border1&signal$times <=border2] <- 1 
ftime  <- seq(0, 31, 0.01) 
sigimp <- approxfun(signal$times, signal$importr3, rule = 2) 
Sigimp <- approx(signal$times, signal$importr3, xout=ftime ,rule = 2)$y 
forcingsr3 <- cbind(ftime, Sigimp) 

signal <- as.data.frame(list(times = times, 
importr4 = rep(0, length(times)))) 
border3 <-ifelse(rtimes[sum(r2,r3,r4)]==max(rtimes), 
rtimes[sum(r2,r3,r4)], 
rtimes[sum(r2,r3,r4)+1]) 
signal$importr4[signal$times >border2&signal$times <=border3] <- 1 
ftime  <- seq(0, 31, 0.01) 
sigimp <- approxfun(signal$times, signal$importr4, rule = 2) 
Sigimp <- approx(signal$times, signal$importr4, 
xout=ftime ,rule = 2)$y 
forcingsr4 <- cbind(ftime, Sigimp) 

signal <- as.data.frame(list(times = times, 
importr5 = rep(0, length(times)))) 
border4 <-ifelse(rtimes[sum(r2,r3,r4,r5)]==max(rtimes), 
rtimes[sum(r2,r3,r4,r5)], 
rtimes[sum(r2,r3,r4,r5)+1]) 
signal$importr5[signal$times >border3&signal$times <=border4] <- 1 
ftime  <- seq(0, 31, 0.01) 
sigimp <- approxfun(signal$times, signal$importr5, rule = 2) 
Sigimp <- approx(signal$times, signal$importr5, 
xout=ftime ,rule = 2)$y 
forcingsr5 <- cbind(ftime, Sigimp) 

signal <- as.data.frame(list(times = times, 
importr6 = rep(0, length(times)))) 
border5 <-ifelse(rtimes[sum(r2,r3,r4,r5,r6)]==max(rtimes), 
rtimes[sum(r2,r3,r4,r5,r6)], 
rtimes[sum(r2,r3,r4,r5,r6)+1]) 
signal$importr6[signal$times >border4&signal$times <=border5] <- 1 
ftime  <- seq(0, 31, 0.01) 
sigimp <- approxfun(signal$times, signal$importr6, rule = 2) 
Sigimp <- approx(signal$times, signal$importr6, 
xout=ftime ,rule = 2)$y 
forcingsr6 <- cbind(ftime, Sigimp) 

forcings<-list(forcingsr2,forcingsr3,forcingsr4, 
forcingsr5,forcingsr6) 


setwd("C:...\\fme") 
system("R CMD SHLIB modelFME3.c") 
dyn.load(paste("modelFME3", .Platform$dynlib.ext, sep = "")) 
#dyn.unload(paste("modelFME3", .Platform$dynlib.ext, sep = "")) 

model <- function(pars, times=seq(0,31, by = 0.01)) { 

state <- c(a = 0,b=0,c=0,d=0,e=0,f=0,g=0) 
return(ode(y = state, 
times = times, 
func = "derivs", 
parms = pars, 
dllname = "modelFME3", 
initforc = "forcc", 
forcings=forcings, 
method="lsoda", 
rtol = 1e-8, 
atol =1e-8, 
hmax=0.01, 
initfunc = "initmod", 
nout=6, 
outnames="cp.a","p2.c","p3.d","p4.e","p5.f","p6.g")) 
} 

data <-model(pars)

once we run the code an error message will be returned as follows:
Error in lsodar(y, times, func, parms, rtol, atol, jacfunc, jactype, 
rootfunc,  : ,jactype' must be one of 'fullint', 'fullusr', 'bandusr' or 'bandint....

1. if familiar with the error, could you please explain why the error
 message is prompted using the code above?

2. what specific suggestions would you have to correct the code above 
to allow the estimation of the variables of interest? 

Appreciate your input, 
Andras Farkas