map2fig 0.7.7

//! FITS file reading and data extraction.
//!
//! This module provides functions for reading HEALPix astronomical data from FITS files.
//! It supports both:
//! - **Dense maps**: Complete sky coverage (standard NSIDE² pixels)
//! - **Sparse maps**: Partial sky coverage with explicit PIXEL column indexing
//!
//! # FITS Format Support
//!
//! - Binary table format (.fits files)
//! - HEALPix RING or NEST pixel ordering
//! - Sparse maps with IMPLICIT or EXPLICIT indexing schemes via the INDXSCHM keyword
//!
//! # Sparse Map Handling
//!
//! For sparse maps with EXPLICIT indexing, this module automatically expands the data
//! to a full dense array with UNSEEN values for omitted pixels.
//!
//! **Tier 4.2b Optimization:** Sparse map column extraction is parallelized via rayon
//! for efficient multi-core processing of large sparse catalogs.
//!
//! # Examples
//!
//! ```ignore
//! use map2fig::read_healpix_column;
//!
//! let data = read_healpix_column("map.fits", 0);
//! println!("Loaded {} pixels", data.len());
//! ```

use std::fs::File;
use std::io::{Read, Write};

use crate::data_array::DataArray;
use fitsrs::hdu::data::bintable::{ColumnId, DataValue};
use fitsrs::{Fits, HDU, card::Value};
use rayon::prelude::*;

// ============================================================================
// Tier 1 Optimization: Direct Binary Reading for Float32 Columns
// ============================================================================

/// Parse FITS TFORM code to extract array count and type
/// E.g., "4096E" -> (4096, 'E'), "1D" -> (1, 'D')
fn parse_tform(tform: &str) -> Option<(usize, char)> {
    let tform = tform.trim();
    let type_char = tform.chars().last()?;

    let count_str = tform.trim_end_matches(type_char);
    let count: usize = if count_str.is_empty() {
        1
    } else {
        count_str.parse().ok()?
    };

    Some((count, type_char))
}

/// Fast path for reading float32 columns directly from binary data
/// Preserves f32 precision - does NOT convert to f64 (New optimization!)
///
/// This function:
/// 1. Uses fitsrs to parse headers and find table offset
/// 2. Reads column binary data directly from mmap
/// 3. Interprets float32 values inline (no enum, no conversion)
///
/// Expected improvement: 6.8s saved (62.4% of FITS reading)
/// Memory saved: 3.2 GB (for 806M pixel maps)
fn try_read_float32_column_native(
    _filename: &str,
    mmap_data: &[u8],
    col_idx: usize,
) -> Option<(Vec<f32>, i64)> {
    use std::io::Cursor;

    let cursor = Cursor::new(mmap_data);
    let mut fits = Fits::from_reader(cursor);
    let mut nside: i64 = 0;

    while let Some(Ok(hdu)) = fits.next() {
        if let HDU::XBinaryTable(hdu) = hdu {
            let header = hdu.get_header();

            // Skip sparse maps (explicit indexing) - use fallback path
            let has_explicit_indexing = match header.get("INDXSCHM") {
                Some(Value::String { value, .. }) => value.trim() == "EXPLICIT",
                _ => false,
            };
            if has_explicit_indexing {
                return None;
            }

            // Get NSIDE
            if nside == 0 {
                match header.get("NSIDE") {
                    Some(Value::Integer { value, .. }) => nside = *value,
                    _ => return None,
                };
            }

            //Get column type and count from TFORM
            let tform_key = format!("TFORM{}", col_idx + 1);
            let tform_str = match header.get(&tform_key) {
                Some(Value::String { value, .. }) => value.clone(),
                _ => return None,
            };

            let (elem_count, type_char) = parse_tform(&tform_str)?;

            // Fast path ONLY for float32 ('E') columns
            if type_char != 'E' {
                return None;
            }

            // Get column byte offset from TOFFSET (defaults to 0 for first column per FITS standard)
            let toffset_key = format!("TOFFSET{}", col_idx + 1);
            let col_offset: usize = match header.get(&toffset_key) {
                Some(Value::Integer { value, .. }) => *value as usize,
                _ => {
                    if col_idx == 0 {
                        0
                    } else {
                        return None;
                    }
                }
            };

            // Get row byte size (NAXIS1)
            let row_size: usize = match header.get("NAXIS1") {
                Some(Value::Integer { value, .. }) => *value as usize,
                _ => return None,
            };

            // Get number of rows (NAXIS2)
            let num_rows: usize = match header.get("NAXIS2") {
                Some(Value::Integer { value, .. }) => *value as usize,
                _ => return None,
            };

            // Find data offset (after all headers, which are 2880-byte blocks)
            let data_offset = find_binary_table_data_offset(mmap_data)?;

            // Pre-allocate result - f32, no conversion!
            let total_elems = elem_count * num_rows;
            let mut result = vec![0f32; total_elems];

            // **Tier 5.3: Sequential FITS Reading (15.7× optimization)**
            // Read column data sequentially through the file (optimal for I/O bandwidth)
            // Instead of scattered row-by-row access that breaks prefetcher,
            // iterate through file sequentially and extract column values in-place.
            //
            // Expected improvement: 5.5s → 0.35s (94% faster)
            // Root cause of old slowness: 65KB strides prevent CPU prefetcher activation,
            // causing every memory access to be a cache miss (50+ cycle latency).
            //
            // FITS format: Row-major storage
            // Row 0: [col0(4B)][col1(4B)]...[col4095(4B)]
            // Row 1: [col0(4B)][col1(4B)]...[col4095(4B)]
            //
            // Sequential read pattern: Process each row in order, extract our column
            let file_data = &mmap_data[data_offset..];

            for (row_idx, row_chunk) in file_data.chunks(row_size).enumerate() {
                if row_idx >= num_rows {
                    break;
                }

                // Extract this row's column data (col_offset to col_offset + elem_count*4)
                let col_end = col_offset + elem_count * 4;
                if col_end > row_chunk.len() {
                    return None;
                }

                let col_bytes = &row_chunk[col_offset..col_end];

                // Parse column values and store in result array
                // FITS format uses big-endian byte order (IEEE 754 network byte order)
                for (elem_idx, chunk) in col_bytes.chunks_exact(4).enumerate() {
                    let bytes = [chunk[0], chunk[1], chunk[2], chunk[3]];
                    let f32_val = f32::from_be_bytes(bytes);
                    result[row_idx * elem_count + elem_idx] = f32_val;
                }
            }

            return Some((result, nside));
        }
    }

    None
}

