daisy <-
function(x, metric = "euclidean", stand = F, type = list())
{
	meanabsdev <- function(y)
	{
		mean(abs(y - mean(y, na.rm = T)), na.rm = T)
	}
	levs <- function(y)
	{
		levels(as.factor(y))
	}
#check type of input matrix 
	if(!is.data.frame(x) && !is.numeric(x))
		stop(message = "x is not a dataframe or a numeric matrix.")
	if(!is.null(type$asymm)) if(!all(sapply(lapply(as.data.frame(x[, type$
			asymm]), levs), length) == 2))
			stop(message =
				"asymmetric binary variable has more than 2 levels."
				)	
#transform variables and construct `type' vector
	type2 <- sapply(x, data.class)
	x <- data.matrix(x)
	if (length(n <- names(type2[type$ordratio])))
		x[, n] <- codes(as.ordered(x[, n]))
	if (length(n <- names(type2[type$logratio])))
		x[, n] <- log10(x[, n])
	type2[type$asymm] <- "A"
	type2[type$ordratio] <- "O"
	type2[type2 == "numeric"] <- "I"
	type2[type2 == "ordered"] <- "O"
	type2[type2 == "factor"] <- "N"	
#standardize, if necessary
	if(all(type2 == "I")) {
		if(stand) {
			x <- scale(x, scale = apply(x, 2, meanabsdev))
		}
		jdat <- "2"
		if(metric == "manhattan")
			ndyst <- 2
		else ndyst <- 1
	}
	else {
		colmin <- apply(x, 2, min, na.rm = T)
		colextr <- apply(x, 2, max, na.rm = T) - colmin
		x <- scale(x, center = colmin, scale = colextr)
		jdat <- "1"
		ndyst <- 0
	}
	type2 <- paste(type2, collapse = "")
#put info about NAs in arguments for the Fortran call
	jtmd <- ifelse(is.na(rep(1, nrow(x)) %*% x), -1, 1)
	valmisdat <- min(x, na.rm = T) - 0.5
	x[is.na(x)] <- valmisdat
	valmd <- rep(valmisdat, ncol(x))	
#call Fortran routine
	storage.mode(x) <- "double"
	storage.mode(valmd) <- "double"
	storage.mode(jtmd) <- "integer"
	storage.mode(type2) <- "character"
	res <- .Fortran("daisy",
		as.integer(nrow(x)),
		as.integer(ncol(x)),
		x,
		valmd,
		jtmd,
		as.character(jdat),
		type2,
		as.integer(ndyst),
		dis = double(1 + (nrow(x) * (nrow(x) - 1))/2))	
#adapt Fortran output to S-Plus:
#convert lower matrix, read by rows, to upper matrix, read by rows.
	disv <- res$dis[-1]
	disv[disv == -1] <- NA
	full <- matrix(0, nrow(x), nrow(x))
	full[!lower.tri(full, diag = T)] <- disv
	disv <- t(full)[lower.tri(full)]	
#give warning if some dissimilarities are missimg
	if(is.na(min(disv))) attr(disv, "NA.message") <- 
			"NA-values in the dissimilarity matrix !"	
#construct S-Plus object
	class(disv) <- "dissimilarity"
	attr(disv, "Labels") <- dimnames(x)[[1]]
	attr(disv, "Size") <- nrow(x)
	attr(disv, "Metric") <- ifelse(ndyst == 0, "mixed", metric)
	disv
}

print.dissimilarity <- 
function(x, ...)
{
	cat("Dissimilarities :\n")
	print(as.vector(x), ...)
	cat("\n")
	if(!is.null(attr(x, "na.message")))
		cat("Warning : ", attr(x, "NA.message"), "\n")
	cat("Metric : ", attr(x, "Metric"), "\n")
	cat("Number of objects : ", attr(x, "Size"), "\n")
	invisible(x)
}

