########### ACP e CLUSTERING 
########### dados santarém 
#### 3 aula 


load("/Users/pedrosilva/Documents/ISA25_26/MMA/Aulas/aulas/dadosMulti.rdata")


### variáveis não normalizadas
### 
head(santarem)
round(cov(santarem),2)

D<-dist(santarem) ## matriz de distâncias


### hierarquical agglomerative (ascending) clustering (HAC)

hc.S=hclust(D,method="single")
hc.C=hclust(D,method="complete")
hc.A=hclust(D,method="average")
hc.W=hclust(D,method="ward.D2")


par(mfrow=c(2,2)) ### 4 plots 2x2 
plot(hc.S)
plot(hc.C)
plot(hc.A)
plot(hc.W)

### corte do dendrograma no numéro 
### de clusters sugerido pelo dendrograma 

cls.S=cutree(hc.S,4)
cls.C=cutree(hc.C,3)
cls.A=cutree(hc.A,2)
cls.W=cutree(hc.W,2)

X_ACP<-prcomp(santarem) ### ACP sobre a matriz de covariancias
summary(X_ACP) ### O PFP explica mais de 98% da variância


### single 
plot(X_ACP$x[,1:2],pch=16,col=cls.S,cex=1.5,asp=TRUE,main="Single")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP$x[,1:2],labels=substr(rownames(santarem),1,3),pos=3)
color_biplot(santarem,col.ind =cls.S,
             labels.ind=substr(rownames(santarem),1,3),
             cex.ind=.75,main="single")
#biplot(X_ACP)


plot(X_ACP$x[,1:2],pch=16,col=cls.C,cex=1.5,main="Complete")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP$x[,1:2],labels=substr(rownames(santarem),1,1),pos=3)
color_biplot(santarem,col.ind =cls.C,labels.ind=substr(rownames(santarem),1,3),cex.ind=.75,main="complete")



plot(X_ACP$x[,1:2],pch=16,col=cls.A,cex=1.5,asp=TRUE,main="Average")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP$x[,1:2],labels=substr(rownames(santarem),1,1),pos=3)
color_biplot(santarem,col.ind =cls.A,labels.ind=substr(rownames(santarem),1,3),cex.ind=.75,main="average")


plot(X_ACP$x[,1:2],pch=16,col=cls.W,cex=1.5,asp=TRUE,main="Ward")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP$x[,1:2],labels=substr(rownames(santarem),1,1),pos=3)
color_biplot(santarem,col.ind =cls.W,labels.ind=substr(rownames(santarem),1,3),cex.ind=.75,main="ward")







########### CLUSTERING sobre os dados  normalizadas

santarem.R<-scale(santarem,scale=TRUE) ### variable standardization
D.R=dist(santarem.R)
X_ACP_R<-prcomp(santarem.R) ## <=> prcomp(santarem, scale=TRUE) 
summary(X_ACP_R)
s.R<-summary(X_ACP_R)

## pela regra de KAISER devímaos reter 4 PCs (var >= 1), 
## que explicam mais de 86% da variabilidade total do conjunto
## de dados

par(mfrow=c(1,1))
plot(s.R$sdev^2,type="b",lty=2)
csum<-cumsum(s.R$sdev^2)
plot(csum,type="b",lty=2) ## we have an elbow point at 4 PCs

scores.R<-X_ACP_R$x

hc.S=hclust(D.R,method="single")
hc.C=hclust(D.R,method="complete")
hc.A=hclust(D.R,method="average")
hc.W=hclust(D.R,method="ward.D2")

par(mfrow=c(2,2))
plot(hc.S,h=-1)
plot(hc.C,h=-1)
plot(hc.A,h=-1)
plot(hc.W,h=-1)

cls.S=cutree(hc.S,3)
cls.C=cutree(hc.C,6)
cls.A=cutree(hc.A,3)
cls.W=cutree(hc.W,4)

par(mfrow=c(1,1))
plot(X_ACP_R$x[,1:2],pch=16,col=cls.S,cex=1.5,asp=FALSE,main="Single")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP_R$x[,1:2],labels=substr(rownames(santarem),1,3),pos=3)

plot(X_ACP_R$x[,1:2],pch=16,col=cls.C,cex=1.5,asp=TRUE,main="Complete")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP_R$x[,1:2],labels=substr(rownames(santarem),1,3),pos=3)

plot(X_ACP_R$x[,1:2],pch=16,col=cls.A,cex=1.5,asp=TRUE,main="Average")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP_R$x[,1:2],labels=substr(rownames(santarem),1,3),pos=3)

plot(X_ACP_R$x[,1:2],pch=16,col=cls.W,cex=1.5,asp=TRUE,main="Ward")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP_R$x[,1:2],labels=substr(rownames(santarem),1,3),pos=3)

color_biplot(santarem.R,col.ind =cls.W+2,labels.ind=substr(rownames(santarem),1,3),cex.ind=.75)

###qualidade de representação de Rio Maior 
sum(scores.R["Rio Maior",1:2]^2)/sum(scores.R["Rio Maior",]^2)
##
acos(sqrt(sum(scores.R["Coruche",1:2]^2)/sum(scores.R["Coruche",]^2)))/pi*2





par(mfrow=c(1,1))

### qualidade de representação de Rio Maior 
cos2<-rowSums(scores.R[,1:2]^2)/sum(scores.R^2)

