wasm4pm 26.6.10

//! Western Electric SPC Rules and Process Capability (Cp, Cpk, DPMO, Six-Sigma).
//!
//! Ported from knhk/rust/knhk-dflss. Pure Rust stdlib math — trivially WASM-compatible.
//!
//! ## Control Chart Monitoring
//!
//! [`check_western_electric_rules`] detects four classes of special causes against a
//! trailing window of [`ChartData`] observations:
//!
//! | Rule | Signal |
//! |------|--------|
//! | Rule 1 | Single point beyond UCL or LCL (beyond 3σ) |
//! | Rule 2 | 9 consecutive points on the same side of the center line |
//! | Rule 3 | 6 consecutive strictly monotone points |
//! | Rule 4 | 2 of 3 consecutive points beyond 2σ on the same side |
//!
//! ## Process Capability
//!
//! [`ProcessCapability::calculate`] computes Cp, Cpk, DPMO, and sigma level from
//! a sample and specification limits. [`dpmo_to_sigma`] applies the standard 1.5σ
//! long-term shift: 3.4 DPMO → 6.0σ.

use tracing::{debug, info};

// ---------------------------------------------------------------------------
// Types (ported from knhk internal/chart.rs)
// ---------------------------------------------------------------------------

/// A single observation on a control chart.
#[derive(Debug, Clone)]
#[allow(dead_code)]
pub struct ChartData {
    /// ISO-8601 timestamp or label.
    pub timestamp: String,
    /// Measured value.
    pub value: f64,
    /// Upper Control Limit.
    pub ucl: f64,
    /// Center Line (process mean).
    pub cl: f64,
    /// Lower Control Limit.
    pub lcl: f64,
    /// Optional subgroup raw values.
    pub subgroup_data: Option<Vec<f64>>,
}

/// Direction of a shift relative to the center line.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[allow(dead_code)]
pub enum ShiftDirection {
    Above,
    Below,
}

/// Direction of a trend.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[allow(dead_code)]
pub enum TrendDirection {
    Increasing,
    Decreasing,
}

/// A detected special-cause signal.
#[derive(Debug, Clone, PartialEq)]
#[allow(dead_code)]
pub enum SpecialCause {
    /// Rule 1: Point beyond UCL or LCL (beyond 3σ).
    OutOfControl { value: f64, ucl: f64, lcl: f64 },
    /// Rule 2: N consecutive points on same side of CL.
    Shift {
        direction: ShiftDirection,
        count: usize,
    },
    /// Rule 3: N consecutive points increasing or decreasing.
    Trend {
        direction: TrendDirection,
        count: usize,
    },
    /// Rule 4: 2 of 3 consecutive points beyond 2σ on the same side of CL.
    TwoOfThree { direction: ShiftDirection },
}

// ---------------------------------------------------------------------------
// Western Electric Rules (ported from knhk rules.rs)
// ---------------------------------------------------------------------------

