# Numerical study of AQuA (preprint at: https://arxiv.org/abs/1801.09124)
# Script in the language of R (https://www.r-project.org/)
# Requires packages gurobi, plyr, matrixcalc, Matrix, MASS
# To get the gurobi solver, see http://www.iam.fmph.uniba.sk/ospm/Harman/design/gurobi_inst.txt 
# Before running this study, load functions in AQuA_OptimalDesign.txt
# available at http://www.iam.fmph.uniba.sk/ospm/Harman/design/


##################
# Subsection 6.1 #
##################

library(gurobi)
library(MASS)

# Compute the matrix of all 64 possible regressors

sbw6 <- Fx.cube(~x1+x2+x3+x4+x5+x6-1, lower=rep(0, 6))
Ns <- 6:24; nN <- length(Ns)

# Compute the D-optimum exact values by the KL algorithm

opt.RCKL6D <- rep(0, nN)
des.RCKL6D <- matrix(0, ncol=nrow(sbw6), nrow=nN)
for(iN in 1:nN) {
   res <- od.KL(sbw6, Ns[iN], crit="D", t.max=60)
   opt.RCKL6D[iN] <- res$Phi.best
   des.RCKL6D[iN,] <- res$w.best
}

# Note: Analytic formulas for the optimal designs are in Neubauer et al. (2000)

# Apply the REX algorithm to compute the D-optimal approximate design
app6D <- od.AA(sbw6, crit="D", eff=0.99999999999)

effs <- c(0.95, 1.0); neff <- length(effs); nruns <- 3
t.max <- 1000; mat <- array(0, dim=c(nN, neff, nruns)) 
t.act.IQP6Dp <- eff.IQP6Dp <- mat
t.act.IQP6Dm <- eff.IQP6Dm <- mat

pert.D <- function(Fx, w.opt, eff) {
  # Find a random design with D-efficiency eff (approximately)
  # w.opt is the optimal design

  n <- nrow(Fx); m <- ncol(Fx)
  D.opt <- optcrit(Fx, w.opt, crit="D")

  # Generate a random design w.dir with efficiency 0
  ind.dir <- sample(1:n, m-1)
  w.dir <- rep(0, n); w.dir[ind.dir] <- 1/(m-1)

  # Slowly move in the direction of w.dir to get a design with required efficiency
  for(i in 0:10000) {
     del <- i/10000
     w.pert <- (1-del)*w.opt + del*w.dir        
     eff.pert <- optcrit(Fx, w.pert, crit="D") / D.opt 
     if(eff.pert < eff + 1e-12) break
  }

  return(list(w.pert=w.pert, eff.pert=eff.pert))
}

for(irun in 1:nruns) {

   # For the plot of a compact information during the computation
   plot(c(0.90, 1.05), c(5, 25), type="n", xlab="efficiency", ylab="N", main=paste("Model: sbw6; criterion: D; ", irun))

   for(iN in 1:nN) {
      for(ieff in 1:neff) { 
         res.pert <- pert.D(sbw6, app6D$w.best, effs[ieff])

         res.comp <- od.AQuA(sbw6, Ns[iN], M.anchor=infmat(sbw6, Ns[iN]*res.pert$w.pert), crit="D", ver="+", conic=FALSE, t.max=t.max)
         eff.IQP6Dp[iN, ieff, irun] <- res.comp$Phi.best/opt.RCKL6D[iN]
         t.act.IQP6Dp[iN, ieff, irun] <- res.comp$t.act
         text(effs[ieff]-0.01, Ns[iN], paste(round(eff.IQP6Dp[iN, ieff, irun], 4), "(+)"), cex=0.7, col="darkblue") 

         res.comp <- od.AQuA(sbw6, Ns[iN], M.anchor=infmat(sbw6, Ns[iN]*res.pert$w.pert), crit="D", ver="-", conic=FALSE, t.max=t.max)
         eff.IQP6Dm[iN, ieff, irun] <- res.comp$Phi.best/opt.RCKL6D[iN]
         t.act.IQP6Dm[iN, ieff, irun] <- res.comp$t.act
         text(effs[ieff]+0.01, Ns[iN], paste(round(eff.IQP6Dm[iN, ieff, irun], 4), "(-)"), cex=0.7, col="darkblue") 
      }
   }
}

# Compute the medians

med <- matrix(0, nrow=nN, ncol=neff) 
t.act.IQP6Dp.med <- eff.IQP6Dp.med <- med
t.act.IQP6Dm.med <- eff.IQP6Dm.med <- med

for(iN in 1:nN) {
   for(ieff in 1:neff) { 
      t.act.IQP6Dp.med[iN, ieff] <- median(t.act.IQP6Dp[iN, ieff,])
      t.act.IQP6Dm.med[iN, ieff] <- median(t.act.IQP6Dm[iN, ieff,])
      eff.IQP6Dp.med[iN, ieff] <- median(eff.IQP6Dp[iN, ieff,])
      eff.IQP6Dm.med[iN, ieff] <- median(eff.IQP6Dm[iN, ieff,])
   }
}

# Compute the whiskers

whi <- matrix(0, nrow=nN, ncol=neff) 
lt.act.IQP6D.whil <- lt.act.IQP6D.whiu <- whi
eff.IQP6D.whil <- eff.IQP6D.whiu <- whi

for(iN in 1:nN) {
   for(ieff in 1:neff) { 
      lt.act.IQP6D.whil[iN, ieff] <- log(min(c(t.act.IQP6Dp[iN, ieff,], t.act.IQP6Dm[iN, ieff,])), base=10)
      lt.act.IQP6D.whiu[iN, ieff] <- log(max(c(t.act.IQP6Dp[iN, ieff,], t.act.IQP6Dm[iN, ieff,])), base=10)
      eff.IQP6D.whil[iN, ieff] <- min(c(eff.IQP6Dp[iN, ieff,], eff.IQP6Dm[iN, ieff,]))
      eff.IQP6D.whiu[iN, ieff] <- max(c(eff.IQP6Dp[iN, ieff,], eff.IQP6Dm[iN, ieff,]))
   }
}

# Plot the graph of the efficiencies

