pub fn run(

    args: StandardArgs,

    haplotype_path: PathBuf,

    decoy_samples: Option<PathBuf>,

    mark_shorter_alleles: bool,

    want_png: bool,

    graph_args: GraphArgs,

) -> Result<()> {

    if args.selection == Selection::Unphased {

        return Err(eyre!("Running with unphased data is not supported."))

    }

    let decoy_samples = read_sample_ids(&decoy_samples)?;

    let ht = read_haplotype_file(haplotype_path)?;

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

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

    let (contig, variant_pos) = parse_snp_coord(&args.coords)?;

    let img_type = match want_png {

        false => "svg",

    let mut img_output = args.output.clone();

    let name = match mark_shorter_alleles {

        true => "differences_shorter_alleles_marked",

        false => "differences"

    push_to_output(&args, &mut img_output, name, img_type);

    svg::save(&img_output, &svg::Document::new())

        .wrap_err(eyre!("Failed writing to {:?}", img_output))?;

    let mut output = args.output.clone();

    push_to_output(&args, &mut output, "differences", "csv");

    let mut writer = open_csv_writer(output)?;

    let mut only_longest: Option<Vec<usize>> = None;

    tracing::debug!("Reading into a sample-variant matrix");

    let vcf = match args.selection {

            let vcf = read_vcf_to_matrix(

                &args,

                contig,

                variant_pos,

                Some((start.pos as u64, end.pos as u64)),

            )?;

            select_carrier_haplotypes(vcf, variant_pos, &args.coords, &args.selection)

            let vcf = read_vcf_to_matrix(&args, contig, variant_pos, None)?;

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

            let only_longest_lengths = select_only_longest_alleles(&shared_lengths);

            only_longest = Some(only_longest_lengths.iter().map(|s| s.index).collect());

            let mut vcf = select_only_longest_haplotypes(&shared_lengths, vcf);

            // Use start and end from the haplotype to select columns from the matrix by range

            vcf.select_columns_by_range(

                vcf.idx_by_hapvariant(start).unwrap()..vcf.idx_by_hapvariant(end).unwrap() + 1,

            );

            vcf

        Selection::All => read_vcf_to_matrix(

            &args,

            contig,

            variant_pos,

            Some((start.pos as u64, end.pos as u64)),

        )?,

    tracing::debug!("Finished reading the matrix");

    let vcf = transform_gt_matrix_to_match_matrix(vcf, &ht);

    tracing::debug!("Finished transforming to match matrix");

    let mut npy_output = args.output.clone();

    push_to_output(&args, &mut npy_output, "differences", "npy");

    ndarray_npy::write_npy(npy_output, &vcf.matrix)?;

    let index_order = sort_indexes_for_diff_graph(&vcf, &decoy_samples, mark_shorter_alleles);

    tracing::debug!("Finished sorting");

    let shared_ranges = find_shared_haplotype_ranges(&vcf);

    print_ranges_to_csv(

        &vcf,

        &shared_ranges,

        &decoy_samples,

        &only_longest,

        &mut writer,

    )?;

    tracing::debug!(

    if want_png {

        let mut dg = MatrixGraph::new(&vcf, graph_args, decoy_samples, mark_shorter_alleles);

        dg.draw_graph(&index_order);

        svg::save(img_output, &dg.document)?;

    tracing::debug!("Finished drawing graph");

    Ok(())

}

fn sort_indexes_for_diff_graph(

    vcf: &PhasedMatrix,

    decoy_samples: &Option<Vec<String>>,

    mark_shorter_alleles: bool,

) -> Vec<usize> {

    // Take all from start to variant index, reverse and calculate 1 count in parallel

    // to get the amount of markers shared to the left

    let mut values: Vec<(usize, i32)> = vcf

        .matrix

        .slice(s![.., 0..vcf.variant_idx()])

        .axis_iter(Axis(0))

        .into_par_iter()

        .enumerate()

        .map(|(y, row)| {

            let mut count = 0;

            for i in row.iter().rev() {

                match i {

                    0 => break,

                    1 => count += 1,

            (y, count)

        })

        .collect();

    // Sort shortest alleles to the top

    if mark_shorter_alleles {

        let shared_lengths = shared_lengths_by_majority(vcf, vcf.variant_idx());

        let only_longest = select_only_longest_alleles(&shared_lengths);

        let only_longest: Vec<usize> = only_longest.iter().map(|s| s.index).collect();

        values.sort_by(|a, b| {

            only_longest

                .contains(&a.0)

                .cmp(&only_longest.contains(&b.0))

        });

    }

    if let Some(samples) = decoy_samples {

    }

    values.sort_by(|a, b| a.1.cmp(&b.1));

    values.iter().map(|v| v.0).collect::<Vec<usize>>()

}

pub fn transform_gt_matrix_to_match_matrix(

    mut vcf: PhasedMatrix,

    ht: &Vec<HapVariant>,

) -> PhasedMatrix {

    // Info to users

    match ht.len().cmp(&vcf.matrix().shape()[1]) {

        std::cmp::Ordering::Equal => (),

    let coords = std::mem::take(vcf.coords_mut());

    vcf.matrix

        .axis_iter_mut(Axis(0))

        .into_par_iter()

        .for_each_with(&coords, |coords, mut row| {

            for (i, gt) in row.iter_mut().enumerate() {

                let coord = coords.get(i).unwrap();

                if let Some(variant_idx) = ht.iter().position(|c| c == coord) {

                    *gt = match &ht[variant_idx].gt == gt {

                        true => 1,

                        false => 0,

        });

    // Swap back

    vcf.set_coords(coords);

    vcf

}

pub fn find_shared_haplotype_ranges(vcf: &PhasedMatrix) -> Vec<(usize, usize, usize)> {

    vcf.matrix

        .axis_iter(Axis(0))

        .into_par_iter()

        .enumerate()

        .map(|(idx, row)| {

            let (mut last_start, mut last_end) = (0, 0);

            for (i, gt) in row.iter().enumerate() {

                match gt {

                    1 => last_end = i,

                        if vcf.variant_idx() <= last_end && vcf.variant_idx() >= last_start {

                            return (idx, last_start, last_end);

                        }

                        last_start = i;

                        last_end = i;

            (idx, 0, vcf.ncoords() - 1)

        })

        .collect()

}

    writer.write_record(vec![

        "id",

        "start",

        "stop",

        "length",

        "markers",

        "is_decoy",

        "is_longest",

    ])?;

    for (idx, start, stop) in ranges {

        let start_pos = vcf.get_pos(*start);

        let stop_pos = vcf.get_pos(*stop);

        let allele = if *idx >= vcf.nsamples() / 2 { 2 } else { 1 };

        let mut row = vec![

                format!("{}_{allele}",vcf.get_sample_name(*idx)),

                start_pos.to_string(),

                stop_pos.to_string(),

                (stop_pos - start_pos).to_string(),

                (stop - start).to_string(),

            ];

        match decoy_samples {

            None => row.push("na".to_string()),

        match only_longest {

            Some(longest) => match longest.contains(idx) {

                false => row.push("false".to_string()),

            None => row.push("na".to_string()),

        writer.write_record(row)?;

    Ok(())

}