#Hierarchical clustering By Joe Irgon
with(LinearAlgebra):

initT:=proc(d1) global T,TS,n,d:
d:=d1;
n:=RowDimension(d);
T:=[seq(i,i=1..n)]:
TS:=T:
Heirarch():

end:

Heirarch:=proc() local i,j,V,d1,cluster:
global T,TS,d:


while nops(T)<>1 do
	
cluster:=MinLoc():
TS:=subsop(cluster[1]=NULL,cluster[2]=[op(T[cluster[1]]),op(T[cluster[2]])],T);
T:=subsop(cluster[1]=NULL,cluster[2]=[T[cluster[1]],T[cluster[2]]],T);

od:
T:
end:

#Find min cluster distances 
MinLoc:=proc() local i,j,indx,rec,dd,k,m,n:
global T,TS,d:
	n:=nops(T):
	rec:=infinity:	

	for i to n do
   	   for j to i - 1 do
		dd:=0;
		for k in TS[i] do
		   for m in TS[j] do
			dd:=dd+d[k,m];
		   od:
		od:
		dd:=dd/(nops(TS[i])*nops(TS[j]));
	   	if dd < rec then
		 rec:=dd:
		 indx:=[i,j]:
		fi:
	   od:
	od:
indx:
end:

initT:=proc(d1) global T,d:
local n:
d:=d1;
n:=RowDimension(d);
T:=[seq(i,i=1..n)]:

Heirarch():

end:

Heirarch:=proc() local i,j,indx,rec,n:
global T,d:


while nops(T)<>1 do
	
#Find min cluster distances 

n:=nops(T):
rec:=infinity:	

	for i to n do
   	   for j to i - 1 do
		if d[i,j]<rec then
		 rec:=d[i,j]:
		 indx:=[i,j]:
		fi:
	   od:
	od:

i:=indx[1];j:=indx[2];
d[i,1..n]:=1/2*d[i,1..n]+1/2*d[j,1..n];
d[1..n,i]:=d[i,1..n]^%T;
d[i,i]:=0;
d:=DeleteRow(d,j);d:=DeleteColumn(d,j);
T:=[op(1..j-1,T),op(j+1..i-1,T),[T[i],T[j]],op(i+1..n,T)];
od:
T:
end: