use std::collections::{HashMap, HashSet};
use std::fs::File;
use std::io::{BufRead, BufReader, BufWriter, Write};
use std::path::{Path, PathBuf};
use std::sync::Arc;

use anyhow::{anyhow, bail, Context, Result};
use chrono::NaiveDate;
use rand::rngs::StdRng;
use rand::{Rng, SeedableRng};
use rayon::prelude::*;
use rayon::ThreadPoolBuilder;
use serde::Deserialize;

use crate::download::ensure_reference_db;
use crate::stats::{IlluminaProbeRecord, ReferenceVariant, StatsStore};
use crate::{SyntheticArgs, SyntheticFormat};

const HEADER_TEXT: &str = r#"# This data file generated by Dynamic DNA (DDNA) Laboratories at: Thu Nov 7 16:03:14 2024
#						
# This file contains raw genetic data, including data that is not used in DDNA reports. 						
# This data has undergone a general quality review however only a subset of markers have been 						
# individually reviewed for accuracy. As such, this data is suitable only for research,  						
# educational, and informational use and not for any medical or diagnostic use. 						
# 						
# Below is a text version of your data.  Fields are TAB-separated and						
# each line corresponds to a single SNP.  For each SNP, we provide its identifier,  						
# it chromosomal location realtive to Build 38 of the human reference genome, and the 						
# genotype call oriented with respect to the plus strand on the human reference sequence.						
# More information about Dynamic DNA Laboratories can be found at:						
# https://dynamicdnalabs.com						
# 						
# More information on reference human assembly builds:						
# https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.40/						
#						
# rsid	chromosome	position	genotype	gs	baf	lrr
"#;

enum VariantKind {
    Snp,
    Mnv,
    Insertion,
    Deletion,
}

pub fn run_synthetic(args: SyntheticArgs) -> Result<()> {
    if !(0.0..=1.0).contains(&args.alt_frequency) {
        bail!("--alt-frequency must be between 0 and 1");
    }
    if !(0.0..=1.0).contains(&args.no_call_frequency) {
        bail!("--no-call-frequency must be between 0 and 1");
    }
    if args.no_call_token.trim().is_empty() {
        bail!("--no-call-token must be non-empty");
    }
    if args.count == 0 {
        bail!("--count must be at least 1");
    }
    let output_template = args.output.to_string_lossy().to_string();
    if args.count > 1 && !(output_template.contains("{index}") || output_template.contains("{id}"))
    {
        bail!("When --count > 1, --output must include an '{{index}}' or '{{id}}' placeholder");
    }
    if args.id_min > args.id_max {
        bail!("--id-min must be <= --id-max");
    }
    if args.month_min < 1 || args.month_max > 12 || args.month_min > args.month_max {
        bail!("Month range must be between 1 and 12");
    }
    if args.day_min < 1 || args.day_max > 31 || args.day_min > args.day_max {
        bail!("Day range must be between 1 and 31");
    }

    let sqlite_path = resolve_sqlite_path(&args)?;
    let store = if read_only_db_requested() {
        StatsStore::connect_read_only(&sqlite_path)?
    } else {
        StatsStore::connect(&sqlite_path)?
    };
    let references = store.all_references(args.limit)?;
    if references.is_empty() {
        bail!(
            "No reference rows found in {}",
            sqlite_path.to_string_lossy()
        );
    }
    let references = Arc::new(references);

    let overlays = load_overlay_specs(&args)?
        .map(Arc::new)
        .unwrap_or_else(|| Arc::new(Vec::new()));
    let illumina_manifest = load_illumina_probe_manifest(&args, &store)?.map(Arc::new);
    let biallelic_choices = if args.biallelic {
        Some(Arc::new(build_biallelic_choices(
            references.as_ref(),
            overlays.as_ref(),
        )?))
    } else {
        None
    };

    let plans = build_file_plans(&output_template, &args)?;

    let thread_count = args
        .threads
        .unwrap_or_else(|| rayon::current_num_threads().max(1));
    let pool = ThreadPoolBuilder::new()
        .num_threads(thread_count)
        .build()
        .context("build synthetic thread pool")?;

    let results: Vec<usize> = pool.install(|| {
        plans
            .par_iter()
            .map(|plan| {
                write_single_file(
                    &plan.path,
                    references.as_ref(),
                    overlays.as_ref(),
                    args.alt_frequency,
                    args.no_call_frequency,
                    args.no_call_token.as_str(),
                    plan.seed,
                    args.format,
                    biallelic_choices.as_ref().map(|choices| choices.as_ref()),
                    illumina_manifest.as_ref().map(|manifest| manifest.as_ref()),
                    args.qc_edge_cases,
                )
            })
            .collect::<Result<Vec<_>>>()
    })?;
    if args.biallelic {
        validate_biallelic_outputs(&plans, args.format)?;
    }

    let total_rows: usize = results.iter().sum();
    println!(
        "🧪 Generated {} file(s), {} total rows (alt freq {:.2}%)",
        results.len(),
        total_rows,
        args.alt_frequency * 100.0
    );
    Ok(())
}

fn read_only_db_requested() -> bool {
    match std::env::var("BVS_READ_ONLY_DB") {
        Ok(value) => matches!(
            value.to_ascii_lowercase().as_str(),
            "1" | "true" | "yes" | "on"
        ),
        Err(_) => false,
    }
}

fn resolve_sqlite_path(args: &SyntheticArgs) -> Result<PathBuf> {
    if read_only_db_requested() {
        let baked_in = PathBuf::from("/app/data/genostats.sqlite");
        if baked_in.exists() {
            return Ok(baked_in);
        }
    }
    ensure_reference_db(Some(&args.sqlite), args.force_download)
}