/// Check the four classic Western Electric special-cause rules against the
/// trailing window of chart data.
///
/// Rules evaluated:
/// 1. Point beyond UCL or LCL (beyond 3σ).
/// 2. 9 consecutive points on same side of center line.
/// 3. 6 consecutive points increasing or decreasing.
/// 4. 2 of 3 consecutive points beyond 2σ on the same side of CL.
///
/// Returns all signals found (may be empty).
///
/// # Example
///
/// ```
/// use wasm4pm::spc::{ChartData, check_western_electric_rules, SpecialCause};
///
/// // Rule 1: a point beyond UCL fires OutOfControl
/// let data = vec![ChartData {
///     timestamp: "t1".to_string(),
///     value: 100.0,
///     ucl: 13.0,
///     cl: 10.0,
///     lcl: 7.0,
///     subgroup_data: None,
/// }];
/// let alerts = check_western_electric_rules(&data);
/// assert_eq!(alerts.len(), 1);
/// assert!(matches!(alerts[0], SpecialCause::OutOfControl { .. }));
///
/// // In-control point → no alerts
/// let ok = vec![ChartData {
///     timestamp: "t2".to_string(),
///     value: 10.5,
///     ucl: 13.0,
///     cl: 10.0,
///     lcl: 7.0,
///     subgroup_data: None,
/// }];
/// assert!(check_western_electric_rules(&ok).is_empty());
/// ```
#[allow(dead_code)]
pub fn check_western_electric_rules(data: &[ChartData]) -> Vec<SpecialCause> {
    let mut alerts = Vec::new();
    let buffer_size = data.len();

    // Emit warm-up event when buffer reaches critical size for rules 2-4
    if buffer_size == 9 {
        debug!(
            buffer_size = buffer_size,
            status = "ok",
            service_name = "wpm",
            "spc.warm_up_complete"
        );
    }

    // Rule 1: Point beyond UCL or LCL (applies to any buffer size, any single point).
    //
    // Defect-class SPC-NaN (this iteration): `NaN > x` and `NaN < x` are both
    // `false`, so a corrupt `NaN` data point used to slip past Rule 1 silently
    // — the chart appeared "in control" even though the input was invalid.
    // Van der Aalst process-mining doctrine treats corrupt evidence as a
    // first-class defect, not a discrepancy. Treat any non-finite value as an
    // out-of-control signal so the autonomic loop notices and degrades safely.
    if let Some(latest) = data.last() {
        if !latest.value.is_finite() || latest.value > latest.ucl || latest.value < latest.lcl {
            // z_score: (value - cl) / sigma, where sigma = (ucl - cl) / 3.0.
            // For non-finite values, z_score is reported as NaN (auditable signal).
            let sigma = (latest.ucl - latest.cl) / 3.0;
            let z_score = if sigma > 0.0 {
                (latest.value - latest.cl) / sigma
            } else {
                f64::NAN
            };
            debug!(
                rule = 1,
                value = latest.value,
                ucl = latest.ucl,
                lcl = latest.lcl,
                alert_level = 3,
                status = "error",
                service_name = "wpm",
                rule_fired = "rule_1",
                "SPC Rule 1 fired: point beyond control limits"
            );
            // GAP-2 IMPLEMENTATION: classified span with canonical rule_violated string,
            // rule_number, and z_score for independent Rank-1 oracle validation.
            // Auditors can verify |z_score| > 3.0 independently from raw values.
            info!(
                target: "autonomic.spc",
                rule_violated = "rule_1_outlier",
                rule_number = 1u32,
                spc_z_score = z_score,
                spc_metric_value = latest.value,
                spc_control_limit_mean = latest.cl,
                spc_control_limit_ucl = latest.ucl,
                spc_control_limit_lcl = latest.lcl,
                spc_consecutive_count = 1u32,
                status = "error",
                service_name = "wpm",
                "spc.rule_violation_classified"
            );
            alerts.push(SpecialCause::OutOfControl {
                value: latest.value,
                ucl: latest.ucl,
                lcl: latest.lcl,
            });
        }
    }

    // Rules 2, 3, and 4 require at least 9 points (Rule 4 needs 3, but we guard uniformly)
    if data.len() < 9 {
        debug!(
            buffer_size,
            "SPC rules 2-4 not yet active, insufficient buffer"
        );
        return alerts;
    }

    let recent = &data[data.len().saturating_sub(9)..];

    // Rule 2: 9 consecutive points on same side of center line.
    // Branchless direction selection via bitmask key.
    if recent.len() >= 9 {
        let above = recent.iter().filter(|d| d.value > d.cl).count() == 9;
        let below = recent.iter().filter(|d| d.value < d.cl).count() == 9;
        const SHIFT_DIR: [Option<ShiftDirection>; 4] = [
            None,
            Some(ShiftDirection::Above),
            Some(ShiftDirection::Below),
            None, // conflict — both true is impossible
        ];
        let key = (above as usize) | ((below as usize) << 1);
        if let Some(dir) = SHIFT_DIR[key] {
            info!(
                rule = 2,
                shift_direction = ?dir,
                count = 9,
                alert_level = 2,
                status = "error",
                service_name = "wpm",
                rule_fired = "rule_2",
                "SPC Rule 2 fired: 9 consecutive points on same side of center line"
            );
            // GAP-2 IMPLEMENTATION: classified span with canonical rule_violated,
            // rule_number, and consecutive_count for Rank-1 oracle (count >= 9).
            let direction_str = match dir {
                ShiftDirection::Above => "above",
                ShiftDirection::Below => "below",
            };
            info!(
                target: "autonomic.spc",
                rule_violated = "rule_2_shift",
                rule_number = 2u32,
                spc_shift_direction = direction_str,
                spc_consecutive_count = 9u32,
                spc_z_score = 0.0f64,  // Rule 2 is count-based; z_score not applicable
                status = "error",
                service_name = "wpm",
                "spc.rule_violation_classified"
            );
            alerts.push(SpecialCause::Shift {
                direction: dir,
                count: 9,
            });
        }
    }

    // Rule 3: 6 consecutive points increasing or decreasing.
    // Branchless direction selection via bitmask key.
    if recent.len() >= 6 {
        let last_6 = &recent[recent.len() - 6..];
        let values: Vec<f64> = last_6.iter().map(|d| d.value).collect();
        let incr = values.windows(2).all(|w| w[1] > w[0]);
        let decr = values.windows(2).all(|w| w[1] < w[0]);
        const TREND_DIR: [Option<TrendDirection>; 4] = [
            None,
            Some(TrendDirection::Increasing),
            Some(TrendDirection::Decreasing),
            None, // conflict — impossible
        ];
        let key = (incr as usize) | ((decr as usize) << 1);
        if let Some(dir) = TREND_DIR[key] {
            info!(
                rule = 3,
                trend_direction = ?dir,
                count = 6,
                alert_level = 2,
                status = "error",
                service_name = "wpm",
                rule_fired = "rule_3",
                "SPC Rule 3 fired: 6 consecutive monotone points"
            );
            // GAP-2 IMPLEMENTATION: classified span with canonical rule_violated,
            // rule_number, and monotonic_sequence_length for Rank-1 oracle (length >= 6).
            let trend_str = match dir {
                TrendDirection::Increasing => "increasing",
                TrendDirection::Decreasing => "decreasing",
            };
            info!(
                target: "autonomic.spc",
                rule_violated = "rule_3_trend",
                rule_number = 3u32,
                spc_trend_direction = trend_str,
                spc_consecutive_count = 6u32,
                spc_monotonic_sequence_length = 6u32,
                spc_z_score = 0.0f64,  // Rule 3 is sequence-based; z_score not applicable
                status = "error",
                service_name = "wpm",
                "spc.rule_violation_classified"
            );
            alerts.push(SpecialCause::Trend {
                direction: dir,
                count: 6,
            });
        }
    }

    // Rule 4: 2 of 3 consecutive points beyond 2σ on the same side of CL.
    //
    // The 2σ boundary is derived from the control limits stored in ChartData:
    //   sigma = (ucl - cl) / 3.0
    //   upper_2sigma = cl + 2.0 * sigma  (i.e. cl + (ucl - cl) * 2/3)
    //   lower_2sigma = cl - 2.0 * sigma  (i.e. cl - (ucl - cl) * 2/3)
    //
    // We scan the last 3 points. If ≥2 are beyond the 2σ line above CL, or ≥2 are
    // beyond the 2σ line below CL, the rule fires.
    if recent.len() >= 3 {
        let last_3 = &recent[recent.len() - 3..];

        // Count how many of the last 3 points are beyond 2σ above or below CL.
        let beyond_2sigma_above = last_3
            .iter()
            .filter(|d| {
                let sigma = (d.ucl - d.cl) / 3.0;
                sigma > 0.0 && d.value > d.cl + 2.0 * sigma
            })
            .count();

        let beyond_2sigma_below = last_3
            .iter()
            .filter(|d| {
                let sigma = (d.ucl - d.cl) / 3.0;
                sigma > 0.0 && d.value < d.cl - 2.0 * sigma
            })
            .count();

        if beyond_2sigma_above >= 2 {
            info!(
                rule = 4,
                direction = ?ShiftDirection::Above,
                count = beyond_2sigma_above,
                alert_level = 2,
                status = "error",
                service_name = "wpm",
                rule_fired = "rule_4",
                "SPC Rule 4 fired: 2+ of 3 points beyond 2σ above center line"
            );
            // GAP-2 IMPLEMENTATION: classified span with canonical rule_violated and consecutive_count.
            info!(
                target: "autonomic.spc",
                rule_violated = "rule_4_two_of_three",
                rule_number = 4u32,
                spc_shift_direction = "above",
                spc_consecutive_count = beyond_2sigma_above as u32,
                spc_z_score = 0.0f64,  // Rule 4 uses 2σ threshold, not 3σ; z_score is min 2.0
                status = "error",
                service_name = "wpm",
                "spc.rule_violation_classified"
            );
            alerts.push(SpecialCause::TwoOfThree {
                direction: ShiftDirection::Above,
            });
        } else if beyond_2sigma_below >= 2 {
            info!(
                rule = 4,
                direction = ?ShiftDirection::Below,
                count = beyond_2sigma_below,
                alert_level = 2,
                status = "error",
                service_name = "wpm",
                rule_fired = "rule_4",
                "SPC Rule 4 fired: 2+ of 3 points beyond 2σ below center line"
            );
            // GAP-2 IMPLEMENTATION: classified span with canonical rule_violated and consecutive_count.
            info!(
                target: "autonomic.spc",
                rule_violated = "rule_4_two_of_three",
                rule_number = 4u32,
                spc_shift_direction = "below",
                spc_consecutive_count = beyond_2sigma_below as u32,
                spc_z_score = 0.0f64,  // Rule 4 uses 2σ threshold, not 3σ; z_score is min 2.0
                status = "error",
                service_name = "wpm",
                "spc.rule_violation_classified"
            );
            alerts.push(SpecialCause::TwoOfThree {
                direction: ShiftDirection::Below,
            });
        }
    }

    // Log summary of alerts detected
    if !alerts.is_empty() {
        debug!(
            alert_count = alerts.len(),
            rule_fired = ?alerts
                .iter()
                .map(|a| match a {
                    SpecialCause::OutOfControl { .. } => "rule_1",
                    SpecialCause::Shift { .. } => "rule_2",
                    SpecialCause::Trend { .. } => "rule_3",
                    SpecialCause::TwoOfThree { .. } => "rule_4",
                })
                .collect::<Vec<_>>(),
            status = "error",
            service_name = "wpm",
            "SPC check completed with alerts"
        );
    } else {
        debug!(
            alert_count = 0,
            status = "ok",
            service_name = "wpm",
            "SPC check completed: no alerts"
        );
    }

    alerts
}

