diff --git a/R/display-xy.r b/R/display-xy.r
index 3cd3963b..f2e7cccd 100644
--- 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
index 11625b26..0271b93d 100644
--- 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