diff --git a/benchmark/scripts/collect_results.py b/benchmark/scripts/collect_results.py
new file mode 100644
index 00000000000..c15ce63aa84
--- /dev/null
+++ b/benchmark/scripts/collect_results.py
@@ -0,0 +1,41 @@
+"""
+Parse per-run accuracy files into a single results.csv.
+
+Output columns: label, epsilon, private, accuracy
+"""
+import pathlib, csv, re
+
+RESULTS = pathlib.Path("benchmark/results")
+
+rows = []
+
+def parse_acc(path: pathlib.Path) -> float:
+ txt = path.read_text().strip()
+ # SystemDS writes a bare float.
+ return float(txt)
+
+# Non-private baseline.
+baseline_path = RESULTS / "acc_baseline.txt"
+if baseline_path.exists():
+ rows.append(dict(label="baseline", epsilon="inf",
+ private=0, accuracy=parse_acc(baseline_path)))
+
+# DP runs.
+for eps in [0.5, 1, 4, 8]:
+ p = RESULTS / f"acc_eps_{eps}.txt"
+ if p.exists():
+ rows.append(dict(label=f"ε={eps}", epsilon=eps,
+ private=1, accuracy=parse_acc(p)))
+ else:
+ print(f"Warning: {p} not found — skipping")
+
+out = RESULTS / "results.csv"
+with open(out, "w", newline="") as f:
+ w = csv.DictWriter(f, fieldnames=["label","epsilon","private","accuracy"])
+ w.writeheader()
+ w.writerows(rows)
+
+print(f"Wrote {out}")
+for r in rows:
+ print(f" {r['label']:12s} acc={r['accuracy']:.4f}")
+
diff --git a/benchmark/scripts/eval.dml b/benchmark/scripts/eval.dml
new file mode 100644
index 00000000000..7f0c89fbaa4
--- /dev/null
+++ b/benchmark/scripts/eval.dml
@@ -0,0 +1,22 @@
+# eval.dml — compute binary classification accuracy on held-out test set.
+# Arguments: data_dir, model_path, out_acc
+data_dir = $data_dir;
+model_path = $model_path;
+out_acc = $out_acc;
+
+X_test = read(data_dir + "/X_test.csv",
+ data_type="matrix", value_type="double", format="csv");
+y_test = read(data_dir + "/y_test.csv",
+ data_type="matrix", value_type="double", format="csv");
+w = read(model_path,
+ data_type="matrix", value_type="double", format="csv");
+
+scores = X_test %*% w;
+preds = (scores > 0.0); # threshold at 0 (log-odds)
+correct = sum(preds == y_test);
+n_test = nrow(y_test);
+accuracy = correct / n_test;
+
+print("Accuracy: " + accuracy);
+write(accuracy, out_acc, format="csv");
+
diff --git a/benchmark/scripts/fedavg_dp.dml b/benchmark/scripts/fedavg_dp.dml
new file mode 100644
index 00000000000..d4722fdb2a6
--- /dev/null
+++ b/benchmark/scripts/fedavg_dp.dml
@@ -0,0 +1,130 @@
+# ── fedavg_dp.dml ────────────────────────────────────────────────────────────
+# Arguments (passed via -nvargs):
+# data_dir : path to benchmark/data/
+# epsilon : DP privacy budget ε (use 9999 for non-private baseline)
+# delta : DP delta (1e-5 can be used)
+# clip_norm : per-example gradient L2-norm clip bound (default 4.0) —
+# sensitivity = clip_norm / n follows from this, since
+# clipping each record's gradient contribution to clip_norm
+# is what makes that bound actually hold.
+# n_rounds : number of FedAvg rounds (default 50)
+# lr : learning rate (default 0.1)
+# w1_rows : rows in worker 1 shard
+# w2_rows : rows in worker 2 shard
+# w3_rows : rows in worker 3 shard
+# w4_rows : rows in worker 4 shard
+# n_features : number of features
+# out : path to write final weights
+# private : 1 = apply DP noise (default), 0 = non-private baseline
+
+data_dir = $data_dir;
+epsilon = $epsilon;
+delta = $delta;
+clip_norm = ifdef($clip_norm, 4.0);
+n_rounds = ifdef($n_rounds, 50);
+lr = ifdef($lr, 0.1);
+private = ifdef($private, 1);
+out_path = $out;
+
+# Per-round epsilon: dp_gaussian is called once per round, and
+# DPBudgetAccountant composes the cost of every release against the single
+# total budget set by dp_set_budget(). Spending the full epsilon on every
+# round would exhaust the budget almost immediately, so split it evenly
+# across rounds instead.
+round_epsilon = epsilon / n_rounds;
+
+w1r = $w1_rows;
+w2r = $w2_rows;
+w3r = $w3_rows;
+w4r = $w4_rows;
+d = $n_features;
+n = w1r + w2r + w3r + w4r;
+
+# Sensitivity of the released (mean) gradient to a single record changing:
+# with each record's contribution clipped to L2-norm <= clip_norm below,
+# the sum can move by at most clip_norm, so the average moves by clip_norm/n.
+sensitivity = clip_norm / n;
+
+# ── Build federated matrix from 4 local workers ───────────────────────────
+# Row ranges are 0-based [begin, end) for each partition.
+r1s=0; r1e=w1r;
+r2s=w1r; r2e=w1r+w2r;
+r3s=w1r+w2r; r3e=w1r+w2r+w3r;
+r4s=w1r+w2r+w3r; r4e=n;
+
+X = federated(
+ addresses=list(
+ "localhost:8301/" + data_dir + "/worker1/X_train.csv",
+ "localhost:8302/" + data_dir + "/worker2/X_train.csv",
+ "localhost:8303/" + data_dir + "/worker3/X_train.csv",
+ "localhost:8304/" + data_dir + "/worker4/X_train.csv"),
+ ranges=list(
+ list(r1s, 0), list(r1e, d),
+ list(r2s, 0), list(r2e, d),
+ list(r3s, 0), list(r3e, d),
+ list(r4s, 0), list(r4e, d)));
+
+y = federated(
+ addresses=list(
+ "localhost:8301/" + data_dir + "/worker1/y_train.csv",
+ "localhost:8302/" + data_dir + "/worker2/y_train.csv",
+ "localhost:8303/" + data_dir + "/worker3/y_train.csv",
+ "localhost:8304/" + data_dir + "/worker4/y_train.csv"),
+ ranges=list(
+ list(r1s, 0), list(r1e, 1),
+ list(r2s, 0), list(r2e, 1),
+ list(r3s, 0), list(r3e, 1),
+ list(r4s, 0), list(r4e, 1)));
+
+# ── Initialise weights ────────────────────────────────────────────────────
+w = matrix(0, rows=d, cols=1);
+
+if (private == 1) {
+ eps = dp_set_budget($epsilon, $delta);
+}
+
+# ── Training rounds ───────────────────────────────────────────────────────
+for (round in 1:n_rounds) {
+
+ # Forward pass — executes on federated workers.
+ scores = X %*% w; # (n × 1), federated
+ probs = 1.0 / (1.0 + exp(-scores)); # (n × 1), federated
+
+ residuals = probs - y; # (n × 1), federated
+
+ # DP noise injection — only when private=1.
+ if (private == 1) {
+ # Per-example L2-norm clipping: each record's raw gradient
+ # contribution is X[i,:] * residual_i, with norm
+ # ||X[i,:]||_2 * |residual_i|. Scaling it down to clip_norm
+ # whenever it exceeds that bound is what makes
+ # sensitivity = clip_norm / n (set above) an actual, provable
+ # bound instead of an arbitrary constant.
+ row_norms = sqrt(rowSums(X^2)); # (n × 1) ||X[i,:]||_2
+ contrib_norms = abs(residuals) * row_norms; # (n × 1) ||X[i,:]*residual_i||_2
+ clip_scale = clip_norm / pmax(contrib_norms, clip_norm); # (n × 1), in (0,1]
+ clipped_residuals = residuals * clip_scale; # (n × 1)
+
+ # Gradient aggregation — t(X) %*% clipped_residuals is a (d × 1)
+ # local sum that the coordinator collects in one federated
+ # aggregate instruction.
+ grad = t(X) %*% clipped_residuals / n; # (d × 1), LOCAL after agg, clipped
+
+ noisy_grad = dp_gaussian(grad,
+ "identity",
+ sensitivity=sensitivity,
+ epsilon=round_epsilon,
+ delta=delta);
+ w = w - lr * noisy_grad;
+ } else {
+ # Gradient aggregation — t(X) %*% residuals is an (d × 1) local sum
+ # that the coordinator collects in one federated aggregate instruction.
+ grad = t(X) %*% residuals / n; # (d × 1), LOCAL after agg
+ w = w - lr * grad;
+ }
+}
+
+# ── Write model weights ───────────────────────────────────────────────────
+write(w, out_path, format="csv");
+print("FedAvg done. epsilon=" + epsilon + " rounds=" + n_rounds);
+
diff --git a/benchmark/scripts/plot.py b/benchmark/scripts/plot.py
new file mode 100644
index 00000000000..ac0a82850e3
--- /dev/null
+++ b/benchmark/scripts/plot.py
@@ -0,0 +1,114 @@
+"""
+Read results.csv and produce two figures:
+
+1. accuracy_vs_epsilon.png
+ Line plot: x = ε, y = accuracy.
+ Horizontal dashed line = non-private baseline.
+ Points labelled with accuracy values.
+
+2. privacy_cost.png
+ Bar chart showing accuracy loss relative to baseline (utility cost of DP).
+"""
+import pathlib
+import csv
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
+RESULTS = pathlib.Path("benchmark/results")
+
+# ── Load ──────────────────────────────────────────────────────────────────
+rows = []
+with open(RESULTS / "results.csv") as f:
+ for r in csv.DictReader(f):
+ rows.append({
+ "label": r["label"],
+ "epsilon": float(r["epsilon"]) if r["epsilon"] != "inf" else None,
+ "private": int(r["private"]),
+ "accuracy": float(r["accuracy"]),
+ })
+
+baseline = next(r for r in rows if r["private"] == 0)
+dp_rows = sorted([r for r in rows if r["private"] == 1],
+ key=lambda r: r["epsilon"])
+
+eps_vals = [r["epsilon"] for r in dp_rows]
+acc_vals = [r["accuracy"] for r in dp_rows]
+baseline_acc = baseline["accuracy"]
+
+# ── Figure 1: Accuracy vs ε ───────────────────────────────────────────────
+fig, ax = plt.subplots(figsize=(7, 4.5))
+
+ax.plot(eps_vals, acc_vals, marker="o", linewidth=2,
+ color="#028090", label="DP-FedAvg (Gaussian)")
+ax.axhline(baseline_acc, linestyle="--", color="#1C3A5E",
+ linewidth=1.5, label=f"Non-private baseline ({baseline_acc:.3f})")
+
+# Annotate each DP point.
+for eps, acc in zip(eps_vals, acc_vals):
+ ax.annotate(f"{acc:.3f}", xy=(eps, acc),
+ xytext=(0, 8), textcoords="offset points",
+ ha="center", fontsize=9, color="#028090")
+
+ax.set_xscale("log")
+ax.set_xticks(eps_vals)
+ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
+ax.set_xlabel("Privacy budget ε (smaller = stronger privacy)", fontsize=11)
+ax.set_ylabel("Test accuracy", fontsize=11)
+ax.set_title("Accuracy vs. Privacy Budget — DP-FedAvg on Adult (4 workers)",
+ fontsize=12)
+ax.legend(fontsize=9)
+ax.set_ylim(max(0, min(acc_vals) - 0.05), min(1.0, baseline_acc + 0.05))
+ax.grid(True, which="both", linestyle=":", alpha=0.5)
+
+plt.tight_layout()
+out1 = RESULTS / "accuracy_vs_epsilon.png"
+fig.savefig(out1, dpi=150)
+print(f"Saved {out1}")
+plt.close()
+
+# ── Figure 2: Utility cost (accuracy drop) ────────────────────────────────
+fig, ax = plt.subplots(figsize=(6, 4))
+
+drops = [baseline_acc - acc for acc in acc_vals]
+colors = ["#B91C1C" if d > 0.02 else "#028090" for d in drops]
+# Position bars at their true ε value on a log-scaled x-axis (rather than
+# evenly-spaced categorical slots) so the visual spacing between 0.5→1 and
+# 4→8 reflects the same 2x ratio. Bar widths scale with x so they stay a
+# constant fraction of their slot in log space instead of shrinking/growing.
+widths = [e * 0.4 for e in eps_vals]
+bars = ax.bar(eps_vals, drops, color=colors, width=widths,
+ edgecolor="white")
+ax.set_xscale("log")
+ax.set_xticks(eps_vals)
+ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
+
+for bar, drop in zip(bars, drops):
+ ax.text(bar.get_x() + bar.get_width() / 2,
+ bar.get_height() + 0.001,
+ f"{drop:.3f}", ha="center", va="bottom", fontsize=9)
+
+ax.legend(["drop > 0.02", "drop <= 0.02"])
+
+ax.axhline(0, color="black", linewidth=0.8)
+ax.set_xlabel("Privacy budget ε", fontsize=11)
+ax.set_ylabel("Accuracy drop vs. baseline", fontsize=11)
+ax.set_title("Utility Cost of Differential Privacy — DP-FedAvg on Adult",
+ fontsize=11)
+ax.grid(True, axis="y", linestyle=":", alpha=0.5)
+
+plt.tight_layout()
+out2 = RESULTS / "privacy_cost.png"
+fig.savefig(out2, dpi=150)
+print(f"Saved {out2}")
+plt.close()
+
+# ── Console summary table ─────────────────────────────────────────────────
+print()
+print(f"{'ε':>8} {'accuracy':>10} {'drop':>8}")
+print("-" * 34)
+print(f"{'baseline':>8} {baseline_acc:10.4f} {'—':>8}")
+for eps, acc, drop in zip(eps_vals, acc_vals, drops):
+ print(f"{eps:>8.1f} {acc:10.4f} {drop:8.4f}")
+
diff --git a/benchmark/scripts/prepare_data.py b/benchmark/scripts/prepare_data.py
new file mode 100644
index 00000000000..dce8e8e8006
--- /dev/null
+++ b/benchmark/scripts/prepare_data.py
@@ -0,0 +1,123 @@
+"""
+Download the UCI Adult dataset, binarise labels, standardise features,
+split into 4 equal horizontal partitions for federated workers, and write
+SystemDS .mtd metadata files alongside each CSV shard.
+
+Outputs
+-------
+benchmark/data/worker{1..4}/X_train.csv + X_train.csv.mtd
+benchmark/data/worker{1..4}/y_train.csv + y_train.csv.mtd
+benchmark/data/X_test.csv + X_test.csv.mtd
+benchmark/data/y_test.csv + y_test.csv.mtd
+benchmark/data/meta.txt # n_train, n_test, n_features
+"""
+
+import json, os, pathlib
+import numpy as np
+import pandas as pd
+from sklearn.model_selection import train_test_split
+from sklearn.preprocessing import StandardScaler
+
+ADULT_TRAIN_URL = (
+ "https://archive.ics.uci.edu/ml/machine-learning-databases"
+ "/adult/adult.data"
+)
+ADULT_TEST_URL = (
+ "https://archive.ics.uci.edu/ml/machine-learning-databases"
+ "/adult/adult.test"
+)
+
+COLS = [
+ "age","workclass","fnlwgt","education","education_num","marital_status",
+ "occupation","relationship","race","sex","capital_gain","capital_loss",
+ "hours_per_week","native_country","label",
+]
+NUMERIC = ["age","fnlwgt","education_num","capital_gain",
+ "capital_loss","hours_per_week"]
+
+DATA_DIR = pathlib.Path("benchmark/data")
+N_WORKERS = 4
+
+def download(url, dest):
+ import urllib.request
+ if not dest.exists():
+ print(f"Downloading {url}")
+ urllib.request.urlretrieve(url, dest)
+
+def load_adult(path, skip_rows=0):
+ df = pd.read_csv(path, names=COLS, skipinitialspace=True,
+ skiprows=skip_rows, na_values="?").dropna()
+ # binarise label: >50K → 1, else 0
+ df["label"] = (df["label"].str.strip().str.rstrip(".") == ">50K").astype(float)
+ # one-hot encode categoricals
+ cats = [c for c in COLS[:-1] if c not in NUMERIC]
+ df = pd.get_dummies(df, columns=cats, drop_first=True)
+ return df
+
+def write_csv_and_mtd(arr: np.ndarray, path: pathlib.Path, description: str):
+ """Write a CSV and a matching SystemDS .mtd metadata file."""
+ path.parent.mkdir(parents=True, exist_ok=True)
+ np.savetxt(path, arr, delimiter=",", fmt="%.8f")
+ rows, cols = arr.shape
+ mtd = {
+ "data_type": "matrix",
+ "value_type": "double",
+ "rows": rows,
+ "cols": cols,
+ "format": "csv",
+ "header": False,
+ "description": description,
+ }
+ with open(str(path) + ".mtd", "w") as f:
+ json.dump(mtd, f, indent=2)
+ print(f" {path} ({rows} × {cols})")
+
+# ── Download ─────────────────────────────────────────────────────────────────
+download(ADULT_TRAIN_URL, DATA_DIR / "raw" / "adult.data")
+download(ADULT_TEST_URL, DATA_DIR / "raw" / "adult.test")
+
+train_df = load_adult(DATA_DIR / "raw" / "adult.data")
+test_df = load_adult(DATA_DIR / "raw" / "adult.test", skip_rows=1)
+
+# Align columns (test may have different dummies after get_dummies).
+train_df, test_df = train_df.align(test_df, join="left", axis=1, fill_value=0)
+
+# ── Feature / label split ────────────────────────────────────────────────────
+feature_cols = [c for c in train_df.columns if c != "label"]
+X_train = train_df[feature_cols].values.astype(float)
+y_train = train_df["label"].values.reshape(-1, 1).astype(float)
+X_test = test_df[feature_cols].values.astype(float)
+y_test = test_df["label"].values.reshape(-1, 1).astype(float)
+
+# ── Standardise (fit on train only) ─────────────────────────────────────────
+scaler = StandardScaler()
+X_train = scaler.fit_transform(X_train)
+X_test = scaler.transform(X_test)
+
+n_train, n_features = X_train.shape
+n_test = X_test.shape[0]
+print(f"Train: {n_train} × {n_features} | Test: {n_test} × {n_features}")
+
+# ── Partition across workers (equal horizontal splits) ──────────────────────
+indices = np.array_split(np.arange(n_train), N_WORKERS)
+for i, idx in enumerate(indices, start=1):
+ wdir = DATA_DIR / f"worker{i}"
+ write_csv_and_mtd(X_train[idx], wdir / "X_train.csv",
+ f"Adult features, worker {i}")
+ write_csv_and_mtd(y_train[idx], wdir / "y_train.csv",
+ f"Adult labels, worker {i}")
+
+# ── Test set (coordinator-local) ────────────────────────────────────────────
+write_csv_and_mtd(X_test, DATA_DIR / "X_test.csv", "Adult test features")
+write_csv_and_mtd(y_test, DATA_DIR / "y_test.csv", "Adult test labels")
+
+# ── Metadata for DML scripts ─────────────────────────────────────────────────
+worker_rows = [len(idx) for idx in indices]
+with open(DATA_DIR / "meta.txt", "w") as f:
+ f.write(f"n_train={n_train}\n")
+ f.write(f"n_test={n_test}\n")
+ f.write(f"n_features={n_features}\n")
+ for i, r in enumerate(worker_rows, start=1):
+ f.write(f"worker{i}_rows={r}\n")
+print("Wrote benchmark/data/meta.txt")
+
diff --git a/benchmark/scripts/run_benchmark.sh b/benchmark/scripts/run_benchmark.sh
new file mode 100755
index 00000000000..70966fa30ad
--- /dev/null
+++ b/benchmark/scripts/run_benchmark.sh
@@ -0,0 +1,15 @@
+# 1. Prepare data (once).
+python benchmark/scripts/prepare_data.py
+
+# 2. Run the sweep (starts workers, trains, evaluates, stops workers).
+bash benchmark/scripts/run_sweep.sh
+
+# 3. Collect results into CSV.
+python benchmark/scripts/collect_results.py
+
+# 4. Generate plots.
+python benchmark/scripts/plot.py
+
+# 5. Confirm outputs exist.
+ls -lh benchmark/results/accuracy_vs_epsilon.png benchmark/results/privacy_cost.png
+
diff --git a/benchmark/scripts/run_sweep.sh b/benchmark/scripts/run_sweep.sh
new file mode 100755
index 00000000000..4696ea5240a
--- /dev/null
+++ b/benchmark/scripts/run_sweep.sh
@@ -0,0 +1,69 @@
+#!/usr/bin/env bash
+# Runs FedAvg for each epsilon value and the non-private baseline,
+# then evaluates accuracy. Results are appended to results/results.csv.
+set -euo pipefail
+
+REPO_ROOT="$(cd "$(dirname "$0")/../.." && pwd)"
+JAR="$REPO_ROOT/target/SystemDS.jar"
+SCRIPTS="$REPO_ROOT/benchmark/scripts"
+DATA="$REPO_ROOT/benchmark/data"
+RESULTS="$REPO_ROOT/benchmark/results"
+mkdir -p "$RESULTS"
+
+# Read dataset metadata written by prepare_data.py.
+source <(grep -E '^(n_train|n_test|n_features|worker[1-4]_rows)=' \
+ "$DATA/meta.txt" | sed 's/=/="/;s/$/"/')
+
+COMMON_ARGS="\
+ data_dir=$DATA \
+ n_features=$n_features \
+ w1_rows=$worker1_rows \
+ w2_rows=$worker2_rows \
+ w3_rows=$worker3_rows \
+ w4_rows=$worker4_rows \
+ n_rounds=300 \
+ lr=1.0 \
+ clip_norm=4.0 \
+ delta=1e-5"
+
+# ── Start workers ─────────────────────────────────────────────────────────
+echo "=== Starting federated workers ==="
+bash "$SCRIPTS/start_workers.sh"
+
+run_one() {
+ local label="$1" # e.g. "eps_0.5" or "baseline"
+ local extra="$2" # extra -nvargs for this run
+ local model="$RESULTS/model_${label}.csv"
+ local acc_file="$RESULTS/acc_${label}.txt"
+
+ echo ""
+ echo "--- Training: $label ---"
+ java --add-modules=jdk.incubator.vector -jar "$JAR" \
+ -f "$SCRIPTS/fedavg_dp.dml" \
+ -nvargs $COMMON_ARGS $extra out="$model" \
+ 2>&1 | tee "$RESULTS/train_${label}.log" | grep -E "FedAvg|error|Error" || true
+
+ echo " Evaluating …"
+ java --add-modules=jdk.incubator.vector -jar "$JAR" \
+ -f "$SCRIPTS/eval.dml" \
+ -nvargs data_dir="$DATA" model_path="$model" out_acc="$acc_file" \
+ 2>&1 | grep "Accuracy:"
+}
+
+# ── Non-private baseline ──────────────────────────────────────────────────
+run_one "baseline" "private=0 epsilon=9999"
+
+# ── DP runs ───────────────────────────────────────────────────────────────
+for EPS in 0.5 1 4 8; do
+ LABEL="eps_${EPS}"
+ run_one "$LABEL" "private=1 epsilon=${EPS}"
+done
+
+# ── Stop workers ──────────────────────────────────────────────────────────
+echo ""
+echo "=== Stopping federated workers ==="
+bash "$SCRIPTS/stop_workers.sh"
+
+echo ""
+echo "All runs complete. Logs and model files in $RESULTS/"
+
diff --git a/benchmark/scripts/start_workers.sh b/benchmark/scripts/start_workers.sh
new file mode 100755
index 00000000000..23b5634230c
--- /dev/null
+++ b/benchmark/scripts/start_workers.sh
@@ -0,0 +1,24 @@
+#!/usr/bin/env bash
+# Start 4 local SystemDS federated workers on ports 8301-8304.
+# Each worker is given the absolute path to its data shard directory.
+set -e
+REPO_ROOT="$(cd "$(dirname "$0")/../.." && pwd)"
+JAR="$REPO_ROOT/target/SystemDS.jar"
+DATA_DIR="$REPO_ROOT/benchmark/data"
+LOG_DIR="$REPO_ROOT/benchmark/results"
+mkdir -p "$LOG_DIR"
+
+for i in 1 2 3 4; do
+ PORT=$((8300 + i))
+ echo "Starting worker $i on port $PORT …"
+ java --add-modules=jdk.incubator.vector -jar "$JAR" \
+ -w "$PORT" \
+ > "$LOG_DIR/worker${i}.log" 2>&1 &
+ echo $! > "$LOG_DIR/worker${i}.pid"
+done
+
+# Give workers time to bind their ports.
+sleep 3
+echo "Workers running. PIDs:"
+for i in 1 2 3 4; do cat "$LOG_DIR/worker${i}.pid"; done
+
diff --git a/benchmark/scripts/stop_workers.sh b/benchmark/scripts/stop_workers.sh
new file mode 100755
index 00000000000..d1ff3cfde41
--- /dev/null
+++ b/benchmark/scripts/stop_workers.sh
@@ -0,0 +1,11 @@
+#!/usr/bin/env bash
+LOG_DIR="$(cd "$(dirname "$0")/../../benchmark/results" && pwd)"
+for i in 1 2 3 4; do
+ PID_FILE="$LOG_DIR/worker${i}.pid"
+ if [ -f "$PID_FILE" ]; then
+ PID=$(cat "$PID_FILE")
+ kill "$PID" 2>/dev/null && echo "Stopped worker $i (PID $PID)" || true
+ rm -f "$PID_FILE"
+ fi
+done
+
diff --git a/pom.xml b/pom.xml
index 068bed2e8ea..724ee4b1f1d 100644
--- a/pom.xml
+++ b/pom.xml
@@ -408,6 +408,7 @@
maven-surefire-plugin
${maven-surefire-plugin.version}
+ plain
${maven.test.skip}
${test-parallel}
${test-threadCount}
diff --git a/src/main/java/org/apache/sysds/common/Builtins.java b/src/main/java/org/apache/sysds/common/Builtins.java
index f5719641df7..6dce5fb30eb 100644
--- a/src/main/java/org/apache/sysds/common/Builtins.java
+++ b/src/main/java/org/apache/sysds/common/Builtins.java
@@ -116,6 +116,9 @@ public enum Builtins {
DECISIONTREEPREDICT("decisionTreePredict", true),
DECOMPRESS("decompress", false),
DEDUP("dedup", true),
+ DP_LAPLACE("dp_laplace", false),
+ DP_GAUSSIAN("dp_gaussian", false),
+ DP_SET_BUDGET("dp_set_budget", false),
DEEPWALK("deepWalk", true),
DET("det", false),
DETECTSCHEMA("detectSchema", false),
diff --git a/src/main/java/org/apache/sysds/common/InstructionType.java b/src/main/java/org/apache/sysds/common/InstructionType.java
index e0e77c46c59..ed795d038e2 100644
--- a/src/main/java/org/apache/sysds/common/InstructionType.java
+++ b/src/main/java/org/apache/sysds/common/InstructionType.java
@@ -63,6 +63,7 @@ public enum InstructionType {
MMChain,
Union,
EINSUM,
+ DPBuiltin,
//SP Types
MAPMM,
diff --git a/src/main/java/org/apache/sysds/common/Opcodes.java b/src/main/java/org/apache/sysds/common/Opcodes.java
index 9a894dde13b..f9ed57eae06 100644
--- a/src/main/java/org/apache/sysds/common/Opcodes.java
+++ b/src/main/java/org/apache/sysds/common/Opcodes.java
@@ -194,6 +194,10 @@ public enum Opcodes {
LIST("list", InstructionType.BuiltinNary),
EINSUM("einsum", InstructionType.BuiltinNary),
+ //DP built-in functions
+ DP_LAPLACE("dp_laplace", InstructionType.DPBuiltin),
+ DP_GAUSSIAN("dp_gaussian", InstructionType.DPBuiltin),
+
//Parametrized builtin functions
AUTODIFF("autoDiff", InstructionType.ParameterizedBuiltin),
CONTAINS("contains", InstructionType.ParameterizedBuiltin),
diff --git a/src/main/java/org/apache/sysds/common/Types.java b/src/main/java/org/apache/sysds/common/Types.java
index 624c9eed3c6..611d5011fbb 100644
--- a/src/main/java/org/apache/sysds/common/Types.java
+++ b/src/main/java/org/apache/sysds/common/Types.java
@@ -806,7 +806,8 @@ public static ReOrgOp valueOfByOpcode(String opcode) {
/** Parameterized operations that require named variable arguments */
public enum ParamBuiltinOp {
- AUTODIFF, CDF, CONTAINS, INVALID, INVCDF, GROUPEDAGG, RMEMPTY, REPLACE, REXPAND,
+ AUTODIFF, CDF, CONTAINS, DP_LAPLACE, DP_GAUSSIAN, INVALID, INVCDF,
+ GROUPEDAGG, RMEMPTY, REPLACE, REXPAND,
LOWER_TRI, UPPER_TRI,
TRANSFORMAPPLY, TRANSFORMDECODE, TRANSFORMCOLMAP, TRANSFORMMETA,
TOKENIZE, TOSTRING, LIST, PARAMSERV
diff --git a/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java b/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java
index 61a4b8b8f91..3997d3402c8 100644
--- a/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java
+++ b/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java
@@ -182,7 +182,7 @@ public Lop constructLops()
}
case CONTAINS:
case CDF:
- case INVCDF:
+ case INVCDF:
case REPLACE:
case LOWER_TRI:
case UPPER_TRI:
@@ -194,7 +194,9 @@ public Lop constructLops()
case TOSTRING:
case PARAMSERV:
case LIST:
- case AUTODIFF:{
+ case AUTODIFF:
+ case DP_LAPLACE:
+ case DP_GAUSSIAN:{
ParameterizedBuiltin pbilop = new ParameterizedBuiltin(
inputlops, _op, getDataType(), getValueType(), et);
if( isMultiThreadedOpType() )
@@ -688,7 +690,17 @@ else if( _op == ParamBuiltinOp.TRANSFORMAPPLY ) {
return new MatrixCharacteristics(dc.getRows(), dc.getCols(), -1, dc.getLength());
}
}
-
+ else if( _op == ParamBuiltinOp.DP_LAPLACE || _op == ParamBuiltinOp.DP_GAUSSIAN ) {
+ if( dc.dimsKnown() ) {
+ Hop query = getParameterHop("query");
+ String queryVal = (query instanceof LiteralOp) ? ((LiteralOp)query).getStringValue() : null;
+ if( "colMeans".equals(queryVal) || "colSums".equals(queryVal) )
+ ret = new MatrixCharacteristics(1, dc.getCols(), -1, dc.getCols());
+ else if( "identity".equals(queryVal) )
+ ret = new MatrixCharacteristics(dc.getRows(), dc.getCols(), -1, dc.getLength());
+ }
+ }
+
return ret;
}
@Override
@@ -758,7 +770,8 @@ && getTargetHop().areDimsBelowThreshold() ) {
if (_op == ParamBuiltinOp.TRANSFORMCOLMAP || _op == ParamBuiltinOp.TRANSFORMMETA
|| _op == ParamBuiltinOp.TOSTRING || _op == ParamBuiltinOp.LIST
|| _op == ParamBuiltinOp.CDF || _op == ParamBuiltinOp.INVCDF
- || _op == ParamBuiltinOp.PARAMSERV) {
+ || _op == ParamBuiltinOp.PARAMSERV
+ || _op == ParamBuiltinOp.DP_LAPLACE || _op == ParamBuiltinOp.DP_GAUSSIAN) {
_etype = ExecType.CP;
}
diff --git a/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java b/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java
index 3604121aac8..48ecbaa72df 100644
--- a/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java
+++ b/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java
@@ -204,7 +204,19 @@ public String getInstructions(String output)
compileGenericParamMap(sb, _inputParams);
break;
}
-
+ case DP_LAPLACE: {
+ sb.append(Opcodes.DP_LAPLACE);
+ sb.append(OPERAND_DELIMITOR);
+ compileGenericParamMap(sb, _inputParams);
+ break;
+ }
+ case DP_GAUSSIAN: {
+ sb.append(Opcodes.DP_GAUSSIAN);
+ sb.append(OPERAND_DELIMITOR);
+ compileGenericParamMap(sb, _inputParams);
+ break;
+ }
+
default:
throw new LopsException(this.printErrorLocation() + "In ParameterizedBuiltin Lop, Unknown operation: " + _operation);
}
diff --git a/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java b/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java
index ab0c7993b4e..1c2c6ba8110 100644
--- a/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java
+++ b/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java
@@ -2005,6 +2005,59 @@ else if(this.getOpCode() == Builtins.MAX_POOL || this.getOpCode() == Builtins.AV
}
else
raiseValidateError("Local instruction not allowed in dml script");
+ case DP_LAPLACE: {
+ checkNumParameters(4);
+ checkMatrixParam(getFirstExpr());
+ checkScalarParam(getSecondExpr());
+ checkValueTypeParam(getSecondExpr(), ValueType.STRING);
+ checkScalarParam(getThirdExpr());
+ checkScalarParam(getFourthExpr());
+ String dpLaplaceQuery = getDPQueryLiteral(getSecondExpr());
+ long[] dpLaplaceDims = getDPOutputDims(dpLaplaceQuery,
+ getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2());
+ output.setDataType(DataType.MATRIX);
+ output.setValueType(ValueType.FP64);
+ output.setDimensions(dpLaplaceDims[0], dpLaplaceDims[1]);
+ break;
+ }
+ case DP_GAUSSIAN: {
+ checkNumParameters(5);
+ checkMatrixParam(getFirstExpr());
+ checkScalarParam(getSecondExpr());
+ checkValueTypeParam(getSecondExpr(), ValueType.STRING);
+ checkScalarParam(getThirdExpr());
+ checkScalarParam(getFourthExpr());
+ checkScalarParam(getFifthExpr());
+ String dpGaussianQuery = getDPQueryLiteral(getSecondExpr());
+ long[] dpGaussianDims = getDPOutputDims(dpGaussianQuery,
+ getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2());
+ output.setDataType(DataType.MATRIX);
+ output.setValueType(ValueType.FP64);
+ output.setDimensions(dpGaussianDims[0], dpGaussianDims[1]);
+ break;
+ }
+ case DP_SET_BUDGET: {
+ checkNumParameters(2);
+ checkScalarParam(getFirstExpr());
+ checkScalarParam(getSecondExpr());
+ // resolved entirely at compile time (see DMLTranslator), so both
+ // arguments must be known before the HOP DAG is built.
+ if (!isConstant(getFirstExpr()) || !isConstant(getSecondExpr()))
+ raiseValidateError(getOpCode() + ": 'epsilon' and 'delta' must be compile-time numeric literals",
+ false, LanguageErrorCodes.INVALID_PARAMETERS);
+ double dpSetBudgetEpsilon = getDoubleValue(getFirstExpr());
+ double dpSetBudgetDelta = getDoubleValue(getSecondExpr());
+ if (dpSetBudgetEpsilon <= 0)
+ raiseValidateError(getOpCode() + ": epsilon must be > 0, got " + dpSetBudgetEpsilon,
+ false, LanguageErrorCodes.INVALID_PARAMETERS);
+ if ((dpSetBudgetDelta <= 0) || (dpSetBudgetDelta >= 1))
+ raiseValidateError(getOpCode() + ": delta must be in (0,1), got " + dpSetBudgetDelta,
+ false, LanguageErrorCodes.INVALID_PARAMETERS);
+ output.setDataType(DataType.SCALAR);
+ output.setValueType(ValueType.FP64);
+ output.setDimensions(0, 0);
+ break;
+ }
case COMPRESS:
case DECOMPRESS:
if(OptimizerUtils.ALLOW_SCRIPT_LEVEL_COMPRESS_COMMAND){
@@ -2111,6 +2164,34 @@ else if(this.getOpCode() == Builtins.MAX_POOL || this.getOpCode() == Builtins.AV
}
}
+ /**
+ * dp_laplace/dp_gaussian require the "query" parameter to be a compile-time
+ * string literal so that the output shape (and thus the transformation
+ * matrix T built at runtime) is known during validation.
+ */
+ private String getDPQueryLiteral(Expression queryExpr) {
+ if (!(queryExpr instanceof StringIdentifier))
+ raiseValidateError(getOpCode() + ": 'query' must be a string literal", false,
+ LanguageErrorCodes.INVALID_PARAMETERS);
+ return ((StringIdentifier) queryExpr).getValue();
+ }
+
+ /** Output dimensions of T %*% X for the given named query, X being n x d. */
+ private long[] getDPOutputDims(String query, long n, long d) {
+ switch (query) {
+ case "colMeans":
+ case "colSums":
+ return new long[] {1, d};
+ case "identity":
+ return new long[] {n, d};
+ default:
+ raiseValidateError(getOpCode() + ": unknown query type '" + query
+ + "' (expected colMeans, colSums, or identity)", false,
+ LanguageErrorCodes.INVALID_PARAMETERS);
+ return null; // unreachable
+ }
+ }
+
private void validateEinsum(DataIdentifier output){
if(getSecondExpr() == null)
raiseValidateError("Einsum: at least one input matrix required", false,
diff --git a/src/main/java/org/apache/sysds/parser/DMLProgram.java b/src/main/java/org/apache/sysds/parser/DMLProgram.java
index 2f69cb7ea0d..d2c0b9ae7d3 100644
--- a/src/main/java/org/apache/sysds/parser/DMLProgram.java
+++ b/src/main/java/org/apache/sysds/parser/DMLProgram.java
@@ -36,7 +36,21 @@ public class DMLProgram
private ArrayList _blocks;
private Map> _namespaces;
private boolean _containsRemoteParfor;
-
+
+ /**
+ * Session-wide differential privacy budget resolved at compile time from a
+ * {@code dp_set_budget(epsilon, delta)} call (its arguments must be numeric
+ * literals — see {@code BuiltinFunctionExpression}). Null until such a call
+ * is encountered during HOP construction; consulted by
+ * {@code ExecutionContext#getDPBudgetAccountant()} in place of its hardcoded
+ * default. This is deliberately a plain field (not a Hop/Lop/Instruction):
+ * since {@code Program} holds a reference back to this {@code DMLProgram}
+ * (see {@code Program#getDMLProg()}), the value survives from compile time
+ * through to runtime without needing a runtime instruction at all.
+ */
+ private Double _dpBudgetEpsilon;
+ private Double _dpBudgetDelta;
+
public DMLProgram(){
_blocks = new ArrayList<>();
_namespaces = new HashMap<>();
@@ -67,6 +81,23 @@ public void setContainsRemoteParfor(boolean flag) {
public boolean containsRemoteParfor() {
return _containsRemoteParfor;
}
+
+ public void setDPBudget(double epsilon, double delta) {
+ _dpBudgetEpsilon = epsilon;
+ _dpBudgetDelta = delta;
+ }
+
+ public boolean hasDPBudget() {
+ return _dpBudgetEpsilon != null;
+ }
+
+ public double getDPBudgetEpsilon() {
+ return _dpBudgetEpsilon;
+ }
+
+ public double getDPBudgetDelta() {
+ return _dpBudgetDelta;
+ }
public static boolean isInternalNamespace(String namespace) {
return DEFAULT_NAMESPACE.equals(namespace)
diff --git a/src/main/java/org/apache/sysds/parser/DMLTranslator.java b/src/main/java/org/apache/sysds/parser/DMLTranslator.java
index a8e1667d049..a9d5279c26b 100644
--- a/src/main/java/org/apache/sysds/parser/DMLTranslator.java
+++ b/src/main/java/org/apache/sysds/parser/DMLTranslator.java
@@ -2589,6 +2589,53 @@ else if ( sop.equalsIgnoreCase(Opcodes.NOTEQUAL.toString()) )
case DECOMPRESS:
currBuiltinOp = new UnaryOp(target.getName(), target.getDataType(), ValueType.FP64, OpOp1.DECOMPRESS, expr);
break;
+ case DP_LAPLACE: {
+ String[] dpLaplaceParamNames = {"target", "query", "sensitivity", "epsilon"};
+ LinkedHashMap dpLaplaceParams = new LinkedHashMap<>();
+ dpLaplaceParams.put(dpLaplaceParamNames[0], expr);
+ dpLaplaceParams.put(dpLaplaceParamNames[1], expr2);
+ dpLaplaceParams.put(dpLaplaceParamNames[2], expr3);
+ for (int i = 3; i < dpLaplaceParamNames.length; i++) {
+ dpLaplaceParams.put(dpLaplaceParamNames[i],
+ source.getExpr(i) != null ? processExpression(source.getExpr(i), null, hops) : null);
+ }
+ currBuiltinOp = new ParameterizedBuiltinOp(target.getName(), DataType.MATRIX, ValueType.FP64,
+ ParamBuiltinOp.DP_LAPLACE, dpLaplaceParams);
+ break;
+ }
+ case DP_GAUSSIAN: {
+ String[] dpGaussianParamNames = {"target", "query", "sensitivity", "epsilon", "delta"};
+ LinkedHashMap dpGaussianParams = new LinkedHashMap<>();
+ dpGaussianParams.put(dpGaussianParamNames[0], expr);
+ dpGaussianParams.put(dpGaussianParamNames[1], expr2);
+ dpGaussianParams.put(dpGaussianParamNames[2], expr3);
+ for (int i = 3; i < dpGaussianParamNames.length; i++) {
+ dpGaussianParams.put(dpGaussianParamNames[i],
+ source.getExpr(i) != null ? processExpression(source.getExpr(i), null, hops) : null);
+ }
+ currBuiltinOp = new ParameterizedBuiltinOp(target.getName(), DataType.MATRIX, ValueType.FP64,
+ ParamBuiltinOp.DP_GAUSSIAN, dpGaussianParams);
+ break;
+ }
+ case DP_SET_BUDGET: {
+ // Resolved entirely at compile time: BuiltinFunctionExpression.validateExpression
+ // already enforced that both arguments are numeric literals, so 'expr'/'expr2' are
+ // guaranteed LiteralOps here. There is deliberately no runtime Hop/Lop/Instruction for
+ // this call — the budget is applied directly to the DMLProgram (reachable later from
+ // ExecutionContext via Program.getDMLProg(), see ExecutionContext.getDPBudgetAccountant())
+ // before any instruction executes, so there is nothing for the DAG linearizer to reorder
+ // or drop as dead code.
+ if (_dmlProg.hasDPBudget())
+ throw new LanguageException(source.getOpCode() + ": dp_set_budget may only be called once per "
+ + "script (already set to epsilon=" + _dmlProg.getDPBudgetEpsilon()
+ + ", delta=" + _dmlProg.getDPBudgetDelta() + ")");
+ if (!(expr instanceof LiteralOp) || !(expr2 instanceof LiteralOp))
+ throw new LanguageException(source.getOpCode()
+ + ": epsilon and delta must be compile-time numeric literals");
+ _dmlProg.setDPBudget(((LiteralOp) expr).getDoubleValue(), ((LiteralOp) expr2).getDoubleValue());
+ currBuiltinOp = expr; // echo epsilon back as confirmation
+ break;
+ }
case QUANTIZE_COMPRESS:
currBuiltinOp = new BinaryOp(target.getName(), target.getDataType(), target.getValueType(), OpOp2.valueOf(source.getOpCode().name()), expr, expr2);
break;
diff --git a/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java b/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java
index 67cda352a73..7f25b2364fe 100644
--- a/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java
+++ b/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java
@@ -28,6 +28,7 @@
import org.apache.sysds.conf.ConfigurationManager;
import org.apache.sysds.hops.OptimizerUtils;
import org.apache.sysds.hops.fedplanner.FTypes.FType;
+import org.apache.sysds.parser.DMLProgram;
import org.apache.sysds.runtime.DMLRuntimeException;
import org.apache.sysds.runtime.controlprogram.LocalVariableMap;
import org.apache.sysds.runtime.controlprogram.Program;
@@ -62,6 +63,7 @@
import org.apache.sysds.runtime.meta.MetaData;
import org.apache.sysds.runtime.meta.MetaDataFormat;
import org.apache.sysds.runtime.util.HDFSTool;
+import org.apache.sysds.runtime.privacy.dp.DPBudgetAccountant;
import org.apache.sysds.utils.Statistics;
import java.util.ArrayList;
@@ -90,6 +92,8 @@ public class ExecutionContext {
protected SEALClient _seal_client;
+ private DPBudgetAccountant _dpBudgetAccountant = null;
+
//parfor temporary functions (created by eval)
protected Set _fnNames;
@@ -144,6 +148,23 @@ public void setLineage(Lineage lineage) {
_lineage = lineage;
}
+ /**
+ * Returns the session-scoped {@link DPBudgetAccountant}, lazily initialised on
+ * first use. If the DML script called {@code dp_set_budget(epsilon, delta)}
+ * with compile-time literal arguments, that value (resolved onto the
+ * {@code DMLProgram} during HOP construction — see {@code DMLTranslator}'s
+ * {@code DP_SET_BUDGET} case) is used instead of the hardcoded defaults.
+ */
+ public DPBudgetAccountant getDPBudgetAccountant() {
+ if (_dpBudgetAccountant == null) {
+ DMLProgram dmlProg = (_prog != null) ? _prog.getDMLProg() : null;
+ _dpBudgetAccountant = (dmlProg != null && dmlProg.hasDPBudget())
+ ? new DPBudgetAccountant(dmlProg.getDPBudgetEpsilon(), dmlProg.getDPBudgetDelta())
+ : new DPBudgetAccountant();
+ }
+ return _dpBudgetAccountant;
+ }
+
public boolean isAutoCreateVars() {
return _autoCreateVars;
}
diff --git a/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java b/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java
index 92e11b425dd..ca6baf058b4 100644
--- a/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java
+++ b/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java
@@ -65,6 +65,7 @@
import org.apache.sysds.runtime.instructions.cp.UnaryCPInstruction;
import org.apache.sysds.runtime.instructions.cp.VariableCPInstruction;
import org.apache.sysds.runtime.instructions.cp.UnionCPInstruction;
+import org.apache.sysds.runtime.instructions.cp.DPBuiltinCPInstruction;
import org.apache.sysds.runtime.instructions.cp.EinsumCPInstruction;
import org.apache.sysds.runtime.instructions.cpfile.MatrixIndexingCPFileInstruction;
@@ -226,7 +227,10 @@ public static CPInstruction parseSingleInstruction ( InstructionType cptype, Str
case EINSUM:
return EinsumCPInstruction.parseInstruction(str);
-
+
+ case DPBuiltin:
+ return DPBuiltinCPInstruction.parseInstruction(str);
+
default:
throw new DMLRuntimeException("Invalid CP Instruction Type: " + cptype );
}
diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java
index b8d84ca3898..668d2f36978 100644
--- a/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java
+++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java
@@ -48,7 +48,8 @@ public enum CPType {
EvictLineageCache, EINSUM,
NoOp,
Union,
- QuantizeCompression
+ QuantizeCompression,
+ DPBuiltin
}
protected final CPType _cptype;
diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/DPBuiltinCPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/DPBuiltinCPInstruction.java
new file mode 100755
index 00000000000..1ca7d75bf86
--- /dev/null
+++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/DPBuiltinCPInstruction.java
@@ -0,0 +1,428 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied. See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.runtime.instructions.cp;
+
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.controlprogram.context.ExecutionContext;
+import org.apache.sysds.runtime.instructions.InstructionUtils;
+import org.apache.sysds.runtime.matrix.data.LibMatrixMult;
+import org.apache.sysds.runtime.matrix.data.LibMatrixReorg;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.runtime.privacy.dp.DPBudgetAccountant;
+
+import java.util.LinkedHashMap;
+import java.util.concurrent.ThreadLocalRandom;
+
+/**
+ * CP instruction for differential-privacy release of a linear query over the
+ * original matrix.
+ *
+ * DML syntax (raw-matrix form):
+ * result = dp_laplace(X, query="colMeans", sensitivity=1.0, epsilon=0.5)
+ * result = dp_gaussian(X, query="colMeans", sensitivity=1.0, epsilon=0.5, delta=1e-5)
+ *
+ * The instruction receives the original {@code n x d} matrix {@code X},
+ * builds a transformation matrix {@code T} ({@code k x n}) from the named
+ * {@code query} (see {@link #buildTransform}), and returns a noisy release of
+ * {@code T %*% X}. The noise is not added as a separate elementwise
+ * pass over a materialised aggregate: it is injected by augmenting {@code T}
+ * with an identity block and {@code X} with the noise matrix, so that the
+ * noisy release is the result of a single {@link LibMatrixMult#matrixMult}
+ * call (see {@link #processInstruction} for the derivation).
+ *
+ * Sensitivity norm: {@code sensitivity} is not interchangeable
+ * between the two builtins. {@code dp_laplace} calibrates its noise scale
+ * to the L1 sensitivity of {@code T %*% X} to a single-record
+ * change; {@code dp_gaussian} calibrates its σ to the L2 sensitivity.
+ * For a scalar release (e.g. {@code query="colMeans"} on single-column
+ * {@code X}) the two norms coincide, but for a vector- or matrix-valued
+ * release they generally differ — the caller is responsible for supplying
+ * the norm matching the builtin invoked (see {@link #sensitivityOf}).
+ *
+ * The {@link #sensitivityOf} method is deliberately separated from the
+ * noise-scale computation. It currently returns the caller-supplied
+ * constant. A future rewrite pass could replace the body of this single
+ * method with a static analysis that derives sensitivity from {@code T}'s
+ * column norms and a declared per-record bound on {@code X}; every other
+ * line in this class would stay unchanged.
+ */
+public class DPBuiltinCPInstruction extends ComputationCPInstruction {
+
+ // -----------------------------------------------------------------------
+ // Constants
+ // -----------------------------------------------------------------------
+
+ /** Opcode registered in Builtins and CPInstructionParser. */
+ public static final String OPCODE_LAPLACE = "dp_laplace";
+ public static final String OPCODE_GAUSSIAN = "dp_gaussian";
+
+ // -----------------------------------------------------------------------
+ // Fields
+ // -----------------------------------------------------------------------
+
+ /**
+ * Named parameters extracted from the serialised instruction string.
+ * Keys: "target", "query", "sensitivity", "epsilon", "delta" (Gaussian only).
+ *
+ * Using the same LinkedHashMap convention as
+ * ParameterizedBuiltinCPInstruction so that CPInstructionParser can
+ * call the shared constructParameterMap() helper unchanged.
+ */
+ private final LinkedHashMap _params;
+
+ // -----------------------------------------------------------------------
+ // Constructor (private – use parseInstruction)
+ // -----------------------------------------------------------------------
+
+ private DPBuiltinCPInstruction(
+ CPOperand input,
+ CPOperand output,
+ String opcode,
+ String istr,
+ LinkedHashMap params) {
+ super(CPType.DPBuiltin, null, input, null, output, opcode, istr);
+ _params = params;
+ }
+
+ // -----------------------------------------------------------------------
+ // Static factory / parser
+ // -----------------------------------------------------------------------
+
+ /**
+ * Reconstructs a {@code DPBuiltinCPInstruction} from its serialised
+ * instruction string produced by the LOP layer.
+ *
+ * Expected format (OPERAND_DELIM = '\u00b0'):
+ * dp_gaussian°target=mVar1·MATRIX·FP64°query=colMeans·SCALAR·STRING·true
+ * °sensitivity=1.0·SCALAR·FP64·true°epsilon=0.5·SCALAR·FP64·true
+ * °delta=1e-5·SCALAR·FP64·true°_mVar2·MATRIX·FP64
+ *
+ * The first token is always the opcode; the last token is always the
+ * output operand; the tokens in between are key=value pairs. This matches
+ * the convention used by ParameterizedBuiltinCPInstruction exactly.
+ */
+ public static DPBuiltinCPInstruction parseInstruction(String str) {
+ String[] parts = InstructionUtils.getInstructionPartsWithValueType(str);
+ InstructionUtils.checkNumFields(parts, 5, 6); // laplace=5, gaussian=6
+ String opcode = parts[0];
+
+ // Output operand is always the last token.
+ CPOperand output = new CPOperand(parts[parts.length - 1]);
+
+ // The "target" parameter holds the variable name of the input matrix.
+ // ParameterizedBuiltinCPInstruction.constructParameterMap strips the
+ // type suffixes and returns bare key=value pairs.
+ LinkedHashMap params =
+ ParameterizedBuiltinCPInstruction.constructParameterMap(parts);
+
+ // The target CPOperand is needed by ComputationCPInstruction's
+ // getInputs() / getLineageItem() machinery.
+ CPOperand input = new CPOperand(params.get("target"),
+ org.apache.sysds.common.Types.ValueType.FP64,
+ org.apache.sysds.common.Types.DataType.MATRIX);
+
+ // Validate required keys.
+ if (!params.containsKey("query"))
+ throw new DMLRuntimeException(opcode + ": missing 'query'");
+ if (!params.containsKey("sensitivity"))
+ throw new DMLRuntimeException(opcode + ": missing 'sensitivity'");
+ if (!params.containsKey("epsilon"))
+ throw new DMLRuntimeException(opcode + ": missing 'epsilon'");
+ if (opcode.equals(OPCODE_GAUSSIAN) && !params.containsKey("delta"))
+ throw new DMLRuntimeException(opcode + ": missing 'delta'");
+
+ return new DPBuiltinCPInstruction(input, output, opcode, str, params);
+ }
+
+ // -----------------------------------------------------------------------
+ // Core execution
+ // -----------------------------------------------------------------------
+
+ /**
+ * Executes the DP release.
+ *
+ * - Read the original {@link MatrixBlock} {@code X} from the variable
+ * table.
+ * - Build the transformation matrix {@code T} ({@code k x n}) from
+ * {@code query} (see {@link #buildTransform}).
+ * - Determine sensitivity via {@link #sensitivityOf}.
+ * - Generate a noise {@link MatrixBlock} shaped {@code k x d}.
+ * - Fuse {@code T %*% X + noise} into a single
+ * {@link LibMatrixMult#matrixMult} call (see below).
+ * - Record the release with the session-scoped
+ * {@link DPBudgetAccountant}; throw if budget is exhausted.
+ * - Write the noisy block back to the variable table and release
+ * the input pin.
+ *
+ * Fusion derivation: for {@code T} ({@code k x n}), {@code X}
+ * ({@code n x d}) and noise {@code N} ({@code k x d}), let
+ * {@code T' = [T | I_k]} ({@code k x (n+k)}) and
+ * {@code X' = [X ; N]} ({@code (n+k) x d}). Then
+ * {@code T' %*% X' = T %*% X + I_k %*% N = T %*% X + N}, computed as one
+ * matrix multiply instead of a multiply followed by a separate
+ * elementwise add.
+ */
+ @Override
+ public void processInstruction(ExecutionContext ec) {
+
+ // ── 1. Read original input matrix X ─────────────────────────────────
+ // getMatrixInput pins the block in memory and increments the
+ // reference count; we must call releaseMatrixInput afterwards.
+ MatrixBlock X = ec.getMatrixInput(_params.get("target"));
+
+ // ── 2. Parse DP parameters ──────────────────────────────────────────
+ double epsilon = parsePositiveDouble("epsilon");
+ double delta = instOpcode.equals(OPCODE_GAUSSIAN)
+ ? parsePositiveDouble("delta") : 0.0;
+ String query = _params.get("query");
+
+ // ── 3. Build the transformation matrix T (k x n) ────────────────────
+ MatrixBlock T = buildTransform(query, X.getNumRows());
+
+ // ── 4. Determine sensitivity (caller-supplied constant) ─────────────
+ double sensitivity = sensitivityOf(T);
+
+ // ── 5. Generate noise shaped like the release T %*% X (k x d) ───────
+ MatrixBlock noiseBlock = generateNoise(T.getNumRows(), X.getNumColumns(),
+ sensitivity, epsilon, delta);
+
+ // ── 6. Fuse T %*% X + noise into a single matrix multiply ───────────
+ MatrixBlock Ik = identity(T.getNumRows());
+ MatrixBlock Tp = T.append(Ik, null, true); // [T | I_k]
+ MatrixBlock Xp = X.append(noiseBlock, null, false); // [X ; noise]
+ MatrixBlock outBlock = LibMatrixMult.matrixMult(Tp, Xp);
+
+ // ── 7. Record release and enforce budget ────────────────────────────
+ // getDPBudgetAccountant() returns a lazy-initialised DPBudgetAccountant that is
+ // owned by this ExecutionContext (added in a companion EC patch).
+ DPBudgetAccountant accountant = ec.getDPBudgetAccountant();
+ accountant.compose(epsilon, delta, sensitivity); // throws on exhaustion
+
+ // ── 8. Write output and release input pin ───────────────────────────
+ ec.releaseMatrixInput(_params.get("target"));
+ ec.setMatrixOutput(output.getName(), outBlock);
+ }
+
+ // -----------------------------------------------------------------------
+ // Transformation matrix construction
+ // -----------------------------------------------------------------------
+
+ /**
+ * Builds the {@code k x n} transformation matrix {@code T} for the given
+ * named query, to be left-multiplied against the {@code n x d} input
+ * {@code X} as {@code T %*% X}.
+ *
+ * - {@code "colMeans"}: {@code T} is {@code 1 x n}, filled with
+ * {@code 1/n} — {@code T %*% X} is the column-mean row vector.
+ * - {@code "colSums"}: {@code T} is {@code 1 x n}, filled with
+ * {@code 1.0} — {@code T %*% X} is the column-sum row vector.
+ * - {@code "identity"}: {@code T} is the {@code n x n} identity
+ * (built sparsely via {@link #identity}) — {@code T %*% X} is
+ * {@code X} itself, i.e. a noisy release of the raw matrix.
+ *
+ * Row-wise aggregates ({@code rowMeans}/{@code rowSums}) reduce across
+ * the feature axis of {@code X}, i.e. they are naturally
+ * {@code X %*% T'} (right-multiply), not {@code T %*% X}, so they are
+ * intentionally not supported here.
+ */
+ private static MatrixBlock buildTransform(String query, int n) {
+ switch (query) {
+ case "colMeans": {
+ MatrixBlock T = new MatrixBlock(1, n, false);
+ T.allocateDenseBlock();
+ double v = 1.0 / n;
+ for (int c = 0; c < n; c++)
+ T.set(0, c, v);
+ T.recomputeNonZeros();
+ return T;
+ }
+ case "colSums": {
+ MatrixBlock T = new MatrixBlock(1, n, false);
+ T.allocateDenseBlock();
+ for (int c = 0; c < n; c++)
+ T.set(0, c, 1.0);
+ T.recomputeNonZeros();
+ return T;
+ }
+ case "identity":
+ return identity(n);
+ default:
+ throw new DMLRuntimeException(
+ "dp_laplace/dp_gaussian: unknown query type '" + query
+ + "' (expected colMeans, colSums, or identity)");
+ }
+ }
+
+ /**
+ * Builds a {@code k x k} identity matrix, sparsely, by reusing the
+ * existing {@link LibMatrixReorg#diag} reorg operator (the same runtime
+ * path DML's {@code diag()} builtin uses to expand a vector into a
+ * diagonal matrix). Keeps memory {@code O(k)} rather than {@code O(k^2)},
+ * which matters for the {@code query="identity"} case where {@code k}
+ * equals the number of rows of {@code X}.
+ */
+ private static MatrixBlock identity(int k) {
+ MatrixBlock ones = new MatrixBlock(k, 1, false);
+ ones.allocateDenseBlock();
+ for (int i = 0; i < k; i++)
+ ones.set(i, 0, 1.0);
+ ones.recomputeNonZeros();
+ return LibMatrixReorg.diag(ones, new MatrixBlock(k, k, true));
+ }
+
+ // -----------------------------------------------------------------------
+ // Sensitivity seam
+ // -----------------------------------------------------------------------
+
+ /**
+ * Returns the sensitivity of the release {@code T %*% X} to a
+ * single-record change, in the norm required by the mechanism actually
+ * invoked: L1 for {@code dp_laplace}, L2 for
+ * {@code dp_gaussian} (see the class Javadoc). The two only coincide
+ * when the release is scalar.
+ *
+ * Returns the caller-supplied literal from the DML script as-is, with
+ * no norm conversion or validation — the DML author must compute the
+ * sensitivity in the correct norm for the builtin they call. A future
+ * rewrite pass could replace this body with an analysis that derives
+ * sensitivity from {@code T}'s column norms and a declared per-record
+ * bound on {@code X}; no other line in this class would need to change.
+ *
+ * @param T the transformation matrix (unused for now; kept as the seam
+ * for a future sensitivity-derivation pass)
+ * @return caller-supplied sensitivity constant, expected to already be
+ * in the L1 norm (Laplace) or L2 norm (Gaussian)
+ */
+ private double sensitivityOf(MatrixBlock T) {
+ return parsePositiveDouble("sensitivity");
+ }
+
+ // -----------------------------------------------------------------------
+ // Noise generation
+ // -----------------------------------------------------------------------
+
+ /**
+ * Generates a {@code rows x cols} noise {@link MatrixBlock} — matching
+ * the shape of the release {@code T %*% X} — filled with samples from the
+ * mechanism-appropriate distribution calibrated to ({@code sensitivity},
+ * {@code epsilon}, {@code delta}).
+ *
+ * Both mechanisms produce a dense block. Sparsity exploitation is
+ * left for future work; for the releases targeted here (e.g. column
+ * means, column sums) the noise is dense regardless.
+ */
+ private MatrixBlock generateNoise(
+ int rows,
+ int cols,
+ double sensitivity,
+ double epsilon,
+ double delta) {
+
+ MatrixBlock noise = new MatrixBlock(rows, cols, false); // dense
+ noise.allocateDenseBlock();
+
+ if (instOpcode.equals(OPCODE_LAPLACE)) {
+ // Laplace mechanism
+ // For a given epsilon, noise is drawn from the Laplace distribution at
+ // scale b = sensitivity / epsilon
+ fillLaplaceNoise(noise, sensitivity / epsilon);
+ } else {
+ // Gaussian mechanism: calibrate sigma for (epsilon, delta)-DP.
+ // For a given epsilon, noise is drawn from the normal distribution at
+ // sigma^2 = 2 * sensitivity^2 * log(1.25/delta) / epsilon^2
+ double sigma = sensitivity
+ * Math.sqrt(2.0 * Math.log(1.25 / delta))
+ / epsilon;
+ fillGaussianNoise(noise, sigma);
+ }
+
+ noise.recomputeNonZeros();
+ return noise;
+ }
+
+ /**
+ * Fills {@code block} with i.i.d. Laplace(0, scale) samples using the
+ * inverse-CDF method.
+ *
+ * For u ~ Uniform(0, 1): X = -scale * sign(u - 0.5) * ln(1 - 2|u - 0.5|)
+ */
+ private static void fillLaplaceNoise(MatrixBlock block, double scale) {
+ ThreadLocalRandom rng = ThreadLocalRandom.current();
+ int rows = block.getNumRows();
+ int cols = block.getNumColumns();
+ for (int r = 0; r < rows; r++) {
+ for (int c = 0; c < cols; c++) {
+ double u = rng.nextDouble(); // u in (0, 1)
+ double v = u - 0.5;
+ // Guard against the degenerate u == 0.5 case (ln(0) = -inf).
+ if (v == 0.0) v = 1e-15;
+ double sample = -scale * Math.signum(v) * Math.log(1.0 - 2.0 * Math.abs(v));
+ block.set(r, c, sample);
+ }
+ }
+ }
+
+ /**
+ * Fills {@code block} with i.i.d. N(0, sigma²) samples.
+ *
+ * Uses {@link ThreadLocalRandom#nextGaussian()} which is thread-safe
+ * and does not require external libraries.
+ */
+ private static void fillGaussianNoise(MatrixBlock block, double sigma) {
+ ThreadLocalRandom rng = ThreadLocalRandom.current();
+ int rows = block.getNumRows();
+ int cols = block.getNumColumns();
+ for (int r = 0; r < rows; r++) {
+ for (int c = 0; c < cols; c++) {
+ block.set(r, c, sigma * rng.nextGaussian());
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Helpers
+ // -----------------------------------------------------------------------
+
+ /**
+ * Parses a parameter value as a positive {@code double}.
+ *
+ * @throws DMLRuntimeException if the key is absent, unparseable, or
+ * non-positive
+ */
+ private double parsePositiveDouble(String key) {
+ String raw = _params.get(key);
+ if (raw == null)
+ throw new DMLRuntimeException(
+ instOpcode + ": parameter '" + key + "' is missing");
+ double v;
+ try {
+ v = Double.parseDouble(raw);
+ } catch (NumberFormatException e) {
+ throw new DMLRuntimeException(
+ instOpcode + ": parameter '" + key
+ + "' is not a valid number: " + raw);
+ }
+ if (!(v > 0.0))
+ throw new DMLRuntimeException(
+ instOpcode + ": parameter '" + key
+ + "' must be strictly positive, got " + v);
+ return v;
+ }
+}
diff --git a/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java b/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java
new file mode 100644
index 00000000000..11ddc3ac79f
--- /dev/null
+++ b/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java
@@ -0,0 +1,277 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied. See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.runtime.privacy.dp;
+
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.instructions.cp.DPBuiltinCPInstruction;
+
+/**
+ * Session-scoped differential privacy budget accountant.
+ *
+ * Tracks composition of DP releases across the lifetime of a DML script
+ * execution. Each call to {@link #compose} records one release and checks
+ * whether the cumulative privacy cost has exceeded the user-specified budget.
+ *
+ * The mechanism type (Laplace vs Gaussian) is inferred from the {@code delta}
+ * argument passed to {@link #compose}:
+ *
+ * - Laplace (delta == 0): pure ε-DP. The budget cost is tracked via
+ * basic composition — each release contributes exactly its ε to a running
+ * sum. This is the tightest possible bound for pure DP and avoids the
+ * looser estimate that results from routing Laplace through the RDP
+ * conversion path (which would introduce an unnecessary δ). Noise scale
+ * is calibrated to L1 sensitivity (see {@link #compose}).
+ * - Gaussian (delta > 0): (ε, δ)-DP via Rényi DP composition.
+ * Rényi divergences at a discrete set of orders α compose additively;
+ * the accumulated sum is converted to (ε, δ) at query time using the
+ * formula from Mironov 2017. This is substantially tighter than basic
+ * composition for repeated Gaussian releases, which is the common case
+ * in federated learning.
+ *
+ * When both mechanisms are used in the same script the total cost is:
+ * ε_total = ε_Laplace_sum + ε_Gaussian_RDP
+ * This follows from basic composition of a pure-DP mechanism with an
+ * approximate-DP mechanism, which is additive in ε.
+ *
+ * Rényi orders tracked (Gaussian path)
+ * α ∈ {2, 4, 8, 16, 32, 64, 128, 256, 512, 1024}. At query time the minimum
+ * converted ε across all orders is taken as the tightest available bound.
+ *
+ * Gaussian RDP divergence
+ * For the Gaussian mechanism with noise scale σ and L2 sensitivity Δf:
+ * D_α = α · Δf² / (2σ²)
+ * σ is back-derived from the caller's (ε, δ) via the standard calibration
+ * formula (see {@link #gaussianSigma}). Note that sensitivity cancels in the
+ * final expression, so the RDP cost depends only on the (ε, δ) parameters.
+ *
+ * RDP → (ε, δ) conversion (Mironov 2017, Proposition 3)
+ * ε(α) = R[α] + log(1 − 1/α) − log(δ·(α−1)) / α
+ *
+ * One instance is created per {@code ExecutionContext} (lazy init). It is
+ * garbage-collected with the context when the script finishes; no state
+ * leaks between script executions or between concurrent scripts.
+ *
+ * Not thread-safe. A single DML script executes instructions sequentially
+ * on one thread, so no synchronisation is needed.
+ *
+ * @see DPBuiltinCPInstruction
+ */
+public class DPBudgetAccountant {
+
+ // -----------------------------------------------------------------------
+ // Rényi orders used for Gaussian composition
+ // -----------------------------------------------------------------------
+
+ private static final double DEFAULT_EPSILON_BUDGET = 1.0;
+
+ private static final double DEFAULT_DELTA = 1e-5;
+
+ /**
+ * Discrete set of Rényi orders α. All must be > 1.
+ * Finer grids give tighter bounds; this set covers the range relevant
+ * for typical ML workloads.
+ */
+ private static final double[] ORDERS = {
+ 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024
+ };
+
+ // -----------------------------------------------------------------------
+ // State
+ // -----------------------------------------------------------------------
+
+ /** Accumulated Rényi divergence at each order (Gaussian releases only). */
+ private final double[] _rdpSum = new double[ORDERS.length];
+
+ /**
+ * Running sum of pure ε from Laplace releases.
+ *
+ * Laplace gives pure ε-DP (no δ). Basic composition is exact and
+ * tighter than the RDP conversion path for Laplace (which would introduce
+ * an unnecessary δ and produce a looser bound). Each Laplace release adds
+ * its ε here; the total is added directly in {@link #totalEpsilonSpent()}.
+ */
+ private double _pureEpsilonSum = 0.0;
+
+ /** Total privacy budget (ε) for the script execution. */
+ private final double _epsilonBudget;
+
+ /** δ used for the Gaussian RDP-to-(ε,δ) conversion. */
+ private final double _delta;
+
+ /** Number of releases recorded so far (for error messages). */
+ private int _releaseCount = 0;
+
+ // -----------------------------------------------------------------------
+ // Constructors
+ // -----------------------------------------------------------------------
+
+ /**
+ * Creates an accountant with the given global budget.
+ *
+ * Typical usage: the DML script sets the budget once at the top
+ * (future work: a {@code dp_set_budget(epsilon, delta)} built-in),
+ * or the accountant is created with defaults and the budget is checked
+ * after each release.
+ *
+ * @param epsilonBudget total ε budget for the script execution (must be > 0)
+ * @param delta δ used for the Gaussian RDP-to-(ε,δ) conversion (must be in (0,1))
+ */
+ public DPBudgetAccountant(double epsilonBudget, double delta) {
+ if (!(epsilonBudget > 0))
+ throw new DMLRuntimeException(
+ "DPBudgetAccountant: epsilonBudget must be > 0, got " + epsilonBudget);
+ if (!(delta > 0 && delta < 1))
+ throw new DMLRuntimeException(
+ "DPBudgetAccountant: delta must be in (0,1), got " + delta);
+ _epsilonBudget = epsilonBudget;
+ _delta = delta;
+ }
+
+ /**
+ * Convenience constructor using a liberal default δ = 1e-5.
+ * Suitable when the calling script does not specify δ explicitly.
+ */
+ public DPBudgetAccountant(double epsilonBudget) {
+ this(epsilonBudget, 1e-5);
+ }
+
+ /**
+ * Default constructor using defaults.
+ * Suitable when the calling script does not specify ε, δ explicitly.
+ */
+ public DPBudgetAccountant() {
+ this(DEFAULT_EPSILON_BUDGET, DEFAULT_DELTA);
+ }
+
+ // -----------------------------------------------------------------------
+ // Core API
+ // -----------------------------------------------------------------------
+
+ /**
+ * Records one DP release and checks the budget.
+ *
+ * This method must be called before the result is written to
+ * the variable table. If the budget is exhausted it throws and the
+ * caller's result is discarded, preventing an unaccounted release.
+ *
+ * Mechanism selection (see class-level Javadoc for details):
+ * - {@code delta == 0} → Laplace, pure ε-DP basic composition
+ * - {@code delta > 0} → Gaussian, Rényi DP composition
+ *
+ * @param epsilon per-release ε parameter (must be > 0)
+ * @param delta per-release δ parameter (0 for Laplace, >0 for Gaussian)
+ * @param sensitivity sensitivity Δf of the released quantity (must be > 0).
+ * The norm depends on the mechanism selected by
+ * {@code delta}: callers must supply the
+ * L1 sensitivity ‖f(D) − f(D′)‖₁ when
+ * {@code delta == 0} (Laplace), and the L2
+ * sensitivity ‖f(D) − f(D′)‖₂ when {@code delta > 0}
+ * (Gaussian). The two coincide for scalar-valued
+ * releases but diverge for vector-valued ones, so
+ * passing the wrong norm silently under- or
+ * over-calibrates the noise.
+ * @throws DMLRuntimeException if the cumulative ε after this release
+ * would exceed the budget
+ */
+ public void compose(double epsilon, double delta, double sensitivity) {
+ _releaseCount++;
+
+ if (delta == 0.0) {
+ // Laplace: pure ε-DP, basic composition — cost is exactly epsilon.
+ _pureEpsilonSum += epsilon;
+ } else {
+ // Gaussian: accumulate Rényi divergence at each order, then convert.
+ for (int i = 0; i < ORDERS.length; i++) {
+ double sigma = gaussianSigma(sensitivity, epsilon, delta);
+ _rdpSum[i] += rdpGaussian(ORDERS[i], sensitivity, sigma);
+ }
+ }
+
+ double spentEpsilon = totalEpsilonSpent();
+ if (spentEpsilon > _epsilonBudget) {
+ throw new DMLRuntimeException(String.format(
+ "Privacy budget exhausted after %d release(s): "
+ + "spent ε ≈ %.6f exceeds budget ε = %.6f (δ = %.2e). "
+ + "Reduce the number of releases or widen the budget.",
+ _releaseCount, spentEpsilon, _epsilonBudget, _delta));
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Inspection
+ // -----------------------------------------------------------------------
+
+ /**
+ * Returns the current total privacy cost as an ε value.
+ *
+ * Total = Laplace pure-ε sum + Gaussian RDP-converted ε (clamped to
+ * zero when no Gaussian releases have been recorded).
+ */
+ public double totalEpsilonSpent() {
+ // Take min_α(ε_α) as the current total privacy cost
+ double gaussianEps = Double.MAX_VALUE;
+ for (int i = 0; i < ORDERS.length; i++) {
+ double alpha = ORDERS[i];
+ double eps = _rdpSum[i]
+ + Math.log(1.0 - 1.0 / alpha)
+ - Math.log(_delta * (alpha - 1.0)) / alpha;
+ if (eps < gaussianEps)
+ gaussianEps = eps;
+ }
+ // Clamp: with no Gaussian releases the RDP sum is 0 and the log-delta
+ // term alone drives gaussianEps to a small positive value; clamp to 0
+ // so Laplace-only scripts are not penalised by δ they never requested.
+ if (gaussianEps < 0) gaussianEps = 0.0;
+ return _pureEpsilonSum + gaussianEps;
+ }
+
+ /** Returns the remaining ε budget (negative if the budget is exceeded). */
+ public double remainingBudget() {
+ return _epsilonBudget - totalEpsilonSpent();
+ }
+
+ /** Returns the number of DP releases recorded so far. */
+ public int releaseCount() {
+ return _releaseCount;
+ }
+
+ // -----------------------------------------------------------------------
+ // Private helpers
+ // -----------------------------------------------------------------------
+
+ /**
+ * Rényi divergence of order α for the Gaussian mechanism (Mironov 2017,
+ * Proposition 3, example 2):
+ * D_α = α · Δf² / (2σ²)
+ */
+ private static double rdpGaussian(double alpha, double sensitivity, double sigma) {
+ return alpha * (sensitivity * sensitivity) / (2.0 * sigma * sigma);
+ }
+
+ /**
+ * Gaussian noise scale σ calibrated to (ε, δ)-DP:
+ * σ = Δf · sqrt(2 · log(1.25 / δ)) / ε
+ * Must match the formula used in {@link DPBuiltinCPInstruction} so that
+ * the RDP cost recorded here is consistent with the noise actually injected.
+ */
+ private static double gaussianSigma(double sensitivity, double epsilon, double delta) {
+ return sensitivity * Math.sqrt(2.0 * Math.log(1.25 / delta)) / epsilon;
+ }
+}
diff --git a/src/main/python/requirements.txt b/src/main/python/requirements.txt
new file mode 100644
index 00000000000..5d17a102c29
--- /dev/null
+++ b/src/main/python/requirements.txt
@@ -0,0 +1,31 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements. See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership. The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License. You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied. See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+numpy
+pandas
+scipy
+py4j
+wheel
+requests
+setuptools
+
+scikit-learn
+matplotlib
diff --git a/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java b/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java
new file mode 100755
index 00000000000..c10dd065fd8
--- /dev/null
+++ b/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java
@@ -0,0 +1,418 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied. See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.test.component.cp;
+
+import org.apache.sysds.parser.DMLProgram;
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.controlprogram.Program;
+import org.apache.sysds.runtime.controlprogram.context.ExecutionContext;
+import org.apache.sysds.runtime.controlprogram.context.ExecutionContextFactory;
+import org.apache.sysds.runtime.privacy.dp.DPBudgetAccountant;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+
+/**
+ * Tests for {@code DPBuiltinCPInstruction} and {@code DPBudgetAccountant}.
+ *
+ * The tests are grouped into three levels:
+ * - Unit tests on DPBudgetAccountant — verify composition, conversion,
+ * and budget enforcement in isolation, with no dependency on the full
+ * SystemDS runtime.
+ * - Noise distribution tests — verify that the noise blocks
+ * generated by the Laplace and Gaussian mechanisms have statistically
+ * correct means and variances (Kolmogorov-Smirnov style sanity checks).
+ * - DML integration tests — run complete DML scripts and verify
+ * end-to-end correctness via the existing AutomatedTestBase machinery.
+ *
+ * The DML integration tests require a built SystemDS jar and are separated
+ * into a companion class {@link org.apache.sysds.test.functions.privacy.dp.DPBuiltinDMLTest}.
+ */
+public class DPBuiltinCPInstructionTest {
+
+ private static final double EPS = 1e-9;
+
+ // =======================================================================
+ // 1. DPBudgetAccountant unit tests
+ // =======================================================================
+
+ @Test
+ public void testAccountantInitialisesAtZeroCost() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(1.0, 1e-5);
+ // No releases yet: total cost should be a large negative number
+ // (conversion formula gives -∞ when rdpSum = 0 for all orders),
+ // so remainingBudget() should exceed the budget.
+ assertTrue("No releases should leave budget intact",
+ acc.remainingBudget() > 0);
+ assertEquals(0, acc.releaseCount());
+ }
+
+ @Test
+ public void testSingleLaplaceReleaseDoesNotExceedBudget() {
+ // epsilon=0.5, budget=1.0: one release should consume < budget.
+ DPBudgetAccountant acc = new DPBudgetAccountant(1.0, 1e-5);
+ acc.compose(0.5, 0.0, 1.0); // Laplace, sensitivity=1
+ assertEquals(1, acc.releaseCount());
+ assertTrue("Single release within budget",
+ acc.totalEpsilonSpent() <= 1.0);
+ }
+
+ @Test
+ public void testSingleGaussianReleaseDoesNotExceedBudget() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(1.0, 1e-5);
+ acc.compose(0.5, 1e-5, 1.0); // Gaussian
+ assertEquals(1, acc.releaseCount());
+ assertTrue("Single Gaussian release within budget",
+ acc.totalEpsilonSpent() <= 1.0);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testBudgetExhaustionThrows() {
+ // Budget = 0.1, but we try to make 10 releases at epsilon=0.5 each.
+ // After enough releases the budget must be exceeded.
+ DPBudgetAccountant acc = new DPBudgetAccountant(0.1, 1e-5);
+ for (int i = 0; i < 10; i++) {
+ acc.compose(0.5, 0.0, 1.0); // will throw before the 10th
+ }
+ }
+
+ @Test
+ public void testCompositionIsMonotonicallyIncreasing() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(100.0, 1e-5); // large budget
+ double prev = acc.totalEpsilonSpent();
+ for (int i = 0; i < 5; i++) {
+ acc.compose(0.3, 1e-5, 1.0);
+ double current = acc.totalEpsilonSpent();
+ assertTrue("Epsilon spent must increase with each release",
+ current > prev);
+ prev = current;
+ }
+ }
+
+ @Test
+ public void testGaussianTighterThanLaplaceForSameEpsilon() {
+ // For the same nominal (ε, δ), Gaussian uses RDP composition which
+ // is tighter than Laplace with basic composition. After 5 releases:
+ // Laplace (basic, worst-case): 5ε
+ // Gaussian (RDP) : something < 5ε
+ double eps = 0.5;
+ double delta = 1e-5;
+
+ DPBudgetAccountant gaussian = new DPBudgetAccountant(100.0, delta);
+ DPBudgetAccountant laplace = new DPBudgetAccountant(100.0, delta);
+
+ for (int i = 0; i < 5; i++) {
+ gaussian.compose(eps, delta, 1.0);
+ laplace.compose(eps, 0.0, 1.0);
+ }
+
+ // After 5 releases, Gaussian RDP bound should be tighter.
+ // (Both may be < 5*eps; the point is Gaussian <= Laplace.)
+ assertTrue("Gaussian RDP bound should be <= Laplace bound after 5 releases",
+ gaussian.totalEpsilonSpent() <= laplace.totalEpsilonSpent() + 1e-6);
+ }
+
+ @Test
+ public void testRemainingBudgetDecreasesMonotonically() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(2.0, 1e-5);
+ double prev = acc.remainingBudget();
+ for (int i = 0; i < 3; i++) {
+ acc.compose(0.2, 1e-5, 1.0);
+ double current = acc.remainingBudget();
+ assertTrue("Remaining budget must decrease", current < prev);
+ prev = current;
+ }
+ }
+
+ @Test
+ public void testHigherEpsilonCostMoreForLaplace() {
+ // For Laplace, the accountant uses basic (pure ε-DP) composition: cost = epsilon.
+ // Sensitivity determines noise scale but NOT the budget consumed — that is set
+ // entirely by the caller's epsilon parameter.
+ // A release at epsilon=1.0 costs more budget than one at epsilon=0.5.
+ DPBudgetAccountant acc1 = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant acc2 = new DPBudgetAccountant(100.0, 1e-5);
+ acc1.compose(0.5, 0.0, 1.0); // epsilon=0.5, Laplace
+ acc2.compose(1.0, 0.0, 1.0); // epsilon=1.0, same sensitivity
+
+ assertTrue("Higher epsilon costs more budget (Laplace basic composition)",
+ acc1.totalEpsilonSpent() < acc2.totalEpsilonSpent());
+ }
+
+ // --- Constructor error paths ------------------------------------
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsZeroEpsilonBudget() {
+ new DPBudgetAccountant(0.0, 1e-5);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsNegativeEpsilonBudget() {
+ new DPBudgetAccountant(-0.5, 1e-5);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsDeltaZero() {
+ new DPBudgetAccountant(1.0, 0.0);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsDeltaOne() {
+ new DPBudgetAccountant(1.0, 1.0);
+ }
+
+ // =======================================================================
+ // 1b. DMLProgram / ExecutionContext.getDPBudgetAccountant() (dp_set_budget)
+ // =======================================================================
+ //
+ // dp_set_budget(epsilon, delta) is resolved entirely at compile time onto
+ // DMLProgram (DMLTranslator's DP_SET_BUDGET case) rather than through a
+ // runtime instruction; ExecutionContext.getDPBudgetAccountant() consults
+ // Program.getDMLProg() on first (lazy) access. These tests exercise that
+ // plumbing directly, without going through the DML compiler.
+
+ @Test
+ public void testDMLProgramHasDPBudgetTracksSetState() {
+ DMLProgram dmlProg = new DMLProgram();
+ assertFalse("No dp_set_budget call yet", dmlProg.hasDPBudget());
+ dmlProg.setDPBudget(2.0, 1e-6);
+ assertTrue("dp_set_budget was called", dmlProg.hasDPBudget());
+ assertEquals(2.0, dmlProg.getDPBudgetEpsilon(), EPS);
+ assertEquals(1e-6, dmlProg.getDPBudgetDelta(), EPS);
+ }
+
+ @Test
+ public void testGetDPBudgetAccountantUsesCompileTimeResolvedBudget() {
+ DMLProgram dmlProg = new DMLProgram();
+ dmlProg.setDPBudget(5.0, 1e-6);
+ ExecutionContext ec = ExecutionContextFactory.createContext(new Program(dmlProg));
+
+ DPBudgetAccountant acc = ec.getDPBudgetAccountant();
+ acc.compose(2.0, 0.0, 1.0); // would exceed the hardcoded default budget of 1.0
+ assertTrue("Compile-time-resolved budget should be used instead of the hardcoded default",
+ acc.remainingBudget() > 0);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testGetDPBudgetAccountantFallsBackToDefaultWithoutDPSetBudget() {
+ // No dp_set_budget call: the hardcoded default budget of epsilon=1.0 applies,
+ // so a release at epsilon=1.5 must be rejected.
+ DMLProgram dmlProg = new DMLProgram();
+ ExecutionContext ec = ExecutionContextFactory.createContext(new Program(dmlProg));
+ ec.getDPBudgetAccountant().compose(1.5, 0.0, 1.0);
+ }
+
+ @Test
+ public void testGetDPBudgetAccountantIsLazyAndCachedPerContext() {
+ // The accountant must be created once and reused across calls on the
+ // same ExecutionContext, not rebuilt (which would reset releaseCount()).
+ DMLProgram dmlProg = new DMLProgram();
+ dmlProg.setDPBudget(10.0, 1e-6);
+ ExecutionContext ec = ExecutionContextFactory.createContext(new Program(dmlProg));
+
+ ec.getDPBudgetAccountant().compose(1.0, 0.0, 1.0);
+ assertEquals("Same accountant instance must be reused across calls",
+ 1, ec.getDPBudgetAccountant().releaseCount());
+ }
+
+ // --- Single-argument convenience constructor -------------------
+
+ @Test
+ public void testConvenienceConstructorDefaultsDeltaTo1e5() {
+ // The one-arg form delegates to (epsilonBudget, 1e-5). A Gaussian
+ // release whose per-release delta matches that default must produce
+ // identical totalEpsilonSpent() from both construction paths.
+ DPBudgetAccountant oneArg = new DPBudgetAccountant(10.0);
+ DPBudgetAccountant twoArg = new DPBudgetAccountant(10.0, 1e-5);
+ oneArg.compose(0.5, 1e-5, 1.0);
+ twoArg.compose(0.5, 1e-5, 1.0);
+ assertEquals("Convenience constructor must default to delta=1e-5",
+ twoArg.totalEpsilonSpent(), oneArg.totalEpsilonSpent(), EPS);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testGaussianBudgetExhaustionThrows() {
+ // Budget = 0.1. Each Gaussian release costs more than 0.005, so 20
+ // releases must exceed the budget well before the loop ends.
+ DPBudgetAccountant acc = new DPBudgetAccountant(0.1, 1e-5);
+ for (int i = 0; i < 20; i++) {
+ acc.compose(0.3, 1e-5, 1.0);
+ }
+ }
+
+ @Test
+ public void testMixedCompositionExceedsEitherAlone() {
+ // Compose one Laplace and one Gaussian release. The total cost must
+ // exceed what either mechanism contributes alone, exercising the
+ // _pureEpsilonSum + gaussianEps addition path in totalEpsilonSpent().
+ DPBudgetAccountant mixed = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant lapOnly = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant gauOnly = new DPBudgetAccountant(100.0, 1e-5);
+
+ mixed.compose(0.5, 0.0, 1.0); // Laplace
+ mixed.compose(0.5, 1e-5, 1.0); // Gaussian
+
+ lapOnly.compose(0.5, 0.0, 1.0);
+ gauOnly.compose(0.5, 1e-5, 1.0);
+
+ assertTrue("Mixed cost must exceed Laplace-only cost",
+ mixed.totalEpsilonSpent() > lapOnly.totalEpsilonSpent());
+ assertTrue("Mixed cost must exceed Gaussian-only cost",
+ mixed.totalEpsilonSpent() > gauOnly.totalEpsilonSpent());
+ }
+
+ // --- Release count across multiple mixed releases --------------
+
+ @Test
+ public void testReleaseCountTracksAllReleases() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(100.0, 1e-5);
+ assertEquals(0, acc.releaseCount());
+ acc.compose(0.1, 0.0, 1.0); // Laplace
+ assertEquals(1, acc.releaseCount());
+ acc.compose(0.1, 1e-5, 1.0); // Gaussian
+ assertEquals(2, acc.releaseCount());
+ acc.compose(0.1, 0.0, 1.0); // Laplace
+ acc.compose(0.1, 0.0, 1.0); // Laplace
+ acc.compose(0.1, 1e-5, 1.0); // Gaussian
+ assertEquals(5, acc.releaseCount());
+ }
+
+ // --- Edge-case inputs for rdpGaussian / gaussianSigma ----------
+
+ @Test
+ public void testGaussianSensitivityCancelsInRDP() {
+ // For the Gaussian mechanism: σ = Δf·sqrt(2·ln(1.25/δ))/ε, so
+ // D_α = α·Δf²/(2σ²) = α·ε²/(4·ln(1.25/δ)).
+ // Sensitivity cancels. Two accountants with the same (ε,δ) but
+ // different sensitivity must report identical totalEpsilonSpent().
+ DPBudgetAccountant acc1 = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant acc2 = new DPBudgetAccountant(100.0, 1e-5);
+ acc1.compose(0.5, 1e-5, 1.0); // sensitivity = 1
+ acc2.compose(0.5, 1e-5, 100.0); // sensitivity = 100, same (ε,δ)
+ assertEquals("Gaussian RDP cost must be independent of sensitivity when (ε,δ) are fixed",
+ acc1.totalEpsilonSpent(), acc2.totalEpsilonSpent(), EPS);
+ }
+
+ @Test
+ public void testGaussianLargerEpsilonCostsMoreBudget() {
+ // D_α ∝ ε², so a release declared at a higher ε (less noise, more
+ // privacy loss) must cost more budget than one at a lower ε.
+ DPBudgetAccountant lowEps = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant highEps = new DPBudgetAccountant(100.0, 1e-5);
+ lowEps.compose(0.1, 1e-5, 1.0);
+ highEps.compose(0.5, 1e-5, 1.0);
+ assertTrue("Larger epsilon per Gaussian release must cost more budget",
+ highEps.totalEpsilonSpent() > lowEps.totalEpsilonSpent());
+ }
+
+ // =======================================================================
+ // 2. Noise distribution tests (statistical sanity checks)
+ // =======================================================================
+ // These tests generate many samples and verify that the empirical mean
+ // is near zero and the empirical variance matches the theoretical value
+ // within a reasonable tolerance.
+ //
+
+ @Test
+ public void testLaplaceNoiseMeanNearZero() {
+ // For 10000 samples the empirical mean should be within 3σ/√n of 0.
+ int n = 10_000;
+ double scale = 2.0;
+ double[] samples = sampleLaplace(n, scale);
+ double mean = mean(samples);
+ double theoreticalStdErr = scale * Math.sqrt(2.0) / Math.sqrt(n);
+ assertTrue("Laplace mean should be near 0",
+ Math.abs(mean) < 5 * theoreticalStdErr);
+ }
+
+ @Test
+ public void testLaplaceNoiseVarianceCorrect() {
+ // Var[Laplace(0, b)] = 2b². Allow 10% relative error for n=10000.
+ int n = 10_000;
+ double scale = 1.5;
+ double[] samples = sampleLaplace(n, scale);
+ double variance = variance(samples);
+ double expected = 2.0 * scale * scale;
+ assertEquals("Laplace variance", expected, variance, 0.1 * expected);
+ }
+
+ @Test
+ public void testGaussianNoiseMeanNearZero() {
+ int n = 10_000;
+ double sigma = 3.0;
+ double[] samples = sampleGaussian(n, sigma);
+ double mean = mean(samples);
+ double theoreticalStdErr = sigma / Math.sqrt(n);
+ assertTrue("Gaussian mean should be near 0",
+ Math.abs(mean) < 5 * theoreticalStdErr);
+ }
+
+ @Test
+ public void testGaussianNoiseVarianceCorrect() {
+ int n = 10_000;
+ double sigma = 2.0;
+ double[] samples = sampleGaussian(n, sigma);
+ double variance = variance(samples);
+ double expected = sigma * sigma;
+ assertEquals("Gaussian variance", expected, variance, 0.1 * expected);
+ }
+
+ // -----------------------------------------------------------------------
+ // Helpers for noise distribution tests
+ // -----------------------------------------------------------------------
+
+ /** Sample n Laplace(0, scale) values using the inverse-CDF method. */
+ private static double[] sampleLaplace(int n, double scale) {
+ java.util.concurrent.ThreadLocalRandom rng =
+ java.util.concurrent.ThreadLocalRandom.current();
+ double[] out = new double[n];
+ for (int i = 0; i < n; i++) {
+ double u = rng.nextDouble();
+ double v = u - 0.5;
+ if (v == 0.0) v = 1e-15;
+ out[i] = -scale * Math.signum(v) * Math.log(1.0 - 2.0 * Math.abs(v));
+ }
+ return out;
+ }
+
+ /** Sample n N(0, sigma²) values. */
+ private static double[] sampleGaussian(int n, double sigma) {
+ java.util.concurrent.ThreadLocalRandom rng =
+ java.util.concurrent.ThreadLocalRandom.current();
+ double[] out = new double[n];
+ for (int i = 0; i < n; i++) {
+ out[i] = sigma * rng.nextGaussian();
+ }
+ return out;
+ }
+
+ private static double mean(double[] xs) {
+ double s = 0;
+ for (double x : xs) s += x;
+ return s / xs.length;
+ }
+
+ private static double variance(double[] xs) {
+ double m = mean(xs);
+ double s = 0;
+ for (double x : xs) s += (x - m) * (x - m);
+ return s / (xs.length - 1);
+ }
+}
diff --git a/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java b/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java
new file mode 100644
index 00000000000..9b2ef30e0af
--- /dev/null
+++ b/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java
@@ -0,0 +1,260 @@
+// ==========================================================================
+// DML integration test
+// ==========================================================================
+//
+// Full integration tests extend AutomatedTestBase and drive the DML runner.
+// Each test:
+// (a) Writes a DML script to a temp file.
+// (b) Provides input matrices via TestUtils.
+// (c) Calls runTest() and reads the output MatrixBlock.
+// (d) Verifies that the noisy result differs from the clean result by a
+// statistically plausible amount (not zero, not astronomically large).
+//
+
+
+package org.apache.sysds.test.functions.privacy.dp;
+
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
+
+import java.io.File;
+import java.io.IOException;
+import java.nio.file.Files;
+import java.util.HashMap;
+
+import org.apache.sysds.parser.LanguageException;
+import org.apache.sysds.parser.ParseException;
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex;
+import org.apache.sysds.test.AutomatedTestBase;
+import org.apache.sysds.test.TestConfiguration;
+import org.apache.sysds.test.TestUtils;
+import org.junit.Test;
+
+public class DPBuiltinDMLTest extends AutomatedTestBase {
+
+ private static final String TEST_DIR = "functions/privacy/dp/";
+ private static final String TEST_CLASS = TEST_DIR + DPBuiltinDMLTest.class.getSimpleName() + "/";
+ private static final int ROWS = 100;
+ private static final int COLS = 10;
+
+ private static final String DML_LAPLACE_TEMPLATE =
+ "X = read($1);\n"
+ + "result = dp_laplace(X, query=\"%s\", sensitivity=1.0, epsilon=$2);\n"
+ + "write(result, $3, format=\"text\");\n";
+ private static final String DML_GAUSSIAN_TEMPLATE =
+ "X = read($1);\n"
+ + "result = dp_gaussian(X, query=\"%s\", sensitivity=1.0, epsilon=$2, delta=1e-5);\n"
+ + "write(result, $3, format=\"text\");\n";
+
+ private static final String DML_LAPLACE = String.format(DML_LAPLACE_TEMPLATE, "colMeans");
+ private static final String DML_GAUSSIAN = String.format(DML_GAUSSIAN_TEMPLATE, "colMeans");
+
+ // dp_set_budget(epsilon, delta) is resolved entirely at compile time (its arguments
+ // must be literals), called via a dummy assignment, then a single dp_laplace release
+ // at $2 records a cost of exactly $2 (Laplace basic composition).
+ private static final String DML_SET_BUDGET_TEMPLATE =
+ "eps = dp_set_budget(%s, 1e-6);\n"
+ + "X = read($1);\n"
+ + "result = dp_laplace(X, query=\"colMeans\", sensitivity=1.0, epsilon=$2);\n"
+ + "write(result, $3, format=\"text\");\n";
+
+ // Two dp_set_budget calls in the same script; DMLTranslator must reject this at
+ // compile time (DMLProgram.hasDPBudget()).
+ private static final String DML_SET_BUDGET_TWICE =
+ "eps = dp_set_budget(3.0, 1e-6);\n"
+ + "eps2 = dp_set_budget(5.0, 1e-6);\n"
+ + "X = read($1);\n"
+ + "result = dp_laplace(X, query=\"colMeans\", sensitivity=1.0, epsilon=$2);\n"
+ + "write(result, $3, format=\"text\");\n";
+
+ // A budget argument computed at runtime (not a literal); must be rejected at
+ // compile time by BuiltinFunctionExpression's isConstant() check.
+ private static final String DML_SET_BUDGET_NON_LITERAL =
+ "X = read($1);\n"
+ + "computed = sum(X) / nrow(X);\n"
+ + "eps = dp_set_budget(computed, 1e-6);\n"
+ + "result = dp_laplace(X, query=\"colMeans\", sensitivity=1.0, epsilon=$2);\n"
+ + "write(result, $3, format=\"text\");\n";
+
+ @Override
+ public void setUp() {
+ addTestConfiguration("DPLaplace", new TestConfiguration(TEST_CLASS, "DPLaplace"));
+ addTestConfiguration("DPGaussian", new TestConfiguration(TEST_CLASS, "DPGaussian"));
+ addTestConfiguration("DPSetBudget", new TestConfiguration(TEST_CLASS, "DPSetBudget"));
+ }
+
+ @Test
+ public void testLaplaceOutputDiffersFromCleanMean() {
+ runColMeansDPTest("DPLaplace", DML_LAPLACE, "0.5");
+ }
+
+ @Test
+ public void testGaussianOutputDiffersFromCleanMean() {
+ runColMeansDPTest("DPGaussian", DML_GAUSSIAN, "0.5");
+ }
+
+ @Test
+ public void testLaplaceColSums() {
+ // query="colSums": T is 1 x n filled with 1.0, output is the noisy column-sum row vector.
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ HashMap result = runAndGetResult("DPLaplace",
+ String.format(DML_LAPLACE_TEMPLATE, "colSums"), "0.5", data);
+ assertShape(result, 1, COLS);
+ double maxDiff = maxAbsDiffFromClean(data, result, DPBuiltinDMLTest::colSum);
+ assertTrue("Result should differ from the clean column sums", maxDiff > 0);
+ }
+
+ @Test
+ public void testGaussianIdentity() {
+ // query="identity": T is the n x n identity, output is a noisy release of X itself.
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ HashMap result = runAndGetResult("DPGaussian",
+ String.format(DML_GAUSSIAN_TEMPLATE, "identity"), "0.5", data);
+ assertShape(result, ROWS, COLS);
+ // identity releases X row-by-row, so compare cell-by-cell rather than via a per-column reduction.
+ double maxCellDiff = 0;
+ for (int r = 0; r < ROWS; r++) {
+ for (int c = 0; c < COLS; c++) {
+ double noisy = result.get(new CellIndex(r + 1, c + 1));
+ maxCellDiff = Math.max(maxCellDiff, Math.abs(noisy - data[r][c]));
+ }
+ }
+ assertTrue("Result should differ from the clean matrix", maxCellDiff > 0);
+ }
+
+ @Test
+ public void testHighEpsilonIsCloserToTruth() {
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ // Higher ε → less noise → result closer to the true mean.
+ // NOTE: the DPBudgetAccountant caps total spend at the default budget
+ // (ε = 1.0) regardless of the per-release ε requested, so ε values
+ // here must stay well under that cap or the release is rejected.
+ double noisyLow = runAndGetMaxAbsColMeansDiffFromClean(data, "DPGaussian", DML_GAUSSIAN, "0.1");
+ double noisyHigh = runAndGetMaxAbsColMeansDiffFromClean(data, "DPGaussian", DML_GAUSSIAN, "0.5");
+ assertTrue("ε=0.5 should give less noise than ε=0.1", noisyHigh < noisyLow);
+ }
+
+ @Test
+ public void testSetBudgetLiteralAllowsExceedingDefaultBudget() {
+ // Default budget is epsilon=1.0; a single release at epsilon=1.5 would be
+ // rejected unless dp_set_budget(3.0, ...) widens it first.
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ HashMap result = runAndGetResult("DPSetBudget",
+ String.format(DML_SET_BUDGET_TEMPLATE, "3.0"), "1.5", data);
+ assertShape(result, 1, COLS);
+ }
+
+ @Test
+ public void testSetBudgetNarrowBudgetStillEnforced() {
+ // An explicit narrow budget must still be enforced: epsilon=0.8 exceeds
+ // the explicit budget of 0.5.
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ runExpectingException("DPSetBudget", String.format(DML_SET_BUDGET_TEMPLATE, "0.5"), "0.8", data,
+ DMLRuntimeException.class);
+ }
+
+ @Test
+ public void testSetBudgetCalledTwiceFailsAtCompileTime() {
+ // Thrown from DMLTranslator.processBuiltinFunctionExpression (HOP construction),
+ // which wraps all case-block exceptions in ParseException (see processExpression's
+ // catch-all) — unlike the non-literal check below, which runs during validation
+ // and so surfaces as a bare LanguageException.
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ runExpectingException("DPSetBudget", DML_SET_BUDGET_TWICE, "0.5", data, ParseException.class);
+ }
+
+ @Test
+ public void testSetBudgetRejectsNonLiteralArgs() {
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ runExpectingException("DPSetBudget", DML_SET_BUDGET_NON_LITERAL, "0.5", data, LanguageException.class);
+ }
+
+ private void runColMeansDPTest(String testName, String dml, String epsilonStr) {
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ HashMap result = runAndGetResult(testName, dml, epsilonStr, data);
+ assertShape(result, 1, COLS);
+ // Must differ from the exact (clean) mean by a non-trivial amount.
+ // (A single-seed exact-equality check is fragile; use range check.)
+ double maxDiff = maxAbsDiffFromClean(data, result, DPBuiltinDMLTest::colMean);
+ assertTrue("Result should differ from the clean mean", maxDiff > 0);
+ }
+
+ private double runAndGetMaxAbsColMeansDiffFromClean(double[][] data, String testName, String dml, String epsilonStr) {
+ HashMap result = runAndGetResult(testName, dml, epsilonStr, data);
+ return maxAbsDiffFromClean(data, result, DPBuiltinDMLTest::colMean);
+ }
+
+ private static void assertShape(HashMap result, int expectedRows, int expectedCols) {
+ int maxRow = 0, maxCol = 0;
+ for (CellIndex ci : result.keySet()) {
+ maxRow = Math.max(maxRow, ci.row);
+ maxCol = Math.max(maxCol, ci.column);
+ }
+ assertEquals("Result should have " + expectedRows + " row(s)", expectedRows, maxRow);
+ assertEquals("Result should have " + expectedCols + " column(s)", expectedCols, maxCol);
+ }
+
+ @FunctionalInterface
+ private interface CleanColumnFn {
+ double apply(double[][] data, int col);
+ }
+
+ /** Computes max|noisy(1,c) - clean(data,c)| across the (1 x COLS) row-vector releases. */
+ private static double maxAbsDiffFromClean(double[][] data, HashMap result,
+ CleanColumnFn cleanFn) {
+ double maxDiff = 0;
+ for (int c = 0; c < COLS; c++) {
+ double clean = cleanFn.apply(data, c);
+ double noisy = result.get(new CellIndex(1, c + 1));
+ maxDiff = Math.max(maxDiff, Math.abs(noisy - clean));
+ }
+ return maxDiff;
+ }
+
+ private static double colMean(double[][] data, int c) {
+ double sum = 0;
+ for (int r = 0; r < ROWS; r++)
+ sum += data[r][c];
+ return sum / ROWS;
+ }
+
+ private static double colSum(double[][] data, int c) {
+ double sum = 0;
+ for (int r = 0; r < ROWS; r++)
+ sum += data[r][c];
+ return sum;
+ }
+
+ private HashMap runAndGetResult(String testName, String dml, String epsilonStr,
+ double[][] data)
+ {
+ prepareScript(testName, dml, epsilonStr, data);
+ runTest(true, false, null, -1);
+ return readDMLMatrixFromOutputDir("result");
+ }
+
+ private void runExpectingException(String testName, String dml, String epsilonStr, double[][] data,
+ Class expectedException)
+ {
+ prepareScript(testName, dml, epsilonStr, data);
+ runTest(true, true, expectedException, -1);
+ }
+
+ private void prepareScript(String testName, String dml, String epsilonStr, double[][] data) {
+ getAndLoadTestConfiguration(testName);
+ writeInputMatrixWithMTD("X", data, false);
+
+ fullDMLScriptName = getScript();
+ try {
+ File scriptFile = new File(fullDMLScriptName);
+ scriptFile.getParentFile().mkdirs();
+ Files.write(scriptFile.toPath(), dml.getBytes());
+ }
+ catch (IOException e) {
+ throw new RuntimeException(e);
+ }
+
+ programArgs = new String[]{ "-args", input("X"), epsilonStr, output("result") };
+ }
+}