// ---------------------------------------------------------------------------
// Process Capability (ported from knhk capability.rs)
// ---------------------------------------------------------------------------

/// Computed process capability indices.
#[derive(Debug, Clone, PartialEq)]
#[allow(dead_code)]
pub struct ProcessCapability {
    pub cp: f64,
    pub cpk: f64,
    pub sigma_level: f64,
    pub dpmo: f64,
    pub mean: f64,
    pub std_dev: f64,
    pub usl: f64,
    pub lsl: f64,
}

/// Errors returned by capability calculations.
#[derive(Debug, Clone, PartialEq)]
#[allow(dead_code)]
pub enum CapabilityError {
    EmptyData,
    InvalidLimits,
}

impl std::fmt::Display for CapabilityError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            CapabilityError::EmptyData => write!(f, "Cannot calculate capability with empty data"),
            CapabilityError::InvalidLimits => {
                write!(
                    f,
                    "Invalid specification limits: USL must be greater than LSL"
                )
            }
        }
    }
}

impl std::error::Error for CapabilityError {}

impl ProcessCapability {
    /// Calculate Cp, Cpk, sigma level, and DPMO from observed data and
    /// specification limits.
    #[allow(dead_code)]
    pub fn calculate(data: &[f64], usl: f64, lsl: f64) -> Result<Self, CapabilityError> {
        if data.is_empty() {
            debug!(
                status = "error",
                service_name = "wpm",
                reason = "empty_data",
                "process_capability.calculate failed"
            );
            return Err(CapabilityError::EmptyData);
        }

        if usl <= lsl {
            debug!(
                status = "error",
                service_name = "wpm",
                reason = "invalid_limits",
                "process_capability.calculate failed"
            );
            return Err(CapabilityError::InvalidLimits);
        }

        let mean = spc_mean(data);
        let std_dev = spc_std_dev(data);

        if std_dev == 0.0 {
            // All data points are identical.
            let is_within_limits = data.iter().all(|&x| x >= lsl && x <= usl);
            debug!(
                mean = mean,
                cp = if is_within_limits { f64::INFINITY } else { 0.0 },
                cpk = if is_within_limits { f64::INFINITY } else { 0.0 },
                sigma_level = if is_within_limits { 6.0 } else { 0.0 },
                status = if is_within_limits { "ok" } else { "error" },
                service_name = "wpm",
                "process_capability.calculate zero variance"
            );
            return Ok(Self {
                cp: if is_within_limits { f64::INFINITY } else { 0.0 },
                cpk: if is_within_limits { f64::INFINITY } else { 0.0 },
                sigma_level: if is_within_limits { 6.0 } else { 0.0 },
                dpmo: if is_within_limits { 0.0 } else { 1_000_000.0 },
                mean,
                std_dev,
                usl,
                lsl,
            });
        }

        // Cp: Process Potential Capability
        let cp = (usl - lsl) / (6.0 * std_dev);

        // Cpk: Process Performance Capability
        let cpk_usl = (usl - mean) / (3.0 * std_dev);
        let cpk_lsl = (mean - lsl) / (3.0 * std_dev);
        let cpk = cpk_usl.min(cpk_lsl);

        // DPMO and Sigma Level
        let z_usl = (usl - mean) / std_dev;
        let z_lsl = (lsl - mean) / std_dev;

        let p_usl = 1.0 - normal_cdf(z_usl);
        let p_lsl = normal_cdf(z_lsl);
        let p_defective = p_usl + p_lsl;

        let dpmo = p_defective * 1_000_000.0;
        let sigma_level = dpmo_to_sigma(dpmo);

        debug!(
            cp = cp,
            cpk = cpk,
            sigma_level = sigma_level,
            dpmo = dpmo,
            mean = mean,
            std_dev = std_dev,
            status = "ok",
            service_name = "wpm",
            "process_capability.calculate completed"
        );

        Ok(Self {
            cp,
            cpk,
            sigma_level,
            dpmo,
            mean,
            std_dev,
            usl,
            lsl,
        })
    }
}

