.packageName <- "hopach"
bootplot<-function(bootobj,hopachobj,ord="bootp",main=NULL,labels=NULL,showclusters=TRUE,...){
	p<-nrow(bootobj)
	k<-ncol(bootobj)
	shownames<-(p<50)
	ordering<-hopachobj$clustering$ord
	if(ord=="bootp"){
		start<-1
		stop<-hopachobj$clust$sizes[1]
		set<-ordering[start:stop]
		ordering[start:stop]<-set[rev(order(bootobj[set,1]))]
		for(i in 2:hopachobj$clust$k){
			start<-stop+1
			stop<-cumsum(hopachobj$clust$sizes)[i]
			set<-ordering[start:stop]
			ordering[start:stop]<-set[rev(order(bootobj[set,i]))]	
		}
	}
	if(ord=="final")
		ordering<-hopachobj$final$ord
	if(ord=="none"){
		ordering<-1:p
		showclusters=FALSE
	}
	bootobj<-bootobj[ordering,]
	colors<-rainbow(k)
	colors<-c(colors[seq(1,k,by=2)],colors[seq(2,k,by=2)])
	main<-paste(main,"Barplot of Bootstrap Reappearance Proportions",sep="\n")
	if(is.null(labels))
		labels<-dimnames(bootobj)[[1]]
	par(oma=c(0,0,0,2))
	barplot(t(bootobj),ylim=c(1,p),border=FALSE,space=0,horiz=TRUE,names.arg=labels[ordering],las=1,main=main,cex.names=0.75,legend.text=FALSE,col=colors,axisnames=shownames,xlab="Proportion",...)
	if(showclusters){
		abline(h=cumsum(hopachobj$clust$sizes))
		mtext(colnames(bootobj),outer=TRUE,side=4,at=cumsum(hopachobj$clust$sizes)/p*0.7+0.16,line=-2,col=colors,las=1,cex=0.6)
	}
}
bootmedoids<-function(data,medoids,d="cosangle",I=1000){
	if(inherits(data,"exprSet")) 
		data<-exprs(data)
	data<-as.matrix(data)
	p<-length(data[,1])
	n<-length(data[1,])
	k<-length(medoids)
	blabs<-matrix(0,nrow=p,ncol=I) #holds the bootstrap cluster labels
	bdist<-matrix(0,nrow=p,ncol=k) #holds distance to medoids (recycled)
	props<-matrix(0,nrow=p,ncol=k) #holds the cluster probabilities
	for(i in 1:I){
		samp<-sample(1:n,replace=TRUE)
		for(j in 1:k){
			bdist[,j]<-distancevector(data[,samp],data[medoids[j],samp],d)
		}
       	blabs[,i]<-apply(bdist,1,order)[1,]
	}
	for(i in 1:k)
		props[,i]<-apply(blabs==i,1,mean,na.rm=TRUE)
	props[medoids,]<-diag(k)
	dimnames(props)<-list(dimnames(data)[[1]],paste("Cluster",0:(k-1),sep=""))
	return(props)
}

boothopach<-function(data,hopachobj,I=1000,hopachlabels=FALSE){
	if(inherits(data,"exprSet")) 
		data<-exprs(data)
	data<-as.matrix(data)
	p<-length(data[,1])
	n<-length(data[1,])
	medoids<-hopachobj$clust$medoids
	d<-hopachobj$metric
	k<-length(medoids)
	blabs<-matrix(0,nrow=p,ncol=I) #holds the bootstrap cluster labels
	bdist<-matrix(0,nrow=p,ncol=k) #holds distance to medoids (recycled)
	props<-matrix(0,nrow=p,ncol=k) #holds the cluster probabilities
	for(i in 1:I){
		samp<-sample(1:n,replace=TRUE)
		for(j in 1:k){
			bdist[,j]<-distancevector(data[,samp],data[medoids[j],samp],d)
		}
       	blabs[,i]<-apply(bdist,1,order)[1,]
	}
	for(i in 1:k)
		props[,i]<-apply(blabs==i,1,mean,na.rm=TRUE)
	props[medoids,]<-diag(k)
	if(hopachlabels)
		labs<-as.character(hopachobj$clust$lab[medoids])
	else
		labs<-paste("Cluster",0:(k-1),sep="")
	dimnames(props)<-list(dimnames(data)[[1]],labs)
	return(props)
}
#Distance matrix related functions

#a. compute distances between all rows of a matrix#

#cosine-angle
disscosangle<-function(X,na.rm=TRUE){
	if(!is.matrix(X))
		stop("arg to disscosangle() must be a matrix")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(na.rm){
		X[is.na(X)]<-0
		N<-apply(X!=0,1,sum)
		N<-sqrt(N%*%t(N))/((X!=0)%*%t(X!=0))
	}
	else
		N<-1
	out<-rowSums(X^2)
	out<-matrix(1,nrow=p,ncol=p)-N*(X%*%t(X))/sqrt(out%*%t(out))
	diag(out)<-0
	suppressWarnings(out<-sqrt(out))
	out[out=="NaN"]<-0
	return(out)
}

#absolute cosine-angle
dissabscosangle<-function(X,na.rm=TRUE){
	if(!is.matrix(X))
		stop("arg to dissabscosangle() must be a matrix")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(na.rm){
		X[is.na(X)]<-0
		N<-apply(X!=0,1,sum)
		N<-sqrt(N%*%t(N))/((X!=0)%*%t(X!=0))
	}
	else
		N<-1
	out<-rowSums(X^2)
	out<-matrix(1,nrow=p,ncol=p)-abs(N*(X%*%t(X))/sqrt(out%*%t(out)))
	diag(out)<-0
	suppressWarnings(out<-sqrt(out))
	out[out=="NaN"]<-0
	return(out)
}

#euclidean
#note: disseuclid(X)=daisy(X)/sqrt(dim(X)[2])
disseuclid<-function(X,na.rm=TRUE){
	if(!is.matrix(X))
		stop("arg to disseuclid() must be a matrix")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(na.rm){
		X[is.na(X)]<-0
		N<-apply(X!=0,1,sum)
		N2<-(X!=0)%*%t(X!=0)
		out<-matrix(rep(rowSums(X^2)/N,p),ncol=p)+t(matrix(rep(rowSums(X^2)/N,p),ncol=p))-2*X%*%t(X)/N2	
	}
        else
		out<-matrix(rep(rowMeans(X^2),p),ncol=p)+t(matrix(rep(rowMeans(X^2),p),ncol=p))-2*X%*%t(X)/n
	diag(out)<-0
	suppressWarnings(out<-sqrt(out))
	out[out=="NaN"]<-0
	return(out)
}

#absolute euclidean
dissabseuclid<-function(X,na.rm=TRUE){    
	if(!is.matrix(X))
		stop("arg to dissabseuclid() must be a matrix")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(na.rm){
		X[is.na(X)]<-0
		N<-apply(X!=0,1,sum)
		N2<-(X!=0)%*%t(X!=0)
		out1<-matrix(rep(rowSums(X^2)/N,p),ncol=p)+t(matrix(rep(rowSums(X^2)/N,p),ncol=p))-2*X%*%t(X)/N2	
		out2<-matrix(rep(rowSums(X^2)/N,p),ncol=p)+t(matrix(rep(rowSums(X^2)/N,p),ncol=p))+2*X%*%t(X)/N2
	}
	else{
	        out1<-matrix(rep(rowMeans(X^2),p),ncol=p)+t(matrix(rep(rowMeans(X^2),p),ncol=p))-2*X%*%t(X)/n
	        out2<-matrix(rep(rowMeans(X^2),p),ncol=p)+t(matrix(rep(rowMeans(X^2),p),ncol=p))+2*X%*%t(X)/n
	}
        out1<-pmin(out1,out2)
        diag(out1)<-0
        suppressWarnings(out1<-sqrt(out1))
 	out1[out1=="NaN"]<-0
	return(out1)
}

