debruijn/

kmer.rs

// Copyright 2017 10x Genomics

//! Represent kmers with statically know length in compact integer types
//!
//! A kmer is a DNA sequence with statically-known length K. The sequence is stored the smallest possible integer type
//! Efficient methods for reverse complement and shifting new bases into the left or right of a sequence are provided.
//! Kmers implement `Eq` to test if two kmers represent the same string.
//! Kmers implement `Ord`, which corresponds to the natural lexicographic ordering.
//!
//! ```
//! use debruijn::*;
//! use debruijn::kmer::*;
//!
//! let k1 = Kmer16::from_ascii(b"ACGTACGTACGTACGT");
//!
//! // Reverse complement
//! let rc_k1 = k1.rc();
//!
//! // Double reverse complement
//! let k1_copy = rc_k1.rc();
//! assert_eq!(k1, k1_copy);
//!
//! // Push one base onto the left
//! assert_eq!(k1.extend_left(base_to_bits(b'T')), Kmer16::from_ascii(b"TACGTACGTACGTACG"));
//!
//! // Generate a set of kmers from a string, then sort
//! let mut all_kmers = Kmer16::kmers_from_ascii(b"TACGTACGTACGTACGTT");
//! all_kmers.sort();
//! assert_eq!(all_kmers,
//!     vec![
//!         Kmer16::from_ascii(b"ACGTACGTACGTACGT"),
//!         Kmer16::from_ascii(b"CGTACGTACGTACGTT"),
//!         Kmer16::from_ascii(b"TACGTACGTACGTACG")
//!     ]);

use kmersize_derive::KmerSize;
use num_traits::FromPrimitive;
use num_traits::PrimInt;
use serde::Deserialize;
use serde::Serialize;
use serde::de::DeserializeOwned;
use std;
use std::fmt;
use std::fmt::Binary;
use std::fmt::Debug;
use std::fmt::Display;
use std::hash::Hash;
use std::marker::PhantomData;

use crate::bits_to_base;
use crate::Kmer;
use crate::Mer;

use serde_big_array::BigArray;

// Pre-defined kmer types


/// 128-base kmer, backed by four u64s
pub type Kmer128 = VarLenKmer<u64, 4, K128>;
/// 112-base kmer, backed by seven u32s
pub type Kmer112 = VarLenKmer<u32, 7, K112>;
/// 96-base kmer, backed by three u64s
pub type Kmer96 = VarLenKmer<u64, 3, K96>;
/// 80-base kmer, backed by five u32s
pub type Kmer80 = VarLenKmer<u32, 5, K80>;
/// 64-base kmer, backed by a single u128
pub type Kmer64 = IntKmer<u128>;
/// 63-base kmer, backed by a single u128
pub type Kmer63 = VarIntKmer<u128, K63>;
/// 62-base kmer, backed by a single u128
pub type Kmer62 = VarIntKmer<u128, K62>;
/// 61-base kmer, backed by a single u128
pub type Kmer61 = VarIntKmer<u128, K61>;
/// 60-base kmer, backed by a single u128
pub type Kmer60 = VarIntKmer<u128, K60>;
/// 59-base kmer, backed by a single u128
pub type Kmer59 = VarIntKmer<u128, K59>;
/// 58-base kmer, backed by a single u128
pub type Kmer58 = VarIntKmer<u128, K58>;
/// 57-base kmer, backed by a single u128
pub type Kmer57 = VarIntKmer<u128, K57>;
/// 56-base kmer, backed by a single u128
pub type Kmer56 = VarIntKmer<u128, K56>;
/// 55-base kmer, backed by a single u128
pub type Kmer55 = VarIntKmer<u128, K55>;
/// 54-base kmer, backed by a single u128
pub type Kmer54 = VarIntKmer<u128, K54>;
/// 53-base kmer, backed by a single u128
pub type Kmer53 = VarIntKmer<u128, K53>;
/// 52-base kmer, backed by a single u128
pub type Kmer52 = VarIntKmer<u128, K52>;
/// 51-base kmer, backed by a single u128
pub type Kmer51 = VarIntKmer<u128, K51>;
/// 50-base kmer, backed by a single u128
pub type Kmer50 = VarIntKmer<u128, K50>;
/// 49-base kmer, backed by a single u128
pub type Kmer49 = VarIntKmer<u128, K49>;
/// 48-base kmer, backed by a single u128
pub type Kmer48 = VarIntKmer<u128, K48>;
/// 47-base kmer, backed by a single u128
pub type Kmer47 = VarIntKmer<u128, K47>;
/// 46-base kmer, backed by a single u128
pub type Kmer46 = VarIntKmer<u128, K46>;
/// 45-base kmer, backed by a single u128
pub type Kmer45 = VarIntKmer<u128, K45>;
/// 44-base kmer, backed by a single u128
pub type Kmer44 = VarIntKmer<u128, K44>;
/// 43-base kmer, backed by a single u128
pub type Kmer43 = VarIntKmer<u128, K43>;
/// 42-base kmer, backed by a single u128
pub type Kmer42 = VarIntKmer<u128, K42>;
/// 41-base kmer, backed by a single u128
pub type Kmer41 = VarIntKmer<u128, K41>;
/// 40-base kmer, backed by a single u128
pub type Kmer40 = VarIntKmer<u128, K40>;
/// 39-base kmer, backed by a single u128
pub type Kmer39 = VarIntKmer<u128, K39>;
/// 38-base kmer, backed by a single u128
pub type Kmer38 = VarIntKmer<u128, K38>;
/// 37-base kmer, backed by a single u128
pub type Kmer37 = VarIntKmer<u128, K37>;
/// 36-base kmer, backed by a single u128
pub type Kmer36 = VarIntKmer<u128, K36>;
/// 35-base kmer, backed by a single u128
pub type Kmer35 = VarIntKmer<u128, K35>;
/// 34-base kmer, backed by a single u128
pub type Kmer34 = VarIntKmer<u128, K34>;
/// 33-base kmer, backed by a single u128
pub type Kmer33 = VarIntKmer<u128, K33>;
/// 32-base kmer, backed by a single u64
pub type Kmer32 = IntKmer<u64>;
/// 31-base kmer, backed by a single u64
pub type Kmer31 = VarIntKmer<u64, K31>;
/// 30-base kmer, backed by a single u64
pub type Kmer30 = VarIntKmer<u64, K30>;
/// 29-base kmer, backed by a single u64
pub type Kmer29 = VarIntKmer<u64, K29>;
/// 28-base kmer, backed by a single u64
pub type Kmer28 = VarIntKmer<u64, K28>;
/// 27-base kmer, backed by a single u64
pub type Kmer27 = VarIntKmer<u64, K27>;
/// 25-base kmer, backed by a single u64
pub type Kmer26 = VarIntKmer<u64, K26>;
/// 25-base kmer, backed by a single u64
pub type Kmer25 = VarIntKmer<u64, K25>;
/// 24-base kmer, backed by a single u64
pub type Kmer24 = VarIntKmer<u64, K24>;
/// 23-base kmer, backed by a single u64
pub type Kmer23 = VarIntKmer<u64, K23>;
/// 22-base kmer, backed by a single u64
pub type Kmer22 = VarIntKmer<u64, K22>;
/// 21-base kmer, backed by a single u64
pub type Kmer21 = VarIntKmer<u64, K21>;
/// 20-base kmer, backed by a single u64
pub type Kmer20 = VarIntKmer<u64, K20>;
/// 19-base kmer, backed by a single u64
pub type Kmer19 = VarIntKmer<u64, K19>;
/// 18-base kmer, backed by a single u64
pub type Kmer18 = VarIntKmer<u64, K18>;
/// 17-base kmer, backed by a single u64
pub type Kmer17 = VarIntKmer<u64, K17>;
/// 16-base kmer, backed by a single u32
pub type Kmer16 = IntKmer<u32>;
/// 15-base kmer, backed by a single u32
pub type Kmer15 = VarIntKmer<u32, K15>;
/// 14-base kmer, backed by a single u32
pub type Kmer14 = VarIntKmer<u32, K14>;
/// 13-base kmer, backed by a single u32
pub type Kmer13 = VarIntKmer<u32, K13>;
/// 12-base kmer, backed by a single u32
pub type Kmer12 = VarIntKmer<u32, K12>;
/// 11-base kmer, backed by a single u32
pub type Kmer11 = VarIntKmer<u32, K11>;
/// 10-base kmer, backed by a single u32
pub type Kmer10 = VarIntKmer<u32, K10>;
/// 9-base kmer, backed by a single u32
pub type Kmer9 = VarIntKmer<u32, K9>;
/// 8-base kmer, backed by a single u16
pub type Kmer8 = IntKmer<u16>;
/// 7-base kmer, backed by a single u16
pub type Kmer7 = VarIntKmer<u16, K7>;
/// 6-base kmer, backed by a single u16
pub type Kmer6 = VarIntKmer<u16, K6>;
/// 5-base kmer, backed by a single u16
pub type Kmer5 = VarIntKmer<u16, K5>;
/// 4-base kmer, backed by a single u8
pub type Kmer4 = IntKmer<u8>;
/// 3-base kmer, backed by a single u8
pub type Kmer3 = VarIntKmer<u8, K3>;
/// 2-base kmer, backed by a single u8
pub type Kmer2 = VarIntKmer<u8, K2>;



