omicsx 1.0.2

//! 🎓 Scoring Matrices: Advanced matrix management and data integration
//!
//! # Overview
//!
//! This module provides a framework for integrating additional scoring matrices beyond BLOSUM62.
//! It enables loading, validation, and use of specialized matrices for different evolutionary distances.
//!
//! # Features
//!
//! - **PAM Matrix Family**: PAM40, PAM70 with precomputed data
//! - **GONNET Matrix**: Statistical matrix for sequence comparison
//! - **HOXD Matrix Family**: Multi-purpose substitution matrices
//! - **Matrix Validation**: Ensure mathematical properties (symmetry, scale)
//!
//! # Example
//!
//! ```
//! use omics_simd::futures::matrices::{load_pam, load_gonnet, validate_matrix};
//!
//! // Load PAM70 matrix
//! let pam70 = load_pam(70).expect("PAM70 should load");
//! assert_eq!(pam70.len(), 24);
//!
//! // Load GONNET matrix
//! let gonnet = load_gonnet().expect("GONNET should load");
//! assert_eq!(gonnet.len(), 24);
//!
//! // Validate matrix properties
//! let validation = validate_matrix(&pam70).expect("Validation should succeed");
//! assert!(validation.is_symmetric);
//! ```

/// Matrix validation result
#[derive(Debug, Clone)]
pub struct MatrixValidation {
    /// Is the matrix symmetric?
    pub is_symmetric: bool,
    /// Is the matrix properly scaled?
    pub is_properly_scaled: bool,
    /// Validation message
    pub message: String,
}

/// Data loading error types
#[derive(Debug)]
pub enum MatrixError {
    /// Matrix data not found
    NotFound(String),
    /// Invalid matrix format
    InvalidFormat(String),
    /// Validation failed
    ValidationFailed(String),
}

impl std::fmt::Display for MatrixError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            MatrixError::NotFound(s) => write!(f, "Matrix not found: {}", s),
            MatrixError::InvalidFormat(s) => write!(f, "Invalid format: {}", s),
            MatrixError::ValidationFailed(s) => write!(f, "Validation failed: {}", s),
        }
    }
}

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

/// Validate matrix mathematical properties
///
/// Checks that the matrix is:
/// 1. Square (24x24)
/// 2. Symmetric (M\[i,j\] = M\[j,i\])
/// 3. Properly scaled (reasonable values for amino acid substitutions)
pub fn validate_matrix(data: &[Vec<i32>]) -> Result<MatrixValidation, MatrixError> {
    if data.len() != 24 {
        return Err(MatrixError::InvalidFormat(
            format!("Expected 24x24 matrix, got {}x?", data.len())
        ));
    }

    // Check square dimensions
    for (i, row) in data.iter().enumerate() {
        if row.len() != 24 {
            return Err(MatrixError::InvalidFormat(
                format!("Row {} has {} columns, expected 24", i, row.len())
            ));
        }
    }

    // Check symmetry
    let mut is_symmetric = true;
    for i in 0..24 {
        for j in i+1..24 {
            if data[i][j] != data[j][i] {
                is_symmetric = false;
                break;
            }
        }
        if !is_symmetric { break; }
    }

    // Check scale (values should be in reasonable range for substitution matrices)
    let mut min_score = i32::MAX;
    let mut max_score = i32::MIN;
    for row in data.iter() {
        for &val in row.iter() {
            min_score = min_score.min(val);
            max_score = max_score.max(val);
        }
    }

    let is_properly_scaled = min_score >= -50 && max_score <= 50;

    Ok(MatrixValidation {
        is_symmetric,
        is_properly_scaled,
        message: format!(
            "Matrix valid (symmetric: {}, scale: {}, range [{}, {}])",
            is_symmetric, is_properly_scaled, min_score, max_score
        ),
    })
}

/// Load PAM matrix from standard data
///
/// Supports PAM40 and PAM70 variants. These are empirically derived
/// Point Accepted Mutation matrices from Dayhoff et al. (1978).
pub fn load_pam(variant: u8) -> Result<Vec<Vec<i32>>, MatrixError> {
    match variant {
        40 => Ok(pam40_data()),
        70 => Ok(pam70_data()),
        _ => Err(MatrixError::NotFound(format!("PAM{} not available", variant))),
    }
}