#[allow(clippy::too_many_arguments)]
fn write_single_file(
    path: &PathBuf,
    references: &[ReferenceVariant],
    overlays: &[OverlaySpec],
    alt_frequency: f64,
    no_call_frequency: f64,
    no_call_token: &str,
    seed: Option<u64>,
    format: SyntheticFormat,
    biallelic_choices: Option<&HashMap<i64, BiallelicChoice>>,
    illumina_manifest: Option<&IlluminaProbeManifest>,
    qc_edge_cases: bool,
) -> Result<usize> {
    if let Some(parent) = path.parent() {
        if !parent.as_os_str().is_empty() {
            std::fs::create_dir_all(parent)
                .with_context(|| format!("Create directory {:?}", parent))?;
        }
    }
    let mut rng = match seed {
        Some(seed) => StdRng::seed_from_u64(seed),
        None => StdRng::from_entropy(),
    };
    let mut forced_no_call_emitted = false;
    let default_frequencies = GenotypeFrequencies::from_alt_frequency(alt_frequency)?;

    let mut overlay_assignments =
        prepare_overlay_assignments(overlays, default_frequencies, biallelic_choices, &mut rng)?;

    let sample_id = path
        .file_stem()
        .map(|s| s.to_string_lossy().to_string())
        .unwrap_or_else(|| "SAMPLE".to_string());

    let file = File::create(path).with_context(|| format!("Create {:?}", path))?;
    let mut writer = BufWriter::new(file);
    match format {
        SyntheticFormat::Ddna => {
            writer
                .write_all(HEADER_TEXT.as_bytes())
                .context("write header")?;
        }
        SyntheticFormat::Illumina => {
            let num_snps = if let Some(manifest) = illumina_manifest {
                count_illumina_manifest_outputs(manifest, references, overlays)
            } else {
                let ref_rsids: HashSet<i64> = references.iter().map(|r| r.rsid).collect();
                let extra = overlay_assignments
                    .keys()
                    .filter(|rsid| !ref_rsids.contains(rsid))
                    .count();
                references.len() + extra
            };
            let num_snps = if qc_edge_cases {
                num_snps + ILLUMINA_QC_EDGE_CASES.len()
            } else {
                num_snps
            };
            writer
                .write_all(illumina_header(num_snps).as_bytes())
                .context("write header")?;
        }
    }

    if format == SyntheticFormat::Illumina {
        if let Some(manifest) = illumina_manifest {
            return write_illumina_manifest_file(
                &mut writer,
                manifest,
                references,
                overlays,
                &mut overlay_assignments,
                &default_frequencies,
                no_call_frequency,
                no_call_token,
                biallelic_choices,
                &sample_id,
                &mut rng,
                &mut forced_no_call_emitted,
                qc_edge_cases,
            );
        }
    }

    let mut written = 0usize;
    for reference in references {
        if let Some(assignment) = overlay_assignments.remove(&reference.rsid) {
            let genotype = apply_no_call(
                &assignment.genotype,
                no_call_frequency,
                no_call_token,
                &mut rng,
            );
            let genotype =
                force_first_no_call(genotype, no_call_token, &mut forced_no_call_emitted);
            write_row(
                &mut writer,
                format,
                &sample_id,
                &format!("rs{}", assignment.spec.rsid),
                &format!("rs{}", assignment.spec.rsid),
                &assignment.spec.chromosome,
                assignment.spec.position,
                &genotype,
                &design_for_assignment(&assignment, biallelic_choices, &genotype),
                &mut rng,
            )?;
            written += 1;
            continue;
        }

        let genotype = synthesize_reference_genotype(
            reference,
            &default_frequencies,
            biallelic_choices,
            &mut rng,
        );
        let genotype = apply_no_call(&genotype, no_call_frequency, no_call_token, &mut rng);
        let genotype = force_first_no_call(genotype, no_call_token, &mut forced_no_call_emitted);
        let (design_reference, design_alts) =
            design_parts_for_reference(reference, biallelic_choices);
        write_row(
            &mut writer,
            format,
            &sample_id,
            &format!("rs{}", reference.rsid),
            &format!("rs{}", reference.rsid),
            &reference.chromosome,
            reference.position,
            &genotype,
            &design_for(Some(&design_reference), &design_alts, &genotype),
            &mut rng,
        )
        .with_context(|| format!("write row for rs{}", reference.rsid))?;
        written += 1;
    }

    for assignment in overlay_assignments.into_values() {
        let genotype = apply_no_call(
            &assignment.genotype,
            no_call_frequency,
            no_call_token,
            &mut rng,
        );
        let genotype = force_first_no_call(genotype, no_call_token, &mut forced_no_call_emitted);
        write_row(
            &mut writer,
            format,
            &sample_id,
            &format!("rs{}", assignment.spec.rsid),
            &format!("rs{}", assignment.spec.rsid),
            &assignment.spec.chromosome,
            assignment.spec.position,
            &genotype,
            &design_for_assignment(&assignment, biallelic_choices, &genotype),
            &mut rng,
        )?;
        written += 1;
    }

    if qc_edge_cases {
        written += write_qc_edge_cases(&mut writer, format, &sample_id)?;
    }

    writer.flush()?;
    Ok(written)
}

fn illumina_header(num_snps: usize) -> String {
    format!(
        "[Header]\n\
         GSGT Version\t2.0.5\n\
         Processing Date\tBioVault Synthetic\n\
         Content\t\tBioVaultSynthetic.bpm\n\
         Num SNPs\t{n}\n\
         Total SNPs\t{n}\n\
         Num Samples\t1\n\
         Total Samples\t1\n\
         File \t1 of 1\n\
         [Data]\n\
         Sample ID\tSample Name\tSNP Name\tSNP\tChr\tPosition\t\
         Allele1 - Top\tAllele2 - Top\tAllele1 - Forward\tAllele2 - Forward\t\
         Allele1 - Plus\tAllele2 - Plus\tAllele1 - Design\tAllele2 - Design\t\
         Allele1 - AB\tAllele2 - AB\tPlus/Minus Strand\n",
        n = num_snps
    )
}

