-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathTRPallminopt.R
98 lines (80 loc) · 5.49 KB
/
TRPallminopt.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
taskid <- as.numeric(commandArgs(trailingOnly=TRUE))[1]
ntasks <- as.numeric(commandArgs(trailingOnly=TRUE))[2]
tname <- commandArgs(trailingOnly=TRUE)[3]
targetpt <- commandArgs(trailingOnly=TRUE)[4]
DST <- commandArgs(trailingOnly=TRUE)[5]
resume <- commandArgs(trailingOnly=TRUE)[6]
noneq <- commandArgs(trailingOnly=TRUE)[7]
if (noneq) paramvary <- list(commandArgs(trailingOnly=TRUE)[8], as.numeric(commandArgs(trailingOnly=TRUE)[9]))
location<-"../scratch/"
tag <- "20160419p"
currenttag <- paste0(tname,"_",tag)
if (targetpt=="DS") rDSTall <- TRUE else rDSTall <- FALSE #commandArgs(trailingOnly=TRUE)[5]
drtag <- ifelse(rDSTall == TRUE, paste0("rDSTall.",currenttag), currenttag)
tasktag <- paste0(currenttag,".",taskid)
source("TPPmat.R")
dssetup <- setup.model(DRera=FALSE, treatSL=FALSE, treatnovel=FALSE)
drsetup <- setup.model(DRera=TRUE, treatSL=TRUE, treatnovel=FALSE)
novelsetup <- setup.model(DRera=TRUE, treatSL=TRUE, treatnovel=TRUE)
values <- set.values()
genericvalues <- mergedvalues <- append(append(values[[1]], values[[2]]), append(values[[3]], values[[4]]))
tallynames <- colnames(equilib()$log)[-(1:(length(dssetup$statenames)+1))]
elementnames <- c("all",set.novelvalues()$elementnames)
alldrout <- numeric(0)
i <- 1; while(file.exists(paste0(location,"DRcalibration_",drtag,".",i,".csv")))
{alldrout <- rbind(alldrout, read.csv(paste0(location,"DRcalibration_",drtag,".",i,".csv"), header = TRUE)); i <- i+1} #saved results from dr sampling runs at time 0
tolerance <- 1.5
drout <- alldrout[alldrout[,"rrinc"]/alldrout[,"inc"] > 1/tolerance*alldrout[,"targetdr"] & alldrout[,"rrinc"]/alldrout[,"inc"] < tolerance*alldrout[,"targetdr"], ]
header <- c("inew", "ids","idr","targetprev","targetcoprev", "targetdr", "targetpt","DST", "rDSTall", names(unlist(genericvalues)))
header <- append(header, paste0( rep(tallynames, times=2*11), rep(0:10,each=length(tallynames)), rep(c("allmin", "allopt"), each=11*length(tallynames)) ) )
if (noneq)
{ if(!file.exists(paste0(location,"Allminopt.",paramvary[[1]],paramvary[[2]],"_", targetpt,DST,"_",tasktag,".csv"))) { write(header, file=paste0(location,"Allminopt.",paramvary[[1]],paramvary[[2]],"_", targetpt,DST,"_",tasktag,".csv"), sep=",", ncol=length(header)) }
}else
if(!file.exists(paste0(location,"Allminopt","_", targetpt,DST,"_",tasktag,".csv"))) { write(header, file=paste0(location,"Allminopt","_", targetpt,DST,"_",tasktag,".csv"), sep=",", ncol=length(header)) }
ilimits <- ceiling(seq(0,nrow(drout), length=ntasks+1))
istart <- ilimits[taskid]+1
if (resume) istart <- max(read.csv(paste0(location,"Allminopt","_", targetpt,DST,"_",tasktag,".csv"))$inew)+1
if (istart > ilimits[taskid+1]) stop("Already finished this taskid")
for (inew in istart:ilimits[taskid+1])
{
iter <- unlist(c(inew,unlist(drout[inew,c("ids", "idr", "targetprev","targetcoprev","targetdr")]), targetpt, DST, rDSTall)) #will include these labels as part of returned output
valuevect <- drout[inew, 5+(1:length(unlist(genericvalues)))] ; names(valuevect) <- names(unlist(genericvalues))
v <- 0; for (pname in names(genericvalues))
{ genericvalues[[pname]] <- unlist(valuevect[v+(1:length(genericvalues[[pname]]))]); #names(genericvalues[[pname]]) <- names()
v <- v + length(genericvalues[[pname]]) # move forward to start of next par vector in sampled values
}
state <- drout[inew,drsetup$statenames]
newstate <- numeric(length(novelsetup$statenames)); names(newstate) <- novelsetup$statenames
newstate[drsetup$statenames] <- unlist(state)
iresult <- numeric(0) #will become full row of results for this dr row (inew), for all elements and levels
for (level in c("minimal","optimal"))
{
print(paste0("Evaluating TRP optimal for all ",level," for Simulation #", inew," of ",nrow(drout)," for ",targetpt,DST,tasktag))
m <- o <- NA
if (level=="minimal") m <- elementnames[2:7] else m <- NA
if (level=="optimal") o <- elementnames[2:7] else o <- NA
valueset <- sampleTRP(mergedvalues = genericvalues, targetpt = targetpt, DST = DST,
minimals=m, optimals=o , HIV="nonHIV")
s_cr <- valueset$cres[1]; r_cr <- valueset$cres[2]
novelstate <- newstate
for (name in novelsetup$statenames) if (length(grep("^S", name))==0 & length(grep("^C", name))==0 )
{
if (length(grep("+c+",name))==1 )
{
if (length(grep("+.Rr+", name))==1 )
{ novelstate[name] <- r_cr * newstate[str_replace(string = name, pattern = "c", replacement = "")]
} else novelstate[name] <- s_cr * max(newstate[str_replace(string = name, pattern = "c", replacement = "")], newstate[str_replace(string = name, pattern = "c", replacement = "0")], na.rm = TRUE)
} else
{
if (length(grep("+.Rr+", name))==1 )
{ novelstate[name] <- (1-r_cr) * newstate[name]
} else novelstate[name] <- (1-s_cr) * newstate[name]
}
}
parset <- create.pars(setup = novelsetup, values = valueset, T, T, T)
outset <- ode(y=unlist(novelstate), times=0:10, func=dxdt, parms=parset$fullpars, do.tally=TRUE, paramvary=paramvary, method="adams")
iresult <- append(iresult, as.vector(t(outset[,tallynames])))
}
if(noneq) write(unlist(c(iter, valuevect, iresult)), file=paste0(location,"Allminopt.",paramvary[[1]],paramvary[[2]],"_", targetpt,DST,"_",tasktag,".csv"), sep=",", append=TRUE, ncol=length(header)) else
write(unlist(c(iter, valuevect, iresult)), file=paste0(location,"Allminopt","_", targetpt,DST,"_",tasktag,".csv"), sep=",", append=TRUE, ncol=length(header))
}