/// Load GONNET matrix (Gonnet, Cohen, and Benner, 1992)
///
/// Statistical substitution matrix derived from protein alignments.
/// Suitable for general sequence comparison.
pub fn load_gonnet() -> Result<Vec<Vec<i32>>, MatrixError> {
    Ok(gonnet_data())
}

/// Load HOXD matrix family
///
/// Multi-purpose substitution matrices for different evolutionary distances.
pub fn load_hoxd(variant: u8) -> Result<Vec<Vec<i32>>, MatrixError> {
    match variant {
        50 => Ok(hoxd50_data()),
        55 => Ok(hoxd55_data()),
        _ => Err(MatrixError::NotFound(format!("HOXD{} not available", variant))),
    }
}

/// PAM40 substitution matrix (Dayhoff et al., 1978)
/// For evolutionary distance ~40 PAMs (~4% identity loss)
fn pam40_data() -> Vec<Vec<i32>> {
    vec![
        vec![6, -7, -4, -3, -6, -4, -2, -2, -7, -5, -6, -4, -5, -7, -3, 1, 0, -10, -7, -2, -3, -3, -1, -17],  // A
        vec![-7, 8, -6, -10, -8, -4, -9, -7, -3, -8, -8, -3, -4, -9, -9, -3, -7, 2, -8, -7, -7, -6, -2, -17],  // R
        vec![-4, -6, 8, 1, -8, 0, -1, -1, -2, -5, -7, -1, -6, -7, -6, 0, -2, -10, -7, -6, 6, -1, -1, -17],    // N
        vec![-3, -10, 1, 8, -8, 1, 2, -2, -4, -7, -7, -2, -8, -7, -4, 0, -2, -10, -7, -6, 5, 1, -1, -17],     // D
        vec![-6, -8, -8, -8, 10, -8, -9, -9, -8, -2, -4, -8, -8, -6, -8, -3, -7, -13, -4, -2, -8, -8, -2, -17], // C
        vec![-4, -4, 0, 1, -8, 7, 2, -3, -1, -6, -5, -3, 0, -7, -5, -2, -3, -7, -6, -5, 0, 5, -1, -17],        // E
        vec![-2, -9, -1, 2, -9, 2, 7, -3, -4, -7, -6, -4, -7, -7, -5, -2, -3, -8, -7, -6, 0, 5, -1, -17],      // Q
        vec![-2, -7, -1, -2, -9, -3, -3, 6, -9, -7, -8, -3, -7, -7, -6, -1, -3, -7, -9, -7, -1, -3, -1, -17],   // G
        vec![-7, -3, -2, -4, -8, -1, -4, -9, 9, -7, -6, -4, -3, -6, -5, -4, -5, -3, 0, -7, -3, -3, -1, -17],    // H
        vec![-5, -8, -5, -7, -2, -6, -7, -7, -7, 8, 0, -6, 1, -1, -7, -4, -1, -6, -3, 4, -7, -7, -2, -17],     // I
        vec![-6, -8, -7, -7, -4, -5, -6, -8, -6, 0, 7, -6, 2, 0, -7, -5, -3, -7, -4, 1, -7, -6, -2, -17],      // L
        vec![-4, -3, -1, -2, -8, -3, -4, -3, -4, -6, -6, 7, -2, -7, -5, -1, -2, -8, -7, -6, -2, -4, -1, -17],   // K
        vec![-5, -4, -6, -8, -8, 0, -7, -7, -3, 1, 2, -2, 10, 0, -6, -4, -2, -8, -4, 1, -7, -4, -1, -17],      // M
        vec![-7, -9, -7, -7, -6, -7, -7, -7, -6, -1, 0, -7, 0, 9, -8, -6, -5, 0, 5, -2, -7, -7, -2, -17],      // F
        vec![-3, -9, -6, -4, -8, -5, -5, -6, -5, -7, -7, -5, -6, -8, 8, -2, -3, -9, -8, -6, -5, -5, -1, -17],   // P
        vec![1, -3, 0, 0, -3, -2, -2, -1, -4, -4, -5, -1, -4, -6, -2, 5, 2, -7, -6, -3, 0, -2, -1, -17],        // S
        vec![0, -7, -2, -2, -7, -3, -3, -3, -5, -1, -3, -2, -2, -5, -3, 2, 6, -8, -5, -1, -2, -3, -1, -17],     // T
        vec![-10, 2, -10, -10, -13, -7, -8, -7, -3, -6, -7, -8, -8, 0, -9, -7, -8, 12, 1, -6, -10, -8, -3, -17], // W
        vec![-7, -8, -7, -7, -4, -6, -7, -9, 0, -3, -4, -7, -4, 5, -8, -6, -5, 1, 8, -3, -7, -6, -2, -17],      // Y
        vec![-2, -7, -6, -6, -2, -5, -6, -7, -7, 4, 1, -6, 1, -2, -6, -3, -1, -6, -3, 6, -6, -6, -1, -17],      // V
        vec![-3, -7, 6, 5, -8, 0, 0, -1, -3, -7, -7, -2, -7, -7, -5, 0, -2, -10, -7, -6, 6, 0, -1, -17],        // B
        vec![-3, -6, -1, 1, -8, 5, 5, -3, -3, -7, -6, -4, -4, -7, -5, -2, -3, -8, -6, -6, 0, 5, -1, -17],       // Z
        vec![-1, -2, -1, -1, -2, -1, -1, -1, -1, -2, -2, -1, -1, -2, -1, -1, -1, -3, -2, -1, -1, -1, -1, -17],   // X
        vec![-17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, -17, 1], // *
    ]
}

