radivot — Oct 31, 2012, 5:41 PM
library(myelo)
Loading required package: deSolve
tpts=c(0, 5, 20, 60, 120, 180, 240, 300) # the time points at which we have data, i.e. tpts=unique(epo$time)
(f=file.path(system.file(paste("libs",Sys.getenv("R_ARCH"),sep=""), package = "myelo"),
paste("myelo",.Platform$dynlib.ext,sep="")))
[1] "/usr/local/lib/R/site-library/myelo/libs/myelo.so"
dyn.load(f)
sig=c(summary(lm(epo.e~as.factor(time),data=epo))$coef[2,2],
summary(lm(epo.m~as.factor(time),data=epo))$coef[2,2],
summary(lm(epo.i~as.factor(time),data=epo))$coef[2,2])
(Sig=rep(1,24)%*%t(sig)) # make the sd's match the dims of the data
[,1] [,2] [,3]
[1,] 5804 1236 2152
[2,] 5804 1236 2152
[3,] 5804 1236 2152
[4,] 5804 1236 2152
[5,] 5804 1236 2152
[6,] 5804 1236 2152
[7,] 5804 1236 2152
[8,] 5804 1236 2152
[9,] 5804 1236 2152
[10,] 5804 1236 2152
[11,] 5804 1236 2152
[12,] 5804 1236 2152
[13,] 5804 1236 2152
[14,] 5804 1236 2152
[15,] 5804 1236 2152
[16,] 5804 1236 2152
[17,] 5804 1236 2152
[18,] 5804 1236 2152
[19,] 5804 1236 2152
[20,] 5804 1236 2152
[21,] 5804 1236 2152
[22,] 5804 1236 2152
[23,] 5804 1236 2152
[24,] 5804 1236 2152
require(bbmle)
Loading required package: bbmle
Loading required package: stats4
nLLHC<-function(kde,kdi,ke,kex,kon,kt,scale,d,ts,IW)
{ p0=10^c(kde,kdi,ke,kex,kon,kt,scale)
yout <- ode(y = c(Epo=2100,EpoR=516,EpoEpoR=0,EpoEpoRi=0,dEpoi=0,dEpoe=0), times = ts,
func = "derivscHC", parms = p0, dllname = "myelo",initfunc = "parmscHC",
nout = 3, outnames = c("y1", "y2","y3"))
E=cbind(rep(yout[,8],each=3),rep(yout[,9],each=3),rep(yout[,10],each=3))
res=(E-d[,2:4])/IW
sum(res^2)
}
thetahatHC=c(kde = -2,kdi=-3, ke= -1, kex=-2,kon=-4,kt=-2,scale= 2)
fit0 <- mle2(nLLHC,start=as.list(thetahatHC),
method="L-BFGS-B",
data=list(d=epo,ts=tpts,IW=Sig))
(s=summary(fit0))
Maximum likelihood estimation
Call:
mle2(minuslogl = nLLHC, start = as.list(thetahatHC), method = "L-BFGS-B",
data = list(d = epo, ts = tpts, IW = Sig))
Coefficients:
Estimate Std. Error z value Pr(z)
kde -1.84462 0.01191 -154.9 <2e-16 ***
kdi -2.67768 0.04268 -62.7 <2e-16 ***
ke -1.19591 0.00856 -139.7 <2e-16 ***
kex -2.60420 0.08262 -31.5 <2e-16 ***
kon -3.93587 0.00938 -419.5 <2e-16 ***
kt -1.62612 0.01160 -140.1 <2e-16 ***
scale 2.20740 0.00115 1923.2 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-2 log L: 1236
myplotHC=function(p) {
p0=10^p
yout <-ode(y = c(Epo=2100,EpoR=516,EpoEpoR=0,EpoEpoRi=0,dEpoi=0,dEpoe=0),
times = tpts, func = "derivscHC", parms = p0,
dllname = "myelo",initfunc = "parmscHC",
nout = 3, outnames = c("y1", "y2","y3"))
with(epo,matplot(time,epo[,-1]))
matlines(yout[,1],yout[,-(1:7)])
}
myplotHC(coef(s)[,1]) # fit still looks good
myprof=function(k) profile(fit0,which=k,alpha=0.05,maxsteps=20)
library(parallel)
cl <- makeCluster(4)
clusterCall(cl, function() getwd())
[[1]]
[1] "/home/radivot"
[[2]]
[1] "/home/radivot"
[[3]]
[1] "/home/radivot"
[[4]]
[1] "/home/radivot"
clusterCall(cl, function() Sys.info()[c("nodename","machine")])
[[1]]
nodename machine
"CML" "x86_64"
[[2]]
nodename machine
"CML" "x86_64"
[[3]]
nodename machine
"CML" "x86_64"
[[4]]
nodename machine
"CML" "x86_64"
clusterEvalQ(cl, require(myelo))
[[1]]
[1] TRUE
[[2]]
[1] TRUE
[[3]]
[1] TRUE
[[4]]
[1] TRUE
clusterEvalQ(cl, require(bbmle))
[[1]]
[1] TRUE
[[2]]
[1] TRUE
[[3]]
[1] TRUE
[[4]]
[1] TRUE
clusterExport(cl, ls())
clusterCall(cl, function() ls(.GlobalEnv))
[[1]]
[1] "Sig" "cl" "f" "fit0" "myplotHC"
[6] "myprof" "nLLHC" "s" "sig" "thetahatHC"
[11] "tpts"
[[2]]
[1] "Sig" "cl" "f" "fit0" "myplotHC"
[6] "myprof" "nLLHC" "s" "sig" "thetahatHC"
[11] "tpts"
[[3]]
[1] "Sig" "cl" "f" "fit0" "myplotHC"
[6] "myprof" "nLLHC" "s" "sig" "thetahatHC"
[11] "tpts"
[[4]]
[1] "Sig" "cl" "f" "fit0" "myplotHC"
[6] "myprof" "nLLHC" "s" "sig" "thetahatHC"
[11] "tpts"
clusterCall(cl, function() myprof)
[[1]]
function (k)
profile(fit0, which = k, alpha = 0.05, maxsteps = 20)
[[2]]
function (k)
profile(fit0, which = k, alpha = 0.05, maxsteps = 20)
[[3]]
function (k)
profile(fit0, which = k, alpha = 0.05, maxsteps = 20)
[[4]]
function (k)
profile(fit0, which = k, alpha = 0.05, maxsteps = 20)
clusterCall(cl, function() dyn.load(f))
[[1]]
DLL name: myelo
Filename: /usr/local/lib/R/site-library/myelo/libs/myelo.so
Dynamic lookup: TRUE
[[2]]
DLL name: myelo
Filename: /usr/local/lib/R/site-library/myelo/libs/myelo.so
Dynamic lookup: TRUE
[[3]]
DLL name: myelo
Filename: /usr/local/lib/R/site-library/myelo/libs/myelo.so
Dynamic lookup: TRUE
[[4]]
DLL name: myelo
Filename: /usr/local/lib/R/site-library/myelo/libs/myelo.so
Dynamic lookup: TRUE
system.time(PL<-clusterApply(cl, 1:7, myprof)) # takes ~95 secs on two, so not quite cut in half from 140
user system elapsed
0.004 0.000 54.230
par(mfrow=c(2,4)) # 4 jobs to one cpu and 3 to the other => expect this
junk<-lapply(PL,plot) # on a quad proc the time was 46 secs, i.e. 3 cpus with 2 to do, 1 with 1.
stopCluster(cl)