Librerias requeridas

setwd("C:/Users/Alumnos/Desktop")

#fastcluster requerido para TDAmapper
install.packages("fastcluster")

#TDAmapper para mapper1D, mapper2D
install.packages("TDAmapper") 

#TDA para homologia persistente y sampleo uniforme circulo,esfera, toro
install.packages("TDA")

#igraph para mostrar grafos
install.packages("igraph")                  

#Para visualizacion 3D
install.packages("ggplot2")
install.packages("rgl")

#Para grabar grafos en formato GEFX
install.packages("rgexf")
install.packages("plyr")

#paquete para hacer PCA
install.packages("FactoMineR")

#paquete para algoritmo TDAmapper
require(fastcluster) 
library("TDAmapper")

#paquete para pers.hom, muestreo S1,S2,S1xS1
library("TDA")

#para manejo de grafos (y grabar GEXF)
library("igraph")
library("rgexf")
library("plyr")

#para visualizaciones de PCA y 3D
library("ggplot2")
library("rgl")

#para PCA
library("FactoMineR")

m_XXX <- mapper1D(
  distance_matrix = dist(d_XXX),
  filter_values =filtro,
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

dist() calcula la distancia (euclideana) de un conjunto de puntos en R^n
filtro es un vector de valores (funcion en el conjunto de puntos)

mapper1d retorna una lista de 6 con:

adjacency : matriz de adjacency del grafo
num_vertices : numero de vertices
level_of_vertex : rango de la funcion en cada vertice
points_in_vertex : datos en cada vertices()
points_in_level : datos en cada rango
vertices_in_level : vertices en cada rango

mapper1d retorna una lista de 6 con:

adjacency : matriz de adjacency del grafo
num_vertices : numero de vertices
level_of_vertex : rango de la funcion en cada vertice
points_in_vertex : datos en cada vertices()
points_in_level : datos en cada rango
vertices_in_level : vertices en cada rango

numpoints<-1000
epsilon<-0.1
r<-1

#generar datos uniformes en el circulo

datos_circulo_uniforme<-as.data.frame(circleUnif(numpoints))

mapper_circulo_uniforme <- mapper1D(
  distance_matrix = dist(datos_circulo_uniforme),
  filter_values =datos_circulo_uniforme[,1],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

grafo_circulo_uniforme <- graph.adjacency(mapper_circulo_uniforme$adjacency, mode="undirected")

#generar datos circulo
numpoints<-100
t<-runif(numpoints)

datos_circulo<-data.frame(x=cos(2*pi*t),y=sin(2*pi*t))

mapper_circulo <- mapper1D(
  distance_matrix = dist(datos_circulo),
  filter_values =datos_circulo[,1],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

grafo_circulo <- graph.adjacency(mapper_circulo$adjacency, mode="undirected")

numpoints<-1000
#generar datos circulo con ruido
datos_circulo_con_ruido<-datos_circulo_uniforme+epsilon*matrix(rnorm(2*numpoints,mean=0,sd=0.1),ncol=2)

mapper_circulo_con_ruido <- mapper1D(
  distance_matrix = dist(datos_circulo_con_ruido),
  filter_values =datos_circulo_con_ruido[,1],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

grafo_circulo_con_ruido <- graph.adjacency(mapper_circulo_con_ruido$adjacency, mode="undirected")

numpoints<-100
datos_ocho=data.frame( x=2*r*cos(0.5*(1:numpoints)), y=r*sin(1:numpoints) )

mapper_ocho<- mapper1D(
  distance_matrix = dist(datos_ocho),
  filter_values =datos_ocho[,1],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)



grafo_ocho <- graph.adjacency(mapper_ocho$adjacency, mode="undirected")

#elipse
datos_elipse=data.frame( x=2*r*cos(1:numpoints), y=sin(1:numpoints) )

mapper_elipse<- mapper1D(
  distance_matrix = dist(datos_elipse),
  filter_values =datos_elipse[,1],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)




grafo_elipse <- graph.adjacency(mapper_elipse$adjacency, mode="undirected")

#dos espirales
set.seed("1")
t <- runif(numpoints, min=1, max=6.3) # theta

datos_espirales <- data.frame( x = c( t*cos(t), -t*cos(t) ), y = c( t*sin(t), -t*sin(t) ) )

mapper_espirales<- mapper1D(
  distance_matrix = dist(datos_espirales),
  filter_values =datos_espirales[,1],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)



grafo_espirales <- graph.adjacency(mapper_espirales$adjacency, mode="undirected")

#transforma datos en R^2
theta<-60/180*pi

Rotation_Matrix<-matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), nrow=2, ncol=2)

Otra_Matrix <- matrix(c(1, 0, 0, 1), nrow=2, ncol=2) 

transformacion_elipse <- (Rotation_Matrix %*% ( Otra_Matrix   %*%  t(datos_elipse)))   

datos_transformados_elipse<-as.data.frame(t(transformacion_elipse))

mapper_elipse_transformada<- mapper1D(
  distance_matrix = dist(datos_transformados_elipse),
  filter_values =datos_transformados_elipse[,1],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)


grafo_elipse_transformada <- graph.adjacency(mapper_elipse_transformada$adjacency, mode="undirected")

maxdimension <- 1
maxscale <- 5
Diag <- ripsDiag(X = datos_transformados_elipse, maxdimension, maxscale,library = "GUDHI", printProgress = FALSE)

#generar datos curva de momento trigonometrica en R4
datos_curva_momento<-data.frame(x1=cos((1:numpoints)),y1=sin((1:numpoints)),x2=cos(2*(1:numpoints)),y2=sin(2*(1:numpoints)))

X<-(datos_curva_momento$x1+2)*datos_curva_momento$x2
Y<-(datos_curva_momento$x1+2)*datos_curva_momento$y2
Z<-datos_curva_momento$y1

plot3d(X,Y,Z)

mapper_curva_momento<- mapper1D(
  distance_matrix = dist(datos_curva_momento),
  filter_values =X,
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

grafo_curva_momento <- graph.adjacency(mapper_curva_momento$adjacency, mode="undirected")

#generar datos curva de momento trigonometrica en R4 de nuevo
numpoints<-1000

datos_curva_momento<-data.frame(x1=cos((1:numpoints)),y1=sin((1:numpoints)),x2=cos(2*(1:numpoints)),y2=sin(2*(1:numpoints)))

X<-(datos_curva_momento$x1+2)*datos_curva_momento$x2
Y<-(datos_curva_momento$x1+2)*datos_curva_momento$y2
Z<-datos_curva_momento$y1
plot3d(X,Y,Z)

mapper_curva_momento<- mapper1D(
  distance_matrix = dist(datos_curva_momento),
  filter_values =X,
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

grafo_curva_momento <- graph.adjacency(mapper_curva_momento$adjacency, mode="undirected")

#generar datos curva de momento trigonometrica en R6
datos_curva_momento2<-data.frame(x1=cos((1:numpoints)),y1=sin((1:numpoints)),x2=cos(2*(1:numpoints)),y2=sin(2*(1:numpoints)),x3=cos(3*(1:numpoints)),y3=sin(3*(1:numpoints)))

mapper_curva_momento2<- mapper1D(
  distance_matrix = dist(datos_curva_momento2),
  filter_values =X,
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

mapper_XXX<-mapper2D(
  distance_matrix=dist(datos_XXX), 
  filter_values = list(filtro1,filtro2),    
  c(num_intervals,num_intervals), # use default
  percent_overlap = 50,
  num_bins_when_clustering = 10
)

numpoints<-1000
num_intervals = 10
percent_overlap = 50
num_bins_when_clustering = 10

r<-1

#esfera en R^3 uniforme

datos_esfera_uniforme=as.data.frame(sphereUnif(numpoints, 2, r))

plot3d(datos_esfera_uniforme$V1, datos_esfera_uniforme$V2, datos_esfera_uniforme$V3)

mapper_esfera_uniforme <- mapper1D(
  distance_matrix = dist(datos_esfera_uniforme),
  filter_values =datos_esfera_uniforme[,3],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10
  )

grafo_esfera_uniforme <- graph.adjacency(mapper_esfera_uniforme$adjacency, mode="undirected")

mapper_esfera2d<-mapper2D(
  distance_matrix=dist(datos_esfera_uniforme), 
  filter_values = list(datos_esfera_uniforme[,1],datos_esfera_uniforme[,2]),    
  c(num_intervals,num_intervals), # use default
  percent_overlap = 50,
  num_bins_when_clustering = 10
)

grafo_esfera2d <- graph.adjacency(mapper_esfera2d$adjacency, mode="undirected")

#necesita librerias rxgef,plyr
saveAsGEXF(grafo_esfera2d)

#Miremos como luce los datos del grafo con MDS
distMatrixgrafoEsfera2D <- shortest.paths(grafo_esfera2d, v=V(grafo_esfera2d), to=V(grafo_esfera2d))
fit<-cmdscale(distMatrixgrafoEsfera2D,3)
plot3d(fit)

#Homologia persistente del grafo???
#maxdimension <- 1
#maxscale <- 5
#Diag <- ripsDiag(distMatrixgrafoEsfera2D, maxdimension, maxscale,dist="arbitrary",library = "GUDHI", printProgress = FALSE)
#plot(Diag[["diagram"]],barcode=TRUE)

#persistence homology de la esfera se queda...
#maxdimension <- 1
#maxscale <- 2
#Diag <- ripsDiag(X = datos_esfera_uniforme, maxdimension, maxscale,library = "GUDHI", printProgress = FALSE)
#plot(Diag[["diagram"]],barcode=TRUE)

#toro uniforme
R<-2

datos_toro_uniforme=as.data.frame(torusUnif(numpoints, r, R))

plot3d(datos_toro_uniforme$x, datos_toro_uniforme$y, datos_toro_uniforme$z)

mapper_toro_uniforme <- mapper1D(
  distance_matrix = dist(datos_toro_uniforme),
  filter_values =datos_toro_uniforme[,3],
  num_intervals = 10,
  percent_overlap = 50,
  num_bins_when_clustering = 10)

grafo_toro_uniforme <- graph.adjacency(mapper_toro_uniforme$adjacency, mode="undirected")

mapper_toro_uniforme_2d <- mapper2D(
  distance_matrix=dist(datos_toro_uniforme),
  filter_values = list(datos_toro_uniforme[,1],datos_toro_uniforme[,2]),    
  num_intervals =c(num_intervals,num_intervals), # use default
  percent_overlap=50, # use default
  num_bins_when_clustering=10  # use default
)

[1] "Level set has only one point"

grafo_toro_uniforme2d <- graph.adjacency(mapper_toro_uniforme_2d$adjacency, mode="undirected")

#necesita librerias rxgef,plyr
saveAsGEXF(grafo_toro_uniforme2d)

#Miremos como luce los datos del grafo con MDS
distMatrix_grafo_torouniforme2D <- shortest.paths(grafo_toro_uniforme2d, v=V(grafo_toro_uniforme2d), to=V(grafo_toro_uniforme2d))
fit<-cmdscale(distMatrix_grafo_torouniforme2D,3)
plot3d(fit)

#Homologia persistente del grafo???
#maxdimension <- 1
#maxscale <- 5
#Diag <- ripsDiag(distMatrix_grafo_torouniforme2D, maxdimension, #maxscale,dist="arbitrary",library = "GUDHI", printProgress = FALSE)
#plot(Diag[["diagram"]],barcode=TRUE)


#Homologia persitente se queda
#maxdimension <- 1
#maxscale <- 5
#Diag <- ripsDiag(X = datos_toro_uniforme, maxdimension, maxscale,library = "GUDHI", printProgress = FALSE)
#plot(Diag[["diagram"]],barcode=TRUE)

#puntos en el toro plano
set.seed("1")
U <- runif(numpoints, min=0, max=1) 
V <- runif(numpoints, min=0, max=1) 

X<-(r*cos(2*pi*U)+R)*cos(2*pi*V)
Y<-(r*cos(2*pi*U)+R)*sin(2*pi*V)
Z<-r*sin(2*pi*U)

datos_toro_plano=as.data.frame(cbind(X,Y,Z))

plot3d(X,Y,Z)

#distancia toro plano
distance_matrix1<-matrix(0,nrow=numpoints,ncol=numpoints)
for (i in 1:numpoints){
  for (j in i:numpoints){
    distance_matrix1[i,j]<-sqrt((min(abs(U[i]-U[j]),1-abs(U[i]-U[j])))^2+(min(abs(V[i]-V[j]),1-abs(V[i]-V[j])))^2)
    distance_matrix1[j,i]<-distance_matrix1[i,j]
  }
}

fit<-cmdscale(distance_matrix1,3)
plot3d(fit)

mapper_toro_plano <- mapper1D(
  distance_matrix=dist(distance_matrix1),
  filter_values = cos(U),   
  num_intervals =10, # use default
  percent_overlap=50, # use default
  num_bins_when_clustering=10  # use default
)

grafo_toro_plano <- graph.adjacency(mapper_toro_plano$adjacency, mode="undirected")

mapper_toro_plano2d <- mapper2D(
  distance_matrix=dist(distance_matrix1),
  filter_values = list(cos(U),sin(U)),  
  num_intervals = c(10,10), # use default
  percent_overlap=50, # use default
  num_bins_when_clustering=10  # use default
)

[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set has only one point"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"
[1] "Level set is empty"

grafo_toro_plano2d <- graph.adjacency(mapper_toro_plano2d$adjacency, mode="undirected")

install.packages("locfit")
library(locfit)

#OJO!!!!!! el paquete ks y el paquete TDA tienen ambos funciones llamadas kde con distintos parametros!!!!
install.packages("ks")
library(ks)

data(chemdiab)

summary(chemdiab)

normdiab <- chemdiab
normdiab[,1:5] <- scale(normdiab[,1:5],center=FALSE)
normdiab.dist = dist(normdiab[,1:5])

filter.kde <- ks::kde(normdiab[,1:5],H=diag(1,nrow = 5),eval.points = normdiab[,1:5])$estimate

diabetes_pca <- PCA(normdiab,graph = FALSE)
qplot(diabetes_pca$ind$coord[,1],diabetes_pca$ind$coord[,2])
plot3d(diabetes_pca$ind$coord[,1],diabetes_pca$ind$coord[,2],diabetes_pca$ind$coord[,3])

diab.mapper <- mapper1D(
  normdiab.dist,
  filter.kde,
  num_intervals = 4,
  percent_overlap = 50,
  num_bins_when_clustering = 20)

diab.graph <- graph.adjacency(diab.mapper$adjacency, mode="undirected")

plot(diab.graph )