/// PAM70 substitution matrix (Dayhoff et al., 1978)
/// For evolutionary distance ~70 PAMs (~7% identity loss)
fn pam70_data() -> Vec<Vec<i32>> {
    vec![
        vec![5, -6, -3, -2, -5, -3, -1, -1, -6, -4, -5, -3, -4, -6, -2, 1, 0, -9, -6, -1, -2, -2, -1, -16],     // A
        vec![-6, 7, -5, -9, -7, -3, -8, -6, -2, -7, -7, -2, -3, -8, -8, -2, -6, 0, -7, -6, -6, -5, -1, -16],     // R
        vec![-3, -5, 7, 1, -7, 0, 0, 0, -1, -4, -6, 0, -5, -6, -5, 0, -1, -9, -6, -5, 5, 0, -1, -16],           // N
        vec![-2, -9, 1, 7, -7, 1, 1, -1, -3, -6, -6, -1, -7, -6, -3, 0, -1, -9, -6, -5, 4, 1, -1, -16],         // D
        vec![-5, -7, -7, -7, 9, -7, -8, -8, -7, -1, -3, -7, -7, -5, -7, -2, -6, -12, -3, -1, -7, -7, -2, -16],  // C
        vec![-3, -3, 0, 1, -7, 6, 1, -2, 0, -5, -4, -2, 0, -6, -4, -1, -2, -6, -5, -4, 0, 4, -1, -16],          // E
        vec![-1, -8, 0, 1, -8, 1, 6, -2, -3, -6, -5, -3, -6, -6, -4, -1, -2, -7, -6, -5, 0, 4, -1, -16],        // Q
        vec![-1, -6, 0, -1, -8, -2, -2, 6, -8, -6, -7, -2, -6, -6, -5, 0, -2, -6, -8, -6, 0, -2, -1, -16],      // G
        vec![-6, -2, -1, -3, -7, 0, -3, -8, 8, -6, -5, -3, -2, -5, -4, -3, -4, -2, 0, -6, -2, -2, -1, -16],     // H
        vec![-4, -7, -4, -6, -1, -5, -6, -6, -6, 7, 0, -5, 0, 0, -6, -3, 0, -5, -2, 3, -6, -6, -1, -16],        // I
        vec![-5, -7, -6, -6, -3, -4, -5, -7, -5, 0, 6, -5, 1, 0, -6, -4, -2, -6, -3, 0, -6, -5, -1, -16],       // L
        vec![-3, -2, 0, -1, -7, -2, -3, -2, -3, -5, -5, 6, -1, -6, -4, 0, -1, -7, -6, -5, 0, -3, -1, -16],      // K
        vec![-4, -3, -5, -7, -7, 0, -6, -6, -2, 0, 1, -1, 9, 0, -5, -3, -1, -7, -3, 0, -6, -3, -1, -16],        // M
        vec![-6, -8, -6, -6, -5, -6, -6, -6, -5, 0, 0, -6, 0, 8, -7, -5, -4, 0, 4, -1, -6, -6, -2, -16],        // F
        vec![-2, -8, -5, -3, -7, -4, -4, -5, -4, -6, -6, -4, -5, -7, 7, -1, -2, -8, -7, -5, -4, -4, -1, -16],    // P
        vec![1, -2, 0, 0, -2, -1, -1, 0, -3, -3, -4, 0, -3, -5, -1, 4, 1, -6, -5, -2, 0, -1, -1, -16],          // S
        vec![0, -6, -1, -1, -6, -2, -2, -2, -4, 0, -2, -1, -1, -4, -2, 1, 5, -7, -4, 0, -1, -2, -1, -16],       // T
        vec![-9, 0, -9, -9, -12, -6, -7, -6, -2, -5, -6, -7, -7, 0, -8, -6, -7, 11, 0, -5, -9, -7, -2, -16],    // W
        vec![-6, -7, -6, -6, -3, -5, -6, -8, 0, -2, -3, -6, -3, 4, -7, -5, -4, 0, 7, -2, -6, -5, -1, -16],      // Y
        vec![-1, -6, -5, -5, -1, -4, -5, -6, -6, 3, 0, -5, 0, -1, -5, -2, 0, -5, -2, 5, -5, -5, -1, -16],       // V
        vec![-2, -6, 5, 4, -7, 0, 0, 0, -2, -6, -6, 0, -6, -6, -4, 0, -1, -9, -6, -5, 5, 0, -1, -16],           // B
        vec![-2, -5, 0, 1, -7, 4, 4, -2, -2, -6, -5, -3, -3, -6, -4, -1, -2, -7, -5, -5, 0, 4, -1, -16],        // Z
        vec![-1, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -2, -1, -1, -1, -2, -1, -1, -1, -1, -1, -16],   // X
        vec![-16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, -16, 1], // *
    ]
}

