# DDOA

m = c() # True mean
  
mbest = which.max(m) # Index of the true best alternative

M = length(m) # Number of alternatives

r = 1000 # Number of replications

n = 25000 # Total budget

sn = seq(0, 0, length = M) # Count the number of allocated samples over r replications

gxstar = seq(-0.5, -0.5, length = M) # \gamma^n_{x^{*, n}, x}

gx = seq(0.5, 0.5, length = M) # \gamma^n_{x, x^{*, n}}

u = seq(0, 0, length = M) # u^n_{x^{*, n}, x}

ixstar = seq(0, 0, length = M) # I^n_{x^{*, n}} (\gamma^n_{x^{*, n}, x})

ix = seq(0, 0, length = M) # I^n_x (\gamma^n_{x, x^{*, n}})

v = seq(0, 0, length = M) # \Gamma^n_{x^{*, n}, x}

ics = seq(0, 0, length = n) # Count the number of incorrect selection

initialsize = 2 # Naive initialization with two constants 0 and 1

warmup = 10 # Warmup size

for(i in 1:r){
  
  N = seq(initialsize, initialsize, length = M) # Count the number of allocated samples in the current replication
  
  wmatrix = matrix(0, nrow = M, ncol = n+initialsize) # To save allocated samples
  
  wmatrix[, 1] = 1
  
  wmatrix[, 2] = 0
  
  t = seq(0.5, 0.5, length = M) # Sample mean
  
  csampled = 1 # 1 if the current best is sampled, 0 otherwise
  
  for(j in 1:(M*warmup)){ # Warmup
    
    c = (j-1)%%M + 1
    
    N[c] = N[c] + 1
    
    w = samplegenerator(m[c], ...) # Generate a sample for alternative c
    t[c] = (t[c]*(N[c]-1)+w)/N[c] # Update sample mean
    
    wmatrix[c, N[c]] = w
    
    d = which.max(t) # The current best alternative
    
    if(d != mbest) # Check if incorrect selection is made
      ics[j] = ics[j] + 1
  }
  
  clast = 0 # The current best at the last stage
  
  for(j in (M*warmup+1):n){
    
    c = which.max(t) # The current best
    
    uc = 1
    
    for(k in 1:M){
      if(k != c & t[k] == t[c])
        uc = 0
    }
    
    if(uc == 0){
      for(k in 1:M){
        if(k != c & t[k] == t[c] & N[k] < N[c])
          c = k
      }
      
      N[c] = N[c] + 1
      
      w = samplegenerator(m[c], ...)
      t[c] = (t[c]*(N[c]-1)+w)/N[c]
      
      wmatrix[c, N[c]] = w
      
      d = which.max(t)
      
      if(d != mbest)
        ics[j] = ics[j] + 1
      
      csampled = 1
      clast = c
    }
    
    if(uc == 1){
      
      if(csampled == 0 & c == clast){ # Only one alternative needs rate update
        ixstar = ixstarlast
        ix = ixlast
        v = vlast
        
        k = c1last
        
        rn = N[c]/N[k]
        
        wmaxc = max(wmatrix[c, 1:N[c]])
        wmaxk = max(wmatrix[k, 1:N[k]])
        
        f = function(x){
          ordc1 = sum(wmatrix[c, 1:N[c]]*exp(x*(wmatrix[c, 1:N[c]]-wmaxc)))
          ordc0 = sum(exp(x*(wmatrix[c, 1:N[c]]-wmaxc)))
          
          ordk1 = sum(wmatrix[k, 1:N[k]]*exp(-rn*x*(wmatrix[k, 1:N[k]]-wmaxk)))
          ordk0 = sum(exp(-rn*x*(wmatrix[k, 1:N[k]]-wmaxk)))
          
          return(ordc1/ordc0 - ordk1/ordk0)
        }
        
        lb = 2*gxstar[k]
        
        gxstar[k] = uniroot(f, c(lb, 0), extendInt = "upX")$root
        
        if(gxstar[k] == 0){
          gxstar[k] = 1e-6
        }
        
        gx[k] = -rn*gxstar[k]
        
        u[k] = sum(wmatrix[c, 1:N[c]]*exp(gxstar[k]*wmatrix[c, 1:N[c]]))/sum(exp(gxstar[k]*wmatrix[c, 1:N[c]]))
        
        ixstar[k] = gxstar[k]*u[k] - log(1/N[c]*sum(exp(gxstar[k]*wmatrix[c, 1:N[c]])))
        
        ix[k] = gx[k]*u[k] - log(1/N[k]*sum(exp(gx[k]*wmatrix[k, 1:N[k]])))
        
        v[k] = N[c]*ixstar[k] + N[k]*ix[k]
      }
      
      else{ # All alternatives need rate update
        for(k in 1:M){
          if(k != c){
            rn = N[c]/N[k]
            
            wmaxc = max(wmatrix[c, 1:N[c]])
            wmaxk = max(wmatrix[k, 1:N[k]])
            
            f = function(x){
              ordc1 = sum(wmatrix[c, 1:N[c]]*exp(x*(wmatrix[c, 1:N[c]]-wmaxc)))
              ordc0 = sum(exp(x*(wmatrix[c, 1:N[c]]-wmaxc)))
              
              ordk1 = sum(wmatrix[k, 1:N[k]]*exp(-rn*x*(wmatrix[k, 1:N[k]]-wmaxk)))
              ordk0 = sum(exp(-rn*x*(wmatrix[k, 1:N[k]]-wmaxk)))
              
              return(ordc1/ordc0 - ordk1/ordk0)
            }
            
            lb = 2*gxstar[k]
            
            gxstar[k] = uniroot(f, c(lb, 0), extendInt = "upX")$root
            
            if(gxstar[k] == 0){
              gxstar[k] = 1e-6
            }
            
            gx[k] = -rn*gxstar[k]
            
            u[k] = sum(wmatrix[c, 1:N[c]]*exp(gxstar[k]*wmatrix[c, 1:N[c]]))/sum(exp(gxstar[k]*wmatrix[c, 1:N[c]]))
            
            ixstar[k] = gxstar[k]*u[k] - log(1/N[c]*sum(exp(gxstar[k]*wmatrix[c, 1:N[c]])))
            
            ix[k] = gx[k]*u[k] - log(1/N[k]*sum(exp(gx[k]*wmatrix[k, 1:N[k]])))
            
            v[k] = N[c]*ixstar[k] + N[k]*ix[k]
          }
        }
      }
      
      ss = 0
      
      for(k in 1:M){
        if(k != c)
          ss = ss + ixstar[k]/ix[k]
      }
      
      if(ss > 1){
        N[c] = N[c] + 1
        
        w = samplegenerator(m[c], ...)
        t[c] = (t[c]*(N[c]-1)+w)/N[c]
        
        wmatrix[c, N[c]] = w
        
        d = which.max(t)
        
        if(d != mbest)
          ics[j] = ics[j] + 1
        
        csampled = 1
        clast = c
      }
      
      else{
        v[c] = max(abs(v)) + 1
        
        c1 = which.min(v) # The alternative with the lowest rate
        
        N[c1] = N[c1] + 1
        
        w = samplegenerator(m[c1], ...)
        t[c1] = (t[c1]*(N[c1]-1)+w)/N[c1]
        
        wmatrix[c1, N[c1]] = w
        
        d = which.max(t)
        
        if(d != mbest)
          ics[j] = ics[j] + 1
        
        csampled = 0
        clast = c
        c1last = c1 # The alternative sampled at the last stage
        
        # Rates at the last stage
        ixstarlast = ixstar
        ixlast = ix
        vlast = v
      }
    }
  }
  
  sn = sn + N
  
  print(paste(100*i/r, "%")) # Report progress
}

prop = sn/r/(n+M*initialsize) # Allocation proportion over r replications

ics = ics/r # 1-PCS over r replications


# Plot PCS

par(ps=12)

par(mar = c(4, 4, 1, 1))

plot(ics, type="l", xlab="n", ylab="1-PCS", 
     col="seagreen4", ylim=c(0, 1), cex.axis=1, lwd=0.5)

paxis = seq(2500, n, by = 5000) # Positions to add point symbols

points(paxis, ics[paxis], pch=1, cex=1.2, col="seagreen4") # Add point symbols
