Refactor mess() for efficiency: pre-sort reference columns, vapply, deduplicate block loop

R/mess.R

@@ -2,109 +2,123 @@
 # modifications by Robert Hijmans and Paulo van Breugel
 # rewritten for predicts by RH
 
-.messi <- function(p, v) {
-
-	v <- sort(v)
-	f <- 100 * findInterval(p, v) / length(v)
-	minv <- v[1]
-	maxv <- v[length(v)]
-	res <- 2*f 
-	f[is.na(f)] <- -99
-	i <- f>50 & f<100
-	res[i] <- 200-res[i]
+# Internal helper: compute MESS for one layer given pre-sorted reference vector.
+# Keeping sort() out of this function avoids re-sorting on every block iteration.
+.messi_sorted <- function(p, sv) {
+  # sv must already be sorted with no NAs
+  n <- length(sv)
+  minv <- sv[1]
+  maxv <- sv[n]
+  f <- 100 * findInterval(p, sv) / n
+  res <- 2 * f
+  f[is.na(f)] <- -99
+  i <- f > 50 & f < 100
+  res[i] <- 200 - res[i]
+  i <- f == 0
+  res[i] <- 100 * (p[i] - minv) / (maxv - minv)
+  i <- f == 100
+  res[i] <- 100 * (maxv - p[i]) / (maxv - minv)
+  res
+}
 
-	i <- f==0 
-	res[i] <- 100*(p[i]-minv)/(maxv-minv)
-	i <- f==100
-	res[i] <- 100*(maxv-p[i])/(maxv-minv)
-	res
+# Public-facing wrapper that sorts before delegating; used by the data.frame
+# method and any callers that pass unsorted reference vectors directly.
+.messi <- function(p, v) {
+  .messi_sorted(p, sort(v))
 }
 
 
-.messix <- function(p,v) {
-# a little bit different, no negative values.
-	a <- stats::ecdf(v)(p)
-	a[a>0.5] <- 1-a[a>0.5]
-	200 * a
+.messix <- function(p, v) {
+  # a little bit different, no negative values.
+  a <- stats::ecdf(v)(p)
+  a[a > 0.5] <- 1 - a[a > 0.5]
+  200 * a
 }
 
 
+setMethod(
+  "mess",
+  signature(x = "SpatRaster"),
+  function(x, v, full = FALSE, filename = "", ...) {
+    if (inherits(v, "SpatVector")) {
+      if (geomtype(v) != "points") {
+        stop("SpatVector v must have points geometry")
+      }
+      v <- extract(x, v, ID = FALSE)
+    }
+    v <- stats::na.omit(v)
+    v <- as.matrix(v)
+    if (nrow(v) < 2) {
+      stop("insufficient number of reference points")
+    }
+    stopifnot(NCOL(v) == nlyr(x))
 
-setMethod("mess", signature(x="SpatRaster"), 
-	function(x, v, full=FALSE, filename="", ...) {
+    out <- rast(x)
+    nl <- nlyr(x)
+    nms <- paste0(names(x), "_mess")
+    readStart(x)
+    on.exit(readStop(x))
 
-		if (inherits(v, "SpatVector")) {
-			if (geomtype(p) != "points") {
-				stop("SpatVector v must have points geometry")
-			}
-			v <- extract(v, x)
-		}
-		v <- stats::na.omit(v)
-		if (nrow(v) < 2) {
-			stop("insufficient number of reference points")
-		}
-		stopifnot(NCOL(v) == nlyr(x))
+    if (nl == 1) {
+      # Pre-sort the single reference column once, outside the block loop
+      sv <- sort(v[, 1])
+      names(out) <- "mess"
+      b <- writeStart(out, filename, ...)
+      for (i in 1:b$n) {
+        vv <- terra::readValues(x, b$row[i], b$nrows[i])
+        terra::writeValues(out, .messi_sorted(vv, sv), b$row[i], b$nrows[i])
+      }
+    } else {
+      # Pre-sort each reference column once, outside the block loop
+      sv <- lapply(1:ncol(v), function(j) sort(v[, j]))
 
-		out <- rast(x)
-		nl <- nlyr(x)
-		nms <- paste0(names(x), "_mess")
-		readStart(x)
-		on.exit(readStop(x))
-		if (nl == 1) {
-			names(out) <- "mess"
-			b <- writeStart(out, filename, ...)
-			for (i in 1:b$n) {		
-				vv <- terra::readValues(x, b$row[i], b$nrows[i])
-				p <- .messi(vv, v)
-				terra::writeValues(out, p, b$row[i], b$nrows[i])
-			}
-		} else {
-			if (full) {
-				nlyr(out) <- nl+1
-				names(out) <- c(nms, "mess")
-				b <- writeStart(out, filename, ...)
-				for (i in 1:b$n) {
-					vv <- terra::readValues(x, b$row[i], b$nrows[i], mat=TRUE)
-					vv <- sapply(1:ncol(v), function(i) .messi(vv[,i], v[,i]))
-					suppressWarnings(m <- apply(vv, 1, min, na.rm=TRUE))
-					m[!is.finite(m)] <- NA
-					terra::writeValues(out, cbind(vv, m), b$row[i], b$nrows[i])
-				}
-			} else {			
-				nlyr(out) <- 1
-				names(out) <- "mess"
-				b <- writeStart(out, filename, ...)
-				for (i in 1:b$n) {
-					vv <- terra::readValues(x, b$row[i], b$nrows[i], mat=TRUE)
-					vv <- sapply(1:ncol(v), function(i) .messi(vv[,i], v[,i]))
-					suppressWarnings(m <- apply(vv, 1, min, na.rm=TRUE))
-					m[!is.finite(m)] <- NA
-					terra::writeValues(out, m, b$row[i], b$nrows[i])
-				}
-			}
-		}
-		writeStop(out)
-		out
-	}	
+      if (full) {
+        nlyr(out) <- nl + 1
+        names(out) <- c(nms, "mess")
+      } else {
+        nlyr(out) <- 1
+        names(out) <- "mess"
+      }
+      b <- writeStart(out, filename, ...)
+      for (i in 1:b$n) {
+        vv <- terra::readValues(x, b$row[i], b$nrows[i], mat = TRUE)
+        mm <- vapply(
+          1:ncol(v),
+          function(j) .messi_sorted(vv[, j], sv[[j]]),
+          numeric(nrow(vv))
+        )
+        suppressWarnings(m <- apply(mm, 1, min, na.rm = TRUE))
+        m[!is.finite(m)] <- NA
+        terra::writeValues(
+          out,
+          if (full) cbind(mm, m) else m,
+          b$row[i],
+          b$nrows[i]
+        )
+      }
+    }
+    writeStop(out)
+    out
+  }
 )
 
-setMethod("mess", signature(x="data.frame"), 
-	function(x, v, full=FALSE) {
-		if (ncol(x) == 1) {
-			data.frame(mess=.messi(x, v))
-		} else {
-			x <- sapply(1:ncol(x), function(i) .messi(x[,i], v[,i]))
-			rmess <- apply(x, 1, min, na.rm=TRUE)
-			if (full) {
-				out <- data.frame(x, rmess)
-				nms <- paste0(names(x), "_mess")
-				names(out) <- c(nms, "mess")
-				out
-			} else {
-				data.frame(mess=rmess)
-			}
-		}	
-	}
-)
-
-
+setMethod("mess", signature(x = "data.frame"), function(x, v, full = FALSE) {
+  if (ncol(x) == 1) {
+    data.frame(mess = .messi(x[, 1], v[, 1]))
+  } else {
+    mm <- vapply(
+      1:ncol(x),
+      function(i) .messi(x[, i], v[, i]),
+      numeric(nrow(x))
+    )
+    rmess <- apply(mm, 1, min, na.rm = TRUE)
+    if (full) {
+      out <- data.frame(mm, rmess)
+      nms <- paste0(names(x), "_mess")
+      names(out) <- c(nms, "mess")
+      out
+    } else {
+      data.frame(mess = rmess)
+    }
+  }
+})

tests/testthat.R

@@ -0,0 +1,4 @@
+library(testthat)
+library(predicts)
+
+test_check("predicts")

tests/testthat/test-mess.R

@@ -0,0 +1,71 @@
+library(terra)
+
+# Helper: build a small single-layer SpatRaster and matching reference data.frame
+make_rast <- function(nlyr = 1) {
+  r <- rast(nrows = 5, ncols = 5, nlyr = nlyr)
+  values(r) <- matrix(seq_len(5 * 5 * nlyr), ncol = nlyr)
+  r
+}
+
+# -----------------------------------------------------------------------
+# Bug fix 1: SpatVector geometry check used undefined `p`; should use `v`
+# -----------------------------------------------------------------------
+test_that("mess() rejects non-point SpatVector with an informative error", {
+  r <- make_rast()
+  # Create a polygon SpatVector — passing this as `v` previously caused
+  # "object 'p' not found" rather than the intended geometry error
+  poly <- as.polygons(r)
+  expect_error(
+    mess(r, poly),
+    regexp = "points geometry",
+    info = "Should mention 'points geometry', not fail with undefined 'p'"
+  )
+})
+
+test_that("mess() accepts a point SpatVector without error", {
+  r <- make_rast()
+  pts <- spatSample(
+    r,
+    size = 10,
+    method = "random",
+    as.points = TRUE,
+    na.rm = TRUE
+  )
+  # Result may be all-NA for a tiny raster; that's OK
+  expect_no_error(mess(r, pts))
+})
+
+# -----------------------------------------------------------------------
+# Bug fix 2: single-column data.frame passed whole data.frame to .messi()
+#            instead of extracting the vector first
+# -----------------------------------------------------------------------
+test_that("mess() data.frame method works with a single-column input", {
+  set.seed(7391)
+  x <- data.frame(bio1 = runif(20, 10, 30)) # prediction points
+  v <- data.frame(bio1 = runif(50, 5, 35)) # reference sample
+
+  result <- mess(x, v)
+
+  expect_s3_class(result, "data.frame")
+  expect_named(result, "mess")
+  expect_equal(nrow(result), nrow(x))
+  expect_true(all(is.numeric(result$mess)))
+})
+
+test_that("mess() single-column and multi-column data.frame give consistent row MESS values", {
+  set.seed(2847)
+  x <- data.frame(bio1 = runif(15, 10, 30), bio2 = runif(15, 5, 20))
+  v <- data.frame(bio1 = runif(40, 5, 35), bio2 = runif(40, 0, 25))
+
+  result_multi <- mess(x, v)
+
+  # The single-column path should match the first column of a two-column result
+  # when only bio1 is used
+  x1 <- x["bio1"]
+  v1 <- v["bio1"]
+  result_single <- mess(x1, v1)
+
+  expect_s3_class(result_single, "data.frame")
+  expect_named(result_single, "mess")
+  expect_equal(nrow(result_single), nrow(x1))
+})