use anyhow::{Context, Result};
use clap::{Parser, Subcommand};
use std::collections::HashSet;
use std::path::{Path, PathBuf};

use skope::{derive_sample_name, find_fastx_files, find_fastx_files_recursive, validate_k_s};

const DEFAULT_KMER_LENGTH: u8 = 31;
const DEFAULT_SMER_LENGTH: u8 = 9;

/// Validate the FracMinHash fraction is in (0, 1]
fn validate_fraction(fraction: f64) -> Result<()> {
    if !(fraction > 0.0 && fraction <= 1.0) {
        return Err(anyhow::anyhow!(
            "Invalid --fraction {fraction}: must be in (0, 1] (1 = retain all)"
        ));
    }
    Ok(())
}

fn validate_sample_names(names: &[String]) -> Result<()> {
    let mut seen = HashSet::new();
    let mut duplicates = Vec::new();

    for name in names {
        if !seen.insert(name) && !duplicates.contains(name) {
            duplicates.push(name.clone());
        }
    }

    if !duplicates.is_empty() {
        return Err(anyhow::anyhow!(
            "Duplicate sample names detected: {}. Please rename files or use --names to specify unique names.",
            duplicates.join(", ")
        ));
    }

    Ok(())
}

/// Check if a path is a special input (FIFO / process substitution)
#[cfg(unix)]
fn is_special_input_path(path: &Path) -> bool {
    use std::os::unix::fs::FileTypeExt;
    let path_str = path.to_string_lossy();
    if path_str.starts_with("/dev/fd/") || path_str.starts_with("/proc/self/fd/") {
        return true;
    }
    if let Ok(metadata) = std::fs::metadata(path) {
        return metadata.file_type().is_fifo();
    }
    false
}

#[cfg(not(unix))]
fn is_special_input_path(_path: &Path) -> bool {
    false
}

/// Expand sample inputs (files and directories) into lists of files per sample
fn expand_sample_inputs(inputs: &[PathBuf]) -> Result<(Vec<Vec<PathBuf>>, Vec<bool>)> {
    let mut expanded_samples = Vec::new();
    let mut is_directory = Vec::new();

    for input in inputs {
        // Check stdin special case
        if input.to_string_lossy() == "-" {
            expanded_samples.push(vec![input.clone()]);
            is_directory.push(false);
            continue;
        }

        // Check for special input paths (FIFOs / process substitution)
        if is_special_input_path(input) {
            expanded_samples.push(vec![input.clone()]);
            is_directory.push(false);
            continue;
        }

        // Check path exists
        if !input.exists() {
            return Err(anyhow::anyhow!("Path does not exist: {}", input.display()));
        }

        let metadata = std::fs::metadata(input)
            .with_context(|| format!("Failed to access path: {}", input.display()))?;

        if metadata.is_file() {
            expanded_samples.push(vec![input.clone()]);
            is_directory.push(false);
        } else if metadata.is_dir() {
            let files = find_fastx_files(input)?;
            expanded_samples.push(files);
            is_directory.push(true);
        } else {
            return Err(anyhow::anyhow!(
                "Path is neither a regular file nor directory: {}",
                input.display()
            ));
        }
    }

    Ok((expanded_samples, is_directory))
}

/// Expand background inputs into a flat list of fastx files (directories searched recursively)
fn expand_background_inputs(inputs: &[PathBuf]) -> Result<Vec<PathBuf>> {
    let mut files = Vec::new();
    for input in inputs {
        if input.to_string_lossy() == "-" || is_special_input_path(input) {
            files.push(input.clone());
            continue;
        }
        if !input.exists() {
            return Err(anyhow::anyhow!(
                "Background path does not exist: {}",
                input.display()
            ));
        }
        if std::fs::metadata(input)?.is_dir() {
            files.extend(find_fastx_files_recursive(input)?);
        } else {
            files.push(input.clone());
        }
    }
    Ok(files)
}

/// Parse a base-count string with K/M/G/T suffix into a bp count
fn parse_bases(s: &str) -> Result<u64> {
    let s = s.trim().to_uppercase();
    let (num_str, multiplier) = if s.ends_with('T') {
        (&s[..s.len() - 1], 1_000_000_000_000u64)
    } else if s.ends_with('G') {
        (&s[..s.len() - 1], 1_000_000_000u64)
    } else if s.ends_with('M') {
        (&s[..s.len() - 1], 1_000_000u64)
    } else if s.ends_with('K') {
        (&s[..s.len() - 1], 1_000u64)
    } else {
        (s.as_str(), 1u64)
    };

    let num: u64 = num_str
        .parse()
        .with_context(|| format!("Invalid base count: {}", s))?;

    Ok(num * multiplier)
}