par(mfrow=c(2, 2))
plot(Ns, eff.IQP6Dp.med[, 2], xlab="N", ylab="efficiency", ylim=c(0.90, 1), type="n", main="(a) D_sbw_eff100", cex.axis=0.8); grid()
for(iN in 1:nN) lines(c(Ns[iN], Ns[iN]), c(eff.IQP6D.whil[iN, 2], eff.IQP6D.whiu[iN, 2]), lwd=2, col="darkgray")
lines(Ns, eff.IQP6Dp.med[, 2], type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, eff.IQP6Dm.med[, 2], type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")

# Plot the graph of the computation times (vertical axis is log of the time)

plot(Ns, log(t.act.IQP6Dp.med[, 2], base=10), xlab="N", ylab="log(time)", ylim=c(-2, 3), type="n", main="(b) D_sbw_time100", cex.axis=0.8); grid()
for(iN in 1:nN) lines(c(Ns[iN], Ns[iN]), c(lt.act.IQP6D.whil[iN, 2], lt.act.IQP6D.whiu[iN, 2]), lwd=2, col="darkgray")
lines(Ns, log(t.act.IQP6Dp.med[, 2], base=10), type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, log(t.act.IQP6Dm.med[, 2], base=10), type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")

# Plot the graph of the efficiencies

plot(Ns, eff.IQP6Dp.med[, 1], xlab="N", ylab="efficiency", ylim=c(0.90, 1), type="n", main="(c) D_sbw_eff095", cex.axis=0.8); grid()
for(iN in 1:nN) lines(c(Ns[iN], Ns[iN]), c(eff.IQP6D.whil[iN, 1], eff.IQP6D.whiu[iN, 1]), lwd=2, col="darkgray")
lines(Ns, eff.IQP6Dp.med[, 1], type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, eff.IQP6Dm.med[, 1], type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")

# Plot the graph of the computation times (vertical axis is log of the time)

plot(Ns, log(t.act.IQP6Dp.med[, 1], base=10), xlab="N", ylab="log(time)", ylim=c(-2, 3), type="n", main="(d) D_sbw_time095", cex.axis=0.8); grid()
for(iN in 1:nN) lines(c(Ns[iN], Ns[iN]), c(lt.act.IQP6D.whil[iN, 1], lt.act.IQP6D.whiu[iN, 1]), lwd=2, col="darkgray") 
lines(Ns, log(t.act.IQP6Dp.med[, 1], base=10), type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, log(t.act.IQP6Dm.med[, 1], base=10), type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")


# Compute the A-optimum exact values by the KL heuristic

opt.RCKL6A <- rep(0, nN)
des.RCKL6A <- matrix(0, ncol=nrow(sbw6), nrow=nN)
for(iN in 1:nN) {
   res <- od.KL(sbw6, Ns[iN], crit="A", t.max=60)
   opt.RCKL6A[iN] <- res$Phi.best
   des.RCKL6A[iN,] <- res$w.best
}

# Apply the REX algorithm to compute the A-optimal approximate design
app6A <- od.AA(sbw6, crit="A", eff=0.99999999999)

effs <- c(0.95, 1.0); neff <- length(effs); nruns <- 3
t.max <- 1000; mat <- array(0, dim=c(nN, neff, nruns)) 
t.act.IQP6Ap <- eff.IQP6Ap <- mat
t.act.IQP6Am <- eff.IQP6Am <- mat

pert.A <- function(Fx, w.opt, eff) {
  # Find a random design with A-efficiency eff (approximately)
  # w.opt is the optimal design

  n <- nrow(Fx); m <- ncol(Fx)
  A.opt <- optcrit(Fx, w.opt, crit="A")

  # Generate a random design w.dir with efficiency 0
  ind.dir <- sample(1:n, m-1)
  w.dir <- rep(0, n); w.dir[ind.dir] <- 1/(m-1)

  # Slowly move in the direction of w.dir to get a design with required efficiency
  for(i in 0:10000) {
     del <- i/10000
     w.pert <- (1-del)*w.opt + del*w.dir        
     eff.pert <- optcrit(Fx, w.pert, crit="A") / A.opt 
     if(eff.pert < eff + 1e-12) break
  }

  return(list(w.pert=w.pert, eff.pert=eff.pert))
}

par(mfrow=c(1, 1))
for(irun in 1:nruns) {

   # For the plot of a compact information during the computation
   plot(c(0.90, 1.05), c(5, 25), type="n", xlab="efficiency", ylab="N", main=paste("Model: sbw6; criterion: A; ", irun))

   for(iN in 1:nN) {
      for(ieff in 1:neff) { 
         res.pert <- pert.A(sbw6, app6A$w.best, effs[ieff])

         res.comp <- od.AQuA(sbw6, Ns[iN], M.anchor=infmat(sbw6, Ns[iN]*res.pert$w.pert), crit="A", ver="+", conic=FALSE, t.max=t.max)
         eff.IQP6Ap[iN, ieff, irun] <- res.comp$Phi.best/opt.RCKL6A[iN]
         t.act.IQP6Ap[iN, ieff, irun] <- res.comp$t.act
         text(effs[ieff]-0.01, Ns[iN], paste(round(eff.IQP6Ap[iN, ieff, irun], 4), "(+)"), cex=0.7, col="darkblue") 

         res.comp <- od.AQuA(sbw6, Ns[iN], M.anchor=infmat(sbw6, Ns[iN]*res.pert$w.pert), crit="A", ver="-", conic=FALSE, t.max=t.max)
         eff.IQP6Am[iN, ieff, irun] <- res.comp$Phi.best/opt.RCKL6A[iN]
         t.act.IQP6Am[iN, ieff, irun] <- res.comp$t.act
         text(effs[ieff]+0.01, Ns[iN], paste(round(eff.IQP6Am[iN, ieff, irun], 4), "(-)"), cex=0.7, col="darkblue") 
      }
   }
}

# Compute the medians

med <- matrix(0, nrow=nN, ncol=neff) 
t.act.IQP6Ap.med <- eff.IQP6Ap.med <- med
t.act.IQP6Am.med <- eff.IQP6Am.med <- med

for(iN in 1:nN) {
   for(ieff in 1:neff) { 
      t.act.IQP6Ap.med[iN, ieff] <- median(t.act.IQP6Ap[iN, ieff,])
      t.act.IQP6Am.med[iN, ieff] <- median(t.act.IQP6Am[iN, ieff,])
      eff.IQP6Ap.med[iN, ieff] <- median(eff.IQP6Ap[iN, ieff,])
      eff.IQP6Am.med[iN, ieff] <- median(eff.IQP6Am[iN, ieff,])
   }
}

# Compute the whiskers

whi <- matrix(0, nrow=nN, ncol=neff) 
lt.act.IQP6A.whil <- lt.act.IQP6A.whiu <- whi
eff.IQP6A.whil <- eff.IQP6A.whiu <- whi

for(iN in 1:nN) {
   for(ieff in 1:neff) { 
      lt.act.IQP6A.whil[iN, ieff] <- log(min(c(t.act.IQP6Ap[iN, ieff,], t.act.IQP6Am[iN, ieff,])), base=10)
      lt.act.IQP6A.whiu[iN, ieff] <- log(max(c(t.act.IQP6Ap[iN, ieff,], t.act.IQP6Am[iN, ieff,])), base=10)
      eff.IQP6A.whil[iN, ieff] <- min(c(eff.IQP6Ap[iN, ieff,], eff.IQP6Am[iN, ieff,]))
      eff.IQP6A.whiu[iN, ieff] <- max(c(eff.IQP6Ap[iN, ieff,], eff.IQP6Am[iN, ieff,]))
   }
}

# Plot the graph of the efficiencies

par(mfrow=c(2, 2))
plot(Ns, eff.IQP6Ap.med[, 2], xlab="N", ylab="efficiency", ylim=c(0.90, 1), type="n", main="(a) A_sbw_eff100", cex.axis=0.8); grid()
for(iN in 1:nN) {
   lines(c(Ns[iN], Ns[iN]), c(eff.IQP6A.whil[iN, 2], eff.IQP6A.whiu[iN, 2]), lwd=2, col="darkgray")
} 
lines(Ns, eff.IQP6Ap.med[, 2], type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, eff.IQP6Am.med[, 2], type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")

# Plot the graph of the computation times (vertical axis is log of the time)

plot(Ns, log(t.act.IQP6Ap.med[, 2], base=10), xlab="N", ylab="log(time)", ylim=c(-2, 3), type="n", main="(b) A_sbw_time100", cex.axis=0.8); grid()
for(iN in 1:nN) {
   lines(c(Ns[iN], Ns[iN]), c(lt.act.IQP6A.whil[iN, 2], lt.act.IQP6A.whiu[iN, 2]), lwd=2, col="darkgray")
} 
lines(Ns, log(t.act.IQP6Ap.med[, 2], base=10), type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, log(t.act.IQP6Am.med[, 2], base=10), type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")

# Plot the graph of the efficiencies

plot(Ns, eff.IQP6Ap.med[, 1], xlab="N", ylab="efficiency", ylim=c(0.90, 1), type="n", main="(c) A_sbw_eff095", cex.axis=0.8); grid()
for(iN in 1:nN) {
   lines(c(Ns[iN], Ns[iN]), c(eff.IQP6A.whil[iN, 1], eff.IQP6A.whiu[iN, 1]), lwd=2, col="darkgray")
} 
lines(Ns, eff.IQP6Ap.med[, 1], type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, eff.IQP6Am.med[, 1], type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")

# Plot the graph of the computation times (vertical axis is log of the time)

plot(Ns, log(t.act.IQP6Ap.med[, 1], base=10), xlab="N", ylab="log(time)", ylim=c(-2, 3), type="n", main="(d) A_sbw_time095", cex.axis=0.8); grid()
for(iN in 1:nN) {
   lines(c(Ns[iN], Ns[iN]), c(lt.act.IQP6A.whil[iN, 1], lt.act.IQP6A.whiu[iN, 1]), lwd=2, col="darkgray")
} 
lines(Ns, log(t.act.IQP6Ap.med[, 1], base=10), type="b", pch=24, cex=0.7, bg="darkblue", col="darkblue")
lines(Ns, log(t.act.IQP6Am.med[, 1], base=10), type="b", pch=25, cex=0.7, bg="darkblue", col="darkblue")


##################
# Subsection 6.2 #
##################

plot.mix <- function (w, k, pch, col, bg, main) {
   # Plot of mixture designs with 3 components on a discretized simplex
   # w ... design vector
   # k ... numer of (each) factor levels
   # pch ... plot character
   # col ... design points color
   # main ... main title of the figure

   trans <- function(a) {
      x1 <- a[2] + a[3]/2
      x2 <- sqrt(3)/2 * a[3]
      return(c(x1, x2))
   }

   plot(c(0, 1), c(0, sqrt(3)/2), type="n", xlab="", ylab="", main=main, axes=FALSE, asp=1)
   A <- trans(c(1,0,0)); B <- trans(c(0,1,0)); C <- trans(c(0,0,1)) 
   lines(c(A[1], B[1]), c(A[2], B[2]))
   lines(c(B[1], C[1]), c(B[2], C[2]))
   lines(c(C[1], A[1]), c(C[2], A[2]))

   Fx.aux <- Fx.simp(~x1+x2+x3-1, n.levels=k); n <- nrow(Fx.aux)

   for(i in 1:n) {
      x <- trans(Fx.aux[i,])
      points(x[1], x[2], pch=19, cex=0.1, col=rgb(0.8,0.8,0.8))
   }

   for(i in 1:n) {
      x <- trans(Fx.aux[i,])
      if(1.5*sqrt(w[i])>0.1) {
         points(x[1], x[2], pch=pch, cex=1.5*sqrt(w[i]), col=col, bg=bg)
      }
   }
}


# Procedure for generating the model and the constraints

bridge.Scheffe <- function(N, d, type="quad") {

   if(type=="quad") {
      Fx <- Fx.simp(~x1+x2+x3+I(x1*x2)+I(x1*x3)+I(x2*x3)-1, d)
   } else {
      Fx <- Fx.simp(~x1+x2+x3-1, d)
   }
   m <- ncol(Fx); n <- nrow(Fx)
  
   b1 <- rep(1, d+1)  
   A1 <- matrix(0, nrow=d+1, ncol=n)  
   I1 <- c(0, 0, 1)
   for (j in 3:d) I1 <- c(I1, sum(1:(j-1)))
   for (i in 1:d) A1[i, 1:(d-i+1) + (i-1)*d-I1[i]] <- 1

   A1[d+1, ] <- 1; b1[d+1] <- N
 
   At1 <- matrix(0, nrow=2*nrow(Fx), ncol=3)
   for(i in 1:nrow(Fx)){
      riadky <- i
      p <- Fx[i, c(1,2,3)]
      for(j in 1:nrow(Fx)){
         if(all(Fx[j,1:3]==p) || all(Fx[j,1:3]==c(p[2],p[3],p[1])) || all(Fx[j,1:3]==c(p[3],p[1],p[2]))) riadky <- c(riadky, j)
      }
      l <- length(unique(riadky))
      At1[i, 1:l] <- sort(unique(riadky))
   }

   for(i in 1:nrow(Fx)){
      riadky <- i; p <- Fx[i, c(1,2,3)]
      for(j in 1:nrow(Fx)){
	 if(all(Fx[j,1:3]==c(p[1],p[3],p[2])) || all(Fx[j,1:3]==c(p[2],p[1],p[3])) || all(Fx[j,1:3]==c(p[3],p[2],p[1]))) riadky <- c(riadky, j)
      }
      riadky <- riadky[-1]
      l <- length(unique(riadky))
      At1[i+nrow(Fx), 1:l] <- sort(unique(riadky))
   }

   At1 <- At1[!duplicated(At1), ]
   row_sub <- apply(At1, 1, function(row) all(row !=0 ))
   At1 <- At1[row_sub,]

   A3 <- matrix(0, nrow=2*nrow(At1), ncol=nrow(Fx))
   for(i in 1:nrow(At1)){
      A3[(i-1)*2+1, At1[i, 1]] <- 1
      A3[(i-1)*2+1, At1[i, 2]] <- -1
      A3[(i-1)*2+2, At1[i, 1]] <- 1
      A3[(i-1)*2+2, At1[i, 3]] <- -1
   }
   b3 <- rep(0, nrow(A3))

   return(list(Fx=Fx, b1=b1, A1=A1, b3=b3, A3=A3))
}

library(matrixcalc); library(Matrix)

N <- 18; d <- 41
br.s <- bridge.Scheffe(N, d)

# Compute all required mixture designs

br1841appD <- od.SOCP(br.s$Fx, br.s$b1, br.s$A1, b3=br.s$b3, A3=br.s$A3, t.max=100)
br1841iqpDcon <- od.AQuA(br.s$Fx, br.s$b1, br.s$A1, b3=br.s$b3, A3=br.s$A3, bin=TRUE, conic=TRUE, t.max=100)
br1841iqpD <- od.AQuA(br.s$Fx, br.s$b1, br.s$A1, b3=br.s$b3, A3=br.s$A3, bin=TRUE, conic=FALSE, t.max=100)
br1841misocpD <- od.MISOCP(br.s$Fx, br.s$b1, br.s$A1, b3=br.s$b3, A3=br.s$A3, bin=TRUE, t.max=100)
br1841appI <- od.SOCP(br.s$Fx, br.s$b1, br.s$A1, b3=br.s$b3, A3=br.s$A3, crit="I", t.max=100)
br1841iqpIcon <- od.AQuA(br.s$Fx, br.s$b1, br.s$A1, bin=TRUE, conic=TRUE, crit="I", t.max=100, b3=br.s$b3, A3=br.s$A3)
br1841iqpI <- od.AQuA(br.s$Fx, br.s$b1, br.s$A1, bin=TRUE, conic=FALSE, crit="I", t.max=100, b3=br.s$b3, A3=br.s$A3)
br1841misocpI <- od.MISOCP(br.s$Fx, br.s$b1, br.s$A1, bin=TRUE, crit="I", t.max=100, b3=br.s$b3, A3=br.s$A3)

# Plot the designs

windows()
par(mfrow=c(2, 2))
plot.mix(br1841appD$w.best, d, 21, "black", "black", "(a) D_mix_SOCP")
plot.mix(br1841misocpD$w.best, d, 21, "black", "red", "(b) D_mix_MISOCP")
plot.mix(br1841iqpD$w.best, d, 21, "black", "darkblue", "(c) D_mix_AQuA_IQP")
plot.mix(br1841iqpDcon$w.best, d, 21, "black", "lightblue", "(d) D_mix_AQuA_MICQP")

windows()
par(mfrow=c(2, 2))
plot.mix(br1841appI$w.best, d, 21, "black", "black", "(a) I_mix_SOCP")
plot.mix(br1841misocpI$w.best, d, 21, "black", "red", "(b) I_mix_MISOCP")
plot.mix(br1841iqpI$w.best, d, 21, "black", "darkblue", "(c) I_mix_AQuA_IQP")
plot.mix(br1841iqpIcon$w.best, d, 21, "black", "lightblue", "(d) I_mix_AQuA_MICQP")


##################
# Subsection 6.3 #
##################

# Original dataset wine https://www.kaggle.com/zynicide/wine-reviews#winemag-data-130k-v2.csv
# Cleaning of wine: remove all row that have missing or incorrect entries in any of the following variables:
# country, description, points, price, province, title, variety, winary. Then remove duplicate rows.

# Note: The following file has almost 40MB; it make take a while to download

wine <- read.csv("http://www.iam.fmph.uniba.sk/ospm/Harman/data/wine_clean_semicolon.csv", sep=";")
wine$points <- wine$points - 90
wine$cind <- as.numeric(wine$country)
attach(wine); n <- nrow(wine)


D.val <- matrix(0, nrow=10, ncol=4)

for(k in 1:10) {
  sel <- c();
  for(i in 1:42) {
    cnk <- (1:n)[cind==i]
    ord <- order(-points[cnk])
    sel <- c(sel, cnk[ord[1]])
  }
   
  sel <- c(sel, sample(setdiff(1:n, sel), 1958))
  A3 <- matrix(0, nrow=42, ncol=n); for(j in 1:42) A3[j,] <- as.numeric(cind==j)
  Fx <- cbind(1, log(price, 10), points)

  res.mix <- od.SOCP(Fx[sel,], b1=1000, A1=matrix(price[sel], nrow=1), b2=0, A2=matrix(points[sel], nrow=1), b3=rep(1, 42), A3=A3[,sel], bin=TRUE, t.max=60)
  D.val[k, 1] <- res.mix$Phi.best

  for(j in 2:4) {
    res.mix <- od.AQuA(Fx, b1=1000, A1=matrix(price, nrow=1), b2=0, A2=matrix(points, nrow=1), b3=rep(1,42), A3=A3, M.anchor=res.mix$M.best, conic=TRUE, bin=TRUE, t.max=60)
   D.val[k, j] <- res.mix$Phi.best
  }
}

# Plot the D-optimal subsample

windows()
pp <- cbind(points, price); pp <- unique(pp)
nn <- nrow(pp); nums <- rep(0, nn)
for(i in 1:nn) nums[i] <- sum(points==pp[i,1] & price==pp[i,2])  
plot(pp[,1]+90, log(pp[,2], base=10), pch=19, cex=0.1*sqrt(nums), col="darkgray", xlab="points", ylab="log(price)", main="D_wine_AQuA")
points(points[res.mix$supp]+90, log(price[res.mix$supp], base=10), pch=21, bg="darkblue", col="black")

ingraph <- function (x, Y, xmin, xmax, ymin, ymax, eps) 
{
  n <- length(x); mxy <- max(Y)
  x <- (x-min(x))/(max(x)-min(x))*(xmax-xmin)+xmin
  Y <- (Y-0*min(Y))/(max(Y)-0*min(Y))*(ymax-ymin)+ymin
  for(j in 1:nrow(Y)) 
     for(i in 1:(n-1))
        lines(x[i:(i+1)], Y[j, i:(i+1)], type="b", pch=19)
  lines(c(xmin, xmax), c(ymin, ymin))
  lines(c(xmin, xmin), c(ymin, ymax))
  text(xmin-eps, ymin, "0")
  text(xmin-eps, ymax, round(mxy))
}

ingraph(1:4, D.val[1:5,], 81.5, 86.5, 2.7, 3.5, 0.7)


I.val <- matrix(0, nrow=10, ncol=4)

for(k in 1:10) {
  sel <- c();
  for(i in 1:42) {
    cnk <- (1:n)[cind==i]
    ord <- order(-points[cnk])
    sel <- c(sel, cnk[ord[1]])
  }
   
  sel <- c(sel, sample(setdiff(1:n, sel), 1958))
  A3 <- matrix(0, nrow=42, ncol=n); for(j in 1:42) A3[j,] <- as.numeric(cind==j)
  Fx <- cbind(1, log(price, 10), points)
  FxA <- Fx.ItoA(Fx)

  resA.mix <- od.SOCP(FxA[sel,], b1=1000, A1=matrix(price[sel], nrow=1), b2=0, A2=matrix(points[sel], nrow=1), b3=rep(1, 42), A3=A3[,sel], crit="A", bin=TRUE, t.max=60)
  I.val[k, 1] <- optcrit(Fx[sel,], resA.mix$w.best, crit="I", eta=rep(1/length(sel), length(sel)))

  for(j in 2:4) {
    resA.mix <- od.AQuA(FxA, b1=1000, A1=matrix(price, nrow=1), b2=0, A2=matrix(points, nrow=1), b3=rep(1,42), A3=A3, M.anchor=resA.mix$M.best, crit="A", conic=TRUE, bin=TRUE, t.max=60)
    I.val[k, j] <- optcrit(Fx, resA.mix$w.best, crit="I", eta=rep(1/n, n))
  }
}

# Plot the I-optimal subsample

windows()
pp <- cbind(points, price); pp <- unique(pp)
nn <- nrow(pp); nums <- rep(0, nn)
for(i in 1:nn) nums[i] <- sum(points==pp[i,1] & price==pp[i,2])  
plot(pp[,1]+90, log(pp[,2], base=10), pch=19, cex=0.1*sqrt(nums), col="darkgray", xlab="points", ylab="log(price)", main="I_wine_AQuA")
points(points[resA.mix$supp]+90, log(price[resA.mix$supp], base=10), pch=21, bg="darkblue", col="black")

ingraph <- function (x, Y, xmin, xmax, ymin, ymax, eps) 
{
  n <- length(x); mxy <- max(Y)
  x <- (x-min(x))/(max(x)-min(x))*(xmax-xmin)+xmin
  Y <- (Y-0*min(Y))/(max(Y)-0*min(Y))*(ymax-ymin)+ymin
  for(j in 1:nrow(Y)) 
     for(i in 1:(n-1))
        lines(x[i:(i+1)], Y[j, i:(i+1)], type="b", pch=19)
  lines(c(xmin, xmax), c(ymin, ymin))
  lines(c(xmin, xmin), c(ymin, ymax))
  text(xmin-eps, ymin, "0")
  text(xmin-eps, ymax, round(mxy))
}

ingraph(1:4, I.val[1:5,], 81.5, 86.5, 2.7, 3.5, 0.7)

#### END ####