abaddon 0.2.0-rc.1

//! GPU-accelerated holographic tensor reconstruction.
//!
//! This module provides CUDA-accelerated reconstruction for HoloTensor format,
//! enabling progressive loading of neural network weights with quality proportional
//! to fragments loaded.
//!
//! ## Encoding Schemes
//!
//! ### Spectral (SHE)
//! Uses 2D IDCT (Inverse Discrete Cosine Transform) for reconstruction.
//! Each fragment contributes frequency coefficients that are accumulated
//! before final inverse transform.
//!
//! ### Random Projection (RPH)
//! Uses least-squares reconstruction from random projections.
//! Projections are accumulated and normalized for final output.
//!
//! ### Low-Rank Distributed (LRDF)
//! Uses batched outer products for rank-1 component accumulation.
//! Each fragment contributes (u, s, v) triplets that sum to the output.
//!
//! ## Architecture
//!
//! ```text
//! ┌─────────────────────────────────────────────────────────────────────────────┐
//! │                    HoloTensor GPU Reconstruction                             │
//! ├─────────────────────────────────────────────────────────────────────────────┤
//! │                                                                             │
//! │  ┌──────────┐    ┌──────────┐    ┌──────────┐    ┌──────────┐              │
//! │  │Fragment 0│    │Fragment 1│    │Fragment 2│    │Fragment N│              │
//! │  └────┬─────┘    └────┬─────┘    └────┬─────┘    └────┬─────┘              │
//! │       │               │               │               │                     │
//! │       ▼               ▼               ▼               ▼                     │
//! │  ┌─────────────────────────────────────────────────────────────┐           │
//! │  │              Accumulator Buffer (GPU)                        │           │
//! │  │  Spectral: DCT coefficients accumulating                     │           │
//! │  │  RPH: Projection sum accumulating                            │           │
//! │  │  LRDF: Rank-1 matrices accumulating                          │           │
//! │  └─────────────────────────────────────────────────────────────┘           │
//! │                               │                                             │
//! │                               ▼ (when quality threshold reached)            │
//! │  ┌─────────────────────────────────────────────────────────────┐           │
//! │  │              GPU Reconstruction Kernel                       │           │
//! │  │  Spectral: IDCT-2D                                           │           │
//! │  │  RPH: Final normalization                                    │           │
//! │  │  LRDF: No finalization needed                                │           │
//! │  └─────────────────────────────────────────────────────────────┘           │
//! │                               │                                             │
//! │                               ▼                                             │
//! │  ┌─────────────────────────────────────────────────────────────┐           │
//! │  │              Reconstructed Tensor (GPU memory)               │           │
//! │  └─────────────────────────────────────────────────────────────┘           │
//! │                                                                             │
//! └─────────────────────────────────────────────────────────────────────────────┘
//! ```

/// CUDA-accelerated holographic tensor reconstruction.
#[cfg(feature = "cuda")]
pub mod cuda {
    use std::sync::Arc;

    use cudarc::driver::{CudaDevice, CudaSlice, DeviceSlice, LaunchAsync, LaunchConfig};
    use cudarc::nvrtc::Ptx;
    use haagenti::holotensor::{HoloFragment, HoloTensorHeader, HolographicEncoding, QualityCurve};

    // Haagenti's FFT-capable DCT context (when haagenti-gpu feature is enabled)
    #[cfg(feature = "haagenti-gpu")]
    use haagenti_cuda::dct_gpu::GpuDctContext;

    // ==================== Error Types ====================

    /// Errors from GPU holographic operations.
    #[derive(Debug, thiserror::Error)]
    pub enum GpuHoloError {
        /// CUDA device initialization failed.
        #[error("Failed to initialize CUDA device {device_id}: {message}")]
        DeviceInit {
            /// CUDA device ID.
            device_id: usize,
            /// Error message.
            message: String,
        },

        /// Kernel loading failed.
        #[error("Failed to load kernel: {message}")]
        KernelLoad {
            /// Error message.
            message: String,
        },

        /// Kernel not loaded.
        #[error("Kernel not loaded: {kernel}")]
        KernelNotLoaded {
            /// Kernel name.
            kernel: String,
        },

        /// Kernel execution failed.
        #[error("Kernel execution failed: {message}")]
        KernelExec {
            /// Error message.
            message: String,
        },

        /// Memory allocation failed.
        #[error("Memory allocation failed: {message}")]
        MemoryAlloc {
            /// Error message.
            message: String,
        },

        /// Memory copy failed.
        #[error("Memory copy failed: {message}")]
        MemoryCopy {
            /// Error message.
            message: String,
        },

        /// Synchronization failed.
        #[error("Synchronization failed: {message}")]
        Synchronize {
            /// Error message.
            message: String,
        },

        /// Invalid input data.
        #[error("Invalid input: {message}")]
        InvalidInput {
            /// Error message.
            message: String,
        },

        /// Insufficient fragments for reconstruction.
        #[error("Insufficient fragments: need at least {min_required}, have {available}")]
        InsufficientFragments {
            /// Minimum fragments required.
            min_required: u16,
            /// Number of fragments available.
            available: u16,
        },

        /// Unsupported encoding.
        #[error("Unsupported encoding: {encoding:?}")]
        UnsupportedEncoding {
            /// The unsupported encoding type.
            encoding: HolographicEncoding,
        },

        /// Fragment decode error.
        #[error("Fragment decode error: {message}")]
        FragmentDecode {
            /// Error message.
            message: String,
        },

        /// Stream creation failed.
        #[error("Failed to create CUDA stream {stream_id}: {message}")]
        StreamCreate {
            /// Stream index.
            stream_id: usize,
            /// Error message.
            message: String,
        },

        /// Quality target not reached.
        #[error("Quality target {target} not reached, current quality: {current}")]
        QualityNotReached {
            /// Target quality score.
            target: f32,
            /// Current quality score.
            current: f32,
        },
    }

    // ==================== Accumulator State ====================

    /// Accumulator state for progressive reconstruction.
    #[derive(Debug)]
    pub enum AccumulatorState {
        /// Spectral: accumulated DCT coefficients
        Spectral {
            /// Accumulated frequency coefficients
            coefficients: CudaSlice<f32>,
            /// Mask of which coefficients are present
            present_mask: CudaSlice<u8>,
            /// Width of tensor
            width: usize,
            /// Height of tensor
            height: usize,
        },
        /// RPH: accumulated projections
        RandomProjection {
            /// Sum of projections
            projection_sum: CudaSlice<f32>,
            /// Number of projections accumulated
            num_projections: u32,
            /// Projection dimension
            proj_dim: usize,
            /// Output dimension
            output_dim: usize,
            /// Base seed for projection generation
            seed: u64,
        },
        /// LRDF: accumulated rank-1 matrices
        LowRankDistributed {
            /// Accumulated output matrix
            output: CudaSlice<f32>,
            /// Number of components accumulated
            num_components: u32,
            /// Number of matrix rows.
            rows: usize,
            /// Number of matrix columns.
            cols: usize,
        },
    }

    // ==================== Kernel Configuration ====================

    /// Kernel launch configuration for optimal GPU utilization.
    ///
    /// These values are tuned based on device properties (SM count, warp size,
    /// max threads per block, shared memory size).
    #[derive(Debug, Clone, Copy)]
    pub struct KernelConfig {
        /// Block size for 1D kernels (accumulate, finalize).
        pub block_size_1d: u32,
        /// Block size for 2D kernels (IDCT, outer product).
        pub block_size_2d: u32,
        /// Maximum shared memory per block (bytes).
        pub shared_mem_limit: u32,
        /// Number of streaming multiprocessors.
        pub sm_count: u32,
        /// Warp size (typically 32).
        pub warp_size: u32,
    }

    impl Default for KernelConfig {
        fn default() -> Self {
            // Conservative defaults for SM 7.0+ (Volta and newer)
            Self {
                block_size_1d: 256,
                block_size_2d: 16,
                shared_mem_limit: 48 * 1024, // 48 KB
                sm_count: 80,
                warp_size: 32,
            }
        }
    }

    impl KernelConfig {
        /// Create configuration optimized for 1D operations.
        pub fn for_1d_operation(&self, num_elements: usize) -> (u32, u32) {
            let grid = ((num_elements as u32) + self.block_size_1d - 1) / self.block_size_1d;
            (grid.max(1), self.block_size_1d)
        }

        /// Create configuration optimized for 2D operations (matrix ops).
        pub fn for_2d_operation(&self, rows: usize, cols: usize) -> ((u32, u32), (u32, u32)) {
            let grid_x = ((cols as u32) + self.block_size_2d - 1) / self.block_size_2d;
            let grid_y = ((rows as u32) + self.block_size_2d - 1) / self.block_size_2d;
            (
                (grid_x.max(1), grid_y.max(1)),
                (self.block_size_2d, self.block_size_2d),
            )
        }

        /// Calculate optimal block size for given shared memory requirement.
        pub fn optimal_block_for_shared_mem(&self, shared_per_thread: u32) -> u32 {
            let max_threads = self.shared_mem_limit / shared_per_thread.max(1);
            // Round down to nearest warp size multiple
            let threads = max_threads.min(1024);
            (threads / self.warp_size) * self.warp_size
        }
    }

    // ==================== GPU Tensor ====================

    /// A tensor stored in GPU memory.
    ///
    /// Wraps a CUDA buffer with shape information and provides
    /// convenience methods for data transfer.
    pub struct GpuTensor {
        data: CudaSlice<f32>,
        rows: usize,
        cols: usize,
        device: Arc<CudaDevice>,
    }

    impl GpuTensor {
        /// Creates a new GPU tensor.
        pub fn new(
            data: CudaSlice<f32>,
            rows: usize,
            cols: usize,
            device: Arc<CudaDevice>,
        ) -> Self {
            Self {
                data,
                rows,
                cols,
                device,
            }
        }

        /// Returns the number of rows.
        pub fn rows(&self) -> usize {
            self.rows
        }

        /// Returns the number of columns.
        pub fn cols(&self) -> usize {
            self.cols
        }

        /// Returns the total number of elements.
        pub fn len(&self) -> usize {
            self.rows * self.cols
        }

        /// Returns true if the tensor is empty.
        pub fn is_empty(&self) -> bool {
            self.rows == 0 || self.cols == 0
        }

        /// Copies the tensor data to host memory.
        pub fn to_host(&self) -> Result<Vec<f32>, GpuHoloError> {
            let mut host = vec![0.0f32; self.len()];
            self.device
                .dtoh_sync_copy_into(&self.data, &mut host)
                .map_err(|e| GpuHoloError::MemoryCopy {
                    message: e.to_string(),
                })?;
            Ok(host)
        }

        /// Returns a reference to the raw GPU data.
        pub fn raw(&self) -> &CudaSlice<f32> {
            &self.data
        }

        /// Converts to a Vec, copying from GPU to CPU.
        pub fn to_vec(&self) -> Vec<f32> {
            self.to_host().unwrap_or_default()
        }
    }

    // ==================== GPU Context ====================

    /// GPU holographic reconstruction context.
    ///
    /// Manages CUDA kernels and state for holographic tensor reconstruction.
    pub struct GpuHoloContext {
        device: Arc<CudaDevice>,
        device_id: usize,
        spectral_kernel_loaded: bool,
        rph_kernel_loaded: bool,
        lrdf_kernel_loaded: bool,
        /// Kernel configuration based on device properties.
        kernel_config: KernelConfig,
        /// Haagenti's FFT-capable DCT context for large tensor reconstruction.
        /// Uses FFT-based O(n log n) IDCT for dimensions > 4096, providing 40-80x speedup.
        #[cfg(feature = "haagenti-gpu")]
        fft_dct_ctx: Mutex<Option<GpuDctContext>>,
    }

    impl GpuHoloContext {
        /// Threshold for using FFT-based IDCT (dimensions > this use FFT).
        /// For tensors like MLP layers (28672x8192), FFT provides 40-80x speedup.
        /// TEMP: Set to very high value to disable FFT path (GpuDctContext hangs in WSL)
        /// TODO: Debug GpuDctContext::with_device() hang and re-enable FFT for perf gain
        pub const FFT_THRESHOLD: usize = 1_000_000_000;

        /// Creates a new GPU holographic context for the specified device.
        pub fn new(device_id: usize) -> Result<Self, GpuHoloError> {
            let device = CudaDevice::new(device_id).map_err(|e| GpuHoloError::DeviceInit {
                device_id,
                message: e.to_string(),
            })?;

            // Query device properties for optimal kernel configuration
            let kernel_config = KernelConfig::default();

            Ok(Self {
                device,
                device_id,
                spectral_kernel_loaded: false,
                rph_kernel_loaded: false,
                lrdf_kernel_loaded: false,
                kernel_config,
                #[cfg(feature = "haagenti-gpu")]
                fft_dct_ctx: Mutex::new(None), // Lazily initialized on first use
            })
        }

        /// Creates a new context using an existing CUDA device.
        pub fn with_device(device: Arc<CudaDevice>, device_id: usize) -> Self {
            Self {
                device,
                device_id,
                spectral_kernel_loaded: false,
                rph_kernel_loaded: false,
                lrdf_kernel_loaded: false,
                kernel_config: KernelConfig::default(),
                #[cfg(feature = "haagenti-gpu")]
                fft_dct_ctx: Mutex::new(None),
            }
        }

        /// Returns the kernel configuration.
        pub fn kernel_config(&self) -> &KernelConfig {
            &self.kernel_config
        }

        /// Sets a custom kernel configuration.
        pub fn with_kernel_config(mut self, config: KernelConfig) -> Self {
            self.kernel_config = config;
            self
        }

        /// Returns the CUDA device ID.
        pub fn device_id(&self) -> usize {
            self.device_id
        }

        /// Returns a reference to the CUDA device.
        pub fn device(&self) -> &Arc<CudaDevice> {
            &self.device
        }

        // ==================== Kernel Loading ====================

        /// Loads the spectral (IDCT) reconstruction kernels.
        pub fn load_spectral_kernel(&mut self) -> Result<(), GpuHoloError> {
            if self.spectral_kernel_loaded {
                return Ok(());
            }

            let ptx = Ptx::from_src(SPECTRAL_KERNEL_PTX);
            self.device
                .load_ptx(
                    ptx,
                    "holo_spectral",
                    &[
                        "holo_spectral_accumulate",
                        "holo_spectral_idct_1d_rows",
                        "holo_spectral_idct_1d_cols",
                        "holo_spectral_idct_2d", // DD-8 stub: retained for future Nihil optimization
                    ],
                )
                .map_err(|e| GpuHoloError::KernelLoad {
                    message: e.to_string(),
                })?;

            self.spectral_kernel_loaded = true;
            Ok(())
        }

        /// Loads the RPH (random projection) reconstruction kernels.
        pub fn load_rph_kernel(&mut self) -> Result<(), GpuHoloError> {
            if self.rph_kernel_loaded {
                return Ok(());
            }

            let ptx = Ptx::from_src(RPH_KERNEL_PTX);
            self.device
                .load_ptx(
                    ptx,
                    "holo_rph",
                    &[
                        "holo_rph_accumulate",
                        "holo_rph_finalize",
                        "holo_rph_generate_projection", // DD-8 stub: retained for future RPH batching
                    ],
                )
                .map_err(|e| GpuHoloError::KernelLoad {
                    message: e.to_string(),
                })?;

            self.rph_kernel_loaded = true;
            Ok(())
        }

        /// Loads the LRDF (low-rank) reconstruction kernels.
        pub fn load_lrdf_kernel(&mut self) -> Result<(), GpuHoloError> {
            if self.lrdf_kernel_loaded {
                return Ok(());
            }

            let ptx = Ptx::from_src(LRDF_KERNEL_PTX);
            self.device
                .load_ptx(
                    ptx,
                    "holo_lrdf",
                    &["holo_lrdf_outer_product", "holo_lrdf_outer_product_batched"],
                )
                .map_err(|e| GpuHoloError::KernelLoad {
                    message: e.to_string(),
                })?;

            self.lrdf_kernel_loaded = true;
            Ok(())
        }

        /// Loads all holographic reconstruction kernels.
        pub fn load_all_kernels(&mut self) -> Result<(), GpuHoloError> {
            self.load_spectral_kernel()?;
            self.load_rph_kernel()?;
            self.load_lrdf_kernel()?;
            Ok(())
        }

        // ==================== Accumulator Creation ====================

        /// Creates an accumulator for progressive reconstruction.
        pub fn create_accumulator(
            &self,
            header: &HoloTensorHeader,
        ) -> Result<AccumulatorState, GpuHoloError> {
            let (width, height) = Self::extract_2d_dims(header)?;
            let total_size = width * height;

            match header.encoding {
                HolographicEncoding::Spectral => {
                    // Allocate coefficient buffer and presence mask
                    let coefficients: CudaSlice<f32> = self
                        .device
                        .alloc_zeros(total_size)
                        .map_err(|e| GpuHoloError::MemoryAlloc {
                            message: e.to_string(),
                        })?;

                    let present_mask: CudaSlice<u8> =
                        self.device.alloc_zeros(total_size).map_err(|e| {
                            GpuHoloError::MemoryAlloc {
                                message: e.to_string(),
                            }
                        })?;

                    Ok(AccumulatorState::Spectral {
                        coefficients,
                        present_mask,
                        width,
                        height,
                    })
                },

                HolographicEncoding::RandomProjection => {
                    // For RPH, we accumulate projected values
                    let proj_dim = Self::compute_projection_dim(total_size);
                    let projection_sum: CudaSlice<f32> = self
                        .device
                        .alloc_zeros(total_size)
                        .map_err(|e| GpuHoloError::MemoryAlloc {
                            message: e.to_string(),
                        })?;

                    Ok(AccumulatorState::RandomProjection {
                        projection_sum,
                        num_projections: 0,
                        proj_dim,
                        output_dim: total_size,
                        seed: header.seed,
                    })
                },

                HolographicEncoding::LowRankDistributed => {
                    // For LRDF, we accumulate the output matrix directly
                    let output: CudaSlice<f32> =
                        self.device.alloc_zeros(total_size).map_err(|e| {
                            GpuHoloError::MemoryAlloc {
                                message: e.to_string(),
                            }
                        })?;

                    Ok(AccumulatorState::LowRankDistributed {
                        output,
                        num_components: 0,
                        rows: height,
                        cols: width,
                    })
                },
            }
        }

        // ==================== Spectral Reconstruction ====================

        /// HCT3 format magic bytes: "HCT3" = 0x48435433
        const HCT3_MAGIC: u32 = 0x48435433;

        /// Accumulates a spectral fragment into the accumulator.
        /// Supports both legacy format and HCT3 format from CompressiveSpectralEncoder.
        pub fn accumulate_spectral(
            &self,
            fragment: &HoloFragment,
            accumulator: &mut AccumulatorState,
        ) -> Result<(), GpuHoloError> {
            if !self.spectral_kernel_loaded {
                return Err(GpuHoloError::KernelNotLoaded {
                    kernel: "holo_spectral".to_string(),
                });
            }

            let (coefficients, present_mask, width, height) = match accumulator {
                AccumulatorState::Spectral {
                    coefficients,
                    present_mask,
                    width,
                    height,
                } => (coefficients, present_mask, *width, *height),
                _ => {
                    return Err(GpuHoloError::InvalidInput {
                        message: "Expected spectral accumulator".to_string(),
                    })
                },
            };

            let fragment_data = &fragment.data;
            if fragment_data.len() < 4 {
                return Err(GpuHoloError::FragmentDecode {
                    message: "Fragment too small".to_string(),
                });
            }

            // Check for HCT3 format magic
            let magic = u32::from_le_bytes([
                fragment_data[0],
                fragment_data[1],
                fragment_data[2],
                fragment_data[3],
            ]);

            if magic == Self::HCT3_MAGIC {
                // HCT3 format from CompressiveSpectralEncoder
                return self.accumulate_spectral_hct3(
                    fragment,
                    coefficients,
                    present_mask,
                    width,
                    height,
                );
            }

            // Legacy format: [num_coeffs: u32][indices: u32...][values: f32...]
            let num_coeffs = magic as usize; // First 4 bytes are num_coeffs in legacy format

            let expected_size = 4 + num_coeffs * 4 + num_coeffs * 4; // header + indices + values
            if fragment_data.len() < expected_size {
                return Err(GpuHoloError::FragmentDecode {
                    message: format!(
                        "Fragment size mismatch: expected {}, got {}",
                        expected_size,
                        fragment_data.len()
                    ),
                });
            }

            // Extract indices and values
            let indices_start = 4;
            let values_start = 4 + num_coeffs * 4;

            let mut indices = Vec::with_capacity(num_coeffs);
            let mut values = Vec::with_capacity(num_coeffs);

            for i in 0..num_coeffs {
                let idx_offset = indices_start + i * 4;
                let idx = u32::from_le_bytes([
                    fragment_data[idx_offset],
                    fragment_data[idx_offset + 1],
                    fragment_data[idx_offset + 2],
                    fragment_data[idx_offset + 3],
                ]);
                indices.push(idx);

                let val_offset = values_start + i * 4;
                let val = f32::from_le_bytes([
                    fragment_data[val_offset],
                    fragment_data[val_offset + 1],
                    fragment_data[val_offset + 2],
                    fragment_data[val_offset + 3],
                ]);
                values.push(val);
            }

            self.accumulate_spectral_kernel(
                &indices,
                &values,
                coefficients,
                present_mask,
                width * height,
            )
        }

        /// Accumulates HCT3 format spectral fragment.
        /// HCT3 Fragment 0: [magic][width][height][retain_count][essential_count][detail_per_frag][bitmap][essential_f16...]
        /// HCT3 Fragment k>0: [detail_f16...]
        fn accumulate_spectral_hct3(
            &self,
            fragment: &HoloFragment,
            coefficients: &mut CudaSlice<f32>,
            present_mask: &mut CudaSlice<u8>,
            width: usize,
            height: usize,
        ) -> Result<(), GpuHoloError> {
            use half::f16;

            let data = &fragment.data;
            let n = width * height;

            if fragment.index == 0 {
                // Parse HCT3 header (24 bytes)
                if data.len() < 24 {
                    return Err(GpuHoloError::FragmentDecode {
                        message: "HCT3 fragment 0 too short for header".to_string(),
                    });
                }

                // Skip magic (already checked), read dimensions
                let hct3_width = u32::from_le_bytes([data[4], data[5], data[6], data[7]]) as usize;
                let hct3_height =
                    u32::from_le_bytes([data[8], data[9], data[10], data[11]]) as usize;
                let retain_count =
                    u32::from_le_bytes([data[12], data[13], data[14], data[15]]) as usize;
                let essential_count =
                    u32::from_le_bytes([data[16], data[17], data[18], data[19]]) as usize;
                let _detail_per_frag =
                    u32::from_le_bytes([data[20], data[21], data[22], data[23]]) as usize;

                // Validate dimensions match
                if hct3_width != width || hct3_height != height {
                    return Err(GpuHoloError::FragmentDecode {
                        message: format!(
                            "HCT3 dimensions mismatch: header ({},{}) vs accumulator ({},{})",
                            hct3_width, hct3_height, width, height
                        ),
                    });
                }

                // Read bitmap and build index map
                let bitmap_bytes = (n + 7) / 8;
                let bitmap_start = 24;
                let bitmap_end = bitmap_start + bitmap_bytes;

                if data.len() < bitmap_end {
                    return Err(GpuHoloError::FragmentDecode {
                        message: "HCT3 fragment 0 truncated (missing bitmap)".to_string(),
                    });
                }

                // Scan bitmap to get indices of retained coefficients
                let mut indices = Vec::with_capacity(retain_count);
                for i in 0..n {
                    let byte_idx = bitmap_start + i / 8;
                    let bit_idx = i % 8;
                    if (data[byte_idx] >> bit_idx) & 1 == 1 {
                        indices.push(i as u32);
                    }
                }

                if indices.len() != retain_count {
                    return Err(GpuHoloError::FragmentDecode {
                        message: format!(
                            "HCT3 bitmap mismatch: {} set bits, expected {}",
                            indices.len(),
                            retain_count
                        ),
                    });
                }

                // Read essential coefficients (f16)
                let coeff_start = bitmap_end;
                let coeff_end = coeff_start + essential_count * 2;

                if data.len() < coeff_end {
                    return Err(GpuHoloError::FragmentDecode {
                        message: "HCT3 fragment 0 truncated (missing essential coefficients)"
                            .to_string(),
                    });
                }

                let mut values = Vec::with_capacity(essential_count);
                for i in 0..essential_count {
                    let offset = coeff_start + i * 2;
                    let f16_val = f16::from_le_bytes([data[offset], data[offset + 1]]);
                    values.push(f16_val.to_f32());
                }

                // Accumulate essential coefficients (first essential_count indices)
                let essential_indices: Vec<u32> = indices[..essential_count].to_vec();
                self.accumulate_spectral_kernel(
                    &essential_indices,
                    &values,
                    coefficients,
                    present_mask,
                    n,
                )?;

                // Store index map for detail fragments (we need this state across fragments)
                // For now, we'll re-read the bitmap for detail fragments (simpler but slightly less efficient)
                Ok(())
            } else {
                // Detail fragment: just f16 coefficients
                // We need to know which indices these correspond to, which requires re-parsing fragment 0
                // For GPU efficiency, we'll fall back to CPU for HCT3 detail fragments
                // This is a limitation - full HCT3 GPU support would require stateful tracking

                // For now, return an error to trigger CPU fallback
                Err(GpuHoloError::FragmentDecode {
                    message:
                        "HCT3 detail fragments not yet supported in GPU path - use CPU fallback"
                            .to_string(),
                })
            }
        }

