fastars 0.1.0

//! Filter criteria definitions.
//!
//! This module defines the criteria used for filtering reads,
//! including quality, length, N-content, and complexity filters.

use crate::trim::LengthConfig;
use super::index::IndexFilterConfig;

/// Decision made by the filter.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FilterDecision {
    /// Read passes all filters.
    Pass,
    /// Read failed quality filter.
    FailQuality,
    /// Read failed length filter.
    FailLength,
    /// Read failed complexity filter.
    FailComplexity,
    /// Read failed N-rate filter.
    FailNRate,
    /// Read failed index barcode filter.
    FailIndex,
}

impl FilterDecision {
    /// Check if the decision is a pass.
    #[inline]
    pub fn is_pass(&self) -> bool {
        matches!(self, FilterDecision::Pass)
    }

    /// Check if the decision is a failure.
    #[inline]
    pub fn is_fail(&self) -> bool {
        !self.is_pass()
    }
}

/// Configuration for read filtering.
#[derive(Debug, Clone)]
pub struct FilterConfig {
    /// Minimum Phred quality for a base to be considered "qualified" (fastp's qualified_quality_phred).
    /// Default: 15 (Q15).
    pub qualified_quality: u8,
    /// Maximum percentage of unqualified (low-quality) bases allowed (fastp's unqualified_percent_limit).
    /// Default: 40.0 (40%).
    pub unqualified_percent_limit: f64,
    /// Maximum number of N bases allowed (fastp's n_base_limit).
    /// Default: 5.
    pub n_base_limit: usize,
    /// Maximum N content as percentage (long mode only, 0.0-100.0).
    pub n_percent_limit: Option<f64>,
    /// Minimum average quality score to pass (optional, fastp default OFF).
    pub min_avg_quality: Option<u8>,
    /// Maximum rate of N bases (0.0-1.0) - optional, for backward compatibility.
    pub max_n_rate: Option<f64>,
    /// Low complexity threshold (0.0-1.0, lower = more complex).
    pub low_complexity_threshold: Option<f64>,
    /// Length filtering configuration.
    pub length_config: LengthConfig,
    /// Whether quality filtering is enabled.
    pub quality_filtering_enabled: bool,
    /// Index barcode filtering configuration.
    pub index_filter: IndexFilterConfig,
}

impl Default for FilterConfig {
    fn default() -> Self {
        Self {
            qualified_quality: 15,
            unqualified_percent_limit: 40.0,
            n_base_limit: 5,
            n_percent_limit: None,
            min_avg_quality: None,
            max_n_rate: None,
            low_complexity_threshold: None,
            length_config: LengthConfig::default(),
            quality_filtering_enabled: true,
            index_filter: IndexFilterConfig::default(),
        }
    }
}

impl FilterConfig {
    /// Create a new filter config with default settings.
    pub fn new() -> Self {
        Self::default()
    }

    /// Create config for short read mode.
    /// Matches fastp defaults: Q15 threshold, 40% unqualified limit, 5 N bases max.
    pub fn short_read() -> Self {
        Self {
            qualified_quality: 15,
            unqualified_percent_limit: 40.0,
            n_base_limit: 5,
            n_percent_limit: None, // Not used in short mode
            min_avg_quality: None, // fastp default OFF
            max_n_rate: None,      // use n_base_limit instead
            low_complexity_threshold: None, // fastp default OFF
            length_config: LengthConfig::short_read(),
            quality_filtering_enabled: true,
            index_filter: IndexFilterConfig::default(),
        }
    }

    /// Create config for long read mode.
    pub fn long_read() -> Self {
        Self {
            qualified_quality: 7, // Lower threshold for long reads
            unqualified_percent_limit: 50.0, // More lenient for long reads
            n_base_limit: 10, // Allow more Ns for long reads
            n_percent_limit: Some(10.0), // Default 10% for long reads
            min_avg_quality: None,
            max_n_rate: None,
            low_complexity_threshold: None, // Complexity less relevant for long reads
            length_config: LengthConfig::long_read(),
            quality_filtering_enabled: true,
            index_filter: IndexFilterConfig::default(),
        }
    }