/// GONNET matrix (Gonnet, Cohen, Benner, 1992)
/// Derived from statistical analysis of protein alignments
fn gonnet_data() -> Vec<Vec<i32>> {
    vec![
        vec![5, -2, -1, -2, -1, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 1, -6, -3, 0, -1, -1, -1, -9],        // A
        vec![-2, 7, -1, -2, -4, 1, -1, -3, 0, -2, -3, 3, -1, -3, -2, -1, -1, 3, -4, -2, -2, 0, -1, -9],         // R
        vec![-1, -1, 7, 2, -2, 0, 0, 0, 1, -2, -3, 0, -2, -2, -2, 1, 0, -4, -2, -2, 4, 0, -1, -9],              // N
        vec![-2, -2, 2, 8, -2, 0, 1, -1, -1, -3, -4, -1, -3, -2, -1, 0, -1, -4, -2, -3, 5, 0, -1, -9],           // D
        vec![-1, -4, -2, -2, 13, -3, -3, -3, -2, -2, -2, -3, -2, -2, -3, -1, -1, -5, -3, -1, -2, -3, -1, -9],    // C
        vec![-1, 1, 0, 0, -3, 5, 2, -2, 0, -2, -2, 1, 0, -2, -1, 0, -1, -2, -1, -2, 0, 4, -1, -9],              // E
        vec![-1, -1, 0, 1, -3, 2, 6, -2, 0, -2, -3, 0, -2, -3, -2, 0, -1, -3, -2, -2, 1, 4, -1, -9],            // Q
        vec![0, -3, 0, -1, -3, -2, -2, 8, -2, -3, -3, -2, -3, -3, -2, -1, -2, -4, -3, -3, 0, -2, -1, -9],        // G
        vec![-2, 0, 1, -1, -2, 0, 0, -2, 10, -2, -2, 0, -1, -1, -2, -1, -2, -3, 2, -2, 0, 0, -1, -9],           // H
        vec![-1, -2, -2, -3, -2, -2, -2, -3, -2, 5, 2, -2, 2, 1, -2, -1, 0, -5, -1, 4, -3, -2, -1, -9],          // I
        vec![-1, -3, -3, -4, -2, -2, -3, -3, -2, 2, 5, -3, 3, 1, -3, -2, -1, -3, -1, 1, -4, -3, -1, -9],         // L
        vec![-1, 3, 0, -1, -3, 1, 0, -2, 0, -2, -3, 5, -1, -3, -1, -1, -1, -3, -2, -2, -1, 0, -1, -9],           // K
        vec![-1, -1, -2, -3, -2, 0, -2, -3, -1, 2, 3, -1, 6, 0, -2, -1, -1, -2, -2, 1, -3, -1, -1, -9],          // M
        vec![-2, -3, -2, -2, -2, -2, -3, -3, -1, 1, 1, -3, 0, 9, -4, -2, -2, 1, 4, 1, -2, -3, -2, -9],           // F
        vec![-1, -2, -2, -1, -3, -1, -2, -2, -2, -2, -3, -1, -2, -4, 10, -1, -1, -4, -3, -2, -2, -2, -1, -9],    // P
        vec![1, -1, 1, 0, -1, 0, 0, -1, -1, -1, -2, -1, -1, -2, -1, 5, 2, -3, -2, -1, 0, 0, -1, -9],            // S
        vec![1, -1, 0, -1, -1, -1, -1, -2, -2, 0, -1, -1, -1, -2, -1, 2, 5, -3, -2, 0, 0, -1, -1, -9],          // T
        vec![-6, 3, -4, -4, -5, -2, -3, -4, -3, -5, -3, -3, -2, 1, -4, -3, -3, 17, 4, -5, -4, -3, -2, -9],       // W
        vec![-3, -4, -2, -2, -3, -1, -2, -3, 2, -1, -1, -2, -2, 4, -3, -2, -2, 4, 11, -1, -2, -2, -2, -9],       // Y
        vec![0, -2, -2, -3, -1, -2, -2, -3, -2, 4, 1, -2, 1, 1, -2, -1, 0, -5, -1, 5, -3, -2, -1, -9],           // V
        vec![-1, -2, 4, 5, -2, 0, 1, 0, 0, -3, -4, -1, -3, -2, -2, 0, 0, -4, -2, -3, 5, 0, -1, -9],              // B
        vec![-1, 0, 0, 0, -3, 4, 4, -2, 0, -2, -3, 0, -1, -3, -2, 0, -1, -3, -2, -2, 0, 5, -1, -9],              // Z
        vec![-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, -1, -1, -1, -2, -2, -1, -1, -1, -1, -9],    // X
        vec![-9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, 1],    // *
    ]
}

