hypercall_vol_oracle/

vol_surface_cache.rs

//! Volatility surface cache for storing and querying implied volatility data.
//!
//! This module provides an efficient cache for volatility surface data with:
//! - O(1) lookups for exact (strike, expiry) matches
//! - Bilinear interpolation for points between grid nodes
//! - Strike quantization to avoid floating-point comparison issues

use serde::{Deserialize, Serialize};
use std::collections::{BTreeSet, HashMap};
use tracing::debug;

use hypercall_margin::black_scholes::black_scholes;
use hypercall_margin::OptionType;

/// A single point on the volatility surface.
#[derive(Debug, Clone, Serialize, Deserialize, utoipa::ToSchema)]
pub struct VolPoint {
    /// Strike price in USD
    pub strike: f64,
    /// Expiry timestamp in Unix seconds
    pub expiry: i64,
    /// Implied volatility as decimal (0.80 = 80%)
    pub iv: f64,
    /// Timestamp when this point was last updated (milliseconds)
    pub timestamp: i64,
}

/// Exported delta-IV pair for API responses.
#[derive(Debug, Clone, Serialize, Deserialize, utoipa::ToSchema)]
pub struct DeltaIvExport {
    pub delta: f64,
    pub iv: f64,
}

/// Exported delta curve for a single expiry.
#[derive(Debug, Clone, Serialize, Deserialize, utoipa::ToSchema)]
pub struct DeltaCurveExport {
    pub expiry: i64,
    pub points: Vec<DeltaIvExport>,
}

/// A delta-IV pair on the vol smile for a single expiry.
#[derive(Debug, Clone, Copy)]
struct DeltaIvPoint {
    /// Call delta (0.0–1.0)
    delta: f64,
    /// Implied volatility as decimal
    iv: f64,
}

/// Volatility surface for a single underlying.
///
/// Supports two data representations:
/// 1. Strike-space: direct (strike, expiry) → IV grid with bilinear interpolation
/// 2. Delta-space: (delta, expiry) → IV curves, queried via strike→delta conversion
///
/// When both are available, strike-space is preferred (exact). When only delta-space
/// data exists (e.g., from Block Scholes delta.iv feed), the caller provides spot
/// price to convert query strikes to deltas for interpolation.
#[derive(Debug, Clone)]
pub struct VolatilitySurface {
    /// Map from quantized_strike -> expiry -> VolPoint
    points: HashMap<i64, HashMap<i64, VolPoint>>,
    /// ATM volatilities by expiry for quick ATM lookups
    atm_vols: HashMap<i64, f64>,
    /// Delta-IV curves by expiry, sorted by delta ascending
    delta_curves: HashMap<i64, Vec<DeltaIvPoint>>,
    /// Strike quantization factor (e.g., 100 = $100 precision)
    strike_precision: f64,
    /// Last update time for the entire surface (milliseconds)
    last_update: Option<i64>,
    /// All known expiries for iteration
    expiries: BTreeSet<i64>,
}

impl Default for VolatilitySurface {
    fn default() -> Self {
        Self::new()
    }
}

impl VolatilitySurface {
    /// Create a new empty volatility surface with default $1 strike precision.
    pub fn new() -> Self {
        Self::with_precision(1.0)
    }

    /// Create a new volatility surface with custom strike precision.
    ///
    /// # Arguments
    /// * `strike_precision` - Strike price bucket size (e.g., 1.0 for $1 buckets)
    pub fn with_precision(strike_precision: f64) -> Self {
        Self {
            points: HashMap::new(),
            atm_vols: HashMap::new(),
            delta_curves: HashMap::new(),
            strike_precision,
            last_update: None,
            expiries: BTreeSet::new(),
        }
    }

    /// Quantize strike to avoid floating point comparison issues.
    fn quantize_strike(&self, strike: f64) -> i64 {
        (strike / self.strike_precision).round() as i64
    }

    /// Insert or update a volatility point.
    pub fn insert(&mut self, strike: f64, expiry: i64, iv: f64) {
        let quantized_strike = self.quantize_strike(strike);
        let now = chrono::Utc::now().timestamp_millis();

        let point = VolPoint {
            strike,
            expiry,
            iv,
            timestamp: now,
        };

        self.points
            .entry(quantized_strike)
            .or_default()
            .insert(expiry, point);

        self.expiries.insert(expiry);
        self.last_update = Some(now);
    }

    /// Set ATM volatility for an expiry.
    pub fn set_atm_vol(&mut self, expiry: i64, iv: f64) {
        self.atm_vols.insert(expiry, iv);
        self.expiries.insert(expiry);
        self.last_update = Some(chrono::Utc::now().timestamp_millis());
    }