/// Trait for specialized integer operations used in DeBruijn Graph
pub trait IntHelp: PrimInt + FromPrimitive + Hash + Serialize + Display + Debug + Binary {
    /// Reverse the order of 2-bit units of the integer
    fn reverse_by_twos(&self) -> Self;

    fn lower_of_two() -> Self;
}

impl IntHelp for u128 {
    #[inline]
    fn reverse_by_twos(&self) -> u128 {
        // swap adjacent pairs
        let mut r = ((self & 0x33333333333333333333333333333333u128) << 2)
            | ((self >> 2) & 0x33333333333333333333333333333333u128);

        // swap nibbles
        r = ((r & 0x0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F0Fu128) << 4)
            | ((r >> 4) & 0x0F0F0F0F0F0F0F0F0F0F0F0F0F0F0F0Fu128);

        // swap bytes
        r = ((r & 0x00FF00FF00FF00FF00FF00FF00FF00FFu128) << 8)
            | ((r >> 8) & 0x00FF00FF00FF00FF00FF00FF00FF00FFu128);

        // swap 2 bytes
        r = ((r & 0x0000FFFF0000FFFF0000FFFF0000FFFFu128) << 16)
            | ((r >> 16) & 0x0000FFFF0000FFFF0000FFFF0000FFFFu128);

        // swap 4 bytes
        r = ((r & 0x00000000FFFFFFFF00000000FFFFFFFFu128) << 32)
            | ((r >> 32) & 0x00000000FFFFFFFF00000000FFFFFFFFu128);

        // swap 8 bytes
        r = ((r & 0x0000000000000000FFFFFFFFFFFFFFFFu128) << 64)
            | ((r >> 64) & 0x0000000000000000FFFFFFFFFFFFFFFFu128);

        r
    }

    #[inline]
    fn lower_of_two() -> u128 {
        0x55555555555555555555555555555555u128
    }
}

impl IntHelp for u64 {
    #[inline]
    fn reverse_by_twos(&self) -> u64 {
        // swap adjacent pairs
        let mut r = ((self & 0x3333333333333333u64) << 2) | ((self >> 2) & 0x3333333333333333u64);

        // swap nibbles
        r = ((r & 0x0F0F0F0F0F0F0F0Fu64) << 4) | ((r >> 4) & 0x0F0F0F0F0F0F0F0Fu64);

        // swap bytes
        r = ((r & 0x00FF00FF00FF00FFu64) << 8) | ((r >> 8) & 0x00FF00FF00FF00FFu64);

        // swap 2 bytes
        r = ((r & 0x0000FFFF0000FFFFu64) << 16) | ((r >> 16) & 0x0000FFFF0000FFFFu64);

        // swap 4 bytes
        r = ((r & 0x00000000FFFFFFFFu64) << 32) | ((r >> 32) & 0x00000000FFFFFFFFu64);

        r
    }

    #[inline]
    fn lower_of_two() -> u64 {
        0x5555555555555555u64
    }
}

impl IntHelp for u32 {
    #[inline]
    fn reverse_by_twos(&self) -> u32 {
        // swap adjacent pairs
        let mut r = ((self & 0x33333333u32) << 2) | ((self >> 2) & 0x33333333u32);

        // swap nibbles
        r = ((r & 0x0F0F0F0Fu32) << 4) | ((r >> 4) & 0x0F0F0F0Fu32);

        // swap bytes
        r = ((r & 0x00FF00FFu32) << 8) | ((r >> 8) & 0x00FF00FFu32);

        // swap 2 bytes
        r = ((r & 0x0000FFFFu32) << 16) | ((r >> 16) & 0x0000FFFFu32);

        r
    }

    #[inline]
    fn lower_of_two() -> u32 {
        0x55555555u32
    }
}

impl IntHelp for u16 {
    #[inline]
    fn reverse_by_twos(&self) -> u16 {
        // swap adjacent pairs
        let mut r = ((self & 0x3333u16) << 2) | ((self >> 2) & 0x3333u16);

        // swap nibbles
        r = ((r & 0x0F0Fu16) << 4) | ((r >> 4) & 0x0F0Fu16);

        // swap bytes
        r = ((r & 0x00FFu16) << 8) | ((r >> 8) & 0x00FFu16);

        r
    }

    #[inline]
    fn lower_of_two() -> u16 {
        0x5555u16
    }
}

impl IntHelp for u8 {
    #[inline]
    fn reverse_by_twos(&self) -> u8 {
        // swap adjacent pairs
        let mut r = ((self & 0x33u8) << 2) | ((self >> 2) & 0x33u8);

        // swap nibbles
        r = ((r & 0x0Fu8) << 4) | ((r >> 4) & 0x0Fu8);

        r
    }

    #[inline]
    fn lower_of_two() -> u8 {
        0x55u8
    }
}