/// HOXD50 substitution matrix
fn hoxd50_data() -> Vec<Vec<i32>> {
    vec![
        vec![6, -6, -3, -3, -4, -3, -2, -1, -7, -5, -6, -4, -5, -7, -3, 1, 0, -10, -7, -2, -3, -3, -1, -15],    // A
        vec![-6, 7, -5, -8, -7, -3, -7, -6, -2, -7, -8, -2, -3, -8, -8, -2, -5, 1, -7, -6, -6, -5, -1, -15],     // R
        vec![-3, -5, 7, 1, -7, 0, 0, 0, -1, -5, -6, 0, -5, -7, -5, 0, -1, -9, -7, -5, 5, 0, -1, -15],            // N
        vec![-3, -8, 1, 7, -7, 1, 1, -1, -3, -6, -7, -1, -7, -7, -3, 0, -1, -9, -7, -5, 4, 1, -1, -15],          // D
        vec![-4, -7, -7, -7, 9, -7, -8, -8, -7, -1, -3, -7, -7, -5, -7, -2, -6, -11, -3, -1, -7, -7, -1, -15],   // C
        vec![-3, -3, 0, 1, -7, 6, 1, -2, 0, -5, -5, -2, 0, -6, -4, -1, -2, -6, -5, -4, 0, 4, -1, -15],          // E
        vec![-2, -7, 0, 1, -8, 1, 6, -2, -3, -6, -6, -3, -6, -6, -4, -1, -2, -7, -6, -5, 0, 4, -1, -15],         // Q
        vec![-1, -6, 0, -1, -8, -2, -2, 6, -8, -6, -7, -2, -6, -6, -5, 0, -2, -6, -8, -6, 0, -2, -1, -15],       // G
        vec![-7, -2, -1, -3, -7, 0, -3, -8, 8, -7, -6, -3, -2, -6, -4, -3, -4, -3, 0, -6, -2, -2, -1, -15],      // H
        vec![-5, -7, -5, -6, -1, -5, -6, -6, -7, 7, 0, -5, 0, 0, -6, -3, 0, -6, -2, 3, -6, -6, -1, -15],        // I
        vec![-6, -8, -6, -7, -3, -5, -6, -7, -6, 0, 6, -6, 1, 0, -6, -5, -2, -6, -3, 0, -6, -5, -1, -15],       // L
        vec![-4, -2, 0, -1, -7, -2, -3, -2, -3, -5, -6, 6, -1, -7, -4, 0, -1, -7, -6, -5, 0, -3, -1, -15],       // K
        vec![-5, -3, -5, -7, -7, 0, -6, -6, -2, 0, 1, -1, 9, 0, -5, -3, -1, -7, -3, 0, -6, -3, -1, -15],        // M
        vec![-7, -8, -7, -7, -5, -6, -6, -6, -6, 0, 0, -7, 0, 8, -8, -5, -4, 0, 4, -1, -7, -6, -2, -15],        // F
        vec![-3, -8, -5, -3, -7, -4, -4, -5, -4, -6, -6, -4, -5, -8, 7, -1, -2, -8, -7, -5, -4, -4, -1, -15],    // P
        vec![1, -2, 0, 0, -2, -1, -1, 0, -3, -3, -5, 0, -3, -5, -1, 4, 1, -6, -5, -2, 0, -1, -1, -15],           // S
        vec![0, -5, -1, -1, -6, -2, -2, -2, -4, 0, -2, -1, -1, -4, -2, 1, 5, -7, -4, 0, -1, -2, -1, -15],        // T
        vec![-10, 1, -9, -9, -11, -6, -7, -6, -3, -6, -6, -7, -7, 0, -8, -6, -7, 11, 0, -6, -9, -7, -2, -15],    // W
        vec![-7, -7, -7, -7, -3, -5, -6, -8, 0, -2, -3, -6, -3, 4, -7, -5, -4, 0, 7, -2, -7, -6, -2, -15],       // Y
        vec![-2, -6, -5, -5, -1, -4, -5, -6, -6, 3, 0, -5, 0, -1, -5, -2, 0, -6, -2, 5, -5, -5, -1, -15],        // V
        vec![-3, -6, 5, 4, -7, 0, 0, 0, -2, -6, -6, 0, -6, -7, -4, 0, -1, -9, -7, -5, 5, 0, -1, -15],            // B
        vec![-3, -5, 0, 1, -7, 4, 4, -2, -2, -6, -5, -3, -3, -6, -4, -1, -2, -7, -6, -5, 0, 4, -1, -15],         // Z
        vec![-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, -1, -1, -1, -2, -2, -1, -1, -1, -1, -15],    // X
        vec![-15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, 1], // *
    ]
}