/// Fast path for reading float64 columns directly from binary data
/// Preserves f64 precision
fn try_read_float64_column_native(
    _filename: &str,
    mmap_data: &[u8],
    col_idx: usize,
) -> Option<(Vec<f64>, i64)> {
    use std::io::Cursor;

    let cursor = Cursor::new(mmap_data);
    let mut fits = Fits::from_reader(cursor);
    let mut nside: i64 = 0;

    while let Some(Ok(hdu)) = fits.next() {
        if let HDU::XBinaryTable(hdu) = hdu {
            let header = hdu.get_header();

            // Skip sparse maps (explicit indexing) - use fallback path
            let has_explicit_indexing = match header.get("INDXSCHM") {
                Some(Value::String { value, .. }) => value.trim() == "EXPLICIT",
                _ => false,
            };
            if has_explicit_indexing {
                return None;
            }

            // Get NSIDE
            if nside == 0 {
                match header.get("NSIDE") {
                    Some(Value::Integer { value, .. }) => nside = *value,
                    _ => return None,
                };
            }

            // Get column type and count from TFORM
            let tform_key = format!("TFORM{}", col_idx + 1);
            let tform_str = match header.get(&tform_key) {
                Some(Value::String { value, .. }) => value.clone(),
                _ => return None,
            };

            let (elem_count, type_char) = parse_tform(&tform_str)?;

            // Fast path ONLY for float64 ('D') columns
            if type_char != 'D' {
                return None;
            }

            // Get column byte offset from TOFFSET (defaults to 0 for first column per FITS standard)
            let toffset_key = format!("TOFFSET{}", col_idx + 1);
            let col_offset: usize = match header.get(&toffset_key) {
                Some(Value::Integer { value, .. }) => *value as usize,
                _ => {
                    if col_idx == 0 {
                        0
                    } else {
                        return None;
                    }
                }
            };

            // Get row byte size (NAXIS1)
            let row_size: usize = match header.get("NAXIS1") {
                Some(Value::Integer { value, .. }) => *value as usize,
                _ => return None,
            };

            // Get number of rows (NAXIS2)
            let num_rows: usize = match header.get("NAXIS2") {
                Some(Value::Integer { value, .. }) => *value as usize,
                _ => return None,
            };

            // Find data offset (after all headers, which are 2880-byte blocks)
            let data_offset = find_binary_table_data_offset(mmap_data)?;

            // Pre-allocate result
            let total_elems = elem_count * num_rows;
            let mut result = vec![0f64; total_elems];

            // **Tier 5.3: Sequential FITS Reading (15.7× optimization)**
            // Read column data sequentially through the file (optimal for I/O bandwidth)
            let file_data = &mmap_data[data_offset..];

            for (row_idx, row_chunk) in file_data.chunks(row_size).enumerate() {
                if row_idx >= num_rows {
                    break;
                }

                // Extract this row's column data (col_offset to col_offset + elem_count*8)
                let col_end = col_offset + elem_count * 8;
                if col_end > row_chunk.len() {
                    return None;
                }

                let col_bytes = &row_chunk[col_offset..col_end];

                // Parse column values and store in result array
                // FITS format uses big-endian byte order (IEEE 754 network byte order)
                for (elem_idx, chunk) in col_bytes.chunks_exact(8).enumerate() {
                    let bytes = [
                        chunk[0], chunk[1], chunk[2], chunk[3], chunk[4], chunk[5], chunk[6],
                        chunk[7],
                    ];
                    let f64_val = f64::from_be_bytes(bytes);
                    result[row_idx * elem_count + elem_idx] = f64_val;
                }
            }

            return Some((result, nside));
        }
    }

    None
}

