Merge pull request #129 from huizezhang-sherry/master

jellyfish optimiser and more indexes

DESCRIPTION

@@ -34,7 +34,9 @@ Suggests:
  rmarkdown,
  tidyr,
  covr,
- plotly
+ plotly,
+ cassowaryr,
+ minerva
 License: MIT + file LICENSE
 LazyData: true
 URL: https://github.com/ggobi/tourr, https://ggobi.github.io/tourr/

NAMESPACE

@@ -12,6 +12,8 @@ S3method(plot,path_index)
 S3method(print,history_array)
 S3method(print,tour_path)
 S3method(str,history_array)
+export(MIC)
+export(TIC)
 export(anchored_orthogonal_distance)
 export(andrews)
 export(angular_breaks)
@@ -41,6 +43,7 @@ export(center)
 export(cmass)
 export(cumulative_radial)
 export(dcor2d)
+export(dcor2d_2)
 export(dependence_tour)
 export(display_andrews)
 export(display_density2d)
@@ -77,6 +80,7 @@ export(lda_pp)
 export(linear_breaks)
 export(little_tour)
 export(local_tour)
+export(loess2d)
 export(mahal_dist)
 export(manual_slice)
 export(mapColors)
@@ -106,12 +110,14 @@ export(save_history)
 export(search_better)
 export(search_better_random)
 export(search_geodesic)
+export(search_jellyfish)
 export(search_polish)
 export(search_posse)
 export(skewness)
 export(slice_index)
 export(sphere_data)
 export(splines2d)
+export(stringy)
 export(thaw)
 importFrom(grDevices,dev.cur)
 importFrom(grDevices,dev.flush)

R/animate.r

@@ -73,9 +73,11 @@ animate <- function(data, tour_path = grand_tour(), display = display_xy(),
  bases <- array(NA, c(ncol(data), ncol(start$target), bs))
 
  # Initialise display
- display$init(data)
- display$render_frame()
- display$render_data(data, start$proj)
+ if (!is.null(start$target)){
+ display$init(data)
+ display$render_frame()
+ display$render_data(data, start$proj)
+ }
 
  b <- 0
  i <- 0
@@ -93,19 +95,23 @@ animate <- function(data, tour_path = grand_tour(), display = display_xy(),
  bs <- 2 * bs
  }
  bases[, , b] <- step$target
+ } else if (is.null(step$proj)){
+ break
  }
 