#[derive(Parser)]
#[command(author, version, about = "Containment and abundance estimation using open syncmers", long_about = None)]
struct Cli {
    #[command(subcommand)]
    command: Commands,
}

#[derive(Subcommand)]
enum IndexCommands {
    /// Build a classification index (.sk) from a directory of fastx files/subdirectories (one group per top-level file or directory)
    #[command(alias = "build")]
    BuildClassify {
        /// Directory of fastx files/subdirectories, one group per top-level file or directory
        groups: PathBuf,

        /// K-mer length (1-61, must be odd)
        #[arg(short = 'k', long = "kmer", value_name = "K", default_value_t = DEFAULT_KMER_LENGTH, value_parser = clap::value_parser!(u8).range(1..=61))]
        kmer_length: u8,

        /// S-mer length used for open syncmer selection (s < k, s must be odd)
        #[arg(short = 's', long = "smer", value_name = "S", default_value_t = DEFAULT_SMER_LENGTH)]
        smer_length: u8,

        /// Number of execution threads (0 = auto)
        #[arg(short = 't', long = "threads", default_value_t = 8)]
        threads: usize,

        /// Path to output index file (.sk) (- for stdout)
        #[arg(short = 'o', long = "output", default_value = "-")]
        output: String,

        /// Suppress progress reporting
        #[arg(short = 'q', long = "quiet", default_value_t = false)]
        quiet: bool,
    },

    /// Build a query index (.sk) from target fastx file(s), optionally masking background syncmers
    BuildQuery {
        /// Path to fastx file (single target unless -i) or directory (one target per child file/subdir)
        targets: PathBuf,

        /// Path to fastx file(s) whose syncmers we wish to drop from our targets
        #[arg(short = 'b', long = "background")]
        background: Vec<PathBuf>,

        /// K-mer length (1-61)
        #[arg(short = 'k', long = "kmer", value_name = "K", default_value_t = DEFAULT_KMER_LENGTH, value_parser = clap::value_parser!(u8).range(1..=61))]
        kmer_length: u8,

        /// S-mer length used for open syncmer selection (s < k, s must be odd)
        #[arg(short = 's', long = "smer", value_name = "S", default_value_t = DEFAULT_SMER_LENGTH)]
        smer_length: u8,

        /// Treat each fastx record as a separate target (default: merge records into one target)
        #[arg(short = 'i', long = "individual", default_value_t = false)]
        individual: bool,

        /// Store syncmer positions (needed for --confidence/--dump-syncmers at query time)
        #[arg(short = 'p', long = "positions", default_value_t = false)]
        positions: bool,

        /// Fraction of target syncmers to keep [0, 1]
        #[arg(short = 'f', long = "fraction", value_name = "FLOAT", default_value_t = 1.0)]
        fraction: f64,

        /// Number of execution threads (0 = auto)
        #[arg(short = 't', long = "threads", default_value_t = 8)]
        threads: usize,

        /// Path to output index file (.sk) (- for stdout)
        #[arg(short = 'o', long = "output", default_value = "-")]
        output: String,

        /// Suppress progress reporting
        #[arg(short = 'q', long = "quiet", default_value_t = false)]
        quiet: bool,
    },

    /// Show metadata for an index (.sk)
    Info {
        /// Path to index file (.sk)
        index: PathBuf,
    },
}