pam <- 
function(x, k, diss = F, metric = "euclidean", stand = F)
{
	meanabsdev <- function(y)
	{
		mean(abs(y - mean(y, na.rm = T)), na.rm = T)
	}
	size <- function(d)
	{
		sqrtdiscr <- sqrt(1 + 8 * length(d))
		if(sqrtdiscr != as.integer(sqrtdiscr))
			return(0)
		(1 + sqrtdiscr)/2
	}
	if(diss) {
#check type of input vector
		if(is.na(min(x))) stop(message = 
				"NA-values in the dissimilarity matrix not allowed."
				)
		if(data.class(x) != "dissimilarity") {
			if(!is.numeric(x) || size(x) == 0) stop(message = 
				  "x is not of class dissimilarity and can not be converted to this class."
				  )	
	#convert input vector to class "dissimilarity"
			class(x) <- "dissimilarity"
			attr(x, "Size") <- size(x)
			attr(x, "Metric") <- "unspecified"
		}
#adapt S-Plus dissimilarities to Fortran:
#convert upper matrix, read by rows, to lower matrix, read by rows.
		n <- attr(x, "Size")
		full <- matrix(0, n, n)
		full[lower.tri(full)] <- x
		dv <- t(full)[!lower.tri(full, diag = T)]	
	#prepare arguments for the Fortran call
		dv <- c(0, dv)
		jp <- 1
		valmd <- double(1)
		jtmd <- integer(1)
		ndyst <- 0
		x2 <- double(n)
		jdyss <- 1
	}
	else {
#check type of input matrix
		if((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x, 
			data.class) == "numeric"))) stop(message = 
				"x is not a numeric dataframe or matrix.")
		x <- data.matrix(x)	#standardize, if necessary
		if(stand) {
			x2 <- scale(x, scale = apply(x, 2, meanabsdev))
		}
		else x2 <- x
		if(metric == "manhattan")
			ndyst <- 2
		else ndyst <- 1
		n <- nrow(x2)
		jp <- ncol(x2)
		jtmd <- ifelse(is.na(rep(1, n) %*% x2), -1, 1)
		valmisdat <- min(x2, na.rm = T) - 0.5
		x2[is.na(x2)] <- valmisdat
		valmd <- rep(valmisdat, jp)
		jdyss <- 0
		dv <- double(1 + (n * (n - 1))/2)
	}
	storage.mode(dv) <- "double"
	storage.mode(x2) <- "double"
	storage.mode(valmd) <- "double"
	storage.mode(jtmd) <- "integer"
	res <- .Fortran("pam",
		as.integer(n),
		as.integer(jp),
		as.integer(k),
		x2,
		dis = dv,
		ok = as.integer(jdyss),
		valmd,
		jtmd,
		as.integer(ndyst),
		integer(n),
		integer(n),
		integer(n),
		double(n),
		double(n),
		avsil = double(n),
		double(n),
		ttsil = as.double(0),
		med = integer(k),
		obj = double(2),
		clu = integer(n),
		clusinf = matrix(0, k, 5),
		silinf = matrix(0, n, 4),
		isol = integer(k))
	sildim <- res$silinf[, 4]
	if(diss) {
		disv <- x	#add labels to Fortran output
		if(length(attr(x, "Labels")) != 0) {
			sildim <- attr(x, "Labels")[sildim]
			names(res$clu) <- attr(x, "Labels")
			res$med <- attr(x, "Labels")[res$med]
		}
	}
	else {
#give warning if some dissimilarities are missing.
		if(res$ok == -1) stop(message = 
				"No clustering performed, NA-values in the dissimilarity matrix.\n"
				)	#adapt Fortran output to S-Plus:
#convert lower matrix, read by rows, to upper matrix, read by rows.
		disv <- res$dis[-1]
		disv[disv == -1] <- NA
		full <- matrix(0, nrow(x), nrow(x))
		full[!lower.tri(full, diag = T)] <- disv
		disv <- t(full)[lower.tri(full)]
		class(disv) <- "dissimilarity"
		attr(disv, "Size") <- nrow(x)
		attr(disv, "Metric") <- metric
		attr(x, "Labels") <- dimnames(x)[[1]]	
	#add labels to Fortran output
		res$med <- x[res$med,  ]
		if(length((dimnames(x)[[1]])) != 0) {
			sildim <- dimnames(x)[[1]][sildim]
			names(res$clu) <- dimnames(x)[[1]]
		}
	}
	names(res$obj) <- c("build", "swap")
	res$isol <- factor(res$isol, levels = c(0, 1, 2), labels = c("no", "L", 
		"L*"))
	names(res$isol) <- 1:k
	dimnames(res$clusinf) <- list(NULL, c("size", "max_diss", "av_diss", 
		"diameter", "separation"))
	if((k != 1) && (k != n)) {
		dimnames(res$silinf) <- list(sildim, c("cluster", "neighbor", 
			"sil_width", ""))	#construct S-Plus object
		clustering <- list(medoids = res$med, clustering = res$clu, 
			objective = res$obj, isolation = res$isol, clusinfo = 
			res$clusinf, silinfo = list(widths = res$silinf[, -4], 
			clus.avg.widths = res$avsil[1:k], avg.width = res$ttsil
			), diss = disv)
	}
	else {
		clustering <- list(medoids = res$med, clustering = res$clu, 
			objective = res$obj, isolation = res$isol, clusinfo = 
			res$clusinf, diss = disv)
	}
	class(clustering) <- c("pam", "partition")
	attr(clustering, "Call") <- sys.call()
	clustering
}

