#Geoff Nicholls 8/6/09 (correction to the #1'st example 12/09)
#parallelise a single serial MCMC process. See for example
#A. Sohn, "Parallel N-ary Speculative Computation of Simulated Annealing"
#IEEE Trans. Parallel Distrib. Syst.,6,997-1005 (1995).
#for the sampling application in a stats context see
#http://www.stats.ox.ac.uk/~nicholls/linkfiles/papers/cui09.pdf

##################################################
#1) Playing with parallel R (snow) and parallel MCMC
#to check it all makes sense
##################################################

#a function to waste some time
rOP<-function(n)
  {
	b<-rnorm(1)
	for (a in 1:n) {b<-b+rnorm(1)}
	return(b)
  }

#the snow parallel library does all the work
library(snow) 
#my laptop has two processors
cl <- makeCluster(2,type="SOCK") 
clusterSetupRNG(cl)

#clusterCall(cl, runif, 3) 
 
system.time(for (k in 1:2){rOP(1e6)})
system.time(clusterCall(cl,rOP,1e6))
#on my system the above makes parallel look very fast
#(more than x 2 speedup) maybe the different RNG's

#do the comparison with a scalar machine
#defined as a parallel system with one machine
#so the software environment is the same
system.time(r<-parLapply(cl,list(1e6,1e6),rOP)) #was as.list() 12/09
stopCluster(cl)
cl <- makeCluster(1,type="SOCK") 
clusterSetupRNG(cl)
system.time(r<-parLapply(cl,1e6,rOP))
stopCluster(cl)


##################################################
#2 parallel MCMC comparison 
 
stopCluster(cl)
nproc<-2
cl <- makeCluster(nproc,type="SOCK") 
clusterSetupRNG(cl)

ProposeAccept<-function(x=0,w=1,Waiting=1)
{
	#example N(0,1)
	xp<-x+w*(2*runif(1)-1)
	a<-(runif(1)<exp( (x^2-xp^2)/2 ))
	#to test the elgorithm, pretend A/R step time-consuming
	for (k in 1:Waiting) {z<-sqrt(1.1)}
	return(list(xp=xp,a=a))
}


Waiting<-10000

N<-1000
count<-1
X<-matrix(NA,N,1)
X[1]<-0
w<-2
#Standard normal, parallel MCMC
system.time(
while (count<N) {    #changed for to while 12/09
	x<-X[count]
	r<-clusterCall(cl,ProposeAccept,x,w,Waiting)
	a<-0
	for (p in 1:nproc) 
	{
		count<-count+1
		if (r[[p]]$a) 
			{X[count]<-r[[p]]$xp; break}
		else 
			{X[count]<-x}
	}
}
)
count
X<-X[1:N] #so get exactly N and dont stop on acceptance etc
plot(X)

#standard normal, serial MCMC
#notice that the number of parallel steps actually
#completed was random in the simulation above
X<-matrix(NA,N,1)
X[1]<-0
w<-2
system.time(
for (j in 2:N) {
	x<-X[j-1]
	r<-ProposeAccept(x,w,Waiting)
	if (r$a) 
		{X[j]<-r$xp}
	else 
		{X[j]<-x}
}
)
plot(X)