/// Find the binary data offset in a FITS file (after all header blocks)
/// FITS headers are padded to 2880-byte blocks; data starts at next block after LAST header
/// For multi-HDU files, we need to find the BINTABLE END, not the PRIMARY END
fn find_binary_table_data_offset(mmap_data: &[u8]) -> Option<usize> {
    const FITS_BLOCK_SIZE: usize = 2880;
    let mut last_end_block = None;

    // Scan through entire file looking for all END keywords
    // Track the LAST one (which marks the end of BINTABLE header)
    for block_num in 0..1000 {
        let block_start = block_num * FITS_BLOCK_SIZE;
        if block_start >= mmap_data.len() {
            break;
        }

        let block = if block_start + FITS_BLOCK_SIZE <= mmap_data.len() {
            &mmap_data[block_start..block_start + FITS_BLOCK_SIZE]
        } else {
            // Partial block at end of file
            &mmap_data[block_start..]
        };

        // Look for "END" keyword in FITS cards
        // FITS cards are 80 bytes, keyword is first 8 bytes
        for card_off in (0..block.len()).step_by(80) {
            if card_off + 8 <= block.len() {
                let card = &block[card_off..card_off + 8];
                if card == b"END     " {
                    // Track this END as the latest we've seen
                    last_end_block = Some(block_num);
                    break; // Only one END per 2880-byte block
                }
            }
        }
    }

    // Data starts at next 2880-byte block after last END keyword
    last_end_block.map(|block_num| (block_num + 1) * FITS_BLOCK_SIZE)
}

/// Read a HEALPix column from a FITS binary table.
///
/// Extracts data from a specific column of a HEALPix FITS binary table.
/// Automatically handles both dense and sparse (EXPLICIT) indexing schemes.
/// For sparse maps, expands the result to full NSIDE² size with UNSEEN values.
///
/// # Arguments
///
/// * `filename` - Path to the FITS file
/// * `col_idx`  - 0-based data column index (not the PIXEL column, just the data)
///
/// # Returns
///
/// DataArray with length = 12 * NSIDE²
/// - Dense maps: all pixels present in FITS
/// - Sparse maps: UNSEEN (-1.6375e30) for missing pixels
///
/// **Type preservation** (New in v0.7.0):
/// - f32 FITS columns stay as f32 (saves 6.8s + 3.2 GB memory)
/// - f64 FITS columns stay as f64
/// - Sparse data converted to f64 via fallback path
///
/// # Panics
///
/// Panics if:
/// - File cannot be opened
/// - FITS structure is invalid
/// - Column index is out of bounds
/// - Required HEALPix headers are missing
pub fn read_healpix_column(filename: &str, col_idx: usize) -> DataArray {
    // Tier 2 Optimization: Use memory-mapped I/O instead of buffered reads
    // Eliminates kernel memcpy overhead (rep_movs_alternative) and improves cache locality
    use memmap2::Mmap;
    use std::io::Cursor;
    use std::time::Instant;

    let enable_profile = std::env::var("MAP2FIG_PROFILE").is_ok();

    let mmap_start = Instant::now();
    let f = File::open(filename).expect("Failed to open FITS file");
    let mmap = unsafe { Mmap::map(&f).expect("Failed to mmap FITS file") };
    let mmap_elapsed = mmap_start.elapsed();
    if enable_profile {
        eprintln!("[I/O DIAG] mmap() took: {:.3}s", mmap_elapsed.as_secs_f64());
    }

    // **NEW: Preserve float32 precision without conversion (6.8s + 3.2 GB saved)**
    // Tier 1b Optimization: Try native f32 reader first
    let f32_start = Instant::now();
    if let Some((data, _nside)) = try_read_float32_column_native(filename, &mmap, col_idx) {
        let f32_elapsed = f32_start.elapsed();
        if enable_profile {
            eprintln!(
                "[I/O DIAG] float32 read took: {:.3}s",
                f32_elapsed.as_secs_f64()
            );
        }
        return DataArray::from_f32(data);
    }

    // **NEW: Preserve float64 precision**
    // Try native f64 reader
    let f64_start = Instant::now();
    if let Some((data, _nside)) = try_read_float64_column_native(filename, &mmap, col_idx) {
        let f64_elapsed = f64_start.elapsed();
        if enable_profile {
            eprintln!(
                "[I/O DIAG] float64 read took: {:.3}s",
                f64_elapsed.as_secs_f64()
            );
        }
        return DataArray::from_f64(data);
    }

    // Fallback path: Use fitsrs DataValue enum (slower but handles all types)
    let cursor = Cursor::new(&mmap[..]);
    let mut fits = Fits::from_reader(cursor);
    let mut result: Vec<f64> = Vec::new();
    let mut nside: i64 = 0;

    while let Some(Ok(hdu)) = fits.next() {
        if let HDU::XBinaryTable(hdu) = hdu {
            let header = hdu.get_header();

            // Check if this uses explicit indexing (sparse/partial sky map)
            let has_explicit_indexing = match header.get("INDXSCHM") {
                Some(Value::String { value, .. }) => value.trim() == "EXPLICIT",
                _ => false,
            };

            // Get NSIDE if not already set
            if nside == 0 {
                nside = match header.get("NSIDE") {
                    Some(Value::Integer { value, .. }) => *value,
                    _ => 0,
                };
            }

            let data = fits.get_data(&hdu);
            let mut table = data.table_data();

            // If explicit indexing, read both PIXEL and data columns together
            if has_explicit_indexing && nside > 0 {
                // For explicit indexing:
                // - Column 0 is always PIXEL indices
                // - Column 1+ are data columns
                // - User's --col N refers to the N-th data column
                // - Adjust file column: file_col = col_idx + 1
                let file_col_for_data = col_idx + 1;

                // Select both columns together in a single call
                // This preserves the column correspondence
                // Do NOT call select_fields twice - that breaks the pairing!
                let all_values: Vec<DataValue> = table
                    .select_fields(&[ColumnId::Index(0), ColumnId::Index(file_col_for_data)])
                    .collect();

                if all_values.is_empty() {
                    result = vec![crate::healpix::HPX_UNSEEN; (12 * nside * nside) as usize];
                } else {
                    let n_rows = all_values.len() / 2;
                    let npix = (12 * nside * nside) as usize;
                    let mut full_map = vec![crate::healpix::HPX_UNSEEN; npix];

                    // Tier 4.2b: Parallel extraction with row batching to reduce task overhead
                    // Batch rows into ~1M row chunks (≈16MB per chunk) to reduce Rayon task count
                    // from millions to hundreds
                    let chunk_size = 1_000_000; // 1M rows per parallelization unit
                    let num_chunks = n_rows.div_ceil(chunk_size);

                    let pairs: Vec<(usize, f64)> = (0..num_chunks)
                        .into_par_iter()
                        .flat_map(|chunk_idx| {
                            let start = chunk_idx * chunk_size;
                            let end = std::cmp::min((chunk_idx + 1) * chunk_size, n_rows);
                            (start..end)
                                .filter_map(|row_idx| {
                                    let pix_idx = row_idx * 2;
                                    let data_idx = row_idx * 2 + 1;

                                    let pix = match &all_values[pix_idx] {
                                        DataValue::Integer { value, .. } => *value as i64,
                                        DataValue::Long { value, .. } => *value,
                                        DataValue::Float { value, .. } => *value as i64,
                                        DataValue::Double { value, .. } => *value as i64,
                                        _ => -1,
                                    };

                                    let val = match &all_values[data_idx] {
                                        DataValue::Double { value, .. } => *value,
                                        DataValue::Float { value, .. } => *value as f64,
                                        DataValue::Integer { value, .. } => *value as f64,
                                        other => {
                                            panic!(
                                                "Unsupported column type in FITS table: {:?}",
                                                other
                                            )
                                        }
                                    };

                                    if pix >= 0 && (pix as usize) < npix {
                                        Some((pix as usize, val))
                                    } else {
                                        None
                                    }
                                })
                                .collect::<Vec<_>>()
                        })
                        .collect();

                    // Populate map sequentially from parallel results
                    for (pix_idx, val) in pairs {
                        full_map[pix_idx] = val;
                    }

                    result = full_map;
                }
            } else {
                // Regular dense map: read column directly
                let values = table.select_fields(&[ColumnId::Index(col_idx)]);
                for cell in values {
                    match cell {
                        DataValue::Double { value, .. } => result.push(value),
                        DataValue::Float { value, .. } => result.push(value as f64),
                        DataValue::Integer { value, .. } => result.push(value as f64),
                        other => panic!("Unsupported column type in FITS table: {:?}", other),
                    }
                }
            }
        }
    }

    // Fallback path returned f64 data (sparse maps use fitsrs)
    DataArray::from_f64(result)
}