        /// Common kernel invocation for spectral accumulation.
        fn accumulate_spectral_kernel(
            &self,
            indices: &[u32],
            values: &[f32],
            coefficients: &mut CudaSlice<f32>,
            present_mask: &mut CudaSlice<u8>,
            total_size: usize,
        ) -> Result<(), GpuHoloError> {
            let num_coeffs = indices.len();
            if num_coeffs == 0 {
                return Ok(());
            }

            // Copy to GPU
            let d_indices =
                self.device
                    .htod_copy(indices.to_vec())
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            let d_values =
                self.device
                    .htod_copy(values.to_vec())
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            // Launch accumulation kernel
            let func = self
                .device
                .get_func("holo_spectral", "holo_spectral_accumulate")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_spectral_accumulate".to_string(),
                })?;

            let threads_per_block = 256u32;
            let num_blocks = ((num_coeffs as u32) + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            unsafe {
                func.launch(
                    cfg,
                    (
                        &d_indices,
                        &d_values,
                        coefficients,
                        present_mask,
                        num_coeffs as u32,
                        total_size as u32,
                    ),
                )
            }
            .map_err(|e| GpuHoloError::KernelExec {
                message: e.to_string(),
            })?;

            Ok(())
        }

        /// Performs 2D IDCT reconstruction from accumulated spectral coefficients.
        ///
        /// For large tensors (dimensions > 4096), uses haagenti's FFT-based O(n log n) IDCT
        /// which provides 40-80x speedup over the direct O(n²) kernel.
        pub fn finalize_spectral(
            &self,
            accumulator: &AccumulatorState,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            let (coefficients, _present_mask, width, height) = match accumulator {
                AccumulatorState::Spectral {
                    coefficients,
                    present_mask,
                    width,
                    height,
                } => (coefficients, present_mask, *width, *height),
                _ => {
                    return Err(GpuHoloError::InvalidInput {
                        message: "Expected spectral accumulator".to_string(),
                    })
                },
            };

            // Use FFT-based IDCT for large tensors (40-80x faster)
            #[cfg(feature = "haagenti-gpu")]
            if width > Self::FFT_THRESHOLD || height > Self::FFT_THRESHOLD {
                return self.finalize_spectral_fft(coefficients, width, height);
            }

            // Fall back to direct PTX kernel for small tensors
            self.finalize_spectral_direct(coefficients, width, height)
        }

        /// FFT-based IDCT for large tensors using haagenti-cuda.
        ///
        /// Uses cuFFT-based O(n log n) algorithm for 40-80x speedup on large tensors.
        /// For a 28672x8192 tensor: ~1.3s direct → ~0.03s FFT.
        #[cfg(feature = "haagenti-gpu")]
        fn finalize_spectral_fft(
            &self,
            coefficients: &CudaSlice<f32>,
            width: usize,
            height: usize,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            let total_size = width * height;

            // Copy coefficients from GPU to host
            let mut host_coeffs = vec![0.0f32; total_size];
            self.device
                .dtoh_sync_copy_into(coefficients, &mut host_coeffs)
                .map_err(|e| GpuHoloError::MemoryCopy {
                    message: format!("FFT IDCT: failed to copy coefficients to host: {}", e),
                })?;

            // Initialize or get the FFT DCT context
            let mut ctx_guard = self.fft_dct_ctx.lock();
            let dct_ctx = ctx_guard.get_or_insert_with(|| {
                tracing::info!(
                    "Initializing haagenti FFT-IDCT context for large tensor reconstruction"
                );
                GpuDctContext::with_device(self.device.clone())
                    .expect("Failed to create GpuDctContext")
            });

            // Perform FFT-based IDCT
            let reconstructed = dct_ctx.idct_2d(&host_coeffs, width, height).map_err(|e| {
                GpuHoloError::KernelExec {
                    message: format!("FFT IDCT failed: {}", e),
                }
            })?;

            // Copy result back to GPU
            let output = self.device.htod_sync_copy(&reconstructed).map_err(|e| {
                GpuHoloError::MemoryCopy {
                    message: format!("FFT IDCT: failed to copy result to GPU: {}", e),
                }
            })?;

            tracing::debug!(width, height, "FFT-based IDCT complete for large tensor");

            Ok(output)
        }

        /// Direct PTX kernel IDCT for small/medium tensors.
        fn finalize_spectral_direct(
            &self,
            coefficients: &CudaSlice<f32>,
            width: usize,
            height: usize,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            if !self.spectral_kernel_loaded {
                return Err(GpuHoloError::KernelNotLoaded {
                    kernel: "holo_spectral".to_string(),
                });
            }

            let total_size = width * height;

            // Allocate intermediate and output buffers
            let temp: CudaSlice<f32> =
                self.device
                    .alloc_zeros(total_size)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            let output: CudaSlice<f32> =
                self.device
                    .alloc_zeros(total_size)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            // IDCT rows
            let func_rows = self
                .device
                .get_func("holo_spectral", "holo_spectral_idct_1d_rows")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_spectral_idct_1d_rows".to_string(),
                })?;

            let threads_per_block = 256u32;
            let num_blocks = ((height as u32) + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: (width * 4) as u32, // Shared memory for one row
            };

            unsafe { func_rows.launch(cfg, (coefficients, &temp, width as u32, height as u32)) }
                .map_err(|e| GpuHoloError::KernelExec {
                    message: e.to_string(),
                })?;

            // IDCT columns
            let func_cols = self
                .device
                .get_func("holo_spectral", "holo_spectral_idct_1d_cols")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_spectral_idct_1d_cols".to_string(),
                })?;

            let num_blocks_cols = ((width as u32) + threads_per_block - 1) / threads_per_block;

            let cfg_cols = LaunchConfig {
                grid_dim: (num_blocks_cols.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: (height * 4) as u32,
            };

            unsafe { func_cols.launch(cfg_cols, (&temp, &output, width as u32, height as u32)) }
                .map_err(|e| GpuHoloError::KernelExec {
                    message: e.to_string(),
                })?;

            self.device
                .synchronize()
                .map_err(|e| GpuHoloError::Synchronize {
                    message: e.to_string(),
                })?;

            Ok(output)
        }

        // ==================== RPH Reconstruction ====================

        /// Accumulates an RPH fragment into the accumulator.
        pub fn accumulate_rph(
            &self,
            fragment: &HoloFragment,
            accumulator: &mut AccumulatorState,
        ) -> Result<(), GpuHoloError> {
            if !self.rph_kernel_loaded {
                return Err(GpuHoloError::KernelNotLoaded {
                    kernel: "holo_rph".to_string(),
                });
            }

            let (projection_sum, _num_projections, _proj_dim, output_dim, seed) = match accumulator
            {
                AccumulatorState::RandomProjection {
                    projection_sum,
                    num_projections,
                    proj_dim,
                    output_dim,
                    seed,
                } => (
                    projection_sum,
                    num_projections,
                    *proj_dim,
                    *output_dim,
                    *seed,
                ),
                _ => {
                    return Err(GpuHoloError::InvalidInput {
                        message: "Expected RPH accumulator".to_string(),
                    })
                },
            };

            // Parse fragment data: [proj_dim: u32][seed_offset: u64][projection: f32...]
            let fragment_data = &fragment.data;
            if fragment_data.len() < 12 {
                return Err(GpuHoloError::FragmentDecode {
                    message: "Fragment too small for RPH".to_string(),
                });
            }

            let frag_proj_dim = u32::from_le_bytes([
                fragment_data[0],
                fragment_data[1],
                fragment_data[2],
                fragment_data[3],
            ]) as usize;

            let frag_seed_offset = u64::from_le_bytes([
                fragment_data[4],
                fragment_data[5],
                fragment_data[6],
                fragment_data[7],
                fragment_data[8],
                fragment_data[9],
                fragment_data[10],
                fragment_data[11],
            ]);

            let projection_start = 12;
            let expected_size = 12 + frag_proj_dim * 4;
            if fragment_data.len() < expected_size {
                return Err(GpuHoloError::FragmentDecode {
                    message: format!(
                        "Fragment size mismatch for RPH: expected {}, got {}",
                        expected_size,
                        fragment_data.len()
                    ),
                });
            }

            // Extract projection values
            let mut projection = Vec::with_capacity(frag_proj_dim);
            for i in 0..frag_proj_dim {
                let offset = projection_start + i * 4;
                let val = f32::from_le_bytes([
                    fragment_data[offset],
                    fragment_data[offset + 1],
                    fragment_data[offset + 2],
                    fragment_data[offset + 3],
                ]);
                projection.push(val);
            }

            // Copy projection to GPU
            let d_projection =
                self.device
                    .htod_copy(projection)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            // Launch RPH accumulation kernel
            // This generates the projection matrix on-the-fly and accumulates
            let func = self
                .device
                .get_func("holo_rph", "holo_rph_accumulate")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_rph_accumulate".to_string(),
                })?;

            let threads_per_block = 256u32;
            let num_blocks = ((output_dim as u32) + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            let fragment_seed = seed.wrapping_add(frag_seed_offset);

            unsafe {
                func.launch(
                    cfg,
                    (
                        &d_projection,
                        projection_sum,
                        frag_proj_dim as u32,
                        output_dim as u32,
                        fragment_seed,
                    ),
                )
            }
            .map_err(|e| GpuHoloError::KernelExec {
                message: e.to_string(),
            })?;

            // Update projection count
            if let AccumulatorState::RandomProjection {
                num_projections, ..
            } = accumulator
            {
                *num_projections += 1;
            }

            Ok(())
        }

        /// Finalizes RPH reconstruction by normalizing accumulated projections.
        pub fn finalize_rph(
            &self,
            accumulator: &AccumulatorState,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            if !self.rph_kernel_loaded {
                return Err(GpuHoloError::KernelNotLoaded {
                    kernel: "holo_rph".to_string(),
                });
            }

            let (projection_sum, num_projections, output_dim) = match accumulator {
                AccumulatorState::RandomProjection {
                    projection_sum,
                    num_projections,
                    output_dim,
                    ..
                } => (projection_sum, *num_projections, *output_dim),
                _ => {
                    return Err(GpuHoloError::InvalidInput {
                        message: "Expected RPH accumulator".to_string(),
                    })
                },
            };

            if num_projections == 0 {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: 1,
                    available: 0,
                });
            }

            // Allocate output buffer
            let output: CudaSlice<f32> =
                self.device
                    .alloc_zeros(output_dim)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            // Launch finalization kernel (normalize by projection count)
            let func = self
                .device
                .get_func("holo_rph", "holo_rph_finalize")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_rph_finalize".to_string(),
                })?;

            let threads_per_block = 256u32;
            let num_blocks = ((output_dim as u32) + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            unsafe {
                func.launch(
                    cfg,
                    (projection_sum, &output, output_dim as u32, num_projections),
                )
            }
            .map_err(|e| GpuHoloError::KernelExec {
                message: e.to_string(),
            })?;

            self.device
                .synchronize()
                .map_err(|e| GpuHoloError::Synchronize {
                    message: e.to_string(),
                })?;

            Ok(output)
        }

        // ==================== LRDF Reconstruction ====================

        /// Accumulates an LRDF fragment into the accumulator.
        pub fn accumulate_lrdf(
            &self,
            fragment: &HoloFragment,
            accumulator: &mut AccumulatorState,
        ) -> Result<(), GpuHoloError> {
            if !self.lrdf_kernel_loaded {
                return Err(GpuHoloError::KernelNotLoaded {
                    kernel: "holo_lrdf".to_string(),
                });
            }

            let (output, rows, cols) = match accumulator {
                AccumulatorState::LowRankDistributed {
                    ref mut output,
                    rows,
                    cols,
                    ..
                } => (output, *rows, *cols),
                _ => {
                    return Err(GpuHoloError::InvalidInput {
                        message: "Expected LRDF accumulator".to_string(),
                    })
                },
            };

            // Parse fragment data (haagenti format):
            // [rows: u32][cols: u32][num_components: u32][components...]
            // Each component: [sigma: f32][u: f32 * rows][v: f32 * cols]
            let fragment_data = &fragment.data;
            if fragment_data.len() < 12 {
                return Err(GpuHoloError::FragmentDecode {
                    message: "Fragment too small for LRDF header".to_string(),
                });
            }

            // Read header: rows, cols, num_components
            let frag_rows = u32::from_le_bytes([
                fragment_data[0],
                fragment_data[1],
                fragment_data[2],
                fragment_data[3],
            ]) as usize;

            let frag_cols = u32::from_le_bytes([
                fragment_data[4],
                fragment_data[5],
                fragment_data[6],
                fragment_data[7],
            ]) as usize;

            let num_comps = u32::from_le_bytes([
                fragment_data[8],
                fragment_data[9],
                fragment_data[10],
                fragment_data[11],
            ]) as usize;

            // Validate dimensions match accumulator
            if frag_rows != rows || frag_cols != cols {
                return Err(GpuHoloError::FragmentDecode {
                    message: format!(
                        "Fragment dimensions {}x{} don't match accumulator {}x{}",
                        frag_rows, frag_cols, rows, cols
                    ),
                });
            }

            // Check for RAW format marker (num_components == 0xFFFFFFFF)
            // RAW format stores raw f32 data directly instead of SVD components
            const RAW_FORMAT_MARKER: usize = 0xFFFFFFFF;
            if num_comps == RAW_FORMAT_MARKER {
                // RAW format: [header: 12 bytes][raw f32 data: rows*cols*4 bytes]
                let expected_size = 12 + rows * cols * 4;
                if fragment_data.len() < expected_size {
                    return Err(GpuHoloError::FragmentDecode {
                        message: format!(
                            "Fragment size mismatch for RAW LRDF: expected {}, got {}",
                            expected_size,
                            fragment_data.len()
                        ),
                    });
                }

                // Parse raw f32 data from fragment
                let raw_data = &fragment_data[12..];
                let total_elements = rows * cols;
                let mut host_data = Vec::with_capacity(total_elements);
                for i in 0..total_elements {
                    let offset = i * 4;
                    let val = f32::from_le_bytes([
                        raw_data[offset],
                        raw_data[offset + 1],
                        raw_data[offset + 2],
                        raw_data[offset + 3],
                    ]);
                    host_data.push(val);
                }

                // Copy directly to GPU output buffer
                self.device
                    .htod_sync_copy_into(&host_data, output)
                    .map_err(|e| GpuHoloError::MemoryCopy {
                        message: format!("Failed to copy RAW data to GPU: {}", e),
                    })?;

                // Mark as having 1 "component" for finalize check
                if let AccumulatorState::LowRankDistributed { num_components, .. } = accumulator {
                    *num_components = 1;
                }

                return Ok(());
            }

            // SVD format: validate expected size
            let component_size = 4 + rows * 4 + cols * 4; // sigma + u + v
            let expected_size = 12 + num_comps * component_size; // 12-byte header
            if fragment_data.len() < expected_size {
                return Err(GpuHoloError::FragmentDecode {
                    message: format!(
                        "Fragment size mismatch for LRDF SVD: expected {}, got {}",
                        expected_size,
                        fragment_data.len()
                    ),
                });
            }

            // Process each SVD component (skip 12-byte header)
            let mut offset = 12;
            for _ in 0..num_comps {
                // Extract sigma
                let sigma = f32::from_le_bytes([
                    fragment_data[offset],
                    fragment_data[offset + 1],
                    fragment_data[offset + 2],
                    fragment_data[offset + 3],
                ]);
                offset += 4;

                // Extract u vector
                let mut u = Vec::with_capacity(rows);
                for _ in 0..rows {
                    let val = f32::from_le_bytes([
                        fragment_data[offset],
                        fragment_data[offset + 1],
                        fragment_data[offset + 2],
                        fragment_data[offset + 3],
                    ]);
                    u.push(val);
                    offset += 4;
                }

                // Extract v vector
                let mut v = Vec::with_capacity(cols);
                for _ in 0..cols {
                    let val = f32::from_le_bytes([
                        fragment_data[offset],
                        fragment_data[offset + 1],
                        fragment_data[offset + 2],
                        fragment_data[offset + 3],
                    ]);
                    v.push(val);
                    offset += 4;
                }

                // Copy to GPU
                let d_u = self
                    .device
                    .htod_copy(u)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

                let d_v = self
                    .device
                    .htod_copy(v)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

                // Launch outer product kernel
                let func = self
                    .device
                    .get_func("holo_lrdf", "holo_lrdf_outer_product")
                    .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                        kernel: "holo_lrdf_outer_product".to_string(),
                    })?;

                // 2D grid for outer product
                let block_size = 16u32;
                let grid_x = ((cols as u32) + block_size - 1) / block_size;
                let grid_y = ((rows as u32) + block_size - 1) / block_size;

                let cfg = LaunchConfig {
                    grid_dim: (grid_x.max(1), grid_y.max(1), 1),
                    block_dim: (block_size, block_size, 1),
                    shared_mem_bytes: 0,
                };

                unsafe {
                    func.clone().launch(
                        cfg,
                        (&d_u, &d_v, &mut *output, sigma, rows as u32, cols as u32),
                    )
                }
                .map_err(|e| GpuHoloError::KernelExec {
                    message: e.to_string(),
                })?;
            }

            // Update component count
            if let AccumulatorState::LowRankDistributed { num_components, .. } = accumulator {
                *num_components += num_comps as u32;
            }

            Ok(())
        }

        /// Batched LRDF accumulation: processes all SVD components in a single kernel launch.
        ///
        /// Packs all components' u, v, sigma into contiguous GPU arrays and calls
        /// `holo_lrdf_outer_product_batched`. More efficient than the per-component
        /// approach for fragments with many components.
        pub fn accumulate_lrdf_batched(
            &self,
            fragment: &HoloFragment,
            accumulator: &mut AccumulatorState,
        ) -> Result<(), GpuHoloError> {
            if !self.lrdf_kernel_loaded {
                return Err(GpuHoloError::KernelNotLoaded {
                    kernel: "holo_lrdf".to_string(),
                });
            }

            let (output, rows, cols) = match accumulator {
                AccumulatorState::LowRankDistributed {
                    ref mut output,
                    rows,
                    cols,
                    ..
                } => (output, *rows, *cols),
                _ => {
                    return Err(GpuHoloError::InvalidInput {
                        message: "Expected LRDF accumulator".to_string(),
                    })
                },
            };

            let fragment_data = &fragment.data;
            if fragment_data.len() < 12 {
                return Err(GpuHoloError::FragmentDecode {
                    message: "Fragment too small for LRDF header".to_string(),
                });
            }

            let frag_rows = u32::from_le_bytes([
                fragment_data[0],
                fragment_data[1],
                fragment_data[2],
                fragment_data[3],
            ]) as usize;
            let frag_cols = u32::from_le_bytes([
                fragment_data[4],
                fragment_data[5],
                fragment_data[6],
                fragment_data[7],
            ]) as usize;
            let num_comps = u32::from_le_bytes([
                fragment_data[8],
                fragment_data[9],
                fragment_data[10],
                fragment_data[11],
            ]) as usize;

            if frag_rows != rows || frag_cols != cols {
                return Err(GpuHoloError::FragmentDecode {
                    message: format!(
                        "Fragment dimensions {}x{} don't match accumulator {}x{}",
                        frag_rows, frag_cols, rows, cols
                    ),
                });
            }

            if num_comps == 0 {
                return Ok(());
            }

            let component_size = 4 + rows * 4 + cols * 4;
            let expected_size = 12 + num_comps * component_size;
            if fragment_data.len() < expected_size {
                return Err(GpuHoloError::FragmentDecode {
                    message: format!(
                        "Fragment size mismatch for batched LRDF: expected {}, got {}",
                        expected_size,
                        fragment_data.len()
                    ),
                });
            }

            // Pack all components into contiguous arrays
            let mut all_sigma = Vec::with_capacity(num_comps);
            let mut all_u = Vec::with_capacity(num_comps * rows);
            let mut all_v = Vec::with_capacity(num_comps * cols);

            let mut offset = 12;
            for _ in 0..num_comps {
                let sigma = f32::from_le_bytes([
                    fragment_data[offset],
                    fragment_data[offset + 1],
                    fragment_data[offset + 2],
                    fragment_data[offset + 3],
                ]);
                all_sigma.push(sigma);
                offset += 4;

                for _ in 0..rows {
                    let val = f32::from_le_bytes([
                        fragment_data[offset],
                        fragment_data[offset + 1],
                        fragment_data[offset + 2],
                        fragment_data[offset + 3],
                    ]);
                    all_u.push(val);
                    offset += 4;
                }

                for _ in 0..cols {
                    let val = f32::from_le_bytes([
                        fragment_data[offset],
                        fragment_data[offset + 1],
                        fragment_data[offset + 2],
                        fragment_data[offset + 3],
                    ]);
                    all_v.push(val);
                    offset += 4;
                }
            }

            // Upload packed arrays to GPU
            let d_sigma =
                self.device
                    .htod_copy(all_sigma)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;
            let d_u = self
                .device
                .htod_copy(all_u)
                .map_err(|e| GpuHoloError::MemoryAlloc {
                    message: e.to_string(),
                })?;
            let d_v = self
                .device
                .htod_copy(all_v)
                .map_err(|e| GpuHoloError::MemoryAlloc {
                    message: e.to_string(),
                })?;

            // Single kernel launch for all components
            let func = self
                .device
                .get_func("holo_lrdf", "holo_lrdf_outer_product_batched")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_lrdf_outer_product_batched".to_string(),
                })?;

            let block_size = 16u32;
            let grid_x = ((cols as u32) + block_size - 1) / block_size;
            let grid_y = ((rows as u32) + block_size - 1) / block_size;

            let cfg = LaunchConfig {
                grid_dim: (grid_x.max(1), grid_y.max(1), 1),
                block_dim: (block_size, block_size, 1),
                shared_mem_bytes: 0,
            };

            unsafe {
                func.launch(
                    cfg,
                    (
                        &d_u,
                        &d_v,
                        &d_sigma,
                        &mut *output,
                        num_comps as u32,
                        rows as u32,
                        cols as u32,
                    ),
                )
            }
            .map_err(|e| GpuHoloError::KernelExec {
                message: e.to_string(),
            })?;

            if let AccumulatorState::LowRankDistributed { num_components, .. } = accumulator {
                *num_components += num_comps as u32;
            }

            Ok(())
        }

        /// Finalizes LRDF reconstruction (no additional processing needed).
        pub fn finalize_lrdf(
            &self,
            accumulator: &AccumulatorState,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            let (output, num_components, rows, cols) = match accumulator {
                AccumulatorState::LowRankDistributed {
                    output,
                    num_components,
                    rows,
                    cols,
                } => (output, *num_components, *rows, *cols),
                _ => {
                    return Err(GpuHoloError::InvalidInput {
                        message: "Expected LRDF accumulator".to_string(),
                    })
                },
            };

            if num_components == 0 {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: 1,
                    available: 0,
                });
            }

            self.device
                .synchronize()
                .map_err(|e| GpuHoloError::Synchronize {
                    message: e.to_string(),
                })?;

            // Clone the output (LRDF accumulates directly into output)
            let mut result: CudaSlice<f32> =
                self.device
                    .alloc_zeros(rows * cols)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            // Copy from accumulator to result
            self.device
                .dtod_copy(output, &mut result)
                .map_err(|e| GpuHoloError::MemoryCopy {
                    message: e.to_string(),
                })?;

            Ok(result)
        }

        // ==================== Unified API ====================

        /// Accumulates a fragment based on the encoding type.
        pub fn accumulate_fragment(
            &self,
            fragment: &HoloFragment,
            accumulator: &mut AccumulatorState,
            encoding: HolographicEncoding,
        ) -> Result<(), GpuHoloError> {
            match encoding {
                HolographicEncoding::Spectral => self.accumulate_spectral(fragment, accumulator),
                HolographicEncoding::RandomProjection => self.accumulate_rph(fragment, accumulator),
                HolographicEncoding::LowRankDistributed => {
                    self.accumulate_lrdf(fragment, accumulator)
                },
            }
        }

        /// Finalizes reconstruction based on the encoding type.
        pub fn finalize_reconstruction(
            &self,
            accumulator: &AccumulatorState,
            encoding: HolographicEncoding,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            match encoding {
                HolographicEncoding::Spectral => self.finalize_spectral(accumulator),
                HolographicEncoding::RandomProjection => self.finalize_rph(accumulator),
                HolographicEncoding::LowRankDistributed => self.finalize_lrdf(accumulator),
            }
        }

        /// Reconstructs a tensor from fragments.
        ///
        /// This is the high-level API for reconstruction.
        pub fn reconstruct(
            &self,
            header: &HoloTensorHeader,
            fragments: &[HoloFragment],
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            if fragments.is_empty() {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: header.quality_curve.min_fragments,
                    available: 0,
                });
            }

            if (fragments.len() as u16) < header.quality_curve.min_fragments {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: header.quality_curve.min_fragments,
                    available: fragments.len() as u16,
                });
            }

            // Create accumulator
            let mut accumulator = self.create_accumulator(header)?;

            // Accumulate all fragments
            for fragment in fragments {
                self.accumulate_fragment(fragment, &mut accumulator, header.encoding)?;
            }

            // Finalize and return
            self.finalize_reconstruction(&accumulator, header.encoding)
        }

        /// Copies reconstructed data to host memory.
        pub fn copy_to_host(&self, gpu_data: &CudaSlice<f32>) -> Result<Vec<f32>, GpuHoloError> {
            let len = gpu_data.len();
            let mut host_data = vec![0.0f32; len];
            self.device
                .dtoh_sync_copy_into(gpu_data, &mut host_data)
                .map_err(|e| GpuHoloError::MemoryCopy {
                    message: e.to_string(),
                })?;
            Ok(host_data)
        }

        // ==================== Convenience Methods ====================

        /// Reconstructs a 2D LRDF tensor from fragments.
        ///
        /// This is a convenience wrapper for LRDF reconstruction that builds
        /// the header automatically.
        pub fn reconstruct_lrdf(
            &self,
            fragments: &[HoloFragment],
            rows: usize,
            cols: usize,
        ) -> Result<GpuTensor, GpuHoloError> {
            if fragments.is_empty() {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: 1,
                    available: 0,
                });
            }

            // Build header for LRDF encoding
            let header = HoloTensorHeader::new(
                HolographicEncoding::LowRankDistributed,
                haagenti::DType::F32,
                vec![rows as u64, cols as u64],
                fragments.len() as u16,
            );

            let gpu_data = self.reconstruct(&header, fragments)?;

            Ok(GpuTensor {
                data: gpu_data,
                rows,
                cols,
                device: self.device.clone(),
            })
        }

        // ==================== Helper Functions ====================

        /// Extracts 2D dimensions from header.
        fn extract_2d_dims(header: &HoloTensorHeader) -> Result<(usize, usize), GpuHoloError> {
            if header.shape.len() < 2 {
                return Err(GpuHoloError::InvalidInput {
                    message: format!("Expected 2D tensor, got {} dimensions", header.shape.len()),
                });
            }
            Ok((header.shape[1] as usize, header.shape[0] as usize))
        }

        /// Computes projection dimension for RPH.
        fn compute_projection_dim(output_dim: usize) -> usize {
            // Use sqrt(n) as default projection dimension
            ((output_dim as f64).sqrt() as usize).max(16)
        }
    }

    // ==================== PTX Kernel Sources ====================

    /// Spectral (IDCT) kernel PTX source.
    const SPECTRAL_KERNEL_PTX: &str = r#"
