R/SwendsenWangAlg.R

In ick003/SpTMixture: Mixture models with spatio-temporal dependencies

#' Simulate a Potts model using the Swendsen-Wang algorithm
#'
#' This function allows you to simulate a Potts model using the Swendsen-Wang algorithm.
#' @param beta parameters of the Potts model.
#' @param Coord Matrix of coordinates.
#' @param col.obs Binary matrix of colors.
#' @param Nrun Number of MCMC runs.
#' @param sampling Choice of sampling for the components.
#' @keywords Swendsen-Wang Algorithm
#' @references Wang, J. and Swendsen, R. (1990) Cluster Monte Carlo algorithms. \emph{Physica A: Statistical Mechanics and its Applications} \bold{167} 565--579.
#' @export
#' @examples
#' beta = c(0.5,0.5,0.5)
#' Obs.loc = expand.grid(1:50,1:50)
#' id.pos = sample(1:3, 2500, replace = T)
#' Res20 = SwendsenWangAlg(beta, Coord = as.matrix(Obs.loc), col.obs = id.pos, Nrun = 20)
#' Res50 = SwendsenWangAlg(beta, Coord = as.matrix(Obs.loc), col.obs = id.pos, Nrun = 50)
#' Res500 = SwendsenWangAlg(beta, Coord = as.matrix(Obs.loc), col.obs = id.pos, Nrun = 500)
#' par(mfrow = c(2,2))
#' col = grey.colors(3)
#' image(matrix(Res20$col.old, ncol = 50), col = col)
#' image(matrix(Res20$col.new, ncol = 50), col = col)
#' image(matrix(Res50$col.new, ncol = 50), col = col)
#' image(matrix(Res500$col.new, ncol = 50), col = col)
SwendsenWangAlg = function(beta, Coord = NULL, Neigh_Bonds = NULL, col.obs, Nrun = 10 ,sampling = "else"){

  n.vert = length(col.obs)

  n.colors = length(beta)

  # Get neighboorhood bonds

  if(is.null(Neigh_Bonds)){Neigh_Bonds = getBonds(Coord, NN = 4, th = 2)}

  # Select bonds
  col.temp = col.obs

  z.p = loopSW(Bds = Neigh_Bonds, Cols = col.obs, ncolors = n.colors, Nrun = Nrun, Betas = beta)

  RET = list(col.new = z.p, col.old = col.obs)
  return(RET)
}