#  R implementation of the cocktail algorithm for D-optimal designs (approximate theory)
#  
#  output=cocktail(X, maxiter, small, alg)
#  
#  required input
#    X: design matrix (n by m)
#    (n: number of design points)
#    (m: dimension of parameter)
#  
#  optional input 
#    maxiter: maximum number of iterations (default=1000)
#    small: convergence threshold (default=1e-5)
#    alg: choice of algorithm (see below; default=4)
#  
#  output: a list of 
#    i) weights (1 by n), ii) bounds on the log determinant, iii) iteration count, and 
#    iv) convergence status (1 means convergence and 0 otherwise) 

cocktail=function(X, maxiter=1000, small=1e-5, alg=4){ 
  if(!is.matrix(X)) 
    X=as.matrix(X)
  n=dim(X)[1]
  m=dim(X)[2]
  if(n<m||m<2)
    stop("Wrong input sizes") 
  if(!any(alg==0:4)){
    warning("Unknown algorithm chosen; default used instead")
    alg=4
  } 
  if(alg==0)
    print("Using the multiplicative algorithm only")
  if(alg==1)
    print("Using the vertex exchange method (VEM)")
  if(alg==2)
    print("Using both the vertex direction method (VDM) and the multiplicative algorithm")
  if(alg==3){
    print("Using the cocktail algorithm:") 
    print("VDM + nearest neighbor exchange + multiplicative algorithm")
    print("The neighborhood structure is pre-specified")
  }
  if(alg==4){
    print("Using the cocktail algorithm:")
    print("VDM + nearest neighbor exchange + multiplicative algorithm")
    print("The nearest neighbors are computed on the fly")
  } 

  # starting values  
  if(alg==0)
    w=rep(1, n)  
  else{
    w=rep(0, n)
    w[sample(1:n, min(2*m, n))]=1
  } 
  w=w/sum(w)  # weights should sum to one 
  index=1:n
  supp=index[w>0]  # support points 
  Y=matrix(nrow=n, ncol=m)
  for(j in 1:m)
    Y[supp, j]=X[supp, j]*w[supp]
  M=t(Y[supp, ])%*%X[supp, ]  # information matrix 
  deri=rep(0, n)

  for(iter in 1:maxiter){
    P=solve(M, t(X))
    for(j in 1:n)
      deri[j]=sum(X[j,]*P[,j])/m 
    max.d=max(deri)
    if(max.d<1+small)  # testing convergence
      break 
      
    if(alg==1){  # VEM step
      j2=max(index[max.d==deri]) 
      min.d=min(deri[w>0]) 
      j1=min(index[w>0&min.d==deri]) 
      d=X[c(j1, j2), ]%*%P[,c(j1, j2)] 
      stepsize=0.5*(d[1,1]-d[2,2])/(d[1,1]*d[2,2]-d[1,2]*d[2,1]+1e-200) 
      stepsize=max(min(stepsize, w[j2]), -w[j1]) 
      w[j1]=w[j1]+stepsize 
      w[j2]=w[j2]-stepsize 
      M=M+stepsize*X[j1,]%*%t(X[j1,])-stepsize*X[j2,]%*%t(X[j2,])
    }

    if(alg>=2){  # VDM step
      j2=max(index[max.d==deri]) 
      stepsize=(max.d-1)/(m*max.d-1) 
      w=(1-stepsize)*w 
      w[j2]=w[j2]+stepsize 
      M=(1-stepsize)*M+stepsize*X[j2,]%*%t(X[j2,]) 
    }

    if(alg>=3){  # nearest neighbor exchanges
      supp=index[w>0] 
      for(j in 1:(length(supp)-1)){
        j1=supp[j] 
        j2=supp[j+1] 
        if(alg==4){  # nearest neighbor computed on the fly
          near=Inf 
          for(k in (j+1):length(supp)){
            diff=sum(abs(X[j1,]-X[supp[k],]))
            if(near>diff){
              near=diff
              j2=supp[k]
            }
          }
        } 
        d=X[c(j1, j2), ]%*%solve(M, t(X[c(j1, j2), ])) 
        stepsize=0.5*(d[1,1]-d[2,2])/(d[1,1]*d[2,2]-d[1,2]*d[2,1]+1e-200)
        stepsize=max(min(stepsize, w[j2]), -w[j1])
        w[j1]=w[j1]+stepsize 
        w[j2]=w[j2]-stepsize 
        M=M+stepsize*X[j1,]%*%t(X[j1,])-stepsize*X[j2,]%*%t(X[j2,])
      }
    }

    if(alg!=1){  # multiplicative step
      if(alg==0)
        w=w*deri 
      if(alg>=2){
        supp=index[w>0]
        P=solve(M, t(X[supp, ]))
        for(j in 1:length(supp))
          w[supp[j]]=w[supp[j]]*sum(X[supp[j],]*P[,j])/m
      }
      w[supp]=w[supp]/sum(w[supp])  # redundant (numerical safeguard)
      for(j in 1:m)
        Y[supp, j]=X[supp, j]*w[supp]
      M=t(Y[supp, ])%*%X[supp, ]
    }
  }  # end for

  if(iter==maxiter) 
    cat("No convergence in", maxiter, "iterations\n")
  else 
    cat("Convergence in", iter, "iterations\n") 
  lower=determinant(M, logarithm=T)[[1]][1]
  upper=lower+m*(max.d-1)
  out=list(w, c(lower, upper), iter, iter!=maxiter) 
  names(out)=c("weights", "bound.criterion", "iteration.count", "converge") 
  out 
}