fixed bug induced by computing the limits once, not each projection

ggobi · Mar 27, 2024 · d888e29 · d888e29
1 parent da2e575
commit d888e29
Show file tree

Hide file tree

Showing 2 changed files with 32 additions and 16 deletions.
diff --git a/R/display-xy.r b/R/display-xy.r
@@ -97,6 +97,21 @@ display_xy <- function(center = TRUE, axes = "center", half_range = NULL,
  init <- function(data) {
  half_range <<- compute_half_range(half_range, data, center)
  labels <<- abbreviate(colnames(data), 3)
+
+ if (!is.null(ellipse)) {
+ if (nrow(ellipse) == ncol(data)) {
+
+ if (is.null(ellc))
+ ellc <<- qchisq(0.95, ncol(data))
+ else
+ stopifnot(ellc > 0) # Needs to be positive
+ if (is.null(ellmu))
+ ellmu <<- rep(0, ncol(data))
+ else
+ stopifnot(length(ellmu) == ncol(data)) # Right dimension
+ message("Using ellc = ", format(ellc, digits = 2))
+ }
+ }
  }
 
  if (!is.null(edges)) {
@@ -157,15 +172,15 @@ display_xy <- function(center = TRUE, axes = "center", half_range = NULL,
  if (!is.null(ellipse)) {
  if (nrow(ellipse) == nrow(proj)) {
 
- if (is.null(ellc))
- ellc <- qchisq(0.95, nrow(proj))
- else
- stopifnot(ellc > 0) # Needs to be positive
- if (is.null(ellmu))
- ellmu <- rep(0, nrow(proj))
- else
- stopifnot(length(ellmu) == nrow(proj)) # Right dimension
- message("Using ellc = ", format(ellc, digits = 2))
+ # if (is.null(ellc))
+ # ellc <- qchisq(0.95, nrow(proj))
+ # else
+ # stopifnot(ellc > 0) # Needs to be positive
+ # if (is.null(ellmu))
+ # ellmu <- rep(0, nrow(proj))
+ # else
+ # stopifnot(length(ellmu) == nrow(proj)) # Right dimension
+ # message("Using ellc = ", format(ellc, digits = 2))
 
  # Project ellipse into 2D
  # Notation in paper: ellipse=A, ellinv=A^(-1),
@@ -178,11 +193,10 @@ display_xy <- function(center = TRUE, axes = "center", half_range = NULL,
  ellinv <- (evc$vectors) %*% diag(evc$values) %*% t(evc$vectors)
  e2 <- t(proj) %*% ellipse %*% proj
  evc2 <- eigen(e2)
- ell2d <- (evc2$vectors) %*% sqrt(diag(evc2$values*ellc)) %*% t(evc2$vectors)
+ ell2d <- as.matrix((evc2$vectors)) %*% diag(sqrt(evc2$values*ellc)) %*% t(as.matrix(evc2$vectors))
  #e3 <- eigen(ell2d)
  #ell2dinv <- (e3$vectors) %*% diag(e3$values) %*% t(e3$vectors)
 
-
  # Compute the points on an ellipse
  # Generate points on a circle
  sph <- geozoo::sphere.hollow(2, 200)$points
@@ -212,6 +226,7 @@ display_xy <- function(center = TRUE, axes = "center", half_range = NULL,
  cov=ellipse)
  #mdst <- mahal_dist(data, ellipse)
  anomalies <- which(mdst > ellc)
+ #cat("1 ", length(anomalies), "\n")
  points(x[anomalies,],
  col = "red",
  pch = 4,

diff --git a/R/tour-guided-anomaly.r b/R/tour-guided-anomaly.r
@@ -51,14 +51,14 @@ guided_anomaly_tour <- function(index_f, d = 2, alpha = 0.5, cooling = 0.99,
  if (nrow(ellipse) == nrow(proj)) {
 
  if (is.null(ellc))
- ellc <- qchisq(0.95, nrow(proj))
+ ellc <<- qchisq(0.95, nrow(proj))
  else
  stopifnot(ellc > 0) # Needs to be positive
  if (is.null(ellmu))
- ellmu <- rep(0, nrow(proj))
+ ellmu <<- rep(0, nrow(proj))
  else
  stopifnot(length(ellmu) == nrow(proj)) # Right dimension
- message("Using ellc = ", format(ellc, digits = 2))
+ #message("Using ellc = ", format(ellc, digits = 2))
 
  # Check which observations are outside pD ellipse
  mdst <- mahalanobis(data,
@@ -67,13 +67,14 @@ guided_anomaly_tour <- function(index_f, d = 2, alpha = 0.5, cooling = 0.99,
  #mdst <- mahal_dist(data, ellipse)
  anomalies <- which(mdst > ellc)
  stopifnot(length(anomalies) > 0)
+ #cat(length(anomalies), "\n")
 
  # Project ellipse into 2D
  evc <- eigen(ellipse) #
- ellinv <- (evc$vectors) %*% diag(evc$values) %*% t(evc$vectors)
+ ellinv <- (evc$vectors) %*% as.matrix(diag(evc$values)) %*% t(evc$vectors)
  e2 <- t(proj) %*% ellipse %*% proj
  evc2 <- eigen(e2)
- ell2d <- (evc2$vectors) %*% sqrt(diag(evc2$values*ellc)) %*% t(evc2$vectors)
+ ell2d <- as.matrix(evc2$vectors) %*% diag(sqrt(evc2$values*ellc)) %*% t(as.matrix(evc2$vectors))
 
  ell2dinv <- (evc2$vectors) %*% diag(evc2$values*ellc) %*% t(evc2$vectors)
  ellmu2d <- t(as.matrix(ellmu)) %*% proj