#correlation
disscor<-function(X,na.rm=TRUE){
	if(!is.matrix(X))
		stop("arg to disscor() must be a matrix")
	p<-dim(X)[1]
	if(na.rm)
		na<-"pairwise.complete.obs"
	else
		na<-"all.obs"
	out<-matrix(rep(1,p*p),nrow=p,ncol=p)-cor(t(X),use=na)
	diag(out)<-0
	suppressWarnings(out<-sqrt(out))
	out[out=="NaN"]<-0
	return(out)
}

#absolute correlation
dissabscor<-function(X,na.rm=TRUE){
	if(!is.matrix(X))
		stop("arg to dissabscor() must be a matrix")
	p<-dim(X)[1]
	if(na.rm)
		na<-"pairwise.complete.obs"
	else
		na<-"all.obs"
	out<-matrix(rep(1,p*p),nrow=p,ncol=p)-abs(cor(t(X),use=na))
	diag(out)<-0
	suppressWarnings(out<-sqrt(out))
	out[out=="NaN"]<-0
	return(out)
}

#b. compute distances between rows of a matrix and a vector#

#cosine-angle
vdisscosangle<-function(X,y,na.rm=TRUE){
	if(!is.matrix(X))
		stop("First arg to vdisscosangle() must be a matrix")
	if(!is.vector(y))
		stop("Second arg to vdisscosangle() must be a vector")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(length(y)!=n)
		stop("Matrix and vector dimensions do not agree in vdisscosangle()")
	if(na.rm){
		X[is.na(X)]<-0
		y[is.na(y)]<-0
		N<-apply(X!=0,1,sum)
		N<-sqrt(N*sum(y!=0))/((X!=0)%*%(y!=0))
	}
	else
		N<-1
	suppressWarnings(out<-sqrt(as.vector(rep(1,p)-N*(X%*%y)/sqrt(rowSums(X^2)*sum(y^2)))))
	out[out=="NaN"]<-0
	return(out)
}

#absolute cosine-angle
vdissabscosangle<-function(X,y,na.rm=TRUE){
	if(!is.matrix(X))
		stop("First arg to vdissabscosangle() must be a matrix")
	if(!is.vector(y))
		stop("Second arg to vdissabscosangle() must be a vector")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(length(y)!=n)
		stop("Matrix and vector dimensions do not agree in vdissabscosangle()")
	if(na.rm){
		X[is.na(X)]<-0
		y[is.na(y)]<-0
		N<-apply(X!=0,1,sum)
		N<-sqrt(N*sum(y!=0))/((X!=0)%*%(y!=0))
	}
	else
		N<-1
	suppressWarnings(out<-sqrt(as.vector(rep(1,p)-abs(N*(X%*%y)/sqrt(rowSums(X^2)*sum(y^2))))))	
	out[out=="NaN"]<-0
	return(out)
}

#euclidean
vdisseuclid<-function(X,y,na.rm=TRUE){
	if(!is.matrix(X))
		stop("First arg to vdisseuclid() must be a matrix")
	if(!is.vector(y))
		stop("Second arg to vdisseuclid() must be a vector")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(length(y)!=n)
		stop("Matrix and vector dimensions do not agree in vdisseuclid()")
	if(na.rm){
		X[is.na(X)]<-0
		y[is.na(y)]<-0
		NX<-apply(X!=0,1,sum)
		Ny<-sum(y!=0)
		N2<-(X!=0)%*%(y!=0)
		suppressWarnings(out<-sqrt(as.vector(rowSums(X^2)/NX+sum(y^2)/Ny-2*X%*%y/N2)))	
	}
	else
		suppressWarnings(out<-sqrt(as.vector(rowMeans(X^2)+mean(y^2)-2*X%*%y/n)))
	out[out=="NaN"]<-0
	return(out)
}

#absolute euclidean
vdissabseuclid<-function(X,y,na.rm=TRUE){
	if(!is.matrix(X))
		stop("First arg to vdissabseuclid() must be a matrix")
	if(!is.vector(y))
		stop("Second arg to vdissabseuclid() must be a vector")
	dX<-dim(X)
	p<-dX[1]
	n<-dX[2]
	if(length(y)!=n)
		stop("Matrix and vector dimensions do not agree in vdissabseuclid()")
	if(na.rm){
		X[is.na(X)]<-0
		y[is.na(y)]<-0
		NX<-apply(X!=0,1,sum)
		Ny<-sum(y!=0)
		N2<-(X!=0)%*%(y!=0)
		out1<-rowSums(X^2)/NX+sum(y^2)/Ny-2*X%*%y/N2
		out2<-rowSums(X^2)/NX+sum(y^2)/Ny+2*X%*%y/N2
	}
	else{
	        out1<-rowMeans(X^2)+mean(y^2)-2*X%*%y/n
       	 	out2<-rowMeans(X^2)+mean(y^2)+2*X%*%y/n
	}
        suppressWarnings(out1<-sqrt(as.vector(pmin(out1,out2))))
	out1[out1=="NaN"]<-0
	return(out1)
}

#correlation
vdisscor<-function(X,y,na.rm=TRUE){
	if(!is.matrix(X))
		stop("First arg to vdisscor() must be a matrix")
	if(!is.vector(y))
		stop("Second arg to vdisscor() must be a vector")
	p<-dim(X)[1]
	if(length(y)!=length(X[1,]))
		stop("Matrix and vector dimensions do not agree in vdisscor()")
	if(na.rm)
		na<-"pairwise.complete.obs"
	else
		na<-"all.obs"
	suppressWarnings(out<-sqrt(as.vector(rep(1,p)-cor(t(X),y,use=na))))
	out[out=="NaN"]<-0
	return(out)
}

#absolute correlation
vdissabscor<-function(X,y,na.rm=TRUE){
	if(!is.matrix(X))
		stop("First arg to vdissabscor() must be a matrix")
	if(!is.vector(y))
		stop("Second arg to vdissabscor() must be a vector")
	p<-dim(X)[1]
	if(length(y)!=length(X[1,]))
		stop("Matrix and vector dimensions do not agree in vdissabscor()")
	if(na.rm)
		na<-"pairwise.complete.obs"
	else
		na<-"all.obs"
	suppressWarnings(out<-sqrt(as.vector(rep(1,p)-abs(cor(t(X),y,use=na)))))
	out[out=="NaN"]<-0
	return(out)
}

#c. wrapper functions#

#makes a distance matrix from X using distance d#
distancematrix<-function(X,d,na.rm=TRUE){
	X<-as.matrix(X)
	if (d=="cosangle") return(disscosangle(X,na.rm))
	if (d=="abscosangle") return(dissabscosangle(X,na.rm))
	if (d=="euclid") return(disseuclid(X,na.rm))
	if (d=="abseuclid") return(dissabseuclid(X,na.rm))
	if (d=="cor") return(disscor(X,na.rm))
	if (d=="abscor") return(dissabscor(X,na.rm))
	#insert your own distance function here
	stop("Distance metric ",d," not available")
}

#makes a distance vector from X and y using distance d#
distancevector<-function(X,y,d,na.rm=TRUE){
	X<-as.matrix(X)
	y<-as.vector(y)
	if (d=="cosangle") return(vdisscosangle(X,y,na.rm))
	if (d=="abscosangle") return(vdissabscosangle(X,y,na.rm))
	if (d=="euclid") return(vdisseuclid(X,y,na.rm))
	if (d=="abseuclid") return(vdissabseuclid(X,y,na.rm))
	if (d=="cor") return(vdisscor(X,y,na.rm))
	if (d=="abscor") return(vdissabscor(X,y,na.rm))
	#insert your own distance function here
	stop("Distance metric ",d," not available")
}

#d. conversions#

#converts distance matrix to a vector#
dissvector<-function(M){
	if(!is.matrix(M))
		stop("arg to dissvector() must be a matrix")
	dM<-dim(M)
	if(dM[1]!=dM[2])
		stop("arg to dissvector() not a sqaure matrix")
	p<-dM[1]
	count<-1
	v<-rep(0,p*(p-1)/2)
	for (i in 1:(p-1)){
		v[count:(count+p-i-1)]<-M[i,(i+1):p]
		count<-count+p-i
	}
	return(v)
}

