library(survival)

source('sensitivity/R/sensitivitySGD.R')


vax<-scan("vax004po.txt",
          list(sid=integer(0), Vx=integer(0), sex=integer(0), age=integer(0), white=integer(0), black=integer(0),
               asian=integer(0), hisp=integer(0), otherrace=integer(0), education=integer(0),
               US=integer(0), Canada=integer(0), Netherlands=integer(0), site=integer(0), region=integer(0),
               baselineriskscore=integer(0), calenddayssinceentry=integer(0),
               HIVinfection=integer(0), daysinfectiongroup=integer(0), daysinfectiondiagnosis=integer(0), treatmentstart=integer(0),
               daystreatmentstart=integer(0), composite=integer(0), dayscomposite=integer(0), dayssinceinfectiondiagnosis=integer(0),
               grouptime=integer(0), viralload=double(0), CD4list=integer(0)))

id<-unique(vax$sid)

#trtstart<-daystrtstart<-comp<-dayscomp<-z<-white<-risk<-s<-NULL
d <- y <- z <- s <- white <- NULL
count<-1


## assuming vax is sorted by id
## finds for each id
##       treatment or placebo
##       length of treatment
##
##       Vx
##       is white
##       baseline risk score
##       was infected
for(i in 1:length(vax$sid)){
  if (vax$sid[i]==id[count]) {
     d[count]<-vax$treatmentstart[i]
    y[count]<-vax$daystreatmentstart[i]/30
#    comp[count]<-vax$composite[i]
#    dayscomp[count]<-vax$dayscomposite[i]
    z[count]<-vax$Vx[i]
    white[count]<-vax$white[i]
#    risk[count]<-vax$baselineriskscore[i]
    s[count]<-vax$HIVinfection[i]
    count<-count+1
  }
}

z <- z[white==0]
s <- s[white==0]
d <- d[white==0]
y <- y[white==0]

## number patients
N<-length(z)
## number vaxinated
N1 <- sum(z)
## number placebo
N0 <- N-N1

## number infected that where not vacinated
n0 <- sum((1-z)*s)
## number infected that where vacinated
n1 <- sum(z*s)

p0 <- n0/N0
p1 <- n1/N1

#d<-trtstart
#y<-daystrtstart/30    # Converting to months

beta0 <- seq(from=-3,to=3,length.out=101)/12
beta1 <- seq(from=-3,to=3,length.out=101)/12
tau <- 24
t.point <- seq(1,24,by=1)

phi <- c(0.99, 0.95, 0.90, 0.8, 0.5, n0/N0)
psi <- sapply(phi, p0=p0, p1=p1,
              function(phi, p0, p1) log((1-p0-p1+phi*p1)*(phi*p1) /
                                        ((p0 - phi*p1)*(p1 - phi*p1))))
#gctorture(TRUE)

#result <- sensitivity2d.calc.mk3(z=z, s=s, d=d, y=y, beta0=beta0, beta1=beta1, psi=psi,
#                                 tau=tau, time.points=time.points)

listans <- list()

for(p in psi) {
  set.seed(4633097)
  listans[[as.character(p)]] <- sensitivitySGD(z, s, d, y, beta0=beta0,
                                               beta1=beta1, psi=p, tau=tau,
                                               time.points=t.point, ci=0.95,
                                               selection=1L,
                                               groupings=c(0L,1L),
                                               empty.principle.stratum=c(0L,1L),
                                               trigger=1L,
                                               ci.method='bootstrap',
                                               N.boot=1000, verbose=TRUE)
  cat("\n")
}

ans <- do.call(mapply, c(listans, list(FUN=function(...) {
  dotlist <- list(...)
  dlen <- length(dim(dotlist[[1]]))
  dims <- replace(dim(dotlist[[1]]), dlen, length(dotlist))
  dimnames <- replace(dimnames(dotlist[[1]]), dlen, list(names(dotlist)))

  ans <- cbind(...)
  if (length(dims) == 0) {
    dim(ans) <- NULL
    dimnames(ans) <- NULL
  } else {
    dim(ans) <- dims
    dimnames(ans) <- dimnames
  }

  ans
})))

str(ans)
save(ans, phi, beta0, beta1, tau, t.point, file='vaxgen_dataset.rda')