# file "Rsimann2" -- a simulated annealing algorithm in one dimension

M = 100000  # run length
temp = 100  # initial value, left constant if no annealing
finaltemp = 0.02
tempfactor = (temp/finaltemp)^(-1/M)  # for exponential cooling
#tempdiff = (temp-finaltemp)/M    # for linearcooling
doannealing = T

pi = function(x, thetemp) {
    return( ( 0.3 * dnorm(x, 0, 1) + 0.7 * dnorm(x, 20, 1) )^(1/thetemp) );
}

h = function(y) { return(y) }

xlist = rep(0,M)  # for keeping track of chain values
templist = rep(0,M)  # for keeping track of temperature values
hlist = rep(0,M)  # for keeping track of h function values
numxaccept = numtempaccept = temptries = 0;
sigma = 1  # proposal scaling
X = runif(1,-10,30)  # overdispersed starting distribution

for (i in 1:M) {
    # UPDATE TEMPERATURE
    if (doannealing)
      # CHOOSE ONE:
      temp = tempfactor * temp;   # exponential cooling
      # temp = temp - tempdiff;   # linear cooling
    # PROPOSED X MOVE
    Y = X + sigma * rnorm(1)  # proposal value
    U = runif(1)  # for accept/reject
    A = pi(Y,temp) / pi(X,temp)  # for accept/reject
    if (U < A) {
	X = Y  # accept proposal
        numxaccept = numxaccept + 1;
    }
    xlist[i] = X;
    templist[i] = temp;
    hlist[i] = h(X);
}

cat("Ran simulated for", M, "iterations.\n");
cat("Final estimate of maximum is", h(X), "achieved at X=", X, ".\n");
cat("acceptance rate =", numxaccept/M, "\n");

# tf = function(x) pi(x,1); plot(tf,-10,30)
# tf = function(x) pi(x,6); plot(tf,-10,30)
# tf = function(x) pi(x,10); plot(tf,-10,30)
plot(xlist, type='l')
# plot(templist, type='l')