#[allow(clippy::too_many_arguments)]
fn write_illumina_manifest_file(
    writer: &mut BufWriter<File>,
    manifest: &IlluminaProbeManifest,
    references: &[ReferenceVariant],
    overlays: &[OverlaySpec],
    overlay_assignments: &mut HashMap<i64, OverlayAssignment>,
    default_frequencies: &GenotypeFrequencies,
    no_call_frequency: f64,
    no_call_token: &str,
    biallelic_choices: Option<&HashMap<i64, BiallelicChoice>>,
    sample_id: &str,
    rng: &mut StdRng,
    forced_no_call_emitted: &mut bool,
    qc_edge_cases: bool,
) -> Result<usize> {
    let references_by_rsid = references
        .iter()
        .map(|reference| (reference.rsid, reference))
        .collect::<HashMap<_, _>>();
    let references_by_locus = references
        .iter()
        .map(|reference| {
            (
                (reference.chromosome.clone(), reference.position),
                reference,
            )
        })
        .collect::<HashMap<_, _>>();
    let overlays_by_locus = overlays
        .iter()
        .map(|overlay| ((overlay.chromosome.clone(), overlay.position), overlay.rsid))
        .collect::<HashMap<_, _>>();
    let mut genotype_by_rsid: HashMap<i64, String> = HashMap::new();
    let mut represented_rsids: HashSet<i64> = HashSet::new();
    let mut explicit_probe_rsids: HashSet<i64> = HashSet::new();
    let mut written = 0usize;

    for probe in &manifest.probes {
        let rsid = resolve_probe_rsid(probe, &references_by_locus, &overlays_by_locus);
        if let Some(explicit_rsid) = probe.rsid {
            explicit_probe_rsids.insert(explicit_rsid);
        }

        let genotype = if let Some(rsid) = rsid {
            represented_rsids.insert(rsid);
            if let Some(existing) = genotype_by_rsid.get(&rsid) {
                existing.clone()
            } else {
                let raw = if let Some(assignment) = overlay_assignments.get(&rsid) {
                    assignment.genotype.clone()
                } else if let Some(reference) = references_by_rsid.get(&rsid) {
                    synthesize_reference_genotype(
                        reference,
                        default_frequencies,
                        biallelic_choices,
                        rng,
                    )
                } else {
                    synthesize_genotype_from_design(&probe.design, default_frequencies, rng)
                };
                genotype_by_rsid.insert(rsid, raw.clone());
                raw
            }
        } else {
            synthesize_genotype_from_design(&probe.design, default_frequencies, rng)
        };
        let genotype = apply_no_call(&genotype, no_call_frequency, no_call_token, rng);
        let genotype = force_first_no_call(genotype, no_call_token, forced_no_call_emitted);

        write_row(
            writer,
            SyntheticFormat::Illumina,
            sample_id,
            &probe.probe_id,
            &probe.probe_id,
            &probe.chromosome,
            probe.position,
            &genotype,
            &probe.design,
            rng,
        )?;
        if let Some(rsid) = rsid {
            overlay_assignments.remove(&rsid);
        }
        written += 1;
    }

    for reference in references {
        if explicit_probe_rsids.contains(&reference.rsid) {
            continue;
        }
        if let Some(clean_rsids) = &manifest.clean_rsids {
            if !clean_rsids.contains(&reference.rsid) {
                continue;
            }
        }
        let genotype = if let Some(assignment) = overlay_assignments.remove(&reference.rsid) {
            assignment.genotype
        } else {
            synthesize_reference_genotype(reference, default_frequencies, biallelic_choices, rng)
        };
        let genotype = apply_no_call(&genotype, no_call_frequency, no_call_token, rng);
        let genotype = force_first_no_call(genotype, no_call_token, forced_no_call_emitted);
        let (design_reference, design_alts) =
            design_parts_for_reference(reference, biallelic_choices);
        let design = design_for(Some(&design_reference), &design_alts, &genotype);
        let rsid_label = format!("rs{}", reference.rsid);
        write_row(
            writer,
            SyntheticFormat::Illumina,
            sample_id,
            &rsid_label,
            &rsid_label,
            &reference.chromosome,
            reference.position,
            &genotype,
            &design,
            rng,
        )?;
        written += 1;
    }

    for assignment in overlay_assignments.values() {
        if represented_rsids.contains(&assignment.spec.rsid) {
            continue;
        }
        let genotype = apply_no_call(&assignment.genotype, no_call_frequency, no_call_token, rng);
        let genotype = force_first_no_call(genotype, no_call_token, forced_no_call_emitted);
        let rsid_label = format!("rs{}", assignment.spec.rsid);
        write_row(
            writer,
            SyntheticFormat::Illumina,
            sample_id,
            &rsid_label,
            &rsid_label,
            &assignment.spec.chromosome,
            assignment.spec.position,
            &genotype,
            &design_for_assignment(assignment, biallelic_choices, &genotype),
            rng,
        )?;
        written += 1;
    }

    if qc_edge_cases {
        written += write_qc_edge_cases(writer, SyntheticFormat::Illumina, sample_id)?;
    }

    writer.flush()?;
    Ok(written)
}

fn count_illumina_manifest_outputs(
    manifest: &IlluminaProbeManifest,
    references: &[ReferenceVariant],
    overlays: &[OverlaySpec],
) -> usize {
    let represented_rsids = manifest
        .probes
        .iter()
        .filter_map(|probe| probe.rsid)
        .collect::<HashSet<_>>();
    let reference_rsids = references
        .iter()
        .map(|reference| reference.rsid)
        .collect::<HashSet<_>>();
    let clean_references = references
        .iter()
        .filter(|reference| !represented_rsids.contains(&reference.rsid))
        .filter(|reference| {
            manifest
                .clean_rsids
                .as_ref()
                .map(|clean_rsids| clean_rsids.contains(&reference.rsid))
                .unwrap_or(true)
        })
        .count();
    let extra_overlays = overlays
        .iter()
        .filter(|overlay| {
            !reference_rsids.contains(&overlay.rsid) && !represented_rsids.contains(&overlay.rsid)
        })
        .count();
    manifest.probes.len() + clean_references + extra_overlays
}

fn resolve_probe_rsid(
    probe: &IlluminaProbe,
    references_by_locus: &HashMap<(String, i64), &ReferenceVariant>,
    overlays_by_locus: &HashMap<(String, i64), i64>,
) -> Option<i64> {
    probe.rsid.or_else(|| {
        references_by_locus
            .get(&(probe.chromosome.clone(), probe.position))
            .map(|reference| reference.rsid)
            .or_else(|| {
                overlays_by_locus
                    .get(&(probe.chromosome.clone(), probe.position))
                    .copied()
            })
    })
}

fn force_first_no_call(genotype: String, no_call_token: &str, emitted: &mut bool) -> String {
    if *emitted {
        return genotype;
    }
    *emitted = true;
    no_call_token.to_string()
}

/// Split a synthetic genotype string into two single-token alleles.
/// No-call / unknown -> ("-", "-"); 2-char -> (first, last); 1-char -> dup.
fn split_alleles(genotype: &str) -> (String, String) {
    let g = genotype.trim();
    if g.is_empty() || g.chars().all(|c| matches!(c, '-' | '.' | 'N' | 'n' | '0')) {
        return ("-".to_string(), "-".to_string());
    }
    let chars: Vec<char> = g.chars().collect();
    if chars.len() == 1 {
        let s = chars[0].to_string();
        return (s.clone(), s);
    }
    (chars[0].to_string(), chars[chars.len() - 1].to_string())
}

fn design_for(reference: Option<&str>, alternates: &[String], genotype: &str) -> String {
    let alt_refs: Vec<&str> = alternates.iter().map(|s| s.as_str()).collect();
    let kind = determine_variant_kind(reference.unwrap_or(""), &alt_refs);
    let (a, b) = match kind {
        VariantKind::Insertion => ("D".to_string(), "I".to_string()),
        VariantKind::Deletion => ("I".to_string(), "D".to_string()),
        VariantKind::Snp | VariantKind::Mnv => {
            let r = reference
                .and_then(|s| s.chars().next())
                .map(|c| c.to_string());
            let alt = alternates
                .first()
                .and_then(|s| s.chars().next())
                .map(|c| c.to_string());
            match (r, alt) {
                (Some(r), Some(alt)) => (r, alt),
                _ => {
                    let (x, y) = split_alleles(genotype);
                    if x == y {
                        (x, "N".to_string())
                    } else {
                        (x, y)
                    }
                }
            }
        }
    };
    format!("[{}/{}]", a, b)
}

fn ab_code(allele: &str, design_ref: char) -> &'static str {
    match allele.chars().next() {
        Some('-') => "-",
        Some(c) if c == design_ref => "A",
        Some(_) => "B",
        None => "-",
    }
}

struct DdnaQcEdgeCase {
    rsid: &'static str,
    chromosome: &'static str,
    position: i64,
    genotype: &'static str,
    gs: &'static str,
    baf: &'static str,
    lrr: &'static str,
}

struct IlluminaQcEdgeCase {
    snp_name: &'static str,
    design: &'static str,
    chromosome: &'static str,
    position: i64,
    allele1: &'static str,
    allele2: &'static str,
}