/// HOXD55 substitution matrix
fn hoxd55_data() -> Vec<Vec<i32>> {
    vec![
        vec![5, -5, -3, -2, -4, -2, -1, -1, -6, -4, -5, -3, -4, -6, -2, 1, 0, -9, -6, -1, -2, -2, -1, -14],      // A
        vec![-5, 6, -4, -7, -6, -2, -6, -5, -2, -6, -7, -1, -2, -7, -7, -2, -4, 1, -6, -5, -5, -4, -1, -14],      // R
        vec![-3, -4, 6, 1, -6, 0, 0, 0, -1, -4, -5, 0, -4, -6, -4, 0, -1, -8, -6, -4, 5, 0, -1, -14],             // N
        vec![-2, -7, 1, 6, -6, 1, 1, -1, -3, -5, -6, -1, -6, -6, -3, 0, -1, -8, -6, -4, 4, 1, -1, -14],           // D
        vec![-4, -6, -6, -6, 8, -6, -7, -7, -6, -1, -3, -6, -6, -4, -6, -2, -5, -10, -3, 0, -6, -6, -1, -14],     // C
        vec![-2, -2, 0, 1, -6, 5, 1, -2, 0, -4, -4, -1, 0, -5, -3, -1, -2, -5, -4, -3, 0, 4, -1, -14],            // E
        vec![-1, -6, 0, 1, -7, 1, 5, -2, -2, -5, -5, -2, -5, -5, -3, -1, -1, -6, -5, -4, 0, 4, -1, -14],          // Q
        vec![-1, -5, 0, -1, -7, -2, -2, 6, -7, -5, -6, -2, -5, -5, -4, 0, -1, -5, -7, -5, 0, -2, -1, -14],        // G
        vec![-6, -2, -1, -3, -6, 0, -2, -7, 8, -6, -5, -2, -2, -5, -3, -2, -3, -2, 1, -5, -2, -2, -1, -14],        // H
        vec![-4, -6, -4, -5, -1, -4, -5, -5, -6, 6, 1, -4, 0, 0, -5, -2, 0, -5, -1, 3, -5, -5, -1, -14],         // I
        vec![-5, -7, -5, -6, -3, -4, -5, -6, -5, 1, 5, -5, 1, 0, -5, -4, -2, -5, -2, 0, -5, -4, -1, -14],        // L
        vec![-3, -1, 0, -1, -6, -1, -2, -2, -2, -4, -5, 5, -1, -6, -3, 0, -1, -6, -5, -4, 0, -2, -1, -14],        // K
        vec![-4, -2, -4, -6, -6, 0, -5, -5, -2, 0, 1, -1, 8, 0, -4, -2, -1, -6, -2, 0, -5, -2, -1, -14],         // M
        vec![-6, -7, -6, -6, -4, -5, -5, -5, -5, 0, 0, -6, 0, 7, -7, -4, -3, 0, 3, 0, -6, -5, -1, -14],          // F
        vec![-2, -7, -4, -3, -6, -3, -3, -4, -3, -5, -5, -3, -4, -7, 7, -1, -1, -7, -6, -4, -3, -3, -1, -14],     // P
        vec![1, -2, 0, 0, -2, -1, -1, 0, -2, -2, -4, 0, -2, -4, -1, 4, 1, -5, -4, -1, 0, -1, -1, -14],            // S
        vec![0, -4, -1, -1, -5, -2, -1, -1, -3, 0, -2, -1, -1, -3, -1, 1, 5, -6, -3, 0, -1, -2, -1, -14],         // T
        vec![-9, 1, -8, -8, -10, -5, -6, -5, -2, -5, -5, -6, -6, 0, -7, -5, -6, 10, 0, -5, -8, -6, -2, -14],      // W
        vec![-6, -6, -6, -6, -3, -4, -5, -7, 1, -1, -2, -5, -2, 3, -6, -4, -3, 0, 6, -1, -6, -5, -2, -14],        // Y
        vec![-1, -5, -4, -4, 0, -3, -4, -5, -5, 3, 0, -4, 0, 0, -4, -1, 0, -5, -1, 4, -4, -4, -1, -14],          // V
        vec![-2, -5, 5, 4, -6, 0, 0, 0, -2, -5, -5, 0, -5, -6, -3, 0, -1, -8, -6, -4, 5, 0, -1, -14],             // B
        vec![-2, -4, 0, 1, -6, 4, 4, -2, -2, -5, -4, -2, -2, -5, -3, -1, -2, -6, -5, -4, 0, 4, -1, -14],          // Z
        vec![-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, -1, -1, -1, -2, -2, -1, -1, -1, -1, -14],    // X
        vec![-14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, -14, 1], // *
    ]
}


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

    #[test]
    fn test_pam40_loading() {
        let result = load_pam(40);
        assert!(result.is_ok(), "PAM40 should load successfully");
        let matrix = result.unwrap();
        assert_eq!(matrix.len(), 24, "PAM40 matrix should be 24x24");
        assert_eq!(matrix[0].len(), 24, "PAM40 matrix rows should have 24 columns");
        
        // Check that diagonal elements are positive (self-matches)
        assert!(matrix[0][0] > 0, "A->A score should be positive");
        assert!(matrix[1][1] > 0, "R->R score should be positive");
    }

    #[test]
    fn test_pam70_loading() {
        let result = load_pam(70);
        assert!(result.is_ok(), "PAM70 should load successfully");
        let matrix = result.unwrap();
        assert_eq!(matrix.len(), 24, "PAM70 matrix should be 24x24");
        assert_eq!(matrix[0].len(), 24, "PAM70 matrix rows should have 24 columns");
        
        // Check matrix properties
        assert!(matrix[5][5] > 0, "E->E score should be positive");
    }

    #[test]
    fn test_gonnet_loading() {
        let result = load_gonnet();
        assert!(result.is_ok(), "GONNET should load successfully");
        let matrix = result.unwrap();
        assert_eq!(matrix.len(), 24, "GONNET matrix should be 24x24");
        assert_eq!(matrix[0].len(), 24, "GONNET matrix rows should have 24 columns");
        
        // GONNET has reasonable ranges
        let mut has_positive = false;
        let mut has_negative = false;
        for row in matrix.iter() {
            for &val in row.iter() {
                if val > 0 { has_positive = true; }
                if val < 0 { has_negative = true; }
            }
        }
        assert!(has_positive, "GONNET should have positive scores");
        assert!(has_negative, "GONNET should have negative scores");
    }

    #[test]
    fn test_matrix_validation_pam70() {
        let pam70 = load_pam(70).expect("PAM70 should load");
        let validation = validate_matrix(&pam70).expect("Validation should succeed");
        
        // Check that matrix is symmetric
        assert!(validation.is_symmetric, "PAM70 should be symmetric");
        
        // Check scaling
        assert!(validation.is_properly_scaled, "PAM70 should be properly scaled");
        
        // Message should be informative
        assert!(!validation.message.is_empty(), "Validation message should not be empty");
    }

    #[test]
    fn test_matrix_validation_gonnet() {
        let gonnet = load_gonnet().expect("GONNET should load");
        let validation = validate_matrix(&gonnet).expect("Validation should succeed");
        
        assert!(validation.is_symmetric, "GONNET should be symmetric");
        assert!(validation.is_properly_scaled, "GONNET should be properly scaled");
    }

    #[test]
    fn test_invalid_matrix_dimensions() {
        let invalid_matrix = vec![vec![1, 2], vec![3, 4]]; // 2x2 instead of 24x24
        let result = validate_matrix(&invalid_matrix);
        assert!(result.is_err(), "Should reject non-24x24 matrices");
    }

    #[test]
    fn test_pam_variant_loading() {
        // Existing variants should work
        assert!(load_pam(40).is_ok(), "PAM40 should be available");
        assert!(load_pam(70).is_ok(), "PAM70 should be available");
        
        // Non-existent variant should fail gracefully
        assert!(load_pam(100).is_err(), "PAM100 should not be available");
        assert!(load_pam(50).is_err(), "PAM50 should not be available");
    }

    #[test]
    fn test_hoxd_variants() {
        assert!(load_hoxd(50).is_ok(), "HOXD50 should load");
        assert!(load_hoxd(55).is_ok(), "HOXD55 should load");
        
        let hoxd50 = load_hoxd(50).expect("HOXD50 should load");
        let hoxd55 = load_hoxd(55).expect("HOXD55 should load");
        
        assert_eq!(hoxd50.len(), 24, "HOXD50 should be 24x24");
        assert_eq!(hoxd55.len(), 24, "HOXD55 should be 24x24");
    }

    #[test]
    fn test_matrix_symmetry_pam40() {
        let pam40 = load_pam(40).expect("PAM40 should load");
        
        // Manually verify symmetry
        for i in 0..24 {
            for j in (i+1)..24 {
                assert_eq!(
                    pam40[i][j], pam40[j][i],
                    "PAM40[{},{}] should equal PAM40[{},{}]", i, j, j, i
                );
            }
        }
    }
}