src/quantization.py

"""
================================================================================
HETMAN-2-NCM 2.0B LOW-LEVEL COMPUTE FABRIC: QUANTIZATION & FWHT INT8 ENGINE
================================================================================
Subsystem: BitNet b1.58 Ternary Core & Fast Walsh-Hadamard Transform (FWHT B=512)
Specification Reference: spec/02_BITNET_TERNARY_AND_FWHT_INT8_ENGINE.md
                         spec/06_PHASE2_BIDIRECTIONAL_NCM_REFINER_AND_DEDUP.md
                         spec/08_TRAINING_RECIPE_LOSSES_AND_CURRICULUM.md

Features:
1. Fast Walsh-Hadamard Transform (FWHT B=512):
   - Vectorized 9-stage butterfly network in O(d log B) add/sub operations.
   - Outlier attenuation: peak amplitude reduced by sqrt(512) ≈ 22.63x.
   - Hadamard Gain Clamping: dynamic scale eta clamped to [eta_min, eta_max]
     to eliminate recurrent amplification across passes.
2. BitNet b1.58 Ternary Quantization:
   - Absmean scaling with Straight-Through Estimator (STE).
   - Weights mapped to discrete ternary alphabet {-1, 0, +1}.
   - Packed 2-bit physical storage (4 weights per byte, 64-byte cache line alignment).
   - Branchless hardware bitmask decoding (Format 0).
3. Dynamic Symmetric INT8 Activation Quantization:
   - Scale = 127.0 / eta with zero clipping loss under FWHT rotation.
4. Production BitLinear Module:
   - Fused Sub-Layer RMSNorm, FWHT B=512, INT8 activation quantization,
     ternary weight contraction, and INT32 accumulation.
5. Outlier Clamping Safeguards:
   - SwiGLUClamped: activation bounded to [-10.0, +10.0] before down-projection.
   - LogitSoftCapping: classifier logits bounded via C_cap * tanh(logits / C_cap) <= 30.0.
6. Bit-Level Hardware GEMV Simulator:
   - C-level integer add/sub accumulation emulator validating bit-exact parity.
================================================================================
"""

import math
import struct
from typing import Optional, Tuple, Union

import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F


# ==============================================================================
# 0. NUMERICAL CONSTANTS & HARDWARE PARAMETERS
# ==============================================================================

FWHT_BLOCK_SIZE_DEFAULT: int = 512
HADAMARD_NORMALIZATION_512: float = 1.0 / math.sqrt(512.0)  # ~0.04419417382415922
DEFAULT_ETA_MIN: float = 0.05
DEFAULT_ETA_MAX: float = 40.0
DEFAULT_SWIGLU_CLAMP: float = 10.0
DEFAULT_LOGIT_SOFT_CAP: float = 30.0
RADA_CACHE_LINE_BYTES: int = 64
_SHIFTS_4: torch.Tensor = torch.tensor([0, 2, 4, 6], dtype=torch.uint8)



# ==============================================================================
# 1. VECTORIZED FAST WALSH-HADAMARD TRANSFORM (FWHT B=512)
# ==============================================================================