    /// Set the qualified quality threshold (fastp's -q).
    pub fn with_qualified_quality(mut self, quality: u8) -> Self {
        self.qualified_quality = quality;
        self
    }

    /// Set the unqualified percent limit (fastp's -u).
    pub fn with_unqualified_percent_limit(mut self, limit: f64) -> Self {
        self.unqualified_percent_limit = limit.clamp(0.0, 100.0);
        self
    }

    /// Set the N base limit (fastp's -n).
    pub fn with_n_base_limit(mut self, limit: usize) -> Self {
        self.n_base_limit = limit;
        self
    }

    /// Set the N percent limit (long mode only).
    pub fn with_n_percent_limit(mut self, percent: f64) -> Self {
        self.n_percent_limit = Some(percent.clamp(0.0, 100.0));
        self
    }

    /// Set minimum average quality (fastp's -e).
    pub fn with_min_avg_quality(mut self, quality: u8) -> Self {
        self.min_avg_quality = Some(quality);
        self
    }

    /// Disable all quality filtering.
    pub fn without_quality_filter(mut self) -> Self {
        self.quality_filtering_enabled = false;
        self.min_avg_quality = None;
        self
    }

    /// Enable quality filtering.
    pub fn with_quality_filter(mut self) -> Self {
        self.quality_filtering_enabled = true;
        self
    }

    /// Set maximum N rate (legacy, for backward compatibility).
    pub fn with_max_n_rate(mut self, rate: f64) -> Self {
        self.max_n_rate = Some(rate.clamp(0.0, 1.0));
        self
    }

    /// Set complexity threshold.
    pub fn with_complexity_threshold(mut self, threshold: f64) -> Self {
        self.low_complexity_threshold = Some(threshold.clamp(0.0, 1.0));
        self
    }

    /// Disable complexity filter.
    pub fn without_complexity_filter(mut self) -> Self {
        self.low_complexity_threshold = None;
        self
    }

    /// Set length configuration.
    pub fn with_length_config(mut self, config: LengthConfig) -> Self {
        self.length_config = config;
        self
    }

    /// Set minimum length.
    pub fn with_min_length(mut self, length: usize) -> Self {
        self.length_config.min_length = length;
        self
    }

    /// Set maximum length.
    pub fn with_max_length(mut self, length: usize) -> Self {
        self.length_config.max_length = Some(length);
        self
    }

    /// Set index 1 barcode filter.
    pub fn with_index1_filter(mut self, barcode: String, max_mismatches: usize) -> Self {
        self.index_filter = self.index_filter.with_index1(barcode, max_mismatches);
        self
    }

    /// Set index 2 barcode filter.
    pub fn with_index2_filter(mut self, barcode: String, max_mismatches: usize) -> Self {
        self.index_filter = self.index_filter.with_index2(barcode, max_mismatches);
        self
    }
}

/// Apply all filters to a read.
///
/// Filter order matches fastp priority:
/// 1. Length check (cheapest)
/// 2. Unqualified base percentage (fastp's primary quality filter)
/// 3. Average quality (optional)
/// 4. N base count (fastp style, absolute count)
/// 5. N rate (optional, legacy/backward compatibility)
/// 6. Complexity (optional, most expensive)
///
/// # Arguments
/// * `name` - The read name (for index filtering)
/// * `seq` - The read sequence
/// * `qual` - Quality scores (Phred+33 encoded)
/// * `config` - Filter configuration
///
/// # Returns
/// The filter decision indicating pass or failure reason.
pub fn apply_filters(name: &[u8], seq: &[u8], qual: &[u8], config: &FilterConfig) -> FilterDecision {
    // 0. Check index barcode filter first (very cheap)
    if config.index_filter.is_enabled() {
        if !super::index::check_index_filter(name, &config.index_filter) {
            return FilterDecision::FailIndex;
        }
    }

    // 1. Check length (cheapest)
    if !check_length(seq.len(), config) {
        return FilterDecision::FailLength;
    }

    // 2. Unqualified base percentage (fastp's primary quality filter)
    if config.quality_filtering_enabled
        && !check_unqualified_percent(qual, config.qualified_quality, config.unqualified_percent_limit)
    {
        return FilterDecision::FailQuality;
    }

    // 3. Average quality (optional)
    if let Some(min_qual) = config.min_avg_quality {
        if !check_avg_quality(qual, min_qual) {
            return FilterDecision::FailQuality;
        }
    }

    // 4. N base count (fastp style, absolute count)
    if !check_n_base_count(seq, config.n_base_limit) {
        return FilterDecision::FailNRate;
    }

    // 4b. N percent limit (long mode, percentage-based)
    if let Some(percent_limit) = config.n_percent_limit {
        if !check_n_percent(seq, percent_limit) {
            return FilterDecision::FailNRate;
        }
    }

    // 5. N rate (optional, for backward compatibility)
    if let Some(max_rate) = config.max_n_rate {
        if !check_n_rate(seq, max_rate) {
            return FilterDecision::FailNRate;
        }
    }

    // 6. Complexity (optional, most expensive)
    if let Some(threshold) = config.low_complexity_threshold {
        if !check_complexity(seq, threshold) {
            return FilterDecision::FailComplexity;
        }
    }

    FilterDecision::Pass
}