const DDNA_QC_EDGE_CASES: &[DdnaQcEdgeCase] = &[
    DdnaQcEdgeCase {
        rsid: "rs900000001",
        chromosome: "1",
        position: 900000001,
        genotype: "II",
        gs: "0.3333",
        baf: "0.0220",
        lrr: "-0.3294",
    },
    DdnaQcEdgeCase {
        rsid: "rs900000002",
        chromosome: "1",
        position: 900000002,
        genotype: "--",
        gs: "0.0000",
        baf: "NaN",
        lrr: "NaN",
    },
    DdnaQcEdgeCase {
        rsid: "rs900000003",
        chromosome: "6",
        position: 900000003,
        genotype: "--",
        gs: "NaN",
        baf: "NaN",
        lrr: "NaN",
    },
    DdnaQcEdgeCase {
        rsid: ".",
        chromosome: "1",
        position: 900000004,
        genotype: "DD",
        gs: "0.4630",
        baf: "0.0000",
        lrr: "-0.0143",
    },
    DdnaQcEdgeCase {
        rsid: ".",
        chromosome: "1",
        position: 900000005,
        genotype: "TT",
        gs: "0.4222",
        baf: "0.9976",
        lrr: "-0.0408",
    },
];

const ILLUMINA_QC_EDGE_CASES: &[IlluminaQcEdgeCase] = &[
    IlluminaQcEdgeCase {
        snp_name: "BV_QC_UNMAPPED_0_0",
        design: "[A/C]",
        chromosome: "0",
        position: 0,
        allele1: "C",
        allele2: "C",
    },
    IlluminaQcEdgeCase {
        snp_name: "BV_QC_INDEL_I_D",
        design: "[I/D]",
        chromosome: "1",
        position: 900000006,
        allele1: "I",
        allele2: "I",
    },
];

fn write_qc_edge_cases(
    writer: &mut BufWriter<File>,
    format: SyntheticFormat,
    sample_id: &str,
) -> Result<usize> {
    match format {
        SyntheticFormat::Ddna => {
            for row in DDNA_QC_EDGE_CASES {
                writeln!(
                    writer,
                    "{}\t{}\t{}\t{}\t{}\t{}\t{}",
                    row.rsid, row.chromosome, row.position, row.genotype, row.gs, row.baf, row.lrr
                )
                .with_context(|| format!("write QC edge-case row for {}", row.rsid))?;
            }
            Ok(DDNA_QC_EDGE_CASES.len())
        }
        SyntheticFormat::Illumina => {
            for row in ILLUMINA_QC_EDGE_CASES {
                let design_ref = row
                    .design
                    .trim_start_matches('[')
                    .chars()
                    .next()
                    .unwrap_or('N');
                let ab1 = ab_code(row.allele1, design_ref);
                let ab2 = ab_code(row.allele2, design_ref);
                writeln!(
                    writer,
                    "{sid}\t\t{snp}\t{design}\t{chr}\t{pos}\t\
                     {a1}\t{a2}\t{a1}\t{a2}\t{a1}\t{a2}\t{a1}\t{a2}\t{ab1}\t{ab2}\t+",
                    sid = sample_id,
                    snp = row.snp_name,
                    design = row.design,
                    chr = row.chromosome,
                    pos = row.position,
                    a1 = row.allele1,
                    a2 = row.allele2,
                    ab1 = ab1,
                    ab2 = ab2,
                )
                .with_context(|| format!("write QC edge-case row for {}", row.snp_name))?;
            }
            Ok(ILLUMINA_QC_EDGE_CASES.len())
        }
    }
}

#[derive(Debug, Clone)]
struct BiallelicChoice {
    reference: String,
    alternates: Vec<String>,
}

fn build_biallelic_choices(
    references: &[ReferenceVariant],
    overlays: &[OverlaySpec],
) -> Result<HashMap<i64, BiallelicChoice>> {
    let mut choices = HashMap::new();
    for reference in references {
        choices.insert(
            reference.rsid,
            BiallelicChoice {
                reference: reference.reference.clone(),
                alternates: first_alternate_only(&reference.alternates),
            },
        );
    }
    for overlay in overlays {
        if overlay.genotype_options.is_empty() {
            let Some(reference) = &overlay.reference else {
                continue;
            };
            choices
                .entry(overlay.rsid)
                .or_insert_with(|| BiallelicChoice {
                    reference: reference.clone(),
                    alternates: overlay.alternates.first().cloned().into_iter().collect(),
                });
        }
    }
    Ok(choices)
}

fn first_alternate_only(alternates: &str) -> Vec<String> {
    parse_alternates(alternates).into_iter().take(1).collect()
}

fn synthesize_reference_genotype(
    reference: &ReferenceVariant,
    default_frequencies: &GenotypeFrequencies,
    biallelic_choices: Option<&HashMap<i64, BiallelicChoice>>,
    rng: &mut StdRng,
) -> String {
    if let Some(choice) = biallelic_choices.and_then(|choices| choices.get(&reference.rsid)) {
        return synthesize_genotype_from_parts(
            &choice.reference,
            &choice.alternates,
            default_frequencies,
            rng,
        );
    }
    synthesize_genotype(reference, default_frequencies, rng)
}

fn design_parts_for_reference(
    reference: &ReferenceVariant,
    biallelic_choices: Option<&HashMap<i64, BiallelicChoice>>,
) -> (String, Vec<String>) {
    if let Some(choice) = biallelic_choices.and_then(|choices| choices.get(&reference.rsid)) {
        return (choice.reference.clone(), choice.alternates.clone());
    }
    (
        reference.reference.clone(),
        parse_alternates(&reference.alternates),
    )
}

fn design_for_assignment(
    assignment: &OverlayAssignment,
    biallelic_choices: Option<&HashMap<i64, BiallelicChoice>>,
    genotype: &str,
) -> String {
    if let Some(choice) = biallelic_choices.and_then(|choices| choices.get(&assignment.spec.rsid)) {
        return design_for(Some(&choice.reference), &choice.alternates, genotype);
    }
    design_for(
        assignment.spec.reference.as_deref(),
        &assignment.spec.alternates,
        genotype,
    )
}

#[allow(clippy::too_many_arguments)]
fn write_row(
    writer: &mut BufWriter<File>,
    format: SyntheticFormat,
    sample_id: &str,
    rsid_label: &str,
    illumina_snp_name: &str,
    chromosome: &str,
    position: i64,
    genotype: &str,
    design: &str,
    rng: &mut StdRng,
) -> Result<()> {
    match format {
        SyntheticFormat::Ddna => {
            let is_no_call = genotype
                .trim()
                .chars()
                .all(|c| matches!(c, '-' | '.' | 'N' | 'n' | '0'));
            if is_no_call {
                let lrr = rng.gen_range(-1.5..=0.5);
                return writeln!(
                    writer,
                    "{}\t{}\t{}\t{}\t0\t0\t{:.4}",
                    rsid_label, chromosome, position, genotype, lrr
                )
                .with_context(|| format!("write row for {}", rsid_label));
            }
            let gs = rng.gen_range(0.2..=1.0);
            let baf = rng.gen_range(0.0..=1.0);
            let lrr = rng.gen_range(-0.5..=0.5);
            writeln!(
                writer,
                "{}\t{}\t{}\t{}\t{:.4}\t{:.3}\t{:.4}",
                rsid_label, chromosome, position, genotype, gs, baf, lrr
            )
            .with_context(|| format!("write row for {}", rsid_label))
        }
        SyntheticFormat::Illumina => {
            let (a1, a2) = split_alleles(genotype);
            let design_ref = design.trim_start_matches('[').chars().next().unwrap_or('N');
            let ab1 = ab_code(&a1, design_ref);
            let ab2 = ab_code(&a2, design_ref);
            // Synthetic data is internally consistent on the plus strand;
            // Top/Forward/Plus/Design all carry the same call (Strand '+').
            writeln!(
                writer,
                "{sid}\t\t{snp}\t{design}\t{chr}\t{pos}\t\
                 {a1}\t{a2}\t{a1}\t{a2}\t{a1}\t{a2}\t{a1}\t{a2}\t{ab1}\t{ab2}\t+",
                sid = sample_id,
                snp = illumina_snp_name,
                design = design,
                chr = chromosome,
                pos = position,
                a1 = a1,
                a2 = a2,
                ab1 = ab1,
                ab2 = ab2,
            )
            .with_context(|| format!("write row for {}", rsid_label))
        }
    }
}