#converts distance vector to a matrix#
dissmatrix<-function(v){
	if(!is.vector(v))
		stop("arg to dissmatrix() must be a vector")
	p<-(1+sqrt(1+8*length(v)))/2
	M<-matrix(0,nrow=p,ncol=p)
	count<-1
	for (i in 1:(p-1)){
		M[i,(i+1):p]<-v[count:(count+p-i-1)]
		count<-count+p-i
	}
	return(M+t(M))
}

#maps index of distance vector into the [i,j] of corresponding pxp distance matrix#
vectmatrix<-function(index,p){
 	count<-1
 	s<-p-1
 	while((index-s)>0){
		s<-s+(p-1-count)
 		count<-count+1
 	}
 	i<-count
 	j<-p-(s-index)
 	return(c(i,j))
}


#e. correlation ordering#

#computes correlation ordering
correlationordering<-function(dist){
	if(!is.matrix(dist))
		stop("arg to correlationordering() must be a matrix")
	p<-length(dist[1,])
	if(p!=length(dist[,1]))
		stop("arg to correlationordering() not a square matrix")
	a<-dissvector(dist)
	b<-dissvector(abs(matrix(1:p,nrow=p,ncol=p,byrow=TRUE)-matrix(1:p,nrow=p,ncol=p)))
	return(cor(a,b))
}

#optimizes correlation ordering
improveordering<-function(dist,echo=FALSE){
	if(!is.matrix(dist))
		stop("arg to improveordering() must be a matrix")
	p<-length(dist[1,])
	if(p!=length(dist[,1]))
		stop("arg to improveordering() not a square matrix")
	v<-correlationordering(dist)
	if(echo)
		cat("Old order:",v,"\n")
	ord<-neword<-(1:p)
	if(!is.na(v)){
		final<-0
		for(gap in (1:(p-1))){
			while(final==0){
				sum<-0
				for(j in (1:(p-gap))){
					temp<-neword[j]
					neword[j]<-neword[j+gap]
					neword[j+gap]<-temp
					vnew<-correlationordering(dist[neword,neword])
					if(vnew<v) 
						neword<-ord
					else{
				 		sum<-sum+1
						v<-vnew
					}
					ord<-neword
				}
				if(sum==0) 
					final<-1
			}
		}
		if(echo)
			cat("New order:",correlationordering(dist[neword,neword]),"\n")
	}
	return(neword)		
}

dplot<-function(dist,hopachobj,ord="final",col=heat.colors(12),main=NULL,xlab=NULL,ylab=NULL,labels=NULL,showclusters=TRUE,...){
	dist<-as.matrix(dist)
	p<-nrow(dist)
	if(ord!="none")
		main<-paste(main,"Ordered Distance Matrix",sep="\n")
	if(is.null(ylab))
		ylab=""
	if(is.null(xlab))
		xlab=""
	if(showclusters)
		boundary<-cumsum(hopachobj$clustering$size)+0.5
	ordering<-1:p
	if(ord=="final")
		ordering<-hopachobj$final$ord
	if(ord=="cluster")
			ordering<-hopachobj$clustering$ord
	distplot<-dist[ordering,ordering]
	diag(distplot)<-min(dist[upper.tri(dist)])
	par(mfrow=c(1,1),pty="s")
	image(1:p,1:p,distplot[,p:1],main=main,xlab=xlab,ylab=ylab,axes=FALSE,col=col,...)
	if(!is.null(labels)){
		labels<-as.character(labels)
		axis(2,labels=rev(labels[ordering]),at=1:p,cex.axis=0.75,col.axis=2,las=2)
		axis(1,labels=labels[ordering],at=1:p,cex.axis=0.75,col.axis=2,las=2)
	}
	box()
	if(showclusters)
		abline(v=boundary,h=p+1-boundary,lty=2)
}
#Hierarchical Ordered Partitioning and Collapsing Hybrid (HOPACH)#

#1. silhouette related calculations#

#a. silhouettes

#from medoids and distance matrix 
medstosil<-function(medoids,dist){
	if(!is.vector(medoids))
		stop("First arg to medstosil() must be a vector")
	if(!is.matrix(dist))
		stop("Second arg to medstosil() must be a matrix")
	p<-length(dist[1,])
	if(p!=length(dist[,1]))
		stop("Second arg to medstosil() not a square matrix")
	k<-length(medoids)
	if(k<2){
		warning("Less than 2 medoids - can't calculate silhouettes!")
		sil<-rep(NA,p)
		clust<-rep(1,p)
	}
	else{
		#get labels
		clust<-apply(dist[,medoids],1,order)[1,]
		#get clussizes (same order as medoids)
		clussize<-table(clust)
		#get sils
		avgdist<-matrix(0,nrow=p,ncol=k)
		a<-b<-NULL
		for(j in (1:k)){
			for(m in (1:k)){
				subdist<-dist[clust==j,clust==m]
				if(clussize[j]>1) 
					avgdist[clust==j,m]<-apply(as.matrix(subdist),1,sum)
				else
					avgdist[clust==j,m]<-sum(subdist)
				if(m==j){
					if(clussize[m]>1) 
						avgdist[clust==j,m]<-avgdist[clust==j,m]/(clussize[m]-1)
					else 
						avgdist[clust==j,m]<-0
				}
				else 
					avgdist[clust==j,m]<-avgdist[clust==j,m]/clussize[m]
			}
			a[clust==j]<-avgdist[clust==j,j]
			if(clussize[j]>1) 
				b[clust==j]<-apply(as.matrix(avgdist[clust==j,-j]),1,min)
			else 
				b[clust==j]<-min(avgdist[clust==j,-j])
		}
		sil<-(b-a)/pmax(a,b)
	}
	list(clust,sil)
}

#from labels and distance matrix
labelstosil<-function(labels,dist){
	if(!is.vector(labels))
		stop("First arg to labelstosil() must be a vector")
	if(!is.matrix(dist))
		stop("Second arg to labelstosil() must be a matrix")
	p<-length(dist[1,])
	if(p!=length(dist[,1]))
		stop("Second arg to labelstosil() not a square matrix")
	if(length(labels)!=p)
		stop("Distance matrix and labels dimensions do not agree in labelstosil()")
	unlabels<-sort(unique(labels))
	k<-length(unlabels)
	if(k<2){
		warning("Only one label - can't calculate silhouettes!")
		sil<-rep(NA,p)
	}
	else{
		clussize<-table(labels)
		#get sils
		avgdist<-matrix(0,nrow=p,ncol=k)
		a<-b<-NULL
		for(j in (1:k)){
			for(m in (1:k)){
				subdist<-dist[labels==unlabels[j],labels==unlabels[m]]
				if(clussize[j]>1) 
					avgdist[labels==unlabels[j],m]<-apply(as.matrix(subdist),1,sum)
				else 
					avgdist[labels==unlabels[j],m]<-sum(subdist)
				if(m==j){
					if(clussize[m]>1) 
						avgdist[labels==unlabels[j],m]<-avgdist[labels==unlabels[j],m]/(clussize[m]-1)
					else 
						avgdist[labels==unlabels[j],m]<-0
				}
				else 
					avgdist[labels==unlabels[j],m]<-avgdist[labels==unlabels[j],m]/clussize[m]
			}
			a[labels==unlabels[j]]<-avgdist[labels==unlabels[j],j]
			if(clussize[j]>1) 
				b[labels==unlabels[j]]<-apply(as.matrix(avgdist[labels==unlabels[j],-j]),1,min)
			else 
				b[labels==unlabels[j]]<-min(avgdist[labels==unlabels[j],-j])
		}
		sil<-(b-a)/pmax(a,b)
	}
	list(labels,sil)
}


#b. mean/median split silhoutte

