pub const TEST_VCF: &str = "tests/test.vcf.gz";

    fn create_test_matrix(nsamples: usize, variants: usize, step_size: usize) -> PhasedMatrix {

        let variant_idx = variants / 2;

        let mut matrix = Array2::<u8>::from_elem((nsamples * 2, variants), 0);

        assert!(variants > nsamples * step_size * 2);

        let pos_coords: Vec<i64> = (1..variants+1).map(|i| i as i64).collect();

        let phase_coords: Vec<bool> = (0..variants).map(|_| true).collect();

        let ref_coords: Vec<String> = (1..variants+1)

            .map(|i| if i % 2 == 0 { "G".to_string() } else { "A".to_string() })

            .collect();

        let alt_coords: Vec<String> = (1..variants+1)

            .map(|i| if i % 2 == 0 { "T".to_string() } else { "C".to_string() })

            .collect();

        let samples: Vec<String> = (1..nsamples+1).cycle().take(nsamples * 2)

            .map(|s| format!("SAMPLE{s}"))

            .collect();

        let coords: Vec<Coord> = (0..variants)

            .map(|i| Coord {

                contig: String::from("chr9"),

                phase: phase_coords[i],

                pos: pos_coords[i],

                reference: ref_coords[i].clone(),

                alt: alt_coords[i].clone(),

            }).collect();

        for i in 0..nsamples * 2 {

            if i >= nsamples {

                matrix[[i, variant_idx]] = 1;

            }

            matrix[[i, variant_idx - (i+2) * step_size]] = 1;

            matrix[[i, variant_idx + (i+2) * step_size]] = 1;

        for i in (0..nsamples).step_by(2) {

            let (row1, row2) = matrix.multi_slice_mut((ndarray::s![i,..], ndarray::s![i+nsamples, ..]));

            ndarray::Zip::from(row1).and(row2).for_each(::std::mem::swap);

        }

        PhasedMatrix::new(variant_idx, matrix, samples, coords)

    }

    #[test]

    fn select_rows() {

        let mut vcf = create_test_matrix(14, 63, 1);

        vcf.select_rows(vec![0, 3, 20]);

        assert_eq!(vcf.samples(), &vec!["SAMPLE1".to_string(), "SAMPLE4".to_string(), "SAMPLE7".to_string()]);

    }

    #[test]

    fn select_columns_by_range() {

        let mut vcf = create_test_matrix(14, 63, 1);

        let pos = vcf.variant_idx_pos();

        let idx = vcf.variant_idx();

        vcf.select_columns_by_range(10..50);

        let pos2 = vcf.variant_idx_pos();

        let idx2 = vcf.variant_idx();

        assert_ne!(idx, idx2);

        assert_eq!(pos, pos2);

    }

    #[test]

    fn select_columns() {

        let mut vcf = create_test_matrix(14, 63, 1);

        vcf.select_columns_by_idx(&mut [0, 3, 20]);

        assert_eq!(vcf.coords().iter().map(|c| c.pos).collect::<Vec<i64>>(), vec![1, 4, 21]);

    }

    #[test]

    fn select_alt_carriers() {

        let nsamples = 14;

        let vcf = create_test_matrix(nsamples, 63, 1);

        let pos = vcf.variant_idx_pos();

        let vcf = hatk::utils::select_carrier_haplotypes(vcf, pos, &format!("chr9:{pos}"), &ALTS);

        let samples = vec![

            "SAMPLE1", "SAMPLE3", "SAMPLE5", "SAMPLE7", "SAMPLE9", "SAMPLE11", "SAMPLE13",

            "SAMPLE2", "SAMPLE4", "SAMPLE6", "SAMPLE8", "SAMPLE10", "SAMPLE12", "SAMPLE14"

        ];

        assert_eq!(&samples, vcf.samples());

    }

    #[test]

    fn select_ref_carriers() {

        let nsamples = 14;

        let vcf = create_test_matrix(nsamples, 63, 1);

        let pos = vcf.variant_idx_pos();

        let vcf = hatk::utils::select_carrier_haplotypes(

            vcf, pos, &format!("chr9:{pos}"), &hatk::args::Selection::OnlyRefs

        );

        let samples = vec![

            "SAMPLE2", "SAMPLE4", "SAMPLE6", "SAMPLE8", "SAMPLE10", "SAMPLE12", "SAMPLE14",

            "SAMPLE1", "SAMPLE3", "SAMPLE5", "SAMPLE7", "SAMPLE9", "SAMPLE11", "SAMPLE13"

        ];

        assert_eq!(&samples, vcf.samples());

    }

    #[test]

    fn get_sample_names_from_phased_matrix() {

        let nsamples = 14;

        let vcf = create_test_matrix(nsamples, 63, 1);

        let names = vcf.get_sample_names(&[0, 3, 20]);

        assert_eq!(names, vec!["SAMPLE1".to_string(), "SAMPLE4".to_string(), "SAMPLE7".to_string()]);

        let name = vcf.get_sample_name(0);

        assert_eq!("SAMPLE1".to_string(), name);

    }

    #[test]

    fn phased_matrix_setters_and_getters() {

        let mut vcf = create_test_matrix(14, 63, 1);

        vcf.set_coords(vec![]);

        assert_eq!(&Vec::<Coord>::new(), vcf.coords());

        vcf.set_variant_idx(99);

        assert_eq!(99, vcf.variant_idx());

        vcf.set_variant_idx(0);

        assert_eq!(0, vcf.variant_idx());

        let coords = vcf.coords_mut();

        coords.push(hatk::structs::Coord { 

            contig: "".to_string(),

            pos: 10,

            reference: "".to_string(),

            alt: "".to_string(),

            phase: false,

        });

        assert_eq!(10, vcf.variant_idx_pos());

        assert_eq!(1, vcf.ncoords());

        assert_eq!(14 * 2, vcf.nsamples());

    }

    #[test]

    fn test_by_pos_idx_getters() {

        let mut vcf = create_test_matrix(14, 63, 1);

        let coords = vcf.coords_mut();

        coords.iter_mut().for_each(|c| c.pos = c.pos.pow(2));

        assert_eq!(4, vcf.get_nearest_idx_by_pos(26));

        assert_eq!(4, vcf.get_nearest_idx_by_pos(23));

        assert_eq!(5, vcf.get_nearest_idx_by_pos(31));

        assert_eq!(vcf.ncoords(), vcf.get_nearest_idx_by_pos(9999999));

        let coord  = hatk::structs::Coord { 

            contig: "chr9".to_string(),

            pos: 25,

            reference: "A".to_string(),

            alt: "C".to_string(),

            phase: true,

        };

        assert_eq!(4, vcf.idx_by_coord(&coord).unwrap());

    }

    #[test]

    fn carrier_vs_only_longest() {

        let nsamples = 14;

        let vcf = create_test_matrix(nsamples, 63, 1);

        let pos = vcf.variant_idx_pos();

        let vcf = hatk::utils::select_carrier_haplotypes(vcf, pos, &format!("chr9:{pos}"), &ALTS);

        let vcf2 = create_test_matrix(nsamples, 63, 1);

        for row in vcf2.matrix.rows() {

            println!(

                "{}", row.to_vec().iter().fold(String::new(), |cur, nxt| format!("{cur}{nxt}") )

            );

        }

        let lengths = hatk::utils::shared_lengths_by_majority(&vcf2, vcf2.variant_idx());

        let vcf2 = hatk::utils::select_only_longest_haplotypes(&lengths, vcf2);

        assert_eq!(vcf.samples(), vcf2.samples())

    }

    #[test]

    fn select_only_longest_haplotypes() {

        let nsamples = 14;

        let vcf = create_test_matrix(nsamples, 63, 1);

        let lengths = hatk::utils::shared_lengths_by_majority(&vcf, vcf.variant_idx());

        let vcf = hatk::utils::select_only_longest_haplotypes(&lengths, vcf);

        let samples = vec![

            "SAMPLE1", "SAMPLE3", "SAMPLE5", "SAMPLE7", "SAMPLE9", "SAMPLE11", "SAMPLE13",

            "SAMPLE2", "SAMPLE4", "SAMPLE6", "SAMPLE8", "SAMPLE10", "SAMPLE12", "SAMPLE14"

        ];

        assert_eq!(&samples, vcf.samples());

    }

    #[test]

    fn select_only_longest_alleles() {

        let nsamples = 14;

        let vcf = create_test_matrix(nsamples, 63, 1);

        let shared_lengths = hatk::utils::shared_lengths_by_majority(&vcf, vcf.variant_idx());

        let only_longest_lengths = hatk::utils::select_only_longest_alleles(&shared_lengths);

        let indexes: Vec<usize> = only_longest_lengths.iter().map(|l| l.index).collect();

        let correct_indexes = vec![

            0, 15, 2, 17, 4, 19, 6,

            21, 8, 23, 10, 25, 12, 27

        ];

        assert_eq!(indexes, correct_indexes)

    }

    #[test]

    fn longest_common_haplotype() {

        let vcf = create_test_matrix(14, 63, 1);

        let pos = vcf.variant_idx_pos();

        let vcf = hatk::utils::select_carrier_haplotypes(vcf, pos, &format!("chr9:{pos}"), &ALTS);

        let range = shared_haplotype::find_shared_haplotype_phased(&vcf);

        assert_eq!((16..47), range);

        let vcf = create_test_matrix(14, 63, 1);

        let range = shared_haplotype::find_shared_haplotype_phased(&vcf);

        assert_eq!((31..31), range);

    }

    #[test]

    fn gt_matrix_to_match_matrix() {

        let file = PathBuf::from("tests/test_haplotype_carriers.csv");

        let ht = hatk::io::read_haplotype_file(file).unwrap();

        let start = ht.first().unwrap();

        let end = ht.last().unwrap();

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, Some((start.pos as u64, end.pos as u64))).unwrap();

        let vcf = hatk::subcommands::ref_ht_comparison::transform_gt_matrix_to_match_matrix(vcf, &ht);

        let homozygous_variants = vec![

            vcf.matrix.slice(ndarray::s![.., 14]).to_vec(),

            vcf.matrix.slice(ndarray::s![.., 16]).to_vec()

        ];

        for col in homozygous_variants {

            assert_eq!(vec![1; 28], col)

        for row in vcf.matrix.rows() {

            println!(

                "{}", row.to_vec().iter().fold(String::new(), |cur, nxt| format!("{cur}{nxt}") )

            );

        }

        let vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, Some((start.pos as u64, end.pos as u64))).unwrap();

        let vcf = hatk::utils::select_carrier_haplotypes(vcf, 32, COORDS, &ALTS);

        let vcf = hatk::subcommands::ref_ht_comparison::transform_gt_matrix_to_match_matrix(vcf, &ht);

        for row in vcf.matrix.rows() {

            assert_eq!((16..47).map(|_| 1).collect::<Vec<u8>>(), row.to_vec())

    }

    #[test]

    fn graph_to_svg() {

        let vcf = create_test_matrix(14, 63, 1);

        let g = hatk::subcommands::recursive_mbah::tree_petgraph(&vcf, &LocDirection::Left, &None, true).unwrap();

        let mut dg = hatk::graphs::DecayGraph::new(&g, &vcf, hatk::args::GraphArgs::default());

        dg.draw_graph().unwrap();

        svg::save("tests/results/test_run.svg", &dg.document).unwrap();

    }

    #[test]

    fn read_vcf_to_phased_matrix() {

        let vcf1 = create_test_matrix(14, 63, 1);

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let vcf2 = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None).unwrap();

        assert_eq!(vcf1.samples(), vcf2.samples());

        assert_eq!(vcf1.variant_idx(), vcf2.variant_idx());

        assert_eq!(vcf1.coords(), vcf2.coords());

    }

    #[test]

    fn vcf_info_filtering() {

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), info_limit: Some(0.9), ..Default::default() };

        let vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None).unwrap();

        assert_eq!(vcf.ncoords(), 31);

    }

    #[test]

    fn read_vcf_to_contra_matrix() {

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let (matrix, _pos_coords) = hatk::read_vcf::read_vcf_to_contra_matrix(&args, "chr9", None).unwrap();

        let mut one_counter = 0;

        for row in matrix.rows() {

            for value in row {

                assert_ne!(&2, value);

                if value == &1 {

                    one_counter += 1;

                }

        assert_eq!(70, one_counter);

    }

    #[test]

    fn read_vcf_to_missing_matrix() {

        let args = StandardArgs { file: PathBuf::from("tests/test_missing.vcf.gz"), ..Default::default() };

        let vcf = hatk::read_vcf::read_vcf_to_matrix_missing_allowed(&args, "chr9", 32, None).unwrap();

        let mut missing_counter = 0;

        for row in vcf.matrix.rows() {

            for value in row {

                if value == &99 {

                    missing_counter += 1;

                }

        assert_eq!(3, missing_counter);

    }

    #[test]

    fn read_vcf_to_phased_matrix_by_coords() {

        let vcf1 = create_test_matrix(14, 63, 1);

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let vcf2 = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 16, Some((1,32))).unwrap();

        assert_eq!(&vcf1.coords().iter().take(32).cloned().collect::<Vec<Coord>>(), vcf2.coords());

    }

    #[test]

    fn read_vcf_to_phased_matrix_filter_info() {

        let args = StandardArgs { 

            file: PathBuf::from(TEST_VCF),  info_limit: Some(0.9), ..Default::default()

        };

        let vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None).unwrap();

        assert_eq!(vcf.ncoords(), 31);

    }

    #[test]

    fn error_on_not_normalized_vcf() {

        let args = StandardArgs { 

            file: PathBuf::from("tests/test_not_normalized.vcf.gz"), ..Default::default()

        };

        let res = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None);

        assert!(res.is_err());

        let res = hatk::read_vcf::read_vcf_to_matrix_missing_allowed(&args, "chr9", 32, None);

        assert!(res.is_err());

        let res = hatk::read_vcf::read_vcf_to_contra_matrix(&args, "chr9", None);

        assert!(res.is_err());

    }

    #[test]

    fn contrahomozygote() {

        let test = PathBuf::from(TEST_VCF);

        let index_vec: Vec<usize> = (0..14).collect();

        let mut reader = hatk::read_vcf::get_reader(&test, "chr9", None).unwrap();

        let records = reader.records();

        let (merged_gt_vec, _) = read_vcf_to_rwc_vec(records, &index_vec, None).unwrap();

        assert_eq!(vec![0; 63], merged_gt_vec);

        let test_contra = PathBuf::from("tests/test_contrahomozygote.vcf.gz");

        let mut reader = hatk::read_vcf::get_reader(&test_contra, "chr9", None).unwrap();

        let records = reader.records();

        let (merged_gt_vec, _) = read_vcf_to_rwc_vec(records, &index_vec, None).unwrap();

        let mut test_vec = vec![0; 63];

        test_vec[2] = 1;

        test_vec[60] = 1;

        assert_eq!(test_vec, merged_gt_vec);

    }

    #[test]

    fn contrahomozygote_sequence_length() {

        let index_vec: Vec<usize> = (0..14).collect();

        let test_contra = PathBuf::from("tests/test_contrahomozygote.vcf.gz");

        let mut reader = hatk::read_vcf::get_reader(&test_contra, "chr9", None).unwrap();

        let records = reader.records();

        let (merged_gt_vec, _) = read_vcf_to_rwc_vec(records, &index_vec, None).unwrap();

        let map = rwc_vec_to_btreemap(merged_gt_vec);

        let value = map.get(&3).unwrap();

        assert_eq!(value.0, 57);

    }

    #[test]

    fn check_for_haplotype() {

        let vcf = create_test_matrix(14, 63, 1);

        let ht = vec![

            HapVariant{ pos: 32, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 0 },

            HapVariant{ pos: 33, contig: "chr9".into(), reference: "A".into(), alt: "C".into(), gt: 0 }

        ];

        let matches = hatk::subcommands::check_for_haplotype::identical_haplotype_count(&vcf, &ht);

        assert_eq!(matches.len(), 14);

        let ht = vec![

            HapVariant{ pos: 32, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 0 },

            HapVariant{ pos: 33, contig: "chr9".into(), reference: "A".into(), alt: "C".into(), gt: 0 },

            HapVariant{ pos: 34, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 1 },

        ];

        let matches = hatk::subcommands::check_for_haplotype::identical_haplotype_count(&vcf, &ht);

        assert_eq!(matches.len(), 1);

        let ht = vec![

            HapVariant{ pos: 32, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 0 },

            HapVariant{ pos: 34, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 1 },

        ];

        let matches = hatk::subcommands::check_for_haplotype::identical_haplotype_count(&vcf, &ht);

        assert_eq!(matches.len(), 1);

        let ht = vec![

            HapVariant{ pos: 32, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 0 },

            HapVariant{ pos: 34, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 1 },

            // This does not exist in the test matrix so it will not be taken into account

            HapVariant{ pos: 35, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 1 }

        ];

        let matches = hatk::subcommands::check_for_haplotype::identical_haplotype_count(&vcf, &ht);

        assert_eq!(matches.len(), 1);

        let ht = vec![

            HapVariant{ pos: 32, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 0 },

            HapVariant{ pos: 34, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 1 },

            HapVariant{ pos: 35, contig: "chr9".into(), reference: "A".into(), alt: "C".into(), gt: 1 }

        ];

        let matches = hatk::subcommands::check_for_haplotype::identical_haplotype_count(&vcf, &ht);

        assert_eq!(matches.len(), 0);

        let ht = vec![

            HapVariant{ pos: 9999, contig: "chr9".into(), reference: "G".into(), alt: "T".into(), gt: 0 },

        ];

        let matches = hatk::subcommands::check_for_haplotype::identical_haplotype_count(&vcf, &ht);

        assert!(matches.is_empty());

    }

    #[test]

    fn set_vcf_variable_data() {

        let path = PathBuf::from("tests/clinical_data.csv");

        let df = hatk::io::read_variable_data_file(path).unwrap();

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let mut vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None).unwrap();

        vcf.set_variable_data(df).unwrap();

        // Make sure NA is calculated right in the average

        let means = vcf.get_variable_data_mean(&[0,1,2,3], &["aoo".to_string(), "dur".to_string()]).unwrap().unwrap();

        assert_eq!(60.75, means[0]);

        assert_eq!(7.0, means[1]);

        let means = vcf.get_variable_data_mean(

            &[0,1,2,3,4,5,6,7,8,9,10,11,12,13],

            &["aoo".to_string(), "dur".to_string()

        ]).unwrap().unwrap();

        assert_eq!("62.69", format!("{:.2}", means[0]));

        assert_eq!("4.00", format!("{:.2}", means[1]));

        let result = vcf.get_variable_data_mean(&[13], &["aoo".to_string()]);

        assert!(result.is_err());

        let vecs = vcf.get_variable_data_vecs(&[0,1], &["aoo".to_string(), "dur".to_string()]).unwrap().unwrap();

        assert_eq!(vec![88.0, 49.0], vecs[0]);

        assert_eq!(vec![2.0, 10.0], vecs[1]);

    }

    fn standard_args(selection: Selection) -> StandardArgs {

        StandardArgs {

            file: PathBuf::from(TEST_VCF),

            output: PathBuf::from("tests/results"),

            coords: String::from(COORDS),

            selection,

            ..Default::default()

        }

    }

    fn clap_standard_args(select: hatk::clap::ClapSelection) -> hatk::clap::ClapStandardArgs {

        hatk::clap::ClapStandardArgs {

            file: PathBuf::from(TEST_VCF),

            outdir: PathBuf::from("tests/results"),

            coords: String::from(COORDS),

            select,

            ..Default::default()

        }

    }

    #[test]

    fn sample_names_subcommand() {

        let path = PathBuf::from(TEST_VCF);

        let ids = get_sample_names(path).unwrap();

        let vec: Vec<_> = (1..15).map(|v| format!("SAMPLE{v}")).collect();

        assert_eq!(vec, ids);

        let path = PathBuf::from("tests/samples.fam");

        let ids = get_sample_names(path).unwrap();

        let vec = vec!["SAMPLE1", "SAMPLE2", "SAMPLE3", "SAMPLE4",];

        assert_eq!(vec, ids);

        let path = PathBuf::from("tests/samples.fam");

        let res = read_sample_names::run(path);

        assert!(res.is_ok());

    }

    #[test]

    fn shared_haplotype() {

        let args = standard_args(Selection::OnlyLongest);

        shared_haplotype::run(args, false).unwrap();

        let res = std::fs::read_to_string("tests/results/shared_haplotype_only_longest.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn shared_haplotype_unphased() {

        let args = standard_args(Selection::Unphased);

        shared_haplotype::run(args, false).unwrap();

        let res = std::fs::read_to_string("tests/results/shared_haplotype_rwc.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn mbah() {

        let args = standard_args(Selection::All);

        mbah::run(args, GraphArgs::default(), None, None).unwrap();

        let res = std::fs::read_to_string("tests/results/mbah.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn rwc() {

        let file = PathBuf::from(TEST_CONTRA);

        let output = PathBuf::from("tests/results");

        rwc::run(file, None, output, None, false, 10, None).unwrap();

        let res = std::fs::read_to_string("tests/results/genome_wide_rwc.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn compare_haplotypes() {

        let args = standard_args(Selection::All);

        mbah::run(args, GraphArgs::default(), None, None).unwrap();

        let args = standard_args(Selection::OnlyLongest);

        shared_haplotype::run(args, false).unwrap();

        let files = vec![

            PathBuf::from("tests/results/mbah.csv"),

            PathBuf::from("tests/results/shared_haplotype_only_longest.csv"),

        ];

        let outdir = PathBuf::from(OUTDIR);

        let cmd = hatk::clap::SubCommand::CompareHaplotypes {

            haplotypes: files,

            outdir,

            prefix: None,

            csv: true,

            hide_missing: false,

            tag_rows: false,

        };

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/ht_comparison.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    fn check_for_haplotype(selection: hatk::clap::ClapSelection) {

        let args = standard_args(Selection::OnlyLongest);

        shared_haplotype::run(args, false).unwrap();

        let args = clap_standard_args(selection);

        let cmd = hatk::clap::SubCommand::CheckForHaplotype {

            args,

            haplotype: PathBuf::from("tests/results/shared_haplotype_only_longest.csv"),

        };

        hatk::clap::run_cmd(cmd).unwrap();

    }

    #[test]

    fn check_for_haplotype_all() {

        check_for_haplotype(hatk::clap::ClapSelection::All);

        let res = std::fs::read_to_string("tests/results/haplotype_check.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn check_for_haplotype_only_longest() {

        check_for_haplotype(hatk::clap::ClapSelection::OnlyLongest);

        let res = std::fs::read_to_string("tests/results/haplotype_check_only_longest.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn check_for_haplotype_only_refs() {

        check_for_haplotype(hatk::clap::ClapSelection::OnlyRefs);

        let res = std::fs::read_to_string("tests/results/haplotype_check_only_refs.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    fn ref_ht_comparison(selection: hatk::clap::ClapSelection, mark: bool, png: bool) {

        let args = standard_args(Selection::OnlyLongest);

        shared_haplotype::run(args, false).unwrap();

        let args = clap_standard_args(selection);

        let cmd = hatk::clap::SubCommand::CompareToHaplotype {

            args,

            haplotype: PathBuf::from("tests/results/shared_haplotype_only_longest.csv"),

            threads: 8,

            decoy_samples: None,

            mark_shorter_alleles: mark,

            png,

            graph_args: hatk::clap::ClapGraphArgs::default(),

        };

        assert_eq!(8, cmd.threads());

        hatk::clap::run_cmd(cmd).unwrap();

    }

    #[test]

    fn ref_ht_comparison_all() {

        ref_ht_comparison(hatk::clap::ClapSelection::All, true, false);

        let res = std::fs::read_to_string("tests/results/differences.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn ref_ht_comparison_only_refs() {

        ref_ht_comparison(hatk::clap::ClapSelection::OnlyRefs, false, false);

        let res = std::fs::read_to_string("tests/results/differences_only_refs.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn ref_ht_comparison_only_alts() {

        ref_ht_comparison(hatk::clap::ClapSelection::OnlyAlts, false, false);

        let res = std::fs::read_to_string("tests/results/differences_only_alts.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn ref_ht_comparison_only_longest() {

        ref_ht_comparison(hatk::clap::ClapSelection::OnlyLongest, false, false);

        let res = std::fs::read_to_string("tests/results/differences_only_longest.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn ref_ht_comparison_png() {

        ref_ht_comparison(hatk::clap::ClapSelection::All, false, false);

        let res = std::fs::read_to_string("tests/results/differences.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn ref_ht_comparison_fail() {

        let args = clap_standard_args(hatk::clap::ClapSelection::Unphased);

        let cmd = hatk::clap::SubCommand::CompareToHaplotype {

            args,

            haplotype: PathBuf::from("tests/results/shared_haplotype_only_longest.csv"),

            threads: 8,

            decoy_samples: None,

            mark_shorter_alleles: false,

            png: false,

            graph_args: hatk::clap::ClapGraphArgs::default(),

        };

        assert!(hatk::clap::run_cmd(cmd).is_err());

    }

    #[test]

    fn mrca() {

        let args = clap_standard_args(hatk::clap::ClapSelection::OnlyLongest);

        let cmd = hatk::clap::SubCommand::Mrca {

            args,

            recombination_rates: PathBuf::from("tests/recombination_rates.tsv"),

        };

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/mrca_gamma_method.txt").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn mrca_scan() {

        let args = clap_standard_args(hatk::clap::ClapSelection::All);

        let cmd = hatk::clap::SubCommand::MrcaScan {

            args,

            recombination_rates: PathBuf::from("tests/recombination_rates.tsv"),

            threads: 8,

            start: None,

            stop: None,

            no_csv: false,

            plot: true,

            graph_args: hatk::clap::ClapGraphArgs::default(),

            mark_centromere: false,

            step_size: 5,

        };

        assert_eq!(8, cmd.threads());

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/mrca_scan.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn mrca_scan_only_longest_mark_centromere() {

        let args = clap_standard_args(hatk::clap::ClapSelection::OnlyLongest);

        let cmd = hatk::clap::SubCommand::MrcaScan {

            args,

            recombination_rates: PathBuf::from("tests/recombination_rates.tsv"),

            threads: 8,

            start: None,

            stop: None,

            no_csv: false,

            plot: true,

            graph_args: hatk::clap::ClapGraphArgs::default(),

            mark_centromere: true,

            step_size: 5,

        };

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/mrca_scan.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn recursive_mbah() {

        let args = clap_standard_args(hatk::clap::ClapSelection::OnlyLongest);

        let cmd = hatk::clap::SubCommand::RecursiveMBAH {

            args,

            graph_args: hatk::clap::ClapGraphArgs::default(),

            decoy_samples: None,

            variable_data: Some(PathBuf::from("tests/clinical_data.csv")),

            variables: Some(vec!["aoo".to_string()]),

            save_branch_haplotypes: true,

        };

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/recursive_mbah_left_only_longest.svg").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn recursive_mbah_all() {

        let args = clap_standard_args(hatk::clap::ClapSelection::All);

        let cmd = hatk::clap::SubCommand::RecursiveMBAH {

            args,

            graph_args: hatk::clap::ClapGraphArgs::default(),

            decoy_samples: None,

            variable_data: Some(PathBuf::from("tests/clinical_data.csv")),

            variables: Some(vec!["aoo".to_string()]),

            save_branch_haplotypes: true,

        };

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/recursive_mbah_left.svg").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn coverage() {

        let cmd = hatk::clap::SubCommand::Coverage {

            file: PathBuf::from(TEST_VCF),

            npipes: 10,

            bp_per_snp: 2,

            threads: 1,

        };

        let res = hatk::clap::run_cmd(cmd);

        // Contig has under 100 records

        assert!(res.is_err());

    }

    #[test]

    fn samples() {

        let cmd = hatk::clap::SubCommand::Samples {

            file: PathBuf::from(TEST_VCF),

        };

        hatk::clap::run_cmd(cmd).unwrap();

    }

    #[test]

    fn local_rwc_alts() {

        let args = clap_standard_args(hatk::clap::ClapSelection::OnlyAlts);

        let cmd = hatk::clap::SubCommand::SampleSharedHaplotype {

            args,

            threads: 8,

            niters: 100,

            sample_size: 4,

        };

        hatk::clap::run_cmd(cmd).unwrap();

    }

    #[test]

    fn local_rwc_unphased() {

        let args = clap_standard_args(hatk::clap::ClapSelection::Unphased);

        let cmd = hatk::clap::SubCommand::SampleSharedHaplotype {

            args,

            threads: 8,

            niters: 100,

            sample_size: 4,

        };

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/random_sampled_shared_haplotype_rwc.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn sampled_rwc() {

        let cmd = hatk::clap::SubCommand::SampleRWC {

            file: PathBuf::from(TEST_VCF),

            merge_rule: hatk::clap::ClapMergeRule::Start,

            prefix: None,

            samples: None,

            outdir: PathBuf::from("tests/results"),

            info_limit: None,

            top_markers_n: 20,

            threads: 8,

            niters: 100,

            sort: true,

            sample_size: 4,

            min_marker_length: 2,

        };

        hatk::clap::run_cmd(cmd).unwrap();

        let res = std::fs::read_to_string("tests/results/random_sampled_rwc_start.csv").unwrap();

        insta::assert_yaml_snapshot!(res);

    }

    #[test]

    fn matrix_to_graph() {

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None).unwrap();

        let g = recursive_mbah::tree_petgraph(&vcf, &LocDirection::Left, &None, true).unwrap();

        let mut f = std::fs::File::create("tests/results/test.dot").unwrap();

        f.write_all(format!("{}", petgraph::dot::Dot::new(&g)).as_bytes()).unwrap();

    }

    #[test]

    fn matrix_graph_png() {

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let gargs = GraphArgs { ..Default::default() };

        let vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None).unwrap();

        let indexes: Vec<_> = (0..vcf.samples().len()).collect();

        hatk::graphs::matrix_graph::matrix_graph_png(&vcf, gargs, false, None, &indexes);

    }

    #[test]

    fn graph_ancestral_haplotype() {

        let args = StandardArgs { file: PathBuf::from(TEST_VCF), ..Default::default() };

        let vcf = hatk::read_vcf::read_vcf_to_matrix(&args, "chr9", 32, None).unwrap();

        let g = recursive_mbah::tree_petgraph(&vcf, &LocDirection::Left, &None, true).unwrap();

        let dg = DecayGraph::new(&g, &vcf, GraphArgs::default());

        let ht = dg.find_ancestral_haplotype();

        let g2 = recursive_mbah::tree_petgraph(&vcf, &LocDirection::Right, &None, true).unwrap();

        let dg2 = DecayGraph::new(&g2, &vcf, GraphArgs::default());

        let ht2 = dg2.find_ancestral_haplotype();

        assert_eq!(ht.len() + ht2.len(), 55);

    }

    #[test]

    fn read_haplotype_file() {

        let file = PathBuf::from("tests/test_haplotype_carriers.csv");

        let ht = hatk::io::read_haplotype_file(file).unwrap();

        let pos: Vec<i64> = (17..48).collect();

        assert_eq!(pos, ht.iter().map(|h| h.pos).collect::<Vec<i64>>());

        let file = PathBuf::from("tests/test_haplotype_unsorted.csv");

        let result = hatk::io::read_haplotype_file(file);

        assert!(result.is_err());

    }

    #[test]

    fn read_recombination_rates() {

        let rates = PathBuf::from("tests/recombination_rates.tsv");

        let btree = hatk::io::read_recombination_file(rates).unwrap();

        let pos: Vec<_> = (1..64).collect();

        for (i, (p, _)) in btree.iter().enumerate() {

            assert_eq!(&pos[i], p);

        let rates = PathBuf::from("tests/recombination_rates_with_header.tsv");

        let res = hatk::io::read_recombination_file(rates);

        assert!(res.is_err());

    }

    #[test]

    fn read_samples_file() {

        let rates = PathBuf::from("tests/samples.ids");

        let wanted = hatk::io::read_sample_ids(&Some(rates)).unwrap().unwrap();

        let vec = vec!["SAMPLE1", "SAMPLE2", "SAMPLE3", "SAMPLE10",

            "#*@#_)i%*#",

            "foo_foosa_23e41_das",

            "!!!!!!!",

            "#######"

        ];

        assert_eq!(vec, wanted);

        let path = PathBuf::from("tests/does_not_exist.ids");

        let res = hatk::io::read_sample_ids(&Some(path));

        assert!(res.is_err());

    }

    #[test]

    fn read_variable_data() {

        let path = PathBuf::from("tests/clinical_data.csv");

        let df = hatk::io::read_variable_data_file(path).unwrap();

        let array = df["id"].utf8().unwrap();

        let vec: Vec<_> = array.into_iter().filter_map(|v| v).collect();

        assert_eq!(vec[0], "SAMPLE1");

        let array = df["aoo"].i64().unwrap();

        let vec: Vec<_> = array.into_iter().filter_map(|v| v).collect();

        assert_eq!(vec[0], 88);

        assert_eq!(vec[2], 58);

        let path = PathBuf::from("tests/does_not_exist.csv");

        let res = hatk::io::read_variable_data_file(path);

        assert!(res.is_err());

    }