    /// Store a delta-IV point for an expiry.
    ///
    /// Delta-IV curves are used when strike-level data is unavailable.
    /// The caller provides spot price at query time to convert strikes to deltas.
    pub fn set_delta_iv(&mut self, expiry: i64, delta: f64, iv: f64) {
        let curve = self.delta_curves.entry(expiry).or_default();
        // Update existing delta or append
        if let Some(point) = curve.iter_mut().find(|p| (p.delta - delta).abs() < 1e-6) {
            point.iv = iv;
        } else {
            curve.push(DeltaIvPoint { delta, iv });
            curve.sort_by(|a, b| a.delta.partial_cmp(&b.delta).unwrap());
        }
        self.expiries.insert(expiry);
        self.last_update = Some(chrono::Utc::now().timestamp_millis());
    }

    /// Get volatility at exact (strike, expiry) or None.
    pub fn get(&self, strike: f64, expiry: i64) -> Option<f64> {
        let quantized = self.quantize_strike(strike);
        self.points
            .get(&quantized)
            .and_then(|expiries| expiries.get(&expiry))
            .map(|p| p.iv)
    }

    /// Get the full VolPoint at exact (strike, expiry) or None.
    pub fn get_point(&self, strike: f64, expiry: i64) -> Option<&VolPoint> {
        let quantized = self.quantize_strike(strike);
        self.points
            .get(&quantized)
            .and_then(|expiries| expiries.get(&expiry))
    }

    /// Get volatility with interpolation if exact match not found.
    ///
    /// Uses bilinear interpolation across strike and expiry dimensions.
    /// If `spot` is provided and no strike-level data exists, falls back to
    /// delta-space interpolation: converts the query strike to a call delta
    /// and interpolates the delta-IV smile, then interpolates across expiries.
    pub fn get_interpolated_with_spot(
        &self,
        strike: f64,
        expiry: i64,
        spot: Option<f64>,
    ) -> Option<f64> {
        // First try exact match
        if let Some(iv) = self.get(strike, expiry) {
            debug!(
                strike,
                expiry,
                iv,
                spot = ?spot,
                lookup_path = "exact_point",
                "Backend vol surface lookup"
            );
            return Some(iv);
        }

        // Find nearest strikes
        let quantized = self.quantize_strike(strike);
        let strike_keys: Vec<i64> = self.points.keys().copied().collect();

        if strike_keys.is_empty() {
            // No strike-level data — try delta-space interpolation
            if let Some(s) = spot {
                if let Some(iv) = self.interpolate_delta_surface(strike, expiry, s) {
                    debug!(
                        strike,
                        expiry,
                        iv,
                        spot = s,
                        lookup_path = "delta_surface",
                        expiries = self.delta_curves.len(),
                        "Backend vol surface lookup"
                    );
                    return Some(iv);
                }
            }
            // Last resort: flat ATM
            let atm_iv = self.get_atm(expiry);
            if let Some(iv) = atm_iv {
                debug!(
                    strike,
                    expiry,
                    iv,
                    spot = ?spot,
                    lookup_path = "atm_fallback",
                    "Backend vol surface lookup"
                );
            }
            return atm_iv;
        }

        // Find bracketing strikes
        let (lower_strike_q, upper_strike_q) =
            self.find_bracketing_values(quantized, &strike_keys)?;

        // Get expiry maps for each strike
        let lower_expiries = self.points.get(&lower_strike_q)?;
        let upper_expiries = self.points.get(&upper_strike_q)?;

        // Collect all expiries from both strikes
        let all_expiries: BTreeSet<i64> = lower_expiries
            .keys()
            .chain(upper_expiries.keys())
            .copied()
            .collect();

        let expiry_vec: Vec<i64> = all_expiries.iter().copied().collect();
        let (lower_expiry, upper_expiry) = self.find_bracketing_values(expiry, &expiry_vec)?;

        // Get the four corner points
        let ll = lower_expiries.get(&lower_expiry).map(|p| p.iv);
        let lu = lower_expiries.get(&upper_expiry).map(|p| p.iv);
        let ul = upper_expiries.get(&lower_expiry).map(|p| p.iv);
        let uu = upper_expiries.get(&upper_expiry).map(|p| p.iv);

        // Convert quantized strikes back to actual values for interpolation
        let lower_strike = lower_strike_q as f64 * self.strike_precision;
        let upper_strike = upper_strike_q as f64 * self.strike_precision;

        // Bilinear interpolation
        let interpolated = self.bilinear_interpolate(BilinearPoint {
            lower_strike,
            upper_strike,
            lower_expiry,
            upper_expiry,
            strike,
            expiry,
            lower_lower: ll,
            lower_upper: lu,
            upper_lower: ul,
            upper_upper: uu,
        });

        if let Some(iv) = interpolated {
            debug!(
                strike,
                expiry,
                iv,
                spot = ?spot,
                lookup_path = "bilinear_strike_expiry",
                lower_strike,
                upper_strike,
                lower_expiry,
                upper_expiry,
                ll = ?ll,
                lu = ?lu,
                ul = ?ul,
                uu = ?uu,
                "Backend vol surface lookup"
            );
        }

        interpolated
    }