// ---------------------------------------------------------------------------
// Statistics helpers (ported from knhk statistics.rs)
// ---------------------------------------------------------------------------

/// Returns the arithmetic mean of `data`, or `0.0` for an empty slice.
///
/// # Examples
///
/// ```
/// use wasm4pm::spc::spc_mean;
///
/// assert_eq!(spc_mean(&[1.0, 2.0, 3.0]), 2.0);
/// assert_eq!(spc_mean(&[]), 0.0);
/// ```
#[allow(dead_code)]
pub fn spc_mean(data: &[f64]) -> f64 {
    if data.is_empty() {
        return 0.0;
    }
    data.iter().sum::<f64>() / data.len() as f64
}

/// Returns the sample standard deviation (Bessel's correction, N−1 denominator).
/// Returns `0.0` for slices shorter than 2 elements.
///
/// # Examples
///
/// ```
/// use wasm4pm::spc::spc_std_dev;
///
/// // [1, 2, 3]: mean=2, sum-sq-dev=2, N-1=2 → variance=1.0 → σ=1.0
/// assert_eq!(spc_std_dev(&[1.0, 2.0, 3.0]), 1.0);
/// assert_eq!(spc_std_dev(&[]), 0.0);
/// assert_eq!(spc_std_dev(&[42.0]), 0.0);
/// ```
#[allow(dead_code)]
pub fn spc_std_dev(data: &[f64]) -> f64 {
    if data.len() < 2 {
        return 0.0;
    }
    let m = spc_mean(data);
    let variance = data.iter().map(|&x| (x - m).powi(2)).sum::<f64>() / (data.len() - 1) as f64;
    variance.sqrt()
}

