### R code from vignette source 'jss741.Rnw' ################################################### ### Create 4 peaks data and scale to [0,1] ################################################### library("msr") data("fourpeaks") d <- fourpeaks(2000) d[ , 1] <- ( d[ , 1] - min(d[ , 1]) ) / (max(d[ , 1]) - min(d[ , 1])) ################################################### ### Create Morse-Smale complex ################################################### ms <- msc.nn(y = d[ , 1], x = d[ , 2:3], pLevel = 0.1, knn = 15) summary(ms) ################################################### ### plot surface ################################################### par(mfcol = c(1, 2)) fancy <- c('#0066CC', '#CCCC00', '#D22905') colormap = colorRamp(fancy, interpolate = "linear", bias = 0.5 ) colors <- rgb(colormap(d[ , 1]), maxColorValue = 255) par(mar = c(8, 5, 5, 4)) plot(d[ , 2], d[ , 3], type = "p", xlab = "", ylab = "", col = colors, pch = 19, cex.axis = 2.5 ) pal <- brewer.pal(9, "Set1") pal <- colorRampPalette(pal) pal <- pal(length(ms$level[[1]]$mins)) colors <- pal[ms$level[[1]]$partition] par(mar=c(8,5,5,4)) plot(d[ , 2], d[ , 3], col = colors, type = "p", xlab = "", ylab = "", pch = 19, cex.axis = 2.5) ################################################### ### plot persistencies ################################################### np <- length(ms$persistence) ms$persistence[np] <- 1 par(mar = c(6, 4, 4, 4) + 3) plot(ms$persistence, np:1 + 1, xlab = "Persistence Percentage", ylab = "Extrema", cex = 2, cex.lab = 2, cex.axis = 2, type = "p", pch = 19) ################################################### ### create Morse-Smale comples with 15 persistence levels ################################################### ms <- msc.nn(y = d[ , 1], x=d[ , 2:3], nLevels = 15, knn = 15) ################################################### ### summary print the Morse-Smale complex at level 5 ################################################### level5 <- ms$level[[5]] summary(level5) ################################################### ### Create predictive Morse-Smale complex with SVM ################################################### ms.svm <- msc.nn.svm(y = d[ , 1], x = d[ , 2:3], nLevels = 15, knn = 15, cost = 1, precompute = TRUE) ################################################### ### Create predictive Morse-Smale complex with kernel density ################################################### ms.kd1 <- msc.nn.kd(y = d[ , 1], x = d[ , 2:3], nLevels=15, knn = 15, bw = 0.05) ms.kd2 <- msc.nn.kd(y = d[ , 1], x = d[ , 2:3], nLevels=15, knn = 15, bw = 0.15) ################################################### ### predict partition on test data ################################################### test <- fourpeaks(2000) ms.svm$predictLevel <- 12 psvm <- predict(ms.svm, newdata = test) ms.kd1$predictLevel <- 12 ms.kd2$predictLevel <- 12 pkd1 <- predict(ms.kd1, newdata = test) pkd2 <- predict(ms.kd2, newdata = test) ################################################### ### plot predicted partitions ################################################### par(mfcol = c(1, 4)) pId <- ms.svm$level[[12]]$partition == 8 colors <- vector(length = length(ms.svm$y)) colors[] = fancy[1] colors[pId] = fancy[3] par(mar=c(6, 4, 4, 4) + 1) plot(ms$x, col = colors, type = "p", xlab = "", ylab = "", pch = 19, cex.axis = 2, cex = 1 ) title("Ground Truth", cex.main = 3.5) colors <- rgb(colormap(psvm[ , 8]), maxColorValue = 255) par(mar=c(6,4,4,4) + 1) plot(test[ , 2], test[ , 3], type = "p", xlab = "", ylab = "", col = colors, pch = 19, cex.axis = 2, cex = 1 ) title("SVM", cex.main = 3.5) colors <- rgb(colormap(pkd1[ , 8]), maxColorValue = 255) par(mar = c(6, 4, 4, 4) + 1) plot(test[ , 2], test[ , 3], type = "p", xlab = "", ylab = "", col = colors, pch = 19, cex.axis = 2, cex = 1 ) title("KD bw = 0.05", cex.main = 3.5) colors <- rgb(colormap(pkd2[ , 8]), maxColorValue = 255) par(mar = c(6, 4, 4, 4) + 1) plot(test[ , 2], test[ , 3], type = "p", xlab = "", ylab = "", col = colors, pch = 19, cex.axis = 2, cex = 1 ) title("KD bw = 0.15", cex.main = 3.5) ################################################### ### Visualize Morse-Smale complex at level 7 and 12 ################################################### ms$predictLevel <- 2 plot(ms, span = 0.9) ms$predictLevel <- 7 plot(ms) ms$predictLevel <- 12 plot(ms) ################################################### ### Linear model for regression at level 1 ################################################### ms.svm$predictLevel <- 1 msr.lm <- msc.lm(ms.svm) summary(msr.lm) pLevel <- msr.lm$ms$predictLevel print(pLevel) print(msr.lm$lms[[pLevel]]$cv) ################################################### ### Print mean squared residual ################################################### p <- predict(msr.lm, newdata=test) print(mean((p - test[ , 1])^2)) ################################################### ### Weighted linear model ################################################### msr.slm <- msc.slm(ms.svm) print(msr.slm$slm[[pLevel]]$cv) ################################################### ### Weighted linear model a t different persitence levels ################################################### msr.slm$ms$predictLevel = 1 p1 <- predict(msr.slm, newdata=test) print(mean((p1 - test[ , 1])^2)) msr.slm$ms$predictLevel = 7 p2 <- predict(msr.slm, newdata=test) print(mean((p2 - test[ , 1])^2)) msr.slm$ms$predictLevel = 12 p3 <- predict(msr.slm, newdata=test) print(mean((p3 - test[ , 1])^2)) ################################################### ### plot predicted regression surfaces ################################################### par(mfcol = c(1, 3)) ps <- p1 ps[ps > 1] <- 1 ps[ps < 0] <- 0 colors <- rgb(colormap(ps), maxColorValue = 255) plot( test[ ,2], test[ ,3], type = "p", xlab = "", ylab = "", col = colors, pch = 19, cex.axis = 2.5, cex = 1 ) ps <- p2 ps[ps > 1] <- 1 ps[ps < 0] <- 0 colors <- rgb(colormap(ps), maxColorValue = 255) plot( test[ ,2], test[ ,3], type = "p", xlab = "", ylab = "", col = colors, pch = 19, cex.axis = 2.5, cex = 1 ) ps <- p3 ps[ps > 1] <- 1 ps[ps < 0] <- 0 colors <- rgb(colormap(ps), maxColorValue=255) plot( test[ ,2], test[ ,3], type = "p", xlab = "", ylab = "", col = colors, pch = 19, cex.axis = 2.5, cex = 1 ) ################################################### ### load UCI crime data set ################################################### data("uci_crime_subset") ################################################### ### Build Morse-Samle complex ################################################### ms <- msc.nn(y = crimes[ , 100], x = crimes[ , 1:99], knn = 120, nLevels = 5) np <- length(ms$persistence) ms$persistence[np] <- 1 par(mar=c(5,4,4,2) + 3) plot(ms$persistence, np:1 + 1, xlab = "Persistence Percentage", ylab = "Extrema", cex = 2.5, cex.lab = 2, cex.axis = 2, type = "p", pch = 19) ################################################### ### Show Morse-Smale visualization ################################################### ms$predictLevel = 3 obj <- plot(ms) ################################################### ### Selcetd subset for visualization ################################################### obj$plotList <- c(1, 12, 27, 20, 32, 33, 97) names <- colnames(crimes)[obj$plotList] print(names) plot(obj) ################################################### ### Load energy data ################################################### data("camera_estimation") energy$E <- (energy$E - min(energy$E)) / (max(energy$E) - min(energy$E)) train <- energy[1:2000, ] test <- energy[8001:10000, ] ################################################### ### Kernel density estimator Morse-Smale complex ################################################### ms <- msc.nn.kd(y = train$E, x = train[ , 2:10], knn = 40, nLevels = 17, bw = 0.3) ################################################### ### Linaer model based on kernel density estimation Morse-Smale complex ################################################### ms.lm <- msc.lm(ms) ################################################### ### show automatically selected prediction persistence level ################################################### pLevel <- ms.lm$ms$predictLevel print(pLevel) ################################################### ### plot persitence levels ################################################### cv <- c() for(i in 1:ms.lm$ms$nLevels){ cv[i] <- ms.lm$lms[[i]]$cv} par(mar = c(5, 4, 4, 2) + 3) plot(1:ms.lm$ms$nLevels, cv, xlab = "Persistence Level", ylab = "CV error", pch = 19, cex = 2.5, cex.lab = 2, cex.axis = 2) ################################################### ### ANOVA decomposition based on Morse-Smale decomposition ################################################### df <- data.frame(test, pID = as.factor(ms$level[[pLevel]]$partition) ) lm1 <- lm(E ~ . * pID, data = df) ################################################### ### select a Morse-SMale complex from the persistence hierarchy and build an elastic net regrssion model ################################################### ms15 <- msc.sublevels(ms, pLevel) ms.l1 <- msc.elnet(ms15) ################################################### ### show coefficents of elastic net at different levels ################################################### b1 <- as.vector(ms.l1$elnet[[1]]$lm[[1]]$beta) b5 <- as.vector(ms.l1$elnet[[1]]$lm[[5]]$beta) print( cbind( b1, b5) ) ################################################### ### stepwise model selection for linear models ################################################### ms.lm1 <- msc.lm(ms15, modelSelect=T) ################################################### ### Weighted linear model ################################################### ms.slm <- msc.slm(ms) ################################################### ### Selected prediction level for weighted linear model ################################################### print(ms.slm$ms$predictLevel) ################################################### ### mean squared residuals of weighted linear model ################################################### penergy <- predict(ms.slm, test) mean((test$E - penergy)^2)