[R-sig-dyn-mod] Fw: forcing parameters for modMCMC in C

[Andras Farkas]

 Andras Farkas, Pharm.D CEO, Optimum Dosing Strategies LLC 63, Reeve Ave, Bloomingdale, NJ, 07403 www.optimum-dosing-strategies.org

     On Friday, April 14, 2017 4:07 PM, Andras Farkas <motyocska at yahoo.com> wrote:
 

  Dear All,
I have this R code:

require(deSolve) 
require(FME) 
pars <- list(k1 = 0.08,k2=0.2,v1=15,v2=50) 
intimes <- c(0,0.5,12,12.5,24) 
input  <- c(800,0,1500,0,0) 
forc <- approxfun(intimes, input, method="constant", rule=2) 
derivs <- function(t, state, pars) { 
kintimes <- c(0,12.01) 
kinput  <- list(pars$k1,pars$k2) 
names(kinput) <-c("k1","k2") 
kforc <- approxfun(kintimes, kinput, method="constant", rule=2) 
vintimes <- c(0,12.01) 
vinput  <- list(pars$v1,pars$v2) 
names(vinput) <-c("v1","v2") 
vforc <- approxfun(vintimes, vinput, method="constant", rule=2) 
k <- kforc(t) 
v <- vforc(t) 
inp <- forc(t) 
dy1 <- - k * state[1] + inp 
dc <- state[1]/v 
return(list(dy1,dc)) 
} 

model <- function(pars, times=seq(0, 25, by = 0.1)) { 
state <- 0 
return(ode(y = state, times = times, func = derivs, parms = pars,method="lsoda",rtol = 1e-8,atol =1e-8)) 
} 

observed <- matrix (nc=2,byrow=2,data=c(0,0,0.5,26.14,11.5,10.84,13.9,17.91,20,5.288)) 
colnames(observed) <- c("time", "2") 

tout <- unique(sort(c(intimes, observed[,"time"]))) 

obj <- function(x, parset = names(x)) { 
pars[parset] <- x 
out <- model(pars, tout) 
return(modCost(obs = observed, model = out)) 
} 

prior <- function(p) 
return( sum(((p-c(0.15,0.15,30,30))/c(0.08,0.08,20,20))^2 )) 

final <- modMCMC(p = c(k1=0.08,k2=0.16,v1=20,v2=45), f = obj, prior=prior,lower=c(0.0001,0.00001,0.0001,0.0001),jump = NULL, niter = 100,updatecov=10, wvar0 = 1,var0=NULL,burninlength = 10) 


summary(final) 

this will allow the estimation of 2 sets of parameters for k and v respectively over time...

would you know of an example that implements these forcings of parameters in C, or thoughts on how to write the code? Got this far, which will work with ode, ie I can get the parameter change implemented but when I try modMCMC, the identification of bestpar fails...

#R code for c run

require(deSolve) 
require(FME) 
pars <- c(k1 = 0.08,k2=0.2,v1=15,v2=50) 
intimes <- c(0,0.5,12,12.5,24) 
input  <- c(800,0,1500,0,0) 
forc <- approxfun(intimes, input, method="constant", rule=2) 
input<-data.frame(t,forc(t)) 
t<-seq(0, 25, by = 0.1) 
kintimes <- c(0,12.01) 
kinput  <- c(pars[1:2]) 
names(kinput) <-c("k1","k2") 
kforc <- approxfun(kintimes, kinput, method="constant", rule=2) 
kin <-data.frame(t,kforc(t)) 

vintimes <- c(0,12.01) 
vinput  <- c(pars[3:4]) 
names(vinput) <-c("v1","v2") 
vforc <- approxfun(vintimes, vinput, method="constant", rule=2) 
vin <-data.frame(t,vforc(t)) 

forcings<-list(input,vin,kin) 

setwd("D:\\ANDRAS\\motyocska\\Backup\\cafb0e06-2f6e-46f8-8b35-79ed0e7c9c81\\20130907_062400_motyocska\\G\\drug info\\antibiotics\\Publications\\ICAAC\\2018\\vancoITCSf") 
system("R CMD SHLIB modelSIM_forcedtest.c") 
dyn.load(paste("modelSIM_forcedtest", .Platform$dynlib.ext, sep = "")) 
#dyn.unload(paste("modelSIM_parsforced", .Platform$dynlib.ext, sep = "")) 


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

state <- c(a=0) 
return(ode(y = state, 
times = times, 
func = "derivs", 
parms = pars, 
dllname = "modelSIM_forcedtest", 
initforc = "forcc", 
forcings=forcings, 
method="lsoda", 
rtol = 1e-8, 
atol =1e-8, 
hmax=0.01, 
initfunc = "initmod", 
nout=1, 
outnames=c("cp"))) 
} 


modelout <-data.frame(model(pars)) 
plot(modelout$time,modelout$cp,type="l")

observed <- matrix (nc=2,byrow=2,data=c(0,0,0.5,26.14,11.5,10.84,13.9,17.91,20,5.288)) 
colnames(observed) <- c("time", "cp") 

tout <- unique(sort(c(intimes, observed[,"time"]))) 

obj <- function(x, parset = names(x)) { 
pars[parset] <- x 
out <- model(pars, tout) 
return(modCost(obs = observed, model = out)) 
} 

prior <- function(p) 
return( sum(((p-c(0.15,0.15,30,30))/c(0.08,0.08,20,20))^2 )) 

final <- modMCMC(p = c(k1=0.08,k2=0.16,v1=20,v2=45), f = obj, prior=prior,lower=c(0.0001,0.00001,0.0001,0.0001),jump = NULL, niter = 100,updatecov=10, wvar0 = 1,var0=NULL,burninlength = 10) 


summary(final)

#C code

/* file modelSIM.c */ 

#include <R.h>

static double parms[4];

static double forc[3];
#define v11 parms[0] 
#define v12 parms[1]
#define k11 parms[2] 
#define k12 parms[3] 
#define import1 forc[0]
#define import2 forc[1] 
#define import3 forc[2]

/* initializer  */

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

{
int N=4;

odeparms(&N, parms);

}

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

/* Derivatives and 1 output variable */ 

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

double *yout, int *ip) 

{    double k,v; 


if (ip[0] <1) error("nout should be at least 1"); 

v=import2; 
k=import3; 
ydot[0]= import1 -k* y[0]; 
yout[0]=y[0]/v;

} 


/* END file modelSIM.c */


goal is to get back final$bestpar with compiled code similar to the results by the non-compiled setup...

appreciate your help...


Andras


Andras Farkas, 


   
	[[alternative HTML version deleted]]