|
| 1 | +#!/usr/bin/env python3 |
| 2 | +""" |
| 3 | +scripts/optimise_vqc_params.py |
| 4 | +Optimise VQC parameters using fully vectorised numpy simulation. |
| 5 | +
|
| 6 | +Runs all 100 events simultaneously using numpy matrix operations — |
| 7 | +no CUDA-Q dependency, completes in ~5-15 minutes. |
| 8 | +
|
| 9 | +Usage: |
| 10 | + python3 scripts/optimise_vqc_params.py [--trials 8] [--maxiter 2000] |
| 11 | +""" |
| 12 | + |
| 13 | +import sys |
| 14 | +import json |
| 15 | +import ast |
| 16 | +import time |
| 17 | +import argparse |
| 18 | +import numpy as np |
| 19 | +from pathlib import Path |
| 20 | +from scipy.optimize import minimize |
| 21 | + |
| 22 | +_ROOT = Path(__file__).resolve().parent.parent |
| 23 | +_SRC = _ROOT / "src" |
| 24 | +sys.path.insert(0, str(_ROOT)) |
| 25 | +sys.path.insert(0, str(_SRC)) |
| 26 | + |
| 27 | +from quantum.circuit import compress_features, FEATURE_DIM, N_QUBITS, N_LAYERS |
| 28 | + |
| 29 | + |
| 30 | +def _apply_ry_qubit(states: np.ndarray, qubit: int, angles: np.ndarray) -> np.ndarray: |
| 31 | + """Apply Ry(angle[i]) to qubit for batch of states. states: (N, 2**N_QUBITS)""" |
| 32 | + n_states = 2 ** N_QUBITS |
| 33 | + new_states = states.copy() |
| 34 | + stride = 2 ** qubit |
| 35 | + c = np.cos(angles / 2) |
| 36 | + s = np.sin(angles / 2) |
| 37 | + for i in range(0, n_states, stride * 2): |
| 38 | + for j in range(stride): |
| 39 | + idx0 = i + j |
| 40 | + idx1 = i + j + stride |
| 41 | + a = states[:, idx0].copy() |
| 42 | + b = states[:, idx1].copy() |
| 43 | + new_states[:, idx0] = c * a - s * b |
| 44 | + new_states[:, idx1] = s * a + c * b |
| 45 | + return new_states |
| 46 | + |
| 47 | + |
| 48 | +def _apply_cx(states: np.ndarray, control: int, target: int) -> np.ndarray: |
| 49 | + """Apply CNOT(control→target) to batch of states.""" |
| 50 | + n_states = 2 ** N_QUBITS |
| 51 | + new_states = states.copy() |
| 52 | + for idx in range(n_states): |
| 53 | + if (idx >> control) & 1: |
| 54 | + flipped = idx ^ (1 << target) |
| 55 | + if idx < flipped: |
| 56 | + new_states[:, idx], new_states[:, flipped] = \ |
| 57 | + states[:, flipped].copy(), states[:, idx].copy() |
| 58 | + return new_states |
| 59 | + |
| 60 | + |
| 61 | +def vqc_batch(feats: np.ndarray, params: np.ndarray) -> np.ndarray: |
| 62 | + """ |
| 63 | + Vectorised VQC for N events simultaneously. |
| 64 | + feats: (N, N_QUBITS) compressed features |
| 65 | + params: (N_PARAMS,) VQC rotation angles |
| 66 | + Returns: (N,) confidence values — high = anomaly |
| 67 | + """ |
| 68 | + N = len(feats) |
| 69 | + n_states = 2 ** N_QUBITS |
| 70 | + states = np.zeros((N, n_states), dtype=complex) |
| 71 | + states[:, 0] = 1.0 |
| 72 | + |
| 73 | + # Feature encoding |
| 74 | + for i in range(N_QUBITS): |
| 75 | + states = _apply_ry_qubit(states, i, feats[:, i] * np.pi) |
| 76 | + |
| 77 | + # Ansatz layers |
| 78 | + param_idx = 0 |
| 79 | + for _ in range(N_LAYERS): |
| 80 | + for i in range(N_QUBITS - 1): |
| 81 | + states = _apply_cx(states, i, i + 1) |
| 82 | + for i in range(N_QUBITS): |
| 83 | + angle = params[param_idx % len(params)] |
| 84 | + states = _apply_ry_qubit(states, i, np.full(N, angle)) |
| 85 | + param_idx += 1 |
| 86 | + |
| 87 | + # P(zero state) → confidence = 1 - P(zero) |
| 88 | + return 1.0 - np.abs(states[:, 0]) ** 2 |
| 89 | + |
| 90 | + |
| 91 | +def main(): |
| 92 | + parser = argparse.ArgumentParser() |
| 93 | + parser.add_argument("--trials", type=int, default=8) |
| 94 | + parser.add_argument("--maxiter", type=int, default=2000) |
| 95 | + parser.add_argument("--seed", type=int, default=42) |
| 96 | + args = parser.parse_args() |
| 97 | + |
| 98 | + print(f"N_QUBITS={N_QUBITS} N_LAYERS={N_LAYERS} FEATURE_DIM={FEATURE_DIM}") |
| 99 | + |
| 100 | + demo_events = json.load(open(_ROOT / "data/samples/demo_events.json")) |
| 101 | + demo_labels = json.load(open(_ROOT / "data/samples/demo_labels.json")) |
| 102 | + |
| 103 | + feats, gt = [], [] |
| 104 | + for e in demo_events: |
| 105 | + fv = e["feature_vector"] |
| 106 | + if isinstance(fv, str): |
| 107 | + fv = ast.literal_eval(fv) |
| 108 | + if len(fv) < FEATURE_DIM: |
| 109 | + fv = fv + [0.0] * (FEATURE_DIM - len(fv)) |
| 110 | + feats.append(compress_features(fv)) |
| 111 | + gt.append(demo_labels[e["id"]]["anomaly"]) |
| 112 | + |
| 113 | + feats = np.array(feats, dtype=np.float64) |
| 114 | + gt = np.array(gt) |
| 115 | + N_PARAMS = N_QUBITS * N_LAYERS * 2 |
| 116 | + print(f"Events: {len(feats)} ({int(gt.sum())} anomalous, {int((1-gt).sum())} benign)") |
| 117 | + print(f"Params: {N_PARAMS}") |
| 118 | + |
| 119 | + # Warm-up timing |
| 120 | + print("\nWarm-up (5 runs)...") |
| 121 | + rng = np.random.default_rng(args.seed) |
| 122 | + x0 = rng.uniform(-np.pi, np.pi, N_PARAMS) |
| 123 | + t0 = time.monotonic() |
| 124 | + for _ in range(5): |
| 125 | + vqc_batch(feats, x0) |
| 126 | + warmup = (time.monotonic() - t0) / 5 |
| 127 | + print(f" {warmup*1000:.0f}ms per evaluation") |
| 128 | + est_total = warmup * args.maxiter * args.trials / 60 |
| 129 | + print(f" Estimated: ~{est_total:.0f} min total ({args.trials} trials × {args.maxiter} iters)") |
| 130 | + |
| 131 | + best_result = None |
| 132 | + best_loss = float("inf") |
| 133 | + call_count = [0] |
| 134 | + start = time.monotonic() |
| 135 | + |
| 136 | + def loss(params): |
| 137 | + call_count[0] += 1 |
| 138 | + confs = vqc_batch(feats, params) |
| 139 | + mean_anom = float(confs[gt == 1].mean()) |
| 140 | + mean_ben = float(confs[gt == 0].mean()) |
| 141 | + sep = mean_anom - mean_ben |
| 142 | + var = float(confs[gt == 1].std() + confs[gt == 0].std()) |
| 143 | + obj = -sep + 0.3 * var |
| 144 | + if call_count[0] % 100 == 0: |
| 145 | + acc = float(((confs >= 0.5).astype(int) == gt).mean()) |
| 146 | + elapsed = time.monotonic() - start |
| 147 | + print(f" eval {call_count[0]:5d} sep={sep:.3f} " |
| 148 | + f"acc={acc:.0%} loss={obj:.4f} {elapsed:.0f}s") |
| 149 | + return float(obj) |
| 150 | + |
| 151 | + for trial in range(args.trials): |
| 152 | + print(f"\nTrial {trial + 1}/{args.trials}") |
| 153 | + x0 = rng.uniform(-np.pi, np.pi, N_PARAMS) |
| 154 | + call_count = [0] |
| 155 | + start = time.monotonic() |
| 156 | + |
| 157 | + result = minimize( |
| 158 | + loss, x0, |
| 159 | + method = "Nelder-Mead", |
| 160 | + options = {"maxiter": args.maxiter, "xatol": 0.003, |
| 161 | + "fatol": 0.001, "disp": False}, |
| 162 | + ) |
| 163 | + elapsed = time.monotonic() - start |
| 164 | + print(f" {result.message} loss={result.fun:.4f} " |
| 165 | + f"evals={result.nfev} {elapsed:.0f}s") |
| 166 | + |
| 167 | + if result.fun < best_loss: |
| 168 | + best_loss = result.fun |
| 169 | + best_result = result |
| 170 | + # Save best params after each improvement |
| 171 | + json.dump(list(result.x.tolist()), |
| 172 | + open(_ROOT / "data/vqc_params_v2.json", "w"), indent=2) |
| 173 | + print(f" ✓ New best — saved") |
| 174 | + |
| 175 | + print(f"\n{'='*60}") |
| 176 | + print(f"Best loss: {best_loss:.4f}") |
| 177 | + |
| 178 | + confs = vqc_batch(feats, best_result.x) |
| 179 | + preds = (confs >= 0.5).astype(int) |
| 180 | + tp = int(((preds==1)&(gt==1)).sum()) |
| 181 | + fp = int(((preds==1)&(gt==0)).sum()) |
| 182 | + fn = int(((preds==0)&(gt==1)).sum()) |
| 183 | + tn = int(((preds==0)&(gt==0)).sum()) |
| 184 | + f1 = 2*tp/(2*tp+fp+fn) if (2*tp+fp+fn)>0 else 0 |
| 185 | + |
| 186 | + print(f"TP={tp} FP={fp} FN={fn} TN={tn} F1={f1:.3f} Acc={(tp+tn)/100:.0%}") |
| 187 | + print(f"Anomaly conf: mean={confs[gt==1].mean():.3f} std={confs[gt==1].std():.3f}") |
| 188 | + print(f"Benign conf: mean={confs[gt==0].mean():.3f} std={confs[gt==0].std():.3f}") |
| 189 | + print(f"Separation: {confs[gt==1].mean() - confs[gt==0].mean():.3f}") |
| 190 | + print(f"\nSaved to data/vqc_params_v2.json") |
| 191 | + |
| 192 | + |
| 193 | +if __name__ == "__main__": |
| 194 | + main() |
0 commit comments