/// A Kmer sequence with a statically know K. K will fill the underlying integer type.
#[derive(Copy, Clone, PartialEq, PartialOrd, Eq, Ord, Hash, Serialize, Deserialize)]
pub struct IntKmer<T: IntHelp> {
    pub storage: T,
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp + Sized> IntKmer<T> {
    fn msk() -> T {
        T::one() << 1 | T::one()
    }

    #[inline(always)]
    fn to_byte(v: T) -> u8 {
        T::to_u8(&v).unwrap()
    }

    #[inline(always)]
    fn t_from_byte(v: u8) -> T {
        T::from_u8(v).unwrap()
    }

    #[inline(always)]
    fn t_from_u64(v: u64) -> T {
        T::from_u64(v).unwrap()
    }

    #[inline]
    fn addr(&self, pos: usize) -> usize {
        let top_base = Self::k() - 1;
        (top_base - pos) * 2
    }

    #[inline(always)]
    fn _k() -> usize {
        // 4 bases per byte
        std::mem::size_of::<T>() * 4
    }

    #[inline(always)]
    fn _bits() -> usize {
        std::mem::size_of::<T>() * 8
    }

    #[inline(always)]
    pub fn top_mask(n_bases: usize) -> T {
        if n_bases > 0 {
            // first pos bases
            let one = T::one();
            ((one << (n_bases * 2)) - one) << (Self::_bits() - n_bases * 2)
        } else {
            T::zero()
        }
    }

    #[inline(always)]
    pub fn bottom_mask(n_bases: usize) -> T {
        if n_bases > 0 {
            // first pos bases
            let one = T::one();
            (one << (n_bases * 2)) - one
        } else {
            T::zero()
        }
    }
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp> Mer for IntKmer<T> {
    #[inline(always)]
    fn len(&self) -> usize {
        Self::_k()
    }

    #[inline(always)]
    fn is_empty(&self) -> bool {
        false
    }

    /// Get the letter at the given position.
    #[inline(always)]
    fn get(&self, pos: usize) -> u8 {
        let bit = self.addr(pos);
        Self::to_byte(self.storage >> bit & Self::msk())
    }

    fn set_mut(&mut self, pos: usize, v: u8) {
        let bit = self.addr(pos);
        let mask = !(Self::msk() << bit);

        self.storage = (self.storage & mask) | (Self::t_from_byte(v) << bit);
    }

    /// Set a slice of bases in the kmer, using the packed representation in value.
    /// Sets n_bases, starting at pos. Bases must always be packed into the upper-most
    /// bits of the value.
    #[inline(always)]
    fn set_slice_mut(&mut self, pos: usize, n_bases: usize, value: u64) {
        debug_assert!(pos + n_bases <= Self::k());

        let v_shift = if Self::_bits() < 64 {
            value >> (64 - Self::_bits())
        } else {
            value
        };

        let v = if Self::_bits() > 64 {
            Self::t_from_u64(v_shift) << (Self::_bits() - 64)
        } else {
            Self::t_from_u64(v_shift)
        };

        let top_mask = Self::top_mask(pos);
        let bottom_mask = Self::bottom_mask(Self::k() - (pos + n_bases));
        let mask = top_mask | bottom_mask;

        let value_slide = v >> (2 * pos);

        self.storage = (self.storage & mask) | (value_slide & !mask);
    }

    /// Return the reverse complement of this kmer
    fn rc(&self) -> Self {
        // not bits to get complement, then reverse order
        let new = !self.storage.reverse_by_twos();

        // NOTE: IntKmer always fills the bits, so we don't need to shift here.
        IntKmer { storage: new }
    }

    fn at_count(&self) -> u32 {
        // A's and T's have upper_bit ^ lower_bit == 0
        // count how many of these are present
        let mix_base_bits = !((self.storage >> 1) ^ self.storage);
        let mask_lower = mix_base_bits & IntHelp::lower_of_two();
        mask_lower.count_ones()
    }

    fn gc_count(&self) -> u32 {
        // A's and T's have upper_bit ^ lower_bit == 1
        // count how many of these are present
        let mix_base_bits = (self.storage >> 1) ^ self.storage;
        let mask_lower = mix_base_bits & IntHelp::lower_of_two();
        mask_lower.count_ones()
    }
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp> Kmer for IntKmer<T> {
    fn empty() -> Self {
        IntKmer { storage: T::zero() }
    }

    #[inline(always)]
    fn k() -> usize {
        Self::_k()
    }

    #[inline(always)]
    fn from_u64(v: u64) -> IntKmer<T> {
        IntKmer {
            storage: Self::t_from_u64(v),
        }
    }

    fn to_u64(&self) -> u64 {
        T::to_u64(&self.storage).unwrap()
    }

    /// Shift the base v into the left end of the kmer
    fn extend_left(&self, v: u8) -> Self {
        let new = self.storage >> 2 | (Self::t_from_byte(v) << ((Self::k() - 1) * 2));
        IntKmer { storage: new }
    }

    fn extend_right(&self, v: u8) -> Self {
        let new = self.storage << 2;
        let mut kmer = IntKmer { storage: new };
        kmer.set_mut(Self::k() - 1, v);
        kmer
    }

    /// Hamming distance between this kmer and another kmer
    fn hamming_dist(&self, other: Self) -> u32 {
        let bit_diffs = self.storage ^ other.storage;
        let two_bit_diffs = (bit_diffs | bit_diffs >> 1) & IntHelp::lower_of_two();
        two_bit_diffs.count_ones()
    }
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp> fmt::Debug for IntKmer<T> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut s = String::new();
        for pos in 0..Self::k() {
            s.push(bits_to_base(self.get(pos)))
        }

        write!(f, "{}", s)
    }
}

/// Helper trait for declaring the K value of a Kmer. Will be removed when const generics are available
pub trait KmerSize: Ord + Hash + Copy + fmt::Debug {
    #[allow(non_snake_case)]
    fn K() -> usize;
}

/// A fixed-length Kmer sequence that may not fill the bits of T
///
/// side:             L           R
/// bases: 0   0   0  A  C  G  T  T
/// bits:             H  ........ L
/// bit :  14  12  10 8  6  4  2  0
///
/// sorting the integer will give a lexicographic sorting of the corresponding string.
///  kmers that don't fill `storage` are always aligned to the least significant bits
#[derive(Copy, Clone, PartialEq, PartialOrd, Eq, Ord, Hash, Serialize, Deserialize)]
pub struct VarIntKmer<T: IntHelp, KS: KmerSize> {
    pub storage: T,
    pub phantom: PhantomData<KS>,
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp, KS: KmerSize> Kmer for VarIntKmer<T, KS> {
    fn empty() -> Self {
        VarIntKmer {
            storage: T::zero(),
            phantom: PhantomData,
        }
    }

    #[inline]
    fn k() -> usize {
        Self::_k()
    }

    fn to_u64(&self) -> u64 {
        T::to_u64(&self.storage).unwrap()
    }

    fn from_u64(v: u64) -> Self {
        VarIntKmer {
            storage: Self::t_from_u64(v),
            phantom: PhantomData,
        }
    }

    /// Shift the base v into the left end of the kmer
    fn extend_left(&self, v: u8) -> Self {
        let new = self.storage >> 2;
        let mut kmer = VarIntKmer {
            storage: new,
            phantom: PhantomData,
        };
        kmer.set_mut(0, v);
        kmer
    }

    fn extend_right(&self, v: u8) -> Self {
        let new = self.storage << 2 & !Self::top_mask(0);
        let mut kmer = VarIntKmer {
            storage: new,
            phantom: PhantomData,
        };
        kmer.set_mut(Self::k() - 1, v);
        kmer
    }

    fn hamming_dist(&self, other: Self) -> u32 {
        let bit_diffs = self.storage ^ other.storage;
        let two_bit_diffs = (bit_diffs | bit_diffs >> 1) & IntHelp::lower_of_two();
        two_bit_diffs.count_ones()
    }
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp, KS: KmerSize> VarIntKmer<T, KS> {
    #[inline(always)]
    fn msk() -> T {
        T::one() << 1 | T::one()
    }

    fn to_byte(v: T) -> u8 {
        T::to_u8(&v).unwrap()
    }

    fn t_from_byte(v: u8) -> T {
        T::from_u8(v).unwrap()
    }

    fn t_from_u64(v: u64) -> T {
        T::from_u64(v).unwrap()
    }

    #[inline(always)]
    fn addr(&self, pos: usize) -> usize {
        let top_base = Self::k() - 1;
        (top_base - pos) * 2
    }

    // K of this kmer
    #[inline(always)]
    fn _k() -> usize {
        KS::K()
    }

    // Bits used by this kmer
    #[inline(always)]
    fn _bits() -> usize {
        Self::_k() * 2
    }

    #[inline(always)]
    fn _total_bits() -> usize {
        std::mem::size_of::<T>() * 8
    }

    // mask the unused bits at the top, plus the requested number of bases
    #[inline(always)]
    pub fn top_mask(n_bases: usize) -> T {
        let unused_bits = Self::_total_bits() - Self::_bits();
        let mask_bits = n_bases * 2 + unused_bits;

        if mask_bits > 0 {
            let one = T::one();
            ((one << mask_bits) - one) << (Self::_total_bits() - mask_bits)
        } else {
            T::zero()
        }
    }

    #[inline(always)]
    pub fn bottom_mask(n_bases: usize) -> T {
        if n_bases > 0 {
            // first pos bases
            let one = T::one();
            (one << (n_bases * 2)) - one
        } else {
            T::zero()
        }
    }
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp, KS: KmerSize> Mer for VarIntKmer<T, KS> {
    #[inline(always)]
    fn len(&self) -> usize {
        Self::_k()
    }

    fn is_empty(&self) -> bool {
        Self::_k() == 0
    }

    /// Get the letter at the given position.
    fn get(&self, pos: usize) -> u8 {
        let bit = self.addr(pos);
        Self::to_byte(self.storage >> bit & Self::msk())
    }

    fn set_mut(&mut self, pos: usize, v: u8) {
        let bit = self.addr(pos);
        let mask = !(Self::msk() << bit);

        self.storage = (self.storage & mask) | (Self::t_from_byte(v) << bit);
    }

    /// Set a slice of bases in the kmer, using the packed representation in value.
    /// Sets n_bases, starting at pos. Incoming bases must always be packed into the upper-most
    /// bits of the value.
    #[inline(always)]
    fn set_slice_mut(&mut self, pos: usize, n_bases: usize, value: u64) {
        debug_assert!(pos + n_bases <= Self::k());

        // Move bases to the top of the smaller type
        let v_shift = if Self::_total_bits() < 64 {
            value >> (64 - Self::_total_bits())
        } else {
            value
        };

        // Move bases up to the top of this type
        let v = if Self::_total_bits() > 64 {
            Self::t_from_u64(v_shift) << (Self::_total_bits() - 64)
        } else {
            Self::t_from_u64(v_shift)
        };

        // Mask for where the bases will sit in the kmer
        let top_mask = Self::top_mask(pos);
        let bottom_mask = Self::bottom_mask(Self::k() - (pos + n_bases));
        let mask = top_mask | bottom_mask;

        // Move the base down to their home: position + unused high-order bits
        let shift = (2 * pos) + (Self::_total_bits() - Self::_bits());
        let value_slide = v >> shift;

        self.storage = (self.storage & mask) | (value_slide & !mask);
    }

    /// Return the reverse complement of this kmer
    fn rc(&self) -> Self {
        // not bits to get complement, then reverse order
        let mut new = !self.storage.reverse_by_twos();

        // deal with case when the kmer doesn't fill the bits
        if Self::k() < std::mem::size_of::<T>() * 4 {
            let up_shift = 2 * (std::mem::size_of::<T>() * 4 - Self::k());
            new = new >> up_shift;
        }

        VarIntKmer {
            storage: new,
            phantom: PhantomData,
        }
    }

    fn at_count(&self) -> u32 {
        // A's and T's have upper_bit ^ lower_bit == 0
        // count how many of these are present
        let mix_base_bits = !((self.storage >> 1) ^ self.storage);
        let mask_lower = mix_base_bits & !Self::top_mask(0) & IntHelp::lower_of_two();
        mask_lower.count_ones()
    }

    fn gc_count(&self) -> u32 {
        // A's and T's have upper_bit ^ lower_bit == 1
        // count how many of these are present
        let mix_base_bits = (self.storage >> 1) ^ self.storage;
        let mask_lower = mix_base_bits & !Self::top_mask(0) & IntHelp::lower_of_two();
        mask_lower.count_ones()
    }
}

impl<T: PrimInt + FromPrimitive + Hash + IntHelp, KS: KmerSize> fmt::Debug for VarIntKmer<T, KS> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut s = String::new();
        for pos in 0..Self::k() {
            s.push(bits_to_base(self.get(pos)))
        }

        write!(f, "{}", s)
    }
}

/// A fixed-length Kmer sequence that may not fill the bits of T
///
/// side:             L           R
/// bases: 0   0   0  A  C  G  T  T
/// bits:             H  ........ L
/// bit :  14  12  10 8  6  4  2  0
///
/// sorting the integer will give a lexicographic sorting of the corresponding string.
///  kmers that don't fill `storage` are always aligned to the least significant bits
#[derive(Copy, Clone, PartialEq, PartialOrd, Eq, Ord, Hash, Serialize, Deserialize)]
pub struct VarLenKmer<T: IntHelp + DeserializeOwned, const LEN: usize, KS: KmerSize> 
{
    #[serde(with="BigArray")]
    pub storage: [T; LEN],
    pub phantom: PhantomData<KS>,
}

impl<T: IntHelp + DeserializeOwned, const LEN: usize, KS: KmerSize> Kmer for VarLenKmer<T, LEN, KS> {
    fn empty() -> Self {
        VarLenKmer {
            storage: [T::zero(); LEN],
            phantom: PhantomData,
        }
    }

    #[inline]
    fn k() -> usize {
        Self::_k()
    }

    fn to_u64(&self) -> u64 {
        unimplemented!() // cant really turn the array into an u64
    }

    fn from_u64(_v: u64) -> Self {
        unimplemented!() // cant really turn the u64 into an array
    }

    /// Shift the base v into the left end of the kmer
    fn extend_left(&self, v: u8) -> Self {
        let mut new_kmer = VarLenKmer::empty();
        // for the back storage block, shift right by two
        let new_back_block = self.storage[LEN-1] >> 2;
        new_kmer.storage[LEN-1] = new_back_block;

        // go over the rest of the blocks
        // should we only have one block, this is skipped
        for step in (0..(LEN-1)).rev() {
            // always get the back base of block and attach it to the front of the previous one
            let moved_v = self.get_by_addr(step, 0);
            new_kmer.set_by_addr(step+1, Self::t_bits()-2, moved_v);

            // then, shift the block by 2 and add it to the new_kmer
            let new_block = self.storage[step] >> 2;
            new_kmer.storage[step] = new_block;
        }

        // finally, add the new base where the k-mer ends in the front block
        new_kmer.set_mut(0, v);
        new_kmer
    }

    /// shift the base v into the right end of the k-mer
    fn extend_right(&self, v: u8) -> Self {
        let mut new_kmer = VarLenKmer::empty();
        // for the front storage block, shift left by two with the mask
        let new_front_block = self.storage[0] << 2 & !Self::unused_bits_mask(0);
        new_kmer.storage[0] = new_front_block;

        // go over the rest of the blocks
        // should we only have one block, this is skipped
        for step in 1..(LEN) {
            // always get the front base of block and attach it to the back of the previous one
            let moved_v = self.get_by_addr(step, Self::t_bits()-2);
            new_kmer.set_by_addr(step-1, 0, moved_v);

            // then, shift the block by 2 and add it to the new_kmer
            let new_block = self.storage[step] << 2;
            new_kmer.storage[step] = new_block;
        }

        // finally, add the new base to the back block
        new_kmer.set_mut(Self::k()-1, v);
        new_kmer
    }

    fn hamming_dist(&self, other: Self) -> u32 {
        let mut dist = 0;
        for block in 0..LEN {
            let bit_diffs = self.storage[block] ^ other.storage[block];
            let two_bit_diffs = (bit_diffs | bit_diffs >> 1) & IntHelp::lower_of_two();
            dist += two_bit_diffs.count_ones();
        }

        dist
    }
}

impl<T: IntHelp + DeserializeOwned, const LEN: usize, KS: KmerSize> VarLenKmer<T, LEN, KS> {
    fn to_byte(v: T) -> u8 {
        T::to_u8(&v).unwrap()
    }

    fn t_from_byte(v: u8) -> T {
        T::from_u8(v).unwrap()
    }

    fn t_from_u64(v: u64) -> T {
        T::from_u64(v).unwrap()
    }

    /// get the (block, bit) for the position in the k-mer
    /// 
    /// bits start counting from the back
    #[inline(always)]
    fn addr(&self, pos: usize) -> (usize, usize) {
        let overall_i = (Self::_k() - 1 - pos) * 2;
        let block = LEN - 1 - (overall_i / Self::t_bits());
        let bit = overall_i % Self::t_bits();
        (block, bit)
    }

    #[inline(always)]
    fn get_by_addr(&self, block: usize, bit: usize) -> u8 {
        let mask = Self::t_from_byte(3);
        Self::to_byte((self.storage[block] >> bit) & mask)
    }

    #[inline(always)]
    fn set_by_addr(&mut self, block: usize, bit: usize, value: u8) {
        let mask = 3;
        // check that value has only one base, cast to T, and mov to correct position
        let value = Self::t_from_byte(value & mask) << bit;
        // remove the old bits from the block
        let block_mask = Self::t_from_byte(3) << bit;
        let new_block = self.storage[block] & !block_mask;
        // add new block with new base into storage
        self.storage[block] = new_block | value;
    }

    /// K of this kmer
    #[inline(always)]
    fn _k() -> usize {
        KS::K()
    }

    /// Bits used by this kmer
    #[inline(always)]
    fn _bits() -> usize {
        Self::_k() * 2
    }

    /// bits in the data type
    #[inline(always)]
    fn _total_bits() -> usize {
        std::mem::size_of::<T>() * 8 * LEN
    }

    /// number of bits in the type T
    #[inline(always)]
    fn t_bits() -> usize {
        std::mem::size_of::<T>() * 8
    }

    // TODO check
    /// mask the unused bits at the top, plus the requested number of bases
    #[inline(always)]
    pub fn unused_bits_mask(n_bases: usize) -> T {
        let unused_bits = Self::_total_bits() - Self::_bits();

        assert!(unused_bits < Self::t_bits(), "please use a smaller k-mer impl for this k");

        let mask_bits = n_bases * 2 + unused_bits;

        if mask_bits > 0 {
            let one = T::one();
            ((one << mask_bits) - one) << (Self::t_bits() - mask_bits)
        } else {
            T::zero()
        }
    }
}

impl<T: IntHelp + DeserializeOwned, const LEN: usize, KS: KmerSize> Mer for VarLenKmer<T, LEN, KS> {
    #[inline(always)]
    fn len(&self) -> usize {
        Self::_k()
    }

    fn is_empty(&self) -> bool {
        Self::_k() == 0
    }

    #[inline(always)]
    /// Get the letter at the given position.
    fn get(&self, pos: usize) -> u8 {
        let (block, bit) = self.addr(pos);
        self.get_by_addr(block, bit)
    }

    #[inline(always)]
    fn set_mut(&mut self, pos: usize, v: u8) {
        let (block, bit) = self.addr(pos);
        self.set_by_addr(block, bit, v);
    }

    /// Set a slice of bases in the kmer, using the packed representation in value.
    /// Sets n_bases, starting at pos. Incoming bases must always be packed into the front
    /// bits of the value.
    #[inline(always)]
    fn set_slice_mut(&mut self, pos: usize, n_bases: usize, value: u64) {
        /* let one = T::one();
        let value_bits = n_bases * 2;
        // get starting address
        let (start_block, start_bit) = self.addr(pos);
        let start_bits = start_bit + 2;
        let rest_bits = value_bits.saturating_sub(start_bits);
        let middle_blocks = rest_bits / Self::t_bits();
        let back_bits = rest_bits.checked_rem(middle_blocks).unwrap_or(rest_bits);

        // if the slice fits into the first block, with some room in the back -> mask the back
        let gap = start_bits.saturating_sub(value_bits);
        let gap_mask = if value_bits < start_bit { (one << gap) - one } else { T::zero() };
        // move first subslice to back, cast to T, then move back up by the gap
        let first_subslice = Self::t_from_u64(value >> (64 - (start_bits - gap))) << gap;

        // create mask to clear old bits in block, ignore potential back bits
        let first_mask = if start_bits == Self::t_bits() {
            !gap_mask
        } else {
            ((one << start_bits) - one) & !gap_mask
        };
        // clear old bits from block, and add starting slice
        self.storage[start_block] = (self.storage[start_block] & !first_mask) | first_subslice;

        if middle_blocks == 0 && back_bits == 0 { return; } // we are done

        // for each full block, extract subslice into T and replace block (in case of T u128 there are no middle blocks)
        let end_middle_block = start_block + middle_blocks;
        let t_mask = u64::MAX.checked_shl(Self::t_bits() as u32).unwrap_or(0);

        for (i, full_block) in ((start_block+1)..(end_middle_block+1)).enumerate() {
            // move relevant bits to back of value  and apply mask
            let subslice = (value >> (64 - start_bits - Self::t_bits() * i)) & !t_mask;
            // cast into T and replace block
            self.storage[full_block] = Self::t_from_u64(subslice);
        }

        // insert back bits
        if back_bits > 0 {
            let back_gap_mask = (one << (Self::t_bits() - back_bits)) - one;
            // shift last bits to end of value, and then back to the start of T
            let subslice = Self::t_from_u64(value >> (64 - value_bits) & !t_mask);
            let subslice = subslice << (Self::t_bits() - back_bits);
            let end_block = end_middle_block + 1;

            self.storage[end_block] = (self.storage[end_block] & back_gap_mask) | subslice;
        } */
        // turn the value slice into a 32 b k-mer
        let slice = IntKmer::<u64>::from_u64(value);

        // iterate over bases which should be transferred
        for i in  0..n_bases {
            let val = slice.get(i);
            self.set_mut(i + pos, val);
        }
    }

    /// Return the reverse complement of this kmer
    fn rc(&self) -> Self {
        let mut new_kmer = VarLenKmer::empty();

        // get the rc of each base from back to front and add into a new k-mer
        for rc_base in (0..Self::k()).rev().map(|pos| 3 - self.get(pos)) {
            new_kmer = new_kmer.extend_right(rc_base);
        }
        new_kmer
    }

    fn at_count(&self) -> u32 {
        // A's and T's have upper_bit ^ lower_bit == 0
        // count how many of these are present
        
        // first block with mask on unused bits
        let mix_base_bits = !((self.storage[0] >> 1) ^ self.storage[0]);
        let mask_lower = mix_base_bits & !Self::unused_bits_mask(0) & IntHelp::lower_of_two();
        let mut count = mask_lower.count_ones();

        for block in 1..LEN {
            let mix_base_bits = !((self.storage[block] >> 1) ^ self.storage[block]);
            let mask_lower = mix_base_bits & IntHelp::lower_of_two();
            count += mask_lower.count_ones();
        }

        count
    }

    fn gc_count(&self) -> u32 {
        // A's and T's have upper_bit ^ lower_bit == 1
        // count how many of these are present
        let mix_base_bits = (self.storage[0] >> 1) ^ self.storage[0];
        let mask_lower = mix_base_bits & !Self::unused_bits_mask(0) & IntHelp::lower_of_two();
        let mut count = mask_lower.count_ones();
    
        for block in 1..LEN {
            let mix_base_bits = (self.storage[block] >> 1) ^ self.storage[block];
            let mask_lower = mix_base_bits & IntHelp::lower_of_two();
            count += mask_lower.count_ones();
        }

        count
    }
}

impl<T: IntHelp + DeserializeOwned, const LEN: usize, KS: KmerSize> fmt::Debug for VarLenKmer<T, LEN, KS> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut s = String::new();
        for pos in 0..Self::k() {
            s.push(bits_to_base(self.get(pos)))
        }

        write!(f, "{}", s)
    }
}

/// Marker struct for generating K=128 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K128;

/// Marker struct for generating K=112 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K112;

/// Marker struct for generating K=96 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K96;

/// Marker struct for generating K=80 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K80;

/// Marker struct for generating K=79 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K79;
/// Marker struct for generating K=63 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K63;

/// Marker struct for generating K=62 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K62;

/// Marker struct for generating K=61 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K61;

/// Marker struct for generating K=60 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K60;

/// Marker struct for generating K=59 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K59;

/// Marker struct for generating K=58 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K58;

/// Marker struct for generating K=57 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K57;

/// Marker struct for generating K=56 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K56;

/// Marker struct for generating K=55 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K55;

/// Marker struct for generating K=54 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K54;

/// Marker struct for generating K=53 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K53;

/// Marker struct for generating K=52 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K52;

/// Marker struct for generating K=51 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K51;

/// Marker struct for generating K=50 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K50;

/// Marker struct for generating K=49 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K49;

/// Marker struct for generating K=48 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K48;

/// Marker trait for generating K=47 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K47;

/// Marker trait for generating K=46 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K46;

/// Marker trait for generating K=45 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K45;

/// Marker trait for generating K=44 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K44;

/// Marker trait for generating K=43 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K43;

/// Marker trait for generating K=42 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K42;

/// Marker trait for generating K=41 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K41;

/// Marker trait for generating K=40 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K40;

/// Marker trait for generating K=39 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K39;

/// Marker trait for generating K=38 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K38;

/// Marker trait for generating K=37 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K37;

/// Marker trait for generating K=36 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K36;

/// Marker trait for generating K=35 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K35;

/// Marker trait for generating K=34 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K34;

/// Marker trait for generating K=33 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K33;

/// Marker trait for generating K=32 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K32;
/// Marker trait for generating K=31 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K31;

/// Marker trait for generating K=30 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K30;

/// Marker trait for generating K=29 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K29;

/// Marker trait for generating K=28 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K28;

/// Marker trait for generating K=27 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K27;

/// Marker trait for generating K=26 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K26;

/// Marker trait for generating K=25 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K25;

/// Marker trait for generating K=24 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K24;

/// Marker trait for generating K=23 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K23;

/// Marker trait for generating K=22 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K22;

/// Marker trait for generating K=21 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K21;

/// Marker trait for generating K=20 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K20;

/// Marker trait for generating K=19 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K19;

/// Marker trait for generating K=18 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K18;

/// Marker trait for generating K=17 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K17;

/// Marker trait for generating K=15 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K15;

/// Marker trait for generating K=14 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K14;

/// Marker trait for generating K=13 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K13;

/// Marker trait for generating K=12 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K12;

/// Marker trait for generating K=11 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K11;

/// Marker trait for generating K=10 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K10;

/// Marker trait for generating K=9 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K9;

/// Marker trait for generating K=7 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K7;

/// Marker trait for generating K=6 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K6;

/// Marker trait for generating K=5 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K5;
/// Marker trait for generating K=4 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K4;
/// Marker trait for generating K=3 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K3;
/// Marker trait for generating K=2 Kmers
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
pub struct K2;

/// Marker trait for generating K=0 Kmers (for testing)
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
struct K0;

/// Marker trait for generating K=300 Kmers (for testing)
#[derive(Debug, Hash, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, KmerSize)]
struct K300;

#[cfg(test)]
mod tests {
use std::fmt::Debug;
use std::time;

use super::*;
    use crate::vmer::Lmer;
    use rand::{self, Rng, RngCore};