labelstomss<-function(labels,dist,khigh=9,within="med",between="med",hierarchical=TRUE){
	if(!is.vector(labels))
		stop("First arg to labelstomss() must be a vector")
	if(!is.matrix(dist))
		stop("Second arg to labelstomss() must be a matrix")
	p<-length(dist[1,])
	if(p!=length(dist[,1]))
		stop("Second arg to labelstomss() not a square matrix")
	if(length(labels)!=p)
		stop("Distance matrix and labels dimensions do not agree in labelstomss()")
	unlabels<-sort(unique(labels))
	k<-length(unlabels)
	ss<-NULL
	for(i in 1:k){
		labs<-(1:p)[labels==unlabels[i]]
		pp<-length(labs)
		if(pp<3)
			ss[i]<-NA
		else{
			dissvec<-dissvector(dist[labs,labs])
			bestk<-silcheck(dissvec,min(khigh,(length(labs)-1),na.rm=TRUE),diss=TRUE)[1]
			if(within=="med") 
				ss[i]<-median(pam(dissvec,bestk,diss=TRUE)$silinfo$widths[,3])
			if(within=="mean") 
				ss[i]<-pam(dissvec,bestk,diss=TRUE)$silinfo$avg.width
		}
	}
	if(sum(is.na(ss))==k)
		out<-NA
	else{
		if(hierarchical==TRUE){
			parentlab<-trunc(labels/10)
			unparentlab<-sort(unique(parentlab))
			pk<-length(unparentlab)
			pss<-NULL
			for(i in 1:pk){
				if(between=="med") 
					pss[i]<-median(ss[trunc(unlabels/10)==unparentlab[i]],na.rm=TRUE)
				if(between=="mean") 
					pss[i]<-mean(ss[trunc(unlabels/10)==unparentlab[i]],na.rm=TRUE)
			}
			if(between=="med") 
				out<-median(pss,na.rm=TRUE)
			if(between=="mean") 
				out<-mean(pss,na.rm=TRUE)
		}
		else{
			if(between=="med") 
				out<-median(ss,na.rm=TRUE)
			if(between=="mean") 
				out<-mean(ss,na.rm=TRUE)
		}
	}
	return(out)
}

#c. optimizing number of clusters with average silhouette or mss

#silcheck (silhouettes)
silcheck<-function(data,kmax=9,diss=FALSE,echo=FALSE,graph=FALSE){
	sil<-NULL
	m<-min(kmax,max((!diss)*(dim(data)[1]-1),(diss)*(-0.5+sqrt(1+8*length(data))),na.rm=TRUE))
	if(m<2)
		out<-c(1,NA)
	else{
		for(i in 1:(m-1))
			sil[i]<-pam(data,k=(i+1),diss=diss)$silinfo$avg.width
		if(echo)
			cat("best k = ",order(sil)[length(sil)]+1,", sil(k) = ",round(max(sil),4),"\n")
		if(graph){
			plot(2:m,sil,type="n",xlab="Number of Clusters",ylab="Average Silhouette")
			text(2:m,sil,2:m)
		}
		out<-c(order(sil)[length(sil)]+1,max(sil))
	}
	return(out)
}

#msscheck (mss)
msscheck<-function(dist,kmax=9,khigh=9,within="med",between="med",force=FALSE,echo=FALSE,graph=FALSE){
	if(!is.matrix(dist))
		stop("First arg to msscheck() must be a matrix")
	p<-length(dist[1,])
	if(p!=length(dist[,1]))
		stop("First arg to msscheck() not a square matrix")
	if(p<3)
		out<-c(1,NA)
	else{
		dvec<-dissvector(dist)
		kmax<-min(kmax,p-1,na.rm=TRUE)
		if(force)
			mss<-0
		else
			mss<-labelstomss(rep(1,p),dist,khigh,within,between)
		for(k in 2:kmax)
			mss[k]<-labelstomss(pam(dvec,k,diss=TRUE)$clust,dist,khigh,within,between)
		shift<-0
		if(force){
			mss<-mss[-1]
			shift<-1
		}
		if(echo)
			cat("best k = ",order(mss)[1]+shift,", mss(k) = ",round(min(mss),4),"\n")
		if(graph){
			kmin<-ifelse(force,2,1)
			plot(kmin:kmax,mss,type="n",xlab="Number of Clusters",ylab="MSS")
			text(kmin:kmax,mss,kmin:kmax)
		}
		out<-c(order(mss)[1]+shift,min(mss))
	}
	return(out)
}

#2. functions for making the tree#

#a. msssplitcluster: splits a cluster# 
	#clust1 is gene by subjects dataframe for cluster1
	#l1 is an integer label of cluster 1: e.g 11,12,13,22, etc describing its path so far in the tree.
	#id1 is id's of cluster1 indicating row numbers in original data frame subdata
	#kmax is the maximum number of groups
	#khigh  is the maximum number of child groups for each group when computing mss
	#medoid1 is the row number (in full data matrix!) indicating  medoid for cluster 1
	#medtodist is the distance from each gene in clust1 to the neighboring cluster medoid,
	#right is 1 if medoid2 is to the right and is 0 if clust1 is the last cluster so that medoid2
	# is the medoid2 to the left of clust1
	#silh is the silhouette of clust1
	#dist1 is the distance matrix for all genes in clust1 
msssplitcluster<-function(clust1,l1,id1,medoid1,med2dist,right,dist1,kmax=9,khigh=9,within="med",between="med"){
	if(!medoid1) 
		warning("Medoid missing - continue to split cluster")
	else{
		if(sum(medoid1==id1)==0 & medoid1) 
			warning("Medoid not in cluster - continue to split cluster")
	}
	if(is.matrix(clust1)){			
		p1<-length(clust1[,1])
		n<-length(clust1[1,])
	}
	else p1<-1
	if(p1<3) 
		k1<-1				
	else{
		l<-length(clust1[,1])
		dissvec<-dissvector(dist1)
		kmax<-min(p1-1,kmax,na.rm=TRUE)
		khigh<-min(p1-1,khigh,na.rm=TRUE)
		k1<-msscheck(dist1,kmax,khigh,within,between)[1]
		if(k1>1){
			pamobj<-pam(dissvec,k1,diss=TRUE)
			newclussizes<-pamobj$clusinfo[,1]
			newmedoids1<-id1[pamobj$medoids]  
			newlabels1<-pamobj$clustering
			distnewmedoids<-NULL 
			for(j in (1:k1)) 
				distnewmedoids[j]<-mean(med2dist[newlabels1==newlabels1[pamobj$medoids[j]]])  
			if(right==1) 
				ord<-rev(order(distnewmedoids))
			else 
				ord<-order(distnewmedoids)
			newmedoids1<-newmedoids1[ord]
			newclussizes<-newclussizes[ord]
			oldlab<-newlabels1
			for(j in (1:k1)) 
				newlabels1[oldlab==ord[j]]<-j
			newlabels1<-rep(10*l1,l)+newlabels1
		}
	}
	if(k1==1){
		newmedoids1<-medoid1
		newlabels1<-rep(10*l1,p1)
		newclussizes<-p1
	}
	for(a in (1:length(newmedoids1))){
		if(sum(newmedoids1[a]==id1)==0) 
			warning("Problem with new medoids after splitting cluster")
	}
	list(k1,newmedoids1,newlabels1,newclussizes)
}

#b. mssnextlevel: calls mssspltitcluster to produce the next level of the tree#
	#data is the data frame
	#prevlevel is the previous level of the tree
	#dmat is the distance matrix
	#kmax is the maximum number of groups
	#khigh is the maximum number of child groups for each group when computing mss
	#within and between are either "med" for median split silhouette or "mean"
	# for mean split silhouette
