Examples and tests

.gitignore

@@ -192,7 +192,7 @@ cython_debug/
 
 # Datafiles
 examples/*/output.csv
-examples/gprat_*/include/
+examples/*/include/
 examples/gpflow_reference/GPflow
 /.vscode
 benchmark_results_*

core/src/gpu/cuda/gp_algorithms.cu

@@ -276,7 +276,7 @@ hpx::shared_future<double *> 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<<<n_blocks, threads_per_block, 0, stream>>>(
-        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));
 

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<std::int64_t>(M),
         static_cast<std::int64_t>(N),
+        static_cast<std::int64_t>(M),
         alpha,
         f_A,
         static_cast<std::int64_t>(M),

core/src/gpu/sycl/sycl_gp_algorithms.cpp

@@ -186,7 +186,7 @@ hpx::shared_future<double *> 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);
             });
 

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()

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,
+        }
+    )
+)

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)})
+)

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,
+        }
+    )
+)

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