.version 8.0
.target sm_89
.address_size 64

// Accumulate coefficients into buffer
// Args: indices, values, coefficients, mask, num_coeffs, buffer_size
.visible .entry holo_spectral_accumulate(
    .param .u64 indices_ptr,
    .param .u64 values_ptr,
    .param .u64 coeffs_ptr,
    .param .u64 mask_ptr,
    .param .u32 num_coeffs,
    .param .u32 buffer_size
)
{
    .reg .u32 %gid, %idx, %coeff_idx;
    .reg .u64 %indices_addr, %values_addr, %coeffs_addr, %mask_addr;
    .reg .f32 %val;
    .reg .pred %p;

    // Get thread ID: gid = blockIdx.x * blockDim.x + threadIdx.x
    mov.u32 %gid, %ctaid.x;
    mov.u32 %idx, %ntid.x;
    mul.lo.u32 %gid, %gid, %idx;
    mov.u32 %idx, %tid.x;
    add.u32 %gid, %gid, %idx;

    // Bounds check
    ld.param.u32 %idx, [num_coeffs];
    setp.ge.u32 %p, %gid, %idx;
    @%p bra EXIT;

    // Load index and value
    ld.param.u64 %indices_addr, [indices_ptr];
    ld.param.u64 %values_addr, [values_ptr];

    mad.wide.u32 %indices_addr, %gid, 4, %indices_addr;
    mad.wide.u32 %values_addr, %gid, 4, %values_addr;

    ld.global.u32 %coeff_idx, [%indices_addr];
    ld.global.f32 %val, [%values_addr];

    // Bounds check on coefficient index
    ld.param.u32 %idx, [buffer_size];
    setp.ge.u32 %p, %coeff_idx, %idx;
    @%p bra EXIT;

    // Store coefficient and mark as present
    ld.param.u64 %coeffs_addr, [coeffs_ptr];
    ld.param.u64 %mask_addr, [mask_ptr];

    mad.wide.u32 %coeffs_addr, %coeff_idx, 4, %coeffs_addr;
    cvt.u64.u32 %mask_addr, %coeff_idx;
    add.u64 %mask_addr, %mask_addr, %mask_addr;

    st.global.f32 [%coeffs_addr], %val;

EXIT:
    ret;
}

// 1D IDCT (Type-III DCT) on rows
// x[n] = sqrt(2/N) * [ X[0]/sqrt(2) + sum_{k=1}^{N-1} X[k] * cos(pi*(2n+1)*k / (2N)) ]
// One thread per row, inner loop over width columns.
// Args: input, output, width, height
.visible .entry holo_spectral_idct_1d_rows(
    .param .u64 input_ptr,
    .param .u64 output_ptr,
    .param .u32 width,
    .param .u32 height
)
{
    .reg .u32 %row, %col, %k, %w, %h, %tmp, %n2p1;
    .reg .u64 %in_base, %out_addr, %coeff_addr;
    .reg .f32 %sum, %val, %cos_val, %scale, %pi_2n, %angle, %w_f, %k_f, %n2p1_f;
    .reg .pred %p;

    // Get row index = blockIdx.x * blockDim.x + threadIdx.x
    mov.u32 %row, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %row, %row, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %row, %row, %tmp;

    ld.param.u32 %h, [height];
    setp.ge.u32 %p, %row, %h;
    @%p bra ROW_EXIT;

    ld.param.u32 %w, [width];
    ld.param.u64 %in_base, [input_ptr];
    ld.param.u64 %out_addr, [output_ptr];

    // Compute row base offsets: in_base += row * width * 4, out_addr += row * width * 4
    mul.lo.u32 %tmp, %row, %w;
    mad.wide.u32 %in_base, %tmp, 4, %in_base;
    mad.wide.u32 %out_addr, %tmp, 4, %out_addr;

    // Pre-compute scale = sqrt(2.0 / width)
    cvt.rn.f32.u32 %w_f, %w;
    mov.f32 %scale, 2.0;
    div.full.f32 %scale, %scale, %w_f;
    sqrt.approx.f32 %scale, %scale;

    // Pre-compute pi_2n = PI / (2 * width)
    mov.f32 %pi_2n, 3.14159265358979;
    mov.f32 %angle, 2.0;
    mul.f32 %angle, %angle, %w_f;
    div.full.f32 %pi_2n, %pi_2n, %angle;

    // For each output column n = 0..width
    mov.u32 %col, 0;
ROW_COL_LOOP:
    setp.ge.u32 %p, %col, %w;
    @%p bra ROW_COL_END;

    // DC term: sum = X[0] * scale * (1/sqrt(2))
    ld.global.f32 %val, [%in_base];
    mul.f32 %sum, %val, %scale;
    mul.f32 %sum, %sum, 0.7071067811865476;

    // Pre-compute (2*col + 1) as float for AC loop
    mul.lo.u32 %n2p1, %col, 2;
    add.u32 %n2p1, %n2p1, 1;
    cvt.rn.f32.u32 %n2p1_f, %n2p1;

    // AC terms: sum += X[k] * scale * cos((2n+1) * k * pi_2n)
    mov.u32 %k, 1;
ROW_AC_LOOP:
    setp.ge.u32 %p, %k, %w;
    @%p bra ROW_AC_END;

    // Load X[k] from in_base + k * 4
    mad.wide.u32 %coeff_addr, %k, 4, %in_base;
    ld.global.f32 %val, [%coeff_addr];

    // angle = (2n+1) * k * pi_2n
    cvt.rn.f32.u32 %k_f, %k;
    mul.f32 %angle, %n2p1_f, %k_f;
    mul.f32 %angle, %angle, %pi_2n;

    cos.approx.f32 %cos_val, %angle;

    // sum += X[k] * scale * cos(angle)
    mul.f32 %val, %val, %scale;
    fma.rn.f32 %sum, %val, %cos_val, %sum;

    add.u32 %k, %k, 1;
    bra ROW_AC_LOOP;

ROW_AC_END:
    // Store result at output[row * width + col]
    st.global.f32 [%out_addr], %sum;

    add.u64 %out_addr, %out_addr, 4;
    add.u32 %col, %col, 1;
    bra ROW_COL_LOOP;

ROW_COL_END:
ROW_EXIT:
    ret;
}

// 1D IDCT (Type-III DCT) on columns
// x[n] = sqrt(2/H) * [ X[0]/sqrt(2) + sum_{k=1}^{H-1} X[k] * cos(pi*(2n+1)*k / (2H)) ]
// One thread per column, inner loop over height rows.
// Column access pattern: input[k * width + col], output[n * width + col]
// Args: input, output, width, height
.visible .entry holo_spectral_idct_1d_cols(
    .param .u64 input_ptr,
    .param .u64 output_ptr,
    .param .u32 width,
    .param .u32 height
)
{
    .reg .u32 %col, %row, %k, %w, %h, %tmp, %n2p1;
    .reg .u64 %in_base, %out_base, %coeff_addr, %out_addr;
    .reg .f32 %sum, %val, %cos_val, %scale, %pi_2n, %angle, %h_f, %k_f, %n2p1_f;
    .reg .pred %p;

    // Get column index = blockIdx.x * blockDim.x + threadIdx.x
    mov.u32 %col, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %col, %col, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %col, %col, %tmp;

    ld.param.u32 %w, [width];
    setp.ge.u32 %p, %col, %w;
    @%p bra COL_EXIT;

    ld.param.u32 %h, [height];
    ld.param.u64 %in_base, [input_ptr];
    ld.param.u64 %out_base, [output_ptr];

    // Pre-compute scale = sqrt(2.0 / height)
    cvt.rn.f32.u32 %h_f, %h;
    mov.f32 %scale, 2.0;
    div.full.f32 %scale, %scale, %h_f;
    sqrt.approx.f32 %scale, %scale;

    // Pre-compute pi_2n = PI / (2 * height)
    mov.f32 %pi_2n, 3.14159265358979;
    mov.f32 %angle, 2.0;
    mul.f32 %angle, %angle, %h_f;
    div.full.f32 %pi_2n, %pi_2n, %angle;

    // For each output row n = 0..height
    mov.u32 %row, 0;
COL_ROW_LOOP:
    setp.ge.u32 %p, %row, %h;
    @%p bra COL_ROW_END;

    // DC term: sum = input[0 * width + col] * scale * (1/sqrt(2))
    mad.wide.u32 %coeff_addr, %col, 4, %in_base;
    ld.global.f32 %val, [%coeff_addr];
    mul.f32 %sum, %val, %scale;
    mul.f32 %sum, %sum, 0.7071067811865476;

    // Pre-compute (2*row + 1) as float for AC loop
    mul.lo.u32 %n2p1, %row, 2;
    add.u32 %n2p1, %n2p1, 1;
    cvt.rn.f32.u32 %n2p1_f, %n2p1;

    // AC terms: sum += input[k*width+col] * scale * cos((2n+1) * k * pi_2n)
    mov.u32 %k, 1;
COL_AC_LOOP:
    setp.ge.u32 %p, %k, %h;
    @%p bra COL_AC_END;

    // Load coefficient: input[k * width + col]
    mul.lo.u32 %tmp, %k, %w;
    add.u32 %tmp, %tmp, %col;
    mad.wide.u32 %coeff_addr, %tmp, 4, %in_base;
    ld.global.f32 %val, [%coeff_addr];

    // angle = (2*row+1) * k * pi_2n
    cvt.rn.f32.u32 %k_f, %k;
    mul.f32 %angle, %n2p1_f, %k_f;
    mul.f32 %angle, %angle, %pi_2n;

    cos.approx.f32 %cos_val, %angle;

    // sum += input[k*w+col] * scale * cos(angle)
    mul.f32 %val, %val, %scale;
    fma.rn.f32 %sum, %val, %cos_val, %sum;

    add.u32 %k, %k, 1;
    bra COL_AC_LOOP;

COL_AC_END:
    // Store result at output[row * width + col]
    mul.lo.u32 %tmp, %row, %w;
    add.u32 %tmp, %tmp, %col;
    mad.wide.u32 %out_addr, %tmp, 4, %out_base;
    st.global.f32 [%out_addr], %sum;

    add.u32 %row, %row, 1;
    bra COL_ROW_LOOP;

COL_ROW_END:
COL_EXIT:
    ret;
}

// DD-8 STUB: Fused 2D IDCT kernel.
// Currently unused - separable row + col IDCT handles 2D reconstruction.
// Potential future optimization: a single-pass 2D IDCT could reduce kernel
// launch overhead and improve cache locality for large tensors. Profile
// separable path under Nihil before implementing.
.visible .entry holo_spectral_idct_2d(
    .param .u64 input_ptr,
    .param .u64 output_ptr,
    .param .u32 width,
    .param .u32 height
)
{
    .reg .u32 %tmp;
    mov.u32 %tmp, 0;
    ret;
}
"#;

    /// RPH (random projection) kernel PTX source.
    const RPH_KERNEL_PTX: &str = r#"
.version 8.0
.target sm_89
.address_size 64

// Accumulate projection into output buffer
// Uses on-the-fly random weight generation via XORShift64 PRNG.
// Each thread handles one output element, iterating over projection dims.
// Reference: haagenti-hct SeededRng (xorshift64: <<13, >>7, <<17)
.visible .entry holo_rph_accumulate(
    .param .u64 projection_ptr,
    .param .u64 output_ptr,
    .param .u32 proj_dim,
    .param .u32 output_dim,
    .param .u64 seed
)
{
    .reg .u32 %gid, %out_idx, %proj_d, %out_d, %tmp, %i;
    .reg .u64 %proj_base, %proj_addr, %out_addr, %rng_state, %tmp64;
    .reg .f32 %sum, %proj_val, %rand_val, %out_val, %scale;
    .reg .pred %p;

    // Get output index = blockIdx.x * blockDim.x + threadIdx.x
    mov.u32 %gid, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %gid, %gid, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %out_idx, %gid, %tmp;

    ld.param.u32 %out_d, [output_dim];
    setp.ge.u32 %p, %out_idx, %out_d;
    @%p bra RPH_EXIT;

    // Save projection base pointer (Bug fix: don't modify in loop)
    ld.param.u64 %proj_base, [projection_ptr];
    ld.param.u32 %proj_d, [proj_dim];

    // Per-thread PRNG initialization:
    // state = seed + 1 (matching SeededRng::new wrapping_add(1))
    // then XOR with out_idx spread across both halves for decorrelation
    ld.param.u64 %rng_state, [seed];
    add.u64 %rng_state, %rng_state, 1;
    cvt.u64.u32 %tmp64, %out_idx;
    shl.b64 %rng_state, %rng_state, 1;
    xor.b64 %rng_state, %rng_state, %tmp64;
    shl.b64 %tmp64, %tmp64, 32;
    xor.b64 %rng_state, %rng_state, %tmp64;
    // Ensure non-zero (XORShift64 produces 0 forever from state=0)
    setp.ne.u64 %p, %rng_state, 0;
    @%p bra RPH_SEED_OK;
    mov.u64 %rng_state, 1;
RPH_SEED_OK:

    // scale = 1.0 / sqrt(proj_dim)
    cvt.rn.f32.u32 %scale, %proj_d;
    sqrt.approx.f32 %scale, %scale;
    rcp.approx.f32 %scale, %scale;

    mov.f32 %sum, 0.0;

    // Sum projection[i] * random_weight for each projection dimension
    mov.u32 %i, 0;
RPH_PROJ_LOOP:
    setp.ge.u32 %p, %i, %proj_d;
    @%p bra RPH_END_PROJ;

    // Load projection value: proj_base + i * 4
    mad.wide.u32 %proj_addr, %i, 4, %proj_base;
    ld.global.f32 %proj_val, [%proj_addr];

    // XORShift64: x ^= x << 13; x ^= x >> 7; x ^= x << 17
    shl.b64 %tmp64, %rng_state, 13;
    xor.b64 %rng_state, %rng_state, %tmp64;
    shr.u64 %tmp64, %rng_state, 7;
    xor.b64 %rng_state, %rng_state, %tmp64;
    shl.b64 %tmp64, %rng_state, 17;
    xor.b64 %rng_state, %rng_state, %tmp64;

    // Convert to float [-1, 1):
    // Take bits 40-63 (24 bits) -> [0, 2^24), divide by 2^23 -> [0, 2), sub 1 -> [-1, 1)
    shr.u64 %tmp64, %rng_state, 40;
    cvt.rn.f32.u64 %rand_val, %tmp64;
    div.full.f32 %rand_val, %rand_val, 8388608.0;
    sub.f32 %rand_val, %rand_val, 1.0;

    // Apply projection scale
    mul.f32 %rand_val, %rand_val, %scale;

    // Accumulate: sum += proj_val * rand_val
    fma.rn.f32 %sum, %proj_val, %rand_val, %sum;

    add.u32 %i, %i, 1;
    bra RPH_PROJ_LOOP;

RPH_END_PROJ:
    // Store to output: output[out_idx] += sum
    ld.param.u64 %out_addr, [output_ptr];
    mad.wide.u32 %out_addr, %out_idx, 4, %out_addr;

    ld.global.f32 %out_val, [%out_addr];
    add.f32 %out_val, %out_val, %sum;
    st.global.f32 [%out_addr], %out_val;

RPH_EXIT:
    ret;
}

// Finalize RPH by dividing by projection count
.visible .entry holo_rph_finalize(
    .param .u64 input_ptr,
    .param .u64 output_ptr,
    .param .u32 size,
    .param .u32 num_projections
)
{
    .reg .u32 %gid, %sz, %np, %tmp;
    .reg .u64 %in_addr, %out_addr;
    .reg .f32 %val, %divisor;
    .reg .pred %p;

    mov.u32 %gid, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %gid, %gid, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %gid, %gid, %tmp;

    ld.param.u32 %sz, [size];
    setp.ge.u32 %p, %gid, %sz;
    @%p bra EXIT;

    ld.param.u64 %in_addr, [input_ptr];
    ld.param.u64 %out_addr, [output_ptr];
    ld.param.u32 %np, [num_projections];

    mad.wide.u32 %in_addr, %gid, 4, %in_addr;
    mad.wide.u32 %out_addr, %gid, 4, %out_addr;

    ld.global.f32 %val, [%in_addr];
    cvt.rn.f32.u32 %divisor, %np;
    div.rn.f32 %val, %val, %divisor;
    st.global.f32 [%out_addr], %val;

EXIT:
    ret;
}

// DD-8 STUB: On-GPU projection matrix generation.
// Currently unused - projections arrive pre-computed in fragment data and
// random weights are generated on-the-fly in holo_rph_accumulate.
// Future use: pre-generating the full projection matrix on-GPU could enable
// batched RPH accumulation and avoid redundant PRNG work across fragments
// sharing the same seed. Relevant once Nihil replaces Candle and RPH
// reconstruction becomes latency-critical.
.visible .entry holo_rph_generate_projection(
    .param .u64 output_ptr,
    .param .u32 row,
    .param .u32 col,
    .param .u32 rows,
    .param .u32 cols,
    .param .u64 seed
)
{
    ret;
}
"#;

    /// LRDF (low-rank distributed) kernel PTX source.
    const LRDF_KERNEL_PTX: &str = r#"
.version 8.0
.target sm_89
.address_size 64

// Outer product: output += sigma * u * v^T
.visible .entry holo_lrdf_outer_product(
    .param .u64 u_ptr,
    .param .u64 v_ptr,
    .param .u64 output_ptr,
    .param .f32 sigma,
    .param .u32 rows,
    .param .u32 cols
)
{
    .reg .u32 %row, %col, %r, %c, %idx, %tmp;
    .reg .u64 %u_addr, %v_addr, %out_addr;
    .reg .f32 %u_val, %v_val, %prod, %out_val, %sig;
    .reg .pred %p;

    // Get 2D thread position: col = blockIdx.x * blockDim.x + threadIdx.x
    mov.u32 %col, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %col, %col, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %col, %col, %tmp;

    // row = blockIdx.y * blockDim.y + threadIdx.y
    mov.u32 %row, %ctaid.y;
    mov.u32 %tmp, %ntid.y;
    mul.lo.u32 %row, %row, %tmp;
    mov.u32 %tmp, %tid.y;
    add.u32 %row, %row, %tmp;

    // Bounds check
    ld.param.u32 %r, [rows];
    ld.param.u32 %c, [cols];
    setp.ge.u32 %p, %row, %r;
    @%p bra EXIT;
    setp.ge.u32 %p, %col, %c;
    @%p bra EXIT;

    // Load u[row] and v[col]
    ld.param.u64 %u_addr, [u_ptr];
    ld.param.u64 %v_addr, [v_ptr];

    mad.wide.u32 %u_addr, %row, 4, %u_addr;
    mad.wide.u32 %v_addr, %col, 4, %v_addr;

    ld.global.f32 %u_val, [%u_addr];
    ld.global.f32 %v_val, [%v_addr];

    // Compute sigma * u * v
    ld.param.f32 %sig, [sigma];
    mul.f32 %prod, %u_val, %v_val;
    mul.f32 %prod, %prod, %sig;

    // Add to output[row * cols + col]
    ld.param.u64 %out_addr, [output_ptr];
    mul.lo.u32 %idx, %row, %c;
    add.u32 %idx, %idx, %col;
    mad.wide.u32 %out_addr, %idx, 4, %out_addr;

    ld.global.f32 %out_val, [%out_addr];
    add.f32 %out_val, %out_val, %prod;
    st.global.f32 [%out_addr], %out_val;

EXIT:
    ret;
}

