|
| 1 | +### For automatically reducing K methods. |
| 2 | + |
| 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"))){ |
| 7 | + stop("kmeans/pkmeans is not supported in reduceK.") |
| 8 | + } |
| 9 | + |
| 10 | + if(algorithm[1] %in% .PMC.CT$algorithm.gbd){ |
| 11 | + ret <- pmclust.reduceK.spmd(K = K, algorithm = algorithm) |
| 12 | + } else if(algorithm[1] %in% .PMC.CT$algorithm.dmat){ |
| 13 | + ret <- pmclust.reduceK.dmat(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.reduceK.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), silent = TRUE) |
| 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), silent = TRUE) |
| 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), silent = TRUE) |
| 48 | + em.update.class() |
| 49 | + N.CLASS <- get.N.CLASS(K) |
| 50 | + |
| 51 | + |
| 52 | + # Reduce K if error occurs. |
| 53 | + repeat{ |
| 54 | + if((class(PARAM.new) == "try-error" || |
| 55 | + .pmclustEnv$CHECK$convergence == 99) && |
| 56 | + K > 1){ |
| 57 | + # Drop specific i.k if available or |
| 58 | + # drop the smallest class or |
| 59 | + # drop the class with the smallest eta among all small classes or |
| 60 | + # drop all classes with 0 elements. |
| 61 | + if(.pmclustEnv$CONTROL$stop.at.fail && .pmclustEnv$FAIL.i.k > 0){ |
| 62 | + i.k <- .pmclustEnv$FAIL.i.k |
| 63 | + } else{ |
| 64 | + i.k <- which(N.CLASS == min(N.CLASS)) |
| 65 | + } |
| 66 | + if(i.k > 1 && min(N.CLASS) > 0){ |
| 67 | + i.k <- i.k[which.min(PARAM.new$ETA[i.k])] |
| 68 | + } |
| 69 | + K <- K - length(i.k) |
| 70 | + comm.cat("- Reduce: ", K, "\n") |
| 71 | + |
| 72 | + # Initial global storage. |
| 73 | + PARAM.org <- set.global(K = K) |
| 74 | + |
| 75 | + # Replacing PARAM.org by previous PARAM.new. |
| 76 | + PARAM.org$ETA <- PARAM.new$ETA[-i.k] / sum(PARAM.new$ETA[-i.k]) |
| 77 | + PARAM.org$log.ETA <- log(PARAM.org$ETA) |
| 78 | + PARAM.org$MU <- matrix(PARAM.new$MU[, -i.k], ncol = K) |
| 79 | + PARAM.org$SIGMA <- PARAM.new$SIGMA[-i.k] |
| 80 | + |
| 81 | + # Update steps. |
| 82 | + e.step.spmd(PARAM.org) |
| 83 | + PARAM.new <- try(method.step(PARAM.org), silent = TRUE) |
| 84 | + em.update.class() |
| 85 | + N.CLASS <- get.N.CLASS(K) |
| 86 | + } else{ |
| 87 | + break |
| 88 | + } |
| 89 | + } |
| 90 | + |
| 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 | + |
| 98 | + ret |
| 99 | +} # End of pmclust.reduceK.spmd(). |
| 100 | + |
0 commit comments