Rishit/wf512

benchmark_driver.py

@@ -0,0 +1,304 @@
+import json
+import subprocess
+import argparse
+import sys
+import os
+import time
+import csv
+from datetime import datetime
+
+def load_config(config_path):
+    with open(config_path, 'r') as f:
+        return json.load(f)
+
+def compile_rebound(profile=False):
+    print(f"Compiling REBOUND (Profile mode: {profile})...")
+    opt_flags = "-march=native -O3"
+    if profile:
+        opt_flags += " -DPROF"
+
+    os.environ['OPT'] = opt_flags
+
+    try:
+        subprocess.check_call(['make', 'clean'], cwd='examples/whfast512_solar_system_jac', stdout=subprocess.DEVNULL)
+        subprocess.check_call(['make'], cwd='examples/whfast512_solar_system_jac', stdout=subprocess.DEVNULL)
+        print("Compilation successful.")
+        return True
+    except subprocess.CalledProcessError:
+        print("Compilation failed!")
+        return False
+
+def run_benchmark(config, executable='./rebound'):
+    cmd = [executable]
+    for key, value in config['args'].items():
+        cmd.append(f"--{key}")
+        cmd.append(str(value))
+
+    try:
+        result = subprocess.run(
+            cmd, 
+            cwd='examples/whfast512_solar_system_jac',
+            capture_output=True,
+            text=True,
+            check=True
+        )
+
+        profile_data = None
+        if "PROFILING_START" in result.stdout:
+            # Extract and echo profiling block
+            print("\n" + "-"*40)
+            print(f"Profiling Output for {config['name']}:")
+            start = result.stdout.find("PROFILING_START")
+            end = result.stdout.find("PROFILING_END") + 13
+            profiling_block = result.stdout[start:end]
+            print(profiling_block)
+            print("-"*40 + "\n")
+
+            # Parse profiling numbers for CSV when running in profiled mode
+            profile_data = {}
+            lines = profiling_block.strip().splitlines()
+            for line in lines:
+                line = line.strip()
+                if not line or line.startswith("PROFILING_"):
+                    continue
+                # Lines are of the form "Label:  value"
+                if ":" not in line:
+                    continue
+                label, val_str = line.split(":", 1)
+                label = label.strip()
+                val_str = val_str.strip()
+                try:
+                    val = float(val_str)
+                except ValueError:
+                    continue
+
+                # Map human-readable labels to stable CSV keys
+                if label.startswith("Total Walltime"):
+                    profile_data["prof_total_walltime"] = val
+                elif label.startswith("Kepler"):
+                    profile_data["prof_kepler"] = val
+                elif "Stiefel" in label:
+                    profile_data["prof_kepler_stiefel"] = val
+                elif "f/g func" in label:
+                    profile_data["prof_kepler_fg"] = val
+                elif label.startswith("Interaction"):
+                    profile_data["prof_interaction"] = val
+                elif label.startswith("Forces"):
+                    profile_data["prof_forces"] = val
+                elif "Loop 1" in label:
+                    profile_data["prof_loop1"] = val
+                elif "Loop 2" in label:
+                    profile_data["prof_loop2"] = val
+                elif "Loop 3" in label:
+                    profile_data["prof_loop3"] = val
+                elif "Loop 4" in label:
+                    profile_data["prof_loop4"] = val
+                elif "Reduce" in label:
+                    profile_data["prof_reduce"] = val
+                elif "Stellar" in label:
+                    profile_data["prof_stellar"] = val
+                elif label.startswith("GR"):
+                    profile_data["prof_gr"] = val
+                elif label.startswith("Transform"):
+                    profile_data["prof_transform"] = val
+                elif label.startswith("Jump"):
+                    profile_data["prof_jump"] = val
+                elif label.startswith("Sync"):
+                    profile_data["prof_sync"] = val
+
+        lines = result.stdout.strip().split('\n')
+        last_line = lines[-1]
+
+        if "PROFILING" in last_line: 
+            pass
+
+        parts = last_line.split(',')
+        walltime = float(parts[0])
+        verification = float(parts[1])
+        energy_err = float(parts[2]) if len(parts) > 2 else None
+        return walltime, verification, energy_err, profile_data
+    except subprocess.CalledProcessError as e:
+        print(f"Error running configuration {config['name']}: {e}")
+        return None, None, None, None
+    except ValueError:
+        print(f"Error parsing output for {config['name']}: {result.stdout}")
+        return None, None, None, None
+
+def write_csv(results, output_file, profile_mode):
+    """Write results to CSV file"""
+    # Use WHFast512 Jacobi (Baseline) as reference for speedup
+    baseline = next((r for r in results if 'Jacobi (Baseline)' in r['name']), results[0])
+    reference_time = baseline['time'] if baseline['time'] is not None else 0
+    reference_x = results[0]['verify'] if results else 0
+
+    with open(output_file, 'w', newline='') as csvfile:
+        fieldnames = ['configuration', 'time_s', 'speedup', 'position_error', 'energy_error',
+                      'profile_mode', 'timestamp']
+
+        # For profiled runs, add per-component timing columns
+        if profile_mode == "profiled":
+            profile_fields = [
+                'prof_total_walltime',
+                'prof_kepler',
+                'prof_kepler_stiefel',
+                'prof_kepler_fg',
+                'prof_interaction',
+                'prof_forces',
+                'prof_loop1',
+                'prof_loop2',
+                'prof_loop3',
+                'prof_loop4',
+                'prof_reduce',
+                'prof_stellar',
+                'prof_gr',
+                'prof_transform',
+                'prof_jump',
+                'prof_sync',
+            ]
+            fieldnames.extend(profile_fields)
+
+        writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
+
+        writer.writeheader()
+
+        timestamp = datetime.now().strftime('%Y-%m-%d %H:%M:%S')
+
+        for res in results:
+            name = res['name']
+            t = res['time']
+
+            if t is None:
+                row = {
+                    'configuration': name,
+                    'time_s': 'FAILED',
+                    'speedup': 'N/A',
+                    'position_error': 'N/A',
+                    'energy_error': 'N/A',
+                    'profile_mode': profile_mode,
+                    'timestamp': timestamp
+                }
+                if profile_mode == "profiled":
+                    for pf in profile_fields:
+                        row[pf] = ''
+                writer.writerow(row)
+                continue
+
+            speedup = reference_time / t if t > 0 else 0
+
+            # Position Error vs reference
+            pos_error = 0.0
+            if reference_x != 0:
+                pos_error = abs(res['verify'] - reference_x) / abs(reference_x)
+
+            # Energy Error (from simulation)
+            energy_error = res.get('energy_err', 0.0)
+            if energy_error is None: energy_error = 0.0
+
+            row = {
+                'configuration': name,
+                'time_s': f"{t:.6f}",
+                'speedup': f"{speedup:.4f}",
+                'position_error': f"{pos_error:.6e}",
+                'energy_error': f"{energy_error:.6e}",
+                'profile_mode': profile_mode,
+                'timestamp': timestamp
+            }
+
+            # Attach profiling numbers if available
+            if profile_mode == "profiled":
+                profile_data = res.get('profile', {}) or {}
+                for pf in profile_fields:
+                    val = profile_data.get(pf, '')
+                    # Format floats nicely, leave missing as empty
+                    if isinstance(val, float):
+                        row[pf] = f"{val:.6f}"
+                    else:
+                        row[pf] = val
+
+            writer.writerow(row)
+
+    print(f"Results written to: {output_file}")
+
+def print_table(results):
+    # Use WHFast512 Jacobi (Baseline) as reference for speedup
+    baseline = next((r for r in results if 'Jacobi (Baseline)' in r['name']), results[0])
+    reference_time = baseline['time'] if baseline['time'] is not None else 0
+    reference_x = results[0]['verify'] if results else 0
+
+    print(f"\n" + "="*100)
+    print(f"{'Method/Configuration':<30} | {'Time (s)':<10} | {'Speedup':<10} | {'Pos Error':<12} | {'Energy Err':<12}")
+    print("-" * 100)
+
+    for res in results:
+        name = res['name']
+        t = res['time']
+
+        if t is None:
+            print(f"{name:<30} | {'FAILED':<10} | {'-':<10} | {'-':<12} | {'-':<12}")
+            continue
+
+        speedup = reference_time / t if t > 0 else 0
+
+        # Position Error vs reference
+        pos_error = 0.0
+        if reference_x != 0:
+            pos_error = abs(res['verify'] - reference_x) / abs(reference_x)
+
+        # Energy Error (from simulation)
+        energy_error = res.get('energy_err', 0.0)
+        if energy_error is None: energy_error = 0.0
+
+        print(f"{name:<30} | {t:<10.4f} | {speedup:<10.2f}x | {pos_error:<12.2e} | {energy_error:<12.2e}")
+    print("=" * 100 + "\n")
+
+def main():
+    parser = argparse.ArgumentParser(description='Run REBOUND Ablation Benchmarks')
+    parser.add_argument('--config', default='benchmarks.json', help='Path to configuration JSON')
+    parser.add_argument('--profile', action='store_true', help='Enable profiling (-DPROF)')
+    parser.add_argument('--filter', help='Run only configs containing this string')
+    parser.add_argument('--output', default=None, help='Custom output CSV filename')
+    args = parser.parse_args()
+
+    configs = load_config(args.config)
+
+    if args.filter:
+        configs = [c for c in configs if args.filter in c['name']]
+
+    if not compile_rebound(args.profile):
+        sys.exit(1)
+
+    results = []
+    print(f"\nRunning {len(configs)} configurations...\n")
+
+    for i, cfg in enumerate(configs):
+        print(f"[{i+1}/{len(configs)}] Running: {cfg['name']}...", end='', flush=True)
+        t, v, e, prof = run_benchmark(cfg)
+        if t is not None:
+            print(f" Done ({t:.4f}s)")
+        else:
+            print(f" FAILED")
+
+        results.append({
+            'name': cfg['name'],
+            'time': t,
+            'verify': v,
+            'energy_err': e,
+            'profile': prof,
+        })
+
+    print_table(results)
+
+    # Determine output filename
+    if args.output:
+        output_file = args.output
+    else:
+        profile_suffix = "_profiled" if args.profile else "_noprofile"
+        timestamp_suffix = datetime.now().strftime('%Y%m%d_%H%M%S')
+        output_file = f"benchmark_results{profile_suffix}_{timestamp_suffix}.csv"
+
+    # Write CSV output
+    profile_mode = "profiled" if args.profile else "no_profile"
+    write_csv(results, output_file, profile_mode)
+
+if __name__ == "__main__":
+    main()

