1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
#' Functions for listing all possible partitions or just counting the number of them
#'
#' The function \code{nkpartitions} lists all possible partitions of n objects in to k clusters.
#'
# #' @usage nkpartitions(n, k, exact = TRUE, print = FALSE)
# #' nkpar(n, k)
#'
#' @param n Number of units/objects.
#' @param k Number of clusters/groups.
#' @param exact Search for partitions with exactly \code{k} or at most \code{k} clusters.
#' @param print Print results as they are found.
#'
#' @return The matrix or number of possible partitions.
#'
#' @author Chris Andrews
#'
#' @examples
#' n <- 8 # If larger, the number of partitions increases dramatically,
#' # as does if we increase the number of clusters
#' net <- matrix(NA, ncol = n, nrow = n)
#' clu <- rep(1:2, times = c(3, 5))
#' tclu <- table(clu)
#' net[clu == 1, clu == 1] <- rnorm(n = tclu[1] * tclu[1], mean = 0, sd = 1)
#' net[clu == 1, clu == 2] <- rnorm(n = tclu[1] * tclu[2], mean = 4, sd = 1)
#' net[clu == 2, clu == 1] <- rnorm(n = tclu[2] * tclu[1], mean = 0, sd = 1)
#' net[clu == 2, clu == 2] <- rnorm(n = tclu[2] * tclu[2], mean = 0, sd = 1)
#' # Computation of criterion function with the correct partition
#' nkpar(n = n, k = length(tclu)) # Computing the number of partitions
#' all.par <- nkpartitions(n = n, k = length(tclu)) # Forming the partitions
#' all.par <- lapply(apply(all.par, 1, list), function(x) x[[1]])
#' # to make a list out of the matrix
#' res <- critFunC(M = net, clu = clu, approaches = "val",
#' blocks = c("nul", "com"), preSpecM = 4)
#' plot(res) # We get the original partition
#'
#' @keywords cluster
#'
#' @export
"nkpartitions" <-
function(n, k, exact=TRUE, print=FALSE) {
# n objects
# k subgroups
# exactly k or at most k?
# print results as they are found?
# Author: Chris Andrews <candrews@buffalo.edu>
if (n != floor(n) | n<=0) stop("n must be positive integer")
if (k != floor(k) | k<=0) stop("k must be positive integer")
if (print) {
printnkp <- function(a) {
for (j in seq(max(a))) cat("{", seq(along=a)[a==j], "} ")
cat("\n")
}
}
# How many?
Stirling2nd <- function(n, k) {
sum((-1)^seq(0,k-1) * choose(k, seq(0,k-1)) * (k-seq(0,k-1))^n) /
factorial(k)
}
rows <- Stirling2nd(n,k)
if (!exact & k>1) {
for (i in seq(k-1,1)) {
rows <- rows + Stirling2nd(n,i)
}
}
if (print) cat("rows =",rows,"\n")
# Allocate space
theparts <- matrix(NA, nrow=rows, ncol=n)
# begin counting
howmany <- 0
# all in one group
a <- rep(1,n)
# does this count?
if (!exact | (k==1)) {
# increase count, store, and print
howmany <- howmany + 1
theparts[howmany,] <- a
if (print) printnkp(a)
}
# search for others
repeat {
# start at high end
last <- n
repeat {
# increment it if possible
if ((a[last] <= max(a[1:(last-1)])) & (a[last] < k)) {
a[last] <- a[last]+1
# does this count?
if (!exact | max(a)==k) {
# increase count, store, and print
howmany <- howmany + 1
theparts[howmany,] <- a
if (print) printnkp(a)
}
# start again at high end.
break
}
# otherwise set to 1 and move to a different object
a[last] <- 1
if (last>2) {
last <- last-1
next
}
# report the partitions
return(theparts)
}
}
}
|