// Batched outer products: output += sum_i sigma[i] * u[i] * v[i]^T
// Memory layout: u = [u0(rows), u1(rows), ...], v = [v0(cols), v1(cols), ...], sigma = [s0, s1, ...]
// Each thread (row, col) accumulates all components in registers, single write to output.
.visible .entry holo_lrdf_outer_product_batched(
    .param .u64 u_ptr,
    .param .u64 v_ptr,
    .param .u64 sigma_ptr,
    .param .u64 output_ptr,
    .param .u32 num_components,
    .param .u32 rows,
    .param .u32 cols
)
{
    .reg .u32 %row, %col, %r, %c, %nc, %idx, %tmp, %i, %off;
    .reg .u64 %u_base, %v_base, %sig_base, %out_addr, %addr;
    .reg .f32 %u_val, %v_val, %sig_val, %prod, %out_val, %sum;
    .reg .pred %p;

    // Get 2D thread position: col = blockIdx.x * blockDim.x + threadIdx.x
    mov.u32 %col, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %col, %col, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %col, %col, %tmp;

    // row = blockIdx.y * blockDim.y + threadIdx.y
    mov.u32 %row, %ctaid.y;
    mov.u32 %tmp, %ntid.y;
    mul.lo.u32 %row, %row, %tmp;
    mov.u32 %tmp, %tid.y;
    add.u32 %row, %row, %tmp;

    // Bounds check
    ld.param.u32 %r, [rows];
    ld.param.u32 %c, [cols];
    setp.ge.u32 %p, %row, %r;
    @%p bra BATCH_EXIT;
    setp.ge.u32 %p, %col, %c;
    @%p bra BATCH_EXIT;

    // Load base pointers and component count
    ld.param.u64 %u_base, [u_ptr];
    ld.param.u64 %v_base, [v_ptr];
    ld.param.u64 %sig_base, [sigma_ptr];
    ld.param.u32 %nc, [num_components];

    // Accumulate all components in registers
    mov.f32 %sum, 0.0;
    mov.u32 %i, 0;

BATCH_COMP_LOOP:
    setp.ge.u32 %p, %i, %nc;
    @%p bra BATCH_COMP_END;

    // sigma[i]: sig_base + i * 4
    mad.wide.u32 %addr, %i, 4, %sig_base;
    ld.global.f32 %sig_val, [%addr];

    // u[i * rows + row]: u_base + (i * rows + row) * 4
    mul.lo.u32 %off, %i, %r;
    add.u32 %off, %off, %row;
    mad.wide.u32 %addr, %off, 4, %u_base;
    ld.global.f32 %u_val, [%addr];

    // v[i * cols + col]: v_base + (i * cols + col) * 4
    mul.lo.u32 %off, %i, %c;
    add.u32 %off, %off, %col;
    mad.wide.u32 %addr, %off, 4, %v_base;
    ld.global.f32 %v_val, [%addr];

    // sum += sigma[i] * u[i][row] * v[i][col]
    mul.f32 %prod, %u_val, %v_val;
    fma.rn.f32 %sum, %sig_val, %prod, %sum;

    add.u32 %i, %i, 1;
    bra BATCH_COMP_LOOP;

BATCH_COMP_END:
    // output[row * cols + col] += sum
    ld.param.u64 %out_addr, [output_ptr];
    mul.lo.u32 %idx, %row, %c;
    add.u32 %idx, %idx, %col;
    mad.wide.u32 %out_addr, %idx, 4, %out_addr;

    ld.global.f32 %out_val, [%out_addr];
    add.f32 %out_val, %out_val, %sum;
    st.global.f32 [%out_addr], %out_val;

BATCH_EXIT:
    ret;
}
"#;

    /// Fused reconstruction + dequantization kernel PTX source.
    ///
    /// Combines holographic reconstruction output with INT4/INT8 dequantization
    /// in a single pass to avoid extra memory bandwidth.
    const FUSED_KERNEL_PTX: &str = r#"
.version 8.0
.target sm_89
.address_size 64

// Fused holographic reconstruction + F32 to F16 conversion
// For use after IDCT/RPH/LRDF reconstruction when output needs to be F16
.visible .entry holo_fused_f32_to_f16(
    .param .u64 input_ptr,     // F32 reconstructed data
    .param .u64 output_ptr,    // F16 output
    .param .u32 size
)
{
    .reg .u32 %gid, %sz, %tmp;
    .reg .u64 %in_addr, %out_addr;
    .reg .f32 %val;
    .reg .b16 %h_val;
    .reg .pred %p;

    // Get thread index
    mov.u32 %gid, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %gid, %gid, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %gid, %gid, %tmp;

    ld.param.u32 %sz, [size];
    setp.ge.u32 %p, %gid, %sz;
    @%p bra EXIT;

    ld.param.u64 %in_addr, [input_ptr];
    ld.param.u64 %out_addr, [output_ptr];

    // Load F32, convert to F16, store
    mad.wide.u32 %in_addr, %gid, 4, %in_addr;
    ld.global.f32 %val, [%in_addr];

    cvt.rn.f16.f32 %h_val, %val;

    mad.wide.u32 %out_addr, %gid, 2, %out_addr;
    st.global.u16 [%out_addr], %h_val;

EXIT:
    ret;
}

// Fused: apply scale and zero-point correction to reconstructed values
// For quantized weight reconstruction
// output = (reconstructed - zero_point) * scale
.visible .entry holo_fused_dequant_f32(
    .param .u64 input_ptr,         // F32 reconstructed (may represent quantized values)
    .param .u64 scales_ptr,        // Per-block scales (F32)
    .param .u64 zeros_ptr,         // Per-block zero points (F32)
    .param .u64 output_ptr,        // F32 dequantized output
    .param .u32 size,
    .param .u32 block_size
)
{
    .reg .u32 %gid, %sz, %blk_sz, %blk_idx, %tmp;
    .reg .u64 %in_addr, %scale_addr, %zero_addr, %out_addr;
    .reg .f32 %val, %scale, %zero, %result;
    .reg .pred %p;

    mov.u32 %gid, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %gid, %gid, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %gid, %gid, %tmp;

    ld.param.u32 %sz, [size];
    setp.ge.u32 %p, %gid, %sz;
    @%p bra EXIT2;

    ld.param.u64 %in_addr, [input_ptr];
    ld.param.u64 %scale_addr, [scales_ptr];
    ld.param.u64 %zero_addr, [zeros_ptr];
    ld.param.u64 %out_addr, [output_ptr];
    ld.param.u32 %blk_sz, [block_size];

    // Calculate block index
    div.u32 %blk_idx, %gid, %blk_sz;

    // Load reconstructed value
    mad.wide.u32 %in_addr, %gid, 4, %in_addr;
    ld.global.f32 %val, [%in_addr];

    // Load scale and zero for this block
    mad.wide.u32 %scale_addr, %blk_idx, 4, %scale_addr;
    mad.wide.u32 %zero_addr, %blk_idx, 4, %zero_addr;
    ld.global.f32 %scale, [%scale_addr];
    ld.global.f32 %zero, [%zero_addr];

    // Dequantize: (val - zero) * scale
    sub.f32 %result, %val, %zero;
    mul.f32 %result, %result, %scale;

    // Store
    mad.wide.u32 %out_addr, %gid, 4, %out_addr;
    st.global.f32 [%out_addr], %result;

EXIT2:
    ret;
}

// Fused IDCT + F16 output in single kernel
// Performs 1D IDCT on a row and outputs as F16
.visible .entry holo_spectral_idct_f16(
    .param .u64 coeffs_ptr,    // F32 frequency coefficients
    .param .u64 output_ptr,    // F16 output
    .param .u32 width,
    .param .u32 height
)
{
    .reg .u32 %row, %col, %w, %h, %tmp;
    .reg .u64 %coeff_addr, %out_addr;
    .reg .f32 %sum;
    .reg .b16 %h_val;
    .reg .pred %p;

    mov.u32 %row, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %row, %row, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %row, %row, %tmp;

    ld.param.u32 %h, [height];
    setp.ge.u32 %p, %row, %h;
    @%p bra EXIT3;

    ld.param.u32 %w, [width];
    ld.param.u64 %coeff_addr, [coeffs_ptr];
    ld.param.u64 %out_addr, [output_ptr];

    // Process each column in this row
    mov.u32 %col, 0;
COL_LOOP:
    setp.ge.u32 %p, %col, %w;
    @%p bra END_COL;

    // Simplified IDCT placeholder - just copy for now
    // Real implementation would use shared memory for efficiency
    mov.f32 %sum, 0.0;

    // Convert to F16 and store
    cvt.rn.f16.f32 %h_val, %sum;

    // Calculate output offset: (row * width + col) * 2
    mul.lo.u32 %row, %row, %w;
    add.u32 %row, %row, %col;
    mad.wide.u32 %out_addr, %row, 2, %out_addr;
    st.global.u16 [%out_addr], %h_val;

    add.u32 %col, %col, 1;
    bra COL_LOOP;

END_COL:
EXIT3:
    ret;
}

// Scale reconstructed values by a constant factor
.visible .entry holo_scale_values(
    .param .u64 data_ptr,
    .param .f32 scale,
    .param .u32 size
)
{
    .reg .u32 %gid, %sz, %tmp;
    .reg .u64 %addr;
    .reg .f32 %val, %scl;
    .reg .pred %p;

    mov.u32 %gid, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %gid, %gid, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %gid, %gid, %tmp;

    ld.param.u32 %sz, [size];
    setp.ge.u32 %p, %gid, %sz;
    @%p bra EXIT4;

    ld.param.u64 %addr, [data_ptr];
    ld.param.f32 %scl, [scale];

    mad.wide.u32 %addr, %gid, 4, %addr;
    ld.global.f32 %val, [%addr];
    mul.f32 %val, %val, %scl;
    st.global.f32 [%addr], %val;

EXIT4:
    ret;
}
"#;

    impl GpuHoloContext {
        // ==================== Fused Operations ====================

        /// Loads the fused reconstruction + dequantization kernels.
        pub fn load_fused_kernel(&mut self) -> Result<(), GpuHoloError> {
            let ptx = Ptx::from_src(FUSED_KERNEL_PTX);
            self.device
                .load_ptx(
                    ptx,
                    "holo_fused",
                    &[
                        "holo_fused_f32_to_f16",
                        "holo_fused_dequant_f32",
                        "holo_spectral_idct_f16",
                        "holo_scale_values",
                    ],
                )
                .map_err(|e| GpuHoloError::KernelLoad {
                    message: e.to_string(),
                })?;

            Ok(())
        }

        /// Converts reconstructed F32 values to F16.
        ///
        /// Useful when inference expects F16 weights.
        pub fn convert_f32_to_f16(
            &self,
            input: &CudaSlice<f32>,
        ) -> Result<CudaSlice<half::f16>, GpuHoloError> {
            let size = input.len();

            let output: CudaSlice<half::f16> =
                self.device
                    .alloc_zeros(size)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            let func = self
                .device
                .get_func("holo_fused", "holo_fused_f32_to_f16")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_fused_f32_to_f16".to_string(),
                })?;

            let threads_per_block = 256u32;
            let num_blocks = ((size as u32) + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            unsafe { func.launch(cfg, (input, &output, size as u32)) }.map_err(|e| {
                GpuHoloError::KernelExec {
                    message: e.to_string(),
                }
            })?;

            self.device
                .synchronize()
                .map_err(|e| GpuHoloError::Synchronize {
                    message: e.to_string(),
                })?;

            Ok(output)
        }

        /// Applies per-block dequantization to reconstructed values.
        ///
        /// Formula: output = (input - zero_point) * scale
        ///
        /// # Arguments
        ///
        /// * `input` - Reconstructed F32 values (may represent quantized data)
        /// * `scales` - Per-block scale factors
        /// * `zeros` - Per-block zero points
        /// * `block_size` - Number of values per quantization block
        pub fn dequantize_reconstructed(
            &self,
            input: &CudaSlice<f32>,
            scales: &[f32],
            zeros: &[f32],
            block_size: usize,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            let size = input.len();

            // Copy scales and zeros to GPU
            let d_scales =
                self.device
                    .htod_copy(scales.to_vec())
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            let d_zeros =
                self.device
                    .htod_copy(zeros.to_vec())
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            let output: CudaSlice<f32> =
                self.device
                    .alloc_zeros(size)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            let func = self
                .device
                .get_func("holo_fused", "holo_fused_dequant_f32")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_fused_dequant_f32".to_string(),
                })?;

            let threads_per_block = 256u32;
            let num_blocks = ((size as u32) + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            unsafe {
                func.launch(
                    cfg,
                    (
                        input,
                        &d_scales,
                        &d_zeros,
                        &output,
                        size as u32,
                        block_size as u32,
                    ),
                )
            }
            .map_err(|e| GpuHoloError::KernelExec {
                message: e.to_string(),
            })?;

            self.device
                .synchronize()
                .map_err(|e| GpuHoloError::Synchronize {
                    message: e.to_string(),
                })?;

            Ok(output)
        }

        /// Scales all values in a buffer by a constant factor.
        pub fn scale_values(
            &self,
            data: &mut CudaSlice<f32>,
            scale: f32,
        ) -> Result<(), GpuHoloError> {
            let size = data.len();

            let func = self
                .device
                .get_func("holo_fused", "holo_scale_values")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_scale_values".to_string(),
                })?;

            let threads_per_block = 256u32;
            let num_blocks = ((size as u32) + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            unsafe { func.launch(cfg, (data as &CudaSlice<f32>, scale, size as u32)) }.map_err(
                |e| GpuHoloError::KernelExec {
                    message: e.to_string(),
                },
            )?;

            Ok(())
        }

        /// Reconstructs and dequantizes in a single pipeline.
        ///
        /// This is the high-level fused API combining:
        /// 1. Holographic reconstruction from fragments
        /// 2. Dequantization with per-block scales/zeros
        /// 3. Optional F32→F16 conversion
        pub fn reconstruct_and_dequantize(
            &self,
            header: &HoloTensorHeader,
            fragments: &[HoloFragment],
            scales: &[f32],
            zeros: &[f32],
            block_size: usize,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            // Step 1: Reconstruct
            let reconstructed = self.reconstruct(header, fragments)?;

            // Step 2: Dequantize
            let dequantized =
                self.dequantize_reconstructed(&reconstructed, scales, zeros, block_size)?;

            Ok(dequantized)
        }

        /// Reconstructs, dequantizes, and converts to F16.
        pub fn reconstruct_dequantize_f16(
            &self,
            header: &HoloTensorHeader,
            fragments: &[HoloFragment],
            scales: &[f32],
            zeros: &[f32],
            block_size: usize,
        ) -> Result<CudaSlice<half::f16>, GpuHoloError> {
            let dequantized =
                self.reconstruct_and_dequantize(header, fragments, scales, zeros, block_size)?;

            self.convert_f32_to_f16(&dequantized)
        }

        // ==================== Quality Calibration ====================

        /// Measures the actual reconstruction quality by comparing to original data.
        ///
        /// Returns a quality score between 0.0 and 1.0 based on normalized MSE.
        pub fn measure_reconstruction_quality(
            &self,
            reconstructed: &CudaSlice<f32>,
            original: &[f32],
        ) -> Result<f32, GpuHoloError> {
            // Copy reconstructed data to host
            let mut reconstructed_host = vec![0.0f32; original.len()];
            self.device
                .dtoh_sync_copy_into(reconstructed, &mut reconstructed_host)
                .map_err(|e| GpuHoloError::MemoryCopy {
                    message: e.to_string(),
                })?;

            // Calculate normalized MSE
            let mse: f32 = original
                .iter()
                .zip(reconstructed_host.iter())
                .map(|(o, r)| (o - r).powi(2))
                .sum::<f32>()
                / original.len() as f32;

            // Calculate variance of original for normalization
            let mean: f32 = original.iter().sum::<f32>() / original.len() as f32;
            let variance: f32 =
                original.iter().map(|x| (x - mean).powi(2)).sum::<f32>() / original.len() as f32;

            // Quality = 1 - normalized_mse, clamped to [0, 1]
            let nmse = if variance > 1e-10 {
                mse / variance
            } else {
                mse
            };
            let quality = (1.0 - nmse).clamp(0.0, 1.0);

            Ok(quality)
        }

        /// Calibrates quality curve by testing reconstruction at various fragment counts.
        ///
        /// Returns a new `QualityCurve` fitted to actual performance on test data.
        pub fn calibrate_quality_curve(
            &mut self,
            original: &[f32],
            width: usize,
            height: usize,
            encoding: HolographicEncoding,
            num_fragments: u16,
        ) -> Result<QualityCurve, GpuHoloError> {
            use haagenti::holotensor::HoloTensorEncoder;

            // Encode the test data
            let encoder = HoloTensorEncoder::new(encoding).with_fragments(num_fragments);
            let (header, fragments) = encoder.encode_2d(original, width, height).map_err(|e| {
                GpuHoloError::InvalidInput {
                    message: format!("Failed to encode test data: {}", e),
                }
            })?;

            // Measure quality at each fragment count
            let mut quality_points = Vec::with_capacity(num_fragments as usize);

            for k in 1..=num_fragments {
                // Reconstruct with k fragments
                let partial_frags: Vec<_> = fragments.iter().take(k as usize).cloned().collect();
                let reconstructed = self.reconstruct(&header, &partial_frags)?;
                let quality = self.measure_reconstruction_quality(&reconstructed, original)?;
                quality_points.push((k as f32 / num_fragments as f32, quality));
            }

            // Fit polynomial curve to quality points
            // Using simple least squares for [a0, a1, a2, a3]
            let coefficients = Self::fit_polynomial(&quality_points);

            // Find min and sufficient fragments
            let min_fragments = quality_points
                .iter()
                .position(|&(_, q)| q > 0.1)
                .map(|i| (i + 1) as u16)
                .unwrap_or(1);

            let sufficient_fragments = quality_points
                .iter()
                .position(|&(_, q)| q > 0.99)
                .map(|i| (i + 1) as u16)
                .unwrap_or(num_fragments);

            Ok(QualityCurve {
                coefficients,
                min_fragments,
                sufficient_fragments,
            })
        }

        /// Fits a 3rd-degree polynomial to quality data points.
        fn fit_polynomial(points: &[(f32, f32)]) -> [f32; 4] {
            // Simple polynomial regression using normal equations
            // y = a0 + a1*x + a2*x^2 + a3*x^3
            let n = points.len() as f32;
            if n < 4.0 {
                // Not enough points, return linear approximation
                return [0.0, 1.0, 0.0, 0.0];
            }

            // Build sums for normal equations
            let mut sx = [0.0f32; 7]; // x^0 to x^6
            let mut sy = [0.0f32; 4]; // y*x^0 to y*x^3

            for &(x, y) in points {
                let mut xp = 1.0f32;
                for i in 0..7 {
                    sx[i] += xp;
                    if i < 4 {
                        sy[i] += y * xp;
                    }
                    xp *= x;
                }
            }

            // Simplified solution: use least squares approximation
            // For robustness, just compute a reasonable fit
            let a0 = points.first().map(|&(_, y)| y).unwrap_or(0.0);
            let a1 = if n > 1.0 {
                points.last().map(|&(_, y)| y).unwrap_or(1.0) - a0
            } else {
                1.0
            };

            [a0, a1, 0.0, 0.0]
        }
    }

    // ==================== Streaming Holographic Context ====================

    /// Progressive holographic loader for streaming reconstruction.
    ///
    /// Supports incremental fragment feeding with quality tracking.
    pub struct ProgressiveHoloLoader {
        context: GpuHoloContext,
        header: HoloTensorHeader,
        accumulator: AccumulatorState,
        fragments_loaded: u16,
        current_quality: f32,
    }

    impl ProgressiveHoloLoader {
        /// Creates a new progressive loader.
        pub fn new(
            context: GpuHoloContext,
            header: HoloTensorHeader,
        ) -> Result<Self, GpuHoloError> {
            let accumulator = context.create_accumulator(&header)?;

            Ok(Self {
                context,
                header,
                accumulator,
                fragments_loaded: 0,
                current_quality: 0.0,
            })
        }

        /// Feeds a fragment into the loader, updating the reconstruction.
        ///
        /// Returns the new quality estimate.
        pub fn feed(&mut self, fragment: &HoloFragment) -> Result<f32, GpuHoloError> {
            self.context.accumulate_fragment(
                fragment,
                &mut self.accumulator,
                self.header.encoding,
            )?;

            self.fragments_loaded += 1;
            self.current_quality = self
                .header
                .quality_curve
                .predict(self.fragments_loaded, self.header.total_fragments);

            Ok(self.current_quality)
        }

        /// Returns the current quality estimate.
        pub fn quality(&self) -> f32 {
            self.current_quality
        }

        /// Returns the number of fragments loaded.
        pub fn fragments_loaded(&self) -> u16 {
            self.fragments_loaded
        }

        /// Checks if minimum fragments are available for reconstruction.
        pub fn can_reconstruct(&self) -> bool {
            self.fragments_loaded >= self.header.quality_curve.min_fragments
        }

        /// Checks if quality target has been reached.
        pub fn is_sufficient(&self, target: f32) -> bool {
            self.current_quality >= target
        }

        /// Finalizes reconstruction with current fragments.
        pub fn finalize(&self) -> Result<CudaSlice<f32>, GpuHoloError> {
            if !self.can_reconstruct() {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: self.header.quality_curve.min_fragments,
                    available: self.fragments_loaded,
                });
            }

            self.context
                .finalize_reconstruction(&self.accumulator, self.header.encoding)
        }

        /// Returns the header.
        pub fn header(&self) -> &HoloTensorHeader {
            &self.header
        }
    }

    // ==================== Streaming Context ====================

    use cudarc::driver::CudaStream;

    /// Stream pool for async holographic operations.
    ///
    /// Enables pipelining: while one fragment is being transferred,
    /// another is being accumulated on the GPU.
    pub struct HoloStreamPool {
        device: Arc<CudaDevice>,
        streams: Vec<CudaStream>,
        num_streams: usize,
    }

    impl HoloStreamPool {
        /// Creates a new stream pool.
        ///
        /// # Arguments
        ///
        /// * `device` - CUDA device reference
        /// * `num_streams` - Number of concurrent streams (2-4 recommended)
        pub fn new(device: Arc<CudaDevice>, num_streams: usize) -> Result<Self, GpuHoloError> {
            let mut streams = Vec::with_capacity(num_streams);

            for i in 0..num_streams {
                let stream =
                    device
                        .fork_default_stream()
                        .map_err(|e| GpuHoloError::StreamCreate {
                            stream_id: i,
                            message: e.to_string(),
                        })?;
                streams.push(stream);
            }

            Ok(Self {
                device,
                streams,
                num_streams,
            })
        }

        /// Returns a reference to stream at the given index (wraps around).
        pub fn get_stream(&self, index: usize) -> &CudaStream {
            &self.streams[index % self.num_streams]
        }

        /// Returns the number of streams in the pool.
        pub fn num_streams(&self) -> usize {
            self.num_streams
        }

        /// Synchronizes all streams in the pool.
        pub fn synchronize_all(&self) -> Result<(), GpuHoloError> {
            self.device
                .synchronize()
                .map_err(|e| GpuHoloError::Synchronize {
                    message: e.to_string(),
                })?;
            Ok(())
        }
    }

    /// Streaming holographic reconstruction context.
    ///
    /// Provides pipelined fragment loading with overlapped H2D transfers
    /// and GPU accumulation for maximum throughput.
    ///
    /// ## Pipeline Architecture
    ///
    /// ```text
    /// Fragment 0: [Transfer] [Accumulate]
    /// Fragment 1:            [Transfer  ] [Accumulate]
    /// Fragment 2:                         [Transfer  ] [Accumulate]
    /// ```
    ///
    /// ## Usage
    ///
    /// ```ignore
    /// let ctx = StreamingHoloContext::new(0, 4)?;
    /// let result = ctx.reconstruct_streaming(&header, fragments.iter(), 0.95)?;
    /// ```
    pub struct StreamingHoloContext {
        ctx: GpuHoloContext,
        stream_pool: HoloStreamPool,
        pipeline_depth: usize,
    }

    impl StreamingHoloContext {
        /// Creates a new streaming context.
        ///
        /// # Arguments
        ///
        /// * `device_id` - CUDA device ID
        /// * `pipeline_depth` - Number of concurrent operations (2-4 recommended)
        pub fn new(device_id: usize, pipeline_depth: usize) -> Result<Self, GpuHoloError> {
            let mut ctx = GpuHoloContext::new(device_id)?;

            // Load all kernels upfront with better error tracking
            ctx.load_spectral_kernel().map_err(|e| {
                eprintln!("StreamingHoloContext: spectral kernel failed: {:?}", e);
                e
            })?;
            ctx.load_rph_kernel().map_err(|e| {
                eprintln!("StreamingHoloContext: rph kernel failed: {:?}", e);
                e
            })?;
            ctx.load_lrdf_kernel().map_err(|e| {
                eprintln!("StreamingHoloContext: lrdf kernel failed: {:?}", e);
                e
            })?;
            ctx.load_fused_kernel().map_err(|e| {
                eprintln!("StreamingHoloContext: fused kernel failed: {:?}", e);
                e
            })?;

            let stream_pool = HoloStreamPool::new(Arc::clone(&ctx.device), pipeline_depth)?;

            Ok(Self {
                ctx,
                stream_pool,
                pipeline_depth,
            })
        }

        /// Returns the underlying context for synchronous operations.
        pub fn context(&self) -> &GpuHoloContext {
            &self.ctx
        }

        /// Returns mutable reference to the underlying context.
        pub fn context_mut(&mut self) -> &mut GpuHoloContext {
            &mut self.ctx
        }

        /// Returns the stream pool for advanced async operations.
        pub fn stream_pool(&self) -> &HoloStreamPool {
            &self.stream_pool
        }

        /// Reconstructs from fragments with pipelined streaming.
        ///
        /// This method overlaps H2D transfers with accumulation to maximize
        /// throughput. Fragments are processed in batches according to pipeline depth.
        ///
        /// # Arguments
        ///
        /// * `header` - HoloTensor header with encoding info
        /// * `fragments` - Iterator of fragments to process
        /// * `min_quality` - Stop early if this quality is reached (0.0 to process all)
        ///
        /// # Returns
        ///
        /// Reconstructed tensor on GPU
        pub fn reconstruct_streaming<'a, I>(
            &self,
            header: &HoloTensorHeader,
            fragments: I,
            min_quality: f32,
        ) -> Result<CudaSlice<f32>, GpuHoloError>
        where
            I: Iterator<Item = &'a HoloFragment>,
        {
            let mut accumulator = self.ctx.create_accumulator(header)?;
            let mut fragments_loaded: u16 = 0;

            // Process fragments with pipelining
            for chunk in fragments.collect::<Vec<_>>().chunks(self.pipeline_depth) {
                for (i, fragment) in chunk.iter().enumerate() {
                    let _stream = self.stream_pool.get_stream(i);

                    // Accumulate fragment (transfers data and runs kernel)
                    self.ctx
                        .accumulate_fragment(fragment, &mut accumulator, header.encoding)?;
                    fragments_loaded += 1;
                }

                // Synchronize after each batch
                self.stream_pool.synchronize_all()?;

                // Check if quality target reached
                if min_quality > 0.0 {
                    let quality = header
                        .quality_curve
                        .predict(fragments_loaded, header.total_fragments);
                    if quality >= min_quality {
                        break;
                    }
                }
            }

            // Check minimum fragment requirement
            if fragments_loaded < header.quality_curve.min_fragments {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: header.quality_curve.min_fragments,
                    available: fragments_loaded,
                });
            }

            // Finalize reconstruction
            self.ctx
                .finalize_reconstruction(&accumulator, header.encoding)
        }

        /// Reconstructs with a callback after each fragment.
        ///
        /// Useful for progress reporting or early termination based on quality.
        ///
        /// # Arguments
        ///
        /// * `header` - HoloTensor header
        /// * `fragments` - Iterator of fragments
        /// * `callback` - Called after each fragment with (fragments_loaded, quality)
        ///                Return `false` to stop early.
        pub fn reconstruct_with_callback<'a, I, F>(
            &self,
            header: &HoloTensorHeader,
            fragments: I,
            mut callback: F,
        ) -> Result<CudaSlice<f32>, GpuHoloError>
        where
            I: Iterator<Item = &'a HoloFragment>,
            F: FnMut(u16, f32) -> bool,
        {
            let mut accumulator = self.ctx.create_accumulator(header)?;
            let mut fragments_loaded: u16 = 0;

            for fragment in fragments {
                self.ctx
                    .accumulate_fragment(fragment, &mut accumulator, header.encoding)?;
                fragments_loaded += 1;

                let quality = header
                    .quality_curve
                    .predict(fragments_loaded, header.total_fragments);

                // Call user callback - return false to stop
                if !callback(fragments_loaded, quality) {
                    break;
                }
            }

            // Check minimum fragment requirement
            if fragments_loaded < header.quality_curve.min_fragments {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: header.quality_curve.min_fragments,
                    available: fragments_loaded,
                });
            }

            self.ctx
                .finalize_reconstruction(&accumulator, header.encoding)
        }

        /// Reconstructs, dequantizes, and converts to F16 in a streaming pipeline.
        ///
        /// This is the high-level API for model loading that combines:
        /// 1. Streaming holographic reconstruction
        /// 2. Dequantization with per-block scales/zeros
        /// 3. F32 to F16 conversion
        pub fn reconstruct_dequantize_f16_streaming<'a, I>(
            &self,
            header: &HoloTensorHeader,
            fragments: I,
            scales: &[f32],
            zeros: &[f32],
            block_size: usize,
            min_quality: f32,
        ) -> Result<CudaSlice<half::f16>, GpuHoloError>
        where
            I: Iterator<Item = &'a HoloFragment>,
        {
            // Step 1: Reconstruct with streaming
            let reconstructed = self.reconstruct_streaming(header, fragments, min_quality)?;

            // Step 2: Dequantize
            let dequantized =
                self.ctx
                    .dequantize_reconstructed(&reconstructed, scales, zeros, block_size)?;

            // Step 3: Convert to F16
            self.ctx.convert_f32_to_f16(&dequantized)
        }

        /// Returns streaming statistics for the last operation.
        pub fn stats(&self) -> StreamingHoloStats {
            StreamingHoloStats {
                pipeline_depth: self.pipeline_depth,
                num_streams: self.stream_pool.num_streams(),
            }
        }
    }

    /// Statistics for streaming operations.
    #[derive(Debug, Clone)]
    pub struct StreamingHoloStats {
        /// Pipeline depth used
        pub pipeline_depth: usize,
        /// Number of CUDA streams
        pub num_streams: usize,
    }

    // ==================== Memory Coalescing Optimization ====================

    /// Optimized kernel PTX with memory coalescing.
    ///
    /// These kernels use vectorized loads/stores (ld.global.v4, st.global.v4)
    /// for coalesced memory access when accessing consecutive elements.
    const COALESCED_KERNEL_PTX: &str = r#"