- dev.hold()
- on.exit(dev.flush())
- if (plat$os == "win" || plat$iface == "rstudio") {
- display$render_frame()
- } else {
- display$render_transition()
+ if (!is.null(start$target)){
+ dev.hold()
+ on.exit(dev.flush())
+ if (plat$os == "win" || plat$iface == "rstudio") {
+ display$render_frame()
+ } else {
+ display$render_transition()
+ }
+ display$render_data(data, step$proj, step$target)
+ dev.flush()
+ if (step$step < 0) break # break after rendering final projection
+ Sys.sleep(1 / fps)
  }
- display$render_data(data, step$proj, step$target)
- dev.flush()
- if (step$step < 0) break # break after rendering final projection
- Sys.sleep(1 / fps)
  }
  },
  interrupt = function(cond) {

R/geodesic-path.r

@@ -35,6 +35,7 @@ new_geodesic_path <- function(name, generator, frozen = NULL, ...) {
 
  gen <- generator(current, data, tries, ...)
  target <- gen$target
+ if (inherits(target, "multi-bases")) return(list(target = target))
 
  # generator has run out, so give up
  if (is.null(target)) {
@@ -46,7 +47,6 @@ new_geodesic_path <- function(name, generator, frozen = NULL, ...) {
  return(NULL)
  }
 
- #cat("generation: dist = ", dist, "\n")
  }
  list(ingred = geodesic_path(current, target, frozen, ...), index = gen$index, tries = tries)
  }
@@ -61,7 +61,7 @@ new_geodesic_path <- function(name, generator, frozen = NULL, ...) {
 
 #' @export
 "print.tour_path" <- function(x, ...) {
- cat("Tour path:", attr(x, "name"), "\n")
+ message("Tour path: ", attr(x, "name"))
 
  # params <- as.list(environment(get("generator", environment(g))))
  # str(params)

R/geodesic.r

@@ -99,7 +99,7 @@ geodesic_info <- function(Fa, Fz, epsilon = 1e-6) {
  }
 
  # if (Fa.equivalent(Fz)) return();
- # cat("dim Fa",nrow(Fa),ncol(Fa),"dim Fz",nrow(Fz),ncol(Fz),"\n")
+ # message("dim Fa",nrow(Fa),ncol(Fa),"dim Fz",nrow(Fz),ncol(Fz))
 
  # Compute the SVD: Fa'Fz = Va lambda Vz' --------------------------------
  sv <- svd(t(Fa) %*% Fz)

R/interesting-indices.r

@@ -1,11 +1,45 @@
+#' Scagnostic indexes.
+#'
+#' Compute the scagnostic measures from the cassowaryr package
+#' @export
+stringy <- function(){
+ function(mat){
+ cassowaryr::sc_stringy(mat[,1], mat[,2])
+ }
+}
+
+
+#' Maximum and total information coefficient index.
+#'
+#' Compute the maximum and total information coefficient indexes,
+#' see \code{minerva::mine}.
+#'
+#' @rdname indexes
+#' @export
+MIC <- function(){
+ function(mat){
+ minerva::mine(mat[,1], mat[,2], alpha = 0.3, est = "mic_e")$MIC
+ }
+}
+
+#' @rdname indexes
+#' @export
+TIC <- function(){
+ function(mat){
+ minerva::mine(mat[,1], mat[,2], est = "mic_e", alpha = 0.3)$TIC
+ }
+}
+
 #' Distance correlation index.
 #'
-#' Computes the distance correlation based index on
-#' 2D projections of the data.
+#' Computes the distance correlation based index on 2D projections of the data.
+#' \code{dcor2d_2} uses the faster implementation of the distance correlation
+#' for bivariate data, see \code{energy::dcor2d}.
 #'
 #' @keywords hplot
 #' @importFrom stats na.omit
 #' @export
+#' @rdname dcor
 dcor2d <- function() {
  function(mat) {
  xy <- na.omit(data.frame(x = mat[, 1], y = mat[, 2]))
@@ -14,15 +48,26 @@ dcor2d <- function() {
  }
 }
 
-#' Spline based index.
+#' @rdname dcor
+#' @export
+dcor2d_2 <- function() {
+ function(mat) {
+ xy <- na.omit(data.frame(x = mat[, 1], y = mat[, 2]))
+ measure <- with(xy, energy::dcor2d(x, y, type = "U"))
+ return(measure)
+ }
+}
+
+#' Spline/loess based index.
 #'
 #' Compares the variance in residuals of a fitted
-#' spline model to the overall variance to find
+#' spline/loess model to the overall variance to find
 #' functional dependence in 2D projections
 #' of the data.
 #'
 #' @keywords hplot
 #' @importFrom stats residuals var
+#' @rdname spline-loess
 #' @export
 splines2d <- function() {
  function(mat) {
@@ -38,6 +83,20 @@ splines2d <- function() {
  }
 }
 
+#' @rdname spline-loess
+#' @export
+loess2d <- function() {
+ function(mat) {
+ mat <- as.data.frame(mat)
+ colnames(mat) <- c("x", "y")
+ loess_fit <- loess(y ~ x, data = mat, span = 0.05)
+ loess_fit2 <- loess(x ~ y, data = mat, span = 0.05)
+ measure <- max(1 - var(residuals(loess_fit), na.rm = T) / var(mat$y, na.rm = T),
+ 1 - var(residuals(loess_fit2), na.rm = T) / var(mat$y, na.rm = T)
+ )
+ return(measure)
+ }
+}
 
 
 #' Normality index.

R/search-better.r

@@ -48,7 +48,6 @@ search_better <- function(current, alpha = 0.5, index, tries, max.tries = Inf,..
  warning("cur_index is zero!")
  }
 
- cat("Old", cur_index, "\n")
  try <- 1
 
  while (try < max.tries) {
@@ -70,7 +69,7 @@ search_better <- function(current, alpha = 0.5, index, tries, max.tries = Inf,..
 
 
  if (new_index > cur_index) {
- cat("New", new_index, "try", try, "\n")
+ message("Target: ", sprintf("%.3f", new_index), ", try: ", try)
 
  nr <- nrow(rcd_env[["record"]])
  rcd_env[["record"]][nr, "info"] <- "new_basis"
@@ -83,24 +82,8 @@ search_better <- function(current, alpha = 0.5, index, tries, max.tries = Inf,..
  try <- try + 1
  }
 
- cat("No better bases found after ", max.tries, " tries. Giving up.\n",
- sep = ""
- )
- cat("Final projection: \n")
- if (ncol(current) == 1) {
- for (i in 1:length(current)) {
- cat(sprintf("%.3f", current[i]), " ")
- }
- cat("\n")
- }
- else {
- for (i in 1:nrow(current)) {
- for (j in 1:ncol(current)) {
- cat(sprintf("%.3f", current[i, j]), " ")
- }
- cat("\n")
- }
- }
+ message("No better bases found after ", max.tries, " tries. Giving up.")
+ print_final_proj(current)
 
  rcd_env[["record"]] <- dplyr::mutate(
  rcd_env[["record"]],
@@ -139,8 +122,6 @@ search_better_random <- function(current, alpha = 0.5, index, tries,
  if (cur_index == 0) {
  warning("cur_index is zero!")
  }
-
- cat("Old", cur_index, "\n")
  try <- 1
  while (try < max.tries) {
  new_basis <- basis_nearby(current, alpha, method)
@@ -160,8 +141,7 @@ search_better_random <- function(current, alpha = 0.5, index, tries,
  )
 
  if (new_index > cur_index) {
- cat("New", new_index, "try", try, "\n")
- cat("Accept \n")
+ message("Target: ", sprintf("%.3f", new_index), ", try: ", try, ", accept")
 
  nr <- nrow(rcd_env[["record"]])
  rcd_env[["record"]][nr, "info"] <- "new_basis"
@@ -176,8 +156,8 @@ search_better_random <- function(current, alpha = 0.5, index, tries,
  rand <- stats::runif(1)
 
  if (prob > rand) {
- cat("New", new_index, "try", try, "\n")
- cat("Accept with probability, prob =", prob, "\n")
+ message("Target: ", sprintf("%.3f", new_index), ", try: ", try,
+  ", probabilistic accept p = ", sprintf("%.3f", prob))
 
  nr <- nrow(rcd_env[["record"]])
  rcd_env[["record"]][nr, "info"] <- "new_basis"
@@ -193,24 +173,8 @@ search_better_random <- function(current, alpha = 0.5, index, tries,
  try <- try + 1
  }
 
- cat("No better bases found after ", max.tries, " tries. Giving up.\n",
- sep = ""
- )
- cat("Final projection: \n")
- if (ncol(current) == 1) {
- for (i in 1:length(current)) {
- cat(sprintf("%.3f", current[i]), " ")
- }
- cat("\n")
- }
- else {
- for (i in 1:nrow(current)) {
- for (j in 1:ncol(current)) {
- cat(sprintf("%.3f", current[i, j]), " ")
- }
- cat("\n")
- }
- }
+ message("No better bases found after ", max.tries, " tries. Giving up.")
+ print_final_proj(current)
 
  rcd_env[["record"]] <- dplyr::mutate(
  rcd_env[["record"]],

R/search-frozen-geodesic.r

@@ -44,13 +44,13 @@ search_frozen_geodesic <- function(current, index, tries, max.tries = 5, n = 5,
  dig3 <- function(x) sprintf("%.3f", x)
  pdiff <- (new_index - cur_index) / cur_index
  if (pdiff > 0.001) {
- cat("New index: ", dig3(new_index), " (", dig3(peak$alpha$maximum), " away)\n", sep = "")
+ message("New index: ", dig3(new_index), " (", dig3(peak$alpha$maximum), " away)", sep = "")
  current <<- new_basis
  cur_index <<- new_index
 
  return(list(target = new_basis[[1]]))
  }
- cat("Best was: ", dig3(new_index), " (", dig3(peak$alpha$maximum), " away). Trying again...\n", sep = "")
+ message("Best was: ", dig3(new_index), " (", dig3(peak$alpha$maximum), " away). Try again...", sep = "")
 
  try <- try + 1
  }

R/search-geodesic.r

@@ -58,44 +58,21 @@ search_geodesic <- function(current, alpha = 1, index, tries, max.tries = 5, ...
  warning("either the cur_index or the new_index is zero!")
  } else {
  pdiff <- (new_index - cur_index) / cur_index
-
- dig3 <- function(x) sprintf("%.3f", x)
-
- cat(
- "Value ", dig3(new_index), " ",
- sprintf("%.1f", pdiff * 100), "% better "
- )
  if (pdiff > 0.001 & proj_dist(current, new_basis[[1]]) > 1e-2) { # FIXME: pdiff should pbly be a changeable parameter
- cat(" - NEW BASIS\n")
-
+ message("Target: ", sprintf("%.3f", new_index), ", ",
+ sprintf("%.1f", pdiff * 100), "% better ")
  current <- new_basis
  cur_index <- new_index
 
  return(list(target = new_basis[[1]]))
  }
- cat("\n")
  }
 
  try <- try + 1
  }
- cat("No better bases found after ", max.tries, " tries. Giving up.\n",
- sep = ""
- )
- cat("Final projection: \n")
- if (ncol(current) == 1) {
- for (i in 1:length(current)) {
- cat(sprintf("%.3f", current[i]), " ")
- }
- cat("\n")
- }
- else {
- for (i in 1:nrow(current)) {
- for (j in 1:ncol(current)) {
- cat(sprintf("%.3f", current[i, j]), " ")
- }
- cat("\n")
- }
- }
+
+ message("No better bases found after ", max.tries, " tries. Giving up.")
+ print_final_proj(current)
 
  NULL
 }