\HeaderA{clara}{Clustering Large Applications}{clara}
\keyword{cluster}{clara}
\begin{Description}\relax
Computes a \code{"clara"} object, a list representing a clustering of
the data into \code{k} clusters.
\end{Description}
\begin{Usage}
\begin{verbatim}
clara(x, k, metric = "euclidean", stand = FALSE, samples = 5,
      sampsize = min(n, 40 + 2 * k), trace = 0, medoids.x = TRUE,
      keep.data = medoids.x, rngR = FALSE)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{x}] data matrix or data frame, each row corresponds to an observation,
and each column corresponds to a variable.  All variables must be numeric.
Missing values (NAs) are allowed.
\item[\code{k}] integer, the number of clusters.
It is required that \eqn{0 < k < n}{} where \eqn{n}{} is the number of
observations (i.e., n = \code{nrow(x)}).
\item[\code{metric}] character string specifying the metric to be used for calculating
dissimilarities between observations.
The currently available options are "euclidean" and "manhattan".
Euclidean distances are root sum-of-squares of differences, and
manhattan distances are the sum of absolute differences.

\item[\code{stand}] logical, indicating if the measurements in \code{x} are
standardized before calculating the dissimilarities.  Measurements
are standardized for each variable (column), by subtracting the
variable's mean value and dividing by the variable's mean absolute
deviation.

\item[\code{samples}] integer, number of samples to be drawn from the dataset.
\item[\code{sampsize}] integer, number of observations in each
sample. \code{sampsize} should be higher than the number of clusters
(\code{k}) and at most the number of observations (n = \code{nrow(x)}).
\item[\code{trace}] integer indicating a \emph{trace level} for diagnostic
output during the algorithm.
\item[\code{medoids.x}] logical indicating if the medoids should be
returned, identically to some rows of the input data \code{x}.  If
\code{FALSE}, \code{keep.data} must be false as well, and the medoid
indices, i.e., row numbers of the medoids will still be returned
(\code{i.med} component), and the algorithm saves space by needing
one copy less of \code{x}.
\item[\code{keep.data}] logical indicating if the (\emph{scaled} if
\code{stand} is true) data should be kept in the result.
Setting this to \code{FALSE} saves memory (and hence time), but
disables \code{\LinkA{clusplot}{clusplot}()}ing of the result.  Use
\code{medoids.x = FALSE} to save even more memory.
\item[\code{rngR}] logical indicating if \R's random number generator should
be used instead of the primitive clara()-builtin one.  If true, this
also means that each call to \code{clara()} returns a different result
-- though only slightly different in good situations.
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
\code{clara} is fully described in chapter 3 of Kaufman and Rousseeuw (1990).
Compared to other partitioning methods such as \code{pam}, it can deal with
much larger datasets.  Internally, this is achieved by considering
sub-datasets of fixed size (\code{sampsize}) such that the time and
storage requirements become linear in \eqn{n}{} rather than quadratic.

Each sub-dataset is partitioned into \code{k} clusters using the same
algorithm as in \code{\LinkA{pam}{pam}}.\\
Once \code{k} representative objects have been selected from the
sub-dataset, each observation of the entire dataset is assigned
to the nearest medoid.

The sum of the dissimilarities of the observations to their closest
medoid is used as a measure of the quality of the clustering.  The
sub-dataset for which the sum is minimal, is retained.  A further
analysis is carried out on the final partition.

Each sub-dataset is forced to contain the medoids obtained from the
best sub-dataset until then.  Randomly drawn observations are added to
this set until \code{sampsize} has been reached.
\end{Details}
\begin{Value}
an object of class \code{"clara"} representing the clustering.  See
\code{\LinkA{clara.object}{clara.object}} for details.
\end{Value}
\begin{Note}\relax
By default, the random sampling is implemented with a \emph{very}
simple scheme (with period \eqn{2^{16} = 65536}{}) inside the Fortran
code, independently of \R's random number generation, and as a matter
of fact, deterministically.  Alternatively, we recommend setting
\code{rngR = TRUE} which uses \R's random number generators.  Then,
\code{clara()} results are made reproducible typically by using
\code{\LinkA{set.seed}{set.seed}()} before calling \code{clara}.

The storage requirement of \code{clara} computation (for small
\code{k}) is about
\eqn{O(n \times p) + O(j^2)}{O(n * p) + O(j^2)} where
\eqn{j = \code{sampsize}}{}, and \eqn{(n,p) = \code{dim(x)}}{}.
The CPU computing time (again assuming small \code{k}) is about
\eqn{O(n \times p \times j^2 \times N)}{O(n * p * j^2 * N)}, where
\eqn{N = \code{samples}}{}.

For ``small'' datasets, the function \code{\LinkA{pam}{pam}} can be used
directly.  What can be considered \emph{small}, is really a function
of available computing power, both memory (RAM) and speed.
Originally (1990), ``small'' meant less than 100 observations;
later, the authors said \emph{``small (say with fewer than 200
observations)''}..
\end{Note}
\begin{Author}\relax
Kaufman and Rousseuw, originally.
All arguments from \code{trace} on, and most \R{} documentation and all
tests by Martin Maechler.
\end{Author}
\begin{SeeAlso}\relax
\code{\LinkA{agnes}{agnes}} for background and references;
\code{\LinkA{clara.object}{clara.object}}, \code{\LinkA{pam}{pam}},
\code{\LinkA{partition.object}{partition.object}}, \code{\LinkA{plot.partition}{plot.partition}}.
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}
## generate 500 objects, divided into 2 clusters.
x <- rbind(cbind(rnorm(200,0,8), rnorm(200,0,8)),
           cbind(rnorm(300,50,8), rnorm(300,50,8)))
clarax <- clara(x, 2)
clarax
clarax$clusinfo
plot(clarax)

## `xclara' is an artificial data set with 3 clusters of 1000 bivariate
## objects each.
data(xclara)
(clx3 <- clara(xclara, 3))
## Plot similar to Figure 5 in Struyf et al (1996)
## Not run: plot(clx3, ask = TRUE)


## Try 100 times *different* random samples -- for reliability:
nSim <- 100
nCl <- 3 # = no.classes
set.seed(421)# (reproducibility)
cl <- matrix(NA,nrow(xclara), nSim)
for(i in 1:nSim)
   cl[,i] <- clara(xclara, nCl, medoids.x = FALSE, rngR = TRUE)$cluster
tcl <- apply(cl,1, tabulate, nbins = nCl)
## those that are not always in same cluster (5 out of 3000 for this seed):
(iDoubt <- which(apply(tcl,2, function(n) all(n < nSim))))
if(length(iDoubt)) { # (not for all seeds)
  tabD <- tcl[,iDoubt, drop=FALSE]
  dimnames(tabD) <- list(cluster = paste(1:nCl), obs = format(iDoubt))
  t(tabD) # how many times in which clusters
}

\end{ExampleCode}
\end{Examples}

