Last data update: 2014.03.03

 R: Find disjoint spatial blocks of sites
calcBlockWR Documentation

Find disjoint spatial blocks of sites

Description

Partition the space into disjoint spatial blocks of sites. Call the C++ function calcOrderSite_hpp. For internal use.

Usage

calcBlockW(W, site_order = NULL, dist.center = NULL, dist.max = Inf, 
         verbose = optionsMRIaggr("verbose"))

Arguments

W

the neighbourhood matrix. dgCMatrix. REQUIRED.

site_order

a specific order to go all over the sites. integer vector.

dist.center

the distance between each point and a reference point. numeric vector.

dist.max

the neighbourhood range. numeric vector.

verbose

Should the process be verbose over iterations ? logical.

Details

This function requires to have installed the Matrix and the spam package to work. If no specific order is set, sites are visitating from the first to the last, according to the neighbourhood matrix.

Value

An list containing :

  • [[ls_groups]] : a list containing the index of the sites for each independant group.

  • [[size_groups]] : a vector containing the size of each independant group.

  • [[n_groups]] : an integer giving the number of independant groups.

Examples

#### spatial field
## Not run: 
n <- 100

## End(Not run)

coords <- data.frame(which(matrix(0, nrow = n, ncol = n) == 0,arr.ind = TRUE), 1)
optionsMRIaggr(quantiles.legend = FALSE, axes = FALSE, num.main = FALSE, bg = "white")

#### 1- neighbourhood matrix (king) ####
W_king <- calcW(coords, range = 1.001, row.norm = TRUE)$W

#### find independant groups
Block_king <- calcBlockW(W_king)

## check groups
# diagonal : percent of neighborhing sites whithin group
# extra-diagonal : percent of neighborhing sites between groups
sapply(1:Block_king$n_groups, function(x){
  sapply(1:Block_king$n_groups, function(y){
    sum(spam::rowSums(W_king[Block_king$ls_groups[[x]], Block_king$ls_groups[[y]]] > 0) > 0)
  }) / length(Block_king$ls_groups[[x]])
}
)

## diplay sparse matrix
spam::image(W_king)
spam::image(W_king[unlist(Block_king$ls_groups), unlist(Block_king$ls_groups)])

## display site blocks
col_sites <- unlist(lapply(1:Block_king$n_groups, function(x){
	rep(rainbow(Block_king$n_groups)[x], Block_king$size_groups[x])
}))

multiplot(coords[unlist(Block_king$ls_groups),],
          xlim = c(0,30),ylim = c(0,30),
          col = col_sites, legend = FALSE)


#### 2- neighbourhood matrix (Queen) ####
W_queen <- calcW(coords, range = sqrt(2) + 0.001, row.norm = TRUE)$W

#### find independant groups
Block_queen <- calcBlockW(W_queen)

## check groups
# diagonal : percent of neighborhing sites whithin group
# extra-diagonal : percent of neighborhing sites between groups
sapply(1:Block_queen$n_groups, function(x){
  sapply(1:Block_queen$n_groups, function(y){
    sum(spam::rowSums(W_queen[Block_queen$ls_groups[[x]], Block_queen$ls_groups[[y]]] > 0) > 0)
  }) / length(Block_queen$ls_groups[[x]])
}
)

## diplay sparse matrix
spam::image(W_queen)
spam::image(W_queen[unlist(Block_queen$ls_groups), unlist(Block_queen$ls_groups)])

## display site blocks
col_sites <- unlist(lapply(1:Block_queen$n_groups, function(x){
	rep(rainbow(Block_queen$n_groups)[x], Block_queen$size_groups[x])
}))

multiplot(coords[unlist(Block_queen$ls_groups),],
          xlim = c(0,30), ylim = c(0,30),
          col = col_sites, legend = FALSE)

#### 3- neighbourhood matrix (Regional) ####
W_Regional <- calcW(coords, range = 3, row.norm = TRUE)$W

#### find independant groups
system.time(
  Block_Regional <- calcBlockW(W_Regional)
)

system.time(
Block_Regional_test1 <- calcBlockW(W_Regional, 
     dist.center = sqrt(spam::rowSums(sweep(coords, MARGIN = 2, 
	                                  STATS = apply(coords, 2, median), FUN = "-")^2))
     )
)
system.time(
  Block_Regional_test2 <- calcBlockW(W_Regional, 
     dist.center = sqrt(spam::rowSums(sweep(coords, MARGIN = 2,
                            	      STATS = apply(coords, 2, median), FUN = "-")^2)),
     dist.max = 3
  )
)
# all(unlist(Block_Regional_test1$ls_groups) == unlist(Block_Regional_test2$ls_groups))


## check groups
# diagonal : percent of neighborhing sites whithin group
# extra-diagonal : percent of neighborhing sites between groups
sapply(1:Block_Regional$n_groups,function(x){
   sapply(1:Block_Regional$n_groups,function(y){
    if(length(Block_Regional$ls_groups[[x]]) > 1){
      sum(spam::rowSums(as.matrix(W_Regional[Block_Regional$ls_groups[[x]],
	                                   Block_Regional$ls_groups[[y]]]) > 0) > 0)
    }else{
      sum(W_Regional[Block_Regional$ls_groups[[x]],
	                 Block_Regional$ls_groups[[y]]] > 0) > 0
    }
  }) / length(Block_Regional$ls_groups[[x]])
}
)
# clustering could be improved