    /// Get volatility with interpolation (no spot price for delta conversion).
    pub fn get_interpolated(&self, strike: f64, expiry: i64) -> Option<f64> {
        self.get_interpolated_with_spot(strike, expiry, None)
    }

    /// Interpolate the delta-IV surface for a given strike and expiry.
    ///
    /// Converts strike to call delta using approximate ATM vol, then
    /// linearly interpolates the delta-IV smile at that delta.
    /// Finally interpolates across expiries.
    fn interpolate_delta_surface(&self, strike: f64, expiry: i64, spot: f64) -> Option<f64> {
        if self.delta_curves.is_empty() {
            return None;
        }

        let expiry_vec: Vec<i64> = self.delta_curves.keys().copied().collect();
        let (lower_exp, upper_exp) = find_bracketing_sorted(expiry, &expiry_vec)?;

        let iv_lower = self.interpolate_delta_smile_at_expiry(strike, lower_exp, spot)?;

        if lower_exp == upper_exp {
            return Some(iv_lower);
        }

        let iv_upper = self.interpolate_delta_smile_at_expiry(strike, upper_exp, spot)?;

        // Linear interpolation across expiries
        let t = (expiry - lower_exp) as f64 / (upper_exp - lower_exp) as f64;
        Some(iv_lower + t * (iv_upper - iv_lower))
    }

    /// Interpolate the delta-IV smile at a single expiry for a given strike.
    fn interpolate_delta_smile_at_expiry(
        &self,
        strike: f64,
        expiry: i64,
        spot: f64,
    ) -> Option<f64> {
        let curve = self.delta_curves.get(&expiry)?;
        if curve.is_empty() {
            return None;
        }

        // Use ATM vol as initial guess for delta computation
        let atm_iv = self.atm_vols.get(&expiry).copied().unwrap_or(0.50);
        let now = chrono::Utc::now().timestamp();
        let tte = (expiry - now) as f64 / (365.25 * 86400.0);
        if tte <= 0.0 {
            // Expired — return ATM vol for this expiry if available
            return self.atm_vols.get(&expiry).copied();
        }

        // Convert strike to call delta: δ = N(d1)
        let sqrt_t = tte.sqrt();
        let d1 = ((spot / strike).ln() + 0.5 * atm_iv * atm_iv * tte) / (atm_iv * sqrt_t);
        let delta = normal_cdf(d1);

        // Interpolate the delta-IV curve
        interpolate_curve(curve, delta)
    }

    /// Get ATM volatility for a given expiry.
    pub fn get_atm(&self, expiry: i64) -> Option<f64> {
        // First check explicit ATM vols
        if let Some(&iv) = self.atm_vols.get(&expiry) {
            return Some(iv);
        }

        // Fallback: try to interpolate ATM vol from nearby expiries
        if self.atm_vols.is_empty() {
            return None;
        }

        let expiry_vec: Vec<i64> = self.atm_vols.keys().copied().collect();
        let (lower, upper) = self.find_bracketing_values(expiry, &expiry_vec)?;

        let lower_iv = self.atm_vols.get(&lower)?;
        let upper_iv = self.atm_vols.get(&upper)?;

        if lower == upper {
            return Some(*lower_iv);
        }

        // Linear interpolation between expiries
        let t = (expiry - lower) as f64 / (upper - lower) as f64;
        Some(lower_iv + t * (upper_iv - lower_iv))
    }

    /// Get the last update timestamp.
    pub fn last_update(&self) -> Option<i64> {
        self.last_update
    }

    /// Get all known expiries.
    pub fn expiries(&self) -> &BTreeSet<i64> {
        &self.expiries
    }

    /// Get the number of points in the surface.
    pub fn len(&self) -> usize {
        self.points.values().map(|m| m.len()).sum::<usize>()
            + self.atm_vols.len()
            + self.delta_curves.values().map(|c| c.len()).sum::<usize>()
    }

    /// Check if the surface is empty.
    pub fn is_empty(&self) -> bool {
        self.points.is_empty() && self.atm_vols.is_empty() && self.delta_curves.is_empty()
    }

