R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Draw a sample from a Potts model

simulPotts

R Documentation

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.

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.

Value

A numeric vector containing the regional potential.

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)
)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(MRIaggr)
Loading required package: Rcpp
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MRIaggr/sfMM-simulPotts.Rd_%03d_medium.png", width=480, height=480)
> ### Name: simulPotts
> ### Title: Draw a sample from a Potts model
> ### Aliases: simulPotts
> 
> ### ** Examples
> 
> # spatial field
> ## Don't show: 
> n <- 10
> iter_max <- 100
> ## End(Don't show)
> ## Not run: 
> ##D n <- 30
> ##D 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)
NOTE : specifying argument 'site_order' would speed up the execution of the function 
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
> 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)
NOTE : specifying argument 'site_order' would speed up the execution of the function 
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
> 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)
NOTE : specifying argument 'site_order' would speed up the execution of the function 
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
> 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)
NOTE : specifying argument 'site_order' would speed up the execution of the function 
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
> 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)
NOTE : specifying argument 'site_order' would speed up the execution of the function 
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
> 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)
+ )
   user  system elapsed 
  0.036   0.000   0.037 
> 
> system.time(
+   res <- simulPotts(W, iter_max = 100, G = 3, rho = 6, site_order = site_order, 
+                     fast = TRUE, verbose = FALSE)
+ )
   user  system elapsed 
  0.032   0.000   0.031 
> 
> system.time(
+   res <- simulPotts(W, iter_max = 100, G = 3, rho = 6, site_order = NULL, 
+                     fast = FALSE, verbose = FALSE)
+ )
   user  system elapsed 
  0.056   0.000   0.055 
> 
> system.time(
+   res <- simulPotts(W, iter_max = 100, G = 3, rho = 6, site_order = site_order, 
+                     fast = FALSE, verbose = FALSE)
+ )
   user  system elapsed 
  0.032   0.000   0.034 
> 
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>