#                   NESTED SAMPLING MAIN PROGRAM
# (GNU General Public License software, (C) Sivia and Skilling 2006)
# translated to R by Issac Trotts in 2007
#
# note [1]: R seems to do deep copying of its lists, so the copying part
# of this code is easy.

uniform <- function() runif(1, 0.0,1.0)

# given x=log(a) and y=log(b), return log(a+b)
log.plus <- function(x,y) {
    if(x>y) x+log(1+exp(y-x))
    else    y+log(1+exp(x-y))
}

DBL.MAX <- 1.79769e+308

nested.sampling <- function(n, rprior, explore, max.iter=100) {
  Samples <- list()                  # Weighted samples from posterior
  H    <- 0.0                        # Information, initially 0
  logZ <- -DBL.MAX                   # ln(Evidence Z, initially 0)
  
  ## Set prior objects
  Obj <- list()                         # Collection of n objects
  for(i in 1:n) { Obj[[i]] <- rprior() }
  ## Outermost interval of prior mass
  logwidth <- log(1.0 - exp(-1.0 / n))  # ln(width in prior mass)
  ## NESTED SAMPLING LOOP _____________________________________________
  for(nest in 1:max.iter) {
    ## Find worst object in collection, with Weight <- width * Likelihood
    worst <- 1
    for(i in 2:n) {
      if( Obj[[i]]$logL < Obj[[worst]]$logL )  {
        worst <- i
      }
    }
    Obj[[worst]]$logWt <- logwidth + Obj[[worst]]$logL
    ## Update Evidence Z and Information H
    logZnew <- log.plus(logZ, Obj[[worst]]$logWt)
    H <- exp(Obj[[worst]]$logWt - logZnew) * Obj[[worst]]$logL + 
      exp(logZ - logZnew) * (H + logZ) - logZnew;
    logZ <- logZnew
    ## Posterior Samples (optional)
    Samples[[nest]] <- Obj[[worst]]
    ## Kill worst object in favour of copy of different survivor
    copy <- 1+floor(n * uniform()) %% n   # force 1 <= copy <= n
    while( copy == worst && n > 1 ) {   # don't kill if n is only 1
      copy <- 1+floor(n * uniform()) %% n
    }
    logLstar <- Obj[[worst]]$logL;  # new log likelihood constraint
    Obj[[worst]] <- Obj[[copy]]; # overwrite worst object (see note [1])
    ## Evolve copied object within constraint
    Obj[[worst]] <- explore(Obj[[worst]], logLstar);
    ## Shrink interval
    logwidth <- logwidth - 1.0 / n;
  } # _______ NESTED SAMPLING LOOP (might be ok to terminate early)
  ## Exit with evidence Z, information H, and optional posterior Samples
  list(logZ=logZ, logZ.sd=sqrt(H/n), Hnats=H, Hbits=H/log(2),
       samples=Samples)
}