Skip to content

Commit 04c2320

Browse files
snoweyesnoweye
committed
Need more test
1 parent 554812f commit 04c2320

File tree

2 files changed

+146
-59
lines changed

2 files changed

+146
-59
lines changed

‎R/00_pmclust_reduceK.r‎

Lines changed: 65 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -1,94 +1,100 @@
11
### For automatically reducing K methods.
22

3-
pmclust.reduceK <- function(X = NULL, K = 2, MU = NULL,
4-
algorithm = .PMC.CT$algorithm, RndEM.iter = .PMC.CT$RndEM.iter,
5-
CONTROL = .PMC.CT$CONTROL, method.own.X = .PMC.CT$method.own.X,
6-
rank.own.X = .pbd_env$SPMD.CT$rank.source, comm = .pbd_env$SPMD.CT$comm){
7-
if(algorithm[1] == "kmeans"){
3+
### X should be in spmd, gbd, or dmat and set at .pmclustEnv or so, as used
4+
### in pmclust().
5+
pmclust.reduceK <- function(K = 2, algorithm = .PMC.CT$algorithm){
6+
if(any(algorithm[1] %in% c("kmeans", "kmeans.dmat"))){
87
stop("kmeans/pkmeans is not supported in reduceK.")
98
}
109

11-
# Run through original pmclust().
12-
ret <- pmclust(X = X, K = K, MU = MU, algorithm = algorithm,
13-
RndEM.iter = RndEM.iter, CONTROL = CONTROL,
14-
method.own.X = method.own.X, rank.own.X = rank.own.X,
15-
comm = comm)
10+
if(algorithm[1] %in% .PMC.CT$algorithm.gbd){
11+
ret <- pmclust.reduceK.spmd(X = X, K = K, algorithm = algorithm)
12+
} else if(algorithm[1] %in% .PMC.CT$algorithm.dmat){
13+
ret <- pmclust.reduceK.dmat(X = X, K = K, algorithm = algorithm)
14+
} else{
15+
comm.stop("The algorithm is not found.")
16+
}
17+
18+
ret
19+
} # End of pmclust.reduceK().
20+
21+
22+
pmclust.reduce.spmd <- function(K = 2, algorithm = .PMC.CT$algorithm){
23+
# Get an initial start.
24+
PARAM.org <- set.global(K = K)
25+
PARAM.org <- try(initial.em(PARAM.org))
26+
27+
# Ensure the initial is good. Warning: This may take forever to run!
28+
repeat{
29+
if(class(PARAM.org) == "try-error"){
30+
PARAM.org <- set.global(K = K)
31+
PARAM.org <- try(initial.em(PARAM.org))
32+
} else{
33+
break
34+
}
35+
}
36+
37+
# Update steps.
38+
method.step <- switch(algorithm[1],
39+
"em" = em.step,
40+
"aecm" = aecm.step,
41+
"apecm" = apecm.step,
42+
"apecma" = apecma.step,
43+
NULL)
44+
if(comm.all(is.null(method.step))){
45+
comm.stop("Algorithm is not found.")
46+
}
47+
PARAM.new <- try(method.step(PARAM.org))
48+
em.update.class()
49+
N.CLASS <- get.N.CLASS(K)
1650

17-
# Repeat if error occurs.
51+
52+
# Reduce K if error occurs.
1853
repeat{
19-
if(ret$check$convergence == 99 && K > 1){
54+
if((class(PARAM.new) == "try-error" ||
55+
.pmclustEnv$CHECK$convergence == 99) &&
56+
K > 1){
2057
# Drop specific i.k if available or
2158
# drop the smallest class or
2259
# drop the class with the smallest eta among all small classes or
2360
# drop all classes with 0 elements.
24-
PARAM.new <- ret$param
2561
if(.pmclustEnv$CONTROL$stop.at.fail && .pmclustEnv$FAIL.i.k > 0){
2662
i.k <- .pmclustEnv$FAIL.i.k
2763
} else{
28-
i.k <- which(ret$n.class == min(ret$n.class))
64+
i.k <- which(N.CLASS == min(N.CLASS))
2965
}
30-
if(i.k > 1 && min(ret$n.class) > 0){
66+
if(i.k > 1 && min(N.CLASS) > 0){
3167
i.k <- i.k[which.min(PARAM.new$ETA[i.k])]
3268
}
3369
K <- K - length(i.k)
70+
comm.cat("- Reduce: ", K, "\n")
3471

3572
# Initial global storage.
36-
if(algorithm[1] %in% .PMC.CT$algorithm.gbd){
37-
PARAM.org <- set.global(K = K)
38-
} else if(algorithm[1] %in% .PMC.CT$algorithm.dmat){
39-
PARAM.org <- set.global.dmat(K = K)
40-
} else{
41-
comm.stop("The algorithm is not found.")
42-
}
73+
PARAM.org <- set.global(K = K)
4374

4475
# Replacing PARAM.org by previous PARAM.new.
45-
PARAM.org$ETA <- PARAM.new$ETA[-i.k] / sum(PARAM.org$ETA[-i.k])
76+
PARAM.org$ETA <- PARAM.new$ETA[-i.k] / sum(PARAM.new$ETA[-i.k])
4677
PARAM.org$log.ETA <- log(PARAM.org$ETA)
4778
PARAM.org$MU <- matrix(PARAM.new$MU[, -i.k], ncol = K)
4879
PARAM.org$SIGMA <- PARAM.new$SIGMA[-i.k]
4980

50-
# Need one e-step to initial storage.
51-
if(algorithm[1] %in% .PMC.CT$algorithm.gbd){
52-
e.step.spmd(PARAM.org)
53-
} else if(algorithm[1] %in% .PMC.CT$algorithm.dmat){
54-
e.step.dmat(PARAM.org)
55-
} else{
56-
comm.stop("The algorithm is not found.")
57-
}
58-
5981
# Update steps.
60-
method.step <- switch(algorithm[1],
61-
"em" = em.step,
62-
"aecm" = aecm.step,
63-
"apecm" = apecm.step,
64-
"apecma" = apecma.step,
65-
NULL)
66-
PARAM.new <- method.step(PARAM.org)
67-
68-
# Obtain classifications.
82+
e.step.spmd(PARAM.org)
83+
PARAM.new <- try(method.step(PARAM.org))
6984
em.update.class()
70-
71-
# Get class numbers.
72-
if(algorithm[1] %in% .PMC.CT$algorithm.gbd){
73-
N.CLASS <- get.N.CLASS(K)
74-
} else if(algorithm[1] %in% .PMC.CT$algorithm.dmat){
75-
N.CLASS <- get.N.CLASS.dmat(K)
76-
} else{
77-
comm.stop("The algorithm is not found.")
78-
}
79-
80-
81-
# For return.
82-
ret <- list(algorithm = algorithm[1],
83-
param = PARAM.new,
84-
class = .pmclustEnv$CLASS.spmd,
85-
n.class = N.CLASS,
86-
check = .pmclustEnv$CHECK)
85+
N.CLASS <- get.N.CLASS(K)
8786
} else{
8887
break
8988
}
9089
}
9190

91+
# For return.
92+
ret <- list(algorithm = algorithm[1],
93+
param = PARAM.new,
94+
class = .pmclustEnv$CLASS.spmd,
95+
n.class = N.CLASS,
96+
check = .pmclustEnv$CHECK)
97+
9298
ret
93-
} # end of pmclust.reduceK().
99+
} # End of pmclust.reduceK.spmd().
94100

‎R/00_pmclust_reduceK_dmat.r‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
### A dmat version for automatically reducing K methods.
2+
3+
pmclust.reduce.dmat <- function(K = 2, algorithm = .PMC.CT$algorithm){
4+
# Get an initial start.
5+
PARAM.org <- set.global.dmat(K = K)
6+
PARAM.org <- try(initial.em.dmat(PARAM.org))
7+
8+
# Ensure the initial is good. Warning: This may take forever to run!
9+
repeat{
10+
if(class(PARAM.org) == "try-error"){
11+
PARAM.org <- set.global.dmat(K = K)
12+
PARAM.org <- try(initial.em.dmat(PARAM.org))
13+
} else{
14+
break
15+
}
16+
}
17+
18+
# Update steps.
19+
method.step <- switch(algorithm[1],
20+
"em.dmat" = em.step,
21+
# "aecm.dmat" = aecm.step,
22+
# "apecm.dmat" = apecm.step,
23+
# "apecma.dmat" = apecma.step,
24+
NULL)
25+
if(comm.all(is.null(method.step))){
26+
comm.stop("Algorithm is not found.")
27+
}
28+
PARAM.new <- try(method.step(PARAM.org))
29+
em.update.class.dmat()
30+
N.CLASS <- get.N.CLASS.dmat(K)
31+
32+
33+
# Reduce K if error occurs.
34+
repeat{
35+
if((class(PARAM.new) == "try-error" ||
36+
.pmclustEnv$CHECK$convergence == 99) &&
37+
K > 1){
38+
# Drop specific i.k if available or
39+
# drop the smallest class or
40+
# drop the class with the smallest eta among all small classes or
41+
# drop all classes with 0 elements.
42+
if(.pmclustEnv$CONTROL$stop.at.fail && .pmclustEnv$FAIL.i.k > 0){
43+
i.k <- .pmclustEnv$FAIL.i.k
44+
} else{
45+
i.k <- which(N.CLASS == min(N.CLASS))
46+
}
47+
if(i.k > 1 && min(N.CLASS) > 0){
48+
i.k <- i.k[which.min(PARAM.new$ETA[i.k])]
49+
}
50+
K <- K - length(i.k)
51+
comm.cat("- Reduce: ", K, "\n")
52+
53+
# Initial global storage.
54+
PARAM.org <- set.global.dmat(K = K)
55+
56+
# Replacing PARAM.org by previous PARAM.new.
57+
PARAM.org$ETA <- PARAM.new$ETA[-i.k] / sum(PARAM.new$ETA[-i.k])
58+
PARAM.org$log.ETA <- log(PARAM.org$ETA)
59+
PARAM.org$MU <- matrix(PARAM.new$MU[, -i.k], ncol = K)
60+
PARAM.org$SIGMA <- PARAM.new$SIGMA[-i.k]
61+
62+
# Update steps.
63+
e.step.dmat(PARAM.org)
64+
PARAM.new <- try(method.step(PARAM.org))
65+
em.update.class.dmat()
66+
N.CLASS <- get.N.CLASS.dmat(K)
67+
} else{
68+
break
69+
}
70+
}
71+
72+
# For return.
73+
ret <- list(algorithm = algorithm[1],
74+
param = PARAM.new,
75+
class = .pmclustEnv$CLASS.spmd,
76+
n.class = N.CLASS,
77+
check = .pmclustEnv$CHECK)
78+
79+
ret
80+
} # End of pmclust.reduceK.dmat().
81+

0 commit comments

Comments
 (0)