    use crate::MerImmut;
    use crate::Vmer;

    fn check_hd<T: Kmer + Debug>(k1: T, k2: T) {
        let mut n = 0;
        for i in 0..T::k() {
            if k1.get(i) != k2.get(i) {
                n += 1;
            }
        }

        assert_eq!(k1.hamming_dist(k2), n);
        assert_eq!(k2.hamming_dist(k1), n);
    }

    // Generate random kmers & test the methods for manipulating them
    fn check_kmer<T: Kmer + Debug>() {
        #[allow(non_snake_case)]
        let K = T::k();

        // Random kmer
        let km = random_kmer::<T>();

        // check HD calc
        let km2 = random_kmer::<T>();
        check_hd(km, km2);

        // Reverse complementing
        let rc = km.rc();
        let double_rc = rc.rc();
        assert!(km == double_rc);

        for i in 0..K {
            assert!(km.get(i) == (3 - rc.get(K - 1 - i)))
        }

        // Get and set bases
        for i in 0..K {
            let km2 = km.set(i, 0);
            assert!(km2.get(i) == 0);
        }
        let mut copy_kmer = T::empty();
        for i in 0..K {
            copy_kmer = copy_kmer.set(i, km.get(i));
        }
        assert!(km == copy_kmer);

        // Extend right
        let nb = random_base();
        let ext_r = km.extend_right(nb);
        assert!(ext_r.get(K - 1) == nb);
        assert!(km.get(1) == ext_r.get(0));

        // Extend Left
        let nb = random_base();
        let ext_l = km.extend_left(nb);
        assert!(ext_l.get(0) == nb);
        assert!(ext_l.get(1) == km.get(0));

        // Extend left with an A -- should sort lower
        let lt = km.extend_left(0);
        assert!(lt <= km);

        // Extend left with T -- should sort higher
        let gt = km.extend_left(3);
        assert!(gt >= km);

        // Shift twice
        let l_base = random_base();
        let r_base = random_base();
        let ts = km.set(0, l_base).set(K - 1, r_base);

        let double_shift = ts.extend_left(0).extend_right(r_base);
        assert!(ts == double_shift);

        if T::k() <= 32 {
            // Convert to and from u64.
            let u64_1 = km.to_u64();
            let km2 = T::from_u64(u64_1);
            let u64_2 = km2.to_u64();
            assert_eq!(km, km2);
            assert_eq!(u64_1, u64_2);
        }

        // check AT / GC counter
        let mut at_count = 0;
        let mut gc_count = 0;
        for i in 0..km.len() {
            let base = km.get(i);
            if base == 0 || base == 3 {
                at_count += 1;
            } else {
                gc_count += 1;
            }
        }

        assert_eq!(km.at_count(), at_count);
        assert_eq!(km.gc_count(), gc_count);
    }