mssnextlevel<-function(data,prevlevel,dmat,kmax=9,khigh=9,within="med",between="med"){
	if(!is.matrix(data))
		stop("Frist arg to mssnextlevel() must be a matrix")
	if(!is.matrix(dmat))
		stop("Third arg to mssnextlevel() must be a matrix")
	n<-length(data[1,])
	p<-length(data[,1])
	if(length(dmat[1,])!=p)
		stop("Data and distance matrix dimensions do not agree in mssnextlevel()")
	if(length(dmat[,1])!=p)
		stop("Third arg to mssnextlevel() is not a square matrix")
	id<-1:p
	k<-prevlevel[[1]]
	medoids<-prevlevel[[2]]
	labels<-prevlevel[[4]]
	newk<-0
	newlabels<-newmedoids<-newclussizes<-NULL
	count<-1
	ordlabels<-sort(unique(labels))
	if(length(ordlabels)!=k) 
		warning("Number of unique labels not equal number of clusters in mssnextlevel()")
	if(sum(is.na(medoids))){
		warning("Missing values in medoid vector in nextlevel()")
		medoids[is.na(medoids)]<-FALSE
	}
	if(length(unique(medoids))<k && sum(medoids)) 
		warning("Medoids not unique in mssnextlevel()") 
	checkmeans<-FALSE
	if(length(medoids)==1 && !medoids){
		warning("No medoids provided in mssnxtlevel()")
		usemean<-TRUE
	}
	else{
		if(sum(medoids>1)==k) 
			usemean<-FALSE
		else 
			checkmeans<-TRUE
	}
	for(j in (1:k)){
		clust1<-data[labels==ordlabels[j],]	
		id1<-id[labels==ordlabels[j]]
		if(length(id1)>1) 
			clust1<-as.matrix(clust1)
		l1<-ordlabels[j]
		right<-(j<k)
		medoid1<-ifelse(is.na(medoids[j]),0,medoids[j])
		if (j<k) 
			medoid2<-medoids[j+1]
		else 
			medoid2<-medoids[j-1]
		if(length(id1)>1) 
			med2dist<-apply(as.matrix(dmat[labels==ordlabels[j],labels==labels[medoid2]]),1,mean)
		else 
			med2dist<-mean(dmat[labels==ordlabels[j],labels==labels[medoid2]])
		splitobj<-msssplitcluster(clust1,l1,id1,medoid1,med2dist,right,dmat[labels==l1,labels==l1],kmax,khigh,within,between)
		newlabels[labels==ordlabels[j]]<-splitobj[[3]]
		k1<-splitobj[[1]]
		start<-count
		end<-count+k1-1
		newmedoids[start:end]<-splitobj[[2]]
		newclussizes[start:end]<-splitobj[[4]]
		count<-count+k1		
	}
	count<-newk<-count-1
	newmedoids<-newmedoids[1:newk]
	newclussizes<-newclussizes[1:newk]
	final<-0
	if(count==k)
		final<-1
	if(max(newclussizes)==3) 
		final<-1
	list(newk,newmedoids,newclussizes,newlabels,final,rbind(prevlevel[[6]],cbind(sort(unique(newlabels)),newmedoids)))
}

#c. orderelements: produces an ordering of elements within a set of clusters 
	#level is a level of the tree
	#data is the data frame
	#rel is an indicator of whether to order elements in each cluster with respect 
	# to their own medoid ("own") or the neighboring medoid to the right ("neighbor")
	# or using improveordering() function ("co"). the default is "own"
	#d is an indicator of which distance function to use
	# choices are: "cosangle" (default),"abscosangle","euclid","abseuclid","cor","abscor".
	#dmat is the distance matrix. if this has already been calculated by the user, it can
	# be passed into the function in order to save calculation time
orderelements<-function(level,data,rel="own",d="cosangle",dmat=NULL){
	if(!is.matrix(data))
		stop("Second arg to orderelements() must be a matrix")
	idn<-1:length(data[,1])
	k<-level[[1]]
	labels<-level[[4]]
	medoids<-level[[2]]
	clussizes<-level[[3]]
	ord<-order(labels)
	idnord<-idn[ord]
	subdataord<-data[ord,]
	if(is.matrix(dmat))
		distord<-dmat[ord,]
	labelsord<-labels[ord]
	count<-1
	for(j in (1:k)){
		start<-count
		end<-count+clussizes[j]-1
		if(clussizes[j]>2){
			tempid<-idnord[start:end]
			if(rel=="co"){
				if(is.matrix(dmat)) 
					distj<-distord[,ord][start:end,start:end]
				else 
					distj<-distancematrix(subdataord[start:end,],d)
				idnord[start:end]<-tempid[improveordering(distj)]
			}
			else{
				if(rel=="neighbor"){
					if(j<k) 
						mednext<-medoids[j+1]
					else 
						mednext<-medoids[j-1]
				}
				else 
					mednext<-medoids[j]
				if(is.matrix(dmat)) 
					dmednext<-distord[start:end,mednext]	
				else 
					dmednext<-distancevector(subdataord[start:end,],as.vector(data[mednext,]),d)
				if(rel=="neighbor"){
					if(j<k) 
						ordtemp<-rev(order(dmednext))
					else 
						ordtemp<-order(dmednext)
				}
				else 
					ordtemp<-order(dmednext)
				idnord[start:end]<-tempid[ordtemp]
			}
		}
		else 
			idnord[start:end]<-idnord[start:end]
		count<-count+clussizes[j]
	}
	list(data[idnord,],idnord)
}

#d. mssinitlevel: creates ordered initial level#
	#data is the data matrix
	#kmax is the maximum number of groups
	#khigh is the maximum number of child groups for each group when computing mss
	#d is an indicator of which distance function to use.
	# choices are: "cosangle" (default),"abscosangle","euclid","abseuclid","cor","abscor"	
	#dmat is the distance matrix. if this has already been calculated by the user, it can
	# be passed into the function in order to save calculation time
	#within and between are either "med" for median split silhouette or "mean"
	# for mean split silhouette
	#ord is an indicator of how to order the clusters. choices are to maximize 
	# correlation ordering ("co") or to build a tree of cluster medoids ("clust")
mssinitlevel<-function(data,kmax=9,khigh=9,d="cosangle",dmat=NULL,within="med",between="med",ord="co"){
	if(!is.matrix(data))
		stop("First arg to mssinitlevel() must be a matrix")
	p<-length(data[,1])
	if(!is.matrix(dmat))
		dmat<-distancematrix(data,d)
	if(length(dmat[1,])!=p)
		stop("Data and distance matrix dimensions do not agree in mssinitlevel()")
	if(length(dmat[,1])!=p)
		stop("Distance matrix must be a square matrix in mssinitlevel()")
	m<-msscheck(dmat,kmax,khigh,within,between)
	if(m[1]==1){
		cat("No strong evidence for clusters in the first level - \n continuing to split root node anyway. \n")
		m<-msscheck(dmat,kmax,khigh,within,between,force=TRUE)
	}
	pamobj<-pam(dissvector(dmat),m[1],diss=TRUE)
	rowmedoids<-pamobj$medoids
	final<-ifelse(max(pamobj$clusinfo[,1])<3,1,0)
	if(m[1]>2){			
		medoidsdata<-as.matrix(data[rowmedoids,])
		l<-length(rowmedoids)
		medoidsdist<-dmat[rowmedoids,rowmedoids]
		if(ord=="co")
			medoidsord<-improveordering(medoidsdist)
		if(ord=="clust"){
			mpamobj<-pam(dissvector(medoidsdist),2,diss=TRUE)
			labelsmed<-mpamobj$clustering
			medmed<-mpamobj$medoids
			clussizes<-mpamobj$clusinfo[,1]
			prevlevel<-mssnextlevel(medoidsdata,list(2,medmed,clussizes,labelsmed,0,cbind(c(1,2),medmed),dmat=medoidsdist,kmax,khigh,within,between))
			final<-prevlevel[[5]]
			if(final==0){
				depth<-(l-1)
				for(j in (1:depth)){
					if(final==0){
						clustnext<-mssnextlevel(medoidsdata,prevlevel,dmat=medoidsdist,kmax,khigh,within,between)
						final<-clustnext[[5]]
					}
					if(final==1){ 
						j<-depth
						prevlevel<-clustnext
					}
				}
			}					
			medoidsord<-orderelements(prevlevel,medoidsdata,rel="neighbor",d=d,dmat=medoidsdist)[[2]]
		}
		k<-m[1]
		rowmedoids<-rowmedoids[medoidsord]
		labels<-lab2<-pamobj$clustering
		for(j in (1:k)) 
			lab2[labels==medoidsord[j]]<-j
		output<-list(k,rowmedoids,pamobj$clusinfo[,1][medoidsord],lab2,final,cbind(1:k,rowmedoids))
	}
	else
		output<-list(2,pamobj$medoids,pamobj$clusinfo[,1],pamobj$clustering,final,cbind(1:2,pamobj$medoids))
	return(output)
}

