if(!require("TDA")) {
    install.packages("TDA")
}

if(!require("scatterplot3d")) {
    install.packaged("scatterplot3d")
}

library("TDA")
library("scatterplot3d")

#~ This is for a Mac... might be something else for Linux/Windows
setwd("~/Desktop/")

#~ Get the bone data; exclude first two rows; jump to 4D coordinates (polygon information) after row 5130, so fill will resolve any array alignment problems
bone1 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/001.off",sep = " ", skip = 2,fill = TRUE)
bone2 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/002.off",sep = " ", skip = 2,fill = TRUE)
bone3 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/003.off",sep = " ", skip = 2,fill = TRUE)
bone4 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/004.off",sep = " ", skip = 2,fill = TRUE)
bone5 = read.table("http://www2.stat.duke.edu/~sayan/auto3dgm/Bones/Aligned_Shapes/005.off",sep = " ", skip = 2,fill = TRUE)


#~ Extract the coordinate information
bone1 = bone1[1:5128,]
bone2 = bone2[1:5086,]
bone3 = bone3[1:5091,]
bone4 = bone4[1:5092,]
bone5 = bone5[1:5154,]

#~ Rename the columns
colnames(bone1) = c("Z","Y","Z")
colnames(bone2) = c("Z","Y","Z")
colnames(bone3) = c("Z","Y","Z")
colnames(bone4) = c("Z","Y","Z")
colnames(bone5) = c("Z","Y","Z")

#~ Loop through all datasets
# Programming Note: parse() converts text to a command object; eval() evaluates a given command object
for(i in 1:5) {
    #~ Plot the Data; save plot as a pdf
    eval(parse(text = paste0("pdf('bone_",i,"_scatterplot.pdf'); scatterplot3d(bone",i,",main='3D Bone ",i," Point Cloud');dev.off()")))
    #~ Compute Homology using VR Complex
    eval(parse(text = paste0("boneRips.",i," = ripsDiag(Z = bone",i,",maxdimension = 1, maxscale = 0.015)")))
    #~ Plot Persistence Diagrams and Barcodes; save plots as a pdf
    eval(parse(text = paste0("pdf('bone_rips_",i,".pdf');par(mfrow = c(2,1));plot(boneRips.",i,"[['diagram']],main = 'Bone ",i," Persistence Diagram');plot(boneRips.",i,"[['diagram']],barcode = TRUE,main = 'Bone ",i," Barcode Diagram');dev.off()")))
}

#~ Evaluate over a grid
# Some constants used for evaluation
by = 0.01
buffer = 0.05
seqlen = 10
band = 0.01

# Limit of x coordinates
Xlim = max(bone1[,"X"],bone2[,"X"],bone3[,"X"],bone4[,"X"],bone5[,"X"]) + buffer
Xlim = c(min(bone1[,"X"],bone2[,"X"],bone3[,"X"],bone4[,"X"],bone5[,"X"]) - buffer,Xlim)
Xgrid = seq(Xlim[1],Xlim[2],by = by)
# y
Ylim = max(bone1[,"Y"],bone2[,"Y"],bone3[,"Y"],bone4[,"Y"],bone5[,"Y"]) + buffer
Ylim = c(min(bone1[,"Y"],bone2[,"Y"],bone3[,"Y"],bone4[,"Y"],bone5[,"Y"]) - buffer,Ylim)
Ygrid = seq(Ylim[1],Ylim[2],by = by)
# Z
Zlim = max(bone1[,"Z"],bone2[,"Z"],bone3[,"Z"],bone4[,"Z"],bone5[,"Z"]) + buffer
Zlim = c(min(bone1[,"Z"],bone2[,"Z"],bone3[,"Z"],bone4[,"Z"],bone5[,"Z"]) - buffer,Zlim)
Zgrid = seq(Zlim[1],Zlim[2],by = by)

# construct grid/mesh over the generated sequences
grid = expand.grid(Xgrid,Ygrid,Zgrid)

#~ Compute homology over a grid
for (i in 1:5) {
    #~ Compute bootstrap bandwidth
    eval(parse(text = paste0("bootband = bootstrapBand(X = bone",i,", FUN = kde, Grid = grid, B = 10, parallel = TRUE, alpha=0.05, h=",band,")")))
    #~ Compute Homology over a grid
    eval(parse(text = paste0("boneGrid.",i," = gridDiag(X = bone",i,", FUN = kde, h = ",band,", lim = cbind(Xlim,Ylim,Zlim), by = ",by,", sublevel = FALSE, library = 'Dionysus')")))
    #~ Plot Persistence Diagrams and Barcodes; save plots as a pdf
    eval(parse(text = paste0("pdf('bone_grid_",i,".pdf');par(mfrow = c(2,1));plot(boneGrid.",i,"[['diagram']], band = ",2*bootband$width,",main = 'Bone ",i," Persistence Diagram');plot(boneGrid.",i,"[['diagram']],barcode = TRUE,main = 'Bone ",i," Barcode Diagram');dev.off()")))
}