    fn check_vmer<V: Vmer + MerImmut, T: Kmer>() {
        let vm = random_vmer::<V, T>();
        let l = vm.len();

        let rc = vm.rc();

        for i in 0..l {
            if vm.get(i) != (3 - rc.get(l - 1 - i)) {
                println!("km: {:?}, rc: {:?}", vm, rc);
            }

            assert!(vm.get(i) == (3 - rc.get(l - 1 - i)))
        }

        let double_rc = rc.rc();
        assert!(vm == double_rc);

        for i in 0..l {
            // Get and set
            let vm2 = vm.set(i, 0);
            assert!(vm2.get(i) == 0);
        }

        let mut copy_vmer = V::new(l);
        for i in 0..l {
            copy_vmer = copy_vmer.set(i, vm.get(i));
        }
        assert!(vm == copy_vmer);

        let kmers: Vec<T> = vm.iter_kmers().collect();
        assert_eq!(kmers.len(), vm.len() - T::k() + 1);

        for (pos, kmer) in kmers.iter().enumerate() {
            for i in 0..T::k() {
                assert_eq!(kmer.get(i), vm.get(i + pos))
            }

            if vm.get_kmer::<T>(pos) != *kmer {
                println!(
                    "didn't get same kmer: i:{}, vm: {:?} kmer_iter: {:?}, kmer_get: {:?}",
                    pos,
                    vm,
                    kmer,
                    vm.get_kmer::<T>(pos)
                )
            }
            assert_eq!(*kmer, vm.get_kmer(pos));
        }

        assert!(kmers[0] == vm.first_kmer());
        assert!(kmers[kmers.len() - 1] == vm.last_kmer());
    }

