From e575327b77375f68e60e90799bca042d11dc21c5 Mon Sep 17 00:00:00 2001 From: constracktor <74077030+constracktor@users.noreply.github.com> Date: Thu, 30 Jul 2026 10:52:02 +0200 Subject: [PATCH 1/3] test(gp): add unit test suite with tile-count invariance check No Catch2 unit test suite existed alongside the output_correctness integration test. Add test/src/unit_tests.cpp (wired into test/CMakeLists.txt as GPRat_test_unit) with a first test: predictions, uncertainty, and loss for fixed (untrained) hyperparameters must not depend on n_tiles, since tiling is purely a scheduling/decomposition detail. Compares n_tiles in {2,4,8} against a single-tile baseline for every test point's mean and variance, plus the loss, at machine precision. --- test/CMakeLists.txt | 10 ++++++ test/src/unit_tests.cpp | 80 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+) create mode 100644 test/src/unit_tests.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 94eb9d10..2cf54789 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -55,3 +55,13 @@ add_test( NAME GPRat_test_output_correctness COMMAND GPRat_test_output_correctness WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}") + +add_executable(GPRat_test_unit src/unit_tests.cpp) +target_link_libraries(GPRat_test_unit PRIVATE GPRat::core + Catch2::Catch2WithMain) +target_compile_features(GPRat_test_unit PRIVATE cxx_std_17) + +add_test( + NAME GPRat_test_unit + COMMAND GPRat_test_unit + WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}") diff --git a/test/src/unit_tests.cpp b/test/src/unit_tests.cpp new file mode 100644 index 00000000..0369fbc5 --- /dev/null +++ b/test/src/unit_tests.cpp @@ -0,0 +1,80 @@ +#include "gprat_c.hpp" +#include "target.hpp" +#include "utils_c.hpp" +#include +#include +#include +#include +#include +using Catch::Matchers::WithinRel; + +namespace +{ +// Starts the HPX runtime on construction and stops it on destruction so that +// stop_hpx_runtime() is always called even when a test assertion fails mid-test. +struct hpx_runtime_guard +{ + hpx_runtime_guard() { utils::start_hpx_runtime(0, nullptr); } + + ~hpx_runtime_guard() { utils::stop_hpx_runtime(); } +}; +} // namespace + +namespace gprat::test +{ + +static std::string gprat_data_root() +{ + const char *env = std::getenv("GPRAT_ROOT"); + return env ? env : "../data"; +} + +TEST_CASE("GP: results are tile-count invariant", "[unit][gp][predict][tiling]") +{ + // Tiling is purely a scheduling/decomposition detail: predictions, uncertainty, + // and loss for fixed (untrained) hyperparameters must not depend on n_tiles. + const std::string root = gprat_data_root(); + + constexpr int n = 128, n_reg = 8, n_test = 64; + const double eps = std::numeric_limits::epsilon() * 1'000'000; + + gprat::GP_data train_in(root + "/data_1024/training_input.txt", n, n_reg); + gprat::GP_data train_out(root + "/data_1024/training_output.txt", n, 1); + gprat::GP_data test_in(root + "/data_1024/test_input.txt", n_test, n_reg); + + hpx_runtime_guard hpx_guard; + + // Baseline: a single tile, i.e. no decomposition at all. + const int baseline_tile_size = utils::compute_train_tile_size(n, 1); + const auto [baseline_m_tiles, baseline_m_tile_size] = utils::compute_test_tiles(n_test, 1, baseline_tile_size); + gprat::GP baseline_gp( + train_in.data, train_out.data, 1, baseline_tile_size, n_reg, { 1.0, 1.0, 0.1 }, { false, false, false }); + const auto baseline_pred = + baseline_gp.predict_with_uncertainty(test_in.data, baseline_m_tiles, baseline_m_tile_size); + const double baseline_loss = baseline_gp.calculate_loss(); + + for (const int n_tiles : { 2, 4, 8 }) + { + const int tile_size = utils::compute_train_tile_size(n, n_tiles); + const auto [m_tiles, m_tile_size] = utils::compute_test_tiles(n_test, n_tiles, tile_size); + + gprat::GP gp( + train_in.data, train_out.data, n_tiles, tile_size, n_reg, { 1.0, 1.0, 0.1 }, { false, false, false }); + const auto pred = gp.predict_with_uncertainty(test_in.data, m_tiles, m_tile_size); + + for (int i = 0; i < n_test; ++i) + { + INFO("n_tiles=" << n_tiles << " mean[" << i << "]"); + REQUIRE_THAT(pred[0][static_cast(i)], + WithinRel(baseline_pred[0][static_cast(i)], eps)); + INFO("n_tiles=" << n_tiles << " variance[" << i << "]"); + REQUIRE_THAT(pred[1][static_cast(i)], + WithinRel(baseline_pred[1][static_cast(i)], eps)); + } + + INFO("n_tiles=" << n_tiles << " loss"); + REQUIRE_THAT(gp.calculate_loss(), WithinRel(baseline_loss, eps)); + } +} + +} // namespace gprat::test From 75181a65deddefc9024ef9c2cbba508cca080bbe Mon Sep 17 00:00:00 2001 From: constracktor <74077030+constracktor@users.noreply.github.com> Date: Thu, 30 Jul 2026 10:52:38 +0200 Subject: [PATCH 2/3] fix(gpu): fix predict_with_uncertainty for non-square (test != train) tile sizes Two independent bugs, both only reachable when the test tile size differs from the training tile size (i.e. n_test % n_tile_size != 0, which happens through the ordinary compute_test_tiles fallback path, not just unusual configurations): - CUDA: gen_tile_cross_cov_T's transpose<<<>>> kernel launch passed n_row_tile_size/n_column_tile_size in the wrong order. Harmless when they're equal (the swap is then a no-op), silently wrong otherwise. - SYCL: the same TransposeKernel width/height swap, plus a second, separate bug in the trsm oneMath wrapper that didn't swap its M/N dimension arguments the way the (correct) cuBLAS equivalent does for the row-major- as-column-major convention. The two bugs were masking each other's symptoms until the transpose fix exposed the trsm one. Add a regression test (n_test=48 against n_tile_size=32, which compute_test_tiles cannot round to equal sizes) since every existing GPU test picks sizes that divide evenly -- exactly why this went undetected. Verified the new test fails without its corresponding fix and passes with it, on both CUDA and SYCL hardware. --- core/src/gpu/cuda/gp_algorithms.cu | 2 +- core/src/gpu/sycl/adapter_onemath.cpp | 2 +- core/src/gpu/sycl/sycl_gp_algorithms.cpp | 2 +- test/src/unit_tests.cpp | 52 ++++++++++++++++++++++++ 4 files changed, 55 insertions(+), 3 deletions(-) diff --git a/core/src/gpu/cuda/gp_algorithms.cu b/core/src/gpu/cuda/gp_algorithms.cu index 97ac0c3e..522cfd5f 100644 --- a/core/src/gpu/cuda/gp_algorithms.cu +++ b/core/src/gpu/cuda/gp_algorithms.cu @@ -276,7 +276,7 @@ hpx::shared_future gen_tile_cross_cov_T(std::size_t n_row_tile_size, dim3 n_blocks((n_column_tile_size + BLOCK_SIZE - 1) / BLOCK_SIZE, (n_row_tile_size + BLOCK_SIZE - 1) / BLOCK_SIZE); transpose<<>>( - transposed, d_cross_covariance_tile, n_row_tile_size, n_column_tile_size); + transposed, d_cross_covariance_tile, n_column_tile_size, n_row_tile_size); check_cuda_error(cudaStreamSynchronize(stream)); diff --git a/core/src/gpu/sycl/adapter_onemath.cpp b/core/src/gpu/sycl/adapter_onemath.cpp index a2f0aa39..7d7b7881 100644 --- a/core/src/gpu/sycl/adapter_onemath.cpp +++ b/core/src/gpu/sycl/adapter_onemath.cpp @@ -72,8 +72,8 @@ double *trsm(sycl::queue queue, oneapi::math::uplo::upper, is_transposed, oneapi::math::diag::nonunit, - static_cast(M), static_cast(N), + static_cast(M), alpha, f_A, static_cast(M), diff --git a/core/src/gpu/sycl/sycl_gp_algorithms.cpp b/core/src/gpu/sycl/sycl_gp_algorithms.cpp index 98e8aa09..bbce0806 100644 --- a/core/src/gpu/sycl/sycl_gp_algorithms.cpp +++ b/core/src/gpu/sycl/sycl_gp_algorithms.cpp @@ -186,7 +186,7 @@ hpx::shared_future gen_tile_cross_cov_T(std::size_t n_row_tile_size, [&](sycl::handler &cgh) { auto kernel = - TransposeKernel(transposed, d_cross_covariance_tile, n_row_tile_size, n_column_tile_size, cgh); + TransposeKernel(transposed, d_cross_covariance_tile, n_column_tile_size, n_row_tile_size, cgh); cgh.parallel_for(sycl::nd_range<2>(global_range, local_range), kernel); }); diff --git a/test/src/unit_tests.cpp b/test/src/unit_tests.cpp index 0369fbc5..b362d158 100644 --- a/test/src/unit_tests.cpp +++ b/test/src/unit_tests.cpp @@ -29,6 +29,15 @@ static std::string gprat_data_root() return env ? env : "../data"; } +// Macro that skips the test if no GPU (CUDA or SYCL) is available. +#define GPRAT_SKIP_IF_NO_GPU() \ + do { \ + if (!utils::compiled_with_cuda() && !utils::compiled_with_sycl()) \ + SKIP("GPRat not compiled with GPU support"); \ + if (gprat::gpu_count() == 0) \ + SKIP("No GPU detected"); \ + } while (false) + TEST_CASE("GP: results are tile-count invariant", "[unit][gp][predict][tiling]") { // Tiling is purely a scheduling/decomposition detail: predictions, uncertainty, @@ -77,4 +86,47 @@ TEST_CASE("GP: results are tile-count invariant", "[unit][gp][predict][tiling]") } } +TEST_CASE("GP::predict_with_uncertainty: GPU matches CPU with mismatched tile sizes", "[gpu]") +{ + // Regression test: gen_tile_cross_cov_T's CUDA transpose kernel launch, and separately the + // SYCL TransposeKernel call site plus the SYCL trsm oneMath wrapper, had width/height (resp. + // M/N) arguments swapped -- which only produces wrong results for non-square tiles, i.e. + // whenever the test tile size differs from the training tile size. n_test=48 does not divide + // evenly into n_tile_size=32 here (48 % 32 != 0), so compute_test_tiles falls back to + // m_tile_size = n_test / n_tiles = 12, genuinely different from n_tile_size = 32 -- unlike + // e.g. n_test=64, which compute_test_tiles would keep at m_tile_size == n_tile_size and so + // would not have caught this. + GPRAT_SKIP_IF_NO_GPU(); + + const std::string root = gprat_data_root(); + + constexpr int n = 128, n_tiles = 4, n_reg = 8, n_test = 48; + const int tile_size = utils::compute_train_tile_size(n, n_tiles); + const auto [m_tiles, m_tile_size] = utils::compute_test_tiles(n_test, n_tiles, tile_size); + REQUIRE(m_tile_size != tile_size); + + gprat::GP_data train_in(root + "/data_1024/training_input.txt", n, n_reg); + gprat::GP_data train_out(root + "/data_1024/training_output.txt", n, 1); + gprat::GP_data test_in(root + "/data_1024/test_input.txt", n_test, n_reg); + + gprat::GP gp_cpu(train_in.data, train_out.data, n_tiles, tile_size, n_reg, { 1.0, 1.0, 0.1 }, { true, true, true }); + gprat::GP gp_gpu( + train_in.data, train_out.data, n_tiles, tile_size, n_reg, { 1.0, 1.0, 0.1 }, { true, true, true }, 0, 1); + + hpx_runtime_guard hpx_guard; + const auto cpu_unc = gp_cpu.predict_with_uncertainty(test_in.data, m_tiles, m_tile_size); + const auto gpu_unc = gp_gpu.predict_with_uncertainty(test_in.data, m_tiles, m_tile_size); + + REQUIRE(gpu_unc[0].size() == static_cast(n_test)); + REQUIRE(gpu_unc[1].size() == static_cast(n_test)); + for (int i = 0; i < n_test; ++i) + { + const auto ui = static_cast(i); + INFO("i=" << i); + REQUIRE(gpu_unc[1][ui] >= 0.0); + REQUIRE_THAT(gpu_unc[0][ui], WithinRel(cpu_unc[0][ui], 1e-4)); + REQUIRE_THAT(gpu_unc[1][ui], WithinRel(cpu_unc[1][ui], 1e-4)); + } +} + } // namespace gprat::test From 4065569a76284e477dd5da59df85143da29d37b5 Mon Sep 17 00:00:00 2001 From: constracktor <74077030+constracktor@users.noreply.github.com> Date: Thu, 30 Jul 2026 10:52:52 +0200 Subject: [PATCH 3/3] add examples/prediction_comparison and examples/optimization_comparison Small, non-benchmarking examples that run GPRat, GPflow, and GPyTorch on the same data with matched hyperparameters (each backend under its own interpreter, since they need mutually incompatible environments) and compare results directly, unlike the existing per-backend examples which are benchmarking harnesses. - prediction_comparison: identical kernel hyperparameters, no optimization -- predictions agree across backends to float64 machine precision (~1e-12 to 1e-15). Supports --gpu cuda/sycl to compare GPRat's GPU backend instead of CPU against GPflow/GPyTorch on CPU. - optimization_comparison: 300 Adam iterations from the same starting point with matched Adam settings -- fitted hyperparameters and loss trajectory agree to ~1e-4-1e-7 relative (looser than predictions since it's iterative: independent analytic vs. two different autodiff gradient implementations compound over 300 sequential steps). Widen the examples/*/include/ gitignore pattern (was gprat_*-only) to cover these new directories' cmake --install output too. --- .gitignore | 2 +- examples/optimization_comparison/compare.py | 83 +++++++ .../gpflow_optimize.py | 74 ++++++ .../optimization_comparison/gprat_optimize.py | 56 +++++ .../gpytorch_optimize.py | 76 +++++++ .../optimization_comparison/run_comparison.sh | 53 +++++ examples/prediction_comparison/compare.py | 87 ++++++++ .../prediction_comparison/gpflow_predict.py | 49 ++++ .../prediction_comparison/gprat_predict.py | 62 ++++++ .../prediction_comparison/gpytorch_predict.py | 64 ++++++ .../prediction_comparison/run_comparison.sh | 210 ++++++++++++++++++ 11 files changed, 815 insertions(+), 1 deletion(-) create mode 100644 examples/optimization_comparison/compare.py create mode 100644 examples/optimization_comparison/gpflow_optimize.py create mode 100644 examples/optimization_comparison/gprat_optimize.py create mode 100644 examples/optimization_comparison/gpytorch_optimize.py create mode 100755 examples/optimization_comparison/run_comparison.sh create mode 100644 examples/prediction_comparison/compare.py create mode 100644 examples/prediction_comparison/gpflow_predict.py create mode 100644 examples/prediction_comparison/gprat_predict.py create mode 100644 examples/prediction_comparison/gpytorch_predict.py create mode 100755 examples/prediction_comparison/run_comparison.sh diff --git a/.gitignore b/.gitignore index 2f28fb1b..216ff08f 100644 --- a/.gitignore +++ b/.gitignore @@ -192,7 +192,7 @@ cython_debug/ # Datafiles examples/*/output.csv -examples/gprat_*/include/ +examples/*/include/ examples/gpflow_reference/GPflow /.vscode benchmark_results_* diff --git a/examples/optimization_comparison/compare.py b/examples/optimization_comparison/compare.py new file mode 100644 index 00000000..f229f7c0 --- /dev/null +++ b/examples/optimization_comparison/compare.py @@ -0,0 +1,83 @@ +""" +Run gprat_optimize.py / gpflow_optimize.py / gpytorch_optimize.py, each under +its own interpreter (they depend on mutually incompatible environments), and +compare the fitted kernel hyperparameters and loss trajectory across all +three. + +All three start from the same hyperparameters and use the same Adam settings +(lr=0.1, beta1=0.9, beta2=0.999, epsilon=1e-8) for the same number of +iterations. Their optimizers are otherwise independent implementations +(GPRat: hand-written C++ Adam over analytic gradients; GPflow/GPyTorch: +autodiff + framework Adam), so agreement here demonstrates convergence to the +same optimum, not shared code. +""" + +import json +import subprocess +import sys +from pathlib import Path + +SCRIPT_DIR = Path(__file__).resolve().parent + +BACKENDS = { + "gprat": [sys.executable, str(SCRIPT_DIR / "gprat_optimize.py")], + "gpflow": [ + str(SCRIPT_DIR / "../gpflow_reference/gpflow_cpu_env/bin/python"), + str(SCRIPT_DIR / "gpflow_optimize.py"), + ], + "gpytorch": [ + str(SCRIPT_DIR / "../gpytorch_reference/gpytorch_cpu_env/bin/python"), + str(SCRIPT_DIR / "gpytorch_optimize.py"), + ], +} + +MARKER = "RESULT_JSON:" +CHECKPOINTS = [0, 1, 2, 5, 10, 25, 50, 100, 150, 200, 250, -1] + + +def run_backend(name, cmd): + result = subprocess.run(cmd, capture_output=True, text=True, cwd=SCRIPT_DIR) + for line in result.stdout.splitlines(): + if line.startswith(MARKER): + return json.loads(line[len(MARKER) :]) + raise RuntimeError( + f"{name}: no {MARKER} line found in output.\nstdout:\n{result.stdout}\nstderr:\n{result.stderr}" + ) + + +def main(): + results = {name: run_backend(name, cmd) for name, cmd in BACKENDS.items()} + + n_iter = len(results["gprat"]["losses"]) + for name, r in results.items(): + assert len(r["losses"]) == n_iter, f"{name}: expected {n_iter} loss values, got {len(r['losses'])}" + + print("Loss trajectory (per-sample-normalized negative log marginal likelihood):") + print(f" {'iter':>5} {'gprat':>14} {'gpflow':>14} {'gpytorch':>14}") + for i in CHECKPOINTS: + idx = i if i >= 0 else n_iter + i + print( + f" {idx:>5} " + f"{results['gprat']['losses'][i]:>14.6f} " + f"{results['gpflow']['losses'][i]:>14.6f} " + f"{results['gpytorch']['losses'][i]:>14.6f}" + ) + + print("\nFitted hyperparameters:") + rtol = 1e-2 + ok = True + for field in ("lengthscale", "variance", "noise"): + print(f"{field}:") + for a, b in [("gprat", "gpflow"), ("gprat", "gpytorch"), ("gpflow", "gpytorch")]: + va, vb = results[a][field], results[b][field] + rel_diff = abs(va - vb) / max(abs(va), abs(vb), 1e-12) + print(f" {a:9s} vs {b:9s}: {va:.6f} vs {vb:.6f} (rel_diff={rel_diff:.2e})") + ok &= rel_diff < rtol + + verdict = "PASS: fitted hyperparameters agree" if ok else "FAIL: fitted hyperparameters do not agree" + print(f"\n{verdict} within rtol={rtol:.0e}") + sys.exit(0 if ok else 1) + + +if __name__ == "__main__": + main() diff --git a/examples/optimization_comparison/gpflow_optimize.py b/examples/optimization_comparison/gpflow_optimize.py new file mode 100644 index 00000000..cd270fd4 --- /dev/null +++ b/examples/optimization_comparison/gpflow_optimize.py @@ -0,0 +1,74 @@ +""" +Optimize kernel hyperparameters with GPflow and print the result (fitted +hyperparameters + per-iteration loss) as one JSON line prefixed with +RESULT_JSON:, so compare.py can pick it out of the surrounding TensorFlow log +output. +""" + +import json +import os +from pathlib import Path + +os.environ["CUDA_VISIBLE_DEVICES"] = "" +os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3" + +import numpy as np +import tensorflow as tf +import gpflow + +SCRIPT_DIR = Path(__file__).resolve().parent +DATA_DIR = SCRIPT_DIR / "../../data/data_1024" +TRAIN_SIZE = 512 +N_REG = 8 +OPT_ITER = 300 + +gpflow.config.set_default_float(np.float64) + + +def generate_regressor(x_original, n_regressors): + x_padded = np.pad(x_original, pad_width=(n_regressors - 1, 0), mode="constant") + return np.array([x_padded[i : i + n_regressors] for i in range(len(x_original))]) + + +x_train_in = np.loadtxt(DATA_DIR / "training_input.txt", dtype="d")[:TRAIN_SIZE] +X_train = generate_regressor(x_train_in, N_REG).astype("d") +Y_train = np.loadtxt(DATA_DIR / "training_output.txt", dtype="d")[:TRAIN_SIZE, None] + +model = gpflow.models.GPR( + (X_train, Y_train), + kernel=gpflow.kernels.SquaredExponential(variance=1.0, lengthscales=1.0), + noise_variance=0.1, +) + +opt = tf.keras.optimizers.Adam(learning_rate=0.1, beta_1=0.9, beta_2=0.999, epsilon=1e-08) + + +@tf.function +def optimization_step(): + with tf.GradientTape() as tape: + loss = model.training_loss() + gradients = tape.gradient(loss, model.trainable_variables) + opt.apply_gradients(zip(gradients, model.trainable_variables)) + return loss + + +losses = [] +for _ in range(OPT_ITER): + # Divide by N: GPflow's training_loss() is the raw, unnormalized log + # marginal likelihood, whereas GPRat/GPyTorch both normalize by the + # number of training points. This only rescales the loss for reporting; + # it doesn't change the optimization trajectory, since Adam's update is + # invariant to a constant multiplicative rescaling of the gradient. + losses.append(float(optimization_step().numpy()) / TRAIN_SIZE) + +print( + "RESULT_JSON:" + + json.dumps( + { + "lengthscale": float(model.kernel.lengthscales.numpy()), + "variance": float(model.kernel.variance.numpy()), + "noise": float(model.likelihood.variance.numpy()), + "losses": losses, + } + ) +) diff --git a/examples/optimization_comparison/gprat_optimize.py b/examples/optimization_comparison/gprat_optimize.py new file mode 100644 index 00000000..8b32fa08 --- /dev/null +++ b/examples/optimization_comparison/gprat_optimize.py @@ -0,0 +1,56 @@ +""" +Optimize kernel hyperparameters with GPRat and print the result (fitted +hyperparameters + per-iteration loss) as one JSON line prefixed with +RESULT_JSON:, so compare.py can pick it out of the surrounding HPX/APEX +startup output. +""" + +import json +import re +import sys +from pathlib import Path + +SCRIPT_DIR = Path(__file__).resolve().parent +sys.path.insert(0, str(SCRIPT_DIR / "lib")) +import gprat # noqa: E402 + +DATA_DIR = SCRIPT_DIR / "../../data/data_1024" +TRAIN_SIZE = 512 +N_REG = 8 +N_TILES = 4 +OPT_ITER = 300 +KERNEL_PARAMS = [1.0, 1.0, 0.1] # lengthscale, vertical_lengthscale, noise_variance + +train_in = gprat.GP_data(str(DATA_DIR / "training_input.txt"), TRAIN_SIZE, N_REG) +# n_reg=1 -> no offset: unlike the input, the training output has no +# lookahead padding requirement. +train_out = gprat.GP_data(str(DATA_DIR / "training_output.txt"), TRAIN_SIZE, 1) + +n_tile_size = gprat.compute_train_tile_size(TRAIN_SIZE, N_TILES) + +gp = gprat.GP( + train_in.data, + train_out.data, + N_TILES, + n_tile_size, + kernel_params=KERNEL_PARAMS, + n_reg=N_REG, + trainable=[True, True, True], +) + +gprat.start_hpx(sys.argv, 2) +hpar = gprat.AdamParams(learning_rate=0.1, beta1=0.9, beta2=0.999, epsilon=1e-8, opt_iter=OPT_ITER) +losses = gp.optimize(hpar) +gprat.stop_hpx() + +# kernel_params is a bound C++ struct (SEKParams), not convertible to a +# Python value directly, so pull the fitted values back out of __repr__. +m = re.search( + r"lengthscale=([\d.eE+-]+), vertical_lengthscale=([\d.eE+-]+), noise_variance=([\d.eE+-]+)", repr(gp) +) +lengthscale, variance, noise = (float(g) for g in m.groups()) + +print( + "RESULT_JSON:" + + json.dumps({"lengthscale": lengthscale, "variance": variance, "noise": noise, "losses": list(losses)}) +) diff --git a/examples/optimization_comparison/gpytorch_optimize.py b/examples/optimization_comparison/gpytorch_optimize.py new file mode 100644 index 00000000..1be3c319 --- /dev/null +++ b/examples/optimization_comparison/gpytorch_optimize.py @@ -0,0 +1,76 @@ +""" +Optimize kernel hyperparameters with GPyTorch and print the result (fitted +hyperparameters + per-iteration loss) as one JSON line prefixed with +RESULT_JSON:, so compare.py can pick it out of the surrounding output. +""" + +import json +from pathlib import Path + +import numpy as np +import torch +import gpytorch + +SCRIPT_DIR = Path(__file__).resolve().parent +DATA_DIR = SCRIPT_DIR / "../../data/data_1024" +TRAIN_SIZE = 512 +N_REG = 8 +OPT_ITER = 300 + +torch.set_default_dtype(torch.float64) + + +def generate_regressor(x_original, n_regressors): + x_padded = np.pad(x_original, pad_width=(n_regressors - 1, 0), mode="constant") + return np.array([x_padded[i : i + n_regressors] for i in range(len(x_original))], dtype="d") + + +class ExactGPModel(gpytorch.models.ExactGP): + def __init__(self, train_x, train_y, likelihood): + super().__init__(train_x, train_y, likelihood) + # Zero mean to match GPRat/GPflow, which assume a zero-mean prior. + self.mean_module = gpytorch.means.ZeroMean() + self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel()) + self.covar_module.base_kernel.lengthscale = 1.0 + self.covar_module.outputscale = 1.0 + + def forward(self, x): + return gpytorch.distributions.MultivariateNormal(self.mean_module(x), self.covar_module(x)) + + +x_train_in = np.loadtxt(DATA_DIR / "training_input.txt", dtype="d")[:TRAIN_SIZE] +X_train = torch.from_numpy(generate_regressor(x_train_in, N_REG)) +Y_train = torch.from_numpy(np.loadtxt(DATA_DIR / "training_output.txt", dtype="d")[:TRAIN_SIZE]) + +likelihood = gpytorch.likelihoods.GaussianLikelihood() +likelihood.noise = 0.1 +model = ExactGPModel(X_train, Y_train, likelihood) + +model.train() +likelihood.train() + +optimizer = torch.optim.Adam(model.parameters(), lr=0.1, betas=(0.9, 0.999), eps=1e-8) +# ExactMarginalLogLikelihood already divides by the number of training +# points, matching GPRat's per-sample-normalized loss convention. +mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model) + +losses = [] +for _ in range(OPT_ITER): + optimizer.zero_grad() + output = model(X_train) + loss = -mll(output, Y_train) + loss.backward() + losses.append(loss.item()) + optimizer.step() + +print( + "RESULT_JSON:" + + json.dumps( + { + "lengthscale": model.covar_module.base_kernel.lengthscale.item(), + "variance": model.covar_module.outputscale.item(), + "noise": model.likelihood.noise.item(), + "losses": losses, + } + ) +) diff --git a/examples/optimization_comparison/run_comparison.sh b/examples/optimization_comparison/run_comparison.sh new file mode 100755 index 00000000..eb105112 --- /dev/null +++ b/examples/optimization_comparison/run_comparison.sh @@ -0,0 +1,53 @@ +#!/bin/bash +# Runs Adam hyperparameter optimization with GPRat, GPflow, and GPyTorch from +# the same starting point with matched Adam settings, and compares the +# fitted hyperparameters and loss trajectory. CPU only. +# +# Assumes examples/gpflow_reference/gpflow_cpu_env and +# examples/gpytorch_reference/gpytorch_cpu_env already exist (see +# run_gpflow.sh cpu / run_gpytorch.sh cpu) and that GPRat has been built via +# the release-linux CMake preset. + +SCRIPT_DIR=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) +source "$SCRIPT_DIR/../../site_paths.sh" +cd "$SCRIPT_DIR" + +HOSTNAME=$(hostname -s) + +### ENVIRONMENT SETUP ############################################################################# + +if [[ "$HOSTNAME" == "sven0" || "$HOSTNAME" == "sven1" ]]; then + + export LD_LIBRARY_PATH=$HOME/git_workspace/build-scripts/build/hpx/lib64:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$HOME/git_workspace/build-scripts/build/boost/lib:$LD_LIBRARY_PATH + export LD_PRELOAD=$HOME/git_workspace/build-scripts/build/jemalloc/lib/libjemalloc.so.2 + +elif [[ "$HOSTNAME" == "simcl1n1" || "$HOSTNAME" == "simcl1n2" || \ + "$HOSTNAME" == "simcl1n3" || "$HOSTNAME" == "simcl1n4" ]]; then + + source "$SIMCL1_SPACK_ROOT/spack/share/spack/setup-env.sh" + + if spack env list | grep -q "gprat_cpu_gcc"; then + echo "Found gprat_cpu_gcc environment, activating it." + spack env activate gprat_cpu_gcc + module load gcc/14.1.0 + export LD_LIBRARY_PATH=$(spack location -i hpx)/lib:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$(spack location -i openblas)/lib:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$(spack location -i intel-oneapi-mkl)/lib:$LD_LIBRARY_PATH + fi + +elif [[ "$HOSTNAME" == "pcsgs04" ]]; then + + source "$PCSGS04_SPACK_ROOT/share/spack/setup-env.sh" + +fi + +### INSTALL MATCHING GPRAT BUILD ################################################################## + +GPRAT_ROOT="$SCRIPT_DIR/../.." +cmake --install "$GPRAT_ROOT/build/release-linux" --prefix "$SCRIPT_DIR" +cp "$GPRAT_ROOT/build/release-linux/bindings"/gprat.cpython-*.so "$SCRIPT_DIR/lib/" + +### EXECUTION ##################################################################################### + +python3 compare.py diff --git a/examples/prediction_comparison/compare.py b/examples/prediction_comparison/compare.py new file mode 100644 index 00000000..74bad9e4 --- /dev/null +++ b/examples/prediction_comparison/compare.py @@ -0,0 +1,87 @@ +""" +Run gprat_predict.py / gpflow_predict.py / gpytorch_predict.py, each under its +own interpreter (they depend on mutually incompatible environments), and +compare the predicted mean and variance across all backends. + +All backends are configured with the same kernel hyperparameters and no +optimization, so their predictions should agree to floating-point precision. + +Pass --gpu cuda or --gpu sycl to run GPRat on GPU instead of CPU (matching +whichever GPU build run_comparison.sh installed into lib/) and compare it +against GPflow/GPyTorch on CPU. Only one GPRat variant runs per invocation -- +see run_comparison.sh's comment for why CPU and GPU builds are never loaded +in the same process tree. +""" + +import argparse +import itertools +import json +import subprocess +import sys +from pathlib import Path + +SCRIPT_DIR = Path(__file__).resolve().parent + +parser = argparse.ArgumentParser() +parser.add_argument("--gpu", choices=["cuda", "sycl"], help="Run GPRat on GPU instead of CPU") +args = parser.parse_args() + +gprat_cmd = [sys.executable, str(SCRIPT_DIR / "gprat_predict.py")] +gprat_key = "gprat" +if args.gpu: + gprat_cmd.append("--use-gpu") + gprat_key = "gprat_gpu" + +BACKENDS = { + gprat_key: gprat_cmd, + "gpflow": [ + str(SCRIPT_DIR / "../gpflow_reference/gpflow_cpu_env/bin/python"), + str(SCRIPT_DIR / "gpflow_predict.py"), + ], + "gpytorch": [ + str(SCRIPT_DIR / "../gpytorch_reference/gpytorch_cpu_env/bin/python"), + str(SCRIPT_DIR / "gpytorch_predict.py"), + ], +} + +MARKER = "RESULT_JSON:" + + +def run_backend(name, cmd): + result = subprocess.run(cmd, capture_output=True, text=True, cwd=SCRIPT_DIR) + for line in result.stdout.splitlines(): + if line.startswith(MARKER): + return json.loads(line[len(MARKER) :]) + raise RuntimeError( + f"{name}: no {MARKER} line found in output.\nstdout:\n{result.stdout}\nstderr:\n{result.stderr}" + ) + + +def max_diff(a, b): + assert len(a) == len(b), f"length mismatch: {len(a)} vs {len(b)}" + abs_diff = max(abs(x - y) for x, y in zip(a, b)) + rel_diff = max(abs(x - y) / max(abs(x), abs(y), 1e-12) for x, y in zip(a, b)) + return abs_diff, rel_diff + + +def main(): + predictions = {name: run_backend(name, cmd) for name, cmd in BACKENDS.items()} + + print(f"Test size: {len(predictions[gprat_key]['mean'])}\n") + + rtol = 1e-6 + ok = True + for field in ("mean", "var"): + print(f"{field}:") + for a, b in itertools.combinations(BACKENDS, 2): + abs_diff, rel_diff = max_diff(predictions[a][field], predictions[b][field]) + print(f" {a:9s} vs {b:9s}: max_abs_diff={abs_diff:.3e} max_rel_diff={rel_diff:.3e}") + ok &= rel_diff < rtol + + verdict = "PASS: predictions agree" if ok else "FAIL: predictions do not agree" + print(f"\n{verdict} within rtol={rtol:.0e}") + sys.exit(0 if ok else 1) + + +if __name__ == "__main__": + main() diff --git a/examples/prediction_comparison/gpflow_predict.py b/examples/prediction_comparison/gpflow_predict.py new file mode 100644 index 00000000..2f5e6936 --- /dev/null +++ b/examples/prediction_comparison/gpflow_predict.py @@ -0,0 +1,49 @@ +""" +Predict with GPflow and print the result as one JSON line prefixed with +RESULT_JSON:, so compare.py can pick it out of the surrounding +TensorFlow log output. +""" + +import json +import os +from pathlib import Path + +os.environ["CUDA_VISIBLE_DEVICES"] = "" +os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3" + +import numpy as np +import gpflow + +SCRIPT_DIR = Path(__file__).resolve().parent +DATA_DIR = SCRIPT_DIR / "../../data/data_1024" +TRAIN_SIZE = 512 +TEST_SIZE = 64 +N_REG = 8 + +gpflow.config.set_default_float(np.float64) + + +def generate_regressor(x_original, n_regressors): + x_padded = np.pad(x_original, pad_width=(n_regressors - 1, 0), mode="constant") + return np.array([x_padded[i : i + n_regressors] for i in range(len(x_original))]) + + +x_train_in = np.loadtxt(DATA_DIR / "training_input.txt", dtype="d")[:TRAIN_SIZE] +x_test_in = np.loadtxt(DATA_DIR / "test_input.txt", dtype="d")[:TEST_SIZE] + +X_train = generate_regressor(x_train_in, N_REG).astype("d") +X_test = generate_regressor(x_test_in, N_REG).astype("d") +Y_train = np.loadtxt(DATA_DIR / "training_output.txt", dtype="d")[:TRAIN_SIZE, None] + +model = gpflow.models.GPR( + (X_train, Y_train), + kernel=gpflow.kernels.SquaredExponential(variance=1.0, lengthscales=1.0), + noise_variance=0.1, +) + +mean, var = model.predict_f(X_test) + +print( + "RESULT_JSON:" + + json.dumps({"mean": mean.numpy().flatten().tolist(), "var": var.numpy().flatten().tolist()}) +) diff --git a/examples/prediction_comparison/gprat_predict.py b/examples/prediction_comparison/gprat_predict.py new file mode 100644 index 00000000..516c448e --- /dev/null +++ b/examples/prediction_comparison/gprat_predict.py @@ -0,0 +1,62 @@ +""" +Predict with GPRat and print the result as one JSON line prefixed with +RESULT_JSON:, so compare.py can pick it out of the surrounding HPX/APEX +startup output. +""" + +import argparse +import json +import sys +from pathlib import Path + +parser = argparse.ArgumentParser() +parser.add_argument("--use-gpu", action="store_true", help="Run on GPU (CUDA or SYCL) instead of CPU") +args = parser.parse_args() +sys.argv = [sys.argv[0]] # strip our own flags before they reach gprat.start_hpx below + +SCRIPT_DIR = Path(__file__).resolve().parent +sys.path.insert(0, str(SCRIPT_DIR / "lib")) +import gprat # noqa: E402 + +DATA_DIR = SCRIPT_DIR / "../../data/data_1024" +TRAIN_SIZE = 512 +TEST_SIZE = 64 +N_REG = 8 +N_TILES = 4 +KERNEL_PARAMS = [1.0, 1.0, 0.1] # lengthscale, vertical_lengthscale, noise_variance + +if args.use_gpu: + if not (gprat.compiled_with_cuda() or gprat.compiled_with_sycl()): + print("gprat was not compiled with GPU support but --use-gpu was passed.", file=sys.stderr) + sys.exit(1) + if gprat.gpu_count() == 0: + print("gprat was compiled with GPU support but no GPU was found.", file=sys.stderr) + sys.exit(1) + +train_in = gprat.GP_data(str(DATA_DIR / "training_input.txt"), TRAIN_SIZE, N_REG) +# n_reg=1 -> no offset: unlike the input, the training output has no +# lookahead padding requirement. +train_out = gprat.GP_data(str(DATA_DIR / "training_output.txt"), TRAIN_SIZE, 1) +test_in = gprat.GP_data(str(DATA_DIR / "test_input.txt"), TEST_SIZE, N_REG) + +n_tile_size = gprat.compute_train_tile_size(TRAIN_SIZE, N_TILES) +m_tiles, m_tile_size = gprat.compute_test_tiles(TEST_SIZE, N_TILES, n_tile_size) + +gp_kwargs = dict( + kernel_params=KERNEL_PARAMS, + n_reg=N_REG, + # Predict with the given kernel_params as-is, no optimization: this keeps + # the hyperparameters identical across all backends so predictions are + # directly comparable. + trainable=[False, False, False], +) +if args.use_gpu: + gp_kwargs.update(gpu_id=0, n_units=1) + +gp = gprat.GP(train_in.data, train_out.data, N_TILES, n_tile_size, **gp_kwargs) + +gprat.start_hpx(sys.argv, 2) +mean, var = gp.predict_with_uncertainty(test_in.data, m_tiles, m_tile_size) +gprat.stop_hpx() + +print("RESULT_JSON:" + json.dumps({"mean": list(mean), "var": list(var)})) diff --git a/examples/prediction_comparison/gpytorch_predict.py b/examples/prediction_comparison/gpytorch_predict.py new file mode 100644 index 00000000..ae55a0ea --- /dev/null +++ b/examples/prediction_comparison/gpytorch_predict.py @@ -0,0 +1,64 @@ +""" +Predict with GPyTorch and print the result as one JSON line prefixed with +RESULT_JSON:, so compare.py can pick it out of the surrounding output. +""" + +import json +from pathlib import Path + +import numpy as np +import torch +import gpytorch + +SCRIPT_DIR = Path(__file__).resolve().parent +DATA_DIR = SCRIPT_DIR / "../../data/data_1024" +TRAIN_SIZE = 512 +TEST_SIZE = 64 +N_REG = 8 + +torch.set_default_dtype(torch.float64) + + +def generate_regressor(x_original, n_regressors): + x_padded = np.pad(x_original, pad_width=(n_regressors - 1, 0), mode="constant") + return np.array([x_padded[i : i + n_regressors] for i in range(len(x_original))], dtype="d") + + +class ExactGPModel(gpytorch.models.ExactGP): + def __init__(self, train_x, train_y, likelihood): + super().__init__(train_x, train_y, likelihood) + # Zero mean to match GPRat/GPflow, which assume a zero-mean prior. + self.mean_module = gpytorch.means.ZeroMean() + self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel()) + self.covar_module.base_kernel.lengthscale = 1.0 + self.covar_module.outputscale = 1.0 + + def forward(self, x): + return gpytorch.distributions.MultivariateNormal(self.mean_module(x), self.covar_module(x)) + + +x_train_in = np.loadtxt(DATA_DIR / "training_input.txt", dtype="d")[:TRAIN_SIZE] +x_test_in = np.loadtxt(DATA_DIR / "test_input.txt", dtype="d")[:TEST_SIZE] + +X_train = torch.from_numpy(generate_regressor(x_train_in, N_REG)) +X_test = torch.from_numpy(generate_regressor(x_test_in, N_REG)) +Y_train = torch.from_numpy(np.loadtxt(DATA_DIR / "training_output.txt", dtype="d")[:TRAIN_SIZE]) + +likelihood = gpytorch.likelihoods.GaussianLikelihood() +likelihood.noise = 0.1 +model = ExactGPModel(X_train, Y_train, likelihood) + +model.eval() +likelihood.eval() + +with ( + torch.no_grad(), + gpytorch.settings.fast_pred_var(False), + gpytorch.settings.lazily_evaluate_kernels(False), + gpytorch.settings.fast_computations(covar_root_decomposition=False, log_prob=False, solves=False), +): + f_pred = model(X_test) + mean = f_pred.mean + var = f_pred.variance + +print("RESULT_JSON:" + json.dumps({"mean": mean.tolist(), "var": var.tolist()})) diff --git a/examples/prediction_comparison/run_comparison.sh b/examples/prediction_comparison/run_comparison.sh new file mode 100755 index 00000000..e99088c2 --- /dev/null +++ b/examples/prediction_comparison/run_comparison.sh @@ -0,0 +1,210 @@ +#!/bin/bash +# Input $1 (optional): cpu (default) / cuda / sycl -- which GPRat build to compare. +# Input $2: If $1 is sycl: nvidia/amd/intel. +# +# Runs GPRat, GPflow, and GPyTorch on the same data/hyperparameters (no +# optimization) and compares their predictions. +# +# GPflow/GPyTorch always run on CPU (via their own venvs); only GPRat's +# device varies. GPRat's CPU and GPU builds are never loaded in the same +# process tree -- each invocation of this script activates exactly one +# environment and installs exactly one build into lib/, matching how +# examples/gprat_python/run_gprat_python.sh runs GPRat (mixing two HPX/BLAS +# builds' libraries on LD_LIBRARY_PATH in the same process is untested and +# not worth risking). +# +# Assumes examples/gpflow_reference/gpflow_cpu_env and +# examples/gpytorch_reference/gpytorch_cpu_env already exist (see +# run_gpflow.sh cpu / run_gpytorch.sh cpu) and that GPRat has been built via +# the matching CMake preset (release-linux / release-linux-cuda / +# release-linux-sycl). + +SCRIPT_DIR=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) +source "$SCRIPT_DIR/../../site_paths.sh" +cd "$SCRIPT_DIR" + +HOSTNAME=$(hostname -s) + +if [[ -z "$1" ]]; then + echo "Input parameter is missing. Using default: Run computations on CPU" + GPU="" +elif [[ "$1" == "cuda" || "$1" == "sycl" ]]; then + GPU="--gpu $1" + if [[ \ + "$HOSTNAME" != "simcl1n1" && \ + "$HOSTNAME" != "simcl1n2" && \ + "$HOSTNAME" != "simcl1n3" && \ + "$HOSTNAME" != "simcl1n4" && \ + "$HOSTNAME" != "pcsgs04" ]]; + then + echo "GPU execution with this script is only supported on simcl1n1, simcl1n2, simcl1n3, simcl1n4, and pcsgs04." 1>&2 + exit 1 + fi +elif [[ "$1" != "cpu" ]]; then + echo "Please specify input parameter: cpu/cuda/sycl" + exit 1 +fi + +### ENVIRONMENT SETUP ############################################################################# + +if [[ "$HOSTNAME" == "sven0" || "$HOSTNAME" == "sven1" ]]; then + + export LD_LIBRARY_PATH=$HOME/git_workspace/build-scripts/build/hpx/lib64:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$HOME/git_workspace/build-scripts/build/boost/lib:$LD_LIBRARY_PATH + export LD_PRELOAD=$HOME/git_workspace/build-scripts/build/jemalloc/lib/libjemalloc.so.2 + +elif [[ "$HOSTNAME" == "simcl1n1" || "$HOSTNAME" == "simcl1n2" || \ + "$HOSTNAME" == "simcl1n3" || "$HOSTNAME" == "simcl1n4" ]]; then + + source "$SIMCL1_SPACK_ROOT/spack/share/spack/setup-env.sh" + + if [[ "$1" == "cuda" || "$1" == "sycl" ]]; then + + # simcl1n4 does not have a GPU + if [[ "$HOSTNAME" == "simcl1n4" ]]; then + echo "Machine $HOSTNAME does not have a GPU but you selected GPU execution." 1>&2 + exit 1 + fi + + if spack env list | grep -q "gprat_gpu_clang"; then + echo "Found gprat_gpu_clang environment, activating it." + spack env activate gprat_gpu_clang + export LD_LIBRARY_PATH=$(spack location -i hpx)/lib:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$(spack location -i openblas)/lib:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$(spack location -i intel-oneapi-mkl)/lib:$LD_LIBRARY_PATH + fi + + if [[ "$1" == "cuda" || ( "$1" == "sycl" && "$2" == "nvidia" ) ]]; then + module load cuda/12.0.1 + module load clang/17.0.1 + fi + + if [[ "$1" == "sycl" ]]; then + + if [[ "$2" == "nvidia" ]]; then + + ONEMATH_PATH="${ONEMATH_NVIDIA_ROOT}/lib/" + export LD_LIBRARY_PATH="$ONEMATH_PATH:$LD_LIBRARY_PATH" + + elif [[ "$2" == "amd" ]]; then + + ONEMATH_PATH="${ONEMATH_AMD_ROOT}/lib/" + export LD_LIBRARY_PATH="$ONEMATH_PATH:$LD_LIBRARY_PATH" + + ROCM_PATH=${ROCM_PATH:-/opt/rocm-6.4.0} + if [[ -d "$ROCM_PATH" ]]; then + export LD_LIBRARY_PATH="$ROCM_PATH/lib:$ROCM_PATH/lib64:$ROCM_PATH/hip/lib:$LD_LIBRARY_PATH" + export ROCM_PATH + fi + + COMGR_COMPAT_DIR="/data/scratch-simcl1/breyerml/Programs/.modulefiles/icpx" + if [[ -d "$COMGR_COMPAT_DIR" ]]; then + export LD_LIBRARY_PATH="$COMGR_COMPAT_DIR:$LD_LIBRARY_PATH" + fi + + ONEAPI_SETVARS="/import/sgs.scratch-simcl1/breyerml/Programs/spack/opt/spack/linux-zen4/intel-oneapi-compilers-2025.1.1-5ynklzzqslh265azbglzqdtecdghl7ob/setvars.sh" + if ! command -v icpx &>/dev/null && [[ -f "$ONEAPI_SETVARS" ]]; then + ONEAPI_COMPILER_ROOT="$(dirname $ONEAPI_SETVARS)/compiler/2025.1" + export PATH="$ONEAPI_COMPILER_ROOT/bin:$PATH" + export LD_LIBRARY_PATH="$ONEAPI_COMPILER_ROOT/lib:$LD_LIBRARY_PATH" + elif command -v icpx &>/dev/null; then + ONEAPI_COMPILER_ROOT="$(dirname $(dirname $(which icpx)))" + export LD_LIBRARY_PATH="$ONEAPI_COMPILER_ROOT/lib:$LD_LIBRARY_PATH" + fi + + export HSA_XNACK=1 + + elif [[ "$2" == "intel" ]]; then + + echo "Machine $HOSTNAME does not have an Intel GPU." 1>&2 + exit 1 + + elif [[ "$2" != "nvidia" ]]; then + + echo "Please specify gpu vendor: nvidia/amd/intel" + exit 1 + + fi + + fi + + elif [[ "$1" == "cpu" || -z "$1" ]]; then + + if spack env list | grep -q "gprat_cpu_gcc"; then + echo "Found gprat_cpu_gcc environment, activating it." + spack env activate gprat_cpu_gcc + module load gcc/14.1.0 + export LD_LIBRARY_PATH=$(spack location -i hpx)/lib:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$(spack location -i openblas)/lib:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=$(spack location -i intel-oneapi-mkl)/lib:$LD_LIBRARY_PATH + fi + + fi + +elif [[ "$HOSTNAME" == "pcsgs04" ]]; then + + source "$PCSGS04_SPACK_ROOT/share/spack/setup-env.sh" + + if [[ "$1" == "cuda" || "$1" == "sycl" ]]; then + + if spack env list | grep -q "gprat_gpu_clang"; then + + echo "Found gprat_gpu_clang environment, activating it." + spack env activate gprat_gpu_clang + export LD_LIBRARY_PATH=$(spack location -i hpx)/lib:$LD_LIBRARY_PATH + + if [[ "$1" == "sycl" ]]; then + + if [[ "$2" != "intel" ]]; then + echo "pcsgs04 only has an Intel GPU. Please specify gpu vendor: intel" 1>&2 + exit 1 + fi + + if ! command -v icpx &>/dev/null && [[ -f /opt/intel/oneapi/compiler/2025.3/env/vars.sh ]]; then + source /opt/intel/oneapi/compiler/2025.3/env/vars.sh + fi + + if [[ -f /opt/intel/oneapi/umf/latest/env/vars.sh ]]; then + source /opt/intel/oneapi/umf/latest/env/vars.sh + fi + + if [[ -f /opt/intel/oneapi/mkl/2025.3/env/vars.sh ]]; then + source /opt/intel/oneapi/mkl/2025.3/env/vars.sh + fi + export LD_LIBRARY_PATH="/opt/intel/oneapi/tbb/2022.3/lib/intel64/gcc4.8:$LD_LIBRARY_PATH" + + ONEMATH_PATH="${ONEMATH_INTEL_ROOT}/lib" + export LD_LIBRARY_PATH="$ONEMATH_PATH:$LD_LIBRARY_PATH" + + fi + + else + + echo \ + "Cannot find Spack environment gprat_gpu_clang. Please run spack-repo/environments/setup_gprat_gpu_clang.sh" 1>&2 + exit 1 + + fi + + fi + +fi + +### INSTALL MATCHING GPRAT BUILD ################################################################## + +GPRAT_ROOT="$SCRIPT_DIR/../.." + +if [[ "$1" == "cuda" ]]; then + GPRAT_BUILD_DIR="$GPRAT_ROOT/build/release-linux-cuda" +elif [[ "$1" == "sycl" ]]; then + GPRAT_BUILD_DIR="$GPRAT_ROOT/build/release-linux-sycl" +else + GPRAT_BUILD_DIR="$GPRAT_ROOT/build/release-linux" +fi + +cmake --install "$GPRAT_BUILD_DIR" --prefix "$SCRIPT_DIR" +cp "$GPRAT_BUILD_DIR"/bindings/gprat.cpython-*.so "$SCRIPT_DIR/lib/" + +### EXECUTION ##################################################################################### + +python3 compare.py $GPU