[R-sig-dyn-mod] running ode in C with multiple forcings

[Andras Farkas]

Dear All,
wonder if you could provide input on the following:
#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 codeDOSE <-45TAU <-6pars <-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  <- tsignal <- 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] <- 1ftime  <- seq(0, 31, 0.01)sigimp <- approxfun(signal$times, signal$importr2, rule = 2)Sigimp <- approx(signal$times, signal$importr2, xout=ftime ,rule = 2)$yforcingsr2 <- 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] <- 1ftime  <- seq(0, 31, 0.01)sigimp <- approxfun(signal$times, signal$importr3, rule = 2)Sigimp <- approx(signal$times, signal$importr3, xout=ftime ,rule = 2)$yforcingsr3 <- 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] <- 1ftime  <- seq(0, 31, 0.01)sigimp <- approxfun(signal$times, signal$importr4, rule = 2)Sigimp <- approx(signal$times, signal$importr4, xout=ftime ,rule = 2)$yforcingsr4 <- 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] <- 1ftime  <- seq(0, 31, 0.01)sigimp <- approxfun(signal$times, signal$importr5, rule = 2)Sigimp <- approx(signal$times, signal$importr5, xout=ftime ,rule = 2)$yforcingsr5 <- 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] <- 1ftime  <- seq(0, 31, 0.01)sigimp <- approxfun(signal$times, signal$importr6, rule = 2)Sigimp <- approx(signal$times, signal$importr6, xout=ftime ,rule = 2)$yforcingsr6 <- 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)

this will give me an error message:
 Error in lsodar(y, times, func, parms, rtol, atol, jacfunc, jactype, rootfunc,  :   'jactype' must be one of 'fullint', 'fullusr', 'bandusr' or 'bandint' 
may not be the best approach but my models usually run without Jacobian, may or may not be relevant here... 
Any thoughts why this error message appears and thoughts onhow to get the model to run?
appreciate your insights,
thanks,
Andras

	[[alternative HTML version deleted]]