From d7f6bba216f964a70a9eccf007cbd373efde1ca3 Mon Sep 17 00:00:00 2001 From: Afonso Ferrer <13241113+fonzie42@users.noreply.github.com> Date: Wed, 10 Jun 2026 00:33:59 -0300 Subject: [PATCH] fix(kernels): correct swapped location indices in UnivariateMaternNuggetsStationary GenerateCovarianceMatrix passed (j0, i0) to CalculateDistance instead of (i0, j0), swapping which location set each loop index addresses. Every other kernel in the project passes (i0, j0). The defect is invisible for the symmetric train-train covariance matrix (Euclidean distance is symmetric, so d(loc[i], loc[j]) == d(loc[j], loc[i])) which is why data generation and MLE are unaffected. It does corrupt the rectangular cross-covariance matrix C12 assembled during prediction, where the two location sets (observed and missing) differ: entry (i, j) ends up using d(location1[j], location2[i]) instead of d(location1[i], location2[j]), producing the transposed matrix and therefore wrong kriging predictions. Add an isolated regression test that builds a 3x3 cross-covariance between two distinct location sets and checks every entry against an independently computed reference (Matern nu=0.5 reduces to sigma^2 * exp(-d/beta)). On the old code entry (0,1) returns exp(-1) instead of exp(-3); with the fix all entries match. The existing data-generation test for this kernel continues to pass. --- .../UnivariateMaternNuggetsStationary.cpp | 4 +- ...UnivariateMaternNuggetsCrossCovariance.cpp | 108 ++++++++++++++++++ 2 files changed, 110 insertions(+), 2 deletions(-) create mode 100644 tests/cpp-tests/kernels/concrete/TestUnivariateMaternNuggetsCrossCovariance.cpp diff --git a/src/kernels/concrete/UnivariateMaternNuggetsStationary.cpp b/src/kernels/concrete/UnivariateMaternNuggetsStationary.cpp index ab960043..4b0420cb 100644 --- a/src/kernels/concrete/UnivariateMaternNuggetsStationary.cpp +++ b/src/kernels/concrete/UnivariateMaternNuggetsStationary.cpp @@ -59,7 +59,7 @@ void UnivariateMaternNuggetsStationary::GenerateCovarianceMatrix(T *apMatrixA for (i = 0; i < aRowsNumber; i++) { j0 = aColumnOffset; for (j = 0; j < aColumnsNumber; j++) { - expr = DistanceCalculationHelpers::CalculateDistance(aLocation1, aLocation2, j0, i0, aDistanceMetric, + expr = DistanceCalculationHelpers::CalculateDistance(aLocation1, aLocation2, i0, j0, aDistanceMetric, flag) / aLocalTheta[1]; if (expr == 0) { apMatrixA[i + j * aRowsNumber] = sigma_square + aLocalTheta[3]; @@ -77,7 +77,7 @@ void UnivariateMaternNuggetsStationary::GenerateCovarianceMatrix(T *apMatrixA j0 = aColumnOffset; for (j = 0; j < aColumnsNumber; j++) { flag = 1; - expr = DistanceCalculationHelpers::CalculateDistance(aLocation1, aLocation2, j0, i0, aDistanceMetric, + expr = DistanceCalculationHelpers::CalculateDistance(aLocation1, aLocation2, i0, j0, aDistanceMetric, flag); if (expr == 0) { apMatrixA[i + j * aRowsNumber] = sigma_square + aLocalTheta[3]; diff --git a/tests/cpp-tests/kernels/concrete/TestUnivariateMaternNuggetsCrossCovariance.cpp b/tests/cpp-tests/kernels/concrete/TestUnivariateMaternNuggetsCrossCovariance.cpp new file mode 100644 index 00000000..bf6a029d --- /dev/null +++ b/tests/cpp-tests/kernels/concrete/TestUnivariateMaternNuggetsCrossCovariance.cpp @@ -0,0 +1,108 @@ + +// Copyright (c) 2017-2024 King Abdullah University of Science and Technology, +// All rights reserved. +// ExaGeoStat is a software package, provided by King Abdullah University of Science and Technology (KAUST). + +/** + * @file TestUnivariateMaternNuggetsCrossCovariance.cpp + * @brief Regression test for the location-index order of the UnivariateMaternNuggetsStationary kernel. + * @details The existing kernel test exercises only the symmetric (train-train) covariance matrix, where + * the matrix entry (i, j) is computed from the distance between two points drawn from the SAME location + * set. Because the Euclidean distance is symmetric, an accidental swap of the two location indices is + * invisible there: d(loc[i], loc[j]) == d(loc[j], loc[i]). + * + * This test instead builds a rectangular CROSS-covariance matrix between two DIFFERENT location sets, + * exactly as the prediction path does when assembling C12 between observed and missing points. For two + * distinct location sets, entry (i, j) must use d(location1[i], location2[j]); swapping the indices to + * d(location1[j], location2[i]) produces a different (effectively transposed) matrix and therefore wrong + * predictions. + * + * The expected values are computed independently of the kernel, using the closed-form Matern covariance + * for nu = 0.5, which reduces to the exponential covariance sigma^2 * exp(-d / beta). This makes the + * reference values self-evident and independent of any internal index convention. + * @version 1.1.0 +**/ + +#include +#include + +#include + +#include +#include + +using namespace std; + +using namespace exageostat::common; +using namespace exageostat::dataunits; +using namespace exageostat::kernels; +using namespace exageostat::plugins; + +void TEST_NUGGETS_CROSS_COVARIANCE() { + + SECTION("UnivariateMaternNuggetsStationary cross-covariance respects location-index order") { + + // Matern parameters: sigma^2 = 1, beta (range) = 1, nu = 0.5, nugget = 0.5. + // With nu = 0.5 the Matern covariance reduces to the exponential covariance: + // C(d) = sigma^2 * exp(-d / beta) for d > 0 + // C(0) = sigma^2 + nugget for d = 0 + double theta[4] = {1.0, 1.0, 0.5, 0.5}; + const double sigma_square = theta[0]; + const double beta = theta[1]; + const double nugget = theta[3]; + + const int rows = 3; // number of "row" points, indexed from location set 1 + const int cols = 3; // number of "column" points, indexed from location set 2 + + // Two DISTINCT location sets so the cross-covariance is asymmetric. + // Row points lie on the y-axis, column points lie on the x-axis. + vector row_x = {0.0, 0.0, 0.0}; + vector row_y = {0.0, 1.0, 2.0}; + vector col_x = {0.0, 3.0, 6.0}; + vector col_y = {0.0, 0.0, 0.0}; + + Locations location1(rows, Dimension2D); + location1.SetLocationX(*row_x.data(), rows); + location1.SetLocationY(*row_y.data(), rows); + + Locations location2(cols, Dimension2D); + location2.SetLocationX(*col_x.data(), cols); + location2.SetLocationY(*col_y.data(), cols); + + auto *kernel = PluginRegistry>::Create("UnivariateMaternNuggetsStationary", 1); + REQUIRE(kernel != nullptr); + + // Output matrix in column-major layout: A[i + j * rows]. + vector matrix(rows * cols, 0.0); + + // Distance metric 0 -> Euclidean (metric 1 is reserved for great-circle/Haversine). + const int distance_metric = 0; + kernel->GenerateCovarianceMatrix(matrix.data(), rows, cols, /*row offset*/ 0, /*column offset*/ 0, + location1, location2, location1, theta, distance_metric); + + // Build the reference cross-covariance independently of the kernel. + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + const double dx = row_x[i] - col_x[j]; + const double dy = row_y[i] - col_y[j]; + const double distance = sqrt(dx * dx + dy * dy); + const double expected = (distance == 0.0) ? (sigma_square + nugget) + : sigma_square * exp(-distance / beta); + const double actual = matrix[i + j * rows]; + INFO("entry (" << i << ", " << j << "): distance = " << distance); + REQUIRE(actual == Catch::Approx(expected).margin(1e-9)); + } + } + + // The single most diagnostic entry (0, 1): row point (0,0), column point (3,0). + // Correct : d(location1[0]=(0,0), location2[1]=(3,0)) = 3 -> exp(-3) ~= 0.0498 + // Swapped : d(location1[1]=(0,1), location2[0]=(0,0)) = 1 -> exp(-1) ~= 0.3679 (the bug) + REQUIRE(matrix[0 + 1 * rows] == Catch::Approx(exp(-3.0)).margin(1e-9)); + + delete kernel; + } +} + +TEST_CASE("UnivariateMaternNuggetsStationary cross-covariance index-order regression") { + TEST_NUGGETS_CROSS_COVARIANCE(); +}