/// Check if a length passes the filter.
#[inline]
fn check_length(len: usize, config: &FilterConfig) -> bool {
    crate::trim::check_length(len, &config.length_config)
}

/// Check if average quality meets the minimum.
///
/// # Arguments
/// * `qual` - Quality scores (Phred+33 encoded)
/// * `min` - Minimum average Phred score
pub fn check_avg_quality(qual: &[u8], min: u8) -> bool {
    if qual.is_empty() {
        return false;
    }

    let sum: u64 = qual.iter().map(|&q| q.saturating_sub(33) as u64).sum();
    let avg = sum / qual.len() as u64;
    avg >= min as u64
}

/// Check if N rate is within acceptable limits.
///
/// # Arguments
/// * `seq` - The read sequence
/// * `max_rate` - Maximum allowed N rate (0.0-1.0)
pub fn check_n_rate(seq: &[u8], max_rate: f64) -> bool {
    if seq.is_empty() {
        return true;
    }

    let n_count = seq.iter().filter(|&&b| b == b'N' || b == b'n').count();
    let rate = n_count as f64 / seq.len() as f64;
    rate <= max_rate
}

/// Check if low-quality base percentage is within limit.
///
/// This is fastp's core quality filter. It counts bases with quality below
/// the qualified_quality threshold and fails if the percentage exceeds the limit.
///
/// # Arguments
/// * `qual` - Quality scores (Phred+33 encoded)
/// * `qualified_quality` - Minimum Phred score for a base to be "qualified"
/// * `limit` - Maximum percentage of unqualified bases allowed (0-100)
///
/// # Returns
/// `true` if the read passes (unqualified % <= limit), `false` otherwise.
#[inline]
pub fn check_unqualified_percent(qual: &[u8], qualified_quality: u8, limit: f64) -> bool {
    if qual.is_empty() {
        return false;
    }

    let threshold = qualified_quality.saturating_add(33); // Phred+33 encoded
    let low_qual_count = qual.iter().filter(|&&q| q < threshold).count();
    let percent = (low_qual_count as f64 / qual.len() as f64) * 100.0;
    percent <= limit
}

/// Check if N base count is within limit (fastp style).
///
/// Unlike check_n_rate which uses a percentage, this uses an absolute count.
///
/// # Arguments
/// * `seq` - The read sequence
/// * `limit` - Maximum number of N bases allowed
///
/// # Returns
/// `true` if the read passes (N count <= limit), `false` otherwise.
#[inline]
pub fn check_n_base_count(seq: &[u8], limit: usize) -> bool {
    let n_count = seq.iter().filter(|&&b| b == b'N' || b == b'n').count();
    n_count <= limit
}

/// Check if N content percentage is within limit (long mode).
///
/// # Arguments
/// * `seq` - The read sequence
/// * `percent_limit` - Maximum N content percentage allowed (0-100)
///
/// # Returns
/// `true` if the read passes (N percent <= limit), `false` otherwise.
#[inline]
pub fn check_n_percent(seq: &[u8], percent_limit: f64) -> bool {
    if seq.is_empty() {
        return true;
    }

    let n_count = seq.iter().filter(|&&b| b == b'N' || b == b'n').count();
    let percent = (n_count as f64 / seq.len() as f64) * 100.0;
    percent <= percent_limit
}