## diplay sparse matrix
spam::image(W_Regional)
spam::image(W_Regional[unlist(Block_Regional$ls_groups), unlist(Block_Regional$ls_groups)])

## display site blocks
col_sites <- unlist(lapply(1:Block_Regional$n_groups, function(x){
rep(rainbow(Block_Regional$n_groups)[x], Block_Regional$size_groups[x])
}))

multiplot(coords[unlist(Block_Regional$ls_groups),],
          xlim = c(0,30), ylim = c(0,30),
          col = col_sites, legend = 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/MRIaggr-calcBlockW.Rd_%03d_medium.png", width=480, height=480)
> ### Name: calcBlockW
> ### Title: Find disjoint spatial blocks of sites
> ### Aliases: calcBlockW
> ### Keywords: functions
> 
> ### ** Examples
> 
> #### spatial field
> ## Not run: 
> ##D n <- 100
> ## End(Not run)
> ## Don't show: 
> n <- 10
> ## End(Don't show)
> coords <- data.frame(which(matrix(0, nrow = n, ncol = n) == 0,arr.ind = TRUE), 1)
> optionsMRIaggr(quantiles.legend = FALSE, axes = FALSE, num.main = FALSE, bg = "white")
> 
> #### 1- neighbourhood matrix (king) ####
> W_king <- calcW(coords, range = 1.001, row.norm = TRUE)$W
> 
> #### find independant groups
> Block_king <- calcBlockW(W_king)
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
> 
> ## check groups
> # diagonal : percent of neighborhing sites whithin group
> # extra-diagonal : percent of neighborhing sites between groups
> sapply(1:Block_king$n_groups, function(x){
+   sapply(1:Block_king$n_groups, function(y){
+     sum(spam::rowSums(W_king[Block_king$ls_groups[[x]], Block_king$ls_groups[[y]]] > 0) > 0)
+   }) / length(Block_king$ls_groups[[x]])
+ }
+ )
     [,1] [,2]
[1,]    0    1
[2,]    1    0
> 
> ## diplay sparse matrix
> spam::image(W_king)
> spam::image(W_king[unlist(Block_king$ls_groups), unlist(Block_king$ls_groups)])
> 
> ## display site blocks
> col_sites <- unlist(lapply(1:Block_king$n_groups, function(x){
+ 	rep(rainbow(Block_king$n_groups)[x], Block_king$size_groups[x])
+ }))
> 
> multiplot(coords[unlist(Block_king$ls_groups),],
+           xlim = c(0,30),ylim = c(0,30),
+           col = col_sites, legend = FALSE)
> 
> 
> #### 2- neighbourhood matrix (Queen) ####
> W_queen <- calcW(coords, range = sqrt(2) + 0.001, row.norm = TRUE)$W
> 
> #### find independant groups
> Block_queen <- calcBlockW(W_queen)
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
> 
> ## check groups
> # diagonal : percent of neighborhing sites whithin group
> # extra-diagonal : percent of neighborhing sites between groups
> sapply(1:Block_queen$n_groups, function(x){
+   sapply(1:Block_queen$n_groups, function(y){
+     sum(spam::rowSums(W_queen[Block_queen$ls_groups[[x]], Block_queen$ls_groups[[y]]] > 0) > 0)
+   }) / length(Block_queen$ls_groups[[x]])
+ }
+ )
     [,1] [,2] [,3] [,4]
[1,]    0    1    1    1
[2,]    1    0    1    1
[3,]    1    1    0    1
[4,]    1    1    1    0
> 
> ## diplay sparse matrix
> spam::image(W_queen)
> spam::image(W_queen[unlist(Block_queen$ls_groups), unlist(Block_queen$ls_groups)])
> 
> ## display site blocks
> col_sites <- unlist(lapply(1:Block_queen$n_groups, function(x){
+ 	rep(rainbow(Block_queen$n_groups)[x], Block_queen$size_groups[x])
+ }))
> 
> multiplot(coords[unlist(Block_queen$ls_groups),],
+           xlim = c(0,30), ylim = c(0,30),
+           col = col_sites, legend = FALSE)
> 
> #### 3- neighbourhood matrix (Regional) ####
> W_Regional <- calcW(coords, range = 3, row.norm = TRUE)$W
> 
> #### find independant groups
> system.time(
+   Block_Regional <- calcBlockW(W_Regional)
+ )
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
   user  system elapsed 
  0.000   0.000   0.001 
> 
> system.time(
+ Block_Regional_test1 <- calcBlockW(W_Regional, 
+      dist.center = sqrt(spam::rowSums(sweep(coords, MARGIN = 2, 
+ 	                                  STATS = apply(coords, 2, median), FUN = "-")^2))
+      )
+ )
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
   user  system elapsed 
  0.000   0.000   0.002 
> system.time(
+   Block_Regional_test2 <- calcBlockW(W_Regional, 
+      dist.center = sqrt(spam::rowSums(sweep(coords, MARGIN = 2,
+                             	      STATS = apply(coords, 2, median), FUN = "-")^2)),
+      dist.max = 3
+   )
+ )
0% 10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
**************************************************|
|
   user  system elapsed 
  0.004   0.000   0.001 
> # all(unlist(Block_Regional_test1$ls_groups) == unlist(Block_Regional_test2$ls_groups))
> 
> 
> ## check groups
> # diagonal : percent of neighborhing sites whithin group
> # extra-diagonal : percent of neighborhing sites between groups
> sapply(1:Block_Regional$n_groups,function(x){
+    sapply(1:Block_Regional$n_groups,function(y){
+     if(length(Block_Regional$ls_groups[[x]]) > 1){
+       sum(spam::rowSums(as.matrix(W_Regional[Block_Regional$ls_groups[[x]],
+ 	                                   Block_Regional$ls_groups[[y]]]) > 0) > 0)
+     }else{
+       sum(W_Regional[Block_Regional$ls_groups[[x]],
+ 	                 Block_Regional$ls_groups[[y]]] > 0) > 0
+     }
+   }) / length(Block_Regional$ls_groups[[x]])
+ }
+ )
           [,1]      [,2] [,3]      [,4]  [,5]      [,6]  [,7]  [,8]      [,9]
 [1,] 0.0000000 1.0000000  1.0 1.0000000 1.000 1.0000000 1.000 1.000 1.0000000
 [2,] 1.0000000 0.0000000  1.0 1.0000000 1.000 1.0000000 1.000 1.000 1.0000000
 [3,] 1.0000000 1.0000000  0.0 1.0000000 1.000 1.0000000 1.000 1.000 1.0000000
 [4,] 1.0000000 1.0000000  1.0 0.0000000 1.000 1.0000000 1.000 1.000 1.0000000
 [5,] 1.0000000 1.0000000  1.0 1.0000000 0.000 1.0000000 1.000 1.000 1.0000000
 [6,] 1.0000000 1.0000000  1.0 1.0000000 1.000 0.0000000 1.000 1.000 1.0000000
 [7,] 0.9166667 1.0000000  1.0 1.0000000 1.000 0.8888889 0.000 1.000 1.0000000
 [8,] 1.0000000 0.9166667  1.0 1.0000000 1.000 1.0000000 1.000 0.000 1.0000000
 [9,] 0.8333333 0.9166667  0.9 1.0000000 1.000 1.0000000 1.000 1.000 0.0000000
[10,] 0.7500000 0.7500000  0.9 0.7777778 0.875 0.8888889 1.000 1.000 1.0000000
[11,] 0.6666667 0.6666667  1.0 0.8888889 0.875 1.0000000 0.875 1.000 1.0000000
[12,] 0.6666667 0.5833333  0.9 0.7777778 0.875 0.8888889 0.875 0.875 0.8333333
[13,] 0.6666667 0.5833333  0.5 0.5555556 0.625 0.5555556 0.750 0.625 0.8333333
[14,] 0.2500000 0.2500000  0.2 0.2222222 0.250 0.3333333 0.375 0.250 0.3333333
      [,10] [,11] [,12]     [,13] [,14]
 [1,]   1.0   1.0  1.00 1.0000000     1
 [2,]   1.0   1.0  1.00 1.0000000     1
 [3,]   1.0   1.0  1.00 1.0000000     1
 [4,]   1.0   1.0  1.00 1.0000000     1
 [5,]   1.0   1.0  1.00 1.0000000     1
 [6,]   1.0   1.0  1.00 1.0000000     1
 [7,]   1.0   1.0  1.00 1.0000000     1
 [8,]   1.0   1.0  1.00 1.0000000     1
 [9,]   1.0   1.0  1.00 1.0000000     1
[10,]   0.0   1.0  1.00 1.0000000     1
[11,]   0.8   0.0  1.00 1.0000000     1
[12,]   1.0   1.0  0.00 1.0000000     1
[13,]   0.8   0.8  0.75 0.0000000     1
[14,]   0.4   0.2  0.25 0.3333333     0
> # clustering could be improved
> 
> ## diplay sparse matrix
> spam::image(W_Regional)
> spam::image(W_Regional[unlist(Block_Regional$ls_groups), unlist(Block_Regional$ls_groups)])
> 
> ## display site blocks
> col_sites <- unlist(lapply(1:Block_Regional$n_groups, function(x){
+ rep(rainbow(Block_Regional$n_groups)[x], Block_Regional$size_groups[x])
+ }))
> 
> multiplot(coords[unlist(Block_Regional$ls_groups),],
+           xlim = c(0,30), ylim = c(0,30),
+           col = col_sites, legend = FALSE)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>


Created & Maintained by Osamu Ogasawara (osamu.ogasawara@gmail.com) and IMSBIO