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/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/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 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..b362d158 --- /dev/null +++ b/test/src/unit_tests.cpp @@ -0,0 +1,132 @@ +#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"; +} + +// 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, + // 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)); + } +} + +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