    /// Export all strike-space points for serialisation.
    pub fn export_all_points(&self) -> Vec<VolPoint> {
        self.points
            .values()
            .flat_map(|expiry_map| expiry_map.values())
            .cloned()
            .collect()
    }

    /// Export delta-IV curves for every expiry.
    pub fn export_delta_curves(&self) -> Vec<DeltaCurveExport> {
        self.delta_curves
            .iter()
            .map(|(&expiry, curve)| DeltaCurveExport {
                expiry,
                points: curve
                    .iter()
                    .map(|p| DeltaIvExport {
                        delta: p.delta,
                        iv: p.iv,
                    })
                    .collect(),
            })
            .collect()
    }

    /// Export ATM vols as (expiry, iv) pairs.
    pub fn export_atm_vols(&self) -> Vec<(i64, f64)> {
        self.atm_vols
            .iter()
            .map(|(&expiry, &iv)| (expiry, iv))
            .collect()
    }

    /// Enforce static no-arbitrage constraints on the strike-indexed surface.
    ///
    /// For each expiry, walks strikes ascending and clamps IVs down so that:
    /// 1. Call prices are non-increasing in K (no vertical-spread arb)
    /// 2. Call prices are convex in K (no butterfly arb)
    ///
    /// Both constraints follow from risk-neutral pricing; they must hold for any
    /// arb-free surface. We enforce in price space (via Black-Scholes) then back
    /// out the clamped IV by bisection — this handles non-uniform strike grids
    /// correctly and uses the same pricer the rest of the system trusts.
    ///
    /// Only clamps *down*: if an IV is lower than it should be for arb-freeness,
    /// we do nothing. Lowering an IV cannot create new arbs against neighbors we
    /// have already cleaned.
    ///
    /// Returns the number of IV points modified. Pass `risk_free_rate = 0.0` on
    /// testnet.
    pub fn sanitize_arb_free(&mut self, spot: f64, risk_free_rate: f64) -> u32 {
        if spot <= 0.0 {
            return 0;
        }
        let mut total = 0u32;
        let expiries: Vec<i64> = self.expiries.iter().copied().collect();
        let now_ts = chrono::Utc::now().timestamp();

        for expiry in expiries {
            let tte = (expiry - now_ts) as f64 / (365.25 * 86400.0);
            if tte <= 0.0 {
                continue;
            }

            // Gather (strike, iv, quantized_strike) for this expiry, sorted asc.
            let mut slice: Vec<(f64, f64, i64)> = self
                .points
                .iter()
                .filter_map(|(q, exp_map)| exp_map.get(&expiry).map(|p| (p.strike, p.iv, *q)))
                .collect();
            if slice.len() < 2 {
                continue;
            }
            slice.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));

            // Pass 1: monotonicity — C(K_i) <= C(K_{i-1}).
            for i in 1..slice.len() {
                let (k_prev, iv_prev, _) = slice[i - 1];
                let (k_curr, iv_curr, q_curr) = slice[i];
                let c_prev = bs_call(spot, k_prev, tte, risk_free_rate, iv_prev);
                let c_curr = bs_call(spot, k_curr, tte, risk_free_rate, iv_curr);
                if c_curr > c_prev + PRICE_EPS {
                    let new_iv = solve_call_iv(spot, k_curr, tte, risk_free_rate, c_prev, iv_curr);
                    if (new_iv - iv_curr).abs() > IV_EPS {
                        self.update_iv(q_curr, expiry, new_iv);
                        slice[i].1 = new_iv;
                        total += 1;
                    }
                }
            }