// ---------------------------------------------------------------------------
// Normal CDF / inverse CDF (hand-written, no statrs dependency)
// ---------------------------------------------------------------------------

/// Abramowitz & Stegun approximation of the standard normal CDF Φ(z).
///
/// # Examples
///
/// ```
/// use wasm4pm::spc::normal_cdf;
///
/// // Φ(0) = 0.5 by symmetry of the standard normal
/// assert!((normal_cdf(0.0) - 0.5).abs() < 1e-6);
/// assert!(normal_cdf(1.0) > 0.5);
/// assert!(normal_cdf(-1.0) < 0.5);
/// ```
#[allow(dead_code)]
pub fn normal_cdf(z: f64) -> f64 {
    let t = 1.0 / (1.0 + 0.2316419 * z.abs());
    let d = 0.39894228 * (-z * z / 2.0).exp();
    let prob = 1.0
        - d * t
            * (0.319381530
                + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))));
    if z > 0.0 {
        prob
    } else {
        1.0 - prob
    }
}

/// Rational approximation of the standard normal inverse CDF (Peter Acklam's
/// algorithm).  Handles all three regions: low tail, central, high tail.
/// Maximum absolute error ≈ 1.15 × 10⁻⁹.
#[allow(dead_code)]
pub fn inverse_normal_cdf(p: f64) -> f64 {
    if p <= 0.0 {
        return f64::NEG_INFINITY;
    }
    if p >= 1.0 {
        return f64::INFINITY;
    }

    // Acklam's rational approximation coefficients.
    const A: [f64; 6] = [
        -3.969_683_028_665_376e+01,
        2.209_460_984_245_205e+02,
        -2.759_285_104_469_687e+02,
        1.383_577_518_672_69e2,
        -3.066_479_806_614_716e+01,
        2.506_628_277_459_239e+00,
    ];
    const B: [f64; 5] = [
        -5.447_609_879_822_406e+01,
        1.615_858_368_580_409e+02,
        -1.556_989_798_598_866e+02,
        6.680_131_188_771_972e+01,
        -1.328_068_155_288_572e+01,
    ];
    const C: [f64; 6] = [
        -7.784_894_002_430_293e-03,
        -3.223_964_580_411_365e-01,
        -2.400_758_277_161_838e+00,
        -2.549_732_539_343_734e+00,
        4.374_664_141_464_968e+00,
        2.938_163_982_698_783e+00,
    ];
    const D: [f64; 4] = [
        7.784_695_709_041_462e-03,
        3.224_671_290_700_398e-01,
        2.445_134_137_142_996e+00,
        3.754_408_661_907_416e+00,
    ];

    const P_LOW: f64 = 0.02425;
    const P_HIGH: f64 = 1.0 - P_LOW;

    if p < P_LOW {
        // Lower region.
        let q = (-2.0 * p.ln()).sqrt();
        (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
    } else if p <= P_HIGH {
        // Central region.
        let q = p - 0.5;
        let r = q * q;
        q * (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5])
            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
    } else {
        // Upper region.
        let q = (-2.0 * (1.0 - p).ln()).sqrt();
        -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
    }
}