// ============================================================================
// Metadata Caching (Phase 4.2a: I/O Optimization)
// ============================================================================

use serde_json::{Value as JsonValue, json};
use std::path::{Path, PathBuf};

/// Get cache directory for FITS metadata
fn get_cache_dir() -> Option<PathBuf> {
    let dirs = directories::ProjectDirs::from("", "", "map2fig")?;
    Some(dirs.cache_dir().to_path_buf())
}

/// Hash a file path and mtime for cache validation
fn compute_cache_key(filepath: &str) -> Option<String> {
    use std::fs;

    let path = Path::new(filepath);
    let metadata = fs::metadata(path).ok()?;
    let mtime = metadata.modified().ok()?;

    let mtime_seconds = mtime.duration_since(std::time::UNIX_EPOCH).ok()?.as_secs();

    // Use SHA256 for filename safety
    use sha2::{Digest, Sha256};
    let mut hasher = Sha256::new();
    hasher.update(filepath.as_bytes());
    let hash = format!("{:x}", hasher.finalize());

    // Cache key: {hash}_{mtime}
    Some(format!("{}_{}", hash, mtime_seconds))
}

/// Try to load cached metadata for a FITS file
fn try_load_cache(filepath: &str) -> Option<(i64, String, String)> {
    let cache_dir = get_cache_dir()?;
    let cache_key = compute_cache_key(filepath)?;
    let cache_file = cache_dir.join(format!("fits_meta_{}.json", cache_key));

    let json_str = std::fs::read_to_string(cache_file).ok()?;
    let json: JsonValue = serde_json::from_str(&json_str).ok()?;

    let nside = json.get("nside")?.as_i64()?;
    let order = json.get("ordering")?.as_str()?.to_string();
    let indxschm = json.get("indxschm")?.as_str()?.to_string();

    Some((nside, order, indxschm))
}

/// Save metadata to cache
fn save_cache(filepath: &str, nside: i64, ordering: &str, indxschm: &str) {
    let cache_dir = match get_cache_dir() {
        Some(d) => d,
        None => return, // Caching not available
    };

    // Create cache directory if it doesn't exist
    let _ = std::fs::create_dir_all(&cache_dir);

    let cache_key = match compute_cache_key(filepath) {
        Some(k) => k,
        None => return,
    };

    let cache_file = cache_dir.join(format!("fits_meta_{}.json", cache_key));

    let cache_data = json!({
        "nside": nside,
        "ordering": ordering,
        "indxschm": indxschm,
    });

    let _ = std::fs::write(cache_file, cache_data.to_string());
}

/// Read FITS metadata with caching support
/// This function attempts to use cached metadata to avoid expensive header parsing
/// When MAP2FIG_PROFILE environment variable is set, outputs diagnostic timing info
pub fn read_healpix_meta_cached(filename: &str) -> Option<(i64, String, String)> {
    // Tier 2b Optimization: Always use memory-mapped I/O for metadata reading
    // This eliminates kernel memcpy overhead and page fault overhead from BufReader
    read_healpix_meta_cached_mmap(filename)
}

