import time
import numpy as np
import pandas as pd
from digiqual.pod import (
fit_all_robust_mean_models,
fit_variance_model,
infer_best_distribution,
bootstrap_pod_ci
)
from digiqual.integration import compute_multi_dim_pod
import digiqual.cpp_fallback as cpp_fb
def run_performance_comparison():
# 1. Generate Large Synthetic Physics Dataset (N=300 samples)
np.random.seed(42)
N = 300
length = np.random.uniform(0.5, 10.0, N)
angle = np.random.uniform(-45.0, 45.0, N)
roughness = np.random.uniform(0.0, 1.0, N)
noise_scale = 0.5 + 0.4 * length + 1.0 * roughness
noise = np.random.gumbel(loc=0, scale=noise_scale) - noise_scale * 0.57721
signal = 5.0 + 3.0*length - 0.8*(length**2) + 0.1*(length**3) + 0.1*angle - 5.0*roughness + noise
X = np.column_stack([length, angle, roughness])
y = signal
feature_names = ['Length', 'Angle', 'Roughness']
poi_names = ['Length', 'Angle']
print("Dataset prepared: N=300 samples, 3D inputs")
# 2. Model Fitting & Noise Distribution Inference
models, scores, best_key = fit_all_robust_mean_models(X, y, max_degree=4)
mean_model = models[best_key]
residuals, bw = fit_variance_model(X, y, mean_model)
dist_info = infer_best_distribution(residuals, X, bw)
# 3. Define Evaluation Grid (100 PoI points) & Active Nuisance Bounds
poi_grid = np.column_stack([np.linspace(1, 9, 100), np.linspace(-30, 30, 100)])
nuisance_ranges = {'Roughness': (0.0, 1.0)}
# =========================================================
# WORKFLOW A: C++ Accelerated Backend
# =========================================================
cpp_fb.HAS_CPP = True
print("\n--- Running Workflow A (C++ Accelerated) ---")
t0 = time.perf_counter()
pod_mc_cpp, _ = compute_multi_dim_pod(
poi_grid, nuisance_ranges, mean_model, X, residuals, bw, dist_info,
thresholds=20.0, n_mc_samples=1000, feature_names=feature_names, poi_names=poi_names
)
t_mc_cpp = time.perf_counter() - t0
print(f"[C++] Monte Carlo Marginalisation (100 pts x 1,000 MC): {t_mc_cpp:.4f} s")
t0 = time.perf_counter()
boot_cpp = bootstrap_pod_ci(
X, y, poi_grid, threshold=20.0,
model_type=getattr(mean_model, 'model_type_', 'Polynomial'),
model_params=getattr(mean_model, 'model_params_', 3),
bandwidth=bw, dist_info=dist_info, n_boot=100,
nuisance_ranges=nuisance_ranges, n_jobs=1,
feature_names=feature_names, poi_names=poi_names
)
t_boot_cpp = time.perf_counter() - t0
print(f"[C++] Bootstrap Resampling (100 resamples x 100 pts): {t_boot_cpp:.4f} s")
# =========================================================
# WORKFLOW B: Python-Only Fallback
# =========================================================
cpp_fb.HAS_CPP = False
print("\n--- Running Workflow B (Python-Only Fallback) ---")
t0 = time.perf_counter()
pod_mc_py, _ = compute_multi_dim_pod(
poi_grid, nuisance_ranges, mean_model, X, residuals, bw, dist_info,
thresholds=20.0, n_mc_samples=1000, feature_names=feature_names, poi_names=poi_names
)
t_mc_py = time.perf_counter() - t0
print(f"[Python] Monte Carlo Marginalisation (100 pts x 1,000 MC): {t_mc_py:.4f} s")
t0 = time.perf_counter()
boot_py = bootstrap_pod_ci(
X, y, poi_grid, threshold=20.0,
model_type=getattr(mean_model, 'model_type_', 'Polynomial'),
model_params=getattr(mean_model, 'model_params_', 3),
bandwidth=bw, dist_info=dist_info, n_boot=100,
nuisance_ranges=nuisance_ranges, n_jobs=1,
feature_names=feature_names, poi_names=poi_names
)
t_boot_py = time.perf_counter() - t0
print(f"[Python] Bootstrap Resampling (100 resamples x 100 pts): {t_boot_py:.4f} s")
# Restore C++ status
cpp_fb.HAS_CPP = True
# =========================================================
# SPEEDUP SUMMARY
# =========================================================
print("\n======================================================")
print("SPEEDUP RESULTS")
print("======================================================")
print(f"MC Marginalisation Speedup: {t_mc_py / t_mc_cpp:.2f}x Faster")
print(f"Bootstrap Resampling Speedup: {t_boot_py / t_boot_cpp:.2f}x Faster")
print("======================================================")
if __name__ == "__main__":
run_performance_comparison()