fn load_overlay_specs(args: &SyntheticArgs) -> Result<Option<Vec<OverlaySpec>>> {
    let json_source = match (&args.variants_file, &args.variants_json) {
        (Some(_), Some(_)) => {
            bail!("Use only one of --variants-file or --variants-json");
        }
        (Some(path), None) => Some(
            std::fs::read_to_string(path)
                .with_context(|| format!("Read variants JSON file {:?}", path))?,
        ),
        (None, Some(inline)) => Some(inline.clone()),
        (None, None) => None,
    };

    let Some(raw_json) = json_source else {
        return Ok(None);
    };

    let root: OverlayRoot =
        serde_json::from_str(&raw_json).context("Parse variants JSON payload")?;
    let mut specs = Vec::new();
    for group in root.groups.values() {
        if let Some(alt_frequency) = group.alt_frequency {
            if !(0.0..=1.0).contains(&alt_frequency) {
                bail!("group alt_frequency must be between 0 and 1");
            }
        }
        let group_freqs = if let Some(raw) = group.genotype_frequencies {
            Some(raw)
        } else if let Some(alt_frequency) = group.alt_frequency {
            let frequencies = GenotypeFrequencies::from_alt_frequency(alt_frequency)?;
            Some(GenotypeFrequenciesRaw {
                het: frequencies.het,
                hom_alt: frequencies.hom_alt,
            })
        } else {
            None
        };
        for variant in &group.variants {
            let mut variant = variant.clone();
            if variant.genotype_frequencies.is_none() {
                variant.genotype_frequencies = group_freqs;
            }
            specs.push(OverlaySpec::from_raw(&variant)?);
        }
    }
    Ok(Some(specs))
}

fn load_illumina_probe_manifest(
    args: &SyntheticArgs,
    store: &StatsStore,
) -> Result<Option<IlluminaProbeManifest>> {
    if args.format != SyntheticFormat::Illumina || args.clean_illumina_rsids {
        return Ok(None);
    }
    if let Some(path) = args.illumina_probe_manifest.as_deref() {
        return load_illumina_probe_manifest_tsv(path);
    }
    let probes = store
        .all_illumina_probe_manifest()?
        .into_iter()
        .map(IlluminaProbe::from_record)
        .collect::<Vec<_>>();
    let clean_rsids = store.all_illumina_clean_rsids()?;
    if probes.is_empty() {
        Ok(None)
    } else {
        Ok(Some(IlluminaProbeManifest {
            probes,
            clean_rsids: Some(clean_rsids),
        }))
    }
}

fn load_illumina_probe_manifest_tsv(path: &Path) -> Result<Option<IlluminaProbeManifest>> {
    let mut reader = csv::ReaderBuilder::new()
        .delimiter(b'\t')
        .flexible(true)
        .from_path(path)
        .with_context(|| format!("Open Illumina probe manifest {:?}", path))?;
    let mut probes = Vec::new();
    for row in reader.deserialize::<IlluminaProbeRaw>() {
        let raw = row.with_context(|| format!("Parse Illumina probe manifest {:?}", path))?;
        let probe = IlluminaProbe::from_raw(raw)?;
        probes.push(probe);
    }
    if probes.is_empty() {
        bail!("Illumina probe manifest {:?} has no probes", path);
    }
    Ok(Some(IlluminaProbeManifest {
        probes,
        clean_rsids: None,
    }))
}

fn synthesize_genotype_from_design(
    design: &str,
    default_frequencies: &GenotypeFrequencies,
    rng: &mut StdRng,
) -> String {
    let alleles = design_alleles(design);
    if alleles.is_empty() {
        return "--".to_string();
    }
    let reference = alleles[0].clone();
    let alternates = alleles.into_iter().skip(1).collect::<Vec<_>>();
    synthesize_genotype_from_parts(&reference, &alternates, default_frequencies, rng)
}

fn design_alleles(design: &str) -> Vec<String> {
    design
        .trim()
        .trim_start_matches('[')
        .trim_end_matches(']')
        .split('/')
        .map(|allele| allele.trim())
        .filter(|allele| !allele.is_empty() && *allele != "-")
        .map(|allele| allele.to_string())
        .collect()
}

fn prepare_overlay_assignments(
    overlays: &[OverlaySpec],
    default_frequencies: GenotypeFrequencies,
    biallelic_choices: Option<&HashMap<i64, BiallelicChoice>>,
    rng: &mut StdRng,
) -> Result<HashMap<i64, OverlayAssignment>> {
    let mut assignments: HashMap<i64, OverlayAssignment> = HashMap::new();
    for spec in overlays {
        let genotype = spec.random_genotype(
            default_frequencies,
            biallelic_choices.and_then(|choices| choices.get(&spec.rsid)),
            rng,
        )?;
        if assignments
            .insert(
                spec.rsid,
                OverlayAssignment {
                    spec: spec.clone(),
                    genotype,
                },
            )
            .is_some()
        {
            bail!("Duplicate overlay rsid detected: {}", spec.rsid);
        }
    }
    Ok(assignments)
}

#[derive(Debug, Clone)]
struct OverlaySpec {
    rsid: i64,
    chromosome: String,
    position: i64,
    genotype_options: Vec<String>,
    reference: Option<String>,
    alternates: Vec<String>,
    genotype_frequencies: Option<GenotypeFrequenciesRaw>,
}

impl OverlaySpec {
    fn from_raw(raw: &OverlayVariantRaw) -> Result<Self> {
        let rsid = raw
            .rsid
            .trim()
            .trim_start_matches("rs")
            .parse::<i64>()
            .with_context(|| format!("parse overlay rsid {}", raw.rsid))?;
        let genotype_options = if let Some(options) = &raw.genotypes {
            if options.is_empty() {
                bail!("Variant {} has empty genotype list", raw.rsid);
            }
            options.clone()
        } else {
            Vec::new()
        };

        let (reference, alternates, genotype_frequencies) = if genotype_options.is_empty() {
            let Some(reference) = &raw.reference else {
                bail!(
                    "Variant {} must specify either genotypes or reference/alternates",
                    raw.rsid
                );
            };
            let alternates = raw.alternates.clone().unwrap_or_default();
            (
                Some(reference.clone()),
                alternates,
                raw.genotype_frequencies,
            )
        } else {
            (None, Vec::new(), None)
        };

        Ok(Self {
            rsid,
            chromosome: raw.chromosome.clone(),
            position: raw.position,
            genotype_options,
            reference,
            alternates,
            genotype_frequencies,
        })
    }