.version 8.0
.target sm_89
.address_size 64

// Coalesced vectorized accumulate (4 elements per thread)
// Processes 4 consecutive elements per thread for optimal memory coalescing
.visible .entry holo_coalesced_accumulate_v4(
    .param .u64 src_ptr,
    .param .u64 dst_ptr,
    .param .u32 num_elements
)
{
    .reg .u32 %gid, %base_idx, %sz, %tmp;
    .reg .u64 %src_addr, %dst_addr;
    .reg .v4 .f32 %src_vec, %dst_vec, %result;
    .reg .pred %p;

    // Get thread ID and compute base index (4 elements per thread)
    mov.u32 %gid, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %gid, %gid, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %gid, %gid, %tmp;
    mul.lo.u32 %base_idx, %gid, 4;

    ld.param.u32 %sz, [num_elements];
    setp.ge.u32 %p, %base_idx, %sz;
    @%p bra EXIT;

    ld.param.u64 %src_addr, [src_ptr];
    ld.param.u64 %dst_addr, [dst_ptr];

    // Compute aligned addresses
    mad.wide.u32 %src_addr, %base_idx, 4, %src_addr;
    mad.wide.u32 %dst_addr, %base_idx, 4, %dst_addr;

    // Vectorized load (4 floats = 16 bytes, naturally coalesced)
    ld.global.v4.f32 %src_vec, [%src_addr];
    ld.global.v4.f32 %dst_vec, [%dst_addr];

    // Add vectors
    add.f32 %result.x, %src_vec.x, %dst_vec.x;
    add.f32 %result.y, %src_vec.y, %dst_vec.y;
    add.f32 %result.z, %src_vec.z, %dst_vec.z;
    add.f32 %result.w, %src_vec.w, %dst_vec.w;

    // Vectorized store
    st.global.v4.f32 [%dst_addr], %result;

EXIT:
    ret;
}

// Coalesced IDCT with shared memory (tile-based)
// Each block processes a tile using shared memory for coefficient reuse
.visible .entry holo_coalesced_idct_tile(
    .param .u64 coeffs_ptr,
    .param .u64 output_ptr,
    .param .u32 width,
    .param .u32 height,
    .param .u32 tile_size
)
{
    .shared .align 16 .f32 tile[1024];  // 32x32 tile max
    .reg .u32 %tid_x, %tid_y, %tile_x, %tile_y, %w, %h, %ts;
    .reg .u64 %coeff_addr, %out_addr;
    .reg .f32 %val, %cos_val, %sum;
    .reg .pred %p;

    // Get 2D thread position within block
    mov.u32 %tid_x, %tid.x;
    mov.u32 %tid_y, %tid.y;

    // Get tile position
    mov.u32 %tile_x, %ctaid.x;
    mov.u32 %tile_y, %ctaid.y;

    ld.param.u32 %w, [width];
    ld.param.u32 %h, [height];
    ld.param.u32 %ts, [tile_size];

    // Load coefficients into shared memory (coalesced)
    // Each thread loads one element
    mul.lo.u32 %tile_x, %tile_x, %ts;
    mul.lo.u32 %tile_y, %tile_y, %ts;
    add.u32 %tile_x, %tile_x, %tid_x;
    add.u32 %tile_y, %tile_y, %tid_y;

    setp.ge.u32 %p, %tile_x, %w;
    @%p bra SKIP_LOAD;
    setp.ge.u32 %p, %tile_y, %h;
    @%p bra SKIP_LOAD;

    // Compute global index (row-major)
    mul.lo.u32 %tile_y, %tile_y, %w;
    add.u32 %tile_y, %tile_y, %tile_x;

    ld.param.u64 %coeff_addr, [coeffs_ptr];
    mad.wide.u32 %coeff_addr, %tile_y, 4, %coeff_addr;
    ld.global.f32 %val, [%coeff_addr];

    // Store to shared memory
    mul.lo.u32 %tile_y, %tid_y, %ts;
    add.u32 %tile_y, %tile_y, %tid_x;
    mul.lo.u32 %tile_y, %tile_y, 4;
    mov.u32 %tile_x, tile;
    add.u32 %tile_x, %tile_x, %tile_y;
    st.shared.f32 [%tile_x], %val;

SKIP_LOAD:
    bar.sync 0;

    // IDCT computation using shared memory
    // (Simplified: actual IDCT would iterate over frequencies)
    mov.f32 %sum, 0.0;

    // Store result
    ld.param.u64 %out_addr, [output_ptr];
    st.global.f32 [%out_addr], %sum;

EXIT:
    ret;
}