    pub fn random_vmer<V: Vmer + MerImmut, T: Kmer>() -> V {
        let mut r = rand::thread_rng();
        let len = r.gen_range(T::k(), V::max_len());

        let mut vmer = V::new(len);
        for pos in 0..len {
            let b = (r.next_u64() % 4) as u8;
            vmer = vmer.set(pos, b);
        }
        vmer
    }

    pub fn random_kmer<T: Kmer>() -> T {
        let mut r = rand::thread_rng();
        let mut kmer = T::empty();
        for pos in 0..T::k() {
            let b = (r.next_u64() % 4) as u8;
            kmer = kmer.set(pos, b);
        }
        kmer
    }

    pub fn random_base() -> u8 {
        let mut r = rand::thread_rng();
        (r.next_u64() % 4) as u8
    }

    #[test]
    fn test_lmer_3_kmer_64() {
        for _ in 0..10000 {
            check_vmer::<Lmer<[u64; 3]>, IntKmer<u128>>();
        }
    }

    #[test]
    fn test_lmer_3_kmer_48() {
        for _ in 0..10000 {
            check_vmer::<Lmer<[u64; 3]>, VarIntKmer<u128, K48>>();
        }
    }

    #[test]
    fn test_lmer_3_kmer_32() {
        for _ in 0..10000 {
            check_vmer::<Lmer<[u64; 3]>, IntKmer<u64>>();
        }
    }