#e. collapsing functions#
	#paircoll() collapses a pair of medoids (i,j)
	#collap() calls paircoll() to consider and possibly perform collapsing
	#msscollap() collapses by sequentially calling collap starting with 
	# the closest pair of clusters til there is no more improvement in mss
	#mssmulticollap tries all pairs of clusters and collapses any that improve mss
	###########################################################################################
	#data is the data matrix
	#level is level of the tree
	#d is an indicator of which distance function to use
	# choices are: "cosangle" (default),"abscosangle","euclid","abseuclid","cor","abscor"
	#dmat is the distance matrix. if this has already been calculated by the user, it can
	# be passed into the function in order to save calculation time.
	#newmed is an indicator of which way to find the medoid of the new cluster after collapsing.
	# choices are: "nn" to use the nearest neighbor of the clustersize-weighted 
	# mean of the two medoids as the medoid of a collapsed cluster, "uwnn" to use an unweighted 
	# version of nearest neighbor so that each cluster (rather than each gene) gets equal
	# weight in the mean, "center" to use the cluster center (element with min sum distance 
	# to all others), "medsil" (default) to use the medoid which maximizes the medoid based 
	# silhouette (i.e.: (a-b)/max(a,b), where a=dist(medoid), b=dist(next closest medoid)). 
	#the silhouettes and splits (arguments [[1]] and [[2]] of level) refer to the original
	# splits and loose their meaning if the child cluster(s) are collapsed
	#impr is a margin of improvement required to accept a collapse with msscollap and
	# mssmulticollap. the default is impr=0
paircoll<-function(i,j,data,level,d="cosangle",dmat=NULL,newmed="medsil"){
	if(!is.matrix(data))
		stop("First arg to paircoll() must be a matrix")
	p<-length(data[,1])
	k<-level[[1]]
	labels<-level[[4]]
	medoids<-level[[2]]
	clussizes<-level[[3]]
	N<-length(level[[6]][,1])
	block<-level[[6]][(N-k+1):N,]
	if(N==k)
		prevblock<-NULL
	else
		prevblock<-level[[6]][1:(N-k),]
	if(max(i,j)>k)
		stop("The clusters to collapse do not exist in paircoll()") 
	labeli<-labels[medoids[i]]
	labelj<-labels[medoids[j]]
	oldlabels<-labels
	labels[labels==labelj]<-labeli
	trunclabels<-trunc(oldlabels/10)
	labelparents<-unique(trunclabels)
	parenti<-order(labelparents)[labelparents==trunc(labeli/10)]
	parentj<-order(labelparents)[labelparents==trunc(labelj/10)]
	if(newmed=="nn")
		fakemed<-(data[medoids[i],]*clussizes[i]+data[medoids[j],]*clussizes[j])/(clussizes[i]+clussizes[j])
	if(newmed=="uwnn")
		fakemed<-(data[medoids[i],]+data[medoids[j],])/2
	if(newmed=="nn" || newmed=="uwnn"){ 
		rowsub<-(1:p)[labels==labeli]
		distsfm<-distancevector(data[rowsub,],as.vector(fakemed),d)
		medoids[i]<-rowsub[order(distsfm)[1]]
	}
	else{
		if(is.matrix(dmat)) 
			colldist<-dmat[labels==labeli,labels==labeli]
		else 
			colldist<-distancematrix(data[labels==labeli,],d)
		rowsub<-(1:p)[labels==labeli]
		if(newmed=="center"){
			sumdist<-apply(colldist,1,sum)
			medoids[i]<-rowsub[order(sumdist)==1]
		}
		if(newmed=="medsil"){
			othermed<-medoids[-c(i,j)]
			collp<-length(labels[labels==labeli])
			othern<-length(othermed)
			if(othern==0)
				stop("Not enough medoids to use newmed='medsil' in paircoll()")
			if(is.matrix(dmat)){
			 	if(othern==1) 
					otherdist<-cbind(dmat[labels==labeli,othermed])
				else 
					otherdist<-rbind(dmat[labels==labeli,othermed])
			}
			else{
				if(othern==1)
					otherdist<-distancevector(data[labels==labeli,],data[othermed,],d)
				else{
					othermedmat<-data[othermed,]
					otherdist<-matrix(0,nrow=collp,ncol=othern)
					for(l in 1:othern)
						otherdist[,l]<-distancevector(data[labels==labeli,],othermedmat[l,],d)
				}
			}			
			if(othern==1)
				b<-otherdist
			else
				b<-apply(otherdist,1,min)
			b<-matrix(b,nrow=collp,ncol=collp)
			diag(b)<-0
			b<-abs(b-colldist)/pmax(colldist,b)
			sumdist<-apply(b,1,sum)
			medoids[i]<-rowsub[rev(order(sumdist))==1]
		}
	}	
	k<-k-1	
	clussizes[i]<-clussizes[i]+clussizes[j]
	block[i,2]<-medoids[i]
	if(j<=k){			
		for(l in (j:k)){
			medoids[l]<-medoids[l+1]
 		        clussizes[l]<-clussizes[l+1]
			block[l,]<-block[l+1,]
		}
	}
	medoids<-medoids[1:k]
	clussizes<-clussizes[1:k]
	block<-block[1:k,]
	return(list(k,medoids,clussizes,labels,level[[5]],rbind(prevblock,block)))
}

#note: this version of collap does not have silhbased arg: for use with MSS only (not silhouettes)
collap<-function(data,level,d="cosangle",dmat=NULL,newmed="medsil"){
	if(!is.matrix(data))
		stop("First arg to collap() must be a matrix")
	k<-level[[1]]
	if(k<3){
		warning("Not enough medoids to use newmed='medsil' in collap() - \n using newmed='nn' instead \n") 
		newmed<-"nn"
	}
	medoids<-level[[2]]
	clussizes<-level[[3]]
	if(sum(is.na(clussizes))) 
		warning("NA in clussizes")
	medoidsdata<-data[medoids,]
	if(sum(is.na(medoidsdata))>0) 
		warning("Missing value(s) in medoidsdata in collap()")
	if(is.matrix(dmat)) 
		distmed<-dmat[medoids,medoids]
	else 
		distmed<-distancematrix(medoidsdata,d)
	distv<-dissvector(distmed)
	indexmin<-order(distv)[1]
	best<-vectmatrix(indexmin,k)
	clustfinal<-paircoll(best[1],best[2],data,level,d,dmat,newmed)
	return(clustfinal)
}

msscollap<-function(data,level,khigh,d="cosangle",dmat=NULL,newmed="medsil",within="med",between="med",impr=0){
	if(!is.matrix(data))
		stop("First arg to msscollap() must be a matrix")	
	if(!is.matrix(dmat)) 
		dmat<-distancematrix(data,d)
	newk<-level[[1]]
	mss1<-labelstomss(level[[4]],dmat,khigh,within,between)
	maxncoll<-max(0,newk-2)
        ncoll<-0
	coll<-1
        if(newk<=2) 
		coll<-0
 	while((coll==1) && (ncoll<= maxncoll)){
                levelc<-collap(data,level,d,dmat,newmed)
		mss2<-labelstomss(levelc[[4]],dmat,khigh,within,between)
		if(mss1==0) 
			r<-0
		else 
			r<-(mss1-mss2)/mss1
		if(r<impr) 
			coll<-0 	 
    		else{
			mss1<-mss2
			level<-levelc
                        ncoll<-ncoll+1
                }
	}
	return(level)
} 