// Coalesced F32 to F16 conversion (4 elements per thread)
.visible .entry holo_coalesced_f32_to_f16_v4(
    .param .u64 input_ptr,
    .param .u64 output_ptr,
    .param .u32 size
)
{
    .reg .u32 %gid, %base_idx, %sz, %tmp;
    .reg .u64 %in_addr, %out_addr;
    .reg .v4 .f32 %f32_vec;
    .reg .v4 .b16 %f16_vec;
    .reg .pred %p;

    mov.u32 %gid, %ctaid.x;
    mov.u32 %tmp, %ntid.x;
    mul.lo.u32 %gid, %gid, %tmp;
    mov.u32 %tmp, %tid.x;
    add.u32 %gid, %gid, %tmp;
    mul.lo.u32 %base_idx, %gid, 4;

    ld.param.u32 %sz, [size];
    setp.ge.u32 %p, %base_idx, %sz;
    @%p bra EXIT;

    ld.param.u64 %in_addr, [input_ptr];
    ld.param.u64 %out_addr, [output_ptr];

    mad.wide.u32 %in_addr, %base_idx, 4, %in_addr;
    mad.wide.u32 %out_addr, %base_idx, 2, %out_addr;

    // Vectorized load F32
    ld.global.v4.f32 %f32_vec, [%in_addr];

    // Convert to F16
    cvt.rn.f16.f32 %f16_vec.x, %f32_vec.x;
    cvt.rn.f16.f32 %f16_vec.y, %f32_vec.y;
    cvt.rn.f16.f32 %f16_vec.z, %f32_vec.z;
    cvt.rn.f16.f32 %f16_vec.w, %f32_vec.w;

    // Vectorized store F16
    st.global.v4.u16 [%out_addr], %f16_vec;

EXIT:
    ret;
}
"#;

    /// True CUDA pinned (page-locked) memory buffer.
    ///
    /// Uses cuMemAllocHost_v2 for proper page-locked allocation that provides:
    /// - ~2x higher bandwidth for H2D/D2H transfers
    /// - Ability to overlap transfers with kernel execution
    /// - Required for async memory operations (cudaMemcpyAsync)
    ///
    /// # Safety
    /// The buffer is automatically freed when dropped via cuMemFreeHost.
    pub struct PinnedBuffer {
        ptr: *mut u8,
        size: usize,
        capacity: usize,
    }

    // Safety: PinnedBuffer is just a raw pointer to host memory, safe to send across threads
    unsafe impl Send for PinnedBuffer {}
    unsafe impl Sync for PinnedBuffer {}

    impl PinnedBuffer {
        /// Allocates a new pinned memory buffer of the given size.
        ///
        /// Returns None if CUDA allocation fails.
        pub fn new(size: usize) -> Option<Self> {
            if size == 0 {
                return Some(Self {
                    ptr: std::ptr::NonNull::dangling().as_ptr(),
                    size: 0,
                    capacity: 0,
                });
            }

            // Use cudarc's driver sys module to call cuMemAllocHost_v2
            use std::ffi::c_void;

            let mut host_ptr: *mut c_void = std::ptr::null_mut();
            // Safety: lib() returns the loaded CUDA driver library
            let result =
                unsafe { cudarc::driver::sys::lib().cuMemAllocHost_v2(&mut host_ptr, size) };

            if result != cudarc::driver::sys::CUresult::CUDA_SUCCESS {
                // Fallback: try regular allocation if CUDA pinned fails
                // This can happen if CUDA context isn't active
                return None;
            }

            Some(Self {
                ptr: host_ptr as *mut u8,
                size,
                capacity: size,
            })
        }

        /// Returns a mutable slice to the buffer contents.
        #[inline]
        pub fn as_mut_slice(&mut self) -> &mut [u8] {
            if self.capacity == 0 {
                return &mut [];
            }
            unsafe { std::slice::from_raw_parts_mut(self.ptr, self.size) }
        }

        /// Returns an immutable slice to the buffer contents.
        #[inline]
        pub fn as_slice(&self) -> &[u8] {
            if self.capacity == 0 {
                return &[];
            }
            unsafe { std::slice::from_raw_parts(self.ptr, self.size) }
        }

        /// Returns the raw pointer to the buffer.
        #[inline]
        pub fn as_ptr(&self) -> *const u8 {
            self.ptr
        }

        /// Returns the mutable raw pointer to the buffer.
        #[inline]
        pub fn as_mut_ptr(&mut self) -> *mut u8 {
            self.ptr
        }

        /// Returns the current size (used portion) of the buffer.
        #[inline]
        pub fn len(&self) -> usize {
            self.size
        }

        /// Returns true if the buffer is empty.
        #[inline]
        pub fn is_empty(&self) -> bool {
            self.size == 0
        }

        /// Returns the total capacity of the buffer.
        #[inline]
        pub fn capacity(&self) -> usize {
            self.capacity
        }

        /// Sets the used size of the buffer (must be <= capacity).
        #[inline]
        pub fn set_len(&mut self, len: usize) {
            debug_assert!(len <= self.capacity);
            self.size = len.min(self.capacity);
        }

        /// Copies data from a slice into the buffer.
        ///
        /// Panics if the slice is larger than the buffer capacity.
        pub fn copy_from_slice(&mut self, src: &[u8]) {
            assert!(
                src.len() <= self.capacity,
                "Source slice too large for pinned buffer"
            );
            unsafe {
                std::ptr::copy_nonoverlapping(src.as_ptr(), self.ptr, src.len());
            }
            self.size = src.len();
        }
    }

    impl Drop for PinnedBuffer {
        fn drop(&mut self) {
            if self.capacity > 0 && !self.ptr.is_null() {
                use std::ffi::c_void;

                // Safety: lib() returns the loaded CUDA driver library
                unsafe {
                    let _ = cudarc::driver::sys::lib().cuMemFreeHost(self.ptr as *mut c_void);
                }
            }
        }
    }

    /// Pinned (page-locked) memory pool for fast H2D transfers.
    ///
    /// Uses true CUDA page-locked memory via cuMemAllocHost for maximum
    /// H2D/D2H transfer bandwidth (up to 2x faster than pageable memory).
    ///
    /// Falls back to regular Vec allocation if CUDA pinned allocation fails.
    pub struct PinnedMemoryPool {
        #[allow(dead_code)]
        device: Arc<CudaDevice>,
        /// Pre-allocated pinned buffers by size class
        pinned_pools: std::collections::HashMap<usize, Vec<PinnedBuffer>>,
        /// Fallback regular buffers (when pinned allocation fails)
        fallback_pools: std::collections::HashMap<usize, Vec<Vec<u8>>>,
        /// Size classes (powers of 2 from 4KB to 64MB)
        size_classes: Vec<usize>,
        /// Whether we successfully use pinned memory
        pinned_available: bool,
        /// Stats
        pinned_allocations: std::sync::atomic::AtomicUsize,
        fallback_allocations: std::sync::atomic::AtomicUsize,
    }

    /// Handle to a pinned or fallback buffer from the pool.
    pub enum PooledBuffer {
        /// True CUDA pinned buffer
        Pinned(PinnedBuffer),
        /// Fallback regular Vec buffer
        Fallback(Vec<u8>),
    }

    impl PooledBuffer {
        /// Returns a mutable slice to the buffer.
        pub fn as_mut_slice(&mut self) -> &mut [u8] {
            match self {
                PooledBuffer::Pinned(buf) => buf.as_mut_slice(),
                PooledBuffer::Fallback(buf) => buf.as_mut_slice(),
            }
        }

        /// Returns an immutable slice to the buffer.
        pub fn as_slice(&self) -> &[u8] {
            match self {
                PooledBuffer::Pinned(buf) => buf.as_slice(),
                PooledBuffer::Fallback(buf) => buf.as_slice(),
            }
        }

        /// Returns the capacity of the buffer.
        pub fn capacity(&self) -> usize {
            match self {
                PooledBuffer::Pinned(buf) => buf.capacity(),
                PooledBuffer::Fallback(buf) => buf.capacity(),
            }
        }

        /// Returns true if this is a true pinned buffer.
        pub fn is_pinned(&self) -> bool {
            matches!(self, PooledBuffer::Pinned(_))
        }

        /// Copies data from a slice into the buffer.
        pub fn copy_from_slice(&mut self, src: &[u8]) {
            match self {
                PooledBuffer::Pinned(buf) => buf.copy_from_slice(src),
                PooledBuffer::Fallback(buf) => {
                    buf.clear();
                    buf.extend_from_slice(src);
                },
            }
        }
    }

    impl PinnedMemoryPool {
        /// Creates a new pinned memory pool.
        pub fn new(device: Arc<CudaDevice>) -> Self {
            // Size classes: 4KB, 16KB, 64KB, 256KB, 1MB, 4MB, 16MB, 64MB
            let size_classes: Vec<usize> = (12..=26).step_by(2).map(|exp| 1usize << exp).collect();

            // Test if pinned allocation works
            let pinned_available = PinnedBuffer::new(4096).is_some();

            Self {
                device,
                pinned_pools: std::collections::HashMap::new(),
                fallback_pools: std::collections::HashMap::new(),
                size_classes,
                pinned_available,
                pinned_allocations: std::sync::atomic::AtomicUsize::new(0),
                fallback_allocations: std::sync::atomic::AtomicUsize::new(0),
            }
        }

        /// Returns true if true CUDA pinned memory is available.
        pub fn is_pinned_available(&self) -> bool {
            self.pinned_available
        }

        /// Gets the size class for a given allocation size.
        fn size_class(&self, size: usize) -> usize {
            for &class in &self.size_classes {
                if size <= class {
                    return class;
                }
            }
            // For very large allocations, round up to next power of 2
            size.next_power_of_two()
        }

        /// Allocates a buffer of at least the given size.
        ///
        /// Returns a pinned buffer if available, otherwise falls back to regular Vec.
        pub fn allocate(&mut self, size: usize) -> PooledBuffer {
            use std::sync::atomic::Ordering;
            let class = self.size_class(size);

            // Try pinned pool first
            if self.pinned_available {
                if let Some(pool) = self.pinned_pools.get_mut(&class) {
                    if let Some(mut buf) = pool.pop() {
                        buf.set_len(size);
                        self.pinned_allocations.fetch_add(1, Ordering::Relaxed);
                        return PooledBuffer::Pinned(buf);
                    }
                }

                // Try to allocate new pinned buffer
                if let Some(mut buf) = PinnedBuffer::new(class) {
                    buf.set_len(size);
                    self.pinned_allocations.fetch_add(1, Ordering::Relaxed);
                    return PooledBuffer::Pinned(buf);
                }
            }

            // Fallback to regular pool
            if let Some(pool) = self.fallback_pools.get_mut(&class) {
                if let Some(mut buf) = pool.pop() {
                    buf.resize(size, 0);
                    self.fallback_allocations.fetch_add(1, Ordering::Relaxed);
                    return PooledBuffer::Fallback(buf);
                }
            }

            // Allocate new fallback buffer
            self.fallback_allocations.fetch_add(1, Ordering::Relaxed);
            PooledBuffer::Fallback(vec![0u8; class])
        }

        /// Allocates a buffer of exactly the given size (for legacy Vec<u8> compatibility).
        pub fn allocate_vec(&mut self, size: usize) -> Vec<u8> {
            let class = self.size_class(size);

            // For Vec API, only use fallback pool
            if let Some(pool) = self.fallback_pools.get_mut(&class) {
                if let Some(mut buf) = pool.pop() {
                    buf.resize(size, 0);
                    return buf;
                }
            }

            vec![0u8; class]
        }

        /// Returns a buffer to the pool for reuse.
        pub fn deallocate(&mut self, buf: PooledBuffer) {
            match buf {
                PooledBuffer::Pinned(mut pinned_buf) => {
                    let class = pinned_buf.capacity();
                    pinned_buf.set_len(class); // Reset to full capacity
                    self.pinned_pools
                        .entry(class)
                        .or_insert_with(Vec::new)
                        .push(pinned_buf);
                },
                PooledBuffer::Fallback(buf) => {
                    let class = buf.capacity();
                    self.fallback_pools
                        .entry(class)
                        .or_insert_with(Vec::new)
                        .push(buf);
                },
            }
        }

        /// Returns a Vec buffer to the pool (legacy API).
        pub fn deallocate_vec(&mut self, buf: Vec<u8>) {
            let class = buf.capacity();
            self.fallback_pools
                .entry(class)
                .or_insert_with(Vec::new)
                .push(buf);
        }

        /// Pre-warms the pool by allocating buffers for each size class.
        pub fn prewarm(&mut self, buffers_per_class: usize) {
            for &class in &self.size_classes.clone() {
                for _ in 0..buffers_per_class {
                    if self.pinned_available {
                        if let Some(buf) = PinnedBuffer::new(class) {
                            self.pinned_pools
                                .entry(class)
                                .or_insert_with(Vec::new)
                                .push(buf);
                            continue;
                        }
                    }
                    // Fallback
                    let buf = vec![0u8; class];
                    self.fallback_pools
                        .entry(class)
                        .or_insert_with(Vec::new)
                        .push(buf);
                }
            }
        }

        /// Returns statistics about pool usage.
        pub fn stats(&self) -> PinnedPoolStats {
            use std::sync::atomic::Ordering;

            let mut pinned_buffers = 0;
            let mut pinned_bytes = 0;
            let mut fallback_buffers = 0;
            let mut fallback_bytes = 0;

            for (&class, pool) in &self.pinned_pools {
                pinned_buffers += pool.len();
                pinned_bytes += class * pool.len();
            }

            for (&class, pool) in &self.fallback_pools {
                fallback_buffers += pool.len();
                fallback_bytes += class * pool.len();
            }

            PinnedPoolStats {
                num_size_classes: self.size_classes.len(),
                pinned_buffers,
                pinned_bytes,
                fallback_buffers,
                fallback_bytes,
                pinned_available: self.pinned_available,
                total_pinned_allocations: self.pinned_allocations.load(Ordering::Relaxed),
                total_fallback_allocations: self.fallback_allocations.load(Ordering::Relaxed),
            }
        }
    }

    /// Statistics for pinned memory pool.
    #[derive(Debug, Clone)]
    pub struct PinnedPoolStats {
        /// Number of size classes
        pub num_size_classes: usize,
        /// Pinned buffers currently in pool
        pub pinned_buffers: usize,
        /// Pinned bytes currently in pool
        pub pinned_bytes: usize,
        /// Fallback buffers currently in pool
        pub fallback_buffers: usize,
        /// Fallback bytes currently in pool
        pub fallback_bytes: usize,
        /// Whether true CUDA pinned memory is available
        pub pinned_available: bool,
        /// Total pinned allocations served
        pub total_pinned_allocations: usize,
        /// Total fallback allocations served
        pub total_fallback_allocations: usize,
    }

    impl PinnedPoolStats {
        /// Total buffers in pool (pinned + fallback)
        pub fn total_buffers(&self) -> usize {
            self.pinned_buffers + self.fallback_buffers
        }

        /// Total bytes in pool (pinned + fallback)
        pub fn total_bytes(&self) -> usize {
            self.pinned_bytes + self.fallback_bytes
        }
    }

    // ==================== Multi-GPU Support ====================

    /// Multi-GPU holographic reconstruction context.
    ///
    /// Distributes fragments across multiple GPUs for parallel processing,
    /// then gathers results back to a primary device.
    ///
    /// ## Distribution Strategy
    ///
    /// Fragments are distributed round-robin across devices:
    /// - Device 0: fragments 0, N, 2N, ...
    /// - Device 1: fragments 1, N+1, 2N+1, ...
    /// - ...
    ///
    /// Each device maintains its own accumulator, and results are combined
    /// on the primary device during finalization.
    pub struct MultiGpuHoloContext {
        /// Primary device (where final result is assembled)
        primary_device_id: usize,
        /// Per-device contexts
        contexts: Vec<GpuHoloContext>,
        /// Per-device stream pools for pipelining
        stream_pools: Vec<HoloStreamPool>,
        /// Number of devices
        num_devices: usize,
    }

    impl MultiGpuHoloContext {
        /// Creates a multi-GPU context using specified devices.
        ///
        /// # Arguments
        ///
        /// * `device_ids` - List of CUDA device IDs to use
        /// * `streams_per_device` - Number of streams per device for pipelining
        pub fn new(device_ids: &[usize], streams_per_device: usize) -> Result<Self, GpuHoloError> {
            if device_ids.is_empty() {
                return Err(GpuHoloError::InvalidInput {
                    message: "At least one device required".to_string(),
                });
            }

            let primary_device_id = device_ids[0];
            let mut contexts = Vec::with_capacity(device_ids.len());
            let mut stream_pools = Vec::with_capacity(device_ids.len());

            for &device_id in device_ids {
                let mut ctx = GpuHoloContext::new(device_id)?;
                ctx.load_all_kernels()?;
                ctx.load_fused_kernel()?;

                let pool = HoloStreamPool::new(Arc::clone(ctx.device()), streams_per_device)?;

                contexts.push(ctx);
                stream_pools.push(pool);
            }

            Ok(Self {
                primary_device_id,
                contexts,
                stream_pools,
                num_devices: device_ids.len(),
            })
        }

        /// Creates a context using all available CUDA devices.
        pub fn new_all_devices(streams_per_device: usize) -> Result<Self, GpuHoloError> {
            // Query available devices
            let num_devices = CudaDevice::count().map_err(|e| GpuHoloError::DeviceInit {
                device_id: 0,
                message: format!("Failed to count devices: {}", e),
            })?;

            if num_devices == 0 {
                return Err(GpuHoloError::DeviceInit {
                    device_id: 0,
                    message: "No CUDA devices available".to_string(),
                });
            }

            let device_ids: Vec<usize> = (0..num_devices as usize).collect();
            Self::new(&device_ids, streams_per_device)
        }

        /// Returns the number of devices in use.
        pub fn num_devices(&self) -> usize {
            self.num_devices
        }

        /// Returns the primary device ID.
        pub fn primary_device(&self) -> usize {
            self.primary_device_id
        }

        /// Returns context for a specific device.
        pub fn device_context(&self, device_idx: usize) -> Option<&GpuHoloContext> {
            self.contexts.get(device_idx)
        }

        /// Reconstructs from fragments using all devices in parallel.
        ///
        /// Fragments are distributed across devices, each device accumulates
        /// its portion, and results are combined on the primary device.
        pub fn reconstruct_multi_gpu(
            &self,
            header: &HoloTensorHeader,
            fragments: &[HoloFragment],
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            if fragments.is_empty() {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: header.quality_curve.min_fragments,
                    available: 0,
                });
            }

            // Create accumulators on each device
            let mut accumulators: Vec<AccumulatorState> = self
                .contexts
                .iter()
                .map(|ctx| ctx.create_accumulator(header))
                .collect::<Result<Vec<_>, _>>()?;

            // Distribute fragments round-robin across devices
            for (i, fragment) in fragments.iter().enumerate() {
                let device_idx = i % self.num_devices;
                self.contexts[device_idx].accumulate_fragment(
                    fragment,
                    &mut accumulators[device_idx],
                    header.encoding,
                )?;
            }

            // Finalize reconstruction on each device
            let mut device_results: Vec<CudaSlice<f32>> = Vec::with_capacity(self.num_devices);
            for (ctx, accumulator) in self.contexts.iter().zip(accumulators.iter()) {
                let result = ctx.finalize_reconstruction(accumulator, header.encoding)?;
                device_results.push(result);
            }

            // Combine results on primary device
            // For additive accumulation schemes (RPH, LRDF), we average
            // For spectral, each device may have different coefficients
            self.combine_results(&device_results, header)
        }

        /// Combines results from multiple devices onto the primary device.
        fn combine_results(
            &self,
            device_results: &[CudaSlice<f32>],
            header: &HoloTensorHeader,
        ) -> Result<CudaSlice<f32>, GpuHoloError> {
            if device_results.is_empty() {
                return Err(GpuHoloError::InvalidInput {
                    message: "No device results to combine".to_string(),
                });
            }

            // For single device, just return the result
            if device_results.len() == 1 {
                return Ok(device_results[0].clone());
            }

            let primary_ctx = &self.contexts[0];
            let size = device_results[0].len();

            // Allocate output on primary device
            let output: CudaSlice<f32> =
                primary_ctx
                    .device
                    .alloc_zeros(size)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            // Copy first result to output
            let mut host_buf = vec![0.0f32; size];
            primary_ctx
                .device
                .dtoh_sync_copy_into(&device_results[0], &mut host_buf)
                .map_err(|e| GpuHoloError::MemoryCopy {
                    message: e.to_string(),
                })?;

            // Accumulate other results
            for result in device_results.iter().skip(1) {
                let mut other_buf = vec![0.0f32; size];
                // Note: In production, we'd use P2P transfers between devices
                // For now, go through host memory
                primary_ctx
                    .device
                    .dtoh_sync_copy_into(result, &mut other_buf)
                    .map_err(|e| GpuHoloError::MemoryCopy {
                        message: e.to_string(),
                    })?;

                for (h, o) in host_buf.iter_mut().zip(other_buf.iter()) {
                    *h += *o;
                }
            }

            // For averaging (RPH), normalize by device count
            if header.encoding == HolographicEncoding::RandomProjection {
                let scale = 1.0 / self.num_devices as f32;
                for h in host_buf.iter_mut() {
                    *h *= scale;
                }
            }

            // Copy back to device
            primary_ctx
                .device
                .htod_sync_copy_into(&host_buf, &mut output.clone())
                .map_err(|e| GpuHoloError::MemoryCopy {
                    message: e.to_string(),
                })?;

            Ok(output)
        }

        /// Returns statistics about multi-GPU configuration.
        pub fn stats(&self) -> MultiGpuStats {
            MultiGpuStats {
                num_devices: self.num_devices,
                primary_device: self.primary_device_id,
                streams_per_device: self
                    .stream_pools
                    .first()
                    .map(|p| p.num_streams())
                    .unwrap_or(0),
            }
        }
    }

    /// Statistics for multi-GPU context.
    #[derive(Debug, Clone)]
    pub struct MultiGpuStats {
        /// Number of devices in use
        pub num_devices: usize,
        /// Primary device ID
        pub primary_device: usize,
        /// Streams per device
        pub streams_per_device: usize,
    }

    // ==================== Coalesced Operations ====================

    impl GpuHoloContext {
        /// Loads coalesced (vectorized) kernels for optimized memory access.
        pub fn load_coalesced_kernels(&mut self) -> Result<(), GpuHoloError> {
            let ptx = Ptx::from_src(COALESCED_KERNEL_PTX);
            self.device
                .load_ptx(
                    ptx,
                    "holo_coalesced",
                    &[
                        "holo_coalesced_accumulate_v4",
                        "holo_coalesced_idct_tile",
                        "holo_coalesced_f32_to_f16_v4",
                    ],
                )
                .map_err(|e| GpuHoloError::KernelLoad {
                    message: e.to_string(),
                })?;

            Ok(())
        }

        /// Performs coalesced vectorized accumulation (4 elements per thread).
        ///
        /// This is optimized for large contiguous memory regions where
        /// adjacent threads access adjacent memory locations.
        pub fn accumulate_coalesced_v4(
            &self,
            src: &CudaSlice<f32>,
            dst: &mut CudaSlice<f32>,
        ) -> Result<(), GpuHoloError> {
            let size = src.len().min(dst.len());
            if size == 0 {
                return Ok(());
            }

            let func = self
                .device
                .get_func("holo_coalesced", "holo_coalesced_accumulate_v4")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_coalesced_accumulate_v4".to_string(),
                })?;

            // 4 elements per thread
            let elements_per_thread = 4u32;
            let threads_per_block = 256u32;
            let num_threads = ((size as u32) + elements_per_thread - 1) / elements_per_thread;
            let num_blocks = (num_threads + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            unsafe { func.launch(cfg, (src, dst as &CudaSlice<f32>, size as u32)) }.map_err(
                |e| GpuHoloError::KernelExec {
                    message: e.to_string(),
                },
            )?;

            Ok(())
        }

        /// Converts F32 to F16 using coalesced vectorized operations.
        pub fn convert_f32_to_f16_coalesced(
            &self,
            input: &CudaSlice<f32>,
        ) -> Result<CudaSlice<half::f16>, GpuHoloError> {
            let size = input.len();
            if size == 0 {
                return Err(GpuHoloError::InvalidInput {
                    message: "Empty input buffer".to_string(),
                });
            }

            let output: CudaSlice<half::f16> =
                self.device
                    .alloc_zeros(size)
                    .map_err(|e| GpuHoloError::MemoryAlloc {
                        message: e.to_string(),
                    })?;

            let func = self
                .device
                .get_func("holo_coalesced", "holo_coalesced_f32_to_f16_v4")
                .ok_or_else(|| GpuHoloError::KernelNotLoaded {
                    kernel: "holo_coalesced_f32_to_f16_v4".to_string(),
                })?;

            let elements_per_thread = 4u32;
            let threads_per_block = 256u32;
            let num_threads = ((size as u32) + elements_per_thread - 1) / elements_per_thread;
            let num_blocks = (num_threads + threads_per_block - 1) / threads_per_block;

            let cfg = LaunchConfig {
                grid_dim: (num_blocks.max(1), 1, 1),
                block_dim: (threads_per_block, 1, 1),
                shared_mem_bytes: 0,
            };

            unsafe { func.launch(cfg, (input, &output, size as u32)) }.map_err(|e| {
                GpuHoloError::KernelExec {
                    message: e.to_string(),
                }
            })?;

            self.device
                .synchronize()
                .map_err(|e| GpuHoloError::Synchronize {
                    message: e.to_string(),
                })?;

            Ok(output)
        }
    }

    // ==================== Phase 7: Fault Tolerance ====================

    /// Result of checksum validation.
    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
    pub enum ValidationResult {
        /// Fragment is valid.
        Valid,
        /// Checksum mismatch - fragment is corrupted.
        Corrupted,
        /// Fragment is missing (not yet received).
        Missing,
    }

    /// Fragment validation and recovery options.
    #[derive(Debug, Clone)]
    pub struct FaultToleranceConfig {
        /// Whether to validate checksums on fragment load.
        pub validate_checksums: bool,
        /// Whether to skip corrupted fragments (vs failing).
        pub skip_corrupted: bool,
        /// Minimum quality threshold before failing (0.0 = accept any).
        pub min_quality_threshold: f32,
        /// Whether to replicate essential data across fragments.
        pub essential_redundancy: bool,
        /// Maximum retry count for corrupted fragments.
        pub max_retries: u32,
    }

    impl Default for FaultToleranceConfig {
        fn default() -> Self {
            Self {
                validate_checksums: true,
                skip_corrupted: true,
                min_quality_threshold: 0.5,
                essential_redundancy: true,
                max_retries: 3,
            }
        }
    }

    /// Fault-tolerant holographic decoder.
    ///
    /// Handles:
    /// - Missing fragments (reconstructs with available data)
    /// - Corrupted fragments (skips and adjusts quality)
    /// - Checksum validation
    /// - Automatic quality adjustment
    pub struct FaultTolerantDecoder {
        ctx: GpuHoloContext,
        config: FaultToleranceConfig,
        header: HoloTensorHeader,
        accumulator: AccumulatorState,
        /// Validation status for each fragment.
        fragment_status: Vec<ValidationResult>,
        /// Number of valid fragments loaded.
        valid_count: u16,
        /// Number of corrupted fragments skipped.
        corrupted_count: u16,
    }

    impl FaultTolerantDecoder {
        /// Creates a new fault-tolerant decoder.
        pub fn new(
            ctx: GpuHoloContext,
            header: HoloTensorHeader,
            config: FaultToleranceConfig,
        ) -> Result<Self, GpuHoloError> {
            let accumulator = ctx.create_accumulator(&header)?;
            let fragment_status = vec![ValidationResult::Missing; header.total_fragments as usize];

            Ok(Self {
                ctx,
                config,
                header,
                accumulator,
                fragment_status,
                valid_count: 0,
                corrupted_count: 0,
            })
        }

        /// Validates a fragment's checksum.
        fn validate_fragment(&self, fragment: &HoloFragment) -> ValidationResult {
            if !self.config.validate_checksums {
                return ValidationResult::Valid;
            }

            // Compute XXH3-64 checksum of fragment data
            let computed = Self::compute_checksum(&fragment.data);
            if computed == fragment.checksum {
                ValidationResult::Valid
            } else {
                ValidationResult::Corrupted
            }
        }

        /// Computes XXH3-64 checksum of data.
        fn compute_checksum(data: &[u8]) -> u64 {
            // Simple FNV-1a hash as fallback (production would use xxhash)
            let mut hash: u64 = 0xcbf29ce484222325;
            for &byte in data {
                hash ^= byte as u64;
                hash = hash.wrapping_mul(0x100000001b3);
            }
            hash
        }

        /// Adds a fragment with validation.
        ///
        /// Returns the validation result and current quality estimate.
        pub fn add_fragment(
            &mut self,
            fragment: &HoloFragment,
        ) -> Result<(ValidationResult, f32), GpuHoloError> {
            let idx = fragment.index as usize;
            if idx >= self.fragment_status.len() {
                return Err(GpuHoloError::InvalidInput {
                    message: format!("Fragment index {} out of range", idx),
                });
            }

            // Validate checksum
            let validation = self.validate_fragment(fragment);

            match validation {
                ValidationResult::Valid => {
                    // Accumulate valid fragment
                    self.ctx.accumulate_fragment(
                        fragment,
                        &mut self.accumulator,
                        self.header.encoding,
                    )?;
                    self.fragment_status[idx] = ValidationResult::Valid;
                    self.valid_count += 1;
                },
                ValidationResult::Corrupted => {
                    self.fragment_status[idx] = ValidationResult::Corrupted;
                    self.corrupted_count += 1;

                    if !self.config.skip_corrupted {
                        return Err(GpuHoloError::FragmentDecode {
                            message: format!("Fragment {} checksum validation failed", idx),
                        });
                    }
                },
                ValidationResult::Missing => {
                    // Should not happen for incoming fragment
                },
            }

            // Calculate quality based on valid fragments only
            let quality = self
                .header
                .quality_curve
                .predict(self.valid_count, self.header.total_fragments);

            Ok((validation, quality))
        }

        /// Returns the number of valid fragments loaded.
        pub fn valid_count(&self) -> u16 {
            self.valid_count
        }

        /// Returns the number of corrupted fragments skipped.
        pub fn corrupted_count(&self) -> u16 {
            self.corrupted_count
        }

        /// Returns the number of missing fragments.
        pub fn missing_count(&self) -> u16 {
            self.fragment_status
                .iter()
                .filter(|&&s| s == ValidationResult::Missing)
                .count() as u16
        }

        /// Returns current quality estimate.
        pub fn quality(&self) -> f32 {
            self.header
                .quality_curve
                .predict(self.valid_count, self.header.total_fragments)
        }

        /// Checks if minimum quality threshold is met.
        pub fn meets_threshold(&self) -> bool {
            self.quality() >= self.config.min_quality_threshold
        }

        /// Returns validation status for each fragment.
        pub fn fragment_status(&self) -> &[ValidationResult] {
            &self.fragment_status
        }

        /// Attempts reconstruction with available valid fragments.
        ///
        /// Returns an error if minimum quality threshold is not met.
        pub fn reconstruct(&self) -> Result<CudaSlice<f32>, GpuHoloError> {
            let quality = self.quality();

            if quality < self.config.min_quality_threshold {
                return Err(GpuHoloError::QualityNotReached {
                    target: self.config.min_quality_threshold,
                    current: quality,
                });
            }

            if self.valid_count < self.header.quality_curve.min_fragments {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: self.header.quality_curve.min_fragments,
                    available: self.valid_count,
                });
            }

            self.ctx
                .finalize_reconstruction(&self.accumulator, self.header.encoding)
        }

        /// Returns statistics about fault tolerance.
        pub fn stats(&self) -> FaultToleranceStats {
            FaultToleranceStats {
                total_fragments: self.header.total_fragments,
                valid_count: self.valid_count,
                corrupted_count: self.corrupted_count,
                missing_count: self.missing_count(),
                current_quality: self.quality(),
                meets_threshold: self.meets_threshold(),
            }
        }
    }

    /// Statistics for fault-tolerant decoding.
    #[derive(Debug, Clone)]
    pub struct FaultToleranceStats {
        /// Total fragments expected.
        pub total_fragments: u16,
        /// Valid fragments loaded.
        pub valid_count: u16,
        /// Corrupted fragments skipped.
        pub corrupted_count: u16,
        /// Missing fragments.
        pub missing_count: u16,
        /// Current quality estimate.
        pub current_quality: f32,
        /// Whether minimum threshold is met.
        pub meets_threshold: bool,
    }

    // ==================== Phase 7: Distributed Loading ====================

    /// Source for fragment data.
    ///
    /// Implement this trait to provide fragments from various sources
    /// (local files, HTTP, S3, peer-to-peer, etc.).
    pub trait FragmentSource: Send + Sync {
        /// Fetches a specific fragment by index.
        fn fetch_fragment(&self, index: u16) -> Result<HoloFragment, GpuHoloError>;

        /// Returns the total number of fragments available.
        fn fragment_count(&self) -> u16;

        /// Returns source priority (lower = prefer this source).
        fn priority(&self) -> u32 {
            100
        }

        /// Returns source name for logging.
        fn name(&self) -> &str;
    }

    /// Configuration for distributed loading.
    #[derive(Debug, Clone)]
    pub struct DistributedLoadConfig {
        /// Maximum concurrent fragment fetches.
        pub max_concurrent: usize,
        /// Timeout per fragment fetch (milliseconds).
        pub fetch_timeout_ms: u64,
        /// Whether to try multiple sources for failed fetches.
        pub failover_enabled: bool,
        /// Target quality to stop loading.
        pub target_quality: f32,
        /// Whether to prefer local sources.
        pub prefer_local: bool,
    }

    impl Default for DistributedLoadConfig {
        fn default() -> Self {
            Self {
                max_concurrent: 4,
                fetch_timeout_ms: 5000,
                failover_enabled: true,
                target_quality: 0.95,
                prefer_local: true,
            }
        }
    }

    /// Distributed fragment loader.
    ///
    /// Fetches fragments from multiple sources in parallel with:
    /// - Priority-based source selection
    /// - Failover on source failure
    /// - Early termination on quality target
    pub struct DistributedLoader {
        sources: Vec<Box<dyn FragmentSource>>,
        config: DistributedLoadConfig,
        header: HoloTensorHeader,
    }

    impl DistributedLoader {
        /// Creates a new distributed loader.
        pub fn new(
            sources: Vec<Box<dyn FragmentSource>>,
            header: HoloTensorHeader,
            config: DistributedLoadConfig,
        ) -> Self {
            Self {
                sources,
                config,
                header,
            }
        }

        /// Adds a fragment source.
        pub fn add_source(&mut self, source: Box<dyn FragmentSource>) {
            self.sources.push(source);
        }

        /// Returns sources sorted by priority.
        fn prioritized_sources(&self) -> Vec<&dyn FragmentSource> {
            let mut sources: Vec<_> = self.sources.iter().map(|s| s.as_ref()).collect();
            sources.sort_by_key(|s| s.priority());
            sources
        }

        /// Fetches a fragment with failover.
        fn fetch_with_failover(&self, index: u16) -> Result<HoloFragment, GpuHoloError> {
            let sources = self.prioritized_sources();

            for source in sources {
                match source.fetch_fragment(index) {
                    Ok(fragment) => return Ok(fragment),
                    Err(_) if self.config.failover_enabled => continue,
                    Err(e) => return Err(e),
                }
            }

            Err(GpuHoloError::FragmentDecode {
                message: format!("Failed to fetch fragment {} from any source", index),
            })
        }

        /// Loads fragments until target quality is reached.
        ///
        /// Returns fragments in priority order (essential first).
        pub fn load_to_quality(&self) -> Result<Vec<HoloFragment>, GpuHoloError> {
            if self.sources.is_empty() {
                return Err(GpuHoloError::InvalidInput {
                    message: "No fragment sources configured".to_string(),
                });
            }

            let mut fragments = Vec::new();
            let mut loaded_indices = std::collections::HashSet::new();

            // Calculate how many fragments needed for target quality
            let needed = self
                .header
                .quality_curve
                .fragments_for_quality(self.config.target_quality, self.header.total_fragments);

            // Load fragments in order (essential first if FLAG_ESSENTIAL_FIRST)
            for index in 0..self.header.total_fragments.min(needed) {
                if loaded_indices.contains(&index) {
                    continue;
                }

                match self.fetch_with_failover(index) {
                    Ok(fragment) => {
                        loaded_indices.insert(index);
                        fragments.push(fragment);
                    },
                    Err(e) if self.config.failover_enabled => {
                        // Skip failed fragment, try to continue
                        continue;
                    },
                    Err(e) => return Err(e),
                }
            }

            if fragments.is_empty() {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: self.header.quality_curve.min_fragments,
                    available: 0,
                });
            }

            Ok(fragments)
        }

        /// Loads all fragments from all sources.
        pub fn load_all(&self) -> Result<Vec<HoloFragment>, GpuHoloError> {
            let mut fragments = Vec::with_capacity(self.header.total_fragments as usize);

            for index in 0..self.header.total_fragments {
                match self.fetch_with_failover(index) {
                    Ok(fragment) => fragments.push(fragment),
                    Err(e) if self.config.failover_enabled => continue,
                    Err(e) => return Err(e),
                }
            }

            Ok(fragments)
        }

        /// Returns statistics about sources.
        pub fn stats(&self) -> DistributedLoadStats {
            DistributedLoadStats {
                source_count: self.sources.len(),
                total_fragments: self.header.total_fragments,
                target_quality: self.config.target_quality,
            }
        }
    }

    /// Statistics for distributed loading.
    #[derive(Debug, Clone)]
    pub struct DistributedLoadStats {
        /// Number of fragment sources.
        pub source_count: usize,
        /// Total fragments in tensor.
        pub total_fragments: u16,
        /// Target quality for loading.
        pub target_quality: f32,
    }

    /// In-memory fragment source for testing.
    pub struct MemoryFragmentSource {
        name: String,
        fragments: Vec<HoloFragment>,
        priority: u32,
    }

    impl MemoryFragmentSource {
        /// Creates a new memory source.
        pub fn new(name: impl Into<String>, fragments: Vec<HoloFragment>, priority: u32) -> Self {
            Self {
                name: name.into(),
                fragments,
                priority,
            }
        }
    }

    impl FragmentSource for MemoryFragmentSource {
        fn fetch_fragment(&self, index: u16) -> Result<HoloFragment, GpuHoloError> {
            self.fragments
                .iter()
                .find(|f| f.index == index)
                .cloned()
                .ok_or_else(|| GpuHoloError::FragmentDecode {
                    message: format!("Fragment {} not found in memory source", index),
                })
        }

        fn fragment_count(&self) -> u16 {
            self.fragments.len() as u16
        }

        fn priority(&self) -> u32 {
            self.priority
        }

        fn name(&self) -> &str {
            &self.name
        }
    }

    // ==================== Phase 7: Adaptive Quality ====================

    /// Quality adjustment policy.
    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
    pub enum QualityPolicy {
        /// Fixed quality target.
        Fixed,
        /// Adjust based on memory pressure.
        MemoryAdaptive,
        /// Adjust based on inference latency.
        LatencyAdaptive,
        /// Best quality within constraints.
        BestEffort,
    }

    /// Configuration for adaptive quality management.
    #[derive(Debug, Clone)]
    pub struct AdaptiveQualityConfig {
        /// Quality adjustment policy.
        pub policy: QualityPolicy,
        /// Minimum acceptable quality.
        pub min_quality: f32,
        /// Maximum acceptable quality (1.0 = full).
        pub max_quality: f32,
        /// Memory threshold to trigger quality reduction (bytes).
        pub memory_threshold: usize,
        /// Latency threshold to trigger quality reduction (ms).
        pub latency_threshold_ms: u64,
        /// How aggressively to adjust quality (0.0-1.0).
        pub adjustment_rate: f32,
    }

    impl Default for AdaptiveQualityConfig {
        fn default() -> Self {
            Self {
                policy: QualityPolicy::BestEffort,
                min_quality: 0.7,
                max_quality: 1.0,
                memory_threshold: 8 * 1024 * 1024 * 1024, // 8 GB
                latency_threshold_ms: 100,
                adjustment_rate: 0.1,
            }
        }
    }

    /// Per-layer quality target.
    #[derive(Debug, Clone)]
    pub struct LayerQualityTarget {
        /// Layer name pattern (supports wildcards).
        pub layer_pattern: String,
        /// Quality target for matching layers.
        pub quality: f32,
        /// Priority (higher = apply first).
        pub priority: u32,
    }

    /// Adaptive quality controller.
    ///
    /// Manages quality targets based on:
    /// - System memory pressure
    /// - Inference latency
    /// - Per-layer importance
    /// - User preferences
    pub struct AdaptiveQualityController {
        config: AdaptiveQualityConfig,
        /// Current target quality.
        current_quality: f32,
        /// Per-layer quality overrides.
        layer_targets: Vec<LayerQualityTarget>,
        /// Memory usage samples.
        memory_samples: Vec<usize>,
        /// Latency samples (ms).
        latency_samples: Vec<u64>,
        /// Sample window size.
        sample_window: usize,
    }

    impl AdaptiveQualityController {
        /// Creates a new adaptive quality controller.
        pub fn new(config: AdaptiveQualityConfig) -> Self {
            let initial_quality = match config.policy {
                QualityPolicy::Fixed => config.max_quality,
                QualityPolicy::MemoryAdaptive => config.max_quality,
                QualityPolicy::LatencyAdaptive => config.max_quality,
                QualityPolicy::BestEffort => config.max_quality,
            };

            Self {
                config,
                current_quality: initial_quality,
                layer_targets: Vec::new(),
                memory_samples: Vec::new(),
                latency_samples: Vec::new(),
                sample_window: 10,
            }
        }

        /// Adds a per-layer quality target.
        pub fn add_layer_target(&mut self, target: LayerQualityTarget) {
            self.layer_targets.push(target);
            self.layer_targets
                .sort_by_key(|t| std::cmp::Reverse(t.priority));
        }

        /// Gets quality target for a specific layer.
        pub fn quality_for_layer(&self, layer_name: &str) -> f32 {
            // Check per-layer overrides
            for target in &self.layer_targets {
                if Self::pattern_matches(&target.layer_pattern, layer_name) {
                    return target
                        .quality
                        .clamp(self.config.min_quality, self.config.max_quality);
                }
            }

            // Fall back to current global quality
            self.current_quality
        }

        /// Checks if a pattern matches a layer name.
        fn pattern_matches(pattern: &str, name: &str) -> bool {
            if pattern == "*" {
                return true;
            }
            if pattern.ends_with('*') {
                let prefix = &pattern[..pattern.len() - 1];
                return name.starts_with(prefix);
            }
            if pattern.starts_with('*') {
                let suffix = &pattern[1..];
                return name.ends_with(suffix);
            }
            pattern == name
        }

        /// Records a memory usage sample.
        pub fn record_memory_usage(&mut self, bytes: usize) {
            self.memory_samples.push(bytes);
            if self.memory_samples.len() > self.sample_window {
                self.memory_samples.remove(0);
            }
            self.update_quality();
        }

        /// Records an inference latency sample.
        pub fn record_latency(&mut self, ms: u64) {
            self.latency_samples.push(ms);
            if self.latency_samples.len() > self.sample_window {
                self.latency_samples.remove(0);
            }
            self.update_quality();
        }

        /// Updates quality target based on samples.
        fn update_quality(&mut self) {
            match self.config.policy {
                QualityPolicy::Fixed => {
                    // No adjustment
                },
                QualityPolicy::MemoryAdaptive => {
                    self.adjust_for_memory();
                },
                QualityPolicy::LatencyAdaptive => {
                    self.adjust_for_latency();
                },
                QualityPolicy::BestEffort => {
                    self.adjust_for_memory();
                    self.adjust_for_latency();
                },
            }
        }

        /// Adjusts quality based on memory pressure.
        fn adjust_for_memory(&mut self) {
            if self.memory_samples.is_empty() {
                return;
            }

            let avg_memory: usize =
                self.memory_samples.iter().sum::<usize>() / self.memory_samples.len();

            if avg_memory > self.config.memory_threshold {
                // Reduce quality
                let reduction = self.config.adjustment_rate;
                self.current_quality = (self.current_quality - reduction)
                    .clamp(self.config.min_quality, self.config.max_quality);
            } else if avg_memory < self.config.memory_threshold / 2 {
                // Increase quality if there's headroom
                let increase = self.config.adjustment_rate / 2.0;
                self.current_quality = (self.current_quality + increase)
                    .clamp(self.config.min_quality, self.config.max_quality);
            }
        }

        /// Adjusts quality based on latency.
        fn adjust_for_latency(&mut self) {
            if self.latency_samples.is_empty() {
                return;
            }

            let avg_latency: u64 =
                self.latency_samples.iter().sum::<u64>() / self.latency_samples.len() as u64;

            if avg_latency > self.config.latency_threshold_ms {
                // Reduce quality for faster inference
                let reduction = self.config.adjustment_rate;
                self.current_quality = (self.current_quality - reduction)
                    .clamp(self.config.min_quality, self.config.max_quality);
            } else if avg_latency < self.config.latency_threshold_ms / 2 {
                // Increase quality if there's headroom
                let increase = self.config.adjustment_rate / 2.0;
                self.current_quality = (self.current_quality + increase)
                    .clamp(self.config.min_quality, self.config.max_quality);
            }
        }

        /// Returns current quality target.
        pub fn current_quality(&self) -> f32 {
            self.current_quality
        }

        /// Sets quality target directly.
        pub fn set_quality(&mut self, quality: f32) {
            self.current_quality = quality.clamp(self.config.min_quality, self.config.max_quality);
        }

        /// Returns statistics about quality management.
        pub fn stats(&self) -> AdaptiveQualityStats {
            let avg_memory = if self.memory_samples.is_empty() {
                0
            } else {
                self.memory_samples.iter().sum::<usize>() / self.memory_samples.len()
            };

            let avg_latency = if self.latency_samples.is_empty() {
                0
            } else {
                self.latency_samples.iter().sum::<u64>() / self.latency_samples.len() as u64
            };

            AdaptiveQualityStats {
                policy: self.config.policy,
                current_quality: self.current_quality,
                min_quality: self.config.min_quality,
                max_quality: self.config.max_quality,
                avg_memory_usage: avg_memory,
                avg_latency_ms: avg_latency,
                layer_target_count: self.layer_targets.len(),
            }
        }
    }

    /// Statistics for adaptive quality.
    #[derive(Debug, Clone)]
    pub struct AdaptiveQualityStats {
        /// Current quality policy.
        pub policy: QualityPolicy,
        /// Current quality target.
        pub current_quality: f32,
        /// Minimum allowed quality.
        pub min_quality: f32,
        /// Maximum allowed quality.
        pub max_quality: f32,
        /// Average memory usage (bytes).
        pub avg_memory_usage: usize,
        /// Average latency (ms).
        pub avg_latency_ms: u64,
        /// Number of layer-specific targets.
        pub layer_target_count: usize,
    }

    // ==================== Hot Reload Support ====================

    /// Hot-reload controller for progressive quality improvement.
    ///
    /// Allows loading initial low-quality weights for fast startup,
    /// then progressively improving quality in the background.
    pub struct HotReloadController {
        ctx: GpuHoloContext,
        header: HoloTensorHeader,
        accumulator: AccumulatorState,
        /// Current quality level.
        current_quality: f32,
        /// Fragments loaded so far.
        fragments_loaded: u16,
        /// Whether initial reconstruction is done.
        initial_ready: bool,
        /// Minimum quality for initial readiness.
        initial_threshold: f32,
    }

    impl HotReloadController {
        /// Creates a new hot-reload controller.
        pub fn new(
            ctx: GpuHoloContext,
            header: HoloTensorHeader,
            initial_threshold: f32,
        ) -> Result<Self, GpuHoloError> {
            let accumulator = ctx.create_accumulator(&header)?;

            Ok(Self {
                ctx,
                header,
                accumulator,
                current_quality: 0.0,
                fragments_loaded: 0,
                initial_ready: false,
                initial_threshold,
            })
        }

        /// Adds a fragment and returns whether quality improved.
        ///
        /// Returns (new_quality, quality_improved).
        pub fn add_fragment(
            &mut self,
            fragment: &HoloFragment,
        ) -> Result<(f32, bool), GpuHoloError> {
            self.ctx
                .accumulate_fragment(fragment, &mut self.accumulator, self.header.encoding)?;
            self.fragments_loaded += 1;

            let new_quality = self
                .header
                .quality_curve
                .predict(self.fragments_loaded, self.header.total_fragments);
            let improved = new_quality > self.current_quality;
            self.current_quality = new_quality;

            if !self.initial_ready && new_quality >= self.initial_threshold {
                self.initial_ready = true;
            }

            Ok((new_quality, improved))
        }

        /// Checks if initial quality threshold is met.
        pub fn is_ready(&self) -> bool {
            self.initial_ready
        }

        /// Returns current quality.
        pub fn quality(&self) -> f32 {
            self.current_quality
        }

        /// Returns number of fragments loaded.
        pub fn fragments_loaded(&self) -> u16 {
            self.fragments_loaded
        }

        /// Reconstructs current state.
        ///
        /// Can be called multiple times as quality improves.
        pub fn reconstruct(&self) -> Result<CudaSlice<f32>, GpuHoloError> {
            if self.fragments_loaded < self.header.quality_curve.min_fragments {
                return Err(GpuHoloError::InsufficientFragments {
                    min_required: self.header.quality_curve.min_fragments,
                    available: self.fragments_loaded,
                });
            }

            self.ctx
                .finalize_reconstruction(&self.accumulator, self.header.encoding)
        }

        /// Returns statistics.
        pub fn stats(&self) -> HotReloadStats {
            HotReloadStats {
                fragments_loaded: self.fragments_loaded,
                total_fragments: self.header.total_fragments,
                current_quality: self.current_quality,
                initial_threshold: self.initial_threshold,
                is_ready: self.initial_ready,
            }
        }
    }

    /// Statistics for hot reload.
    #[derive(Debug, Clone)]
    pub struct HotReloadStats {
        /// Fragments loaded so far.
        pub fragments_loaded: u16,
        /// Total fragments available.
        pub total_fragments: u16,
        /// Current quality level.
        pub current_quality: f32,
        /// Threshold for initial readiness.
        pub initial_threshold: f32,
        /// Whether initial threshold is met.
        pub is_ready: bool,
    }
}