            // Pass 2: butterfly / convexity. For convex C(K): at interior K_i,
            // C_i <= w_prev * C_{i-1} + w_next * C_{i+1}, where weights come from
            // linear interpolation between K_{i-1} and K_{i+1}.
            for i in 1..slice.len().saturating_sub(1) {
                let (k_prev, iv_prev, _) = slice[i - 1];
                let (k_curr, iv_curr, q_curr) = slice[i];
                let (k_next, iv_next, _) = slice[i + 1];
                let span = k_next - k_prev;
                if span <= 0.0 {
                    continue;
                }
                let w_prev = (k_next - k_curr) / span;
                let w_next = (k_curr - k_prev) / span;
                let c_prev = bs_call(spot, k_prev, tte, risk_free_rate, iv_prev);
                let c_curr = bs_call(spot, k_curr, tte, risk_free_rate, iv_curr);
                let c_next = bs_call(spot, k_next, tte, risk_free_rate, iv_next);
                let c_max = w_prev * c_prev + w_next * c_next;
                if c_curr > c_max + PRICE_EPS {
                    let new_iv = solve_call_iv(spot, k_curr, tte, risk_free_rate, c_max, iv_curr);
                    if (new_iv - iv_curr).abs() > IV_EPS {
                        self.update_iv(q_curr, expiry, new_iv);
                        slice[i].1 = new_iv;
                        total += 1;
                    }
                }
            }
        }
        total
    }

    /// Mutate the stored IV for an already-present (quantized_strike, expiry).
    fn update_iv(&mut self, quantized_strike: i64, expiry: i64, new_iv: f64) {
        if let Some(exp_map) = self.points.get_mut(&quantized_strike) {
            if let Some(point) = exp_map.get_mut(&expiry) {
                point.iv = new_iv;
                point.timestamp = chrono::Utc::now().timestamp_millis();
            }
        }
    }

    /// Clear all data from the surface.
    pub fn clear(&mut self) {
        self.points.clear();
        self.atm_vols.clear();
        self.delta_curves.clear();
        self.expiries.clear();
        self.last_update = None;
    }

    /// Find bracketing values for interpolation.
    fn find_bracketing_values(&self, target: i64, keys: &[i64]) -> Option<(i64, i64)> {
        if keys.is_empty() {
            return None;
        }

        let mut sorted: Vec<i64> = keys.to_vec();
        sorted.sort();

        match sorted.binary_search(&target) {
            Ok(i) => Some((sorted[i], sorted[i])),
            Err(i) => {
                if i == 0 {
                    // Target is below all values
                    Some((sorted[0], sorted[0]))
                } else if i >= sorted.len() {
                    // Target is above all values
                    let last = *sorted.last()?;
                    Some((last, last))
                } else {
                    Some((sorted[i - 1], sorted[i]))
                }
            }
        }
    }

    /// Perform bilinear interpolation on the four corner values.
    fn bilinear_interpolate(&self, point: BilinearPoint) -> Option<f64> {
        let BilinearPoint {
            lower_strike,
            upper_strike,
            lower_expiry,
            upper_expiry,
            strike,
            expiry,
            lower_lower,
            lower_upper,
            upper_lower,
            upper_upper,
        } = point;

        // Collect available corner values
        let values: Vec<f64> = [lower_lower, lower_upper, upper_lower, upper_upper]
            .iter()
            .filter_map(|v| *v)
            .collect();

        if values.is_empty() {
            return None;
        }

        // If we don't have all four corners, fall back to average
        if values.len() < 4 {
            return Some(values.iter().sum::<f64>() / values.len() as f64);
        }

        let lower_lower = lower_lower?;
        let lower_upper = lower_upper?;
        let upper_lower = upper_lower?;
        let upper_upper = upper_upper?;

        // Handle edge cases where points coincide
        if (upper_strike - lower_strike).abs() < 1e-9 && upper_expiry == lower_expiry {
            return Some(lower_lower);
        }

        // Calculate interpolation weights
        let t = if (upper_strike - lower_strike).abs() < 1e-9 {
            0.0
        } else {
            (strike - lower_strike) / (upper_strike - lower_strike)
        };

        let u = if upper_expiry == lower_expiry {
            0.0
        } else {
            (expiry - lower_expiry) as f64 / (upper_expiry - lower_expiry) as f64
        };

        // Bilinear interpolation formula
        let result = (1.0 - t) * (1.0 - u) * lower_lower
            + t * (1.0 - u) * upper_lower
            + (1.0 - t) * u * lower_upper
            + t * u * upper_upper;

        Some(result)
    }
}

struct BilinearPoint {
    lower_strike: f64,
    upper_strike: f64,
    lower_expiry: i64,
    upper_expiry: i64,
    strike: f64,
    expiry: i64,
    lower_lower: Option<f64>,
    lower_upper: Option<f64>,
    upper_lower: Option<f64>,
    upper_upper: Option<f64>,
}

/// Tolerance for declaring a call-price arb violation (quote-unit dollars).
const PRICE_EPS: f64 = 1e-6;
/// Minimum IV change worth writing back after solving (avoids flapping).
const IV_EPS: f64 = 1e-5;
/// IV floor used during bisection; below this the option is effectively intrinsic.
const IV_FLOOR: f64 = 1e-4;

fn bs_call(spot: f64, strike: f64, tte: f64, r: f64, iv: f64) -> f64 {
    black_scholes(&OptionType::Call, spot, strike, tte, r, iv)
}