    fn random_genotype(
        &self,
        default_frequencies: GenotypeFrequencies,
        biallelic_choice: Option<&BiallelicChoice>,
        rng: &mut StdRng,
    ) -> Result<String> {
        if self.genotype_options.is_empty() {
            let Some(reference) = &self.reference else {
                bail!("Variant {} has no genotype options", self.rsid);
            };
            let frequencies = match &self.genotype_frequencies {
                Some(raw) => parse_frequencies(raw)?,
                None => default_frequencies,
            };
            if let Some(choice) = biallelic_choice {
                return Ok(synthesize_genotype_from_parts(
                    &choice.reference,
                    &choice.alternates,
                    &frequencies,
                    rng,
                ));
            }
            return Ok(synthesize_genotype_from_parts(
                reference,
                &self.alternates,
                &frequencies,
                rng,
            ));
        }
        let idx = rng.gen_range(0..self.genotype_options.len());
        Ok(self.genotype_options[idx].clone())
    }
}

struct OverlayAssignment {
    spec: OverlaySpec,
    genotype: String,
}

#[derive(Debug, Clone)]
struct FilePlan {
    path: PathBuf,
    seed: Option<u64>,
}

fn build_file_plans(template: &str, args: &SyntheticArgs) -> Result<Vec<FilePlan>> {
    let mut rng = StdRng::from_entropy();
    let mut used_ids: HashSet<u32> = HashSet::new();
    let range_size = args.id_max.saturating_sub(args.id_min).saturating_add(1) as usize;
    if args.count > range_size && template.contains("{id}") {
        bail!(
            "ID range ({:?}-{:?}) is too small for {} unique files",
            args.id_min,
            args.id_max,
            args.count
        );
    }

    let mut plans = Vec::with_capacity(args.count);
    for idx in 0..args.count {
        let id = if template.contains("{id}") {
            generate_unique_id(args.id_min, args.id_max, &mut used_ids, &mut rng)?
        } else {
            0
        };

        let (date_str, month, day) = if template.contains("{date}")
            || template.contains("{month}")
            || template.contains("{day}")
            || template.contains("{year}")
        {
            let month = rng.gen_range(args.month_min..=args.month_max);
            let day = rng.gen_range(args.day_min..=args.day_max);
            let date = NaiveDate::from_ymd_opt(args.date_year, month, day).ok_or_else(|| {
                anyhow!(
                    "Invalid date generated: {}/{}/{}",
                    month,
                    day,
                    args.date_year
                )
            })?;
            let formatted = date.format(&args.date_format).to_string();
            (formatted, month, day)
        } else {
            (String::new(), 0, 0)
        };

        let path_str =
            fill_output_template(template, idx, id, &date_str, month, day, args.date_year);
        let seed = args.seed.map(|base| base + idx as u64);
        plans.push(FilePlan {
            path: PathBuf::from(path_str),
            seed,
        });
    }

    Ok(plans)
}

fn generate_unique_id(
    min: u32,
    max: u32,
    used: &mut HashSet<u32>,
    rng: &mut StdRng,
) -> Result<u32> {
    for _ in 0..10_000 {
        let candidate = rng.gen_range(min..=max);
        if used.insert(candidate) {
            return Ok(candidate);
        }
    }
    bail!("Failed to generate unique participant ID within range");
}

fn fill_output_template(
    template: &str,
    index: usize,
    id: u32,
    date_str: &str,
    month: u32,
    day: u32,
    year: i32,
) -> String {
    let mut result = template.to_string();
    if result.contains("{index}") {
        result = result.replace("{index}", &format!("{:04}", index + 1));
    }
    if result.contains("{id}") {
        result = result.replace("{id}", &format!("{:06}", id));
    }
    if result.contains("{date}") {
        result = result.replace("{date}", date_str);
    }
    if result.contains("{month}") {
        result = result.replace("{month}", &format!("{:02}", month));
    }
    if result.contains("{day}") {
        result = result.replace("{day}", &format!("{:02}", day));
    }
    if result.contains("{year}") {
        result = result.replace("{year}", &year.to_string());
    }
    result
}

fn validate_biallelic_outputs(plans: &[FilePlan], format: SyntheticFormat) -> Result<()> {
    let mut observed: HashMap<String, HashSet<String>> = HashMap::new();
    for plan in plans {
        match format {
            SyntheticFormat::Ddna => scan_ddna_biallelic(&plan.path, &mut observed)?,
            SyntheticFormat::Illumina => scan_illumina_biallelic(&plan.path, &mut observed)?,
        }
    }

    let violations = observed
        .iter()
        .filter(|(_, alleles)| alleles.len() > 2)
        .map(|(rsid, alleles)| {
            let mut sorted = alleles.iter().cloned().collect::<Vec<_>>();
            sorted.sort();
            format!("{rsid}: {}", sorted.join(","))
        })
        .collect::<Vec<_>>();
    if !violations.is_empty() {
        bail!(
            "Biallelic validation failed for {} rsID(s): {}",
            violations.len(),
            violations
                .into_iter()
                .take(10)
                .collect::<Vec<_>>()
                .join("; ")
        );
    }
    Ok(())
}

fn scan_ddna_biallelic(
    path: &PathBuf,
    observed: &mut HashMap<String, HashSet<String>>,
) -> Result<()> {
    let file = File::open(path).with_context(|| format!("Open generated DDNA {:?}", path))?;
    for line in BufReader::new(file).lines() {
        let line = line?;
        let trimmed = line.trim();
        if trimmed.is_empty() || trimmed.starts_with('#') {
            continue;
        }
        let fields = trimmed.split_whitespace().collect::<Vec<_>>();
        if fields.len() < 4 {
            continue;
        }
        add_observed_genotype(observed, fields[0], fields[3]);
    }
    Ok(())
}