def fast_walsh_hadamard_transform_512(x: torch.Tensor, normalize: bool = True) -> torch.Tensor:
    """
    Vectorized Fast Walsh-Hadamard Transform (FWHT) for blocks of length B = 512.
    Decomposes H_512 into 9 recursive butterfly stages in FP32 with zero multiplications.

    Mathematical formulation:
        y = (1 / sqrt(512)) * (H_2^{\\otimes 9}) * x

    Properties:
        - Symmetric: H^T = H
        - Orthogonal: H^T * H = I
        - Self-inverse (Involution): H * H = I (when normalized)
        - Condition number: kappa(H) = 1.0

    Args:
        x: Input tensor of shape (..., 512).
        normalize: If True, applies 1/sqrt(512) normalization factor.

    Returns:
        Transformed tensor of identical shape and original dtype.
    """
    assert x.shape[-1] == 512, f"Input last dimension must be 512, got {x.shape[-1]}"
    orig_dtype = x.dtype
    # Perform arithmetic in float32 for absolute numerical stability across butterfly stages
    out = x.to(torch.float32)
    lead_shape = out.shape[:-1]
    out = out.reshape(-1, 512)

    # 9-Stage Butterfly Network: strides 1, 2, 4, 8, 16, 32, 64, 128, 256
    h = 1
    while h < 512:
        out = out.reshape(-1, 512 // (2 * h), 2, h)
        a = out[:, :, 0, :]
        b = out[:, :, 1, :]
        out = torch.stack([a + b, a - b], dim=2).reshape(-1, 512)
        h *= 2

    out = out.reshape(*lead_shape, 512)
    if normalize:
        out = out * HADAMARD_NORMALIZATION_512

    return out.to(orig_dtype)


def blockwise_fwht_forward(
    x: torch.Tensor,
    block_size: int = 512,
    normalize: bool = True
) -> torch.Tensor:
    """
    Applies FWHT block-by-block across arbitrary hidden dimensions divisible by block_size.
    For Hetman-2-NCM d_model = 1536: executes 3 independent 512-point transforms.
    For intermediate_dim = 6144: executes 12 independent 512-point transforms.

    Args:
        x: Input tensor of shape (..., d) where d % block_size == 0.
        block_size: FWHT block size (default 512).
        normalize: Whether to apply isometric 1/sqrt(block_size) scaling.

    Returns:
        Orthogonally rotated tensor of shape (..., d).
    """
    orig_shape = x.shape
    d = orig_shape[-1]
    if d % block_size != 0:
        raise ValueError(f"Dimension {d} is not divisible by FWHT block size {block_size}")

    x_reshaped = x.reshape(-1, d // block_size, block_size)
    x_transformed = fast_walsh_hadamard_transform_512(x_reshaped, normalize=normalize)
    return x_transformed.reshape(orig_shape)


class FastWalshHadamardTransform(nn.Module):
    """
    Online Blockwise Fast Walsh-Hadamard Transform Module with Hadamard Gain Clamping.

    Disperses activation energy across 512-dim blocks, reducing peak outlier magnitude
    by sqrt(512) ≈ 22.63x and enabling stable INT8 quantization without clipping artifacts.

    Hadamard Gain Clamping:
        Dynamic scale eta is bounded by [eta_min, eta_max] to prevent recurrent
        noise amplification or vanishing representation across multi-pass inference.
    """
    def __init__(
        self,
        block_size: int = FWHT_BLOCK_SIZE_DEFAULT,
        eta_min: float = DEFAULT_ETA_MIN,
        eta_max: float = DEFAULT_ETA_MAX,
        simulate_int8: bool = True
    ):
        super().__init__()
        self.block_size = block_size
        self.eta_min = eta_min
        self.eta_max = eta_max
        self.simulate_int8 = simulate_int8

    def forward(self, x: torch.Tensor) -> Tuple[torch.Tensor, float]:
        """
        Applies blockwise FWHT, dynamic INT8 quantization with gain clamping,
        and inverse FWHT (since FWHT is self-inverse).

        Returns:
            Tuple of (reconstructed_tensor, mean_clamped_eta)
        """
        orig_shape = x.shape
        d = orig_shape[-1]
        assert d % self.block_size == 0, f"Last dimension {d} must be a multiple of {self.block_size}"

        # 0. Quarantine / Self-Healing against NaN / Inf contamination
        if torch.isnan(x).any() or torch.isinf(x).any():
            x = torch.nan_to_num(x, nan=0.0, posinf=self.eta_max, neginf=-self.eta_max)

        # 1. Forward FWHT rotation
        x_rot = blockwise_fwht_forward(x, block_size=self.block_size, normalize=True)

        if not self.simulate_int8:
            return x_rot, 0.0

        # 2. Dynamic Outlier scale estimation with Gain Clamping
        # eta = max(|x_rot|) + 1e-5
        eta = torch.amax(torch.abs(x_rot), dim=-1, keepdim=True) + 1e-5
        eta_clamped = torch.clamp(eta, self.eta_min, self.eta_max)

        # 3. Dynamic INT8 quantization & dequantization (STE emulation)
        scale = 127.0 / eta_clamped
        x_scaled = x_rot * scale
        # Clamped strictly to [-127.0, 127.0] to prevent the -128 two's complement negation bug
        x_quant = torch.clamp(torch.round(x_scaled), -127.0, 127.0)
        x_quant_ste = x_scaled + (x_quant - x_scaled).detach()
        x_dequant = x_quant_ste * (eta_clamped / 127.0)

        # 4. Inverse FWHT (FWHT is self-inverse due to orthogonal symmetry)
        x_inv = blockwise_fwht_forward(x_dequant, block_size=self.block_size, normalize=True)
        return x_inv.reshape(orig_shape), eta_clamped.mean().item()


# ==============================================================================
# 2. STRAIGHT-THROUGH ESTIMATORS (STE) FOR TERNARY WEIGHTS & INT8 ACTIVATIONS
# ==============================================================================

class AbsmeanWeightQuantizer(torch.autograd.Function):
    """
    Straight-Through Estimator (STE) for Absmean BitNet b1.58 Ternary Quantization {-1, 0, +1}.

    Forward Pass:
        gamma = mean(|W|).clamp(min=eps)  [Computed in FP32]
        W_tilde = Clip(Round(W / gamma), -1.0, 1.0)

    Backward Pass:
        Straight-Through Estimator with soft boundary saturation threshold tau,
        compensated by (1.0 / gamma) to cancel scale fusion scaling and eliminate
        the 100,000x gradient vanishing trap near zero:
        dL/dW = (dL/dW_tilde) * (1.0 / gamma) * (1 if |W / gamma| <= tau else 0.1)
    """
    @staticmethod
    def forward(
        ctx,
        weight: torch.Tensor,
        eps: float = 1e-6,
        tau: float = 1.0
    ) -> Tuple[torch.Tensor, torch.Tensor]:
        orig_dtype = weight.dtype
        w_f32 = weight.to(torch.float32)
        # Quarantine corrupted weights before absmean computation
        if torch.isnan(w_f32).any() or torch.isinf(w_f32).any():
            w_f32 = torch.nan_to_num(w_f32, nan=0.0, posinf=1.0, neginf=-1.0)
        # Step 1: Compute per-tensor mean absolute scale gamma in FP32 for stability
        gamma = torch.mean(torch.abs(w_f32)).clamp(min=eps)
        # Step 2: Normalize and round to {-1, 0, +1}
        w_scaled = w_f32 / gamma
        w_ternary = torch.clamp(torch.round(w_scaled), -1.0, 1.0).to(orig_dtype)

        ctx.save_for_backward(w_scaled)
        ctx.tau = tau
        ctx.gamma = gamma
        return w_ternary, gamma.to(orig_dtype)

    @staticmethod
    def backward(ctx, grad_w_ternary: torch.Tensor, grad_gamma: Optional[torch.Tensor] = None):
        w_scaled, = ctx.saved_tensors
        tau = ctx.tau
        gamma = ctx.gamma

        # STE: Pass gradient through uninhibited inside [-tau, tau], dampening outside
        mask = (torch.abs(w_scaled) <= tau).to(grad_w_ternary.dtype)
        # Soft boundary gradient prevents saturated weights from permanently deadlocking
        ste_damp = torch.where(mask > 0, torch.ones_like(mask), torch.full_like(mask, 0.1))
        # Compensate for gamma factor from forward dequantization (dL/dW = dL/dW_tilde * (1/gamma))
        grad_weight = (grad_w_ternary / gamma) * ste_damp
        return grad_weight, None, None


class SymmetricActivationQuantizer(torch.autograd.Function):
    """
    Dynamic Symmetric Per-Tensor INT8 Quantizer with Straight-Through Estimator.

    Forward Pass:
        eta = max(|X|)  [Computed in FP32]
        eta_clamped = clamp(eta, eta_min, eta_max)  [Hadamard Gain Clamping]
        scale = 127.0 / max(eta_clamped, eps)
        X_q = Clip(Round(X * scale), -127.0, 127.0)  [Safe symmetric bounds: [-127, +127]]

    Backward Pass:
        dL/dX = dL/dX_q  (Straight-Through)
    """
    @staticmethod
    def forward(
        ctx,
        x: torch.Tensor,
        eps: float = 1e-6,
        eta_min: Optional[float] = None,
        eta_max: Optional[float] = None
    ) -> Tuple[torch.Tensor, torch.Tensor]:
        orig_dtype = x.dtype
        x_f32 = x.to(torch.float32)

        # Compute dynamic per-tensor maximum magnitude in FP32
        eta = torch.amax(torch.abs(x_f32), dim=-1, keepdim=True).clamp(min=eps)
        if eta_min is not None and eta_max is not None:
            eta = torch.clamp(eta, eta_min, eta_max)

        scale = 127.0 / eta
        # STRICT BOUND: [-127.0, 127.0] avoids -128 SIMD two's complement negation bug
        x_q = torch.clamp(torch.round(x_f32 * scale), -127.0, 127.0).to(orig_dtype)
        ctx.save_for_backward(x)
        return x_q, scale.to(orig_dtype)

    @staticmethod
    def backward(ctx, grad_x_q: torch.Tensor, grad_scale: Optional[torch.Tensor] = None):
        # Straight-Through Estimator
        return grad_x_q.clone(), None, None, None


# ==============================================================================
# 3. 2-BIT PACKED STORAGE & BRANCHLESS DECODING ENGINE
# ==============================================================================

class PackedTernaryStorage:
    """
    Hardware-aligned 2-Bit Packed Weight Storage Container.

    Bit-level layouts (4 weights per byte):
        Format 0 (Arithmetic Optimized - Default):
            00_2 =  0
            01_2 = +1
            11_2 = -1 (0x03)
            10_2 =  0 (reserved / NOP)
        Format 1 (Rada v2 Compatibility):
            00_2 =  0
            01_2 = +1
            10_2 = -1 (0x02)
            11_2 =  0 (reserved / NOP)

    Unpacking bitmasks (vectorized via broadcast shifts):
        Bitmask_P = B & 0x55       (01010101_2)
        Bitmask_N = (B >> 1) & 0x55
        Format 0: V_k = P - 2 * (P & N)
        Format 1: V_k = P - N
    """

    @staticmethod
    def is_64byte_aligned(arr: Union[np.ndarray, torch.Tensor]) -> bool:
        """Verifies if buffer address satisfies 64-byte hardware cache line boundary."""
        if isinstance(arr, np.ndarray):
            return (arr.ctypes.data % RADA_CACHE_LINE_BYTES) == 0
        elif isinstance(arr, torch.Tensor):
            return (arr.data_ptr() % RADA_CACHE_LINE_BYTES) == 0
        return False

    @staticmethod
    def allocate_aligned_uint8(shape: Tuple[int, ...]) -> np.ndarray:
        """Allocates a 64-byte aligned numpy uint8 array for cache-line compliance."""
        total_bytes = int(np.prod(shape))
        raw = np.empty(total_bytes + RADA_CACHE_LINE_BYTES, dtype=np.uint8)
        offset = (RADA_CACHE_LINE_BYTES - (raw.ctypes.data % RADA_CACHE_LINE_BYTES)) % RADA_CACHE_LINE_BYTES
        aligned = raw[offset:offset + total_bytes].reshape(shape)
        assert (aligned.ctypes.data % RADA_CACHE_LINE_BYTES) == 0
        return aligned

    @classmethod
    def pack(cls, ternary_weights: torch.Tensor, format_mode: int = 0) -> np.ndarray:
        """
        Packs a 2D ternary weight tensor [N, K] into 64-byte aligned uint8 array [N, K // 4].
        Requires K to be divisible by 4.
        """
        assert ternary_weights.dim() == 2, f"Expected 2D tensor, got {ternary_weights.dim()}D"
        N, K = ternary_weights.shape
        assert K % 4 == 0, f"In_features {K} must be divisible by 4"

        packed_torch = cls.pack_torch(ternary_weights, format_mode=format_mode)
        aligned = cls.allocate_aligned_uint8(packed_torch.shape)
        aligned[:] = packed_torch.cpu().numpy()
        return aligned

    @staticmethod
    def pack_torch(ternary_weights: torch.Tensor, format_mode: int = 0) -> torch.Tensor:
        """
        Pure PyTorch implementation of 2-bit packing (executes on CUDA/MPS/CPU directly).
        Vectorized with zero unnecessary memory allocations:
        Format 0 uses single bitwise AND (w & 3) mapping {-1->3, 0->0, 1->1} in 1 instruction.
        Format 1 maps {-1->2, 0->0, 1->1} for Rada v2 compatibility.
        """
        assert ternary_weights.dim() == 2, f"Expected 2D tensor, got {ternary_weights.dim()}D"
        N, K = ternary_weights.shape
        assert K % 4 == 0, f"In_features {K} must be divisible by 4"

        w = ternary_weights.detach().to(torch.int8)

        if format_mode == 0:
            # Format 0: 0->0, 1->1, -1 (0xFF)->3 (0x03) in single bitwise AND
            encoded = (w & 3).to(torch.uint8).view(N, K // 4, 4)
        else:
            # Format 1 (Rada v2): 0->0, 1->1, -1->2
            encoded = torch.where(w == -1, torch.tensor(2, dtype=torch.uint8, device=w.device),
                      torch.where(w == 1, torch.tensor(1, dtype=torch.uint8, device=w.device),
                                  torch.tensor(0, dtype=torch.uint8, device=w.device))).view(N, K // 4, 4)

        packed = (
            encoded[..., 0] |
            (encoded[..., 1] << 2) |
            (encoded[..., 2] << 4) |
            (encoded[..., 3] << 6)
        )
        return packed.contiguous()

    @staticmethod
    def unpack(packed_bytes: np.ndarray, out_features: int, in_features: int, format_mode: int = 0) -> torch.Tensor:
        """
        Unpacks packed uint8 array [N, K // 4] back into INT8 PyTorch tensor [N, K]
        using branchless bitmask arithmetic.
        """
        packed_t = torch.from_numpy(packed_bytes)
        return PackedTernaryStorage.unpack_torch(packed_t, out_features, in_features, format_mode=format_mode)

    @staticmethod
    def unpack_torch(packed: torch.Tensor, out_features: int, in_features: int, format_mode: int = 0) -> torch.Tensor:
        """
        Pure PyTorch GPU/CPU branchless vectorized unpacking of packed 2-bit weights.
        Uses parallel broadcast shifts without non-contiguous slice writes (2.7x speedup).
        """
        device = packed.device
        shifts = _SHIFTS_4.to(device)

        chunks = (packed.unsqueeze(-1) >> shifts) & 3
        p = (chunks & 1).to(torch.int8)
        n = ((chunks >> 1) & 1).to(torch.int8)

        if format_mode == 0:
            # Format 0: 0->0, 1->1, 3->-1, 2->0
            decoded = p - 2 * (p & n)
        else:
            # Format 1 (Rada v2): 0->0, 1->1, 2->-1, 3->0
            decoded = p - n

        return decoded.view(out_features, in_features).contiguous()


# ==============================================================================
# 4. RMSNORM COMPATIBILITY LAYER
# ==============================================================================

class SubLayerRMSNorm(nn.Module):
    """
    Sub-Layer RMSNorm with learnable scaling weight and epsilon.
    Compatible with all PyTorch releases.
    """
    def __init__(self, dim: int, eps: float = 1e-6):
        super().__init__()
        self.dim = dim
        self.eps = eps
        self.weight = nn.Parameter(torch.ones(dim, dtype=torch.float32))

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        orig_dtype = x.dtype
        x_f32 = x.to(torch.float32)
        variance = x_f32.pow(2).mean(dim=-1, keepdim=True)
        rsqrt = torch.rsqrt(variance + self.eps)
        normalized = x_f32 * rsqrt * self.weight
        return normalized.to(orig_dtype)


# ==============================================================================
# 5. PRODUCTION BITLINEAR MODULE WITH FUSED FWHT & GAIN CLAMPING
# ==============================================================================

class BitLinear(nn.Module):
    """
    Production BitLinear Layer for Hetman-2-NCM.

    Combines:
        1. Sub-Layer RMSNorm (In-register normalization).
        2. Blockwise Fast Walsh-Hadamard Transform (FWHT B=512) for outlier suppression.
        3. Dynamic Symmetric INT8 Activation Quantization with Hadamard Gain Clamping.
        4. Absmean b1.58 Ternary Weight Quantization with Straight-Through Estimator (STE).
        5. Fused INT32 Integer Contraction & Floating-Point Dequantization.
    """
    def __init__(
        self,
        in_features: int,
        out_features: int,
        bias: bool = False,
        use_fwht: bool = True,
        fwht_block_size: int = FWHT_BLOCK_SIZE_DEFAULT,
        eta_min: float = DEFAULT_ETA_MIN,
        eta_max: float = DEFAULT_ETA_MAX,
        rms_norm_eps: float = 1e-6,
        use_rms_norm: bool = True
    ):
        super().__init__()
        self.in_features = in_features
        self.out_features = out_features
        self.use_fwht = use_fwht and (in_features % fwht_block_size == 0)
        self.fwht_block_size = fwht_block_size
        self.eta_min = eta_min
        self.eta_max = eta_max
        self.use_rms_norm = use_rms_norm

        # Latent master weights (maintained in FP32 / BF16 during optimization)
        self.weight = nn.Parameter(torch.empty(out_features, in_features, dtype=torch.float32))
        if bias:
            self.bias = nn.Parameter(torch.zeros(out_features, dtype=torch.float32))
        else:
            self.register_parameter('bias', None)

        if self.use_rms_norm:
            self.rms_norm = SubLayerRMSNorm(in_features, eps=rms_norm_eps)
        else:
            self.rms_norm = None

        self.reset_parameters()

    def reset_parameters(self):
        # Variance-scaled initialization for ternary weight convergence
        std = 1.0 / math.sqrt(self.in_features)
        nn.init.uniform_(self.weight, -math.sqrt(3) * std, math.sqrt(3) * std)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # Step 0: Input numerical quench against NaN / Inf corruption
        if torch.isnan(x).any() or torch.isinf(x).any():
            x = torch.nan_to_num(x, nan=0.0, posinf=self.eta_max, neginf=-self.eta_max)

        # Step 1: Sub-Layer RMSNorm
        if self.rms_norm is not None:
            x_norm = self.rms_norm(x)
        else:
            x_norm = x

        # Step 2: Online Blockwise FWHT (B=512) for outlier attenuation
        if self.use_fwht:
            x_rot = blockwise_fwht_forward(x_norm, block_size=self.fwht_block_size, normalize=True)
        else:
            x_rot = x_norm

        # Step 3: Dynamic Symmetric INT8 Activation Quantization with Gain Clamping
        x_q, scale_x = SymmetricActivationQuantizer.apply(
            x_rot, 1e-6, self.eta_min, self.eta_max
        )

        # Step 4: Absmean Ternary Weight Quantization {-1, 0, +1}
        w_ternary, gamma_w = AbsmeanWeightQuantizer.apply(self.weight)

        # Step 5: Integer Matrix Contraction (simulating integer SIMD add/sub)
        # y_int32 = x_q (INT8) @ w_ternary^T ({-1, 0, +1})
        y_int32 = F.linear(x_q, w_ternary, bias=None)

        # Step 6: Dequantization and Output Scale Fusion (computed safely in FP32)
        scale_fusion = (gamma_w.to(torch.float32) / scale_x.to(torch.float32)).to(x.dtype)
        y_out = y_int32 * scale_fusion

        if self.bias is not None:
            y_out = y_out + self.bias

        return y_out

    def export_packed_weights(
        self,
        format_mode: int = 0,
        return_torch: bool = False
    ) -> Tuple[Union[np.ndarray, torch.Tensor], float]:
        """
        Quantizes latent weights and packs them into 2-bit format.
        Args:
            format_mode: 0 for Format 0 ({-1:3, 0:0, 1:1}), 1 for Format 1 / Rada ({-1:2, 0:0, 1:1}).
            return_torch: If True, returns torch.Tensor (avoids CPU/NumPy copy); else np.ndarray.
        Returns:
            packed_weights: uint8 array/tensor of shape [out_features, in_features // 4]
            gamma_w: float scalar absmean scale factor
        """
        with torch.no_grad():
            w_f32 = self.weight.to(torch.float32)
            gamma_w = torch.mean(torch.abs(w_f32)).clamp(min=1e-6).item()
            w_scaled = w_f32 / gamma_w
            w_ternary = torch.clamp(torch.round(w_scaled), -1.0, 1.0)
            if return_torch:
                packed = PackedTernaryStorage.pack_torch(w_ternary, format_mode=format_mode)
            else:
                packed = PackedTernaryStorage.pack(w_ternary, format_mode=format_mode)
        return packed, gamma_w

    def load_packed_weights(
        self,
        packed_weights: Union[np.ndarray, torch.Tensor],
        gamma_w: float,
        format_mode: int = 0
    ):
        """
        Restores latent weights from packed 2-bit storage and scale factor.
        Supports either NumPy array or PyTorch Tensor directly.
        """
        with torch.no_grad():
            if isinstance(packed_weights, torch.Tensor):
                w_int8 = PackedTernaryStorage.unpack_torch(
                    packed_weights, self.out_features, self.in_features, format_mode=format_mode
                )
            else:
                w_int8 = PackedTernaryStorage.unpack(
                    packed_weights, self.out_features, self.in_features, format_mode=format_mode
                )
            self.weight.copy_(w_int8.to(self.weight.dtype) * gamma_w)


# ==============================================================================
# 6. OUTLIER CLAMPING MODULES (SwiGLU & LOGIT SOFT-CAPPING)
# ==============================================================================

class SwiGLUClamped(nn.Module):
    """
    SwiGLU Feed-Forward Network with Activation Clamping <= 10.0.

    Mathematical Formulation:
        gate = SiLU(W_gate(x))
        up   = W_up(x)
        activated = gate * up
        clamped = Clamp(activated, -clamp_val, +clamp_val)  [Default clamp_val = 10.0]
        output = W_down(clamped)

    Mitigates activation explosion and protects INT8 quantization dynamic range
    during deep multi-layer reasoning.
    """
    def __init__(
        self,
        d_model: int,
        d_ffn: int,
        clamp_val: float = DEFAULT_SWIGLU_CLAMP,
        use_bitlinear: bool = True,
        use_fwht: bool = True
    ):
        super().__init__()
        self.d_model = d_model
        self.d_ffn = d_ffn
        self.clamp_val = clamp_val

        if use_bitlinear:
            self.w_gate = BitLinear(d_model, d_ffn, bias=False, use_fwht=use_fwht)
            self.w_up = BitLinear(d_model, d_ffn, bias=False, use_fwht=use_fwht)
            self.w_down = BitLinear(d_ffn, d_model, bias=False, use_fwht=use_fwht)
        else:
            self.w_gate = nn.Linear(d_model, d_ffn, bias=False)
            self.w_up = nn.Linear(d_model, d_ffn, bias=False)
            self.w_down = nn.Linear(d_ffn, d_model, bias=False)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        gate = F.silu(self.w_gate(x))
        up = self.w_up(x)
        activated = gate * up

        # Outlier clamping: strict bound [-10.0, +10.0]
        clamped = torch.clamp(activated, -self.clamp_val, self.clamp_val)
        return self.w_down(clamped)


def swiglu_clamped_functional(
    x: torch.Tensor,
    w_gate: nn.Module,
    w_up: nn.Module,
    w_down: nn.Module,
    clamp_val: float = DEFAULT_SWIGLU_CLAMP
) -> torch.Tensor:
    """Functional implementation of SwiGLU activation clamping."""
    gate = F.silu(w_gate(x))
    up = w_up(x)
    activated = gate * up
    clamped = torch.clamp(activated, -clamp_val, clamp_val)
    return w_down(clamped)


class LogitSoftCapping(nn.Module):
    """
    Logit Soft-Capping Operator (<= 30.0).

    Mathematical Formulation:
        ell_capped = C_cap * tanh(ell_raw / C_cap)

    Derivative:
        d ell_capped / d ell_raw = 1 - tanh^2(ell_raw / C_cap)

    Properties:
        - When |ell_raw| -> 0: ell_capped ≈ ell_raw (Identity).
        - When |ell_raw| -> infinity: |ell_capped| is strictly bounded by C_cap (30.0).
        - Gradients gracefully saturate to 0 for extreme outliers, preventing NaN loss spikes.
    """
    def __init__(self, cap_val: float = DEFAULT_LOGIT_SOFT_CAP):
        super().__init__()
        self.cap_val = cap_val

    def forward(self, logits: torch.Tensor) -> torch.Tensor:
        return self.cap_val * torch.tanh(logits / self.cap_val)


def logit_soft_capping(logits: torch.Tensor, cap_val: float = DEFAULT_LOGIT_SOFT_CAP) -> torch.Tensor:
    """Functional implementation of logit soft-capping."""
    return cap_val * torch.tanh(logits / cap_val)


# ==============================================================================
# 7. HARDWARE GEMV SIMULATOR (BRANCHLESS INTEGER ADD/SUB EMULATION)
# ==============================================================================

def simulate_hardware_bitnet_gemv(
    x_int8: np.ndarray,              # [K] signed int8
    packed_w_uint8: np.ndarray,      # [N, K // 4] uint8
    scale_x: float,
    gamma_w: float,
    K: int,
    N: int
) -> np.ndarray:
    """
    Simulates the exact AVX-512 VNNI / AVX2 / ARM NEON low-level GEMV kernel.

    Features:
        - In-register branchless bitmask decoding (Format 0).
        - Integer addition and subtraction accumulation (Zero FP multiplications).
        - INT32 master accumulator (safe against INT16 overflow).
        - Single broadcast FP32 dequantization scale fusion.

    Args:
        x_int8: Signed 8-bit integer activation array of length K.
        packed_w_uint8: 2-bit packed weight array of shape [N, K // 4].
        scale_x: Dynamic activation scale factor (127.0 / eta).
        gamma_w: Absmean weight scale factor.
        K: Input dimension (must match length of x_int8).
        N: Output dimension (rows in packed_w_uint8).

    Returns:
        Dequantized float32 output array of length N.
    """
    assert x_int8.dtype == np.int8, f"Expected int8, got {x_int8.dtype}"
    assert packed_w_uint8.dtype == np.uint8, f"Expected uint8, got {packed_w_uint8.dtype}"
    assert x_int8.shape[0] == K, f"Expected K={K}, got {x_int8.shape[0]}"
    assert packed_w_uint8.shape == (N, K // 4), f"Expected shape ({N}, {K // 4})"

    y_out = np.zeros(N, dtype=np.float32)
    scale_fusion = gamma_w / scale_x

    # Row-by-row integer SIMD accumulation
    for row in range(N):
        acc_int32 = 0
        w_row = packed_w_uint8[row]

        for col_byte in range(K // 4):
            b = int(w_row[col_byte])
            x_base = col_byte * 4

            # Branchless bitmask decoding of 4 ternary weights
            p = b & 0x55
            n = (b >> 1) & 0x55

            p0 = p & 1
            n0 = n & 1
            w0 = p0 - 2 * (p0 & n0)

            p1 = (p >> 2) & 1
            n1 = (n >> 2) & 1
            w1 = p1 - 2 * (p1 & n1)

            p2 = (p >> 4) & 1
            n2 = (n >> 4) & 1
            w2 = p2 - 2 * (p2 & n2)

            p3 = (p >> 6) & 1
            n3 = (n >> 6) & 1
            w3 = p3 - 2 * (p3 & n3)

            # Integer add/sub operations only
            acc_int32 += int(x_int8[x_base + 0]) * w0
            acc_int32 += int(x_int8[x_base + 1]) * w1
            acc_int32 += int(x_int8[x_base + 2]) * w2
            acc_int32 += int(x_int8[x_base + 3]) * w3

        # FP32 scale fusion dequantization
        y_out[row] = float(acc_int32) * scale_fusion

    return y_out


# ==============================================================================
# 8. AUTONOMOUS VERIFICATION & COMPLIANCE TEST SUITE
# ==============================================================================

def run_quantization_engine_tests() -> bool:
    """
    Autonomous comprehensive verification suite validating 100% functionality:
    1. FWHT B=512 Orthogonality and Isometry (||H*x|| == ||x||).
    2. FWHT Involution (H * H = I).
    3. FWHT Outlier Attenuation (peak reduced by sqrt(512) ≈ 22.63x).
    4. Hadamard Gain Clamping (bounds eta within [eta_min, eta_max]).
    5. 2-Bit Weight Packing & Unpacking Identity (both NumPy & Torch, 4:1 compression).
    6. Absmean STE Quantizer & Gradient Backpropagation.
    7. Symmetric INT8 Activation Quantizer & STE Flow.
    8. BitLinear End-to-End Execution (Forward, Backward, Optimizer Step).
    9. BitLinear vs. Hardware GEMV Simulator Exact Numerical Parity.
    10. SwiGLU Activation Clamping (bounded <= 10.0).
    11. Logit Soft-Capping (strictly bounded <= 30.0).
    12. 64-Byte Cache-Line Memory Alignment Compliance.
    """
    print("=" * 80)
    print("STARTING HETMAN-2-NCM QUANTIZATION & FWHT ENGINE TEST HARNESS")
    print("=" * 80)

    torch.manual_seed(42)
    np.random.seed(42)

    passed_tests = 0
    total_tests = 15

    # --------------------------------------------------------------------------
    # Test 1: FWHT Orthogonality and Isometry
    # --------------------------------------------------------------------------
    print("\n[Test 1/15] FWHT B=512 Orthogonality and Isometry...")
    x = torch.randn(8, 512, dtype=torch.float32)
    norm_before = torch.norm(x, dim=-1)
    x_rot = fast_walsh_hadamard_transform_512(x, normalize=True)
    norm_after = torch.norm(x_rot, dim=-1)
    norm_err = torch.max(torch.abs(norm_before - norm_after)).item()
    print(f"  Max Euclidean norm discrepancy: {norm_err:.8e}")
    assert norm_err < 1e-5, f"Isometry violated: norm discrepancy {norm_err}"
    print("  -> PASS: Euclidean norm strictly conserved.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 2: FWHT Involution Identity (H * H = I)
    # --------------------------------------------------------------------------
    print("\n[Test 2/15] FWHT B=512 Involution Identity (H * H = I)...")
    x_inv = fast_walsh_hadamard_transform_512(x_rot, normalize=True)
    inv_err = torch.max(torch.abs(x - x_inv)).item()
    print(f"  Involution reconstruction error: {inv_err:.8e}")
    assert inv_err < 1e-5, f"Involution violated: {inv_err}"
    print("  -> PASS: H is self-inverse with zero reconstruction distortion.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 3: Catastrophic Outlier Suppression (22.63x)
    # --------------------------------------------------------------------------
    print("\n[Test 3/15] Outlier Attenuation via FWHT B=512...")
    spike = torch.zeros(1, 512, dtype=torch.float32)
    spike[0, 256] = 100.0  # Massive 100.0 activation spike
    spike_rot = fast_walsh_hadamard_transform_512(spike, normalize=True)
    max_peak = torch.max(torch.abs(spike_rot)).item()
    expected_peak = 100.0 / math.sqrt(512.0)
    print(f"  Input Spike Peak   : 100.0")
    print(f"  Rotated Max Peak   : {max_peak:.4f} (Expected: {expected_peak:.4f})")
    print(f"  Attenuation Factor : {100.0 / max_peak:.2f}x")
    assert abs(max_peak - expected_peak) < 1e-4, f"Outlier attenuation mismatch: {max_peak}"
    assert max_peak <= 4.5, "Outlier magnitude was not suppressed below 4.5"
    print("  -> PASS: Activation spike attenuated by exactly 22.63x.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 4: Hadamard Gain Clamping
    # --------------------------------------------------------------------------
    print("\n[Test 4/15] Hadamard Gain Clamping Stability...")
    fwht_mod = FastWalshHadamardTransform(block_size=512, eta_min=0.05, eta_max=40.0)

    # Subtest A: Giant input that would push eta > 40
    giant_input = torch.randn(2, 512) * 500.0
    _, clamped_eta_giant = fwht_mod(giant_input)
    print(f"  Giant input clamped eta: {clamped_eta_giant:.4f} (Should be <= 40.0)")
    assert clamped_eta_giant <= 40.0 + 1e-5, f"Gain clamping failed for upper bound: {clamped_eta_giant}"

    # Subtest B: Near-zero input that would drop eta < 0.05
    tiny_input = torch.randn(2, 512) * 1e-5
    _, clamped_eta_tiny = fwht_mod(tiny_input)
    print(f"  Tiny input clamped eta: {clamped_eta_tiny:.4f} (Should be >= 0.05)")
    assert clamped_eta_tiny >= 0.05 - 1e-5, f"Gain clamping failed for lower bound: {clamped_eta_tiny}"
    print("  -> PASS: Gain clamping strictly enforces [0.05, 40.0] bounds.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 5: 2-Bit Packing and Unpacking Identity
    # --------------------------------------------------------------------------
    print("\n[Test 5/15] 2-Bit Packing and Hardware Unpacking Identity...")
    N, K = 64, 1536
    raw_ternary = torch.randint(-1, 2, (N, K), dtype=torch.float32)

    # NumPy path
    packed_np = PackedTernaryStorage.pack(raw_ternary)
    unpacked_np = PackedTernaryStorage.unpack(packed_np, N, K).to(torch.float32)
    pack_err_np = torch.max(torch.abs(raw_ternary - unpacked_np)).item()
    print(f"  NumPy packing error: {pack_err_np}")
    assert pack_err_np == 0.0, f"NumPy packing corrupted weights: {pack_err_np}"

    # PyTorch path
    packed_torch = PackedTernaryStorage.pack_torch(raw_ternary)
    unpacked_torch = PackedTernaryStorage.unpack_torch(packed_torch, N, K).to(torch.float32)
    pack_err_torch = torch.max(torch.abs(raw_ternary - unpacked_torch)).item()
    print(f"  Torch packing error: {pack_err_torch}")
    assert pack_err_torch == 0.0, f"Torch packing corrupted weights: {pack_err_torch}"
    print("  -> PASS: Exact 4:1 compression with zero data loss.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 6: Absmean Ternary STE Quantizer & Gradient Flow
    # --------------------------------------------------------------------------
    print("\n[Test 6/15] Absmean Ternary STE Quantizer & Gradient Flow...")
    w_latent = nn.Parameter(torch.randn(32, 128, requires_grad=True))
    w_tern, gamma = AbsmeanWeightQuantizer.apply(w_latent)

    unique_vals = torch.unique(w_tern).tolist()
    print(f"  Ternary quantized unique values: {sorted(unique_vals)}")
    for val in unique_vals:
        assert val in [-1.0, 0.0, 1.0], f"Unexpected quantized value: {val}"

    # Test STE gradient propagation
    loss = (w_tern * 2.0).sum()
    loss.backward()
    assert w_latent.grad is not None, "Gradient did not flow to latent weights"
    assert torch.any(w_latent.grad != 0), "STE produced zero gradients"
    print("  -> PASS: Absmean STE maps strictly to {-1, 0, 1} and flows gradients.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 7: Symmetric Activation Quantizer & STE Flow
    # --------------------------------------------------------------------------
    print("\n[Test 7/15] Symmetric Activation INT8 Quantizer & STE Flow...")
    act = torch.randn(4, 256, requires_grad=True)
    act_q, scale = SymmetricActivationQuantizer.apply(act, 1e-6, 0.05, 40.0)

    assert act_q.min().item() >= -127.0, f"Activation lower bound violated: {act_q.min().item()}"
    assert act_q.max().item() <= 127.0, f"Activation upper bound violated: {act_q.max().item()}"

    loss_act = (act_q * 1.5).sum()
    loss_act.backward()
    assert act.grad is not None and torch.all(act.grad == 1.5), "Activation STE gradient mismatch"
    print("  -> PASS: Activations strictly clamped to [-127, 127] with valid STE flow.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 8: BitLinear Layer Execution & Training Convergence Step
    # --------------------------------------------------------------------------
    print("\n[Test 8/15] BitLinear Forward, Backward & Optimizer Step...")
    layer = BitLinear(in_features=1536, out_features=1536, bias=True, use_fwht=True)
    optimizer = torch.optim.SGD(layer.parameters(), lr=0.01)

    inp = torch.randn(2, 1536, requires_grad=True)
    out = layer(inp)
    assert out.shape == (2, 1536), f"Unexpected output shape: {out.shape}"

    loss_l = out.sum()
    optimizer.zero_grad()
    loss_l.backward()
    optimizer.step()
    assert layer.weight.grad is not None, "BitLinear weight grad is None"
    print("  -> PASS: BitLinear successfully executed forward, backward and optimizer update.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 9: Parity between PyTorch BitLinear and Hardware GEMV Simulator
    # --------------------------------------------------------------------------
    print("\n[Test 9/15] Numerical Parity: BitLinear vs Hardware GEMV Simulator...")
    test_layer = BitLinear(in_features=1536, out_features=1536, bias=False, use_fwht=True)
    test_layer.eval()

    test_input = torch.randn(1, 1536, dtype=torch.float32)

    with torch.no_grad():
        # PyTorch execution
        y_pytorch = test_layer(test_input).squeeze(0).numpy()

        # Manual breakdown for Hardware Simulator
        x_norm = test_layer.rms_norm(test_input)
        x_rot = blockwise_fwht_forward(x_norm, block_size=512, normalize=True)
        eta = torch.amax(torch.abs(x_rot), dim=-1, keepdim=True).clamp(min=1e-6)
        eta_clamped = torch.clamp(eta, test_layer.eta_min, test_layer.eta_max)
        scale_x = (127.0 / eta_clamped).item()
        x_q_int8 = torch.clamp(torch.round(x_rot * scale_x), -127.0, 127.0).squeeze(0).numpy().astype(np.int8)

        gamma_w = torch.mean(torch.abs(test_layer.weight)).clamp(min=1e-6).item()
        w_ternary = torch.clamp(torch.round(test_layer.weight / gamma_w), -1.0, 1.0)
        packed_w = PackedTernaryStorage.pack(w_ternary)

    # Hardware C-level simulation
    y_sim = simulate_hardware_bitnet_gemv(
        x_int8=x_q_int8,
        packed_w_uint8=packed_w,
        scale_x=scale_x,
        gamma_w=gamma_w,
        K=1536,
        N=1536
    )

    max_discrepancy = np.max(np.abs(y_pytorch - y_sim))
    mean_discrepancy = np.mean(np.abs(y_pytorch - y_sim))
    print(f"  Max discrepancy : {max_discrepancy:.6e}")
    print(f"  Mean discrepancy: {mean_discrepancy:.6e}")
    assert max_discrepancy < 1e-4, f"Hardware simulator parity mismatch: {max_discrepancy}"
    print("  -> PASS: 100% bit-level parity with hardware integer GEMV engine.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 10: SwiGLU Activation Clamping (<= 10.0)
    # --------------------------------------------------------------------------
    print("\n[Test 10/15] SwiGLU Activation Clamping Safeguard (<= 10.0)...")
    swiglu = SwiGLUClamped(d_model=1536, d_ffn=6144, clamp_val=10.0, use_bitlinear=True)

    # Inject extreme input to try to blow up activations
    extreme_x = torch.randn(2, 1536) * 100.0
    out_swiglu = swiglu(extreme_x)
    assert not torch.isnan(out_swiglu).any(), "SwiGLU produced NaN output"
    assert not torch.isinf(out_swiglu).any(), "SwiGLU produced Inf output"

    # Functional test: check intermediate activation is bounded
    gate = F.silu(swiglu.w_gate(extreme_x))
    up = swiglu.w_up(extreme_x)
    raw_act = gate * up
    clamped_act = torch.clamp(raw_act, -10.0, 10.0)
    print(f"  Raw intermediate SwiGLU max magnitude: {raw_act.abs().max().item():.2f}")
    print(f"  Clamped intermediate SwiGLU max magnitude: {clamped_act.abs().max().item():.2f}")
    assert clamped_act.abs().max().item() <= 10.0 + 1e-6, "SwiGLU clamp breached!"
    print("  -> PASS: SwiGLU activations strictly clamped within [-10.0, 10.0].")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 11: Logit Soft-Capping (<= 30.0)
    # --------------------------------------------------------------------------
    print("\n[Test 11/15] Logit Soft-Capping Safeguard (<= 30.0)...")
    soft_capper = LogitSoftCapping(cap_val=30.0)

    # Test extreme logits
    extreme_logits = torch.tensor([-10000.0, -100.0, -30.0, 0.0, 30.0, 100.0, 10000.0], dtype=torch.float32)
    capped = soft_capper(extreme_logits)

    print(f"  Input Logits : {extreme_logits.tolist()}")
    print(f"  Capped Logits: {[round(c, 4) for c in capped.tolist()]}")

    assert torch.all(torch.abs(capped) <= 30.0), "Logits exceeded 30.0 soft cap!"
    assert abs(capped[0].item() - (-30.0)) < 1e-4, "Extreme negative logit not bounded at -30.0"
    assert abs(capped[-1].item() - 30.0) < 1e-4, "Extreme positive logit not bounded at +30.0"

    # Check derivative: at 0, derivative is 1.0; at 100, derivative < 0.01; at 300+, derivative < 1e-6
    logits_diff = torch.tensor([0.0, 100.0, 300.0], requires_grad=True)
    capped_diff = soft_capper(logits_diff)
    grad_out = torch.ones_like(capped_diff)
    capped_diff.backward(grad_out)
    grad_at_zero = logits_diff.grad[0].item()
    grad_at_100 = logits_diff.grad[1].item()
    grad_at_300 = logits_diff.grad[2].item()
    print(f"  Gradient at 0.0   : {grad_at_zero:.4f} (Expected: ~1.0000)")
    print(f"  Gradient at 100.0 : {grad_at_100:.4e} (Suppressed by >99%)")
    print(f"  Gradient at 300.0 : {grad_at_300:.4e} (Expected: ~0.0000)")
    assert abs(grad_at_zero - 1.0) < 1e-3, f"Gradient at 0 violated: {grad_at_zero}"
    assert grad_at_100 < 0.01, f"Gradient at 100 not suppressed: {grad_at_100}"
    assert grad_at_300 < 1e-6, f"Gradient at 300 not saturated: {grad_at_300}"
    print("  -> PASS: Logits smoothly capped at 30.0 with gradient saturation.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 12: 64-Byte Cache-Line Memory Alignment Compliance
    # --------------------------------------------------------------------------
    print("\n[Test 12/15] 64-Byte Cache-Line Memory Alignment Verification...")
    aligned_buf = PackedTernaryStorage.allocate_aligned_uint8((1536, 1536 // 4))
    is_aligned = PackedTernaryStorage.is_64byte_aligned(aligned_buf)
    print(f"  Buffer memory pointer: {hex(aligned_buf.ctypes.data)}")
    print(f"  Pointer % 64 bytes  : {aligned_buf.ctypes.data % 64}")
    assert is_aligned, "Buffer is not 64-byte cache line aligned"
    print("  -> PASS: 64-byte cache line alignment strictly enforced.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 13: Format 1 (Rada v2 Compatibility) Packing and Unpacking Identity
    # --------------------------------------------------------------------------
    print("\n[Test 13/15] Format 1 (Rada v2 Compatibility) 2-Bit Packing Identity...")
    raw_ternary_f1 = torch.tensor([[-1.0, 0.0, 1.0, -1.0], [1.0, -1.0, 0.0, 1.0]], dtype=torch.float32)
    packed_f1 = PackedTernaryStorage.pack_torch(raw_ternary_f1, format_mode=1)
    unpacked_f1 = PackedTernaryStorage.unpack_torch(packed_f1, 2, 4, format_mode=1).to(torch.float32)
    err_f1 = torch.max(torch.abs(raw_ternary_f1 - unpacked_f1)).item()
    print(f"  Format 1 roundtrip reconstruction error: {err_f1}")
    assert err_f1 == 0.0, f"Format 1 packing corrupted weights: {err_f1}"
    print("  -> PASS: Format 1 (Rada v2) encodes {-1->2, 0->0, 1->1} with zero error.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 14: FP16 Numerical Stability with Near-Zero Inputs (No NaN / Inf)
    # --------------------------------------------------------------------------
    print("\n[Test 14/15] FP16 Numerical Stability Safeguard (Near-Zero Inputs)...")
    tiny_fp16_x = torch.zeros(2, 512, dtype=torch.float16)
    q_x, scale_fp16 = SymmetricActivationQuantizer.apply(tiny_fp16_x, 1e-6, 0.05, 40.0)
    assert not torch.isnan(q_x).any(), "FP16 SymmetricActivationQuantizer produced NaN"
    assert not torch.isinf(q_x).any(), "FP16 SymmetricActivationQuantizer produced Inf"
    assert not torch.isnan(scale_fp16).any(), "FP16 scale produced NaN"
    assert not torch.isinf(scale_fp16).any(), "FP16 scale produced Inf"
    print(f"  FP16 scale on zeros: {scale_fp16[0, 0].item():.4f} (Finite & valid)")
    print("  -> PASS: FP16 calculations are completely protected against Inf/NaN.")
    passed_tests += 1

    # --------------------------------------------------------------------------
    # Test 15: STE Gradient Vanishing Elimination for Small Weights
    # --------------------------------------------------------------------------
    print("\n[Test 15/15] STE Gradient Vanishing Elimination for Small Weights...")
    w_tiny = nn.Parameter(torch.full((16, 16), 1e-5, requires_grad=True))
    w_tern_tiny, _ = AbsmeanWeightQuantizer.apply(w_tiny)
    loss_tiny = w_tern_tiny.sum()
    loss_tiny.backward()
    grad_tiny_mean = w_tiny.grad.abs().mean().item()
    print(f"  Gradient mean magnitude for weights=1e-5: {grad_tiny_mean:.6f}")
    assert grad_tiny_mean >= 0.5, f"Gradient vanished for small weights: {grad_tiny_mean}"
    print("  -> PASS: Compensated STE prevents gradient vanishing trap (grad magnitude >= 0.5).")
    passed_tests += 1

    print("\n" + "=" * 80)
    print(f"ALL {passed_tests}/{total_tests} QUANTIZATION & FWHT ENGINE TESTS PASSED!")
    print("=" * 80)
    return True


if __name__ == "__main__":
    success = run_quantization_engine_tests()
    if not success:
        raise SystemExit(1)