bench_pricer.cpp

// VolForge — standalone C++ pricing benchmark (Phase 1.4)
// ---------------------------------------------------------------------------
// Pure-C++ latency (no Python boundary) for:
//   1. SoA SIMD batch Greeks   (price + 7 Greeks over a chain)
//   2. the same via the scalar bs_greeks loop  -> SIMD speedup
//   3. Heston whole-chain Carr-Madan FFT
// Reports per-chain latency, throughput, and the active SIMD backend.
//
//   ./bench_pricer [n_strikes=500] [iters=2000]
// ---------------------------------------------------------------------------
#include <chrono>
#include <cstdio>
#include <vector>

#include "greeks.h"
#include "pricer.h"
#include "simd/normal_cdf_simd.h"

namespace {
template <class F>
double time_us(F&& f, int iters) {
    f();  // warm-up
    const auto t0 = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < iters; ++i) f();
    const auto t1 = std::chrono::high_resolution_clock::now();
    return std::chrono::duration<double, std::micro>(t1 - t0).count() / iters;
}
}  // namespace

int main(int argc, char** argv) {
    using namespace volforge;

    const std::size_t n = (argc > 1) ? std::stoul(argv[1]) : 500;
    const int iters = (argc > 2) ? std::stoi(argv[2]) : 2000;

    const double S = 22000.0, r = 0.065, q = 0.0;  // Nifty50-ish
    std::vector<double> K(n), sigma(n), T(n);
    std::vector<double> price(n), delta(n), gamma(n), vega(n), theta(n), rho(n),
        vanna(n), volga(n);
    for (std::size_t i = 0; i < n; ++i) {
        K[i] = S * (0.80 + 0.40 * static_cast<double>(i) / n);  // 80%-120%
        sigma[i] = 0.12 + 0.05 * (static_cast<double>(i) / n);
        T[i] = 7.0 / 365.0;  // weekly expiry
    }

    auto run_batch = [&] {
        bs_greeks_batch(S, r, q, OptionType::Call, K.data(), sigma.data(),
                        T.data(), n, price.data(), delta.data(), gamma.data(),
                        vega.data(), theta.data(), rho.data(), vanna.data(),
                        volga.data());
    };
    auto run_scalar = [&] {
        BSInputs in;
        in.S = S; in.r = r; in.q = q; in.type = OptionType::Call;
        for (std::size_t i = 0; i < n; ++i) {
            in.K = K[i]; in.sigma = sigma[i]; in.T = T[i];
            const Greeks g = bs_greeks(in);
            price[i] = g.price; delta[i] = g.delta;  // sink
        }
    };

    const double simd_us = time_us(run_batch, iters);
    const double scalar_us = time_us(run_scalar, iters);

    // Heston whole-chain FFT.
    HestonParams hp; hp.v0 = 0.04; hp.kappa = 1.5; hp.theta = 0.04; hp.xi = 0.5; hp.rho = -0.6;
    HestonPricer heston(S, r, q, hp, CarrMadanConfig{});
    const double heston_us = time_us([&] { auto c = heston.price_chain_fft(0.25); (void)c; }, 200);

    std::printf("VolForge bench [backend=%s, lanes=%d]\n", simd::simd_backend(),
                simd::simd_lane_width());
    std::printf("  chain size            : %zu strikes\n", n);
    std::printf("  BS batch (SIMD)       : %.3f us/chain   (%.1f M opt/s)  %s\n",
                simd_us, n / simd_us, simd_us < 1000.0 ? "[<1ms OK]" : "[>1ms]");
    std::printf("  BS batch (scalar loop): %.3f us/chain\n", scalar_us);
    std::printf("  SIMD speedup          : %.2fx\n", scalar_us / simd_us);
    std::printf("  Heston chain FFT      : %.3f us  (%d-pt grid)\n", heston_us,
                1 << CarrMadanConfig{}.n_pow2);
    std::printf("  sample ATM price      : %.4f (delta=%.4f vega=%.4f)\n",
                price[n / 2], delta[n / 2], vega[n / 2]);
    return 0;
}