print.pam <-
function(x, ...)
{
	cat("Medoids:\n")
	print(x$medoids, ...)
	cat("Clustering vector:\n")
	print(x$clustering, ...)
	cat("Objective function:\n")
	print(x$objective, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

summary.pam <- 
function(x)
{
	object <- x
	class(object) <- "summary.pam"
	object
}

print.summary.pam <- 
function(x, ...)
{
	cat("Medoids:\n")
	print(x$medoids, ...)
	cat("Clustering vector:\n")
	print(x$clustering, ...)
	cat("Objective function:\n")
	print(x$objective, ...)
	cat("\nNumerical information per cluster:\n")
	print(x$clusinfo, ...)
	cat("\nIsolated clusters:\n")
	cat("L-clusters: ")
	print(names(x$isolation[x$isolation == "L"]), quote = F, ...)
	cat("L*-clusters: ")
	print(names(x$isolation[x$isolation == "L*"]), quote = F, ...)
	cat("\nSilhouette plot information:\n")
	print(x$silinfo[[1]], ...)
	cat("Average silhouette width per cluster:\n")
	print(x$silinfo[[2]], ...)
	cat("Average silhouette width of total data set:\n")
	print(x$silinfo[[3]], ...)
	cat("\n")
	print(x$diss, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

plot.partition <- 
function(x, ...)
{
	if(is.null(x$silinfo))
		stop(message = "No silhouette plot available when the number of clusters equals 1 or n (the number of objects)."
			)
	s <- rev(x$silinfo[[1]][, 3])
	space <- c(0, rev(diff(x$silinfo[[1]][, 1])))
	space[space != 0] <- 0.5
	names <- if(length(s) < 40) rev(dimnames(x$silinfo[[1]])[[1]]) else 
			NULL
	barplot(s, space = space, names = names, xlab = "Silhouette width", 
		ylab = "Object", xlim = c(min(0, min(s)), 1), horiz = T, mgp = 
		c(2.5, 1, 0), ...)
	title(main = paste("Silhouette plot of ", deparse(attr(x, "Call"))), 
		sub = paste("Average silhouette width : ", round(x$silinfo$
		avg.width, digits = 2)), adj = 0)
	invisible()
}

clara <- 
function(x, k, metric = "euclidean", stand = F, samples = 5, sampsize = 40 + 2 * 
	k)
{
	meanabsdev <- function(y)
	{
		mean(abs(y - mean(y, na.rm = T)), na.rm = T)
	}
#check type of input matrix and values of input numbers
	if((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x, data.class) ==
		"numeric")))
		stop(message = "x is not a numeric dataframe or matrix.")
	x <- data.matrix(x)
	n <- nrow(x)
	if(sampsize < k) {
		warning <- paste(c("'sampsize' should be at least", k, 
			"(number of clusters)"), collapse = " ")
		stop(message = warning)
	}
	if(n < sampsize) {
		warning <- paste(c("Number of objects is", n, 
			", should be at least", sampsize, "(sampsize)"), 
			collapse = " ")
		stop(message = warning)
	}
#standardize, if necessary
	if(stand) {
		x2 <- scale(x, scale = apply(x, 2, meanabsdev))
	}
	else x2 <- x
	if(metric == "manhattan")
		ndyst <- 2
	else ndyst <- 1
	n <- nrow(x2)
	jp <- ncol(x2)
	jtmd <- ifelse(is.na(rep(1, n) %*% x2), -1, 1)
	mdata <- ifelse(is.na(min(x2)), 1, 0)
	valmisdat <- min(x2, na.rm = T) - 0.5
	x2[is.na(x2)] <- valmisdat
	valmd <- rep(valmisdat, jp)
	jdyss <- 0
	x2 <- as.vector(t(x2))	#call Fortran routine
	storage.mode(x2) <- "double"
	storage.mode(valmd) <- "double"
	storage.mode(jtmd) <- "integer"
	res <- .Fortran("clara",
		as.integer(n),
		as.integer(jp),
		as.integer(k),
		clu = x2,
		as.integer(samples),
		as.integer(sampsize),
		dis = double(1 + (sampsize * (sampsize - 1))/2),
		as.integer(mdata),
		valmd,
		jtmd,
		as.integer(ndyst),
		integer(sampsize),
		integer(sampsize),
		sample = integer(sampsize),
		integer(k),
		med = integer(k),
		double(k),
		double(k),
		double(k),
		avdis = double(k),
		maxdis = double(k),
		ratdis = double(k),
		size = integer(k),
		obj = as.double(0),
		avsil = double(k),
		ttsil = as.double(0),
		silinf = matrix(0, sampsize, 4),
		stop = as.integer(0),
		double(sampsize),
		double(sampsize),
		double(sampsize),
		integer(sampsize),
		integer(sampsize),
		integer(sampsize),
		integer(sampsize),
		integer(sampsize),
		integer(sampsize))	#give a warning when errors occured
	if(res$stop == 1)
		stop(message = "For each sample at least one object was found which could not be assigned to a cluster (because of missing values)."
			)
	if(res$stop == 2)
		stop(message = "Each of the random samples contains objects between which no distance can be computed."
			)
	sildim <- res$silinf[, 4]	#adapt Fortran output to S-Plus:
#convert lower matrix, read by rows, to upper matrix, read by rows.
	disv <- res$dis[-1]
	disv[disv == -1] <- NA
	full <- matrix(0, sampsize, sampsize)
	full[!lower.tri(full, diag = T)] <- disv
	disv <- t(full)[lower.tri(full)]
	class(disv) <- "dissimilarity"
	attr(disv, "Size") <- sampsize
	attr(disv, "Metric") <- metric	#add labels to Fortran output
	res$med <- x[res$med,  ]
	res$clu <- matrix(res$clu, nrow = n, ncol = jp, byrow = T)[, 1]
	if(length(dimnames(x)[[1]]) != 0) {
		sildim <- dimnames(x)[[1]][sildim]
		res$sample <- dimnames(x)[[1]][res$sample]
		names(res$clu) <- dimnames(x)[[1]]
	}
#add dimnames to Fortran output
	clusinf <- cbind(res$size, res$maxdis, res$avdis, res$ratdis)
	dimnames(clusinf) <- list(NULL, c("size", "max_diss", "av_diss", 
		"isolation"))
	if((k != 1) && (k != n)) {
		dimnames(res$silinf) <- list(sildim, c("cluster", "neighbor", 
			"sil_width", ""))	#construct S-Plus object
		clustering <- list(sample = res$sample, medoids = res$med, 
			clustering = res$clu, objective = res$obj, clusinfo = 
			clusinf, silinfo = list(width = res$silinf[, -4], 
			clus.avg.widths = res$avsil[1:k], avg.width = res$ttsil
			), diss = disv)
	}
	else {
		clustering <- list(sample = res$sample, medoids = res$med, 
			clustering = res$clu, objective = res$obj, clusinfo = 
			clusinf, diss = disv)
	}
	class(clustering) <- c("clara", "partition")
	attr(clustering, "Call") <- sys.call()
	clustering
}

print.clara <-
function(x, ...)
{
	cat("Best sample:\n")
	print(x$sample, quote = F, ...)
	cat("Medoids:\n")
	print(x$medoids, ...)
	cat("Clustering vector:\n")
	print(x$clustering, ...)
	cat("Objective function:\n")
	print(x$objective, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

summary.clara <-
function(x)
{
	object <- x
	class(object) <- "summary.clara"
	object
}

print.summary.clara <-
function(x, ...)
{
	cat("Best sample:\n")
	print(x$sample, quote = F, ...)
	cat("Medoids:\n")
	print(x$medoids, ...)
	cat("Clustering vector:\n")
	print(x$clustering, ...)
	cat("Objective function:\n")
	print(x$objective, ...)
	cat("\nNumerical information per cluster:\n")
	print(x$clusinfo, ...)
	cat("\nSilhouette plot information for best sample:\n")
	print(x$silinfo[[1]], ...)
	cat("Average silhouette width per cluster:\n")
	print(x$silinfo[[2]], ...)
	cat("Average silhouette width of best sample:\n")
	print(x$silinfo[[3]], ...)
	cat("\n")
	print(x$diss, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

fanny <- 
function(x, k, diss = F, metric = "euclidean", stand = F)
{
	meanabsdev <- function(y)
	{
		mean(abs(y - mean(y, na.rm = T)), na.rm = T)
	}
	size <- function(d)
	{
		sqrtdiscr <- sqrt(1 + 8 * length(d))
		if(sqrtdiscr != as.integer(sqrtdiscr))
			return(0)
		(1 + sqrtdiscr)/2
	}
	if(diss) {
#check type of input vector
		if(is.na(min(x))) stop(message = 
				"NA-values in the dissimilarity matrix not allowed."
				)
		if(data.class(x) != "dissimilarity") {
			if(!is.numeric(x) || size(x) == 0) stop(message = 
				  "x is not of class dissimilarity and can not be converted to this class."
				  )	
	#convert input vector to class "dissimilarity"
			class(x) <- "dissimilarity"
			attr(x, "Size") <- size(x)
			attr(x, "Metric") <- "unspecified"
		}
#prepare arguments for the Fortran call
		n <- attr(x, "Size")
		dv <- c(x, 0)
		jp <- 1
		valmd <- double(1)
		jtmd <- integer(1)
		ndyst <- 0
		x2 <- double(n)
		jdyss <- 1
	}
	else {
#check type of input matrix 
		if((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x, 
			data.class) == "numeric"))) stop(message = 
				"x is not a numeric dataframe or matrix.")
		x <- data.matrix(x)	#standardize, if necessary
		if(stand) {
			x2 <- scale(x, scale = apply(x, 2, meanabsdev))
		}
		else x2 <- x
		if(metric == "manhattan")
			ndyst <- 2
		else ndyst <- 1
		n <- nrow(x2)
		jp <- ncol(x2)
		jtmd <- ifelse(is.na(rep(1, n) %*% x2), -1, 1)
		valmisdat <- min(x2, na.rm = T) - 0.5
		x2[is.na(x2)] <- valmisdat
		valmd <- rep(valmisdat, jp)
		jdyss <- 0
		dv <- double(1 + (n * (n - 1))/2)
	}
	storage.mode(dv) <- "double"
	storage.mode(x2) <- "double"
	storage.mode(valmd) <- "double"
	storage.mode(jtmd) <- "integer"
	res <- .Fortran("fanny",
		as.integer(n),
		as.integer(jp),
		as.integer(k),
		x2,
		dis = dv,
		ok = as.integer(jdyss),
		valmd,
		jtmd,
		as.integer(ndyst),
		integer(n),
		integer(n),
		integer(n),
		double(n),
		p = matrix(0, n, k),
		matrix(0, n, k),
		avsil = double(k),
		integer(k),
		double(k),
		double(k),
		double(n),
		ttsil = as.double(0),
		eda = as.double(0),
		edb = as.double(0),
		obj = double(2),
		clu = integer(n),
		silinf = matrix(0, n, 4),
		as.double(1e-15))
	sildim <- res$silinf[, 4]
	if(diss) {
		disv <- x	#add labels to Fortran output
		if(length(attr(x, "Labels")) != 0) {
			sildim <- attr(x, "Labels")[sildim]
			dimnames(res$p) <- list(attr(x, "Labels"), NULL)
			names(res$clu) <- attr(x, "Labels")
		}
	}
	else {
#give warning if some dissimilarities are missing.
		if(res$ok == -1) stop(message = 
				"No clustering performed, NA-values in the dissimilarity matrix.\n"
				)
		disv <- res$dis[ - (1 + (n * (n - 1))/2)]
		disv[disv == -1] <- NA
		class(disv) <- "dissimilarity"
		attr(disv, "Size") <- nrow(x)
		attr(disv, "Metric") <- metric
		attr(x, "Labels") <- dimnames(x)[[1]]	
	#add labels to Fortran output
		if(length(dimnames(x)[[1]]) != 0) {
			sildim <- dimnames(x)[[1]][sildim]
			dimnames(res$p) <- list(dimnames(x)[[1]], NULL)
			names(res$clu) <- dimnames(x)[[1]]
		}
	}
	names(res$obj) <- c("iterations", "objective")
	res$coeff <- c(res$eda, res$edb)
	names(res$coeff) <- c("dunn_coeff", "normalized")
	if((k != 1) && (k != n)) {
		dimnames(res$silinf) <- list(sildim, c("cluster", "neighbor", 
			"sil_width", ""))	#construct S-Plus object
		clustering <- list(membership = res$p, coeff = res$coeff, 
			clustering = res$clu, objective = res$obj, silinfo = 
			list(widths = res$silinf[, -4], clus.avg.widths = res$
			avsil[1:k], avg.width = res$ttsil), diss = disv)
	}
	else {
		clustering <- list(membership = res$p, coeff = res$coeff, 
			clustering = res$clu, objective = res$obj, diss = disv)
	}
	class(clustering) <- c("fanny", "partition")
	attr(clustering, "Call") <- sys.call()
	clustering
}

print.fanny <-
function(x, ...)
{
	print(x$objective, ...)
	cat("Membership coefficients:\n")
	print(x$membership, ...)
	cat("Coefficients:\n")
	print(x$coeff, ...)
	cat("Closest hard clustering:\n")
	print(x$clustering, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

summary.fanny <- 
function(x)
{
	object <- x
	class(object) <- "summary.fanny"
	object
}

print.summary.fanny <- 
function(x, ...)
{
	print(x$objective, ...)
	cat("Membership coefficients:\n")
	print(x$membership, ...)
	cat("Coefficients:\n")
	print(x$coeff, ...)
	cat("Closest hard clustering:\n")
	print(x$clustering, ...)
	cat("\nSilhouette plot information:\n")
	print(x$silinfo[[1]], ...)
	cat("Average silhouette width per cluster:\n")
	print(x$silinfo[[2]], ...)
	cat("Average silhouette width of total data set:\n")
	print(x$silinfo[[3]], ...)
	cat("\n")
	print(x$diss, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

agnes <-
function(x, diss = F, metric = "euclidean", stand = F, method = "average")
{
	meanabsdev <- function(y)
	{
		mean(abs(y - mean(y, na.rm = T)), na.rm = T)
	}
	size <- function(d)
	{
		sqrtdiscr <- sqrt(1 + 8 * length(d))
		if(sqrtdiscr != as.integer(sqrtdiscr))
			return(0)
		(1 + sqrtdiscr)/2
	}
	if(diss) {
#check type of input vector
		if(is.na(min(x))) stop(message =
				"NA-values in the dissimilarity matrix not allowed."
				)
		if(data.class(x) != "dissimilarity") {
			if(!is.numeric(x) || size(x) == 0)
				stop(message =
				  "x is not of class dissimilarity and can not be converted to this class."
				  )
			n <- size(x)
		}
		else n <- attr(x, "Size")
#adapt S-Plus dissimilarities to Fortran:
#convert upper matrix, read by rows, to lower matrix, read by rows.
		full <- matrix(0, n, n)
		full[lower.tri(full)] <- x
		dv <- t(full)[!lower.tri(full, diag = T)]	
#prepare arguments for the Fortran call
		dv <- c(0, dv)
		jp <- 1
		valmd <- double(1)
		jtmd <- integer(1)
		ndyst <- 0
		x2 <- double(n)
		jdyss <- 1
	}
	else {
#check type of input matrix
		if((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x,
			data.class) == "numeric"))) stop(message =
				"x is not a numeric dataframe or matrix.")
		x <- data.matrix(x)	
#standardize, if necessary
		if(stand) {
			x2 <- scale(x, scale = apply(x, 2, meanabsdev))
		}
		else x2 <- x
#put info about metric, size and NAs in arguments for the Fortran call
		if(metric == "manhattan")
			ndyst <- 2
		else ndyst <- 1
		n <- nrow(x2)
		jp <- ncol(x2)
		jtmd <- ifelse(is.na(rep(1, n) %*% x2), -1, 1)
		valmisdat <- min(x2, na.rm = T) - 0.5
		x2[is.na(x2)] <- valmisdat
		valmd <- rep(valmisdat, jp)
		jdyss <- 0
		dv <- double(1 + (n * (n - 1))/2)
	}
        meth <- 1
        if(method == "single")
                meth <- 2
        if(method == "complete")
                meth <- 3
        if(method == "ward")
                meth <- 4
        if(method == "weighted")
                meth <- 5
	jalg <- 1	
#call Fortran routine
	storage.mode(dv) <- "double"
	storage.mode(x2) <- "double"
	storage.mode(valmd) <- "double"
	storage.mode(jtmd) <- "integer"
	merge <- matrix(0, n - 1, 2)
	storage.mode(merge) <- "integer"
	res <- .Fortran("twins",
		as.integer(n),
		as.integer(jp),
		x2,
		dv,
		ok = as.integer(jdyss),
		valmd,
		jtmd,
		as.integer(ndyst),
		as.integer(jalg),
		as.integer(meth),
		integer(n),
		ner = integer(n),
		ban = double(n),
		ac = as.double(0),
		merge = merge)
	if(!diss) {
#give warning if some dissimilarities are missing.
		if(res$ok == -1)
		          stop(message =
				"No clustering performed, NA-values in the dissimilarity matrix.\n"
				)	
#add labels to Fortran output
		if(length(dimnames(x)[[1]]) != 0) {
			ordertree <- res$ner
			res$ner <- dimnames(x)[[1]][res$ner]
		}
	}
	else {
#add labels to Fortran output
		if(length(attr(x, "Labels")) != 0) {
			ordertree <- res$ner
			res$ner <- attr(x, "Labels")[res$ner]
		}
	}
#construct S-Plus object
	clustering <- list(order = res$ner, height = res$ban[-1], ac = res$ac,
		merge = res$merge)
	if(exists("ordertree"))
		clustering$ordertree <- ordertree
	class(clustering) <- "agnes"
	attr(clustering, "Call") <- sys.call()
	clustering
}

print.agnes <- 
function(x, ...)
{
	cat("Merge:\n")
	print(x$merge, ...)
	cat("Order of objects:\n")
	print(x$order, quote = F, ...)
	cat("Height:\n")
	print(x$height, ...)
	cat("Agglomerative coefficient:\n")
	print(x$ac, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

summary.agnes <- 
function(x)
{
	object <- x
	class(object) <- "summary.agnes"
	object
}

print.summary.agnes <- 
function(x, ...)
{
	print.agnes(x, ...)
	invisible()
}

plot.agnes <- 
function(x, ...)
{
	w <- rev(x$height)
	m <- max(x$height)
	w <- rbind(w, m - w)
	barplot(w, xlab = "Height", horiz = T, inside = F, space = 0, axes = F, 
		col = c(0, 2), mgp = c(2.5, 1, 0), ...)
	# mtext("Object", side = 4)
	title(main = paste("Banner of ", deparse(attr(x, "Call"))), sub = paste(
		"Agglomerative Coefficient = ", round(x$ac, digits = 2)), adj
		 = 0)
	flrm <- floor(m)
	axis(1, at = c(0:flrm, m), labels = c(0:flrm, round(m, digits = 2)), 
		...)
	axis(4, at = 0:(length(x$order) - 1), labels = rev(x$order), pos = m, 
		mgp = c(3, 1.25, 0), ...)
	invisible()
}

pltree.agnes <- 
function(x, ...)
{
	if (!exists("plclust", mode = "function"))
		stop("Requires `plclust' which does not exist yet")
	call <- attr(x, "Call")
	attr(x, "Call") <- NULL
	x$height <- sort(x$height)
	labels <- F
	if(!is.null(x$ordertree)) {
		names(x$order) <- names(x$ordertree) <- 1:length(x$ordertree)
		labels <- x$order[names(sort(x$ordertree))]
		x$order <- x$ordertree
	}
	plclust(as.list(x), labels = labels, plot = T, ylab = "Height", ...)
	title(main = paste("Clustering tree of ", deparse(call)), adj = 0)
	invisible()
}

pltree <- 
function(x, ...)
{
	UseMethod("pltree")
}

diana <-
function(x, diss = F, metric = "euclidean", stand = F)
{
	meanabsdev <- function(y)
	{
		mean(abs(y - mean(y, na.rm = T)), na.rm = T)
	}
	size <- function(d)
	{
		sqrtdiscr <- sqrt(1 + 8 * length(d))
		if(sqrtdiscr != as.integer(sqrtdiscr))
			return(0)
		(1 + sqrtdiscr)/2
	}
	if(diss) {
#check type of input vector
		if(is.na(min(x))) stop(message =
				"NA-values in the dissimilarity matrix not allowed."
				)
		if(data.class(x) != "dissimilarity") {
			if(!is.numeric(x) || size(x) == 0)
				stop(message =
				  "x is not of class dissimilarity and can not be converted to this class."
				  )
			n <- size(x)
		}
		else n <- attr(x, "Size")
#adapt S-Plus dissimilarities to Fortran:
#convert upper matrix, read by rows, to lower matrix, read by rows.
		full <- matrix(0, n, n)
		full[lower.tri(full)] <- x
		dv <- t(full)[!lower.tri(full, diag = T)]	
#prepare arguments for the Fortran call
		dv <- c(0, dv)
		jp <- 1
		valmd <- double(1)
		jtmd <- integer(1)
		ndyst <- 0
		x2 <- double(n)
		jdyss <- 1
	}
	else {
#check type of input matrix 
		if((!is.data.frame(x) && !is.numeric(x)) || (!all(sapply(x, 
			data.class) == "numeric"))) stop(message =
				"x is not a numeric dataframe or matrix.")
		x <- data.matrix(x)	
#standardize, if necessary
		if(stand) {
			x2 <- scale(x, scale = apply(x, 2, meanabsdev))
		}
		else x2 <- x
#put info about metric, size and NAs in arguments for the Fortran call
		if(metric == "manhattan")
			ndyst <- 2
		else ndyst <- 1
		n <- nrow(x2)
		jp <- ncol(x2)
		jtmd <- ifelse(is.na(rep(1, n) %*% x2), -1, 1)
		valmisdat <- min(x2, na.rm = T) - 0.5
		x2[is.na(x2)] <- valmisdat
		valmd <- rep(valmisdat, jp)
		jdyss <- 0
		dv <- double(1 + (n * (n - 1))/2)
	}
	jalg <- 2	
#call Fortran routine
	storage.mode(dv) <- "double"
	storage.mode(x2) <- "double"
	storage.mode(valmd) <- "double"
	storage.mode(jtmd) <- "integer"
	merge <- matrix(0, n - 1, 2)
	storage.mode(merge) <- "integer"
	res <- .Fortran("twins",
		as.integer(n),
		as.integer(jp),
		x2,
		dv,
		ok = as.integer(jdyss),
		valmd,
		jtmd,
		as.integer(ndyst),
		as.integer(jalg),
		as.integer(0),
		integer(n),
		ner = integer(n),
		ban = double(n),
		dc = as.double(0),
		merge = merge)
	if(!diss) {
#give warning if some dissimilarities are missing.
		if(res$ok == -1) stop(message =
				"No clustering performed, NA-values in the dissimilarity matrix.\n"
				)	
#add labels to Fortran output
		if(length(dimnames(x)[[1]]) != 0) {
			ordertree <- res$ner
			res$ner <- dimnames(x)[[1]][res$ner]
		}
	}
	else {
		if(length(attr(x, "Labels")) != 0) {
			ordertree <- res$ner
			res$ner <- attr(x, "Labels")[res$ner]
		}
	}
#construct S-Plus object
	clustering <- list(order = res$ner, height = res$ban[-1], dc = res$dc,
		merge = res$merge)
	if(exists("ordertree"))
		clustering$ordertree <- ordertree
	class(clustering) <- "diana"
	attr(clustering, "Call") <- sys.call()
	clustering
}

print.diana <- 
function(x, ...)
{
	cat("Order of objects:\n")
	print(x$order, quote = F, ...)
	cat("Height:\n")
	print(x$height, ...)
	cat("Divisive coefficient:\n")
	print(x$dc, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

summary.diana <- 
function(x)
{
	object <- x
	class(object) <- "summary.diana"
	object
}

print.summary.diana <- 
function(x, ...)
{
	print.diana(x, ...)
	invisible()
}

plot.diana <- 
function(x, ...)
{
	w <- rev(x$height)
	m <- max(x$height)
	w <- rbind(m - w, w)
	barplot(w, xlab = "Height", ylab = "Object", horiz = T, inside = F, 
		space = 0, axes = F, col = c(2, 0), mgp = c(2.5, 1, 0), ...)
	title(main = paste("Banner of ", deparse(attr(x, "Call"))), sub = paste(
		"Divisive Coefficient = ", round(x$dc, digits = 2)), adj = 0)
	flrm <- floor(m)
	axis(1, at = c(0, (0:flrm) + (m - flrm)), labels = c(round(m, digits = 
		2), rev(0:flrm)), ...)
	axis(2, at = 0:(length(x$order) - 1), labels = rev(x$order), pos = 0, 
		mgp = c(3, 1.5, 0), ...)
	invisible()
}

pltree.diana <- 
function(x, ...)
{
	if (!exists("plclust", mode = "function"))
		stop("Requires `plclust' which does not exist yet")
	call <- attr(x, "Call")
	attr(x, "Call") <- NULL
	x$height <- sort(x$height)
	labels <- F
	if(!is.null(x$ordertree)) {
		names(x$order) <- names(x$ordertree) <- 1:length(x$ordertree)
		labels <- x$order[names(sort(x$ordertree))]
		x$order <- x$ordertree
	}
	plclust(as.list(x), labels = labels, plot = T, ylab = "Height", ...)
	title(main = paste("Clustering tree of ", deparse(call)), adj = 0)
	invisible()
}

mona <-
function(x)
{
	levs <- function(y)
	{
		levels(as.factor(y))
	}
#check type of input vector
	if(!(is.matrix(x) || is.data.frame(x)))
		stop(message = "x is not a matrix or a data frame.")
	if(!all(sapply(lapply(as.data.frame(x), levs), length) == 2))
		stop(message = "All variables must be binary (factor with 2 levels)."
			)
	n <- nrow(x)
	jp <- ncol(x)	
#change levels of input matrix
	x2 <- as.character(apply(as.matrix(x), 2, factor))
	x2[x2 == "1"] <- "0"
	x2[x2 == "2"] <- "1"
	x2[x2 == "NA"] <- "2"
	x2 <- paste(x2, collapse = "")	
#call Fortran routine
	storage.mode(x2) <- "character"
	res <- .Fortran("mona",
		as.integer(n),
		as.integer(jp),
		x2 = x2,
		error = as.integer(0),
		nban = integer(n),
		ner = integer(n),
		integer(n),
		lava = integer(n),
		integer(jp))	
#give a warning when errors occured
	if(res$error == 1)
		stop(message = "No clustering performed, an object was found with all values missing."
			)
	if(res$error == 2)
		stop(message = "No clustering performed, a variable was found with at least 50% missing values."
			)
	if(res$error == 3)
		stop(message = "No clustering performed, a variable was found with all non missing values identical."
			)
	if(res$error == 4)
		stop(message = "No clustering performed, all variables have at least one missing value."
			)
	res$x2 <- matrix(as.numeric(substring(res$x2, 1:nchar(res$x2), 1:nchar(
		res$x2))), n, jp)
	dimnames(res$x2) <- dimnames(x)	
#add labels to Fortran output
	if(length(dimnames(x)[[1]]) != 0)
		res$ner <- dimnames(x)[[1]][res$ner]
	if(length(dimnames(x)[[2]]) != 0) {
		lava <- as.character(res$lava)
		lava[lava != "0"] <- dimnames(x)[[2]][res$lava]
		lava[lava == "0"] <- "NULL"
		res$lava <- lava
	}
#construct S-Plus object
	clustering <- list(data = res$x2, order = res$ner, variable = res$lava[
		-1
		], step = res$nban[-1])
	class(clustering) <- "mona"
	attr(clustering, "Call") <- sys.call()
	clustering
}

print.mona <-
function(x, ...)
{
	cat("Revised data:\n")
	print(x$data, quote = F, ...)
	cat("Ordering of objects:\n")
	print(x$order, quote = F, ...)
	cat("Variable used:\n")
	print(x$variable, quote = F, ...)
	cat("Separation step:\n")
	print(x$step, ...)
	cat("\nAvailable arguments:\n")
	print(names(x), ...)
	invisible(x)
}

summary.mona <- 
function(x)
{
	object <- x
	class(object) <- "summary.mona"
	object
}

print.summary.mona <- 
function(x, ...)
{
	print.mona(x, ...)
	invisible()
}

plot.mona <- 
function(x, ...)
{
	w <- rev(x$step)
	w[w == 0] <- max(w) + 1
	m <- max(w)
	barplot(rbind(w, m - w), xlab = "Separation step", ylab = "Object", 
		horiz = T, inside = F, space = 0, axes = F, col = c(2, 0), mgp
		 = c(2.5, 1, 0), ...)
	title(main = paste("Banner of ", deparse(attr(x, "Call"))), adj = 0)
	axis(1, at = 0:m, labels = 0:m, ...)
	axis(2, at = 0:(length(x$order) - 1), labels = rev(x$order), pos = 0, 
		mgp = c(3, 1.5, 0), ...)
	names <- rev(x$variable)
	names[names == "NULL"] <- ""
	text(w, 0:(length(x$order) - 2) + 0.5, labels = paste(" ", names), adj
		 = 0, col = 2, ...)
	invisible()
}

library.dynam("clus", package = "clus")
