#Program 19
# Looks for D-optimal design for Bayesian design for logit(pi) = beta0 + beta1*x, with -10 <= x <= 10

f <- function(x)
{
 f <- c(1,x)
 f
}

infomat <- function(betavec,xvec,deswts)
{
 xmat <- t(matrix(xvec,s,m))
 fxmat <- apply(xmat,2,fx)
 etavec <- as.vector(t(betavec)%*%fxmat)
 expetavec <- exp(etavec)
 modelwtvec <- expetavec/((1+expetavec)^2)
 infomat <- fxmat %*% diag(deswts*modelwtvec) %*% t(fxmat)
 infomat
}

combine <- function(values)
{ #This function calculates the NEGATIVE of the sum of the log(determinant(infomat)) over all the potential beta vectors
 xvec <- 10*cos(pi*values[1:s])
 temp <- (values[(s+1):(2*s)])^2
 deltavec <- temp/sum(temp)
 sum <-  0
 for (m in 1:h)
 {
  betavec <- betamat[m,]
  logdetinfo <- log(det(infomat(betavec,xvec,deltavec)))
  sum <- sum - psivec[m]*logdetinfo
 }
 sum
}

betamat <- matrix(c(-0.2,-0.2,0.2,0.2,0.8,1.2,0.8,1.2),4,2)#matrix(c(0,1),1,2)#
psivec <- c(0.1,0.2,0.3,0.4)#c(1)#
h <- (dim(betamat))[1]
s <- 2
m <- 1
p <- 2

nsims <- 10000
min <- 100 
for (isim in 1:nsims)
{
 values <- runif(2*s)
 out <- optim(values,combine,NULL,method="Nelder-Mead")
 temp <- out$val
 if (temp < min) {min <- temp
 design <- out$par}
}
points <- 10*cos(pi*(design[1:s]))
weights <- (design[(s+1):(2*s)])^2
weights <- weights/sum(weights)
display <- rbind(points,weights)
display
cat("min value of (-logdeterminant)",min,"\n")