/// Find σ in [IV_FLOOR, iv_hint] such that BS call price ≈ target. BS call is
/// monotone increasing in σ, so bisection converges. If the lower bound's price
/// already exceeds target (target below intrinsic), returns IV_FLOOR.
fn solve_call_iv(spot: f64, strike: f64, tte: f64, r: f64, target: f64, iv_hint: f64) -> f64 {
    let mut lo = IV_FLOOR;
    let mut hi = iv_hint.max(IV_FLOOR);
    let c_lo = bs_call(spot, strike, tte, r, lo);
    if target <= c_lo {
        return lo;
    }
    let c_hi = bs_call(spot, strike, tte, r, hi);
    if c_hi <= target {
        return hi; // nothing to lower
    }
    for _ in 0..60 {
        let mid = 0.5 * (lo + hi);
        let c_mid = bs_call(spot, strike, tte, r, mid);
        if c_mid > target {
            hi = mid;
        } else {
            lo = mid;
        }
        if hi - lo < 1e-8 {
            break;
        }
    }
    0.5 * (lo + hi)
}

/// Normal CDF approximation (Abramowitz & Stegun 26.2.17).
fn normal_cdf(x: f64) -> f64 {
    0.5 * (1.0 + libm::erf(x / std::f64::consts::SQRT_2))
}

/// Linear interpolation on a sorted delta-IV curve.
fn interpolate_curve(curve: &[DeltaIvPoint], delta: f64) -> Option<f64> {
    if curve.is_empty() {
        return None;
    }
    if curve.len() == 1 {
        return Some(curve[0].iv);
    }
    // Clamp to curve bounds
    if delta <= curve[0].delta {
        return Some(curve[0].iv);
    }
    if delta >= curve.last().unwrap().delta {
        return Some(curve.last().unwrap().iv);
    }
    // Find bracketing points
    for window in curve.windows(2) {
        let (lo, hi) = (&window[0], &window[1]);
        if delta >= lo.delta && delta <= hi.delta {
            let t = (delta - lo.delta) / (hi.delta - lo.delta);
            return Some(lo.iv + t * (hi.iv - lo.iv));
        }
    }
    None
}