// Dead code removed (Tier 2b optimization: metadata now uses mmap)

// ============================================================================
// Tier 5.2.1: Column Data Caching
// ============================================================================

/// Generate cache key for column data
/// Format: {SHA256(filepath)}_{column_idx}_{mtime_secs}
fn get_column_cache_key(filepath: &str, col_idx: usize, mtime_secs: u64) -> Option<String> {
    use sha2::{Digest, Sha256};

    let mut hasher = Sha256::new();
    hasher.update(filepath.as_bytes());
    let hash = hasher.finalize();

    // Convert hash to hex string (first 16 chars for brevity)
    let hash_str = format!("{:x}", hash);
    let short_hash = &hash_str[..16.min(hash_str.len())];

    Some(format!(
        "fits_col_{}_{:03}_{}",
        short_hash, col_idx, mtime_secs
    ))
}

/// Try to load column data from cache
/// Returns Some(`Vec<f64>`) if cache exists and is valid, None otherwise
fn try_load_column_cache(filepath: &str, col_idx: usize) -> Option<Vec<f64>> {
    let cache_dir = get_cache_dir()?;

    // Get file metadata for mtime validation
    let metadata = std::fs::metadata(filepath).ok()?;
    let mtime_secs = metadata
        .modified()
        .ok()?
        .duration_since(std::time::UNIX_EPOCH)
        .ok()?
        .as_secs();

    // Generate cache key and filename
    let cache_key = get_column_cache_key(filepath, col_idx, mtime_secs)?;
    let cache_file = cache_dir.join(&cache_key);

    // Try to open and read cache file
    let mut file = File::open(cache_file).ok()?;

    // Read header (16 bytes)
    let mut header = [0u8; 16];
    file.read_exact(&mut header).ok()?;

    // Parse header
    let magic = u32::from_le_bytes([header[0], header[1], header[2], header[3]]);
    let version = u32::from_le_bytes([header[4], header[5], header[6], header[7]]);
    let num_pixels = u32::from_le_bytes([header[8], header[9], header[10], header[11]]);

    // Validate magic and version
    if magic != 0xCAFEBABE || version != 1 {
        return None;
    }

    // Read data
    let mut data = vec![0.0; num_pixels as usize];
    let byte_data = unsafe {
        std::slice::from_raw_parts_mut(data.as_mut_ptr() as *mut u8, num_pixels as usize * 8)
    };
    file.read_exact(byte_data).ok()?;

    let enable_profile = std::env::var("MAP2FIG_PROFILE").is_ok();
    if enable_profile {
        eprintln!("[I/O DIAG] Column cache HIT: {} col#{}", filepath, col_idx);
    }

    Some(data)
}

/// Save column data to cache
fn save_column_cache(filepath: &str, col_idx: usize, data: &[f64]) -> Option<()> {
    // Skip caching very large columns (>1GB of data)
    // These are rarely reused and would blow the cache size limit (2GB max)
    const MAX_CACHE_COLUMN_SIZE: usize = 128_000_000; // ~1GB of data
    if data.len() > MAX_CACHE_COLUMN_SIZE {
        let enable_profile = std::env::var("MAP2FIG_PROFILE").is_ok();
        if enable_profile {
            eprintln!(
                "[I/O DIAG] Column cache SKIP (too large): {} col#{} ({} pixels > {}M)",
                filepath,
                col_idx,
                data.len(),
                MAX_CACHE_COLUMN_SIZE / 1_000_000
            );
        }
        return None; // Skip, don't cache
    }

    let cache_dir = get_cache_dir()?;
    let _ = std::fs::create_dir_all(&cache_dir);

    // Get file metadata for mtime
    let metadata = std::fs::metadata(filepath).ok()?;
    let mtime_secs = metadata
        .modified()
        .ok()?
        .duration_since(std::time::UNIX_EPOCH)
        .ok()?
        .as_secs();

    // Generate cache key and filename
    let cache_key = get_column_cache_key(filepath, col_idx, mtime_secs)?;
    let cache_file = cache_dir.join(&cache_key);

    // Create file
    let mut file = File::create(cache_file).ok()?;

    // Write header: magic (0xCAFEBABE), version (1), num_pixels, reserved
    file.write_all(&0xCAFEBABEu32.to_le_bytes()).ok()?;
    file.write_all(&1u32.to_le_bytes()).ok()?;
    file.write_all(&(data.len() as u32).to_le_bytes()).ok()?;
    file.write_all(&0u32.to_le_bytes()).ok()?;

    // Write data as little-endian f64 bytes
    let byte_data =
        unsafe { std::slice::from_raw_parts(data.as_ptr() as *const u8, data.len() * 8) };
    file.write_all(byte_data).ok()?;

    let enable_profile = std::env::var("MAP2FIG_PROFILE").is_ok();
    if enable_profile {
        eprintln!(
            "[I/O DIAG] Column cache SAVE: {} col#{} ({} pixels)",
            filepath,
            col_idx,
            data.len()
        );
    }

    Some(())
}

