"""Does the official QuantLib wheel print the same numbers on every Linux system?

Prices a few options with QuantLib (pip install QuantLib==1.43) and prints every result with 17
significant digits and as a hex float, so comparing two outputs compares bits. It also hashes the
system maths library's exp, log and pow over fixed inputs. All inputs are built with + and * only,
so they are the same bits on every machine.

  python ql_last_digits.py > a.txt               run once in each environment
  python ql_last_digits.py --compare a.txt b.txt count the values that differ

Optional: --dump FILE also writes the raw exp/log/pow results; --compare A.bin B.bin then counts the
inputs on which the two maths libraries disagree.
"""
import array
import hashlib
import math
import os
import platform
import sys

N = 1_000_000


def hexf(v):
    return float(v).hex()


def libm_inputs():
    """Inputs typical of pricing code, built without transcendental functions."""
    xe = [(i - N // 2) * 1.2345678e-5 for i in range(N)]             # exp: about -6.2 to 6.2
    xl = [0.5 + i * 1.5e-6 for i in range(N)]                        # log: 0.5 to 2
    xp = [(1.0 + i * 1e-7, 0.25 + (i % 97) * 0.125) for i in range(N)]  # pow: base 1 to 1.1
    return xe, xl, xp


def libm(dump=None):
    xe, xl, xp = libm_inputs()
    out = {
        "exp": array.array("d", map(math.exp, xe)),
        "log": array.array("d", map(math.log, xl)),
        "pow": array.array("d", (math.pow(b, e) for b, e in xp)),
    }
    for name, a in out.items():
        print("libm %s n=%d sha256=%s" % (name, len(a), hashlib.sha256(a.tobytes()).hexdigest()[:16]))
    if dump:
        with open(dump, "wb") as f:
            for name in ("exp", "log", "pow"):
                out[name].tofile(f)


def show(name, v):
    print("%s = %.17g %s" % (name, v, hexf(v)))


def quantlib():
    import QuantLib as ql

    today = ql.Date(15, ql.May, 2025)
    ql.Settings.instance().evaluationDate = today
    dc, cal = ql.Actual365Fixed(), ql.TARGET()
    spot = ql.QuoteHandle(ql.SimpleQuote(100.0))
    r = ql.YieldTermStructureHandle(ql.FlatForward(today, 0.03, dc))
    q = ql.YieldTermStructureHandle(ql.FlatForward(today, 0.01, dc))
    vol = ql.BlackVolTermStructureHandle(ql.BlackConstantVol(today, cal, 0.25, dc))
    process = ql.BlackScholesMertonProcess(spot, q, r, vol)
    expiry = today + ql.Period(1, ql.Years)

    def option(kind, strike, exercise=None):
        return ql.VanillaOption(ql.PlainVanillaPayoff(kind, strike), exercise or ql.EuropeanExercise(expiry))

    show("discount factor, 1.5 years", r.discount(1.5))
    show("process drift at t=0.5", process.drift(0.5, 100.0))

    o = option(ql.Option.Call, 105.0)
    o.setPricingEngine(ql.AnalyticEuropeanEngine(process))
    show("European call, Black-Scholes formula", o.NPV())
    show("European call, delta", o.delta())

    o.setPricingEngine(ql.MCEuropeanEngine(process, "pseudorandom", timeSteps=52, requiredSamples=20000, seed=42))
    show("European call, Monte Carlo", o.NPV())

    a = ql.AmericanExercise(today, expiry)
    o = option(ql.Option.Put, 95.0, a)
    o.setPricingEngine(ql.BinomialCRRVanillaEngine(process, 801))
    show("American put, binomial tree", o.NPV())
    o.setPricingEngine(ql.MCAmericanEngine(process, "pseudorandom", timeSteps=50, antitheticVariate=True,
                                           requiredSamples=8000, seed=17, nCalibrationSamples=2048))
    show("American put, Longstaff-Schwartz", o.NPV())

    fixings = [cal.advance(today, i, ql.Months) for i in range(1, 13)]
    o = ql.DiscreteAveragingAsianOption(ql.Average.Arithmetic, 0.0, 0, fixings,
                                        ql.PlainVanillaPayoff(ql.Option.Call, 100.0), ql.EuropeanExercise(fixings[-1]))
    o.setPricingEngine(ql.MCDiscreteArithmeticAPEngine(process, "pseudorandom", requiredSamples=20000, seed=11))
    show("Asian call, Monte Carlo", o.NPV())

    o = ql.BarrierOption(ql.Barrier.DownOut, 85.0, 0.0, ql.PlainVanillaPayoff(ql.Option.Call, 100.0),
                         ql.EuropeanExercise(expiry))
    o.setPricingEngine(ql.MCBarrierEngine(process, "pseudorandom", timeSteps=252, requiredSamples=10000, seed=19))
    show("Down-and-out call, Monte Carlo", o.NPV())

    heston = ql.HestonProcess(r, q, spot, 0.0625, 1.5, 0.0625, 0.4, -0.6)
    o = option(ql.Option.Call, 100.0)
    o.setPricingEngine(ql.AnalyticHestonEngine(ql.HestonModel(heston)))
    show("European call, Heston formula", o.NPV())

    # Local volatility from a Black volatility surface with a smile (polynomial in moneyness, no log).
    pillars = [cal.advance(today, p) for p in (ql.Period(1, ql.Months), ql.Period(3, ql.Months),
                                               ql.Period(6, ql.Months), ql.Period(1, ql.Years),
                                               ql.Period(18, ql.Months), ql.Period(2, ql.Years))]
    strikes = [60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0]
    atm = [0.31, 0.29, 0.27, 0.255, 0.245, 0.24]
    vols = ql.Matrix(len(strikes), len(pillars))
    for i, k in enumerate(strikes):
        x = k / 100.0 - 1.0
        for j in range(len(pillars)):
            vols[i][j] = atm[j] - 0.12 * x + 0.30 * x * x / (1.0 + 0.5 * j)
    surface = ql.BlackVolTermStructureHandle(ql.BlackVarianceSurface(today, cal, pillars, strikes, vols, dc))
    local = ql.BlackScholesMertonProcess(spot, q, r, surface)
    o = option(ql.Option.Put, 102.0, ql.EuropeanExercise(cal.advance(today, 6, ql.Months)))
    o.setPricingEngine(ql.MCEuropeanEngine(local, "pseudorandom", timeSteps=52, requiredSamples=8000, seed=42))
    show("European put, local volatility, Monte Carlo", o.NPV())
    o.setPricingEngine(ql.MCEuropeanEngine(local, "lowdiscrepancy", timeSteps=16, brownianBridge=True,
                                           requiredSamples=1023, seed=7))
    show("European put, local volatility, Sobol", o.NPV())


def environment():
    libc = os.confstr("CS_GNU_LIBC_VERSION") if hasattr(os, "confstr") else "?"
    try:
        import QuantLib as ql
        qlv = ql.__version__
    except ImportError:
        qlv = "not installed"
    print("# Python %s | %s | QuantLib %s | %s" % (platform.python_version(), libc, qlv, platform.machine()))


def compare_text(fa, fb):
    def values(f):
        res = {}
        for line in open(f):
            line = line.rstrip("\n")
            if line.startswith("libm "):
                key, _, h = line.partition(" sha256=")
                res[key] = h
            elif " = " in line and not line.startswith("#"):
                key, _, rest = line.partition(" = ")
                res[key] = rest.split()[-1]
        return res

    a, b = values(fa), values(fb)
    diff = [k for k in a if k in b and a[k] != b[k]]
    worst = 0.0
    for k in diff:
        if k.startswith("libm"):
            print("differs: %s (hash)" % k)
            continue
        x, y = float.fromhex(a[k]), float.fromhex(b[k])
        rel = abs(x - y) / max(abs(x), abs(y))
        worst = max(worst, rel)
        print("differs: %-40s %.17g vs %.17g (relative %.1e)" % (k, x, y, rel))
    priced = [k for k in a if k in b and not k.startswith("libm")]
    print("%d of %d printed values differ; largest relative difference %.1e"
          % (sum(1 for k in diff if not k.startswith("libm")), len(priced), worst))


def compare_dumps(fa, fb):
    a, b = array.array("d"), array.array("d")
    with open(fa, "rb") as f:
        a.frombytes(f.read())
    with open(fb, "rb") as f:
        b.frombytes(f.read())
    for i, name in enumerate(("exp", "log", "pow")):
        x, y = a[i * N:(i + 1) * N], b[i * N:(i + 1) * N]
        n = sum(1 for u, v in zip(x, y) if u != v)
        ulps = max((abs(u - v) / math.ulp(max(abs(u), abs(v))) for u, v in zip(x, y) if u != v), default=0)
        print("libm %s: %d of %d results differ (%.3f%%), by at most %d ulp" % (name, n, N, 100.0 * n / N, ulps))


if __name__ == "__main__":
    if len(sys.argv) == 4 and sys.argv[1] == "--compare":
        (compare_dumps if sys.argv[2].endswith(".bin") else compare_text)(sys.argv[2], sys.argv[3])
    else:
        dump = sys.argv[2] if len(sys.argv) == 3 and sys.argv[1] == "--dump" else None
        environment()
        libm(dump)
        quantlib()