benchmark_results_noprofile_20260219_151028.csv

@@ -0,0 +1,23 @@
+configuration,time_s,speedup,position_error,energy_error,profile_mode,timestamp
+WHFast (scalar),11.648404,0.1619,0.000000e+00,6.573279e-10,no_profile,2026-02-19 15:10:28
+WHFast512 Jacobi (Baseline),1.885897,1.0000,2.900910e-05,6.607087e-10,no_profile,2026-02-19 15:10:28
+WHFast512 Jacobi (Fast RSQRT Only),1.919337,0.9826,2.662890e-05,6.609429e-10,no_profile,2026-02-19 15:10:28
+WHFast512 Jacobi (RegBlock Only),1.868693,1.0092,2.900910e-05,6.607087e-10,no_profile,2026-02-19 15:10:28
+WHFast512 Jacobi (Combined),1.923354,0.9805,2.662890e-05,6.609429e-10,no_profile,2026-02-19 15:10:28
+Cache Stumpff: threshold=1e-10,1.834447,1.0280,nan,nan,no_profile,2026-02-19 15:10:28
+WHFast512 DHC,1.824405,1.0337,2.539597e-01,2.560716e-08,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton for 1/gradient,1.966570,0.9590,2.880580e-05,6.600726e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with prior guess,1.812333,1.0406,2.779027e-05,6.609506e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with momentum (beta=2.0),1.810298,1.0418,2.800145e-05,6.616861e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with momentum (beta=1.5),1.811985,1.0408,2.505478e-05,6.609602e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with momentum (beta=1.7),1.809716,1.0421,2.514072e-05,6.607139e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with momentum (beta=1.8),1.812221,1.0407,2.206318e-05,6.602510e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with momentum (beta=1.9),1.815206,1.0389,2.274855e-05,6.604281e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with momentum (beta=2.1),1.811684,1.0410,2.668819e-05,6.603127e-10,no_profile,2026-02-19 15:10:28
+Kepler Solver: Newton start with momentum (beta=2.2),1.812851,1.0403,2.467140e-05,6.604704e-10,no_profile,2026-02-19 15:10:28
+Gravity: Newton for 1/r^3,1.931823,0.9762,2.662890e-05,6.609429e-10,no_profile,2026-02-19 15:10:28
+Best Kepler (beta=2.1) + Fast RSQRT,1.829108,1.0310,2.961162e-05,6.612183e-10,no_profile,2026-02-19 15:10:28
+Best Kepler (beta=1.9) + Fast RSQRT,1.828055,1.0316,2.407835e-05,6.611019e-10,no_profile,2026-02-19 15:10:28
+Fused Stellar-Jacobi (Baseline Kepler),1.864871,1.0113,2.900910e-05,6.607087e-10,no_profile,2026-02-19 15:10:28
+Fused + Momentum Kepler (beta=2.1),1.750748,1.0772,2.668819e-05,6.603127e-10,no_profile,2026-02-19 15:10:28
+Fused + Momentum + Fast RSQRT,1.742217,1.0825,2.961162e-05,6.612183e-10,no_profile,2026-02-19 15:10:28