### biplot com qualidade de representação de cada concelho no PFP
color_biplot(santarem.R,col.var=rep("black",ncol(santarem.R)),
             col.ind =cls.W+1,cex.ind=cos2*25+.5,
             labels.ind=substr(rownames(santarem),1,10),
             main="Ward")




#### PARTITIONAL CLUSTERING 

require(NbClust)
require(cluster)
require(fpc)



NbClust(santarem.R, max.nc = 4,method="kmeans")
km<-kmeans(santarem.R,centers=4,nstart=500)

km$totss ## SSQt
km$withinss ##SSQw
km$betweenss ## SSQb

c1<-colMeans(santarem.R[cls.W==1,])
c2<-colMeans(santarem.R[cls.W==2,])
c3<-colMeans(santarem.R[cls.W==3,])
c4<-colMeans(santarem.R[cls.W==4,])

km<-kmeans(santarem.R,centers=rbind(c1,c2,c3,c4),nstart=500)



km$totss ## SSQt
km$withinss ##SSQw
km$betweenss ## SSQb

sum(km$betweenss,km$withinss)

N <- nrow(santarem.R)
(N-1)*sum(diag(var(santarem.R)))

km.cls<-km$cluster

calinhara(santarem.R,km.cls)


SSQw<-rep(0,9)

for (k in 2:10){
  km<-kmeans(santarem.R,centers=k,nstart=500)
  km.cls<-km$cluster
 # print(calinhara(santarem.R,km.cls))
  sil<-silhouette(km.cls,dist(santarem.R)) ### how well is each el in the class
  s<-summary(sil)
# sil.width <- mean(sil[,3])
  SSQw[k]<-sum(km$withinss)
  print(round(c(k,calinhara(santarem.R,km.cls),s$avg.width),3))
}
plot(SSQw,lty=2,type="b",xlim=c(2,10))

### both CH and silhouette propose 6 clusters
### there is no clear elbow point in the scree plot of SSQw

km<-kmeans(santarem.R,centers=6,nstart=500)
km.cls<-km$cluster
sil<-silhouette(km.cls,dist(santarem.R)) ### how well is each el in the class
s<-summary(sil,ordered=TRUE)
rownames(sil)<-rownames(santarem.R)

plot(sil)
summary(sil)

par(mfrow=c(1,1))

plot(X_ACP$x[,1:2],pch=16,col=km.cls,cex=2,main="k-means")
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
text(X_ACP$x[,1:2],labels=substr(rownames(santarem),1,2),pos=3)

color_biplot(santarem.R,col.var=rep("black",ncol(santarem.R)),
             col.ind =km.cls+1,cex.ind=cos2*25+.5,
             labels.ind=substr(rownames(santarem),1,10),
             main="k-means")




###. classificação das variáveis
###  disregarding the correlation sign 

par(mfrow=c(1,2))

D=as.dist(sqrt(1-cor(santarem)^2))


hc.var<-hclust(D,method="complete")


plot(hc.var)



### usando a correlação 
### accounting for the correlatoion sign


D=as.dist((1-cor(santarem))/2)


hc.var<-hclust(D,method="complete")



plot(hc.var)





#############################################################
### An alternative approach to cluster the variables using
### the kendall-tau ranking distance
#############################################################

require(rankdist) 
santarem.rank<-mapply("rank",santarem,2)
santarem.rank


### construçaõ da matriz de distancias
### usando a distancia de kendall-tau
M=matrix(0,ncol=9,nrow=9)
for (i in 1:9){
  for (j in 1:9){
    r1<-as.vector(santarem.rank[,i])
    r2<-as.vector(santarem.rank[,j])
    M[i,j]<-DistancePair(r1,r2) # computes Kendall distance between the pair of rankings
  }
}

rownames(M)<-colnames<-colnames(santarem)

M[1:5,1:5]
D.k<-as.dist(M)
D.k

par(mfrow=c(1,2))
hc.k<-hclust(D.k)
plot(hc.k)








#### Comparação das diferentes classificações hierárquicas e k-means
#### usando uma diss D definida a partir do índice de RAND: D = 1-RAND

#### TPC: adaptar para o indice de RAND ajustado (ARI): D = (1-ARI)/2
#### ver slides... 

rand <- function(class1,class2){
  n <- length(class1)
  c <- as.dist(outer(class1,class1,"=="))
  d <- as.dist(outer(class2,class2,"=="))
  rand <- sum(c == d)/(n*(n-1)/2)
  return(rand) }


rand(c(1,1,2,3,1,3),c(1,2,1,2,3,3)) # 0.5333333
# 2 random samples of length 1000 with elements extracted from 1,...,10
p1<-sample(1:10,1000,replace=TRUE)
p2<-sample(1:10,1000,replace=TRUE)
rand(p1,p2) # 0.8196997             

cls.S=cutree(hc.S,4)
cls.C=cutree(hc.C,4)
cls.A=cutree(hc.A,4)
cls.W=cutree(hc.W,4)


classif<-cbind(cls.S,cls.C,cls.A,cls.W,km.cls)

M.R=matrix(0,ncol=5,nrow=5)
for (i in 1:5){
  for (j in 1:5){
    r1<-as.vector(classif[,i])
    r2<-as.vector(classif[,j])    
    M.R[i,j]<-rand(r1,r2) # computes Kendall distance between the pair of rankings
  }
}

M.R
rownames(M.R)<-colnames(M.R)<-c("S","C","A","W","K")


D.R <- as.dist(1-M.R) ### converte a semelhanca em diss.

D.R

hc.R<-hclust(D.R)
plot(hc.R)



################################################################