# Programm zur Durchführung einer Clusteranalyse

# V 0.5 vom 20./24.09.2007

# (c) Dr. Alexander Preuß

# Eine vollständige Datenmatrix X muss vorhanden sein!



# Bestimmung von I und J

I <- nrow(X)

J <- ncol(X)

# Bestimmung eines Vektors EinsI und EinsJ

EinsI <- c(rep(1, times=I))

EinsJ <- c(rep(1, times=J))


# Zählvariable für Cluster

cluster <- c(rep(0, times=I))

for (i in 1:I)

	{

	cluster[i]=i

	}

# Bestimmung der Spaltenmittel

smittel <- c(rep(0, times = J))
sstabw <- c(rep(0, times = J))

for (j in 1:J)

	{smittel[j] <- (EinsI%*%X[1:I,j])/I}

# Bestimmung der Spaltenstandardabweichung

for (j in 1:J)

	{

	for (i in 1:I)

		{sstabw[j] <- sstabw[j] + (X[i,j]-smittel[j])^2}

	sstabw[j] <- sqrt(sstabw[j]/(I-1))
	
	}

# Vorprüfung: Alle Fälle gleich?

if (sum(sstabw) == 0)

	{

	stop("Clusteranalyse kann nicht ausgeführt werden, da alle Fälle gleich sind!")

	}


# Bestimmung der Matrix Z (Standardwerte)

Z <- X

for (i in 1:I)

	{

	for (j in 1:J)

		{Z[i,j] <- (X[i,j]-smittel[j])/sstabw[j]}

	}


# Bestimmung der Distanzmatrix Dq (quadrierte euklidische Distanzen
# der standardisierten Werte)

Dq <- matrix(c(rep(0, times = I*I)), nrow=I)


for (i in 1:I)
	
	{
	
	for (j in 1:I)

		{

		for (k in 1:J)

			{

			Dq[i,j] <- Dq[i,j] + (Z[i,k]-Z[j,k])^2
			
			}

		}

	}


##### "Single linkage"




Fusion <- matrix(c(rep(0, times=(I-1)*4)), nrow = (I-1),
		     dimnames=list(seq(I-1),c("Schritt", "Cluster 1", "Cluster 2", "Distanz" )))

Dx <- Dq


Start <- 1



# Vorprüfung: Gibt es bereits identische Fälle?

loeschen <- c(rep(0,times=I))

for (i in 2:I)

	{
		
	j <- 1

	abbruch <- 0

	while (abbruch == 0)

		{ 

		if (Dx[i,j] == 0)

			{

			

			abbruch <- 1

			cmin <- min(i,j)
			cmax <- max(i,j)

			loeschen[cmin] <- cmax 

			Fusion[Start,1] <- Start
			Fusion[Start,2] <- cmin
			Fusion[Start,3] <- cmax
			Fusion[Start,4] <- Dx[i,j]

			Start <- Start + 1

			}

		j <- j + 1

		if (j > (i-1))

			{

			abbruch <- 1

			}

		}

	}
			
for (i in 1:I)

	{

	if (loeschen[i] > 0)

		{
	
		Dxx <- Dx
	
		for (j in 1:I)

			{
			
			Dxx[i,j] <- min(Dx[i,j],Dx[loeschen[i],j],Dx[j,i],Dx[j,loeschen[i]])
			Dxx[j,i] <- min(Dx[i,j],Dx[loeschen[i],j],Dx[j,i],Dx[j,loeschen[i]])

			Dxx[loeschen[i],j] <- 0
			Dxx[j,loeschen[i]] <- 0
				
			
			}

		Dx <- Dxx

		}

	}


# # # Start

for (n in Start:(I-1))

	{

	minimum <- 99999999

	for (i in 2:I)

		{
		
		for (j in 1:(i-1))

			{

			if (Dx[i,j] > 0)

				{
				
					if (Dx[i,j] < minimum)

					{

					minimum <- Dx[i,j]
					C1 <- cluster[i]
					C2 <- cluster[j]
				
					cmin <- min(i,j)
					cmax <- max(i,j)

					}
				}

			}

		}

	
	cluster[cmax] <- cmin
	
	Fusion[n,1] <- n
	Fusion[n,2] <- cmin
	Fusion[n,3] <- cmax
	Fusion[n,4] <- minimum

	
# Erstellen der neuen Distanzmatrix	

Dxx <- Dx

Dxx[cmin,cmax] <- 0
Dxx[cmax,cmin] <- 0

for (i in 1:I)

	{

	Dxx[cmax,i] <- 0
	Dxx[i,cmax] <- 0

	}

for (i in 1:I)

	{

	Dxx[i,cmin] <- min(Dx[i,cmin],Dx[i,cmax],Dx[cmin,i],Dx[cmax,i])
	Dxx[cmin,i] <- min(Dx[i,cmin],Dx[i,cmax],Dx[cmin,i],Dx[cmax,i])

	
	}

	
Dx <- Dxx


}


# Definition der Cluster


# Erzeugen eines Dendrograms



dendro <- c(rep(0,times=I))

dendro[1] <- Fusion[(I-1),3]

dendro[2] <- Fusion[(I-1),2]


for (i in 1:(I-2))

	{

	posc <- 0

	pos <- Fusion[(I-1-i),2]

	for (j in 1:I)

		{

		if (pos == dendro[j])

			{

			posc <- j

			}
		}

	if (posc > 0)

		{

		dendro[(posc+2):I] <- dendro[(posc+1):(I-1)]

		dendro[posc+1] <- Fusion[(I-1-i),3]

		}

	else

		{

		dendro[i+2] <- Fusion[(I-1-i),3]

		}

	}

# Ermittlung der Zahl der Cluster

cls <- sort(Fusion[1:(I-1),2])

clz <- 2

for (i in 2:(I-1))

	{

	if (cls[i] > cls[(i-1)])

		{

		clz <- clz + 1

		}

	}
	

	



# Grafische Ausgabe

windows()

for (i in 1:(I-1))

	{

	xx <- c(0, Fusion[i,4])

	
	# y1 suchen

	for (j in 1:I)

		{

		if (Fusion[i,2] == dendro[j])

			{ 

			yp1 <- j

			}
		
		if (Fusion[i,3] == dendro[j])

			{ 

			yp2 <- j

			}

		}


	y1 <- (yp1/I*I)

	y2 <- (yp2/I*I)

	yy1 <- c(y1, y1)

	yy2 <- c(y2, y2)

	xx3 <- c(Fusion[i,4], Fusion[i,4])

	yy3 <- c(y1, y2)

	
	maxdist <- max(Fusion[1:(I-1),4])


	plot(xx,yy1,type="l", xlim=c(-(maxdist/10),maxdist), ylim=c(0,I), xlab="Distanz", ylab=" ", axes=F)


	par(new=T)

	plot(xx,yy2,type="l", xlim=c(-(maxdist/10),maxdist), ylim=c(0,I), xlab="Distanz", ylab=" ", axes=F)

	par(new=T)
	
	plot(xx3,yy3,type="l", xlim=c(-(maxdist/10),maxdist), ylim=c(0,I), xlab="Distanz", ylab=" ", axes=F)

	par(new=T)

	text(-(maxdist/10),i,rnamen[dendro[i]])

	par(new=T)	

	}

text(-(maxdist/10),I,rnamen[dendro[I]])

title("Dendrogramm")

axis(1)

par(new=T)


# Ermittlung der Clusterzahl

Clusteruebersicht <- matrix(c(rep(0,times=(I*(I+1)))),nrow=(I+1))

for (i in 1:I)

	{

	Clusteruebersicht[(i+1),1] <- i

	}

for (i in 1:(I-1))

	{

	Clusteruebersicht[1,(i+1)] <- Fusion[i,4]

	for (j in 1:I)

		{

		if (Fusion[i,3] == j)

			{

			Clusteruebersicht[(j+1),(i+1)] <- 1

			}

		}

	}



for (i in 1:I)

	{

	for (j in 1:(I-1))

		{

		if (Clusteruebersicht[(i+1),(I+1-j)] == 1)

			{

			Clusteruebersicht[(i+1),(I-j)] <- 1

			}

		}

	}

distanz <- matrix(c(rep(0,(2*(I-1)))), ncol=2)

for (i in 1:(I-1))

	{

	distanz[i,1] <- (Clusteruebersicht[1,i+1]/maxdist)
	distanz[i,2] <- sum(Clusteruebersicht[2:(I+1),(i+1)])

	}