The R Project SVN R-packages

## not run usually:
if(!cluster:::doExtras()) q("no")

library(cluster)

## generate 1500 objects, divided into 2 clusters.
suppressWarnings(RNGversion("3.5.0")) # << as long as we don't have R >= 3.6.0
set.seed(264) ## <-- different than ./clara.R's !
x <- rbind(cbind(rnorm(700, 0,8), rnorm(700, 0,8)),
           cbind(rnorm(800,50,8), rnorm(800,10,8)))

isEq <- function(x,y, epsF = 100)
    is.logical(r <- all.equal(x,y, tol = epsF * .Machine$double.eps)) && r


.proctime00 <- proc.time()

## full size sample {should be = pam()}:
n0 <- length(iSml <- c(1:70, 701:720))
summary(clara0 <- clara(x[iSml,], k = 2, sampsize = n0))
          pam0 <- pam  (x[iSml,], k = 2)
stopifnot(identical(clara0$clustering, pam0$clustering)
        , isEq(clara0$objective, unname(pam0$objective[2]))
          )

summary(clara2 <- clara(x, 2))

clInd <- c("objective", "i.med", "medoids", "clusinfo")
clInS <- c(clInd, "sample")
## clara() {as original code} always draws the *same* random samples !!!!
clara(x, 2, samples = 50)[clInd]

## The for() code below *did* hang (in some 64-bit versions of R / clara / ...):
## but this has been fixed *long* ago
xx <- x[sample(nrow(x)),]
if(FALSE) { ## only run manually
    ff <- "/u/maechler/R/MM/Pkg-ex/cluster/clara2-ex.rds"
    if(file.exists(ff))
        xx <- readRDS(ff)
    else
        saveRDS(xx, file=ff)
}
print(clara(xx, 2, samples= 5, trace = 3)[clInd])

print(clara(xx, 2, samples=50, trace = 2)[clInd])
##                             ^^^^^^^^^ a bit less output ...
## for sample 45
## From that output I gather the sample indices nsel[] and all k{ran} are
ii <- c(194, 1411,
        1430,398,570,72,27,62,1301,368,1390,991,296,431,1019,186,558,258,413,
        647,585,1352,1073,873,377,711,1498,865,1436,1335,189,622,760,226,146,
        145,1349,382,1368,934,204,303,856,1489)
## Trying
pp <- pam(x[ii,], k =2) # gives no problem

for(i in 1:19)
    print(clara(x[sample(nrow(x)),], 2, samples = 50)[clInd])

clara(x, 2, samples = 101)[clInd]
clara(x, 2, samples = 149)[clInd]
clara(x, 2, samples = 200)[clInd]
## Note that this last one is practically identical to the slower  pam() one

x[print(sample(length(x), 20))] <- NA
clara(x, 2, samples = 50)[clInd]

###--- Start a version of  example(clara) -------

## xclara : artificial data with 3 clusters of 1000 bivariate objects each.
data(xclara)
(clx3 <- clara(xclara, 3))
## Plot similar to Figure 5 in Struyf et al (1996)
plot(clx3)

## The  rngR = TRUE case is currently in the non-strict tests
## ./clara-ex.R
## ~~~~~~~~~~~~

## trace output checking
RNGversion("4.3"); RNGkind()
set.seed(144)
xc120 <- xclara[sort(sample.int(nrow(xclara), 120)), ]
cx3.1 <- clara(xc120, 3, trace = 1)
cx3.2 <- clara(xc120, 3, trace = 2)
cx3.3 <- clara(xc120, 3, trace = 3)
cx3.4 <- clara(xc120, 3, trace = 4)
if(FALSE) # no change
cx3.5 <- clara(xc120, 3, trace = 5)

###--- End version of example(clara) -------

## Last Line:
cat('Time elapsed: ', proc.time() - .proctime00,'\n')