/// Stub for non-CUDA builds.
#[cfg(not(feature = "cuda"))]
pub mod cuda {
    /// Placeholder error for non-CUDA builds.
    #[derive(Debug, thiserror::Error)]
    pub enum GpuHoloError {
        /// CUDA support is not enabled.
        #[error("CUDA support not enabled")]
        CudaNotEnabled,
    }

    /// Placeholder context for non-CUDA builds.
    pub struct GpuHoloContext;

    impl GpuHoloContext {
        /// Creates a new context (always fails without CUDA).
        pub fn new(_device_id: usize) -> Result<Self, GpuHoloError> {
            Err(GpuHoloError::CudaNotEnabled)
        }
    }

    /// Placeholder streaming context for non-CUDA builds.
    pub struct StreamingHoloContext;

    impl StreamingHoloContext {
        /// Creates a new streaming context (always fails without CUDA).
        pub fn new(_device_id: usize, _pipeline_depth: usize) -> Result<Self, GpuHoloError> {
            Err(GpuHoloError::CudaNotEnabled)
        }
    }

    /// Placeholder stream pool for non-CUDA builds.
    pub struct HoloStreamPool;

    /// Placeholder statistics for non-CUDA builds.
    #[derive(Debug, Clone)]
    pub struct StreamingHoloStats {
        /// Pipeline depth used
        pub pipeline_depth: usize,
        /// Number of CUDA streams
        pub num_streams: usize,
    }
}

// ==================== Phase 4: gpu_holo Tests ====================
// Trust boundary §6 (HoloTensor Reconstruction) from GPU-CODEC-PIPELINE-TDD.md.

#[cfg(test)]
mod tests {
    use super::cuda::*;
    #[cfg(feature = "cuda")]
    use cudarc::driver::LaunchAsync;

    /// Stub test: non-CUDA build returns CudaNotEnabled.
    #[test]
    #[cfg(not(feature = "cuda"))]
    fn test_holo_context_stub_not_enabled() {
        match GpuHoloContext::new(0) {
            Err(GpuHoloError::CudaNotEnabled) => {},
            other => panic!("Expected CudaNotEnabled, got {:?}", other.err()),
        }
    }

    /// Stub test: StreamingHoloContext non-CUDA build returns CudaNotEnabled.
    #[test]
    #[cfg(not(feature = "cuda"))]
    fn test_streaming_holo_context_stub_not_enabled() {
        match StreamingHoloContext::new(0, 2) {
            Err(GpuHoloError::CudaNotEnabled) => {},
            other => panic!("Expected CudaNotEnabled, got {:?}", other.err()),
        }
    }

