library("TDA") X = read.table("crater.txt") # read data plot(X[,1], X[,2], "p") # grid Xlim <- c(min(X[,1]), max(X[,1])); Ylim <- c(min(X[,2]), max(X[,2])) byx <- (Xlim[2]-Xlim[1])/200; byy <- (Ylim[2]-Ylim[1])/200 Xseq <- seq(Xlim[1], Xlim[2], by = byx); Yseq <- seq(Ylim[1], Ylim[2], by = byy) Grid <- expand.grid(Xseq, Yseq) # kNN density estimator k=200 kNN <- knnDE(X = X, Grid = Grid, k = k) # plot dev.new() #persp(Xseq, Yseq, matrix(kNN, ncol = length(Yseq), nrow = length(Xseq))) # [optional] nice plot persp(Xseq, Yseq, matrix(kNN, ncol = length(Yseq), nrow = length(Xseq)),ylab = "", zlab = "", theta = -70, phi = 65, ltheta = 70,col = 2, border = NA, main = "kNN", expand = 1, shade = 0.9) # compute persistence diagram Diag <- gridDiag(X = X, FUNvalues = matrix(kNN, ncol = length(Yseq), nrow = length(Xseq)), lim = cbind(Xlim, Ylim), by = c(byx, byy), sublevel = FALSE, library = "Dionysus", printProgress = TRUE) # [optional] compute persistence diagram (compute grid and estimator on the fly) #Diag <- gridDiag(X = X, FUN = knnDE, k = 200, lim = cbind(Xlim, Ylim),by = by, sublevel = FALSE, library = "Dionysus", printProgress = TRUE) # plot diagram dev.new() plot(Diag[["diagram"]], main = "kNN Diagram")