/// Public wrapper for normal CDF (used by benchmarks and JTBD validation).
#[allow(dead_code)]
pub fn normal_cdf_public(z: f64) -> f64 {
    normal_cdf(z)
}

/// Public wrapper for inverse normal CDF (used by benchmarks and JTBD validation).
#[allow(dead_code)]
pub fn inverse_normal_cdf_public(p: f64) -> f64 {
    inverse_normal_cdf(p)
}

/// Converts DPMO to a Six-Sigma level (includes the standard 1.5σ long-term shift).
///
/// # Examples
///
/// ```
/// use wasm4pm::spc::dpmo_to_sigma;
///
/// // Six-Sigma quality: 3.4 DPMO → ~6.0σ
/// assert!((dpmo_to_sigma(3.4) - 6.0).abs() < 0.1);
/// // Perfect quality → capped at 6.0σ
/// assert_eq!(dpmo_to_sigma(0.0), 6.0);
/// // Worst case → 0σ
/// assert_eq!(dpmo_to_sigma(1_000_000.0), 0.0);
/// ```
#[allow(dead_code)]
pub fn dpmo_to_sigma(dpmo: f64) -> f64 {
    if dpmo <= 0.0 {
        return 6.0;
    }
    if dpmo >= 1_000_000.0 {
        return 0.0;
    }

    let p_defective = dpmo / 1_000_000.0;
    let z_score = inverse_normal_cdf(1.0 - p_defective);
    z_score + 1.5 // Add 1.5 sigma shift for short-term vs long-term
}

// Tests consolidated in tests/autonomic_tests.rs (spc_tests module)