    /// CUDA: context creation succeeds with valid device.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_context_creation() {
        match GpuHoloContext::new(0) {
            Ok(ctx) => {
                assert_eq!(ctx.device_id(), 0);
            },
            Err(GpuHoloError::DeviceInit { .. }) => {
                eprintln!("Skipping: no CUDA device available");
            },
            Err(e) => panic!("Unexpected error: {:?}", e),
        }
    }

    /// CUDA: spectral kernel loading.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_spectral_kernel_load() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };

        ctx.load_spectral_kernel()
            .expect("spectral kernel should load");
        // Second load should be idempotent
        ctx.load_spectral_kernel()
            .expect("second load should succeed");
    }

    /// CUDA: RPH kernel loading.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_rph_kernel_load() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };

        ctx.load_rph_kernel().expect("RPH kernel should load");
    }

    /// CUDA: LRDF kernel loading.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_lrdf_kernel_load() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };

        ctx.load_lrdf_kernel().expect("LRDF kernel should load");
    }

    /// CUDA: invalid device ID returns error.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_invalid_device() {
        let result = GpuHoloContext::new(999);
        assert!(result.is_err(), "Device 999 should not exist");
    }

    /// CUDA: KernelConfig defaults are sane.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_kernel_config_defaults() {
        let ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };

        let config = ctx.kernel_config();
        // Block dimensions should be reasonable (> 0 and <= 1024)
        assert!(
            config.block_size_1d > 0 && config.block_size_1d <= 1024,
            "1D block size should be 1-1024, got {}",
            config.block_size_1d
        );
        assert!(
            config.block_size_2d > 0 && config.block_size_2d <= 32,
            "2D block size should be 1-32, got {}",
            config.block_size_2d
        );
    }

    // ==================== TDD Phase 4: HoloTensor Reconstruction Tests ====================
    // GPU-CODEC-PIPELINE-TDD.md §6.1-6.5

    /// Helper to build an LRDF fragment with one rank-1 SVD component.
    /// Format: [rows: u32][cols: u32][num_components: u32][sigma: f32][u: f32*rows][v: f32*cols]
    #[cfg(feature = "cuda")]
    fn make_lrdf_fragment(
        index: u16,
        rows: usize,
        cols: usize,
        sigma: f32,
        u: &[f32],
        v: &[f32],
    ) -> haagenti::holotensor::HoloFragment {
        assert_eq!(u.len(), rows);
        assert_eq!(v.len(), cols);
        let mut data = Vec::new();
        data.extend_from_slice(&(rows as u32).to_le_bytes());
        data.extend_from_slice(&(cols as u32).to_le_bytes());
        data.extend_from_slice(&1u32.to_le_bytes()); // 1 component
        data.extend_from_slice(&sigma.to_le_bytes());
        for &val in u {
            data.extend_from_slice(&val.to_le_bytes());
        }
        for &val in v {
            data.extend_from_slice(&val.to_le_bytes());
        }
        haagenti::holotensor::HoloFragment::new(index, data)
    }

    /// §6.2: LRDF rank-1 outer product reconstruction.
    ///
    /// A single fragment with one component: sigma * u * v^T.
    /// Verifies GPU reconstruction matches the expected outer product.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_lrdf_rank1_reconstruction() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_lrdf_kernel().expect("LRDF kernel should load");

        let rows = 3;
        let cols = 2;
        let u = vec![1.0f32, 2.0, 3.0];
        let v = vec![4.0f32, 5.0];
        let sigma = 2.0f32;

        let fragment = make_lrdf_fragment(0, rows, cols, sigma, &u, &v);

        // Build header
        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::LowRankDistributed,
            haagenti::DType::F32,
            vec![rows as u64, cols as u64],
            1, // total fragments
        );

        // Reconstruct
        let gpu_result = ctx.reconstruct(&header, &[fragment]).unwrap();
        let host = ctx.copy_to_host(&gpu_result).unwrap();

        // Expected: sigma * u_i * v_j (row-major)
        let expected = vec![
            2.0 * 1.0 * 4.0,
            2.0 * 1.0 * 5.0, // row 0: 8, 10
            2.0 * 2.0 * 4.0,
            2.0 * 2.0 * 5.0, // row 1: 16, 20
            2.0 * 3.0 * 4.0,
            2.0 * 3.0 * 5.0, // row 2: 24, 30
        ];

        assert_eq!(host.len(), expected.len());
        for (i, (e, g)) in expected.iter().zip(host.iter()).enumerate() {
            assert!(
                (e - g).abs() < 1e-4,
                "LRDF rank-1 mismatch at {}: expected={}, got={}",
                i,
                e,
                g
            );
        }
    }

    /// §6.2: LRDF multi-component reconstruction (2 rank-1 terms).
    ///
    /// Two fragments each contributing one rank-1 component.
    /// Result should be sigma1 * u1 * v1^T + sigma2 * u2 * v2^T.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_lrdf_multi_component_reconstruction() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_lrdf_kernel().expect("LRDF kernel should load");

        let rows = 2;
        let cols = 2;

        // Component 1: 1.0 * [1, 0] * [1, 0]^T = [[1, 0], [0, 0]]
        let frag0 = make_lrdf_fragment(0, rows, cols, 1.0, &[1.0, 0.0], &[1.0, 0.0]);
        // Component 2: 1.0 * [0, 1] * [0, 1]^T = [[0, 0], [0, 1]]
        let frag1 = make_lrdf_fragment(1, rows, cols, 1.0, &[0.0, 1.0], &[0.0, 1.0]);

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::LowRankDistributed,
            haagenti::DType::F32,
            vec![rows as u64, cols as u64],
            2,
        );

        let gpu_result = ctx.reconstruct(&header, &[frag0, frag1]).unwrap();
        let host = ctx.copy_to_host(&gpu_result).unwrap();

        // Expected: identity matrix [[1, 0], [0, 1]]
        let expected = vec![1.0, 0.0, 0.0, 1.0];

        assert_eq!(host.len(), expected.len());
        for (i, (e, g)) in expected.iter().zip(host.iter()).enumerate() {
            assert!(
                (e - g).abs() < 1e-4,
                "LRDF multi-comp mismatch at {}: expected={}, got={}",
                i,
                e,
                g
            );
        }
    }

    /// §6.5: F32 → F16 conversion kernel correctness.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_f32_to_f16_conversion() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_fused_kernel().expect("fused kernel should load");

        let input = vec![1.0f32, 0.5, -1.0, 0.0, 100.0, -0.001];
        let d_input = ctx.device().htod_copy(input.clone()).unwrap();

        let d_output = ctx.convert_f32_to_f16(&d_input).unwrap();

        let mut host_f16 = vec![half::f16::ZERO; input.len()];
        ctx.device()
            .dtoh_sync_copy_into(&d_output, &mut host_f16)
            .unwrap();

        for (i, (&f32_val, &f16_val)) in input.iter().zip(host_f16.iter()).enumerate() {
            let expected = half::f16::from_f32(f32_val);
            assert_eq!(
                f16_val.to_bits(),
                expected.to_bits(),
                "F32→F16 wrong at {}: input={}, got=0x{:04X}, expected=0x{:04X}",
                i,
                f32_val,
                f16_val.to_bits(),
                expected.to_bits()
            );
        }
    }

    /// §6.5: scale_values utility kernel.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_holo_scale_values() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_fused_kernel().expect("fused kernel should load");

        let data = vec![1.0f32, 2.0, 3.0, 4.0];
        let mut d_data = ctx.device().htod_copy(data.clone()).unwrap();

        ctx.scale_values(&mut d_data, 0.5).unwrap();

        let mut host = vec![0.0f32; data.len()];
        ctx.device()
            .dtoh_sync_copy_into(&d_data, &mut host)
            .unwrap();

        let expected = vec![0.5, 1.0, 1.5, 2.0];
        for (i, (e, g)) in expected.iter().zip(host.iter()).enumerate() {
            assert!(
                (e - g).abs() < 1e-6,
                "Scale mismatch at {}: expected={}, got={}",
                i,
                e,
                g
            );
        }
    }

    /// §6.2: LRDF convenience reconstruct_lrdf produces correct GpuTensor.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_lrdf_convenience_reconstruct() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_lrdf_kernel().expect("LRDF kernel should load");

        let rows = 4;
        let cols = 4;
        // Identity-like: sigma=1 for each basis vector pair
        let frag = make_lrdf_fragment(
            0,
            rows,
            cols,
            1.0,
            &[1.0, 1.0, 1.0, 1.0],
            &[1.0, 1.0, 1.0, 1.0],
        );

        let tensor = ctx.reconstruct_lrdf(&[frag], rows, cols).unwrap();

        assert_eq!(tensor.rows(), rows);
        assert_eq!(tensor.cols(), cols);
        assert_eq!(tensor.len(), rows * cols);

        let host = tensor.to_host().unwrap();
        // All 1s outer product: every element should be 1.0
        for (i, &val) in host.iter().enumerate() {
            assert!(
                (val - 1.0).abs() < 1e-4,
                "reconstruct_lrdf all-ones mismatch at {}: got={}",
                i,
                val
            );
        }
    }

    // ==================== Phase 5: §L4 LRDF Batched Outer Product ====================

    /// Helper to build a multi-component LRDF fragment.
    /// Format: [rows: u32][cols: u32][num_components: u32]([sigma: f32][u: f32*rows][v: f32*cols])*
    #[cfg(feature = "cuda")]
    fn make_lrdf_multi_fragment(
        index: u16,
        rows: usize,
        cols: usize,
        components: &[(f32, Vec<f32>, Vec<f32>)], // [(sigma, u, v), ...]
    ) -> haagenti::holotensor::HoloFragment {
        let mut data = Vec::new();
        data.extend_from_slice(&(rows as u32).to_le_bytes());
        data.extend_from_slice(&(cols as u32).to_le_bytes());
        data.extend_from_slice(&(components.len() as u32).to_le_bytes());
        for (sigma, u, v) in components {
            assert_eq!(u.len(), rows);
            assert_eq!(v.len(), cols);
            data.extend_from_slice(&sigma.to_le_bytes());
            for &val in u {
                data.extend_from_slice(&val.to_le_bytes());
            }
            for &val in v {
                data.extend_from_slice(&val.to_le_bytes());
            }
        }
        haagenti::holotensor::HoloFragment::new(index, data)
    }

    /// §L4.1: Single-component batched matches unbatched single outer product.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_lrdf_batched_single_matches_unbatched() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_lrdf_kernel().expect("LRDF kernel should load");

        let rows = 4;
        let cols = 3;
        let sigma = 2.5f32;
        let u = vec![1.0, 0.5, 0.25, 0.125];
        let v = vec![0.3, 0.6, 0.9];

        // Unbatched: single outer product
        let frag_single = make_lrdf_fragment(0, rows, cols, sigma, &u, &v);
        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::LowRankDistributed,
            haagenti::DType::F32,
            vec![rows as u64, cols as u64],
            1,
        );
        let mut acc_single = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_lrdf(&frag_single, &mut acc_single).unwrap();
        let result_single = ctx.finalize_lrdf(&acc_single).unwrap();
        let host_single = ctx.copy_to_host(&result_single).unwrap();

        // Batched: same data through batched path
        let frag_batched = make_lrdf_multi_fragment(0, rows, cols, &[(sigma, u, v)]);
        let mut acc_batched = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_lrdf_batched(&frag_batched, &mut acc_batched)
            .unwrap();
        let result_batched = ctx.finalize_lrdf(&acc_batched).unwrap();
        let host_batched = ctx.copy_to_host(&result_batched).unwrap();

        assert_eq!(host_single.len(), host_batched.len());
        for (i, (&s, &b)) in host_single.iter().zip(host_batched.iter()).enumerate() {
            assert!(
                (s - b).abs() < 1e-5,
                "Batched vs unbatched mismatch at {}: single={}, batched={}",
                i,
                s,
                b
            );
        }
    }

    /// §L4.2: Multi-component batched matches sequential single calls.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_lrdf_batched_multi_matches_sequential() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_lrdf_kernel().expect("LRDF kernel should load");

        let rows = 4;
        let cols = 3;
        let components = vec![
            (2.0f32, vec![1.0, 0.0, 0.0, 0.0], vec![1.0, 0.0, 0.0]),
            (1.5f32, vec![0.0, 1.0, 0.0, 0.0], vec![0.0, 1.0, 0.0]),
            (0.8f32, vec![0.0, 0.0, 1.0, 0.0], vec![0.0, 0.0, 1.0]),
        ];

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::LowRankDistributed,
            haagenti::DType::F32,
            vec![rows as u64, cols as u64],
            1,
        );

        // Sequential: 3 separate single-component fragments
        let mut acc_seq = ctx.create_accumulator(&header).unwrap();
        for (i, (sigma, ref u, ref v)) in components.iter().enumerate() {
            let frag = make_lrdf_fragment(i as u16, rows, cols, *sigma, u, v);
            ctx.accumulate_lrdf(&frag, &mut acc_seq).unwrap();
        }
        let result_seq = ctx.finalize_lrdf(&acc_seq).unwrap();
        let host_seq = ctx.copy_to_host(&result_seq).unwrap();

        // Batched: all 3 components in single fragment
        let frag_batched = make_lrdf_multi_fragment(0, rows, cols, &components);
        let mut acc_batched = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_lrdf_batched(&frag_batched, &mut acc_batched)
            .unwrap();
        let result_batched = ctx.finalize_lrdf(&acc_batched).unwrap();
        let host_batched = ctx.copy_to_host(&result_batched).unwrap();

        assert_eq!(host_seq.len(), host_batched.len());
        for (i, (&s, &b)) in host_seq.iter().zip(host_batched.iter()).enumerate() {
            assert!(
                (s - b).abs() < 1e-5,
                "Batched vs sequential mismatch at {}: seq={}, batched={}",
                i,
                s,
                b
            );
        }

        // Verify expected values: diagonal-like pattern
        // output[0][0] = 2.0, output[1][1] = 1.5, output[2][2] = 0.8, rest ~0
        assert!(
            (host_batched[0] - 2.0).abs() < 1e-5,
            "Expected 2.0 at [0,0]"
        );
        assert!(
            (host_batched[4] - 1.5).abs() < 1e-5,
            "Expected 1.5 at [1,1]"
        );
        assert!(
            (host_batched[8] - 0.8).abs() < 1e-5,
            "Expected 0.8 at [2,2]"
        );
    }

    // ==================== Phase 4: §6.1 Spectral IDCT Reconstruction ====================
    // GPU-CODEC-PIPELINE-TDD.md §6.1: Spectral accumulation and IDCT reconstruction.

    /// Helper to build a legacy spectral fragment.
    /// Format: [num_coeffs: u32][indices: u32...][values: f32...]
    #[cfg(feature = "cuda")]
    fn make_spectral_fragment(
        frag_index: u16,
        coeffs: &[(u32, f32)],
    ) -> haagenti::holotensor::HoloFragment {
        let num_coeffs = coeffs.len() as u32;
        let mut data = Vec::new();
        data.extend_from_slice(&num_coeffs.to_le_bytes());
        for &(idx, _) in coeffs {
            data.extend_from_slice(&idx.to_le_bytes());
        }
        for &(_, val) in coeffs {
            data.extend_from_slice(&val.to_le_bytes());
        }
        haagenti::holotensor::HoloFragment::new(frag_index, data)
    }

    /// §6.1 + §S1.1: DC-only spectral coefficient accumulates correctly and
    /// reconstructs to a constant output.
    ///
    /// A single DC coefficient (index 0) should land at position 0 in the
    /// coefficient buffer with all other positions zero. Full reconstruction
    /// (accumulate → IDCT rows → IDCT cols) should produce a uniform constant
    /// equal to DC / sqrt(width) (for height=1, col IDCT is identity).
    #[test]
    #[cfg(feature = "cuda")]
    fn test_spectral_dc_only_produces_constant_output() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_spectral_kernel()
            .expect("spectral kernel should load");

        let width = 8;
        let height = 1;
        let dc_value = 4.0f32;

        let fragment = make_spectral_fragment(0, &[(0, dc_value)]);

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::Spectral,
            haagenti::DType::F32,
            vec![height as u64, width as u64],
            1,
        );

        // Accumulate and verify coefficients directly
        let mut accumulator = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_spectral(&fragment, &mut accumulator)
            .unwrap();

        if let AccumulatorState::Spectral {
            ref coefficients, ..
        } = accumulator
        {
            let host_coeffs = ctx.copy_to_host(coefficients).unwrap();
            assert_eq!(host_coeffs.len(), width * height);

            // DC coefficient at index 0 should equal the input value
            assert!(
                (host_coeffs[0] - dc_value).abs() < 1e-6,
                "DC coefficient: expected {}, got {}",
                dc_value,
                host_coeffs[0]
            );

            // All other positions should be zero
            for (i, &val) in host_coeffs[1..].iter().enumerate() {
                assert!(
                    val.abs() < 1e-6,
                    "Non-DC position {} should be 0, got {}",
                    i + 1,
                    val
                );
            }
        } else {
            panic!("Expected Spectral accumulator variant");
        }

        // Full reconstruction: DC-only IDCT should produce constant output.
        // For 1D (height=1): x[n] = DC * sqrt(2/N) * 1/sqrt(2) = DC / sqrt(N)
        let result = ctx.finalize_spectral(&accumulator).unwrap();
        let output = ctx.copy_to_host(&result).unwrap();
        let expected = dc_value / (width as f32).sqrt();
        for (i, &val) in output.iter().enumerate() {
            assert!(
                (val - expected).abs() < 1e-3,
                "DC reconstruction at {}: expected {:.6}, got {:.6}",
                i,
                expected,
                val
            );
        }
    }

    /// §6.1: Sparse spectral coefficients accumulate at correct indices.
    ///
    /// Verifies the accumulation kernel places coefficient values at the
    /// correct positions in the frequency buffer, with zeros elsewhere.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_spectral_sparse_accumulation_correct_indices() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_spectral_kernel()
            .expect("spectral kernel should load");

        let width = 4;
        let height = 4;
        let total_size = width * height;

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::Spectral,
            haagenti::DType::F32,
            vec![height as u64, width as u64],
            1,
        );

        // Create accumulator and accumulate a sparse fragment
        let mut accumulator = ctx.create_accumulator(&header).unwrap();
        let coeffs = [(0u32, 1.0f32), (3, 2.0), (7, 3.0), (15, 4.0)];
        let fragment = make_spectral_fragment(0, &coeffs);
        ctx.accumulate_spectral(&fragment, &mut accumulator)
            .unwrap();

        // Read back coefficient buffer from accumulator
        if let AccumulatorState::Spectral {
            ref coefficients, ..
        } = accumulator
        {
            let host_coeffs = ctx.copy_to_host(coefficients).unwrap();
            assert_eq!(host_coeffs.len(), total_size);

            // Indexed positions should have the accumulated values
            assert!(
                (host_coeffs[0] - 1.0).abs() < 1e-6,
                "Index 0: expected 1.0, got {}",
                host_coeffs[0]
            );
            assert!(
                (host_coeffs[3] - 2.0).abs() < 1e-6,
                "Index 3: expected 2.0, got {}",
                host_coeffs[3]
            );
            assert!(
                (host_coeffs[7] - 3.0).abs() < 1e-6,
                "Index 7: expected 3.0, got {}",
                host_coeffs[7]
            );
            assert!(
                (host_coeffs[15] - 4.0).abs() < 1e-6,
                "Index 15: expected 4.0, got {}",
                host_coeffs[15]
            );

            // All non-indexed positions should be zero
            for (i, &val) in host_coeffs.iter().enumerate() {
                if ![0, 3, 7, 15].contains(&i) {
                    assert!(
                        val.abs() < 1e-6,
                        "Non-indexed position {} should be 0, got {}",
                        i,
                        val
                    );
                }
            }
        } else {
            panic!("Expected Spectral accumulator variant");
        }
    }

    // ==================== DD-8 §S1.2: Single AC Coefficient Row IDCT ====================

    /// §S1.2: A single AC coefficient at index k produces a cosine wave.
    /// x[n] = amplitude * sqrt(2/N) * cos(π(2n+1)k / 2N)
    #[test]
    #[cfg(feature = "cuda")]
    fn test_idct_1d_rows_single_ac() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_spectral_kernel()
            .expect("spectral kernel should load");

        let width = 8;
        let height = 1;
        let k = 1usize; // First AC frequency
        let amplitude = 3.0f32;

        // Create coefficient buffer with single AC at index k
        let mut coeffs = vec![0.0f32; width * height];
        coeffs[k] = amplitude;

        // Upload and run row IDCT directly
        let d_input = ctx.device().htod_copy(coeffs).unwrap();
        let d_output: cudarc::driver::CudaSlice<f32> =
            ctx.device().alloc_zeros(width * height).unwrap();

        let func = ctx
            .device()
            .get_func("holo_spectral", "holo_spectral_idct_1d_rows")
            .unwrap();
        let cfg = cudarc::driver::LaunchConfig {
            grid_dim: (1, 1, 1),
            block_dim: (256, 1, 1),
            shared_mem_bytes: (width * 4) as u32,
        };
        unsafe { func.launch(cfg, (&d_input, &d_output, width as u32, height as u32)) }.unwrap();
        ctx.device().synchronize().unwrap();

        let host = ctx.copy_to_host(&d_output).unwrap();
        let scale = (2.0 / width as f32).sqrt();
        let pi_2n = std::f32::consts::PI / (2.0 * width as f32);

        for n in 0..width {
            let expected = amplitude * scale * ((2 * n + 1) as f32 * k as f32 * pi_2n).cos();
            assert!(
                (host[n] - expected).abs() < 1e-3,
                "AC[{}] row IDCT at {}: expected {:.6}, got {:.6}",
                k,
                n,
                expected,
                host[n]
            );
        }
    }

    // ==================== DD-8 §S2.3: End-to-End 2D IDCT ====================

    /// §S2.3: Full 2D IDCT pipeline — accumulate DCT coefficients, finalize,
    /// compare with CPU reference (`haagenti_core::dct::idct_1d_direct`).
    ///
    /// Uses a known 4×4 spatial signal → DCT-II → GPU IDCT → compare.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_spectral_2d_idct_end_to_end() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_spectral_kernel()
            .expect("spectral kernel should load");

        let width = 4;
        let height = 4;

        // Known spatial signal
        let spatial: Vec<f32> = (0..width * height)
            .map(|i| (i as f32 * 0.3).sin())
            .collect();

        // Compute DCT-II coefficients using CPU reference
        let mut dct_coeffs = vec![0.0f32; width * height];
        haagenti_core::dct::dct_2d(&spatial, &mut dct_coeffs, width, height);

        // Build spectral fragments from all non-zero coefficients
        let mut coeff_pairs: Vec<(u32, f32)> = Vec::new();
        for (i, &c) in dct_coeffs.iter().enumerate() {
            if c.abs() > 1e-10 {
                coeff_pairs.push((i as u32, c));
            }
        }
        let fragment = make_spectral_fragment(0, &coeff_pairs);

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::Spectral,
            haagenti::DType::F32,
            vec![height as u64, width as u64],
            1,
        );

        // GPU reconstruction
        let mut acc = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_spectral(&fragment, &mut acc).unwrap();
        let result = ctx.finalize_spectral(&acc).unwrap();
        let gpu_output = ctx.copy_to_host(&result).unwrap();

        // CPU reference IDCT for comparison
        let mut cpu_output = vec![0.0f32; width * height];
        haagenti_core::dct::idct_2d(&dct_coeffs, &mut cpu_output, width, height);

        for (i, (&orig, &recon)) in cpu_output.iter().zip(gpu_output.iter()).enumerate() {
            assert!(
                (orig - recon).abs() < 1e-2,
                "2D IDCT at {}: cpu_reference={:.6}, gpu={:.6}, diff={:.6}",
                i,
                orig,
                recon,
                (orig - recon).abs()
            );
        }
    }

    // ==================== DD-8 §S1.3: GPU vs CPU Proptest ====================

    #[cfg(feature = "cuda")]
    mod idct_proptest {
        use super::*;
        use proptest::prelude::*;

        proptest! {
            #![proptest_config(ProptestConfig::with_cases(20))]
            #[test]
            fn gpu_idct_2d_matches_cpu_reference(
                dim in prop::sample::select(vec![4usize, 8]),
                coeffs_raw in proptest::collection::vec(-5.0f32..5.0f32, 64),
            ) {
                let mut ctx = match GpuHoloContext::new(0) {
                    Ok(ctx) => ctx,
                    Err(_) => { return Ok(()); }
                };
                ctx.load_spectral_kernel().unwrap();

                let width = dim;
                let height = dim;
                let total = width * height;
                let coeffs: Vec<f32> = coeffs_raw[..total].to_vec();

                // CPU reference
                let mut cpu_output = vec![0.0f32; total];
                haagenti_core::dct::idct_2d(&coeffs, &mut cpu_output, width, height);

                // GPU: upload coefficients directly, run finalize_spectral_direct
                let d_coeffs = ctx.device().htod_copy(coeffs).unwrap();

                // Create spectral accumulator and inject coefficients directly
                let acc = AccumulatorState::Spectral {
                    coefficients: d_coeffs,
                    present_mask: ctx.device().alloc_zeros(total).unwrap(),
                    width,
                    height,
                };
                let result = ctx.finalize_spectral(&acc).unwrap();
                let gpu_output = ctx.copy_to_host(&result).unwrap();

                for (i, (&cpu, &gpu)) in cpu_output.iter().zip(gpu_output.iter()).enumerate() {
                    prop_assert!(
                        (cpu - gpu).abs() < 0.05,
                        "IDCT 2D mismatch at {}: cpu={:.6}, gpu={:.6}, diff={:.6}",
                        i, cpu, gpu, (cpu - gpu).abs()
                    );
                }
            }
        }
    }

    // ==================== Phase 4: §6.3 RPH Determinism ====================
    // GPU-CODEC-PIPELINE-TDD.md §6.3: Same seed produces same output.

    /// Helper to build an RPH fragment.
    /// Format: [proj_dim: u32][seed_offset: u64][projection: f32...]
    #[cfg(feature = "cuda")]
    fn make_rph_fragment(
        frag_index: u16,
        proj_dim: usize,
        seed_offset: u64,
        projection: &[f32],
    ) -> haagenti::holotensor::HoloFragment {
        assert_eq!(projection.len(), proj_dim);
        let mut data = Vec::new();
        data.extend_from_slice(&(proj_dim as u32).to_le_bytes());
        data.extend_from_slice(&seed_offset.to_le_bytes());
        for &val in projection {
            data.extend_from_slice(&val.to_le_bytes());
        }
        haagenti::holotensor::HoloFragment::new(frag_index, data)
    }

    /// §R2.2 + §R2.3: RPH reconstruction produces non-zero output and is deterministic.
    ///
    /// Two runs of accumulate+finalize with identical inputs (same header, same seed,
    /// same fragments) must produce bit-identical, non-zero outputs.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_rph_deterministic_same_seed() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_rph_kernel().expect("RPH kernel should load");

        let proj_dim = 8;
        let projection: Vec<f32> = (0..proj_dim).map(|i| (i as f32 + 1.0) * 0.5).collect();

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::RandomProjection,
            haagenti::DType::F32,
            vec![4, 4], // 4x4 = 16 output elements
            1,
        );

        // First run: accumulate + finalize
        let frag1 = make_rph_fragment(0, proj_dim, 42, &projection);
        let mut acc1 = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_rph(&frag1, &mut acc1).unwrap();
        let result1 = ctx.finalize_rph(&acc1).unwrap();
        let host1 = ctx.copy_to_host(&result1).unwrap();

        // Second run with identical inputs
        let frag2 = make_rph_fragment(0, proj_dim, 42, &projection);
        let mut acc2 = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_rph(&frag2, &mut acc2).unwrap();
        let result2 = ctx.finalize_rph(&acc2).unwrap();
        let host2 = ctx.copy_to_host(&result2).unwrap();

        assert_eq!(host1.len(), host2.len(), "Output lengths must match");
        assert_eq!(host1.len(), 16, "Output should be 4x4 = 16 elements");

        // §R2.2: Output must be non-zero (XORShift PRNG fixed)
        assert!(
            host1.iter().any(|v| *v != 0.0),
            "RPH output should be non-zero after XORShift PRNG fix"
        );

        // §R2.3: Determinism — bit-identical across runs
        for (i, (v1, v2)) in host1.iter().zip(host2.iter()).enumerate() {
            assert_eq!(
                v1.to_bits(),
                v2.to_bits(),
                "RPH determinism violation at {}: run1={}, run2={}",
                i,
                v1,
                v2
            );
        }
    }

    /// §R2.3: Different seeds produce different RPH output.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_rph_different_seeds_diverge() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_rph_kernel().expect("RPH kernel should load");

        let proj_dim = 8;
        let projection: Vec<f32> = (0..proj_dim).map(|i| (i as f32 + 1.0) * 0.5).collect();

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::RandomProjection,
            haagenti::DType::F32,
            vec![4, 4],
            1,
        );

        // Run with seed_offset=42
        let frag_a = make_rph_fragment(0, proj_dim, 42, &projection);
        let mut acc_a = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_rph(&frag_a, &mut acc_a).unwrap();
        let result_a = ctx.finalize_rph(&acc_a).unwrap();
        let host_a = ctx.copy_to_host(&result_a).unwrap();

        // Run with seed_offset=999
        let frag_b = make_rph_fragment(0, proj_dim, 999, &projection);
        let mut acc_b = ctx.create_accumulator(&header).unwrap();
        ctx.accumulate_rph(&frag_b, &mut acc_b).unwrap();
        let result_b = ctx.finalize_rph(&acc_b).unwrap();
        let host_b = ctx.copy_to_host(&result_b).unwrap();

        // Different seeds must produce different output
        let differs = host_a
            .iter()
            .zip(host_b.iter())
            .any(|(a, b)| a.to_bits() != b.to_bits());
        assert!(
            differs,
            "Different seed offsets should produce different RPH output"
        );
    }

    /// §6.3: RPH accumulate pipeline exercises fragment parsing and state tracking.
    ///
    /// Verifies the full RPH pipeline: fragment construction, kernel load,
    /// accumulation, and finalization all complete without error.
    /// The accumulator state (num_projections) is correctly incremented.
    #[test]
    #[cfg(feature = "cuda")]
    fn test_rph_accumulate_pipeline_state_tracking() {
        let mut ctx = match GpuHoloContext::new(0) {
            Ok(ctx) => ctx,
            Err(_) => {
                eprintln!("Skipping: no CUDA device available");
                return;
            },
        };
        ctx.load_rph_kernel().expect("RPH kernel should load");

        let frag_proj_dim = 8;
        let projection: Vec<f32> = (0..frag_proj_dim).map(|i| (i as f32 + 1.0) * 0.5).collect();

        let header = haagenti::holotensor::HoloTensorHeader::new(
            haagenti::holotensor::HolographicEncoding::RandomProjection,
            haagenti::DType::F32,
            vec![4, 4],
            1,
        );

        let mut acc = ctx.create_accumulator(&header).unwrap();

        // Verify initial state.
        // Accumulator proj_dim is compute_projection_dim(16) = max(sqrt(16), 16) = 16,
        // which differs from the fragment's proj_dim (8).
        if let AccumulatorState::RandomProjection {
            num_projections,
            proj_dim: pd,
            output_dim: od,
            ..
        } = &acc
        {
            assert_eq!(*num_projections, 0, "Initial num_projections should be 0");
            assert_eq!(*pd, 16, "proj_dim = max(sqrt(output_dim), 16) = 16");
            assert_eq!(*od, 16, "output_dim should be 4*4=16");
        } else {
            panic!("Expected RandomProjection accumulator");
        }

        // Accumulate two fragments with different seed_offsets
        let frag_a = make_rph_fragment(0, frag_proj_dim, 42, &projection);
        ctx.accumulate_rph(&frag_a, &mut acc).unwrap();

        if let AccumulatorState::RandomProjection {
            num_projections, ..
        } = &acc
        {
            assert_eq!(
                *num_projections, 1,
                "num_projections should be 1 after first accumulate"
            );
        }

        let frag_b = make_rph_fragment(1, frag_proj_dim, 999, &projection);
        ctx.accumulate_rph(&frag_b, &mut acc).unwrap();

        if let AccumulatorState::RandomProjection {
            num_projections, ..
        } = &acc
        {
            assert_eq!(
                *num_projections, 2,
                "num_projections should be 2 after second accumulate"
            );
        }

        // Finalize and verify output shape
        let result = ctx.finalize_rph(&acc).unwrap();
        let host = ctx.copy_to_host(&result).unwrap();
        assert_eq!(host.len(), 16, "Finalized output should be 4*4=16 elements");

        // After XORShift PRNG fix, output should be non-zero
        assert!(
            host.iter().any(|v| *v != 0.0),
            "RPH pipeline output should be non-zero after PRNG fix"
        );
    }
}

// Re-exports
pub use cuda::GpuHoloContext;
#[cfg(feature = "cuda")]
pub use cuda::GpuHoloError;
#[cfg(feature = "cuda")]
pub use cuda::HoloStreamPool;
pub use cuda::StreamingHoloContext;
#[cfg(feature = "cuda")]
pub use cuda::StreamingHoloStats;