/// Check if sequence complexity is above threshold.
///
/// Uses a simple dinucleotide diversity measure.
/// A sequence with only one type of dinucleotide has complexity 0.
/// A sequence with even distribution of all 16 dinucleotides has complexity 1.
///
/// # Arguments
/// * `seq` - The read sequence
/// * `threshold` - Minimum complexity (0.0-1.0)
pub fn check_complexity(seq: &[u8], threshold: f64) -> bool {
    if seq.len() < 2 {
        return true; // Too short to assess
    }

    // Count dinucleotides
    let mut counts = [0u32; 16];
    for window in seq.windows(2) {
        if let (Some(i), Some(j)) = (base_to_idx(window[0]), base_to_idx(window[1])) {
            counts[i * 4 + j] += 1;
        }
    }

    // Calculate entropy-based complexity
    let total = (seq.len() - 1) as f64;
    let mut entropy = 0.0;

    for &count in &counts {
        if count > 0 {
            let p = count as f64 / total;
            entropy -= p * p.log2();
        }
    }

    // Normalize to 0-1 range (max entropy for 16 categories is 4 bits)
    let complexity = entropy / 4.0;
    complexity >= threshold
}

/// Convert base to index (A=0, C=1, G=2, T=3).
#[inline]
fn base_to_idx(base: u8) -> Option<usize> {
    match base {
        b'A' | b'a' => Some(0),
        b'C' | b'c' => Some(1),
        b'G' | b'g' => Some(2),
        b'T' | b't' => Some(3),
        _ => None,
    }
}