mssmulticollap<-function(data,level,khigh,d="cosangle",dmat=NULL,newmed="medsil",within="med",between="med",impr=0){
	if(!is.matrix(data))
		stop("First arg to mssmulticollap() must be a matrix")
	if(!is.matrix(dmat)) 
		dmat<-distancematrix(data,d)
	medoids<-level[[2]]
	medoidsdata<-data[medoids,]
	if(sum(is.na(medoidsdata))>0) 
		warning("Missing value(s) in medoidsdata in mssmulticollap()")
	distmed<-dmat[medoids,medoids]
	k<-level[[1]]
	ord<-order(dissvector(distmed))
	mss1<-labelstomss(level[[4]],dmat,khigh,within,between)
	maxncoll<-max(0,k*(k-1)/2)
	ncoll<-0 
	i<-1
	while(i<=maxncoll){
		clusts<-vectmatrix(ord[i],k)
		levelc<-paircoll(clusts[1],clusts[2],data,level,d,dmat,newmed)
		mss2<-labelstomss(levelc[[4]],dmat,khigh,within,between)
		r<-(mss1-mss2)/mss1
		if(r>=impr){
			mss1<-mss2
			level<-levelc
                        ncoll<-ncoll+1
			k<-level[[3]]
			maxncoll<-max(0,k*(k-1)/2)
			i<-0
			medoids<-level[[2]]
			medoidsdata<-data[medoids,]
			if(sum(is.na(medoidsdata))>0) 
				warning("Missing value(s) in medoidsdata in mssmulticollap()")
			distmed<-dmat[medoids,medoids]
			ord<-order(dissvector(distmed))
	        }
		i<-i+1
	}
	return(level)
}

#f. iterating functions to run down the tree#
	#mssrundown() runs down the tree K levels with a stopping rule to 
	# find the main clusters
	#msscomplete() runs down the tree to the final level from any level
	#digits() determines the level of the tree from the labels
	#############################################################################################
	#data is the data matrix
	#K is the maximum number of levels to compute
	#kmax is the maximum number of groups
	#khigh  is the maximum number of child groups for each group when computing mss
	#d is an indicator of which distance function to use
	# choices are: "cosangle" (default),"abscosangle","euclid","abseuclid","cor","abscor"	
	#dmat is the distance matrix. if this has already been calculated by the user, it can
	# be passed into the function in order to save calculation time
	#coll is an indicator of how to collapse. the choices are to begin with the closest 
	# pair of clusters and collapse til there is no more improvement in mss ("seq") 
	# or to try all pairs of clusters and accept any collapse that improves mss ("all").
	#newmed is an indicator of which way to find the medoid of the new cluster after collapsing.
	# choices are: "nn" to use the nearest neighbor of the clustersize-weighted 
	# mean of the two medoids as the medoid of a collapsed cluster, "uwnn" to use an unweighted 
	# version of nearest neighbor so that each cluster (rather than each gene) gets equal
	# weight in the mean, "center" to use the cluster center (element with min sum distance 
	# to all others), "medsil" (default) to use the medoid which maximizes the medoid based 
	# silhouette (i.e.: (a-b)/max(a,b), where a=dist(medoid), b=dist(next closest medoid)). 
	#stop is an indicator that the tree should stop as soon as there is an increase in 
	# mss moving to the next level
	#finish is an indicator that the tree should compute all K levels and return the
	# one with the minimum mss (when finish==FALSE and stop==FALSE, level K is returned)
	#within and between are either "med" for median split silhouette or "mean"
	# for mean split silhouette.
	#impr is a margin of improvement required to accept a collapse with msscollap and
	# mssmulticollap. the default is impr=0.
mssrundown<-function(data,K=16,kmax=9,khigh=9,d="cosangle",dmat=NULL,initord="co",coll="seq",newmed="medsil",stop=TRUE,finish=FALSE,within="med",between="med",impr=0){
	if(!is.matrix(data))
		stop("First arg to mssrundown() must be a matrix")
	if(!is.matrix(dmat))
                dmat<-distancematrix(data,d)
	bestlevel<-level<-mssinitlevel(data,kmax,khigh,d,dmat,within,between,initord)
	bestmss<-mss<-labelstomss(level[[4]],dmat,khigh,within,between)
	bestl<-l<-1
	ind<-0
	cat("Searching for main clusters... \n")
	if(level[[5]]==1)
		return(level)
	while((l<=K) && (ind==0)){
		cat("Level ",l,"\n")
		if(coll=="seq")	
			levelc<-msscollap(data,level,khigh,d,dmat,newmed,within,between,impr)
		if(coll=="all") 
			levelc<-mssmulticollap(data,level,khigh,d,dmat,newmed,within,between,impr)
		mss<-labelstomss(levelc[[4]],dmat,khigh,within,between)
		if(mss>=bestmss & stop==TRUE)
			ind<-1
		else{
			if(mss<bestmss){
				bestlevel<-levelc
				bestmss<-mss
				bestl<-l
			}
		}
		l<-l+1
		if(l<=K){	
			level<-mssnextlevel(data,levelc,dmat,kmax,khigh,within,between)
			if(finish==TRUE){
				if(sum(trunc(level[[4]]/10)*10==level[[4]])==length(level[[4]]) & l<=K){
					ind<-1
					bestlevel<-levelc
					bestmss<-mss
					bestl<-(l-1)
				}
			}
		}
	}
	cat("Identified",bestlevel[[1]]," main clusters in level",bestl,"with MSS =",bestmss,"\n")
	return(bestlevel)
}

msscomplete<-function(level,data,K=16,khigh=9,d="cosangle",dmat=NULL,within="med",between="med"){
	if(!is.matrix(data))
		stop("First arg to msscomplete() must be a matrix")
	if(!is.matrix(dmat))
                dmat<-distancematrix(data,d)
	count<-digits(level[[4]][1])
	cat("Running down without collapsing from Level",count,"\n")
	while((max(level[[3]])>3) & (count<K)){
		level<-newnextlevel(data,level,dmat,2,khigh)
		count<-count+1
		cat("Level",count,"\n")
	}
	return(level)
}

digits<-function(label){
	label<-label[1]
	count<-0
	while(label>=1){
		count<-count+1
		label<-label/10
	}
	return(count)
}

#newnextlevel and newsplitcluster are needed in msscomplete to rundown 
#completely without collapsing back to the main clusters.
	#data is the data matrix
	#prevlevel is the level from which a next level is produced
	#dmat is the distance matrix, as above
	#klow and khigh are the min and max number of children at each node
	#newnextlevel produces the args to newsplitcluster and calls
	# this function to do the splitting of each node