#[derive(Subcommand)]
enum Commands {
    /// Estimate target containment & abundance in fastx file(s) or directories thereof using open syncmers
    Query {
        /// Path to fastx file (treated as single target unless -i set), directory of fastx files/subdirs (one target per child file/subdir), or query index (.sk)
        targets: PathBuf,

        /// Path(s) to fastx files/dirs (- for stdin). Each file/dir is treated as a separate sample
        #[arg(required = true)]
        samples: Vec<PathBuf>,

        /// K-mer length (1-61)
        #[arg(short = 'k', long = "kmer", value_name = "K", default_value_t = DEFAULT_KMER_LENGTH, value_parser = clap::value_parser!(u8).range(1..=61))]
        kmer_length: u8,

        /// S-mer length used for syncmer selection (s < k, s must be odd)
        #[arg(short = 's', long = "smer", value_name = "S", default_value_t = DEFAULT_SMER_LENGTH)]
        smer_length: u8,

        /// Treat each fastx record as separate target (default: merge records into one target)
        #[arg(short = 'i', long = "individual", default_value_t = false)]
        individual: bool,

        /// Report confidence intervals, ANI estimates, and patchiness columns
        #[arg(short = 'c', long = "confidence", default_value_t = false)]
        confidence: bool,

        /// Consider only syncmers unique to each target
        #[arg(short = 'd', long = "discriminatory", default_value_t = false)]
        discriminatory: bool,

        /// Collect syncmer positions (implied by --confidence and --dump-syncmers)
        #[arg(short = 'p', long = "positions", default_value_t = false)]
        positions: bool,

        /// Fraction of target syncmers to keep [0, 1]
        #[arg(short = 'f', long = "fraction", value_name = "FLOAT", default_value_t = 1.0)]
        fraction: f64,

        /// Comma-separated additional abundance thresholds for containment estimation
        #[arg(
            short = 'a',
            long = "abundance-thresholds",
            value_name = "INT,...",
            value_delimiter = ',',
            default_value = "10"
        )]
        abundance_thresholds: Vec<usize>,

        /// Path to fastx file(s) whose syncmers we wish to drop from our targets
        #[arg(short = 'b', long = "background")]
        background: Vec<PathBuf>,

        /// Terminate processing after approximately this many bases (e.g. 50M, 10G)
        #[arg(short = 'l', long = "limit", value_name = "BASES")]
        limit: Option<String>,

        /// Number of execution threads (0 = auto)
        #[arg(short = 't', long = "threads", default_value_t = 8)]
        threads: usize,

        /// Path to output file (- for stdout)
        #[arg(short = 'o', long = "output", default_value = "-")]
        output: String,

        /// Comma-separated sample names (default is file/dir name without extension)
        #[arg(short = 'n', long = "names", value_name = "NAME,...", value_delimiter = ',')]
        sample_names: Option<Vec<String>>,

        /// Sort results
        #[arg(long = "sort", default_value = "containment", value_parser = ["containment", "target", "sample", "input"])]
        sort: String,

        /// Dump selected target syncmers to TSV file (target, position, kmer)
        #[arg(long = "dump-syncmers", value_name = "FILE")]
        dump_syncmers: Option<PathBuf>,

        /// Suppress TOTAL summary rows in output
        #[arg(long = "no-total", default_value_t = false)]
        no_total: bool,

        /// Suppress progress reporting
        #[arg(short = 'q', long = "quiet", default_value_t = false)]
        quiet: bool,
    },

    /// Classify sequences into groups by syncmer content
    Classify {
        /// Path to .sk classification index file or directory of fastx files/subdirectories (one group per top-level file or directory)
        index: PathBuf,

        /// Path(s) to fastx files/dirs (- for stdin)
        #[arg(required = true)]
        samples: Vec<PathBuf>,

        /// K-mer length (only used when index is a directory) (1-61, must be odd)
        #[arg(short = 'k', long = "kmer", value_name = "K", default_value_t = DEFAULT_KMER_LENGTH, value_parser = clap::value_parser!(u8).range(1..=61))]
        kmer_length: u8,

        /// S-mer length (only used when index is a directory)
        #[arg(short = 's', long = "smer", value_name = "S", default_value_t = DEFAULT_SMER_LENGTH)]
        smer_length: u8,

        /// Consider only syncmers unique to each group
        #[arg(short = 'd', long = "discriminatory", default_value_t = false)]
        discriminatory: bool,

        /// Minimum absolute number of syncmer hits for a match
        #[arg(short = 'a', long = "abs-threshold", value_name = "ABS_THRESHOLD", default_value_t = 1)]
        abs_threshold: u64,

        /// Minimum relative proportion (0.0-1.0) of syncmer hits for a match
        #[arg(short = 'r', long = "rel-threshold", value_name = "REL_THRESHOLD", default_value_t = 0.0)]
        rel_threshold: f64,

        /// Terminate processing after approximately this many bases (e.g. 50M, 10G)
        #[arg(short = 'l', long = "limit", value_name = "BASES")]
        limit: Option<String>,

        /// Number of execution threads (0 = auto)
        #[arg(short = 't', long = "threads", default_value_t = 8)]
        threads: usize,

        /// Path to output file (- for stdout)
        #[arg(short = 'o', long = "output", default_value = "-")]
        output: String,

        /// Comma-separated sample names (default is file/dir name without extension)
        #[arg(short = 'n', long = "names", value_name = "NAME,...", value_delimiter = ',')]
        sample_names: Option<Vec<String>>,

        /// Output per-sequence classifications instead of summary
        #[arg(long = "per-seq", default_value_t = false)]
        per_seq: bool,

        /// Suppress progress reporting
        #[arg(short = 'q', long = "quiet", default_value_t = false)]
        quiet: bool,
    },

    /// Generate per-group length histograms based on syncmer classification
    Lenhist {
        /// Path to .sk classification index file, directory of fastx files/subdirectories (one group per top-level entry), or - to disable group filtering (single "all" bucket)
        index: PathBuf,

        /// Path(s) to fastx files/dirs (- for stdin). Each file/dir is treated as a separate sample
        #[arg(required = true)]
        samples: Vec<PathBuf>,

        // Algorithm parameters
        /// K-mer length (only used when index is a directory or -) (1-61, must be odd)
        #[arg(short = 'k', long = "kmer", value_name = "K", default_value_t = DEFAULT_KMER_LENGTH, value_parser = clap::value_parser!(u8).range(1..=61))]
        kmer_length: u8,

        /// S-mer length used for open syncmer selection (only used when index is a directory or -)
        #[arg(short = 's', long = "smer", value_name = "S", default_value_t = DEFAULT_SMER_LENGTH)]
        smer_length: u8,

        /// Consider only syncmers unique to each group
        #[arg(short = 'd', long = "discriminatory", default_value_t = false)]
        discriminatory: bool,

        /// Minimum absolute number of syncmer hits for a match
        #[arg(short = 'a', long = "abs-threshold", value_name = "ABS_THRESHOLD", default_value_t = 1)]
        abs_threshold: u64,

        /// Minimum relative proportion (0.0-1.0) of syncmer hits for a match
        #[arg(short = 'r', long = "rel-threshold", value_name = "REL_THRESHOLD", default_value_t = 0.0)]
        rel_threshold: f64,

        // Processing options
        /// Terminate processing after approximately this many bases (e.g. 50M, 10G)
        #[arg(short = 'l', long = "limit", value_name = "BASES")]
        limit: Option<String>,

        /// Number of execution threads (0 = auto)
        #[arg(short = 't', long = "threads", default_value_t = 8)]
        threads: usize,

        // Output options
        /// Path to output file (- for stdout)
        #[arg(short = 'o', long = "output", default_value = "-")]
        output: String,

        /// Comma-separated sample names (default is file/dir name without extension)
        #[arg(short = 'n', long = "names", value_name = "NAME,...", value_delimiter = ',')]
        sample_names: Option<Vec<String>>,

        /// Suppress progress reporting
        #[arg(short = 'q', long = "quiet", default_value_t = false)]
        quiet: bool,
    },

    /// Build and manage query and classification indexes
    #[command(subcommand)]
    Index(IndexCommands),
}