// Legacy type alias for compatibility
pub type FilterCriteria = FilterConfig;

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

    fn make_qual(scores: &[u8]) -> Vec<u8> {
        scores.iter().map(|&s| s + 33).collect()
    }

    #[test]
    fn test_filter_decision_is_pass() {
        assert!(FilterDecision::Pass.is_pass());
        assert!(!FilterDecision::FailQuality.is_pass());
        assert!(!FilterDecision::FailLength.is_pass());
        assert!(!FilterDecision::FailComplexity.is_pass());
        assert!(!FilterDecision::FailNRate.is_pass());
    }

    #[test]
    fn test_filter_decision_is_fail() {
        assert!(!FilterDecision::Pass.is_fail());
        assert!(FilterDecision::FailQuality.is_fail());
    }

    #[test]
    fn test_filter_config_default() {
        let config = FilterConfig::default();
        // fastp defaults
        assert_eq!(config.qualified_quality, 15);
        assert!((config.unqualified_percent_limit - 40.0).abs() < 0.001);
        assert_eq!(config.n_base_limit, 5);
        assert!(config.min_avg_quality.is_none());
        assert!(config.max_n_rate.is_none());
        assert!(config.low_complexity_threshold.is_none());
        assert!(config.quality_filtering_enabled);
    }

    #[test]
    fn test_filter_config_short_read() {
        let config = FilterConfig::short_read();
        // fastp defaults for short reads
        assert_eq!(config.qualified_quality, 15);
        assert!((config.unqualified_percent_limit - 40.0).abs() < 0.001);
        assert_eq!(config.n_base_limit, 5);
        assert!(config.min_avg_quality.is_none()); // fastp default OFF
        assert!(config.low_complexity_threshold.is_none()); // fastp default OFF
        assert!(config.quality_filtering_enabled);
    }

    #[test]
    fn test_filter_config_long_read() {
        let config = FilterConfig::long_read();
        assert_eq!(config.qualified_quality, 7);
        assert!((config.unqualified_percent_limit - 50.0).abs() < 0.001);
        assert_eq!(config.n_base_limit, 10);
        assert!(config.low_complexity_threshold.is_none());
        assert!(config.quality_filtering_enabled);
    }

    #[test]
    fn test_check_avg_quality_pass() {
        let qual = make_qual(&[30, 30, 30, 30, 30]);
        assert!(check_avg_quality(&qual, 20));
        assert!(check_avg_quality(&qual, 30));
    }

    #[test]
    fn test_check_avg_quality_fail() {
        let qual = make_qual(&[10, 10, 10, 10, 10]);
        assert!(!check_avg_quality(&qual, 20));
    }

    #[test]
    fn test_check_avg_quality_empty() {
        assert!(!check_avg_quality(&[], 10));
    }

    #[test]
    fn test_check_n_rate_pass() {
        let seq = b"ACGTACGTACGT";
        assert!(check_n_rate(seq, 0.05));
    }

    #[test]
    fn test_check_n_rate_fail() {
        let seq = b"ACGTNNNNNNNN";
        assert!(!check_n_rate(seq, 0.05));
    }

    #[test]
    fn test_check_n_rate_boundary() {
        let seq = b"ACGTACGTNN"; // 2/10 = 0.2
        assert!(!check_n_rate(seq, 0.1));
        assert!(check_n_rate(seq, 0.2));
        assert!(check_n_rate(seq, 0.3));
    }

    #[test]
    fn test_check_n_rate_empty() {
        assert!(check_n_rate(&[], 0.05));
    }

    #[test]
    fn test_check_complexity_pass() {
        // Good complexity - mix of bases
        let seq = b"ACGTACGTACGTACGT";
        assert!(check_complexity(seq, 0.3));
    }

    #[test]
    fn test_check_complexity_fail() {
        // Low complexity - poly-A
        let seq = b"AAAAAAAAAAAAAAAA";
        assert!(!check_complexity(seq, 0.3));
    }

    #[test]
    fn test_check_complexity_short() {
        // Very short sequences pass
        assert!(check_complexity(b"A", 0.5));
        assert!(check_complexity(b"", 0.5));
    }

    #[test]
    fn test_apply_filters_pass() {
        let name = b"@test_read";
        let seq = b"ACGTACGTACGTACGTACGTACGTACGTACGT"; // 32bp
        let qual = make_qual(&[30; 32]);
        let config = FilterConfig::short_read();
        let decision = apply_filters(name, seq, &qual, &config);
        assert!(decision.is_pass());
    }

    #[test]
    fn test_apply_filters_fail_length() {
        let name = b"@test_read";
        let seq = b"ACGT"; // Too short
        let qual = make_qual(&[30; 4]);
        let config = FilterConfig::short_read();
        let decision = apply_filters(name, seq, &qual, &config);
        assert_eq!(decision, FilterDecision::FailLength);
    }

    #[test]
    fn test_apply_filters_fail_quality() {
        let name = b"@test_read";
        let seq = b"ACGTACGTACGTACGTACGTACGTACGTACGT";
        let qual = make_qual(&[5; 32]); // Low quality
        let config = FilterConfig::short_read();
        let decision = apply_filters(name, seq, &qual, &config);
        assert_eq!(decision, FilterDecision::FailQuality);
    }

    #[test]
    fn test_apply_filters_fail_n_rate() {
        let name = b"@test_read";
        let seq = b"NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN"; // All N
        let qual = make_qual(&[30; 32]);
        let config = FilterConfig::short_read();
        let decision = apply_filters(name, seq, &qual, &config);
        assert_eq!(decision, FilterDecision::FailNRate);
    }

    #[test]
    fn test_apply_filters_fail_complexity() {
        let name = b"@test_read";
        let seq = b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"; // 40 As
        let qual = make_qual(&[30; 40]);
        // Need to enable complexity filter since fastp default is OFF
        let config = FilterConfig::short_read().with_complexity_threshold(0.3);
        let decision = apply_filters(name, seq, &qual, &config);
        assert_eq!(decision, FilterDecision::FailComplexity);
    }

    #[test]
    fn test_config_builder() {
        let config = FilterConfig::new()
            .with_qualified_quality(20)
            .with_unqualified_percent_limit(30.0)
            .with_n_base_limit(10)
            .with_min_avg_quality(20)
            .with_max_n_rate(0.1)
            .with_complexity_threshold(0.5)
            .with_min_length(50)
            .with_max_length(500);
        assert_eq!(config.qualified_quality, 20);
        assert!((config.unqualified_percent_limit - 30.0).abs() < 0.001);
        assert_eq!(config.n_base_limit, 10);
        assert_eq!(config.min_avg_quality, Some(20));
        assert_eq!(config.max_n_rate, Some(0.1));
        assert_eq!(config.low_complexity_threshold, Some(0.5));
        assert_eq!(config.length_config.min_length, 50);
        assert_eq!(config.length_config.max_length, Some(500));
    }

    #[test]
    fn test_config_disable_filters() {
        let config = FilterConfig::short_read()
            .without_quality_filter()
            .without_complexity_filter();
        assert!(!config.quality_filtering_enabled);
        assert!(config.min_avg_quality.is_none());
        assert!(config.low_complexity_threshold.is_none());
    }

    #[test]
    fn test_n_rate_clamping() {
        let config = FilterConfig::new().with_max_n_rate(1.5);
        assert_eq!(config.max_n_rate, Some(1.0));

        let config = FilterConfig::new().with_max_n_rate(-0.5);
        assert_eq!(config.max_n_rate, Some(0.0));
    }

    // ========================================================================
    // Tests for new fastp-compatible functions
    // ========================================================================

    #[test]
    fn test_check_unqualified_percent_pass() {
        // All bases Q30, threshold Q15 -> 0% unqualified
        let qual = make_qual(&[30, 30, 30, 30, 30]);
        assert!(check_unqualified_percent(&qual, 15, 40.0));
    }

    #[test]
    fn test_check_unqualified_percent_fail() {
        // All bases Q10, threshold Q15 -> 100% unqualified, limit 40%
        let qual = make_qual(&[10, 10, 10, 10, 10]);
        assert!(!check_unqualified_percent(&qual, 15, 40.0));
    }

    #[test]
    fn test_check_unqualified_percent_boundary() {
        // 2/5 = 40% unqualified (exactly at limit)
        let qual = make_qual(&[10, 10, 30, 30, 30]);
        assert!(check_unqualified_percent(&qual, 15, 40.0));
        assert!(!check_unqualified_percent(&qual, 15, 39.0)); // Just below limit
    }

    #[test]
    fn test_check_unqualified_percent_empty() {
        assert!(!check_unqualified_percent(&[], 15, 40.0));
    }

    #[test]
    fn test_check_n_base_count_pass() {
        let seq = b"ACGTACGTNN"; // 2 Ns
        assert!(check_n_base_count(seq, 5));
        assert!(check_n_base_count(seq, 2)); // Exactly at limit
    }

    #[test]
    fn test_check_n_base_count_fail() {
        let seq = b"ACGTNNNNNNN"; // 7 Ns
        assert!(!check_n_base_count(seq, 5));
    }

    #[test]
    fn test_check_n_base_count_empty() {
        assert!(check_n_base_count(b"", 5));
    }

    #[test]
    fn test_fastp_style_filtering() {
        // Test that fastp-style filtering works correctly
        let name = b"@test_read";
        let seq = b"ACGTACGTACGTACGTACGTACGTACGTACGT"; // 32bp, good sequence

        // Good quality - should pass
        let good_qual = make_qual(&[30; 32]);
        let config = FilterConfig::short_read();
        assert!(apply_filters(name, seq, &good_qual, &config).is_pass());

        // 50% low quality (> 40% limit) - should fail
        let mut mixed_qual = vec![30u8 + 33; 16]; // 16 good
        mixed_qual.extend(vec![10u8 + 33; 16]);   // 16 bad (Q10 < Q15)
        assert_eq!(
            apply_filters(name, seq, &mixed_qual, &config),
            FilterDecision::FailQuality
        );
    }

    #[test]
    fn test_quality_filtering_disabled() {
        let name = b"@test_read";
        let seq = b"ACGTACGTACGTACGTACGTACGTACGTACGT";
        let bad_qual = make_qual(&[5; 32]); // Very low quality

        let config = FilterConfig::short_read().without_quality_filter();
        // Should pass since quality filtering is disabled
        // (will only fail n_base_limit check if needed)
        assert!(apply_filters(name, seq, &bad_qual, &config).is_pass());
    }
}