/// Find bracketing values in a sorted slice.
fn find_bracketing_sorted(target: i64, keys: &[i64]) -> Option<(i64, i64)> {
    if keys.is_empty() {
        return None;
    }
    let mut sorted: Vec<i64> = keys.to_vec();
    sorted.sort();
    match sorted.binary_search(&target) {
        Ok(i) => Some((sorted[i], sorted[i])),
        Err(i) => {
            if i == 0 {
                Some((sorted[0], sorted[0]))
            } else if i >= sorted.len() {
                let last = *sorted.last()?;
                Some((last, last))
            } else {
                Some((sorted[i - 1], sorted[i]))
            }
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_insert_and_get() {
        let mut surface = VolatilitySurface::new();
        surface.insert(50000.0, 1735689600, 0.80);

        assert_eq!(surface.get(50000.0, 1735689600), Some(0.80));
        assert_eq!(surface.get(50000.0, 1735689601), None); // Different expiry
        assert_eq!(surface.get(50100.0, 1735689600), None); // Different strike
    }

    #[test]
    fn test_strike_quantization() {
        let mut surface = VolatilitySurface::new();
        surface.insert(50000.0, 1735689600, 0.80);

        // Should find the same point with nearby strike (within ±$0.50 of bucket center)
        // 50000.49 / 1 = 50000.49 rounds to 50000 (same bucket as 50000)
        // 49999.51 / 1 = 49999.51 rounds to 50000 (same bucket as 50000)
        assert_eq!(surface.get(50000.49, 1735689600), Some(0.80));
        assert_eq!(surface.get(49999.51, 1735689600), Some(0.80));

        // Should NOT find with strike in different bucket
        // 50000.51 / 1 = 50000.51 rounds to 50001 (different bucket)
        assert_eq!(surface.get(50000.51, 1735689600), None);
        assert_eq!(surface.get(50001.0, 1735689600), None);
    }

    #[test]
    fn test_atm_vol() {
        let mut surface = VolatilitySurface::new();
        surface.set_atm_vol(1735689600, 0.70);
        surface.set_atm_vol(1735776000, 0.75);

        assert_eq!(surface.get_atm(1735689600), Some(0.70));
        assert_eq!(surface.get_atm(1735776000), Some(0.75));

        // Test interpolation between expiries
        let mid_expiry = (1735689600 + 1735776000) / 2;
        let interpolated = surface.get_atm(mid_expiry);
        assert!(interpolated.is_some());
        let iv = interpolated.unwrap();
        assert!(iv > 0.70 && iv < 0.75, "Expected ~0.725, got {}", iv);
    }

    #[test]
    fn test_interpolation_2x2_grid() {
        let mut surface = VolatilitySurface::new();

        // Create a 2x2 grid
        surface.insert(40000.0, 1735689600, 0.70); // Lower strike, lower expiry
        surface.insert(40000.0, 1735776000, 0.75); // Lower strike, upper expiry
        surface.insert(60000.0, 1735689600, 0.80); // Upper strike, lower expiry
        surface.insert(60000.0, 1735776000, 0.85); // Upper strike, upper expiry

        // Query the center point
        let mid_strike = 50000.0;
        let mid_expiry = (1735689600 + 1735776000) / 2;

        let iv = surface.get_interpolated(mid_strike, mid_expiry);
        assert!(iv.is_some());

        // Expected: average of all four corners = (0.70 + 0.75 + 0.80 + 0.85) / 4 = 0.775
        let expected = 0.775;
        let actual = iv.unwrap();
        assert!(
            (actual - expected).abs() < 0.01,
            "Expected ~{}, got {}",
            expected,
            actual
        );
    }

    #[test]
    fn test_interpolation_exact_match() {
        let mut surface = VolatilitySurface::new();
        surface.insert(50000.0, 1735689600, 0.80);

        // Exact match should return exact value
        let iv = surface.get_interpolated(50000.0, 1735689600);
        assert_eq!(iv, Some(0.80));
    }

    #[test]
    fn test_empty_surface() {
        let surface = VolatilitySurface::new();

        assert!(surface.is_empty());
        assert_eq!(surface.len(), 0);
        assert_eq!(surface.get(50000.0, 1735689600), None);
        assert_eq!(surface.get_interpolated(50000.0, 1735689600), None);
        assert_eq!(surface.get_atm(1735689600), None);
    }

    #[test]
    fn test_clear() {
        let mut surface = VolatilitySurface::new();
        surface.insert(50000.0, 1735689600, 0.80);
        surface.set_atm_vol(1735689600, 0.70);

        assert!(!surface.is_empty());

        surface.clear();

        assert!(surface.is_empty());
        assert_eq!(surface.get(50000.0, 1735689600), None);
        assert_eq!(surface.get_atm(1735689600), None);
    }

    #[test]
    fn test_custom_precision() {
        // Use $1000 strike precision
        let mut surface = VolatilitySurface::with_precision(1000.0);
        surface.insert(50000.0, 1735689600, 0.80);

        // Should find with strikes within ±$500 of bucket center
        // 50499.0 / 1000 = 50.499 rounds to 50 (same bucket as 50000)
        // 49501.0 / 1000 = 49.501 rounds to 50 (same bucket as 50000)
        assert_eq!(surface.get(50499.0, 1735689600), Some(0.80));
        assert_eq!(surface.get(49501.0, 1735689600), Some(0.80));

        // Should NOT find with strike in different bucket
        // 50501.0 / 1000 = 50.501 rounds to 51 (different bucket)
        assert_eq!(surface.get(50501.0, 1735689600), None);
        assert_eq!(surface.get(51500.0, 1735689600), None);
    }

    #[test]
    fn test_expiries_tracking() {
        let mut surface = VolatilitySurface::new();
        surface.insert(50000.0, 1735689600, 0.80);
        surface.insert(50000.0, 1735776000, 0.75);
        surface.set_atm_vol(1735862400, 0.70);

        let expiries = surface.expiries();
        assert_eq!(expiries.len(), 3);
        assert!(expiries.contains(&1735689600));
        assert!(expiries.contains(&1735776000));
        assert!(expiries.contains(&1735862400));
    }

    #[test]
    fn test_partial_grid_interpolation() {
        let mut surface = VolatilitySurface::new();

        // Only have 3 of 4 corners
        surface.insert(40000.0, 1735689600, 0.70);
        surface.insert(40000.0, 1735776000, 0.75);
        surface.insert(60000.0, 1735689600, 0.80);
        // Missing: (60000, 1735776000)

        let mid_strike = 50000.0;
        let mid_expiry = (1735689600 + 1735776000) / 2;

        let iv = surface.get_interpolated(mid_strike, mid_expiry);
        assert!(iv.is_some());

        // Should fall back to average of available points
        let expected = (0.70 + 0.75 + 0.80) / 3.0;
        let actual = iv.unwrap();
        assert!(
            (actual - expected).abs() < 0.01,
            "Expected ~{}, got {}",
            expected,
            actual
        );
    }

    #[test]
    fn test_delta_surface_interpolation() {
        let mut surface = VolatilitySurface::new();
        // Use a far-future expiry so TTE > 0
        let expiry = chrono::Utc::now().timestamp() + 30 * 86400; // 30 days out

        // Set up a vol smile: OTM puts (low delta) have higher IV than ATM
        surface.set_delta_iv(expiry, 0.10, 0.70); // 10-delta: 70%
        surface.set_delta_iv(expiry, 0.25, 0.55); // 25-delta: 55%
        surface.set_delta_iv(expiry, 0.50, 0.45); // ATM: 45%
        surface.set_delta_iv(expiry, 0.75, 0.50); // 75-delta: 50%
        surface.set_delta_iv(expiry, 0.90, 0.60); // 90-delta: 60%
        surface.set_atm_vol(expiry, 0.45);

        let spot = 80000.0;

        // ATM strike (~spot) should give ~ATM IV
        let atm_iv = surface.get_interpolated_with_spot(spot, expiry, Some(spot));
        assert!(atm_iv.is_some());
        let v = atm_iv.unwrap();
        assert!(
            (v - 0.45).abs() < 0.05,
            "ATM query should give ~0.45, got {}",
            v
        );

        // Deep OTM put (very low strike) → low delta → higher IV
        let otm_put_iv = surface.get_interpolated_with_spot(50000.0, expiry, Some(spot));
        assert!(otm_put_iv.is_some());
        assert!(
            otm_put_iv.unwrap() > 0.50,
            "OTM put strike should have higher IV than ATM, got {}",
            otm_put_iv.unwrap()
        );

        // Surface len should include delta points
        assert!(!surface.is_empty());
    }

    #[test]
    fn test_delta_curve_interpolation_helper() {
        let curve = vec![
            DeltaIvPoint {
                delta: 0.10,
                iv: 0.70,
            },
            DeltaIvPoint {
                delta: 0.50,
                iv: 0.45,
            },
            DeltaIvPoint {
                delta: 0.90,
                iv: 0.60,
            },
        ];

        // Exact match
        assert_eq!(interpolate_curve(&curve, 0.50), Some(0.45));

        // Between 0.10 and 0.50 (midpoint delta=0.30)
        let mid = interpolate_curve(&curve, 0.30).unwrap();
        assert!(mid > 0.45 && mid < 0.70, "Expected ~0.575, got {}", mid);

        // Below range → clamp to first
        assert_eq!(interpolate_curve(&curve, 0.01), Some(0.70));

        // Above range → clamp to last
        assert_eq!(interpolate_curve(&curve, 0.99), Some(0.60));
    }

    #[test]
    fn test_normal_cdf() {
        assert!((normal_cdf(0.0) - 0.5).abs() < 1e-6);
        assert!((normal_cdf(1.0) - 0.8413).abs() < 0.01);
        assert!((normal_cdf(-1.0) - 0.1587).abs() < 0.01);
    }

    // --- Export method tests ---

    #[test]
    fn test_export_all_points_empty() {
        let surface = VolatilitySurface::new();
        assert!(surface.export_all_points().is_empty());
    }

    #[test]
    fn test_export_all_points() {
        let mut surface = VolatilitySurface::new();
        surface.insert(50000.0, 1735689600, 0.80);
        surface.insert(60000.0, 1735689600, 0.85);
        surface.insert(50000.0, 1735776000, 0.75);
        assert_eq!(surface.export_all_points().len(), 3);
    }

    #[test]
    fn test_export_delta_curves_empty() {
        let surface = VolatilitySurface::new();
        assert!(surface.export_delta_curves().is_empty());
    }

    #[test]
    fn test_export_delta_curves() {
        let mut surface = VolatilitySurface::new();
        surface.set_delta_iv(1735689600, 0.25, 0.55);
        surface.set_delta_iv(1735689600, 0.50, 0.45);
        surface.set_delta_iv(1735776000, 0.25, 0.60);
        let curves = surface.export_delta_curves();
        assert_eq!(curves.len(), 2);
        for curve in &curves {
            assert!(!curve.points.is_empty());
        }
    }

    #[test]
    fn test_export_atm_vols_empty() {
        let surface = VolatilitySurface::new();
        assert!(surface.export_atm_vols().is_empty());
    }

    #[test]
    fn test_export_atm_vols() {
        let mut surface = VolatilitySurface::new();
        surface.set_atm_vol(1735689600, 0.70);
        surface.set_atm_vol(1735776000, 0.75);
        surface.set_atm_vol(1735862400, 0.80);
        let atm = surface.export_atm_vols();
        assert_eq!(atm.len(), 3);
    }

    // No-arb sanitizer tests live in tests/vol_surface_arb_free_test.rs so they
    // link against the public lib only, bypassing unrelated lib-test build
    // errors in other modules of the crate.
}