fn main() -> Result<()> {
    // Check we have either AVX2 or NEON for SIMD acceleration
    #[cfg(not(any(target_feature = "avx2", target_feature = "neon")))]
    {
        eprintln!(
            "Warning: SIMD acceleration is unavailable. For best performance, compile with `cargo build --release -C target-cpu=native`"
        );
    }

    let cli = Cli::parse();

    match &cli.command {
        Commands::Index(index_cmd) => match index_cmd {
            IndexCommands::BuildClassify {
                groups,
                kmer_length,
                smer_length,
                threads,
                output,
                quiet,
            } => {
                validate_k_s(*kmer_length, *smer_length)?;

                if !groups.is_dir() {
                    return Err(anyhow::anyhow!(
                        "Groups path must be a directory: {}",
                        groups.display()
                    ));
                }

                if *threads > 0 {
                    rayon::ThreadPoolBuilder::new()
                        .num_threads(*threads)
                        .build_global()
                        .context("Failed to initialise thread pool")?;
                }

                let config = skope::BuildClassifyConfig {
                    groups_dir: groups.clone(),
                    kmer_length: *kmer_length,
                    smer_length: *smer_length,
                    threads: *threads,
                    output_path: if output == "-" {
                        None
                    } else {
                        Some(PathBuf::from(output))
                    },
                    quiet: *quiet,
                };

                skope::build_classification_index(&config)
                    .context("Failed to build classification index")?;
            }

            IndexCommands::BuildQuery {
                targets,
                background,
                kmer_length,
                smer_length,
                individual,
                positions,
                fraction,
                threads,
                output,
                quiet,
            } => {
                validate_k_s(*kmer_length, *smer_length)?;
                validate_fraction(*fraction)?;

                if *threads > 0 {
                    rayon::ThreadPoolBuilder::new()
                        .num_threads(*threads)
                        .build_global()
                        .context("Failed to initialise thread pool")?;
                }

                let config = skope::BuildQueryConfig {
                    targets_path: targets.clone(),
                    background_paths: expand_background_inputs(background)?,
                    kmer_length: *kmer_length,
                    smer_length: *smer_length,
                    individual: *individual,
                    positions: *positions,
                    threads: *threads,
                    output_path: if output == "-" {
                        None
                    } else {
                        Some(PathBuf::from(output))
                    },
                    quiet: *quiet,
                    fraction: *fraction,
                };

                skope::build_query_index(&config).context("Failed to build query index")?;
            }

            IndexCommands::Info { index } => {
                skope::run_index_info(index).context("Failed to run index info command")?;
            }
        },

        Commands::Classify {
            index,
            samples,
            sample_names,
            kmer_length,
            smer_length,
            abs_threshold,
            rel_threshold,
            discriminatory,
            threads,
            limit,
            output,
            per_seq,
            quiet,
        } => {
            let (expanded_samples, is_directory) = expand_sample_inputs(samples)?;

            let derived_sample_names: Vec<String> = if let Some(names) = sample_names {
                if names.len() != expanded_samples.len() {
                    return Err(anyhow::anyhow!(
                        "Number of sample names ({}) must match number of samples ({})",
                        names.len(),
                        expanded_samples.len()
                    ));
                }
                names.clone()
            } else {
                samples
                    .iter()
                    .zip(&is_directory)
                    .map(|(p, &is_dir)| derive_sample_name(p, is_dir))
                    .collect()
            };

            validate_sample_names(&derived_sample_names)?;

            // Only validate k/s when index is a directory
            if index.is_dir() {
                validate_k_s(*kmer_length, *smer_length)?;
            }

            if *threads > 0 {
                rayon::ThreadPoolBuilder::new()
                    .num_threads(*threads)
                    .build_global()
                    .context("Failed to initialise thread pool")?;
            }

            let limit_bp = if let Some(s) = limit {
                Some(parse_bases(s)?)
            } else {
                None
            };

            let config = skope::ClassifyConfig {
                index_path: index.clone(),
                sample_paths: expanded_samples,
                sample_names: derived_sample_names,
                kmer_length: *kmer_length,
                smer_length: *smer_length,
                abs_threshold: *abs_threshold,
                rel_threshold: *rel_threshold,
                threads: *threads,
                limit_bp,
                output_path: if output == "-" {
                    None
                } else {
                    Some(PathBuf::from(output))
                },
                per_seq: *per_seq,
                discriminatory: *discriminatory,
                quiet: *quiet,
            };

            skope::run_classification(&config).context("Failed to run classification")?;
        }

        Commands::Query {
            targets,
            samples,
            sample_names,
            kmer_length,
            smer_length,
            threads,
            output,
            quiet,
            abundance_thresholds,
            discriminatory,
            individual,
            limit,
            sort,
            dump_syncmers,
            no_total,
            confidence,
            background,
            positions,
            fraction,
        } => {
            // Expand directories to lists of files
            let (expanded_samples, is_directory) = expand_sample_inputs(samples)?;
            let background_paths = expand_background_inputs(background)?;

            // Derive or validate sample names
            let derived_sample_names: Vec<String> = if let Some(names) = sample_names {
                // User-provided names
                if names.len() != expanded_samples.len() {
                    return Err(anyhow::anyhow!(
                        "Number of sample names ({}) must match number of samples ({})",
                        names.len(),
                        expanded_samples.len()
                    ));
                }
                names.clone()
            } else {
                // Derive from filenames/directory names
                samples
                    .iter()
                    .zip(&is_directory)
                    .map(|(p, &is_dir)| derive_sample_name(p, is_dir))
                    .collect()
            };

            // Validate uniqueness and FracMinHash fraction before resolving target/index parameters
            validate_sample_names(&derived_sample_names)?;
            validate_fraction(*fraction)?;

            // A prebuilt index carries its own k/s and fraction; else validate CLI k/s
            let is_index = !targets.is_dir() && skope::is_query_index(targets);
            let (eff_k, eff_s) = if is_index {
                let (k, s, _fraction) = skope::read_query_index_meta(targets)?;
                (k, s)
            } else {
                validate_k_s(*kmer_length, *smer_length)?;
                (*kmer_length, *smer_length)
            };

            // Configure thread pool if specified (non-zero)
            if *threads > 0 {
                rayon::ThreadPoolBuilder::new()
                    .num_threads(*threads)
                    .build_global()
                    .context("Failed to initialise thread pool")?;
            }

            // Parse limit if provided
            let limit_bp = if let Some(s) = limit {
                Some(parse_bases(s)?)
            } else {
                None
            };

            // Parse sort order
            let sort_order = match sort.as_str() {
                "input" => skope::SortOrder::Original,
                "target" => skope::SortOrder::Target,
                "sample" => skope::SortOrder::Sample,
                "containment" => skope::SortOrder::Containment,
                _ => unreachable!("clap should have validated the sort order"),
            };

            let config = skope::ContainmentConfig {
                targets_path: targets.clone(),
                background_paths,
                sample_paths: expanded_samples,
                sample_names: derived_sample_names,
                kmer_length: eff_k,
                smer_length: eff_s,
                threads: *threads,
                output_path: if output == "-" {
                    None
                } else {
                    Some(PathBuf::from(output))
                },
                quiet: *quiet,
                abundance_thresholds: abundance_thresholds.clone(),
                discriminatory: *discriminatory,
                positions: *positions,
                individual: *individual,
                limit_bp,
                sort_order,
                dump_syncmers_path: dump_syncmers.clone(),
                no_total: *no_total,
                confidence: *confidence,
                fraction: *fraction,
            };

            config
                .execute()
                .context("Failed to run containment analysis")?;
        }
        Commands::Lenhist {
            index,
            samples,
            sample_names,
            kmer_length,
            smer_length,
            abs_threshold,
            rel_threshold,
            discriminatory,
            threads,
            output,
            quiet,
            limit,
        } => {
            // Expand directories to lists of files
            let (expanded_samples, is_directory) = expand_sample_inputs(samples)?;

            // Derive or validate sample names
            let derived_sample_names: Vec<String> = if let Some(names) = sample_names {
                // User-provided names
                if names.len() != expanded_samples.len() {
                    return Err(anyhow::anyhow!(
                        "Number of sample names ({}) must match number of samples ({})",
                        names.len(),
                        expanded_samples.len()
                    ));
                }
                names.clone()
            } else {
                // Derive from filenames/directory names
                samples
                    .iter()
                    .zip(&is_directory)
                    .map(|(p, &is_dir)| derive_sample_name(p, is_dir))
                    .collect()
            };

            // Validate uniqueness
            validate_sample_names(&derived_sample_names)?;

            validate_k_s(*kmer_length, *smer_length)?;

            // Configure thread pool if specified (non-zero)
            if *threads > 0 {
                rayon::ThreadPoolBuilder::new()
                    .num_threads(*threads)
                    .build_global()
                    .context("Failed to initialise thread pool")?;
            }

            // Parse limit if provided
            let limit_bp = if let Some(s) = limit {
                Some(parse_bases(s)?)
            } else {
                None
            };

            // Detect if user wants all seqs (no group filtering)
            let include_all_seqs = index.to_string_lossy() == "-";

            let config = skope::LengthHistogramConfig {
                index_path: index.clone(),
                sample_paths: expanded_samples,
                sample_names: derived_sample_names,
                kmer_length: *kmer_length,
                smer_length: *smer_length,
                abs_threshold: *abs_threshold,
                rel_threshold: *rel_threshold,
                discriminatory: *discriminatory,
                threads: *threads,
                output_path: if output == "-" {
                    None
                } else {
                    Some(PathBuf::from(output))
                },
                quiet: *quiet,
                limit_bp,
                include_all_seqs,
            };

            config
                .execute()
                .context("Failed to run length histogram analysis")?;
        }
    }

    Ok(())
}