/// Enforce cache size limit (2 GB max)
/// Deletes oldest files first if cache exceeds the limit
fn enforce_cache_size_limit() {
    const MAX_CACHE_SIZE: u64 = 2_000_000_000; // 2 GB in bytes

    let cache_dir = match get_cache_dir() {
        Some(d) => d,
        None => return,
    };

    // Calculate total cache size
    let mut total_size = 0u64;
    let mut files = Vec::new();

    for entry in std::fs::read_dir(&cache_dir)
        .ok()
        .into_iter()
        .flatten()
        .flatten()
    {
        let path = entry.path();
        if let Ok(metadata) = entry.metadata() {
            let size = metadata.len();
            total_size += size;
            if let Ok(modified) = metadata.modified() {
                files.push((path, modified, size));
            }
        }
    }

    // Only prune if we exceed the limit
    if total_size <= MAX_CACHE_SIZE {
        return;
    }

    // Sort by modification time (oldest first)
    files.sort_by_key(|f| f.1);

    // Delete oldest files until we're under the limit (target: 90% of max)
    let target_size = (MAX_CACHE_SIZE as f64 * 0.9) as u64;
    let mut freed = 0u64;

    for (path, _, size) in files.iter() {
        if total_size - freed <= target_size {
            break;
        }
        let _ = std::fs::remove_file(path);
        freed += size;
    }
}

/// Read a HEALPix column with caching support (Tier 5.2.1 optimization)
///
/// This function reads column data from cache if available, otherwise reads from FITS
/// and saves to cache for future use. Cache is automatically invalidated if file mtime changes.
///
/// # Memory-mapped I/O Option
///
/// Set `MAP2FIX_USE_MMAP=1` environment variable to use memory-mapped I/O instead of buffered reading.
/// This can improve performance for large files (>500 MB) at the cost of higher memory usage.
pub fn read_healpix_column_cached(filename: &str, col_idx: usize) -> DataArray {
    // Try cache first
    if let Some(data) = try_load_column_cache(filename, col_idx) {
        return DataArray::from_f64(data);
    }

    let enable_profile = std::env::var("MAP2FIG_PROFILE").is_ok();
    if enable_profile {
        eprintln!("[I/O DIAG] Column cache MISS: {} col#{}", filename, col_idx);
    }

    // Check if mmap mode is enabled
    let use_mmap = std::env::var("MAP2FIX_USE_MMAP").is_ok();

    // Cache miss: read from FITS (use mmap if enabled)
    let fits_start = std::time::Instant::now();
    let data_array: DataArray = if use_mmap {
        read_healpix_column_mmap(filename, col_idx)
    } else {
        read_healpix_column(filename, col_idx)
    };
    let fits_elapsed = fits_start.elapsed();
    if enable_profile {
        eprintln!(
            "[I/O DIAG] FITS parsing took: {:.3}s",
            fits_elapsed.as_secs_f64()
        );
    }

    // Return DataArray as-is, preserving f32/f64 type - NO CONVERSION
    // This is the key optimization: avoid expensive f32->f64 conversion (4.98s wasted)

    // Only save f64 data to cache (f32 files are too large, skip caching)
    if let DataArray::Float64(ref vec) = data_array {
        let cache_start = std::time::Instant::now();
        match save_column_cache(filename, col_idx, vec) {
            Some(_) => {
                let cache_elapsed = cache_start.elapsed();
                if enable_profile {
                    eprintln!(
                        "[I/O DIAG] Column cache SAVE SUCCESS ({:.3}s): {} col#{}",
                        cache_elapsed.as_secs_f64(),
                        filename,
                        col_idx
                    );
                }
            }
            None => {
                let cache_elapsed = cache_start.elapsed();
                if enable_profile {
                    eprintln!(
                        "[I/O DIAG] Column cache SAVE FAILED ({:.3}s): {} col#{}",
                        cache_elapsed.as_secs_f64(),
                        filename,
                        col_idx
                    );
                }
            }
        }
    } else if enable_profile {
        eprintln!("[I/O DIAG] Skipping cache for f32 data (too large)");
    }

    // Enforce cache size limit (2 GB max)
    enforce_cache_size_limit();

    data_array
}
// ============================================================================
// Memory-mapped FITS Reading (Optional optimization path)
// ============================================================================