fn scan_illumina_biallelic(
    path: &PathBuf,
    observed: &mut HashMap<String, HashSet<String>>,
) -> Result<()> {
    let file = File::open(path).with_context(|| format!("Open generated Illumina {:?}", path))?;
    let mut in_data = false;
    let mut header: Option<HashMap<String, usize>> = None;
    for line in BufReader::new(file).lines() {
        let line = line?;
        let trimmed = line.trim();
        if trimmed.eq_ignore_ascii_case("[Data]") {
            in_data = true;
            continue;
        }
        if !in_data || trimmed.is_empty() {
            continue;
        }
        let fields = line.split('\t').collect::<Vec<_>>();
        if header.is_none() {
            header = Some(
                fields
                    .iter()
                    .enumerate()
                    .map(|(idx, name)| (name.trim().to_string(), idx))
                    .collect(),
            );
            continue;
        }
        let header = header.as_ref().expect("header must be set");
        let Some(snp_idx) = header.get("SNP Name") else {
            continue;
        };
        let Some(a1_idx) = header.get("Allele1 - Plus") else {
            continue;
        };
        let Some(a2_idx) = header.get("Allele2 - Plus") else {
            continue;
        };
        let Some(rsid) = fields.get(*snp_idx).map(|value| value.trim()) else {
            continue;
        };
        if rsid.is_empty() {
            continue;
        }
        if let Some(a1) = fields.get(*a1_idx) {
            add_observed_allele(observed, rsid, a1);
        }
        if let Some(a2) = fields.get(*a2_idx) {
            add_observed_allele(observed, rsid, a2);
        }
    }
    Ok(())
}

fn add_observed_genotype(
    observed: &mut HashMap<String, HashSet<String>>,
    rsid: &str,
    genotype: &str,
) {
    let (a1, a2) = split_alleles(genotype);
    add_observed_allele(observed, rsid, &a1);
    add_observed_allele(observed, rsid, &a2);
}

fn add_observed_allele(observed: &mut HashMap<String, HashSet<String>>, rsid: &str, allele: &str) {
    let allele = allele.trim();
    if allele.is_empty() || matches!(allele, "-" | "." | "N" | "n" | "0") {
        return;
    }
    observed
        .entry(rsid.to_string())
        .or_default()
        .insert(allele.to_ascii_uppercase());
}

#[derive(Debug, Deserialize)]
struct OverlayRoot {
    #[serde(flatten)]
    groups: HashMap<String, OverlayGroup>,
}

#[derive(Debug, Deserialize)]
struct OverlayGroup {
    #[serde(default)]
    alt_frequency: Option<f64>,
    #[serde(default)]
    genotype_frequencies: Option<GenotypeFrequenciesRaw>,
    variants: Vec<OverlayVariantRaw>,
}

#[derive(Debug, Deserialize, Clone)]
struct OverlayVariantRaw {
    rsid: String,
    chromosome: String,
    position: i64,
    #[serde(default)]
    genotypes: Option<Vec<String>>,
    reference: Option<String>,
    alternates: Option<Vec<String>>,
    #[serde(default)]
    genotype_frequencies: Option<GenotypeFrequenciesRaw>,
}

#[derive(Debug, Clone)]
struct IlluminaProbeManifest {
    probes: Vec<IlluminaProbe>,
    clean_rsids: Option<HashSet<i64>>,
}

#[derive(Debug, Clone)]
struct IlluminaProbe {
    probe_id: String,
    rsid: Option<i64>,
    design: String,
    chromosome: String,
    position: i64,
}

#[derive(Debug, Deserialize)]
struct IlluminaProbeRaw {
    probe_id: String,
    #[serde(default)]
    rsid: String,
    design: String,
    chrom: String,
    pos: i64,
}

impl IlluminaProbe {
    fn from_record(record: IlluminaProbeRecord) -> Self {
        Self {
            probe_id: record.probe_id,
            rsid: record.rsid,
            design: record.design,
            chromosome: record.chromosome,
            position: record.position,
        }
    }

    fn from_raw(raw: IlluminaProbeRaw) -> Result<Self> {
        let rsid = parse_optional_rsid(&raw.rsid)
            .with_context(|| format!("parse manifest rsid {}", raw.rsid))?;
        Ok(Self {
            probe_id: raw.probe_id,
            rsid,
            design: raw.design,
            chromosome: raw.chrom,
            position: raw.pos,
        })
    }
}

fn parse_optional_rsid(value: &str) -> Result<Option<i64>> {
    let value = value.trim();
    if value.is_empty() {
        return Ok(None);
    }
    Ok(Some(value.trim_start_matches("rs").parse::<i64>()?))
}

fn synthesize_genotype(
    reference: &ReferenceVariant,
    default_frequencies: &GenotypeFrequencies,
    rng: &mut StdRng,
) -> String {
    let alt_list = parse_alternates(&reference.alternates);
    synthesize_genotype_from_parts(&reference.reference, &alt_list, default_frequencies, rng)
}

fn parse_frequencies(raw: &GenotypeFrequenciesRaw) -> Result<GenotypeFrequencies> {
    GenotypeFrequencies::from_parts(raw.het, raw.hom_alt)
}

fn parse_alternates(alternates: &str) -> Vec<String> {
    alternates
        .split(',')
        .map(|alt| alt.trim())
        .filter(|alt| !alt.is_empty())
        .map(|alt| alt.to_string())
        .collect()
}

fn synthesize_genotype_from_parts(
    reference: &str,
    alternates: &[String],
    frequencies: &GenotypeFrequencies,
    rng: &mut StdRng,
) -> String {
    let alt_list = alternates
        .iter()
        .map(|alt| alt.as_str())
        .filter(|alt| !alt.is_empty())
        .collect::<Vec<_>>();
    let kind = determine_variant_kind(reference, &alt_list);

    let (genotype_kind, has_alt) = select_genotype_kind(frequencies, !alt_list.is_empty(), rng);
    if !has_alt {
        return homozygous_symbol(reference, None, kind, true);
    }

    match genotype_kind {
        GenotypeKind::HomRef => homozygous_symbol(reference, None, kind, true),
        GenotypeKind::HomAlt => {
            let alt = alt_list[rng.gen_range(0..alt_list.len())];
            homozygous_symbol(reference, Some(alt), kind, false)
        }
        GenotypeKind::Het => {
            let alt = alt_list[rng.gen_range(0..alt_list.len())];
            heterozygous_symbol(reference, alt, kind, rng)
        }
    }
}

fn select_genotype_kind(
    frequencies: &GenotypeFrequencies,
    has_alt: bool,
    rng: &mut StdRng,
) -> (GenotypeKind, bool) {
    if !has_alt {
        return (GenotypeKind::HomRef, false);
    }
    let roll = rng.gen::<f64>();
    let hom_alt_threshold = frequencies.hom_ref + frequencies.het + frequencies.hom_alt;
    let kind = if roll < frequencies.hom_ref {
        GenotypeKind::HomRef
    } else if roll < frequencies.hom_ref + frequencies.het {
        GenotypeKind::Het
    } else if roll < hom_alt_threshold {
        GenotypeKind::HomAlt
    } else {
        GenotypeKind::HomRef
    };
    (kind, true)
}

fn homozygous_symbol(
    reference: &str,
    alternate: Option<&str>,
    kind: VariantKind,
    use_reference: bool,
) -> String {
    let symbol = match kind {
        VariantKind::Snp => {
            if use_reference {
                reference.to_string()
            } else {
                alternate.unwrap_or(reference).to_string()
            }
        }
        VariantKind::Mnv => {
            let allele = if use_reference {
                reference
            } else {
                alternate.unwrap_or(reference)
            };
            allele
                .chars()
                .next()
                .map(|c| c.to_string())
                .unwrap_or_else(|| "?".into())
        }
        VariantKind::Insertion => {
            if use_reference {
                "D".into()
            } else {
                "I".into()
            }
        }
        VariantKind::Deletion => {
            if use_reference {
                "I".into()
            } else {
                "D".into()
            }
        }
    };
    format!("{symbol}{symbol}")
}

