if let Some(prefix) = &args.prefix {

        match args.selection {

            Selection::All => output.push(format!("{prefix}_{name}.{suffix}")),

        match args.selection {

            Selection::All => output.push(format!("{name}.{suffix}")),

            Selection::OnlyAlts => output.push(format!("{name}_only_alts.{suffix}")),

            Selection::OnlyRefs => output.push(format!("{name}_only_refs.{suffix}")),

            Selection::OnlyLongest => output.push(format!("{name}_only_longest.{suffix}")),

            Selection::Unphased => output.push(format!("{name}_rwc.{suffix}")),

}

    if let Some(prefix) = &args.prefix {

        match args.selection {

            Selection::All | Selection::OnlyLongest => output.push(format!("{name}.{suffix}")),

}

    if let Some(wanted) = wanted {

        for i in &wanted {

            if !samples.contains(i) {

                tracing::warn!("Wanted sample {i} is not in the vcf");

            }

        samples

            .iter()

            .enumerate()

            .filter(|(_, s)| wanted.contains(s))

            .map(|(i, _)| i)

            .collect()

        (0..samples.len()).collect()

}

pub fn select_carrier_haplotypes(

    mut vcf: PhasedMatrix,

    variant_pos: i64,

    coords: &str,

    selection: &Selection

) -> PhasedMatrix {

    assert_eq!(

        vcf.variant_idx_pos(),

    let samples_len = vcf.matrix.nrows() / 2;

    let mut haplotype_idxs = vec![];

    for first_allele_idx in 0..samples_len {

        let second_allele_idx = first_allele_idx + samples_len;

        let a = *vcf.matrix.slice(s![first_allele_idx, vcf.variant_idx]).into_scalar();

        let b = *vcf.matrix.slice(s![second_allele_idx, vcf.variant_idx]).into_scalar();

        match selection {

                match (a, b) {

                    (1, 0) => haplotype_idxs.push(first_allele_idx),

                    (0, 1) => haplotype_idxs.push(second_allele_idx),

                match (a, b) {

                    (1, 0) => haplotype_idxs.push(second_allele_idx),

                    (0, 1) => haplotype_idxs.push(first_allele_idx),

    vcf.select_rows(haplotype_idxs);

    vcf

}

pub fn shared_lengths_by_majority(vcf: &PhasedMatrix, variant_idx: usize) -> Vec<SharedLength> {

    let mut lengths: Vec<SharedLength> = (0..vcf.nsamples())

        .map(|i| SharedLength {

            index: i,

            start: None,

            end: None,

            start_pos: 0,

            end_pos: 0,

        })

        .collect();

    // Right side

    let mut mutating_indexes: Vec<usize> = (0..vcf.nsamples()).collect();

    for pos in variant_idx..vcf.matrix.ncols() {

        let submatrix = vcf.matrix.select(ndarray::Axis(0), &mutating_indexes);

        // Returns the majority of indexes

        let majority_idxs = return_majority(&submatrix, pos, &mutating_indexes);

        if majority_idxs.len() == 1 {

            lengths.iter_mut().for_each(|i| {

                if i.end.is_none() {

                    i.end = Some(pos)

                }

            });

            break;

            if mutating_indexes.len() != majority_idxs.len() {

                // Add end positions to all indexes present in the original indexes, but not in the

                // majority of indexes

                mutating_indexes.iter().for_each(|i| {

                    if !majority_idxs.contains(i) {

                        lengths[*i].end = Some(pos)

                    }

                });

            }

        mutating_indexes = majority_idxs;

    let mut mutating_indexes: Vec<usize> = (0..vcf.nsamples()).collect();

    for pos in (0..variant_idx).rev() {

        let submatrix = vcf.matrix.select(ndarray::Axis(0), &mutating_indexes);

        let majority_idxs = return_majority(&submatrix, pos, &mutating_indexes);

        if majority_idxs.len() == 1 {

            lengths.iter_mut().for_each(|i| {

                if i.start.is_none() {

                    i.start = Some(pos)

                }

            });

            break;

        } else if mutating_indexes.len() != majority_idxs.len() {

            mutating_indexes.iter().for_each(|i| {

                if !majority_idxs.contains(i) {

                    lengths[*i].start = Some(pos)

                }

            });

        }

        mutating_indexes = majority_idxs;

    lengths.iter_mut().for_each(|i| {

        if i.start.is_none() {

            i.start = Some(0)

        }

        if i.end.is_none() {

            i.end = Some(vcf.matrix.ncols() - 1)

        }

        i.start_pos = vcf.get_pos(i.start());

        i.end_pos = vcf.get_pos(i.end());

    });

    lengths

}

fn return_majority(matrix: &Array2<u8>, pos: usize, indexes: &[usize]) -> Vec<usize> {

    let zeroes: Vec<usize> = matrix

        .slice(ndarray::s![.., pos])

        .iter()

        .zip(indexes.to_owned())

        .filter_map(|(n, i)| if *n == 0 { Some(i) } else { None })

        .collect();

    let ones: Vec<usize> = matrix

        .slice(ndarray::s![.., pos])

        .iter()

        .zip(indexes.to_owned())

        .filter_map(|(n, i)| if *n == 1 { Some(i) } else { None })

        .collect();

    match zeroes.len().cmp(&ones.len()) {

        std::cmp::Ordering::Less => ones,

        std::cmp::Ordering::Greater => zeroes,

            if ones.len() != 1 {

                tracing::warn!(

            }

            ones

}

pub fn select_only_longest_haplotypes<T: AsRef<SharedLength>>(

    lengths: &Vec<T>,

    mut vcf: PhasedMatrix,

) -> PhasedMatrix {

    let mut only_longest_lengths = vec![];

    for index in 0..lengths.len() / 2 {

        let index2 = index + lengths.len() / 2;

        let allele1 = &lengths[index];

        let allele2 = &lengths[index2];

        let allele1_length = allele1.as_ref().end_pos - allele1.as_ref().start_pos;

        let allele2_length = allele2.as_ref().end_pos - allele2.as_ref().start_pos;

        match allele1_length.cmp(&allele2_length) {

            std::cmp::Ordering::Less => only_longest_lengths.push(index2),

            std::cmp::Ordering::Greater => only_longest_lengths.push(index),

    vcf.select_rows(only_longest_lengths);

    vcf

}

pub fn select_only_longest_alleles(shared_lengths: &Vec<SharedLength>) -> Vec<&SharedLength> {

    let mut only_longest_lengths = vec![];

    for index in 0..&shared_lengths.len() / 2 {

        let allele1 = &shared_lengths[index];

        let allele2 = &shared_lengths[index + shared_lengths.len() / 2];

        let allele1_length = allele1.end_pos - allele1.start_pos;

        let allele2_length = allele2.end_pos - allele2.start_pos;

        match allele1_length.cmp(&allele2_length) {

            std::cmp::Ordering::Less => only_longest_lengths.push(allele2),

            std::cmp::Ordering::Greater => only_longest_lengths.push(allele1),

    only_longest_lengths

}

pub fn two_tail_welch_t_test(x: &[f64], y: &[f64]) -> f64 {

    let n1 = x.len() as f64;

    let n2 = y.len() as f64;

    let mean1 = x.mean();

    let mean2 = y.mean();

    let variance1 = f64::powi(x.std_dev(), 2);

    let variance2 = f64::powi(y.std_dev(), 2);

    let t = (mean1 - mean2) / (variance1 / n1 + variance2 / n2).sqrt();

    let var_by_n1 = variance1 / n1;

    let var_by_n2 = variance2 / n2;

    let upper_part = f64::powi(var_by_n1 + var_by_n2, 2);

    let lower_part = f64::powi(var_by_n1, 2) / (n1 - 1.0) + f64::powi(var_by_n2, 2) / (n2 - 1.0);

    let df = upper_part / lower_part;

    let tdist = StudentsT::new(0.0, 1.0, df).unwrap();

    tdist.cdf(t) * 2.0

}

    #[test]

    fn sample_filtering() {

        let samples = vec!["foo".to_string(), "fii".to_string()];

        let wanted = vec!["foo".to_string(), "fii".to_string()];

        let sample_indexes = filter_samples(&samples, Some(wanted));

        assert_eq!(sample_indexes, vec![0,1]);

        let samples = vec!["foo".to_string(), "fii".to_string()];

        let wanted = vec!["fee".to_string(), "foo".to_string(), "fii".to_string()];

        let sample_indexes = filter_samples(&samples, Some(wanted));

        assert_eq!(sample_indexes, vec![0,1]);

        let samples = vec!["faa".to_string(), "foo".to_string(), "fii".to_string()];

        let wanted = vec!["fee".to_string(), "foo".to_string(), "fii".to_string()];

        let sample_indexes = filter_samples(&samples, Some(wanted));

        assert_eq!(sample_indexes, vec![1,2]);

        let sample_indexes = filter_samples(&samples, None);

        assert_eq!(sample_indexes, vec![0,1,2]);

    }

    #[test]

    fn test_two_tail_welch_t_test() {

        let vec1 = vec![1.0, 2.0, 1.0, 3.0];

        let vec2 = vec![1.0, 2.0, 1.0, 5.0];

        let result = two_tail_welch_t_test(&vec1, &vec2);

        assert_eq!("0.6596", &format!("{result:.4}"));

        let vec1 = vec![1.0, -2.0, 1.0, 3.0];

        let vec2 = vec![-1.0, 2.0, 1.0, 5.0];

        let result = two_tail_welch_t_test(&vec1, &vec2);

        assert_eq!("0.5606", &format!("{result:.4}"));

        let vec1 = vec![1.0, -2.0, 1.0, 3.0];

        let vec2 = vec![-1.0, 2.0, 1.0, 5.0, 5.0, 5.0];

        let result = crate::utils::two_tail_welch_t_test(&vec1, &vec2);

        assert_eq!("0.1959", &format!("{result:.4}"));

    }

    #[test]

    fn test_return_majority() {

        let matrix = ndarray::arr2(&[

                [ 1,  1,  1],

                [ 0,  1,  0]

        ]);

        let maj = return_majority(&matrix, 0, &[5, 8]);

        assert_eq!(maj, vec![5]);

        let maj = return_majority(&matrix, 1, &[5, 8]);

        assert_eq!(maj, vec![5, 8]);

        let maj = return_majority(&matrix, 2, &[5, 8]);

        assert_eq!(maj, vec![5]);

        let matrix = ndarray::arr2(&[

                [ 1,  0,  1, 1],

                [ 0,  0,  1, 1],

                [ 1,  0,  0, 1]

        ]);

        let maj = return_majority(&matrix, 0, &[5, 8, 100]);

        assert_eq!(maj, vec![5, 100]);

        let maj = return_majority(&matrix, 1, &[5, 8, 100]);

        assert_eq!(maj, vec![5, 8, 100]);

        let maj = return_majority(&matrix, 2, &[5, 8, 100]);

        assert_eq!(maj, vec![5, 8]);

        let maj = crate::utils::return_majority(&matrix, 3, &[5, 8, 100]);

        assert_eq!(maj, vec![5, 8, 100]);

    }

    #[test]

    fn test_push_to_output() {

        let mut output = std::path::PathBuf::new();

        let args = crate::args::StandardArgs::default();

        push_to_output(&args, &mut output, "picture", "png");

        assert_eq!(output, std::path::PathBuf::from("picture.png"));

        let mut output = std::path::PathBuf::from("./foo");

        let args = crate::args::StandardArgs::default();

        push_to_output(&args, &mut output, "picture", "png");

        assert_eq!(output, std::path::PathBuf::from("./foo/picture.png"));

        let mut output = std::path::PathBuf::from("./foo");

        let args = crate::args::StandardArgs { 

            prefix: Some("nice".to_string()), 

            ..Default::default()

        };

        push_to_output(&args, &mut output, "picture", "png");

        assert_eq!(output, std::path::PathBuf::from("./foo/nice_picture.png"));

    }