# Program 16
# This first segment of this program searches for a D_s-optimal design for the logistic example
# in Chapter 4
f <- function(xvec)
{ 
 c(1,xvec,xvec[1]*xvec[2])
}
f2 <- function(xvec)
{
 c(1,xvec)
}

ratiodets <- function(variables)
{#Calculates the ratio of the determinants of the information matrices for the full and reduced models
 xvals <- cos(pi*variables[1:lim1])
 xvalsmat <- t(matrix(xvals,s,m))
 wtvals <- (variables[(lim1+1):((m+1)*s)])^2
 deltavec <- wtvals/sum(wtvals)
 fxmat <- apply(xvalsmat,2,f)
 f2xmat <- apply(xvalsmat,2,f2)
 etavec <- as.vector(t(betavec)%*%fxmat)
 expetavec <- exp(etavec)
 modelwtvec <- expetavec/((1+expetavec)^2)
 infomat <- fxmat%*%diag(deltavec*modelwtvec)%*%t(fxmat)
 infomat2 <- f2xmat%*%diag(deltavec*modelwtvec)%*%t(f2xmat)
 -det(infomat)/det(infomat2)
}


m <- 2
p <- 4
p1 <- 1
p2 <- p - p1
s <- 4
lim1 <- m*s
betavec <- c(1,1.1,1,-1)

nsims <- 100
mindet <- 10
#set.seed(123)
for (i in 1:nsims)
{
 initial <- c(runif(lim1),runif(s))
 out <- optim(initial,ratiodets,NULL,method="Nelder-Mead")
 valuenow <- out$val
 if(valuenow < mindet){mindet <- valuenow
 design <- out$par}
}
cat("Min value of ratio\n",mindet,"\n") 
output <- design
out1 <- cos(pi*output[1:lim1])
out2 <- (output[(lim1+1):((m+1)*s)])^2
wts <- out2/sum(out2)
out4 <- cbind(matrix(out1,s,m),wts)
out4 <- out4[order(out4[,1]),]
cat("Design\n")
out5 <- t(out4)
out5


# This segment of the program produces the standardised variance function for the D_s-optimal design
# You need to enter details of the alleged optimal design.  xvalsmat contains the coordinates of
# the support points (1st point, then 2nd point, then ...) and deltavec contains the design weights
# in the same order as the support points are entered in the columns of xvalsmat 

 infomat <- matrix(0,p,p)
 infomat2 <- matrix(0,p2,p2)
 xvalsmat <- matrix(c(-1,-1,-1,1,1,-1,1,1),m,s)
 deltavec <- c(0.255,0.235,0.255,0.255)
 for (i in 1:s)
 {
  xvec <- xvalsmat[,i]
  fvec <- f(xvec)
  f2vec <- f2(xvec)
  eta <- sum(fvec*betavec)
  expeta <- exp(eta)
  wt <- expeta/((1+expeta)^2)
  infomat <- infomat + deltavec[i]*wt*fvec%*%t(fvec)
  infomat2 <- infomat2 + deltavec[i]*wt*f2vec%*%t(f2vec)
 }
invinfo <- solve(infomat)
invinfo2 <- solve(infomat2)

stdvar1 <- function(xvec)
{ 
 fvec <- f(xvec)
 f2vec <- f2(xvec)
 eta <- sum(fvec*betavec)
 expeta <- exp(eta)
 wt <- expeta/((1+expeta)^2)
 wt*(t(fvec)%*%invinfo%*%fvec - t(f2vec)%*%invinfo2%*%f2vec)
}

#Calculate standardised variance at support points
stdvar1(c(-1,-1))
stdvar1(c(-1,1))
stdvar1(c(1,-1))
stdvar1(c(1,1))

#Draw contour plot of stdvar
exp1 <- expression(x[1])
exp2 <- expression(x[2])
x1 <- x2 <- seq(from=-1.1,to=1.1,by=0.005)
lx1 <- length(x1)
y <- rep(0,lx1^2)
kount <- 1
for (i in 1:lx1)
{
 for (j in 1:lx1)
 {
  y[kount] <- stdvar1(c(x1[i],x2[j]))
  kount <- kount + 1
 }
}
ymat <- matrix(y,lx1,lx1,byrow=T)
par(las=1)
contour(x1,x2,ymat,xlab=exp1,ylab=exp2, levels = seq(from=p1-0.6,to=p1+0.2,by=0.2),lwd = c(1,1,1,2,1),
   labcex=1)
lines(c(-1.05,1.05),c(1,1),lty=2)
lines(c(-1.05,1.05),c(-1,-1),lty=2)
lines(c(-1,-1),c(-1.05,1.05),lty=2)
lines(c(1,1),c(-1.05,1.05),lty=2)
points(c(-1,-1,1,1),c(-1,1,-1,1),pch=16,cex=2)