fn heterozygous_symbol(
    reference: &str,
    alternate: &str,
    kind: VariantKind,
    rng: &mut StdRng,
) -> String {
    let (ref_symbol, alt_symbol) = match kind {
        VariantKind::Snp => (reference.to_string(), alternate.to_string()),
        VariantKind::Mnv => (
            reference
                .chars()
                .next()
                .map(|c| c.to_string())
                .unwrap_or_else(|| "?".into()),
            alternate
                .chars()
                .next()
                .map(|c| c.to_string())
                .unwrap_or_else(|| "?".into()),
        ),
        VariantKind::Insertion => ("D".into(), "I".into()),
        VariantKind::Deletion => ("I".into(), "D".into()),
    };

    if rng.gen::<bool>() {
        format!("{ref_symbol}{alt_symbol}")
    } else {
        format!("{alt_symbol}{ref_symbol}")
    }
}

fn apply_no_call(
    genotype: &str,
    no_call_frequency: f64,
    no_call_token: &str,
    rng: &mut StdRng,
) -> String {
    if no_call_frequency > 0.0 && rng.gen::<f64>() < no_call_frequency {
        if no_call_token.len() == 1 {
            format!("{0}{0}", no_call_token)
        } else {
            no_call_token.to_string()
        }
    } else {
        genotype.to_string()
    }
}

fn determine_variant_kind(reference: &str, alts: &[&str]) -> VariantKind {
    let ref_len = reference.len();
    if alts.is_empty() {
        return VariantKind::Snp;
    }
    if ref_len == 1 && alts.iter().all(|alt| alt.len() == 1) {
        return VariantKind::Snp;
    }
    let first_alt_len = alts[0].len();
    let all_alt_same_len = alts.iter().all(|alt| alt.len() == ref_len);
    if all_alt_same_len && ref_len > 1 {
        VariantKind::Mnv
    } else if first_alt_len > ref_len {
        VariantKind::Insertion
    } else if first_alt_len < ref_len {
        VariantKind::Deletion
    } else {
        VariantKind::Mnv
    }
}

#[derive(Debug, Deserialize, Clone, Copy)]
struct GenotypeFrequenciesRaw {
    het: f64,
    hom_alt: f64,
}

#[derive(Debug, Clone, Copy)]
struct GenotypeFrequencies {
    hom_ref: f64,
    het: f64,
    hom_alt: f64,
}

impl GenotypeFrequencies {
    fn from_alt_frequency(alt_frequency: f64) -> Result<Self> {
        if !(0.0..=1.0).contains(&alt_frequency) {
            bail!("--alt-frequency must be between 0 and 1");
        }
        let ref_frequency = 1.0 - alt_frequency;
        Ok(Self {
            hom_ref: ref_frequency * ref_frequency,
            het: 2.0 * ref_frequency * alt_frequency,
            hom_alt: alt_frequency * alt_frequency,
        })
    }

    fn from_parts(het: f64, hom_alt: f64) -> Result<Self> {
        if !(0.0..=1.0).contains(&het) {
            bail!("het frequency must be between 0 and 1");
        }
        if !(0.0..=1.0).contains(&hom_alt) {
            bail!("hom_alt frequency must be between 0 and 1");
        }
        let hom_ref = 1.0 - het - hom_alt;
        if hom_ref < 0.0 {
            bail!("het + hom_alt must be <= 1.0");
        }
        Ok(Self {
            hom_ref,
            het,
            hom_alt,
        })
    }
}

#[derive(Debug, Clone, Copy)]
enum GenotypeKind {
    HomRef,
    Het,
    HomAlt,
}

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

    fn is_het(genotype: &str) -> bool {
        let mut chars = genotype.chars();
        match (chars.next(), chars.next(), chars.next()) {
            (Some(a), Some(b), None) => a != b,
            _ => false,
        }
    }

    #[test]
    fn alt_frequency_is_interpreted_as_hwe_allele_frequency() {
        let frequencies = GenotypeFrequencies::from_alt_frequency(0.2).unwrap();

        assert!((frequencies.hom_ref - 0.64).abs() < 1e-12);
        assert!((frequencies.het - 0.32).abs() < 1e-12);
        assert!((frequencies.hom_alt - 0.04).abs() < 1e-12);
    }

    #[test]
    fn synthesized_genotypes_have_expected_hwe_heterozygote_rate() {
        let frequencies = GenotypeFrequencies::from_alt_frequency(0.5).unwrap();
        let alternates = vec!["G".to_string()];
        let mut rng = StdRng::seed_from_u64(42);
        let sample_count = 20_000usize;

        let mut het_count = 0usize;
        let mut hom_alt_count = 0usize;
        for _ in 0..sample_count {
            let genotype = synthesize_genotype_from_parts("A", &alternates, &frequencies, &mut rng);
            if is_het(&genotype) {
                het_count += 1;
            }
            if genotype == "GG" {
                hom_alt_count += 1;
            }
        }

        let het_rate = het_count as f64 / sample_count as f64;
        let hom_alt_rate = hom_alt_count as f64 / sample_count as f64;
        assert!(
            (het_rate - 0.5).abs() < 0.03,
            "het_rate {het_rate} should be close to 0.5"
        );
        assert!(
            (hom_alt_rate - 0.25).abs() < 0.03,
            "hom_alt_rate {hom_alt_rate} should be close to 0.25"
        );
    }

    #[test]
    fn explicit_overlay_genotypes_override_sampled_frequencies() {
        let overlay = OverlaySpec {
            rsid: 123,
            chromosome: "1".to_string(),
            position: 456,
            genotype_options: vec!["TT".to_string()],
            reference: None,
            alternates: Vec::new(),
            genotype_frequencies: None,
        };
        let default_frequencies = GenotypeFrequencies::from_alt_frequency(0.0).unwrap();
        let mut rng = StdRng::seed_from_u64(7);

        assert_eq!(
            overlay
                .random_genotype(default_frequencies, None, &mut rng)
                .unwrap(),
            "TT"
        );
    }

    #[test]
    fn biallelic_choice_limits_multiallelic_reference_to_first_alt() {
        let reference = ReferenceVariant {
            rsid: 9701872,
            chromosome: "1".to_string(),
            position: 632828,
            reference: "T".to_string(),
            alternates: "C,G".to_string(),
        };
        let choices = build_biallelic_choices(std::slice::from_ref(&reference), &[]).unwrap();
        let frequencies = GenotypeFrequencies::from_alt_frequency(1.0).unwrap();
        let mut rng = StdRng::seed_from_u64(1);

        let genotype =
            synthesize_reference_genotype(&reference, &frequencies, Some(&choices), &mut rng);

        assert_eq!(genotype, "CC");
        let (design_reference, design_alts) =
            design_parts_for_reference(&reference, Some(&choices));
        assert_eq!(
            design_for(Some(&design_reference), &design_alts, &genotype),
            "[T/C]"
        );
    }
}
