#
# Run optimization example from Liang et al JASA 2005; 100:1311-1327
#
liang.example <- function(tune=0.3, othermethod="BFGS") {
  require(lattice)
  H <- function(x) -(x[1]*sin(20*x[2]) + x[2]*sin(20*x[1]))^2*cosh(sin(10*x[1])*x[1])-(x[1]*cos(10*x[2])-x[2]*sin(10*x[1]))^2*cosh(cos(20*x[2])*x[2])
  H1 <- function(x) if(all(x>-1.1 & x<1.1)) { H(x) } else {0}
# function surface
# Run Generalised Wang-Landau
  res1 <- gwl.H(H, tune=tune, tol=1e-4, cuts=c(-1e99, seq(-8, 0, 0.2), 1e99))
  image(x=gp,y=gp,matrix(z, nr=round(sqrt(length(z)))))
  lines(res1$estimates[,1:2])
  bestopt <- 100
# Run other maximizer with 20 different starting values
  for(i in 1:20) {
    sta<-rnorm(2, sd=0.5)
    optres <- optim(sta, H1, method=othermethod)
    if (optres$value<bestopt) bestopt <- optres$value
    arrows(sta[1], sta[2], optres$par[1], optres$par[2], 
          col="blue")
  }
  text(0,0,paste("GWL minimum=",min(res1$lik),"\n",
                 othermethod, " minimum=",bestopt,"\n",
                 "True minimum=", H(c(1.0445, -1.0084)), sep=""))
}

gwl.H <- function(H, delta1=1.718, rho=0.1, tune=0.1, 
                  tol=1e-7, iter1=1000, cuts) {
  g <- rep(1.0, length(cuts)-1)
  x <- runif(2, -0.1,0.1)
  lik <- H(x)
  it <- 0
  delta <- delta1
  iter <- iter1
  while(delta>tol) {
    cat(it," lik=", lik, " x=", x, " delta=",delta," iter=", iter, "\n")
    res <- onegwl.H(H, lik, x, cuts, g, delta, rho, tune, iter)
    x <- res$estimates[iter,]
    lik <- res$lik[iter]
    image(x=gp,y=gp,matrix(z, nr=round(sqrt(length(z)))))
    lines(res$estimates[,1:2])
    title(main=paste("Iteration ",it, " delta=", delta, sep=""))
    delta <- sqrt(1+delta)-1
    iter <- as.integer(1.1*iter)
    it <- it + 1
  }
  res
}

onegwl.H <- function(H, lik1, x1, cuts, g, delta, rho, tune, iter) {
  modelik <- function(x,t) if(all(x>-1.1 & x<1.1)) { exp(-H(x)/t) }else{ 0 }
  accepted <- 0
  it <- 0
  t <- 0.1
  estimates <- matrix(nc=2, nr=iter)
  lik <- double(iter)
  flik1 <- modelik(x1,t)
  cat1 <- max(which(lik1>cuts))
  while(it < iter) {
    it <- it+1
    x2 <- x1 + rnorm(2,0,tune)
    lik2 <- H(x2)
    flik2 <- modelik(x2,t)
    cat2 <- max(which(lik2>cuts))
    qa <- g[cat1]/g[cat2] * flik2/flik1
    if (qa>runif(1)) {
      accepted <- accepted + 1
      x1 <- x2
      lik1 <- lik2
      flik1 <- flik2
      cat1 <- cat2
    }
    upg <- seq(0,length(g)-cat1)
    g[cat1+upg] <- g[cat1+upg] + delta * rho^upg * g[cat1]
    estimates[it,] <- x1
    lik[it] <- lik1
  }
  list(estimates=estimates, lik=lik, weight=g, accepted=accepted/iter)
}    
#
# Likelihood surface for example
#
library(lattice)
H <- function(x) -(x[1]*sin(20*x[2]) + x[2]*sin(20*x[1]))^2*cosh(sin(10*x[1])*x[1])-(x[1]*cos(10*x[2])-x[2]*sin(10*x[1]))^2*cosh(cos(20*x[2])*x[2])
H1 <- function(x) if(all(x>-1.1 & x<1.1)) { H(x) } else {0}
gp <- seq(-1.1,1.1,0.01)
xy <- expand.grid(x=gp,y=gp)
z <- apply(xy,1,H)
print(contourplot(z ~ xy[,1]* xy[,2], region=TRUE))
liang.example()