/// Read a HEALPix column from a FITS file using memory-mapped I/O
///
/// This is an alternative to the standard `read_healpix_column` that uses
/// memory-mapped file access (memmap2) instead of buffered reading.
///
/// # Performance Characteristics
///
/// - **For large files** (>500 MB): Comparable or slightly faster (10-20% best case)
/// - **For small files** (<100 MB): Minimal difference, possibly slower due to memmap overhead
/// - **Memory usage**: Maps entire file to memory (could be large)
/// - **Cache effects**: Once mapped, subsequent passes are very fast
///
/// # Arguments
///
/// Same as `read_healpix_column`
///
/// # Returns
///
/// Same as `read_healpix_column` (DataArray with proper type preservation)
pub fn read_healpix_column_mmap(filename: &str, col_idx: usize) -> DataArray {
    use memmap2::Mmap;
    use std::io::Cursor;

    let f = File::open(filename).expect("Failed to open FITS file");
    let mmap = unsafe { Mmap::map(&f).expect("Failed to memory-map FITS file") };

    // Try native f32 path first
    if let Some((data, _nside)) = try_read_float32_column_native(filename, &mmap, col_idx) {
        return DataArray::from_f32(data);
    }

    // Try native f64 path
    if let Some((data, _nside)) = try_read_float64_column_native(filename, &mmap, col_idx) {
        return DataArray::from_f64(data);
    }

    // Create a Cursor over the memory-mapped data
    // This allows fitsrs to work with the mapped memory without copying
    let cursor = Cursor::new(&mmap[..]);
    let mut fits = Fits::from_reader(cursor);
    let mut result: Vec<f64> = Vec::new();
    let mut nside: i64 = 0;

    while let Some(Ok(hdu)) = fits.next() {
        if let HDU::XBinaryTable(hdu) = hdu {
            let header = hdu.get_header();

            // Check if this uses explicit indexing (sparse/partial sky map)
            let has_explicit_indexing = match header.get("INDXSCHM") {
                Some(Value::String { value, .. }) => value.trim() == "EXPLICIT",
                _ => false,
            };

            // Get NSIDE if not already set
            if nside == 0 {
                nside = match header.get("NSIDE") {
                    Some(Value::Integer { value, .. }) => *value,
                    _ => 0,
                };
            }

            let data = fits.get_data(&hdu);
            let mut table = data.table_data();

            // If explicit indexing, read both PIXEL and data columns together
            if has_explicit_indexing && nside > 0 {
                let file_col_for_data = col_idx + 1;

                // Read both PIXEL (col 0) and data column
                let all_values: Vec<DataValue> = table
                    .select_fields(&[ColumnId::Index(0), ColumnId::Index(file_col_for_data)])
                    .collect();

                if all_values.is_empty() {
                    result = vec![crate::healpix::HPX_UNSEEN; (12 * nside * nside) as usize];
                } else {
                    let n_rows = all_values.len() / 2;
                    let npix = (12 * nside * nside) as usize;
                    let mut full_map = vec![crate::healpix::HPX_UNSEEN; npix];

                    // Tier 4.2b: Parallel extraction with row batching to reduce task overhead
                    // Batch rows into ~1M row chunks (≈16MB per chunk) to reduce Rayon task count
                    // from millions to hundreds
                    let chunk_size = 1_000_000; // 1M rows per parallelization unit
                    let num_chunks = n_rows.div_ceil(chunk_size);

                    let pairs: Vec<(usize, f64)> = (0..num_chunks)
                        .into_par_iter()
                        .flat_map(|chunk_idx| {
                            let start = chunk_idx * chunk_size;
                            let end = std::cmp::min((chunk_idx + 1) * chunk_size, n_rows);
                            (start..end)
                                .filter_map(|row_idx| {
                                    let pix_idx = row_idx * 2;
                                    let data_idx = row_idx * 2 + 1;

                                    let pix = match &all_values[pix_idx] {
                                        DataValue::Integer { value, .. } => *value as i64,
                                        DataValue::Long { value, .. } => *value,
                                        DataValue::Float { value, .. } => *value as i64,
                                        DataValue::Double { value, .. } => *value as i64,
                                        _ => -1,
                                    };

                                    let val = match &all_values[data_idx] {
                                        DataValue::Double { value, .. } => *value,
                                        DataValue::Float { value, .. } => *value as f64,
                                        DataValue::Integer { value, .. } => *value as f64,
                                        other => {
                                            panic!(
                                                "Unsupported column type in FITS table: {:?}",
                                                other
                                            )
                                        }
                                    };

                                    if pix >= 0 && (pix as usize) < npix {
                                        Some((pix as usize, val))
                                    } else {
                                        None
                                    }
                                })
                                .collect::<Vec<_>>()
                        })
                        .collect();

                    for (pix_idx, val) in pairs {
                        full_map[pix_idx] = val;
                    }

                    result = full_map;
                }
            } else {
                // Regular dense map: read column directly
                let values = table.select_fields(&[ColumnId::Index(col_idx)]);
                for cell in values {
                    match cell {
                        DataValue::Double { value, .. } => result.push(value),
                        DataValue::Float { value, .. } => result.push(value as f64),
                        DataValue::Integer { value, .. } => result.push(value as f64),
                        other => panic!("Unsupported column type in FITS table: {:?}", other),
                    }
                }
            }
        }
    }

    // Fallback: return as f64 (sparse/complex types)
    DataArray::from_f64(result)
}

