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

[Daniel Kaschek]

Hi Andras,

nice puzzle ;-)

The "error" in the code has nothing to do with the error message. In
the argument "outname" you just forgot putting c() around your list of
names. Omitting this, every name was interpreted as another argument in
the chain of possible arguments to lsoda().

Best regards,
Daniel

On Sun, 2017-03-12 at 18:51 +0000, Andras Farkas wrote:
> 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
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models