MRIaggr (version 1.1.5)

simulPotts: Draw a sample from a Potts model

Description

Draw a sample from a Potts model. Call the simulPotts_cpp or the simulPottsFast_cpp C++ function.

Usage

simulPotts(W, G, rho, initialization = NULL, site_order = NULL, iter_max = 2500, cv.criterion = 0, fast = FALSE, verbose = TRUE)

Arguments

W

the neighbourhood matrix. dgCMatrix. Should be normalized by row (i.e. rowSums(W)=1). REQUIRED.

G

the number of groups. postive integer. REQUIRED.

rho

the value of the regularization parameter (also called potts parameter or inverse temperature). REQUIRED.

initialization

the probability membership of each voxel to each group used for initialisation. matrix.

site_order

a specific order to go all over the sites. Can be NULL, TRUE or an integer vector.

iter_max

the number of iterations. numeric.

cv.criterion

maximum difference in group probability membership for convergence ? double between 0 and 1.

fast

should simulPottsFast_cpp be used ? logical. Faster but diseable tracing and convergence checking.

verbose

should the simulation be be traced over iterations ? logical.

Value

A numeric vector containing the regional potential.

Details

If no specific order is set, the order is sampled in a uniformed law which increase the computation time. If site_order is set to true, independant spatial block are identifyied using the calcBlockW function. Each independant spatial block is then iteratively considered.

Examples

# spatial field

## Not run: 
# n <- 30
# iter_max <- 500
# ## End(Not run)
G <- 3
coords <- data.frame(which(matrix(0, nrow = n * G, ncol = n * G) == 0, arr.ind = TRUE), 1)
optionsMRIaggr(quantiles.legend = FALSE, axes = FALSE, num.main = FALSE)

# neighbourhood matrix
resW <- calcW(coords, range = sqrt(2), row.norm = TRUE, calcBlockW = TRUE)
W <- resW$W

# with no initial sample
res1 <- simulPotts(W, G = 3, rho = 3, iter_max = iter_max)
multiplot(coords,
          apply(res1$simulation, 1, which.max),
          breaks=seq(0.5, 3.5, 1))

res2 <- simulPotts(W, G = 4, rho = 3, iter_max = iter_max)
multiplot(coords,
          apply(res2$simulation, 1, which.max),
          breaks = seq(0.5, 4.5, 1))

res3 <- simulPotts(W, G = 4, rho = 6, iter_max = iter_max)
multiplot(coords,
          apply(res3$simulation, 1, which.max),
          breaks = seq(0.5, 4.5, 1))

# with specific initialisation
res3.bis <- simulPotts(W, rho = 6, initialization = res3$simulation, iter_max = iter_max / 2)
multiplot(coords,
          apply(res3.bis$simulation, 1, which.max),
          breaks=seq(0.5, 4.5, 1))

res3.ter <- simulPotts(W, rho = 6, initialization = res3$simulation, iter_max = iter_max)
multiplot(coords,
          apply(res3.ter$simulation, 1, which.max),
          breaks=seq(0.5, 4.5, 1))
		  
#### defining site order save time
site_order <- unlist(resW$blocks$ls_groups) - 1

system.time(
  res <- simulPotts(W, iter_max = 100, G = 3, rho = 6, site_order = NULL, 
                    fast = TRUE, verbose = FALSE)
)

system.time(
  res <- simulPotts(W, iter_max = 100, G = 3, rho = 6, site_order = site_order, 
                    fast = TRUE, verbose = FALSE)
)

system.time(
  res <- simulPotts(W, iter_max = 100, G = 3, rho = 6, site_order = NULL, 
                    fast = FALSE, verbose = FALSE)
)

system.time(
  res <- simulPotts(W, iter_max = 100, G = 3, rho = 6, site_order = site_order, 
                    fast = FALSE, verbose = FALSE)
)