################ Script 1 (les scripts sont enchaines les uns sous les autres) ########################
# David NERINI - ACP - Sept 2008
# Ce script permet de générer un nuage de points en 3D,
# de s'y ballader de manière dynamique et
# d'illustrer les notions d'inertie (totale T, intra W et inter B)
# dans un nuage de points. Il effectue egalement les calculs de ces inerties
# et permet de verifier la relation T=B+W.
# Il nécessite le package rgl (OpenGL)
library(rgl)
nx=70
ny=150
nz=100
n=sum(nx,ny,nz)

#On génère un nuage de points consistué de 3 groupes, selon une loi normale
# de paramètres mu et sigme donné. le nuage est représenté en 3D
x=cbind(rnorm(nx,0,1),rnorm(nx,20,2),rnorm(nx,0,2))
y=cbind(rnorm(ny,10,4),rnorm(ny,12,3),rnorm(ny,10,2))
z=cbind(rnorm(nz,-5,1),rnorm(nz,-4,1),rnorm(nz,-3,1))
col=c(rep(1,nx),rep(2,ny),rep(4,nz))
open3d()
spheres3d(rbind(x,y,z),radius=0.2,col=col,alpha=0.8)
rgl.select3d(button="right")

#On représente les centres de gravités des trois classes ...
gx=apply(x,2,mean)
gy=apply(y,2,mean)
gz=apply(z,2,mean)

#... et le centre de gravité total
g=apply(rbind(x,y,z),2,mean)
spheres3d(g,col=3)
rgl.select3d(button="right")


#représente les éléments utilisés pour calculer l'inertie totale
# (chaque point est relié à g)
M=rbind(x,y,z)
Mg=matrix(rep(g,n),n,3,byrow=T)
Mt=matrix(c(t(cbind(M,Mg))),2*n,3,byrow=T)

segments3d(Mt,col=3,alpha=0.5)

M=M-Mg
T=round(1/n*sum(diag(M%*%t(M))),1)
text3d(g[1]+2,g[2],g[3],texts=paste("T=",T))
rgl.select3d(button="right")

#Représente les éléments utilisés pour le calcul de l'inertie inter
# Chaque centre de gravité local est relié à g)

spheres3d(rbind(gx,gy,gz),col=c(1,2,4),radius=0.6)
rgl.select3d(button="right")
segments3d(rbind(g,gx,g,gy,g,gz),col=c(1,1,2,2,4,4),size=10)
rgl.select3d(button="right")
B=rbind(gx,gy,gz)-rbind(g,g,g)
B=round(1/n*sum(diag(B%*%t(B))*c(nx,ny,nz)),1)
text3d(g[1],g[2]+2,g[3],texts=paste("B=",B))
rgl.select3d(button="right")


#Représente les éléments utilisé pour le calcul de l'inertie intra
# Chaque point d'une classe est relié à son centre de gravité local
Mgx=matrix(rep(gx,nx),nx,3,byrow=T)
Mtx=matrix(c(t(cbind(x,Mgx))),2*nx,3,byrow=T)
segments3d(Mtx,col=1,alpha=0.5)

Mgy=matrix(rep(gy,ny),ny,3,byrow=T)
Mty=matrix(c(t(cbind(y,Mgy))),2*ny,3,byrow=T)
segments3d(Mty,col=2,alpha=0.5)

Mgz=matrix(rep(gz,nz),nz,3,byrow=T)
Mtz=matrix(c(t(cbind(z,Mgz))),2*nz,3,byrow=T)
segments3d(Mtz,col=4,alpha=0.5)

Mx=x-Mgx
My=y-Mgy
Mz=z-Mgz
Wx=sum(diag(Mx%*%t(Mx))*1/n)
Wy=sum(diag(My%*%t(My))*1/n)
Wz=sum(diag(Mz%*%t(Mz))*1/n)
W=round(sum(Wx,Wy,Wz),1)
text3d(g[1],g[2],g[3]+2,texts=paste("W=",W))
rgl.select3d(button="right")
text3d(g+10,texts=paste("T=W+B"))





###################### Script 2 #################################
#David NERINI - ACP - Sept 2008
#Ce script permet de montrer les différentes étapes de l'algorithme des nuées dynamiques.
# Il est réalisé à partir d'un nuage de points dans le plan, constitué de 2 classes.
# On cherche, sans information a priori, à retrouver ces classes.

#Creation et représentation du nuage de points (deux gaussiennes)
xA=cbind(rnorm(20),rnorm(20))
xB=cbind(rnorm(10,4,1),rnorm(10,4,1))
x=rbind(xA,xB)
names(x)=paste("p",1:30,sep="")

# Calcul de la matrice d'interdistance entre objets
# (métrique euclidienne usuelle)
#N'hesitez pas a regarder la structure de cet objet

D=as.matrix(dist(x))

#Premiere etape : on tire k représentants des classes (k=2 ici)
#soit aleatoirement avec ind=sample(1:30,2,replace=F)
# soit, dans cet exemple, les points 3 et 5

ind=c(3,5)
n=10
for(k in 1:n){
#dessin du nuage de points
   plot(x,type="n",xlab="",ylab="", main=paste("step ",k))
   text(x,names(x),col=c(rep(3,20),rep(1,10)))
   text(locator(1),"")

   dif=D[ind[1],]-D[ind[2],]
   C1=which(dif<=0)
   C2=which(dif>0)
   C1lig=c(rbind(ind[1],C1))
   C2lig=c(rbind(ind[2],C2))

#etape 1 : dessin des deux représentants tires au hasard,
# association des points du nuage au plus proche représentant
   symbols(x[ind,],circles=c(0.25,0.25),add=T,inches=F)
   text(locator(1),"")
  
   lines(x[C1lig,])
   lines(x[C2lig,])
   text(locator(1),"")

# représentation des centres de gravite
# des classes ainsi definies

   g1=t(as.matrix(apply(x[C1,],2,mean)))
   g2=t(as.matrix(apply(x[C2,],2,mean)))

   points(g1,col="blue")
   points(g2,col="blue")

   text(locator(1),"")
# Etape 2 : Election des nveaux représentants (1er parangon)
   D1=as.matrix(dist(rbind(g1,x[C1,])))[1,-1]
   D2=as.matrix(dist(rbind(g2,x[C2,])))[1,-1]
   ind[1]=C1[which.min(D1)]
   ind[2]=C2[which.min(D2)]
  
   points(x[ind,],col="blue")
   lines(rbind(g1,x[ind[1],]),col="blue")
   lines(rbind(g2,x[ind[2],]),col="blue")
   text(locator(1),"")
  
#les etapes sont reprises n fois ou jusqu a la convergence d'un critere que vous pouvez
#programmer
}