newnextlevel<-function(data,prevlevel,dmat,klow=2,khigh=6){
	if(!is.matrix(data))
		stop("First arg to newnextlevel() must be a matrix")
	p<-length(data[,1])
	n<-length(data[1,])
	if(length(dmat[1,])!=p)
		stop("Distance and data matrix dimensions do not agree in newnextlevel()")
	if(length(dmat[,1])!=p)
		stop("Third arg to newnextlevel() is not a square matrix")
	id<-1:p
	k<-prevlevel[[1]]
	medoids<-prevlevel[[2]]
	labels<-prevlevel[[4]]
	newk<-0
	newlabels<-newmedoids<-newclussizes<-NULL
	count<-1
	ordlabels<-sort(unique(labels))
	if(length(ordlabels)!=k) 
		warning("Number of unique labels not equal number of clusters in newnextlevel()")
	if(sum(is.na(medoids))){
		warning("Missing value(s) in medoid vector in newnextlevel()")
		medoids[is.na(medoids)]<-FALSE
	}
	if(length(unique(medoids))<k && sum(medoids)) 
		warning("Medoids in newnextlevel() are not unique") 
	checkmeans<-FALSE
	if(length(medoids)==1 && !medoids){
		warning("No medoids provided in newnextlevel()")
		usemean<-TRUE
	}
	else{
		if(sum(medoids>1)==k) 
			usemean<-FALSE
		else 
			checkmeans<-TRUE
	}
	for(j in (1:k)){
		clust1<-data[labels==ordlabels[j],]	
		id1<-id[labels==ordlabels[j]]
		if(length(id1)>1){
			kmax<-min(c(khigh,dim(clust1)[1]-1))
			clust1<-as.matrix(clust1)
		}
		l1<-ordlabels[j]
		right<-(j<k)
		medoid1<-ifelse(is.na(medoids[j]),0,medoids[j])
		if (j<k) 
			medoid2<-medoids[j+1]
		else 
			medoid2<-medoids[j-1]
		if(length(id1)>1) 
			med2dist<-apply(as.matrix(dmat[labels==ordlabels[j],labels==labels[medoid2]]),1,mean)
		else 
			med2dist<-mean(dmat[labels==ordlabels[j],labels==labels[medoid2]])
		splitobj<-newsplitcluster(clust1,l1,id1,klow,kmax,medoid1,med2dist,right,dmat[labels==l1,labels==l1]) 
		newlabels[labels==ordlabels[j]]<-splitobj[[3]]
		k1<-splitobj[[1]]
		start<-count
		end<-count+k1-1
		newmedoids[start:end]<-splitobj[[2]]
		newclussizes[start:end]<-splitobj[[4]]
		count<-count+k1-1+1		
	}
	count<-count-1
	newk<-count
	newmedoids<-newmedoids[1:newk]
	newclussizes<-newclussizes[1:newk]
	final<-0
	if(count==k)
		final<-1
	if(max(newclussizes)==3) 
		final<-1
	return(list(newk,newmedoids,newclussizes,newlabels,final,rbind(rbind(prevlevel[[6]],cbind(sort(unique(newlabels)),newmedoids)))))
}

newsplitcluster<-function(clust1,l1,id1,klow=2,khigh=2,medoid1,med2dist,right,dist1){
	if(!medoid1) 
		warning("Medoid missing - continue to split cluster")
	else{
		if(sum(medoid1==id1)==0 & medoid1) 
			warning("Medoid not in cluster - continue to split cluster")
	}
	if(is.matrix(clust1)){			
		p1<-length(clust1[,1])
		n<-length(clust1[1,])
	}
	else p1<-1
	if(p1<3){
		k1<-1				
		newmedoids1<-medoid1
		newlabels1<-rep(10*l1,p1)
		newclussizes<-p1
	}
	else{
		l<-length(clust1[,1])
		dissvec<-dissvector(dist1)
		kmax<-min(p1-1,khigh)
		a<-rep(0,(kmax-klow+2))
		best<-2
		for(j in (klow:kmax)){
			a[j]<-pam(dissvec,j,diss=TRUE)$silinfo$avg.width
			if (a[j]>a[best]) best<-j
		}
		k1<-best
		pamobj<-pam(dissvec,k1,diss=TRUE)
		newclussizes<-pamobj$clusinfo[,1]
		newmedoids1<-id1[pamobj$medoids]
		newlabels1<-pamobj$clustering
		distnewmedoids<-NULL
		for(j in (1:k1)) 
			distnewmedoids[j]<-mean(med2dist[newlabels1==newlabels1[pamobj$medoids[j]]])  
		if(right==1) 
			ord<-rev(order(distnewmedoids))
		else 
			ord<-order(distnewmedoids)
		newmedoids1<-newmedoids1[ord]
		newclussizes<-newclussizes[ord]
		oldlab<-newlabels1
		for(j in (1:k1)) 
			newlabels1[oldlab==ord[j]]<-j
		newlabels1<-rep(10*l1,l)+newlabels1
	}
	for(a in (1:length(newmedoids1))){
		if(sum(newmedoids1[a]==id1)==0) 
			warning("Problem with new medoids after splitting cluster")
	}
	return(list(k1,newmedoids1,newlabels1,newclussizes))
}

#g. wrapper function to build the whole tree with clusters#
	#data is the data matrix.
	#clusters (default="best") tells how to identify the main clusters 
	# clusters="greedy" stops at the first level where MSS increases
	# clusters="none" does not identify main clusters
	# clusters="best" identifies the level<=K with best MSS as main clusters
	#K is the maximum number of levels to compute. right now, still have 
	# computational problem with doing more than 16, which is the default.
	#kmax is the maximum number of groups.
	#khigh  is the maximum number of child groups for each group when computing mss.
	#d is an indicator of which distance function to use.
	# choices are: "cosangle" (default),"abscosangle","euclid","abseuclid","cor","abscor".	
	#dmat is the distance matrix. if this has already been calculated by the user, it can
	# be passed into the function in order to save calculation time.
	#coll is an indicator of how to collapse. the choices are to begin with the closest 
	# pair of clusters and collapse til there is no more improvement in mss ("seq") 
	# or to try all pairs of clusters and accept any collapse that improves mss ("all").
	#newmed is an indicator of which way to find the medoid of the new cluster after collapsing.
	# choices are: "nn" to use the nearest neighbor of the clustersize-weighted 
	# mean of the two medoids as the medoid of a collapsed cluster, "uwnn" to use an unweighted 
	# version of nearest neighbor so that each cluster (rather than each gene) gets equal
	# weight in the mean, "center" to use the cluster center (element with min sum distance 
	# to all others), "medsil" (default) to use the medoid which maximizes the medoid based 
	# silhouette (i.e.: (a-b)/max(a,b), where a=dist(medoid), b=dist(next closest medoid)). 
	#mss is either "med" (default) for median split silhouettes or "mean" for mean 
	# split silhouettes.
	#impr is a margin of improvement required to accept a collapse with msscollap and
	# mssmulticollap. the default is impr=0.
	#initord is "co" (default) if improveordering() is used to order the clusters in 
	# the first level or "clust" if clsutering the medoids is used.
	#ord determines how elements are ordered within clusters: "co" is 
	# using improveordering(), "own" is distance to their own medoid, and "nieghbor"
	# is distance to the neighboring medoid (to the right). 
hopach<-function(data,dmat=NULL,d="cosangle",clusters="best",K=16,kmax=9,khigh=9,coll="seq",newmed="medsil",mss="med",impr=0,initord="co",ord="own"){
	if(inherits(data,"exprSet")) 
		data<-exprs(data)
	data<-as.matrix(data)
	if(K>16){
		K<-16
		warning("K set to 16 - can't do more than 16 levels")
	}
	if(K<1){
		K<-1
		warning("K set to 1 - can't do less than 1 level")
	}
	if(clusters!="none"){
		cuttree<-mssrundown(data,K,kmax,khigh,d,dmat,initord,coll,newmed,stop=(clusters=="greedy"),finish=TRUE,within=mss,between=mss,impr)
		if(cuttree[[1]]>1) 
			cutord<-orderelements(cuttree,data,rel=ord,d,dmat)[[2]]
		else 
			cutord<-NULL
		out1<-list(k=cuttree[[1]],medoids=cuttree[[2]],sizes=cuttree[[3]],labels=cuttree[[4]],order=cutord)
		finaltree<-msscomplete(cuttree,data,K,khigh,d,dmat,within=mss,between=mss)
	}
	else{
		out1<-NULL
		finaltree<-msscomplete(mssinitlevel(as.matrix(data),kmax,khigh,d,dmat,within=mss,between=mss,initord),data,K,khigh,d,dmat,within=mss,between=mss)
	}
	dimnames(finaltree[[6]])<-list(NULL,c("label","medoid"))
	out2<-list(labels=finaltree[[4]],order=orderelements(finaltree,data,rel=ord,d,dmat)[[2]],medoids=finaltree[[6]])
	return(list(clustering=out1,final=out2,call=match.call(),metric=d))
}
.First.lib<-function(libname,pkgname){
    	require(cluster) || stop("can't load without cluster package")
}

