From c77c1c7275e806a6e4008e4dc67496ed6d68e799 Mon Sep 17 00:00:00 2001
From: anshul23102 <anshul23102@iiitd.ac.in>
Date: Sat, 16 May 2026 13:10:21 +0530
Subject: [PATCH] fix(emulator): vectorize calcSpatialCov.matrix, removing
 O(n^3) loop

The inner `for (j in seq_len(nl))` loop was entirely redundant: it
overwrote `H[i, ]` with the same vector on every iteration, doing
nl times more work than necessary. The outer loop over rows was also
unnecessary because R's matrix arithmetic applies the operation to
all elements at once.

Replace both loops with a single vectorized expression:

    tau * exp(-psi * d)

This gives identical results in O(n^2) instead of O(n^3). Add tests
covering correctness, dimensions, symmetry, and agreement with the
already-vectorized calcSpatialCov.list implementation.

Fixes #3964.
---
 modules/emulator/NEWS.md                      |  4 ++
 modules/emulator/R/calcSpatialCov.matrix.R    | 10 +----
 .../tests/testthat/test-calcSpatialCov.R      | 44 +++++++++++++++++++
 3 files changed, 49 insertions(+), 9 deletions(-)
 create mode 100644 modules/emulator/tests/testthat/test-calcSpatialCov.R

diff --git a/modules/emulator/NEWS.md b/modules/emulator/NEWS.md
index 0fd4703c018..19ca3675375 100644
--- a/modules/emulator/NEWS.md
+++ b/modules/emulator/NEWS.md
@@ -2,6 +2,10 @@
 
 * Added keywords and bug reporting URL to DESCRIPTION. No other files changed.
 
+## Fixed
+
+* `calcSpatialCov.matrix()`: removed a redundant nested loop that recomputed each matrix row `n` times unnecessarily, reducing complexity from O(n^3) to O(n^2) (#3964).
+
 
 
 # PEcAn.emulator 1.8.1
diff --git a/modules/emulator/R/calcSpatialCov.matrix.R b/modules/emulator/R/calcSpatialCov.matrix.R
index 7393ff98b9e..224073b76c2 100644
--- a/modules/emulator/R/calcSpatialCov.matrix.R
+++ b/modules/emulator/R/calcSpatialCov.matrix.R
@@ -10,13 +10,5 @@
 ##' @author Michael Dietze
 ##' @export
 calcSpatialCov.matrix <- function(d, psi, tau, ...) {
-  nl <- nrow(d)
-  H  <- matrix(0, nl, nl)
-  for (i in seq_len(nl)) {
-    # for(j in 1:nl){ H[i,j] <- tau*exp(-psi*d[i,j]) }
-    for (j in seq_len(nl)) {
-      H[i, ] <- tau * exp(-psi * d[i, ])
-    }
-  }
-  return(H)
+  tau * exp(-psi * d)
 }
diff --git a/modules/emulator/tests/testthat/test-calcSpatialCov.R b/modules/emulator/tests/testthat/test-calcSpatialCov.R
new file mode 100644
index 00000000000..5bedf75614f
--- /dev/null
+++ b/modules/emulator/tests/testthat/test-calcSpatialCov.R
@@ -0,0 +1,44 @@
+test_that("calcSpatialCov.matrix returns correct exponential covariance", {
+  d <- matrix(c(0, 1, 2,
+                1, 0, 1,
+                2, 1, 0), nrow = 3)
+  psi <- 0.5
+  tau <- 2.0
+
+  result <- calcSpatialCov.matrix(d, psi = psi, tau = tau)
+
+  expected <- tau * exp(-psi * d)
+  expect_equal(result, expected)
+})
+
+test_that("calcSpatialCov.matrix returns a matrix of the same dimensions as d", {
+  d <- matrix(runif(25), nrow = 5)
+  result <- calcSpatialCov.matrix(d, psi = 1, tau = 1)
+  expect_equal(dim(result), dim(d))
+})
+
+test_that("calcSpatialCov.matrix diagonal is tau when d has zeros on diagonal", {
+  n <- 4
+  d <- matrix(1, n, n)
+  diag(d) <- 0
+  tau <- 3.0
+  result <- calcSpatialCov.matrix(d, psi = 1, tau = tau)
+  expect_equal(diag(result), rep(tau, n))
+})
+
+test_that("calcSpatialCov.matrix is symmetric when d is symmetric", {
+  d <- as.matrix(dist(1:5))
+  result <- calcSpatialCov.matrix(d, psi = 0.3, tau = 1.5)
+  expect_equal(result, t(result))
+})
+
+test_that("calcSpatialCov.list and calcSpatialCov.matrix agree on single-component case", {
+  d <- as.matrix(dist(1:4))
+  psi <- 0.4
+  tau <- 2.0
+
+  result_matrix <- calcSpatialCov.matrix(d, psi = psi, tau = tau)
+  result_list   <- calcSpatialCov.list(list(d), psi = psi, tau = tau)
+
+  expect_equal(result_matrix, result_list)
+})