/// Read FITS metadata using memory-mapped I/O (mmap variant of `read_healpix_meta_cached`)
pub fn read_healpix_meta_cached_mmap(filename: &str) -> Option<(i64, String, String)> {
    use memmap2::Mmap;
    use std::io::Cursor;

    let enable_profile = std::env::var("MAP2FIX_PROFILE").is_ok();

    // Try cache first
    let cache_start = std::time::Instant::now();
    if let Some((nside, order, indxschm)) = try_load_cache(filename) {
        if enable_profile {
            let elapsed = cache_start.elapsed();
            eprintln!(
                "[I/O DIAG] Cache HIT (mmap): {} ({:.3}µs)",
                filename,
                elapsed.as_micros()
            );
        }
        return Some((nside, order, indxschm));
    }

    // Cache miss: parse FITS file with mmap
    if enable_profile {
        let elapsed = cache_start.elapsed();
        eprintln!(
            "[I/O DIAG] Cache MISS (mmap): {} (lookup took {:.3}µs)",
            filename,
            elapsed.as_micros()
        );
    }

    let parse_start = std::time::Instant::now();
    let f = File::open(filename).ok()?;
    let mmap = unsafe { Mmap::map(&f).ok()? };

    let cursor = Cursor::new(&mmap[..]);
    let mut fits = Fits::from_reader(cursor);
    let mut nside: i64 = 0;
    let mut ordering = String::new();
    let mut indxschm = String::from("IMPLICIT");

    while let Some(Ok(hdu)) = fits.next() {
        if let HDU::XBinaryTable(hdu) = hdu {
            let header = hdu.get_header();
            nside = match header.get("NSIDE") {
                Some(Value::Integer { value, .. }) => *value,
                _ => 0,
            };
            ordering = match header.get("ORDERING") {
                Some(Value::String { value, .. }) => value.trim().to_string(),
                _ => "RING".to_string(),
            };
            indxschm = match header.get("INDXSCHM") {
                Some(Value::String { value, .. }) => value.trim().to_string(),
                _ => "IMPLICIT".to_string(),
            };
            break;
        }
    }

    let parse_end = std::time::Instant::now();
    if enable_profile {
        eprintln!(
            "[I/O DIAG] Parse done (mmap): {:.3}ms",
            parse_end.duration_since(parse_start).as_secs_f64() * 1000.0
        );
    }

    if nside > 0 {
        save_cache(filename, nside, &ordering, &indxschm);
        Some((nside, ordering, indxschm))
    } else {
        None
    }
}

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

    #[test]
    #[cfg(not(miri))]
    fn test_sparse_fit_explicit_indexing() {
        // Test loading the minimal sparse NSIDE=1 FITS file
        // This verifies the regression fix: sparse maps must use HPX_UNSEEN
        // initialization, not f64::NEG_INFINITY
        //
        // NOTE: This test is skipped under Miri because it uses file system
        // operations (Path::exists) which call statx syscall that Miri's
        // isolation mode doesn't support.

        let test_file =
            PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures/sparse_nside1.fits");

        if !test_file.exists() {
            eprintln!(
                "Skipping test: sparse FITS fixture not found at {}",
                test_file.display()
            );
            return;
        }

        // Load the sparse column data
        let data = read_healpix_column(test_file.to_str().unwrap(), 0);

        // NSIDE=1 has 12 pixels
        assert_eq!(
            data.len(),
            12,
            "NSIDE=1 should have 12 pixels, got {}",
            data.len()
        );

        // The test file has pixels 0, 3, 6, 9 populated with values 1.0, 2.0, 3.0, 4.0
        let expected_values = vec![(0, 1.0), (3, 2.0), (6, 3.0), (9, 4.0)];

        for (idx, expected_val) in expected_values {
            let actual = data.get(idx).expect("pixel out of bounds");
            assert!(
                (actual - expected_val).abs() < 0.01,
                "Pixel {} should be {}, got {}",
                idx,
                expected_val,
                actual
            );
        }

        // Unpopulated pixels should be HPX_UNSEEN
        let unseen_indices = vec![1, 2, 4, 5, 7, 8, 10, 11];
        for idx in unseen_indices {
            let actual = data.get(idx).expect("pixel out of bounds");
            // Use relative tolerance for UNSEEN comparison (healpy may have slight precision differences)
            let diff = (actual - crate::healpix::HPX_UNSEEN).abs();
            let tolerance = (crate::healpix::HPX_UNSEEN.abs() * 1e-6).max(1e-20);
            assert!(
                diff < tolerance,
                "Pixel {} should be HPX_UNSEEN ({}), got {} (diff: {:.2e})",
                idx,
                crate::healpix::HPX_UNSEEN,
                actual,
                diff
            );
        }

        println!("✓ Sparse FITS explicit indexing test passed");
    }

    #[test]
    #[cfg(not(miri))]
    fn test_sparse_map_regression_fix() {
        // Regression test for the v0.5.0 bug where sparse maps were
        // initialized with f64::NEG_INFINITY instead of HPX_UNSEEN.
        //
        // NOTE: This test is skipped under Miri because it uses file system
        // operations (Path::exists) which call statx syscall that Miri's
        // isolation mode doesn't support.
        //
        // The bug manifestation:
        // - Sparse maps initialized with NEG_INFINITY (-inf)
        // - is_seen() filter rejected all pixels (both populated and unpopulated)
        // - Result: "Map contains no valid HEALPix values" panic
        //
        // The fix:
        // - Initialize sparse maps with HPX_UNSEEN (-1.6375e30)
        // - Populated pixels loaded correctly
        // - Unpopulated pixels marked as UNSEEN (filtered out correctly)
        // - is_seen() allows both populated pixels and properly-marked unseen pixels

        let test_file =
            PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures/sparse_nside1.fits");

        if !test_file.exists() {
            eprintln!(
                "Skipping test: sparse FITS fixture not found at {}",
                test_file.display()
            );
            return;
        }

        let data = read_healpix_column(test_file.to_str().unwrap(), 0);

        // Count populated (finite + significant) vs unseen (near HPX_UNSEEN)
        let populated_count = data
            .iter()
            .filter(|&v| v > -1e29) // HPX_UNSEEN is around -1.6375e30
            .count();
        let unseen_count = data
            .iter()
            .filter(|&v| v < -1e29) // Values close to HPX_UNSEEN
            .count();

        // Should have 4 populated distinct values from the test file
        assert_eq!(
            populated_count, 4,
            "Expected 4 populated pixel values, got {}",
            populated_count
        );

        // Should have 8 unseen pixels
        assert_eq!(
            unseen_count, 8,
            "Expected 8 unseen pixels, got {}",
            unseen_count
        );

        // The critical regression check: no NEG_INFINITY values
        let inf_count = data.iter().filter(|v| v == &f64::NEG_INFINITY).count();
        assert_eq!(
            inf_count, 0,
            "No pixels should be NEG_INFINITY, got {}",
            inf_count
        );

        println!("✓ Sparse map regression fix verified");
    }
}