    #[test]
    fn test_lmer_2_kmer_32() {
        for _ in 0..10000 {
            check_vmer::<Lmer<[u64; 3]>, IntKmer<u64>>();
        }
    }

    #[test]
    fn test_lmer_1_kmer_24() {
        for _ in 0..10000 {
            check_vmer::<Lmer<[u64; 1]>, VarIntKmer<u64, K24>>();
        }
    }

    #[test]
    fn test_lmer_1_kmer_20() {
        for _ in 0..10000 {
            check_vmer::<Lmer<[u64; 1]>, VarIntKmer<u64, K20>>();
        }
    }

    #[test]
    fn test_lmer_1_kmer_16() {
        for _ in 0..10000 {
            check_vmer::<Lmer<[u64; 1]>, IntKmer<u32>>();
        }
    }

    #[test]
    fn test_kmer_128() {
        for _ in 0..10000 {
            check_kmer::<VarLenKmer<u64, 4, K128>>();
        }
    }

    #[test]
    fn test_kmer_112() {
        for _ in 0..10000 {
            check_kmer::<VarLenKmer<u32, 7, K112>>();
        }
    }

    #[test]
    fn test_kmer_96() {
        for _ in 0..10000 {
            check_kmer::<VarLenKmer<u64, 3, K96>>();
        }
    }

    #[test]
    fn test_kmer_80() {
        for _ in 0..10000 {
            check_kmer::<VarLenKmer<u32, 5, K80>>();
        }
    }

    #[test]
    fn test_kmer_79() {
        for _ in 0..10000 {
            check_kmer::<VarLenKmer<u32, 5, K79>>();
        }
    }

    #[test]
    fn test_kmer_64() {
        for _ in 0..10000 {
            check_kmer::<IntKmer<u128>>();
        }
    }

    #[test]
    fn test_kmer_63() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K63>>();
        }
    }

    #[test]
    fn test_kmer_62() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K62>>();
        }
    }

    #[test]
    fn test_kmer_61() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K61>>();
        }
    }

    #[test]
    fn test_kmer_60() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K60>>();
        }
    }

    #[test]
    fn test_kmer_59() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K59>>();
        }
    }

    #[test]
    fn test_kmer_58() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K58>>();
        }
    }


    #[test]
    fn test_kmer_57() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K57>>();
        }
    }

    #[test]
    fn test_kmer_56() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K56>>();
        }
    }


    #[test]
    fn test_kmer_55() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K55>>();
        }
    }

    #[test]
    fn test_kmer_54() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K54>>();
        }
    }

    #[test]
    fn test_kmer_53() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K53>>();
        }
    }

    #[test]
    fn test_kmer_52() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K52>>();
        }
    }

    #[test]
    fn test_kmer_51() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K51>>();
        }
    }

    #[test]
    fn test_kmer_50() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K50>>();
        }
    }

    #[test]
    fn test_kmer_49() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K49>>();
        }
    }


    #[test]
    fn test_kmer_48() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K48>>();
        }
    }

    #[test]
    fn test_kmer_47() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K47>>();
        }
    }

    #[test]
    fn test_kmer_46() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K46>>();
        }
    }

    #[test]
    fn test_kmer_45() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K45>>();
        }
    }

    #[test]
    fn test_kmer_44() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K44>>();
        }
    }

    #[test]
    fn test_kmer_43() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K43>>();
        }
    }

    #[test]
    fn test_kmer_42() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K42>>();
        }
    }

    #[test]
    fn test_kmer_41() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K41>>();
        }
    }

    #[test]
    fn test_kmer_40() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K40>>();
        }
    }

    #[test]
    fn test_kmer_39() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K39>>();
        }
    }

    #[test]
    fn test_kmer_38() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K38>>();
        }
    }

    #[test]
    fn test_kmer_37() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K37>>();
        }
    }

    #[test]
    fn test_kmer_36() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K36>>();
        }
    }

    #[test]
    fn test_kmer_35() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K35>>();
        }
    }

    #[test]
    fn test_kmer_34() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K34>>();
        }
    }

    #[test]
    fn test_kmer_33() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u128, K33>>();
        }
    }

    #[test]
    fn test_kmer_32() {
        for _ in 0..10000 {
            check_kmer::<IntKmer<u64>>();
        }
    }

    #[test]
    fn test_kmer_31() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K31>>();
        }
    }

    #[test]
    fn test_kmer_30() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K30>>();
        }
    }

    #[test]
    fn test_kmer_29() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K29>>();
        }
    }

    #[test]
    fn test_kmer_28() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K28>>();
        }
    }

    #[test]
    fn test_kmer_27() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K27>>();
        }
    }

    #[test]
    fn test_kmer_26() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K26>>();
        }
    }

    #[test]
    fn test_kmer_25() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K25>>();
        }
    }

    #[test]
    fn test_kmer_24() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K24>>();
        }
    }

    #[test]
    fn test_kmer_23() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K23>>();
        }
    }

    #[test]
    fn test_kmer_22() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K22>>();
        }
    }

    #[test]
    fn test_kmer_21() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K21>>();
        }
    }

    #[test]
    fn test_kmer_20() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K20>>();
        }
    }

    #[test]
    fn test_kmer_19() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K19>>();
        }
    }

    #[test]
    fn test_kmer_18() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K18>>();
        }
    }

    #[test]
    fn test_kmer_17() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u64, K17>>();
        }
    }

    #[test]
    fn test_kmer_16() {
        for _ in 0..10000 {
            check_kmer::<IntKmer<u32>>();
        }
    }

    #[test]
    fn test_kmer_15() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u32, K15>>();
        }
    }

    #[test]
    fn test_kmer_14() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u32, K14>>();
        }
    }

    #[test]
    fn test_kmer_13() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u32, K13>>();
        }
    }

    #[test]
    fn test_kmer_12() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u32, K12>>();
        }
    }

    #[test]
    fn test_kmer_11() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u32, K11>>();
        }
    }

    #[test]
    fn test_kmer_10() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u32, K10>>();
        }
    }

    #[test]
    fn test_kmer_9() {
        for _ in 0..10000 {
            check_kmer::<VarIntKmer<u32, K9>>();
        }
    }

    #[test]
    fn test_kmer_8() {
        for _ in 0..10000 {
            check_kmer::<IntKmer<u16>>();
        }
    }

    #[test]
    fn test_kmer_7() {
        for _ in 0..10000 {
            check_kmer::<Kmer7>();
        }
    }

    #[test]
    fn test_kmer_6() {
        for _ in 0..10000 {
            check_kmer::<Kmer6>();
        }
    }

    #[test]
    fn test_kmer_5() {
        for _ in 0..10000 {
            check_kmer::<Kmer5>();
        }
    }

    #[test]
    fn test_kmer_4() {
        for _ in 0..10000 {
            check_kmer::<Kmer4>();
        }
    }

    #[test]
    fn test_kmer_3() {
        for _ in 0..10000 {
            check_kmer::<Kmer3>();
        }
    }

    #[test]
    fn test_kmer_2() {
        for _ in 0..10000 {
            check_kmer::<Kmer2>();
        }
    }

    #[test]
    fn test_is_empty() {
        let kmer = IntKmer::<u8>::empty();
        assert!(!kmer.is_empty()); // cant be empty
        let kmer = VarIntKmer::<u8, K0>::empty();
        assert!(kmer.is_empty());
        let kmer = VarLenKmer::<u8, 0, K0>::empty();
        assert!(kmer.is_empty());
    }

    #[test]
    fn test_extend_left() {
        let mut kmer = VarLenKmer::<u32, 5, K80>::empty();

        for base in [1, 1, 1, 3, 0, 1, 2, 0, 3, 0, 1].iter() {
            kmer = kmer.extend_left(*base);
        }

        println!("kmer: {:?}", kmer);

        let mut kmer = VarLenKmer::<u32, 5, K80>::empty();

        for base in [1, 1, 1, 3, 0, 1, 2, 0, 3, 0, 1].iter().rev() {
            kmer = kmer.extend_right(*base);
        }

        println!("kmer: {:?}", kmer);
    }

    #[test]
    fn test_set_mut_slice() {
        let mut kmer = VarLenKmer::<u32, 5, K80>::empty();
        kmer.set_slice_mut(3, 28, 0b01110111010101011111111101010101111111110000000010101010 << 8);
        println!("{:?}", kmer);
    }


    #[test]
    fn test_varlenkmer_speed() {
        let bef = time::Instant::now();
        for _i in 0..100 {
            let mut kmer = VarLenKmer::<u32, 2, K31>::empty();
            kmer.set_slice_mut(6, 15, 0b010101010101010101010101010101);
        }
        let af = bef.elapsed().as_nanos();
        println!("set slice mut, var len 31: {af}");

        let bef = time::Instant::now();
        for _i in 0..100 {
            let mut kmer = Kmer31::empty();
            kmer.set_slice_mut(6, 15, 0b010101010101010101010101010101);
        }
        let af = bef.elapsed().as_nanos();
        println!("set slice mut, con len 31: {af}");

        let bef = time::Instant::now();
        for _i in 0..100 {
            let mut kmer = VarLenKmer::<u32, 2, K31>::empty();
            kmer.set_mut(5, 1);
        }
        let af = bef.elapsed().as_nanos();
        println!("set mut, var len 31: {af}");

        let bef = time::Instant::now();
        for _i in 0..100 {
            let mut kmer = Kmer31::empty();
            kmer.set_mut(5, 1);
        }
        let af = bef.elapsed().as_nanos();
        println!("set mut, con len 31: {af}");

        let bef = time::Instant::now();
        for _i in 0..100 {
            let kmer = Kmer64::from_u64(17);
            for i in 6..21 {
                kmer.get(i);
            }
        }
        let af = bef.elapsed().as_nanos();
        println!("from u64, get f u64: {af}"); 
    }
}