diff --git a/Project.toml b/Project.toml
index 83a9676e..9896a428 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,6 +13,7 @@ CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 CodecZstd = "6b39b394-51ab-5f42-8807-6242bab2b4c2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+GeneticVariantBase = "2447270c-d849-4bf9-ac0d-b5c0b265991c"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
@@ -56,4 +57,4 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "CUDA"]
+test = ["Test"]
diff --git a/src/SnpArrays.jl b/src/SnpArrays.jl
index 4fdfb86f..9b9255fc 100644
--- a/src/SnpArrays.jl
+++ b/src/SnpArrays.jl
@@ -2,6 +2,7 @@ __precompile__()
 
 module SnpArrays
 
+using GeneticVariantBase
 using CodecZlib, CodecXz, CodecBzip2, CodecZstd,  TranscodingStreams
 using Adapt, Glob, LinearAlgebra, LoopVectorization, Missings, Mmap, Printf
 using Requires, SparseArrays, Statistics, StatsBase
@@ -18,7 +19,7 @@ import Tables: table
 export AbstractSnpArray, AbstractSnpBitMatrix, AbstractSnpLinAlg
 export SnpArray, SnpBitMatrix, SnpLinAlg, SnpData, StackedSnpArray
 export compress_plink, decompress_plink, split_plink, merge_plink, write_plink 
-export counts, grm, grm_admixture, maf, mean, minorallele, missingpos, missingrate
+export counts, grm, grm_admixture, maf, maf!, mean, minorallele, missingpos, missingrate
 export std, var, vcf2plink
 export counts, grm, maf, mean, minorallele, missingpos, missingrate, std, var
 export vcf2plink, kinship_pruning
@@ -26,6 +27,8 @@ export ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL
 export CuSnpArray
 import VariantCallFormat: findgenokey, VCF, header
 
+# this is exporting functions necessary 
+
 const ADDITIVE_MODEL = Val(1)
 const DOMINANT_MODEL = Val(2)
 const RECESSIVE_MODEL = Val(3)
@@ -43,6 +46,8 @@ include("linalg_bitmatrix.jl")
 include("reorder.jl")
 include("vcf2plink.jl")
 include("admixture.jl")
+include("iterator.jl")
+
 AbstractSnpArray = Union{SnpArray, SubArray{UInt8, 1, SnpArray}, SubArray{UInt8, 2, SnpArray}, 
     StackedSnpArray, SubArray{UInt8, 1, StackedSnpArray}, SubArray{UInt8, 2, StackedSnpArray}}
 
diff --git a/src/SnpArrays.jl.31595.mem b/src/SnpArrays.jl.31595.mem
new file mode 100644
index 00000000..4c4836c4
--- /dev/null
+++ b/src/SnpArrays.jl.31595.mem
@@ -0,0 +1,60 @@
+        - __precompile__()
+        - 
+ 14681029 module SnpArrays
+        - 
+        - using GeneticVariantBase
+        - using CodecZlib, CodecXz, CodecBzip2, CodecZstd,  TranscodingStreams
+        - using Adapt, Glob, LinearAlgebra, LoopVectorization, Missings, Mmap, Printf
+        - using Requires, SparseArrays, Statistics, StatsBase
+        - import Base: IndexStyle, convert, copyto!, eltype, getindex, setindex!, length, size, wait
+        - import DataFrames: DataFrame, rename!, eachrow
+        - import DelimitedFiles: readdlm, writedlm
+        - import CSV # for CSV.read, to avoid clash with Base.read
+        - import LinearAlgebra: copytri!, mul!
+        - import Statistics: mean, std, var
+        - import StatsBase: counts
+        - import SpecialFunctions: gamma_inc
+        - import VectorizationBase: gesp
+        - import Tables: table
+        - export AbstractSnpArray, AbstractSnpBitMatrix, AbstractSnpLinAlg
+        - export SnpArray, SnpBitMatrix, SnpLinAlg, SnpData, StackedSnpArray
+        - export compress_plink, decompress_plink, split_plink, merge_plink, write_plink 
+        - export counts, grm, grm_admixture, maf, maf!, mean, minorallele, missingpos, missingrate
+        - export std, var, vcf2plink
+        - export counts, grm, maf, mean, minorallele, missingpos, missingrate, std, var
+        - export vcf2plink, kinship_pruning
+        - export ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL
+        - export CuSnpArray
+        - import VariantCallFormat: findgenokey, VCF, header
+        - 
+        - # this is exporting functions necessary 
+        - 
+        - const ADDITIVE_MODEL = Val(1)
+        - const DOMINANT_MODEL = Val(2)
+        - const RECESSIVE_MODEL = Val(3)
+        - 
+        - include("codec.jl")
+        - include("snparray.jl")
+        - include("stackedsnparray.jl")
+        - include("filter.jl")
+        - include("cat.jl")
+        - include("snpdata.jl")
+        - include("grm.jl")
+        - include("kinship_pruning.jl")
+        - include("linalg_direct.jl")
+        - include("linalg_bitmatrix.jl")
+        - include("reorder.jl")
+        - include("vcf2plink.jl")
+        - include("admixture.jl")
+        - include("iterator.jl")
+        - 
+        - AbstractSnpArray = Union{SnpArray, SubArray{UInt8, 1, SnpArray}, SubArray{UInt8, 2, SnpArray}, 
+        -     StackedSnpArray, SubArray{UInt8, 1, StackedSnpArray}, SubArray{UInt8, 2, StackedSnpArray}}
+        - 
+        - datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)
+        - 
+        - function __init__()
+        -     @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
+        - end
+        - 
+        - end # module
diff --git a/src/SnpArrays.jl.31601.mem b/src/SnpArrays.jl.31601.mem
new file mode 100644
index 00000000..50bcb2ea
--- /dev/null
+++ b/src/SnpArrays.jl.31601.mem
@@ -0,0 +1,60 @@
+        - __precompile__()
+        - 
+        - module SnpArrays
+        - 
+        - using GeneticVariantBase
+        - using CodecZlib, CodecXz, CodecBzip2, CodecZstd,  TranscodingStreams
+        - using Adapt, Glob, LinearAlgebra, LoopVectorization, Missings, Mmap, Printf
+        - using Requires, SparseArrays, Statistics, StatsBase
+        - import Base: IndexStyle, convert, copyto!, eltype, getindex, setindex!, length, size, wait
+        - import DataFrames: DataFrame, rename!, eachrow
+        - import DelimitedFiles: readdlm, writedlm
+        - import CSV # for CSV.read, to avoid clash with Base.read
+        - import LinearAlgebra: copytri!, mul!
+        - import Statistics: mean, std, var
+        - import StatsBase: counts
+        - import SpecialFunctions: gamma_inc
+        - import VectorizationBase: gesp
+        - import Tables: table
+        - export AbstractSnpArray, AbstractSnpBitMatrix, AbstractSnpLinAlg
+        - export SnpArray, SnpBitMatrix, SnpLinAlg, SnpData, StackedSnpArray
+        - export compress_plink, decompress_plink, split_plink, merge_plink, write_plink 
+        - export counts, grm, grm_admixture, maf, maf!, mean, minorallele, missingpos, missingrate
+        - export std, var, vcf2plink
+        - export counts, grm, maf, mean, minorallele, missingpos, missingrate, std, var
+        - export vcf2plink, kinship_pruning
+        - export ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL
+        - export CuSnpArray
+        - import VariantCallFormat: findgenokey, VCF, header
+        - 
+        - # this is exporting functions necessary 
+        - 
+        - const ADDITIVE_MODEL = Val(1)
+        - const DOMINANT_MODEL = Val(2)
+        - const RECESSIVE_MODEL = Val(3)
+        - 
+        - include("codec.jl")
+        - include("snparray.jl")
+        - include("stackedsnparray.jl")
+        - include("filter.jl")
+        - include("cat.jl")
+        - include("snpdata.jl")
+        - include("grm.jl")
+        - include("kinship_pruning.jl")
+        - include("linalg_direct.jl")
+        - include("linalg_bitmatrix.jl")
+        - include("reorder.jl")
+        - include("vcf2plink.jl")
+        - include("admixture.jl")
+        - include("iterator.jl")
+        - 
+        - AbstractSnpArray = Union{SnpArray, SubArray{UInt8, 1, SnpArray}, SubArray{UInt8, 2, SnpArray}, 
+        -     StackedSnpArray, SubArray{UInt8, 1, StackedSnpArray}, SubArray{UInt8, 2, StackedSnpArray}}
+        - 
+        - datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)
+        - 
+        - function __init__()
+        0     @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
+        - end
+        - 
+        - end # module
diff --git a/src/SnpArrays.jl.31603.mem b/src/SnpArrays.jl.31603.mem
new file mode 100644
index 00000000..50bcb2ea
--- /dev/null
+++ b/src/SnpArrays.jl.31603.mem
@@ -0,0 +1,60 @@
+        - __precompile__()
+        - 
+        - module SnpArrays
+        - 
+        - using GeneticVariantBase
+        - using CodecZlib, CodecXz, CodecBzip2, CodecZstd,  TranscodingStreams
+        - using Adapt, Glob, LinearAlgebra, LoopVectorization, Missings, Mmap, Printf
+        - using Requires, SparseArrays, Statistics, StatsBase
+        - import Base: IndexStyle, convert, copyto!, eltype, getindex, setindex!, length, size, wait
+        - import DataFrames: DataFrame, rename!, eachrow
+        - import DelimitedFiles: readdlm, writedlm
+        - import CSV # for CSV.read, to avoid clash with Base.read
+        - import LinearAlgebra: copytri!, mul!
+        - import Statistics: mean, std, var
+        - import StatsBase: counts
+        - import SpecialFunctions: gamma_inc
+        - import VectorizationBase: gesp
+        - import Tables: table
+        - export AbstractSnpArray, AbstractSnpBitMatrix, AbstractSnpLinAlg
+        - export SnpArray, SnpBitMatrix, SnpLinAlg, SnpData, StackedSnpArray
+        - export compress_plink, decompress_plink, split_plink, merge_plink, write_plink 
+        - export counts, grm, grm_admixture, maf, maf!, mean, minorallele, missingpos, missingrate
+        - export std, var, vcf2plink
+        - export counts, grm, maf, mean, minorallele, missingpos, missingrate, std, var
+        - export vcf2plink, kinship_pruning
+        - export ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL
+        - export CuSnpArray
+        - import VariantCallFormat: findgenokey, VCF, header
+        - 
+        - # this is exporting functions necessary 
+        - 
+        - const ADDITIVE_MODEL = Val(1)
+        - const DOMINANT_MODEL = Val(2)
+        - const RECESSIVE_MODEL = Val(3)
+        - 
+        - include("codec.jl")
+        - include("snparray.jl")
+        - include("stackedsnparray.jl")
+        - include("filter.jl")
+        - include("cat.jl")
+        - include("snpdata.jl")
+        - include("grm.jl")
+        - include("kinship_pruning.jl")
+        - include("linalg_direct.jl")
+        - include("linalg_bitmatrix.jl")
+        - include("reorder.jl")
+        - include("vcf2plink.jl")
+        - include("admixture.jl")
+        - include("iterator.jl")
+        - 
+        - AbstractSnpArray = Union{SnpArray, SubArray{UInt8, 1, SnpArray}, SubArray{UInt8, 2, SnpArray}, 
+        -     StackedSnpArray, SubArray{UInt8, 1, StackedSnpArray}, SubArray{UInt8, 2, StackedSnpArray}}
+        - 
+        - datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)
+        - 
+        - function __init__()
+        0     @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
+        - end
+        - 
+        - end # module
diff --git a/src/iterator.jl b/src/iterator.jl
new file mode 100644
index 00000000..27659a4c
--- /dev/null
+++ b/src/iterator.jl
@@ -0,0 +1,163 @@
+mutable struct SnpArrayIterator <: GeneticVariantBase.VariantIterator
+   snpdata::SnpData
+end
+
+mutable struct SnpArrayIndex <: GeneticVariantBase.Variant
+    index::Int
+end 
+
+@inline function Base.eltype(::Type{<:VariantIterator})
+    SnpArrayIndex
+end
+
+function Base.iterate(itr::SnpArrayIterator, state=1)
+    if state <= 0
+        throw(BoundsError(itr, state))
+    end
+    if state > size(itr.snpdata.snparray,2)
+        return nothing
+    else
+        index = SnpArrayIndex(state)
+        state = state + 1
+        return (index, state)
+    end
+end
+
+@inline function Base.length(itr::SnpArrayIterator)
+    return size(itr.snpdata.snparray, 2)
+end
+
+function GeneticVariantBase.chrom(s::SnpData, snpindex::SnpArrayIndex)::String
+    result = s.snp_info[snpindex.index,:chromosome]
+    return result
+end
+
+function GeneticVariantBase.pos(s::SnpData, snpindex::SnpArrayIndex)::Int
+    result = s.snp_info[snpindex.index,:position]
+    # println("entered pos function $result $snpindex.index")
+    return result
+end
+
+function GeneticVariantBase.rsid(s::SnpData, snpindex::SnpArrayIndex)::String
+    result = s.snp_info[snpindex.index,:snpid]
+    return result
+end
+
+#SnpData subtype of Genetic Data
+
+function alleles(s::SnpData, snpindex::SnpArrayIndex)::Vector{String}
+    allele1 = s.snp_info[snpindex.index,:allele1]
+    allele2 = s.snp_info[snpindex.index,:allele2]
+    return [allele1, allele2]
+end
+
+function GeneticVariantBase.alt_allele(s::SnpData, snpindex::SnpArrayIndex)::String
+    alt = s.snp_info[snpindex.index,:allele2]
+    return alt
+end
+
+function GeneticVariantBase.ref_allele(s::SnpData, snpindex::SnpArrayIndex)::String
+    ref = s.snp_info[snpindex.index,:allele1]
+    return ref
+end
+
+struct MAFData
+    maf_vector::Vector{Float64}
+end
+
+# fold into GeneticVariantBase.maf function name 
+
+# function calculate_maf_data(s::SnpData)
+#     maf_vector = maf(s.snparray)
+#     result = MAFData(maf_vector)
+#     return result 
+# end
+
+# function maf_index(maf_data::MAFData, snpindex::SnpArrayIndex)
+#     return maf_data.maf_vector[snpindex.index]
+# end
+
+# function GeneticVariantBase.maf(s::SnpData, snpindex::SnpArrayIndex) # This is calculating maf every call 
+#     # maf_vector = calculate_maf_data(s)
+#     maf_vector = maf(s.snparray) # This allocates a lot
+#     # Calculate MAF for one variant if the index is given you might be able to calculate 
+#     # Boolean vector indicating which samples to be used default argument selecting which 
+#     return maf_vector[snpindex.index] 
+#     # return maf_vector[snpindex.index]
+# end 
+
+# Global cache: maps SnpData => MAF vector
+const _maf_cache = IdDict{SnpData, Vector{Float64}}()
+
+"""
+    calculate_maf_data!(s::SnpData)
+
+Precompute and cache MAF values for `s`.
+"""
+function calculate_maf_data!(s::SnpData)
+    _maf_cache[s] = maf(s.snparray)
+    return nothing
+end
+
+"""
+    maf_index(s::SnpData, snpindex::SnpArrayIndex)
+
+Fetches the MAF for `snpindex` from the cache.
+Assumes `calculate_maf_data!` has been called beforehand.
+"""
+function maf_index(s::SnpData, snpindex::SnpArrayIndex)
+    maf_vector = _maf_cache[s]
+    return maf_vector[snpindex.index]
+end
+
+# Override without changing args
+function GeneticVariantBase.maf(s::SnpData, snpindex::SnpArrayIndex)
+    maf_vector = get!(_maf_cache, s) do
+        # compute and store if missing
+        maf(s.snparray)
+    end
+    return maf_vector[snpindex.index]
+end
+
+
+
+function GeneticVariantBase.hwepval(s::SnpData, snpindex::SnpArrayIndex)
+    genotypes = s.snparray[:,snpindex.index]
+
+    n00 = sum(genotypes .== 0x00) 
+    n01 = sum(genotypes .== 0x02) 
+    n11 = sum(genotypes .== 0x03) 
+
+   pval = hwe(n00,n01,n11)
+   return pval
+
+end
+
+ # 0 for homozygous allele 1
+    # 2 Heterozygous
+    # 3 homozygous allele 2
+    # 1 is for missing 
+
+function GeneticVariantBase.alt_dosages!(arr::AbstractArray{T}, s::SnpData, snpindex::SnpArrayIndex; mean_impute=true) where T <: Real
+    GeneticVariantBase.alt_genotypes!(arr, s, snpindex; mean_impute=true)
+    return arr 
+end
+
+# make sure you can read in all genotypes for a sample
+# filtering SNPS
+
+
+function GeneticVariantBase.alt_genotypes!(arr::AbstractArray{T}, s::SnpData, snpindex::SnpArrayIndex; mean_impute=true) where T <: Real
+    Base.copyto!(arr, @view(s.snparray[:, snpindex.index]); impute=mean_impute, center=mean_impute, scale=mean_impute)   # change impute to mean_impute 
+    return arr 
+end
+
+# are we reusing the same arr 
+
+function GeneticVariantBase.n_samples(s::SnpData)::Int
+    return size(s.snparray,1)
+end 
+
+function GeneticVariantBase.n_variants(s::SnpData)::Int
+    return size(s.snparray,2)
+end 
\ No newline at end of file
diff --git a/src/snpdata.jl b/src/snpdata.jl
index 8f2708c7..729c6f19 100644
--- a/src/snpdata.jl
+++ b/src/snpdata.jl
@@ -7,7 +7,7 @@ const PERSON_INFO_KEYS = [:fid, :iid, :father, :mother, :sex, :phenotype]
 
 Type to store SNP and person information along with the SnpArray.
 """
-struct SnpData
+struct SnpData <: GeneticData
     people::Int
     snps::Int
     snparray::SnpArray
diff --git a/src/vcf2plink.jl.31595.mem b/src/vcf2plink.jl.31595.mem
new file mode 100644
index 00000000..55aff01f
--- /dev/null
+++ b/src/vcf2plink.jl.31595.mem
@@ -0,0 +1,112 @@
+        - chrom_map = Dict{String, String}()
+        0 for i in 1:23
+    10016     chrom_map[string(i)] = string(i)
+     4600     chrom_map["chr" * string(i)] = string(i)
+        0 end
+        - chrom_map["chrX"] = "X"
+        - chrom_map["chrY"] = "Y"
+        - chrom_map["chrXY"] = "XY"
+        - chrom_map["chrMT"] = "MT"
+        - chrom_map["chrM"] = "M"
+        - chrom_map["X"] = "X"
+        - chrom_map["Y"] = "Y"
+        - chrom_map["XY"] = "XY"
+        - chrom_map["MT"] = "MT"
+        - chrom_map["M"] = "M"
+        - 
+        - @inline function geno_ismissing(record::VCF.Record, range::UnitRange{Int})
+        -     return record.data[first(range)] == UInt8('.') || record.data[last(range)] == UInt8('.')
+        - end
+        - 
+        - """
+        -     vcf2plink(vcffile, plinkprefix)
+        - 
+        - Convert VCF file to PLINK format. This drops multi-allelic variants and variants with missing ID. Reference alleles are set to A1.
+        - """
+        - function vcf2plink(vcffile, plinkprefix)
+        -     @assert endswith(vcffile, "vcf.gz") || endswith(vcffile, ".vcf")
+        -     nsamples = begin
+        -         vcfio = makestream(vcffile, "r")
+        -         reader = VCF.Reader(vcfio)
+        -         samples = length(header(reader).sampleID)
+        -         makestream(plinkprefix * ".fam", "w") do io
+        -             for (i, id) in enumerate(header(reader).sampleID)
+        -                 println(io, "$i\t$id\t0\t0\t0\t0")
+        -             end
+        -         end
+        -         close(reader)
+        -         samples
+        -     end
+        - 
+        -     nsnps = begin
+        -         vcfio = VCF.Reader(SnpArrays.makestream(vcffile, "r"))
+        -         records = 0 
+        -         for record in vcfio
+        -             length(VCF.alt(record)) > 1 && continue
+        -             try
+        -                 ismissing(VCF.id(record))
+        -             catch e
+        -                 continue
+        -             end
+        -             records += 1
+        -         end
+        -         close(vcfio)
+        -         records
+        -     end
+        - 
+        -     makestream(plinkprefix * ".bed", "w+") do io
+        -         write(io, 0x1b6c)
+        -         write(io, 0x01)
+        -         write(io, Matrix{UInt8}(undef, (nsamples + 3) >> 2, nsnps))
+        -     end
+        -     snparray = SnpArray(plinkprefix * ".bed", nsamples, "r+")
+        -     bimio = makestream(plinkprefix * ".bim", "w")
+        -     vcfio = makestream(vcffile, "r")
+        -     reader = VCF.Reader(vcfio)
+        - 
+        -     j = 0
+        -     for record in reader
+        -         length(VCF.alt(record)) > 1 && continue
+        -         try
+        -             ismissing(VCF.id(record))
+        -         catch e
+        -             continue
+        -         end
+        -         j += 1
+        -         chrom = VCF.chrom(record)
+        -         id = VCF.id(record)[1]
+        -         dist = 0
+        -         pos = VCF.pos(record)
+        -         ref = VCF.ref(record)
+        -         alt = VCF.alt(record)[1]
+        -         println(bimio, "$chrom\t$id\t$dist\t$pos\t$ref\t$alt")
+        -         
+        -         gtkey = findgenokey(record, "GT")
+        -         for (i, _) in enumerate(header(reader).sampleID)
+        -             geno = record.genotype[i]
+        -             # dropped field or "." => 0x2e
+        -             if gtkey > lastindex(geno) || geno_ismissing(record, geno[gtkey])
+        -                 snparray[i,j] = 0x01
+        -             else
+        -                 # "0" => 0x30, "1" => 0x31
+        -                 if record.data[geno[gtkey][1]] == 0x30
+        -                     if record.data[geno[gtkey][3]] == 0x30
+        -                         snparray[i, j] = 0x00
+        -                     else
+        -                         snparray[i, j] = 0x02
+        -                     end
+        -                 elseif record.data[geno[gtkey][1]] == 0x31
+        -                     if record.data[geno[gtkey][3]] == 0x31
+        -                         snparray[i, j] = 0x03
+        -                     else
+        -                         snparray[i, j] = 0x02
+        -                     end
+        -                 end
+        -             end
+        -         end
+        -     end
+        -     close(bimio)
+        -     close(vcfio)
+        -     close(reader)
+        -     snparray
+        - end
diff --git a/test/runtests.jl b/test/runtests.jl
index a4fb7abf..d3477955 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,764 +1,880 @@
 using LinearAlgebra, SnpArrays, SparseArrays, Test
 
 const EUR = SnpArray(SnpArrays.datadir("EUR_subset.bed")) # no missing genotypes
-const mouse = SnpArray(SnpArrays.datadir("mouse.bed")) # has missing genotypes
-
-@testset "size" begin
-@test size(EUR) == (379, 54051)
-@test size(EUR, 1) == 379
-@test size(EUR, 2) == 54051
-@test size(mouse) == (1940, 10150)
-@test size(mouse, 3) == 1
-@test_throws ErrorException size(mouse, 0)
-@test length(mouse) == 1940 * 10150
-@test eltype(mouse) == UInt8
-end
-
-@testset "counts" begin
-cc = counts(EUR, dims=1)
-@test counts(EUR) == sum(cc, dims=2)
-@test size(cc) == (4, size(EUR, 2))
-@test view(cc, :, 1) == [2, 0, 70, 307]
-@test all(sum(cc, dims = 1) .== size(EUR, 1))
-@test all(iszero, view(cc, 2, :))
-rc = counts(EUR, dims=2)
-@test size(rc) == (4, size(EUR, 1))
-@test view(rc, :, 1) == [2997, 0, 13143, 37911]
-@test all(sum(rc, dims = 1) .== size(EUR, 2))
-@test all(iszero, view(rc, 2, :))
-@test_throws ArgumentError counts(EUR, dims=3)
-end
-
-@testset "means" begin
-# ADDITIVE_MODEL
-cmns = mean(EUR, dims = 1)
-@test size(cmns) == (1, size(EUR, 2))
-@test cmns[1] ≈ 1.804749340369393
-@test mean(cmns) ≈ mean(EUR)
-@test minimum(cmns) ≈ 1.0
-@test maximum(cmns) ≈ 1.9788918205804749
-rmns = mean(EUR, dims = 2)
-@test size(rmns) == (size(EUR, 1), 1)
-@test rmns[1] ≈ 1.6459454959205195
-@test mean(rmns) ≈ mean(EUR)
-rmnextrma = extrema(rmns)
-@test rmnextrma[1] ≈ 1.637749532848606
-@test rmnextrma[2] ≈ 1.6627259440158368
-@test_throws ArgumentError mean(EUR, dims=3)
-# DOMINANT_MODEL
-cmns = mean(EUR, dims = 1, model = DOMINANT_MODEL)
-@test size(cmns) == (1, size(EUR, 2))
-@test cmns[1] ≈ 0.9947229551451188
-@test mean(cmns) ≈ mean(EUR, model = DOMINANT_MODEL)
-@test minimum(cmns) ≈ 0.5989445910290238
-@test maximum(cmns) ≈ 1.0
-rmns = mean(EUR, dims = 2, model = DOMINANT_MODEL)
-@test size(rmns) == (size(EUR, 1), 1)
-@test rmns[1] ≈ 0.9445523672087472
-@test mean(rmns) ≈ mean(EUR, model = DOMINANT_MODEL)
-rmnextrma = extrema(rmns)
-@test rmnextrma[1] ≈ 0.9325266877578583
-@test rmnextrma[2] ≈ 0.9530073449149877
-@test_throws ArgumentError mean(EUR, dims=3, model=DOMINANT_MODEL)
-# RECESSIVE_MODEL
-cmns = mean(EUR, dims = 1, model = RECESSIVE_MODEL)
-@test size(cmns) == (1, size(EUR, 2))
-@test cmns[1] ≈ 0.8100263852242744
-@test mean(cmns) ≈ mean(EUR, model = RECESSIVE_MODEL)
-@test minimum(cmns) ≈ 0.0
-@test maximum(cmns) ≈ 0.9868073878627969
-rmns = mean(EUR, dims = 2, model = RECESSIVE_MODEL)
-@test size(rmns) == (size(EUR, 1), 1)
-@test rmns[1] ≈ 0.7013931287117722
-@test mean(rmns) ≈ mean(EUR, model = RECESSIVE_MODEL)
-rmnextrma = extrema(rmns)
-@test rmnextrma[1] ≈ 0.6932341677304767
-@test rmnextrma[2] ≈ 0.7191171301178516
-@test_throws ArgumentError mean(EUR, dims=3, model=RECESSIVE_MODEL)
-end
-
-@testset "var/std" begin
-cvars = var(EUR, dims = 1)
-@test size(cvars) == (1, size(EUR, 2))
-@test cvars[1] ≈ 0.16812553224162724
-@test iszero(minimum(cvars))
-@test maximum(cvars) ≈ 0.8705727966941688
-@test std(EUR, dims = 1) ≈ sqrt.(cvars)
-rvars = var(EUR, dims = 2)
-@test size(rvars) == (size(EUR, 1), 1)
-@test rvars[1] ≈ 0.33960146078503206
-@test minimum(rvars) ≈ 0.31972707244268267
-@test maximum(rvars) ≈ 0.365901927927013
-@test std(EUR, dims = 2) ≈ sqrt.(rvars)
-end
-
-@testset "missingpos" begin
-mp = missingpos(mouse)
-cc = counts(mouse, dims=1)
-@test isa(mp, SparseMatrixCSC)
-@test sum(mp, dims = 1) == view(cc, 2:2, :)
-@test sum(mp, dims = 2) == view(counts(mouse, dims=2), 2:2, :)'
-end
-
-@testset "create bed" begin
-tmpbf = SnpArray("tmp.bed", 5, 3)
-@test isfile("tmp.bed")
-@test all(tmpbf .== 0x00)
-fill!(tmpbf, 0x02)
-tmpbf2 = SnpArray("tmp.bed", 5)
-@test all(tmpbf2 .== 0x02)
-Sys.iswindows() || rm("tmp.bed", force=true)
-tmpbf2 = SnpArray("tmp2.bed", SnpArray(undef, 5, 3))
-fill!(tmpbf2, 0x01)
-@test all(tmpbf2 .== 0x01)
-Sys.iswindows() || rm("tmp2.bed", force=true)
-end
-
-@testset "convert" begin
-# check convert coding
-tmpbf = SnpArray(undef, 4, 1)
-tmpbf[1] = 0x00
-tmpbf[2] = 0x01
-tmpbf[3] = 0x02
-tmpbf[4] = 0x03
-# additive model
-v = convert(Matrix{Float64}, tmpbf, model=ADDITIVE_MODEL)
-@test all([v[1], v[3], v[4]] .== [0.0, 1.0, 2.0])
-@test isnan(v[2])
-# dominant model
-v = convert(Matrix{Float64}, tmpbf, model=DOMINANT_MODEL)
-@test all([v[1], v[3], v[4]] .== [0.0, 1.0, 1.0])
-@test isnan(v[2])
-# recessive model
-v = convert(Matrix{Float64}, tmpbf, model=RECESSIVE_MODEL)
-@test all([v[1], v[3], v[4]] .== [0.0, 0.0, 1.0])
-@test isnan(v[2])
-# mouse test data
-cmns = mean(mouse, dims = 1)
-cvrs = var(mouse, dims = 1)
-cstd = std(mouse, dims = 1)
-v = convert(Vector{Float64}, @view(mouse[:, 1]))
-@test isnan(mean(v))
-@test mean(filter(!isnan, v)) ≈ cmns[1]
-@test isnan(var(v))
-@test var(filter(!isnan, v)) ≈ cvrs[1]
-@test isnan(std(v))
-@test std(filter(!isnan, v)) ≈ cstd[1]
-end
-
-@testset "convert (ADMIXTURE)" begin
-# convert A2 allele frequency to dosage
-@test convert(Float64, 0.5, ADDITIVE_MODEL)  == 1.0
-@test convert(Float64, 0.5, DOMINANT_MODEL)  == 0.75
-@test convert(Float64, 0.5, RECESSIVE_MODEL) == 0.25
-end
-
-@testset "maf" begin
-cc = counts(mouse, dims=1)
-@test all(0 .≤ maf(mouse) .≤ 0.5)
-@test all(minorallele(mouse) .== (cc[1, :] .> cc[4, :]))
-cc = counts(EUR, dims=1)
-@test all(0 .≤ maf(EUR) .≤ 0.5)
-@test all(minorallele(EUR) .== (cc[1, :] .> cc[4, :]))
-end
-
-@testset "grm" begin
-Φgrm = grm(EUR, method=:GRM)
-@test size(Φgrm) == (size(EUR, 1), size(EUR, 1))
-@test issymmetric(Φgrm)
-@test eigmin(Φgrm) > -1e-8
-Φmom = grm(EUR, method=:MoM)
-@test size(Φmom) == (size(EUR, 1), size(EUR, 1))
-@test issymmetric(Φmom)
-Φrbs = grm(EUR, method=:Robust)
-@test size(Φrbs) == (size(EUR, 1), size(EUR, 1))
-@test issymmetric(Φrbs)
-@test eigmin(Φrbs) > -1e-8
-m, n, K = size(EUR, 1), size(EUR, 2), 3
-P = rand(n, K)
-Q = rand(m, K)
-Q ./= sum(Q, dims=2)
-Φreap = grm_admixture(EUR, transpose(P), transpose(Q))
-@test size(Φreap) == (size(EUR, 1), size(EUR, 1))
-@test issymmetric(Φreap)
-@test eigmin(Φrbs) > -1e-8
-end
-
-@testset "filter" begin
-rowmask, colmask =  SnpArrays.filter(mouse, min_success_rate_per_row=0.99,
-    min_success_rate_per_col=0.99)
-SnpArrays.filter(SnpArrays.datadir("mouse"), rowmask, colmask; des="tmp")
-tmpbf = SnpArray("tmp.bed")
-@test size(tmpbf) == (1912, 9997)
-@test all(missingrate(tmpbf, 1) .≤ 0.01)
-@test all(missingrate(tmpbf, 2) .≤ 0.01)
-Sys.iswindows() || rm("tmp.bed", force=true)
-rm("tmp.bim", force=true)
-rm("tmp.fam", force=true)
-rowmask, colmask =  SnpArrays.filter(mouse, min_success_rate_per_row=0.99,
-    min_success_rate_per_col=0.99, min_maf=0.01, min_hwe_pval=1e-8)
-@test (count(rowmask), count(colmask)) == (1911, 9510)
-compress_plink(SnpArrays.datadir("mouse"), "gz")
-SnpArrays.filter(SnpArrays.datadir("mouse.bed.gz"), SnpArrays.datadir("mouse.bim.gz"),
-SnpArrays.datadir("mouse.fam.gz"), 1:5, 1:3; des="tmp")
-tmpbf = SnpArray("tmp.bed")
-@test size(tmpbf) == (5, 3)
-@test isfile("tmp.bed")
-@test isfile("tmp.fam")
-@test isfile("tmp.bim")
-rm(SnpArrays.datadir("mouse.bed.gz"), force=true)
-rm(SnpArrays.datadir("mouse.fam.gz"), force=true)
-rm(SnpArrays.datadir("mouse.bim.gz"), force=true)
-Sys.iswindows() || rm("tmp.bed", force=true)
-rm("tmp.bim", force=true)
-rm("tmp.fam", force=true)
-end
-
-@testset "SnpBitMatrix-vector multiplication" begin
-reltol = 5e-4
-for t in [Float32, Float64]
-    v1 = randn(t, size(EUR, 1))
-    v2 = randn(t, size(EUR, 2))
-    for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
-        @test norm(SnpBitMatrix{t}(EUR, model=model) * v2 -
-            convert(Matrix{t}, EUR, model=model) * v2) /
-            norm(convert(Matrix{t}, EUR, model=model) * v2) < reltol
-        @test norm(SnpBitMatrix{t}(EUR, center=true, model=model) * v2 -
-            convert(Matrix{t}, EUR, center=true, model=model) * v2) /
-            norm(convert(Matrix{t}, EUR, center=true, model=model) * v2) < reltol
-        @test norm(SnpBitMatrix{t}(EUR, scale=true, model=model) * v2 -
-            convert(Matrix{t}, EUR, scale=true, model=model) * v2) /
-            norm(convert(Matrix{t}, EUR, scale=true, model=model) * v2) < reltol
-        @test norm(SnpBitMatrix{t}(EUR, center=true, scale=true, model=model) * v2 -
-            convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) /
-            norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) < reltol
-        @test norm(transpose(SnpBitMatrix{t}(EUR, model=model)) * v1 -
-            transpose(convert(Matrix{t}, EUR, model=model)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, model=model)) * v1) < reltol
-        @test norm(transpose(SnpBitMatrix{t}(EUR, center=true, model=model)) * v1 -
-            transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) < reltol
-        @test norm(transpose(SnpBitMatrix{t}(EUR, scale=true, model=model)) * v1 -
-            transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) < reltol
-        @test norm(transpose(SnpBitMatrix{t}(EUR, center=true, scale=true, model=model)) * v1 -
-            transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) < reltol
-    end
-end
-end
-
-@testset "copyto bitmatrix" begin
-for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
-    mousebm = SnpBitMatrix{t}(mouse, model=ADDITIVE_MODEL, center=false, scale=false)
-    v = copyto!(zeros(1), @view(mousebm[702]))
-    @test v[1] == 0
-    mousebm = SnpBitMatrix{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true)
-    v = copyto!(zeros(1), @view(mousebm[702]))
-    @test v[1] ≈ -1 * mousebm.μ[1] * mousebm.σinv[1]
-    EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
-    EURtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
-    EURtest = copyto!(zeros(t, size(EUR)), EURbm)
-    @test all(EURtest .≈ EURtrue)
-    EURtest = copyto!(zeros(10, 5), @view(EURbm[11:20, 1:2:10]))
-    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
-    EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
-    EURtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
-    EURtest = copyto!(zeros(t, size(EUR)), EURbm)
-    @test all(EURtest .≈ EURtrue)
-    EURtest = copyto!(zeros(10, 5), @view(EURbm[11:20, 1:2:10]))
-    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
-end
-end
-
-@testset "convert (SnpBitMatrix)" begin
-for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
-    EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
-    xtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
-    xtest = convert(Matrix{t}, EURbm)
-    @test all(xtrue .≈ xtest)
-    EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
-    xtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
-    xtest = convert(Matrix{t}, EURbm)
-    @test all(xtrue .≈ xtest)
-    EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
-    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=true, scale=true)
-    xtest = convert(Matrix{t}, @view(EURbm[:, 2:2:10]))
-    @test all(xtrue .≈ xtest)
-    EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
-    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=false, scale=false)
-    xtest = convert(Matrix{t}, @view(EURbm[:, 2:2:10]))
-    @test all(xtrue .≈ xtest)
-end
-end
-
-@testset "copyto SnpLinAlg" begin
-for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
-    # imputing missing data
-    mousela = SnpLinAlg{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true, impute=true)
-    v = copyto!(zeros(1), @view(mousela[702])) # missing entry
-    @test isapprox(v[1], 1.113003134727478, atol=1e-6)
-    # not imputing missing data
-    mousela = SnpLinAlg{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true, impute=false)
-    v = copyto!(zeros(1), @view(mousela[702])) # missing entry
-    @test isnan(v[1])
-    # no missing data
-    EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
-    EURtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
-    EURtest = copyto!(zeros(t, size(EUR)), EURla)
-    @test all(EURtest .≈ EURtrue)
-    EURtest = copyto!(zeros(10, 5), @view(EURla[11:20, 1:2:10]))
-    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
-    EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
-    EURtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
-    EURtest = copyto!(zeros(t, size(EUR)), EURla)
-    @test all(EURtest .≈ EURtrue)
-    EURtest = copyto!(zeros(10, 5), @view(EURla[11:20, 1:2:10]))
-    @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
-end
-end
-
-@testset "convert (SnpLinAlg)" begin
-for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
-    EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
-    xtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
-    xtest = convert(Matrix{t}, EURla)
-    @test all(xtrue .≈ xtest)
-    EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
-    xtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
-    xtest = convert(Matrix{t}, EURla)
-    @test all(xtrue .≈ xtest)
-    EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
-    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=true, scale=true)
-    xtest = convert(Matrix{t}, @view(EURla[:, 2:2:10]))
-    @test all(xtrue .≈ xtest)
-    EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
-    xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=false, scale=false)
-    xtest = convert(Matrix{t}, @view(EURla[:, 2:2:10]))
-    @test all(xtrue .≈ xtest)
-end
-end
-
-@testset "SnpLinAlg-vector multiplication (zeroimpute)" begin
-reltol = 5e-4
-for t in [Float32, Float64]
-    v1 = randn(t, size(EUR, 1))
-    v2 = randn(t, size(EUR, 2))
-    for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
-        @test norm(SnpLinAlg{t}(EUR, model=model, impute=false) * v2 -
-            convert(Matrix{t}, EUR, model=model, impute=false) * v2) /
-            norm(convert(Matrix{t}, EUR, model=model, impute=false) * v2) < reltol
-        @test norm(SnpLinAlg{t}(EUR, center=true, model=model, impute=false) * v2 -
-            convert(Matrix{t}, EUR, center=true, model=model, impute=false) * v2) /
-            norm(convert(Matrix{t}, EUR, center=true, model=model, impute=false) * v2) < reltol
-        @test norm(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false) * v2 -
-            convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * v2) /
-            norm(convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * v2) < reltol
-        @test norm(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * v2 -
-            convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * v2) /
-            norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * v2) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, model=model, impute=false)) * v1 -
-            transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * v1) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=false)) * v1 -
-            transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * v1) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false)) * v1 -
-            transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * v1) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false)) * v1 -
-            transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * v1) < reltol
-    end
-end
-end
-
-@testset "SnpLinAlg-matrix multiplication zeroimpute" begin
-    for t in [Float32, Float64] # Y = A*X
-        X1 = randn(t, size(EUR, 1), 3)
-        X2 = randn(t, size(EUR, 2), 3)
-        for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
-            @test all(isapprox(SnpLinAlg{t}(EUR, model=model, impute=false) * X2,
-                convert(Matrix{t}, EUR, model=model, impute=false) * X2))
-            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, model=model, impute=false) * X2,
-                convert(Matrix{t}, EUR, center=true, model=model, impute=false) * X2))
-            @test all(isapprox(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false) * X2,
-                convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * X2))
-            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * X2,
-                convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * X2))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, model=model, impute=false)) * X1,
-                transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * X1))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=false)) * X1,
-                transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * X1))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false)) * X1,
-                transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * X1))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false)) * X1,
-                transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * X1))
-        end
-    end
-end
-
-@testset "SnpLinAlg-vector multiplication (Miter > 0)" begin
-    EUR11 = [EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR]
-    EUR11la = SnpLinAlg{Float64}(EUR11, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)
-    v = rand(size(EUR11la, 2))
-    vtest = EUR11la * v
-    vtrue = convert(Matrix{Float64}, EUR11, model=ADDITIVE_MODEL, impute=true, center=true, scale=true) * v
-    @test norm(vtest - vtrue) < 5e-4
-end
-
-if get(ENV,"JULIA_SNPARRAYS_TEST_CUDA","") == "true"
-    using CUDA, Adapt
-    @testset "lin. alg. cuda" begin
-    reltol = 5e-4
-    for t in [Float32, Float64]
-        v1 = randn(t, size(EUR, 1))
-        v2 = randn(t, size(EUR, 2))
-        v1_d = adapt(CuVector{t}, v1)
-        v2_d = adapt(CuVector{t}, v2)
-        for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
-            @test norm(collect(CuSnpArray{t}(EUR, model=model) * v2_d) -
-                convert(Matrix{t}, EUR, model=model) * v2) /
-                norm(convert(Matrix{t}, EUR, model=model) * v2) < reltol
-            @test norm(collect(CuSnpArray{t}(EUR, center=true, model=model) * v2_d) -
-                convert(Matrix{t}, EUR, center=true, model=model) * v2) /
-                norm(convert(Matrix{t}, EUR, center=true, model=model) * v2) < reltol
-            @test norm(collect(CuSnpArray{t}(EUR, scale=true, model=model) * v2_d) -
-                convert(Matrix{t}, EUR, scale=true, model=model) * v2) /
-                norm(convert(Matrix{t}, EUR, scale=true, model=model) * v2) < reltol
-            @test norm(collect(CuSnpArray{t}(EUR, center=true, scale=true, model=model) * v2_d) -
-                convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) /
-                norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) < reltol
-            @test norm(collect(transpose(CuSnpArray{t}(EUR, model=model)) * v1_d) -
-                transpose(convert(Matrix{t}, EUR, model=model)) * v1) /
-                norm(transpose(convert(Matrix{t}, EUR, model=model)) * v1) < reltol
-            @test norm(collect(transpose(CuSnpArray{t}(EUR, center=true, model=model)) * v1_d) -
-                transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) /
-                norm(transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) < reltol
-            @test norm(collect(transpose(CuSnpArray{t}(EUR, scale=true, model=model)) * v1_d) -
-                transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) /
-                norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) < reltol
-            @test norm(collect(transpose(CuSnpArray{t}(EUR, center=true, scale=true, model=model)) * v1_d) -
-                transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) /
-                norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) < reltol
-        end
-    end
-    end
-end
-
-@testset "SnpLinAlg-vector multiplication meanimpute" begin
-reltol = 5e-4
-for t in [Float32, Float64]
-    v1 = randn(t, size(EUR, 1))
-    v2 = randn(t, size(EUR, 2))
-    for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
-        @test norm(SnpLinAlg{t}(EUR, model=model, impute=true) * v2 -
-            convert(Matrix{t}, EUR, model=model, impute=true) * v2) /
-            norm(convert(Matrix{t}, EUR, model=model, impute=true) * v2) < reltol
-        @test norm(SnpLinAlg{t}(EUR, center=true, model=model, impute=true) * v2 -
-            convert(Matrix{t}, EUR, center=true, model=model, impute=true) * v2) /
-            norm(convert(Matrix{t}, EUR, center=true, model=model, impute=true) * v2) < reltol
-        @test norm(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true) * v2 -
-            convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * v2) /
-            norm(convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * v2) < reltol
-        @test norm(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true) * v2 -
-            convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * v2) /
-            norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * v2) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, model=model, impute=true)) * v1 -
-            transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * v1) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=true)) * v1 -
-            transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * v1) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true)) * v1 -
-            transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * v1) < reltol
-        @test norm(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true)) * v1 -
-            transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * v1) /
-            norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * v1) < reltol
-    end
-end
-end
-
-@testset "SnpLinAlg-matrix multiplication meanimpute" begin
-    for t in [Float32, Float64] # Y = A*X
-        X1 = randn(t, size(EUR, 1), 3)
-        X2 = randn(t, size(EUR, 2), 3)
-        for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
-            @test all(isapprox(SnpLinAlg{t}(EUR, model=model, impute=true) * X2,
-                convert(Matrix{t}, EUR, model=model, impute=true) * X2))
-            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, model=model, impute=true) * X2,
-                convert(Matrix{t}, EUR, center=true, model=model, impute=true) * X2))
-            @test all(isapprox(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true) * X2,
-                convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * X2))
-            @test all(isapprox(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * X2,
-                convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * X2))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, model=model, impute=true)) * X1,
-                transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * X1))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=true)) * X1,
-                transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * X1))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true)) * X1,
-                transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * X1))
-            @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true)) * X1,
-                transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * X1))
-        end
-    end
-end
-
-@testset "subarrays" begin
-@test all(@view(EUR[1:2:10, 1:2:10]) .==
-[[0x03 0x03 0x02 0x02 0x03];
-[0x03 0x03 0x03 0x03 0x03];
-[0x03 0x03 0x03 0x03 0x03];
-[0x02 0x03 0x02 0x00 0x02];
-[0x03 0x03 0x02 0x02 0x03]
-])
-@test all(convert(Matrix{Float64}, @view(EUR[1:2:10, 1:2:10])) .≈
-[[2.0 2.0 1.0 1.0 2.0];
-[2.0 2.0 2.0 2.0 2.0];
-[2.0 2.0 2.0 2.0 2.0];
-[1.0 2.0 1.0 0.0 1.0];
-[2.0 2.0 1.0 1.0 2.0]
-])
-
-EURsub = @view EUR[1:2:100, 1:2:100]
-EURsubbm = SnpBitMatrix{Float64}(EURsub, model=ADDITIVE_MODEL, center=true, scale=true) # BitMatrix from the SubArray
-EURsubfm = convert(Matrix{Float64}, EURsub, model=ADDITIVE_MODEL, center=true, scale=true) # Float64 Matrix from the SubArray
-
-v1 = randn(size(EURsub, 1))
-v2 = randn(size(EURsub, 2))
-
-@test isapprox(EURsubbm * v2, EURsubfm * v2)
-@test isapprox(EURsubbm' * v1, EURsubfm' * v1)
-end
-
-@testset "split-merge-readwrite" begin
-EUR_data = SnpData(SnpArrays.datadir("EUR_subset"))
-# a small subset of EUR_data
-SnpArrays.filter(EUR_data, collect(1:10), collect(1:20000); des="EUR_subset.small")
-EUR_small = SnpData("EUR_subset.small")
-
-# filter
-chr17 = SnpArrays.filter(SnpArrays.datadir("EUR_subset"); des="tmp.filter.chr17", f_snp = x -> String(x[:chromosome]) == "17")
-@test chr17.snps == 11041
-@test chr17.people == 379
-male = SnpArrays.filter(EUR_data; des="tmp.filter.male", f_person = x -> String(x[:sex]) == "1")
-@test male.snps == 54051
-@test male.people == 178
-chr17_male = SnpArrays.filter(EUR_data; des="tmp.filter.chr17.male", f_snp = x -> String(x[:chromosome]) == "17", f_person = x -> String(x[:sex]) == "1")
-@test chr17_male.snps == 11041
-@test chr17_male.people == 178
-@test size(chr17_male.snparray) == (178, 11041)
-# cleanup
-for ft in ["bim", "fam", "bed"]
-    ft == "bed" && Sys.iswindows() && continue
-    rm("tmp.filter.chr17." * ft, force=true)
-    rm("tmp.filter.male." * ft, force=true)
-    rm("tmp.filter.chr17.male." * ft, force=true)
-end
-
-# split
-splitted = SnpArrays.split_plink("EUR_subset.small"; prefix="tmp.split.chr.")
-piece = splitted["17"]
-@test piece.people == 10
-@test piece.snps == 11041
-@test size(piece.person_info) == (10, 6)
-@test size(piece.snp_info) == (11041, 6)
-@test size(piece.snparray) == (10, 11041)
-@test size(piece.snparray.columncounts) == (4, 11041)
-@test piece.snparray.m == 10
-
-splitted_bysex = SnpArrays.split_plink(EUR_data, :sex; prefix="tmp.split.sex.")
-@test splitted_bysex["1"].people == 178
-@test splitted_bysex["2"].people == 201
-
-# merge
-#@time merged = SnpArrays.merge_plink("tmp.merged", splitted) # write_plink is included here
-#@test EUR_data.people == merged.people
-#@test EUR_data.snps == merged.snps
-#@test EUR_data.person_info == merged.person_info
-#@test EUR_data.snp_info == merged.snp_info # note: the ordering of merged data might be different on other dataset b/c sorted order of chromosomes
-
-#output = SnpData("tmp.merged")
-#@test EUR_data.people == output.people
-#@test EUR_data.snps == output.snps
-#@test EUR_data.person_info == output.person_info
-#@test EUR_data.snp_info == output.snp_info
-
-# cleanup
-#rm("tmp.merged.bim", force=true)
-#rm("tmp.merged.fam", force=true)
-#Sys.iswindows() || rm("tmp.merged.bed", force=true)
-
-# merge from splitted files
-@time merged_from_splitted_files = merge_plink("tmp.split.chr"; des = "tmp2.merged")
-@test EUR_small.people == merged_from_splitted_files.people
-@test EUR_small.snps == merged_from_splitted_files.snps
-@test EUR_small.person_info == merged_from_splitted_files.person_info
-@test EUR_small.snp_info == merged_from_splitted_files.snp_info
-
-# cleanup
-rm("tmp2.merged.bim", force=true)
-rm("tmp2.merged.fam", force=true)
-Sys.iswindows() || rm("tmp2.merged.bed", force=true)
-
-rm("EUR_subset.small.bim", force=true)
-rm("EUR_subset.small.fam", force=true)
-Sys.iswindows() || rm("EUR_subset.small.bed", force=true)
-
-for k in keys(splitted)
-    for ft in ["bim", "fam", "bed"]
-        ft == "bed" && Sys.iswindows() && continue
-        rm("tmp.split.chr.$(k)." * ft, force=true)
-    end
-end
-for k in keys(splitted_bysex)
-    for ft in ["bim", "fam", "bed"]
-        ft == "bed" && Sys.iswindows() && continue
-        rm("tmp.split.sex.$(k)." * ft, force=true)
-    end
-end
-end
-
-@testset "(de)compress" begin
-for format in SnpArrays.ALLOWED_FORMAT
-    # compress mouse Plink files
-    compress_plink(SnpArrays.datadir("mouse"), format)
-    @test isfile(SnpArrays.datadir("mouse.bed." * format))
-    @test isfile(SnpArrays.datadir("mouse.fam." * format))
-    @test isfile(SnpArrays.datadir("mouse.bim." * format))
-    # read in compressed Plink files
-    mouse_zip = SnpArray(SnpArrays.datadir("mouse.bed." * format))
-    fill!(mouse.rowcounts, 0)
-    fill!(mouse.columncounts, 0)
-    @time begin
-        SnpArrays._counts(mouse, 1)
-        SnpArrays._counts(mouse, 2)
-    end
-    fill!(mouse_zip.rowcounts, 0)
-    fill!(mouse_zip.columncounts, 0)
-    @time begin
-        SnpArrays._counts(mouse_zip, 1)
-        SnpArrays._counts(mouse_zip, 2)
-    end
-    @test mouse.m == mouse_zip.m
-    @test norm(mouse.rowcounts - mouse_zip.rowcounts) < 1e-8
-    @test norm(mouse.columncounts - mouse_zip.columncounts) < 1e-8
-    @test norm(mouse.data - mouse_zip.data) < 1e-8
-    # Decompress Plink files
-    decompress_plink(SnpArrays.datadir("mouse"), format, SnpArrays.datadir("mouse2"))
-    @test stat(SnpArrays.datadir("mouse") * ".bed").size ==
-    stat(SnpArrays.datadir("mouse2") * ".bed").size
-    @test stat(SnpArrays.datadir("mouse") * ".fam").size ==
-    stat(SnpArrays.datadir("mouse2") * ".fam").size
-    @test stat(SnpArrays.datadir("mouse") * ".bim").size ==
-    stat(SnpArrays.datadir("mouse2") * ".bim").size
-    # clean up
-    rm(SnpArrays.datadir("mouse.bed." * format), force=true)
-    rm(SnpArrays.datadir("mouse.fam." * format), force=true)
-    rm(SnpArrays.datadir("mouse.bim." * format), force=true)
-    rm(SnpArrays.datadir("mouse2.bed"), force=true)
-    rm(SnpArrays.datadir("mouse2.fam"), force=true)
-    rm(SnpArrays.datadir("mouse2.bim"), force=true)
-end
-@test_throws ArgumentError SnpArray(SnpArrays.datadir("mouse.bed.zip"))
-end
-
-@testset "indexin" begin
-    b = ['1', '2', '3', '4']
-    for a in [['4', '3', '1', '2'], ['3', '2', '6'], ['b', 'c', 'd', 'c']]
-        aind, bmask = SnpArrays.indexin_general(a, b)
-        @test all(a[aind] .== b[bmask])
-    end
-end
-
-@testset "concat" begin
-s = SnpArrays.filter(SnpArrays.datadir("mouse"), 1:2, 1:3; des="mouse.tmp.filtered")
-sd = SnpData("mouse.tmp.filtered")
-sd_hcat = hcat(sd, sd, sd; des="mouse.hcat")
-sd_vcat = vcat(sd, sd, sd; des="mouse.vcat")
-sd_hvcat = hvcat((2, 2), sd, sd, sd, sd; des="mouse.hvcat")
-
-ref = [[0x02 0x02 0x02];
-[0x02 0x02 0x03]]
-
-@test all(sd_hcat.snparray .== [ref ref ref])
-@test all(sd_vcat.snparray .== [ref; ref; ref])
-@test all(sd_hvcat.snparray .== [ref ref; ref ref])
-
-@test sd_hcat.people == 2
-@test sd_hcat.snps == 9
-
-@test sd_vcat.people == 6
-@test sd_vcat.snps == 3
-
-@test sd_hvcat.people == 4
-@test sd_hvcat.snps == 6
-
-for ft in ["bim", "fam", "bed"]
-    ft == "bed" && Sys.iswindows() && continue
-    rm("mouse.tmp.filtered." * ft, force=true)
-    rm("mouse.hcat." * ft, force=true)
-    rm("mouse.vcat." * ft, force=true)
-    rm("mouse.hvcat." * ft, force=true)
-end
-end
-
-@testset "reorder" begin
-mouse_prefix = SnpArrays.datadir("mouse")
-run(`cp $(mouse_prefix * ".bed") mouse_testreorder.bed`)
-run(`cp $(mouse_prefix * ".bim") mouse_testreorder.bim`)
-run(`cp $(mouse_prefix * ".fam") mouse_testreorder.fam`)
-
-mouse_data = SnpData(mouse_prefix)
-mouse_toreorder = SnpData("mouse_testreorder", "r+")
-m, n = size(mouse_toreorder.snparray)
-using Random
-ind = randperm(m);
-SnpArrays.reorder!(mouse_toreorder, ind)
-
-# reread the file
-mouse_toreorder_read = SnpData("mouse_testreorder"; famnm="mouse_testreorder.reordered.fam")
-@test all(mouse_data.snparray[ind, :] .== mouse_toreorder_read.snparray[:, :])
-@test all(mouse_data.person_info[ind, 1] .== mouse_toreorder_read.person_info[:, 1])
-
-# cleanup
-rm("mouse_testreorder.bed", force=true)
-rm("mouse_testreorder.bim", force=true)
-rm("mouse_testreorder.fam", force=true)
-rm("mouse_testreorder.reordered.fam", force=true)
-end
-
-@testset "vcf2plink" begin
-# Download an example VCF file
-isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
-    joinpath(pwd(), "test.08Jun17.d8b.vcf.gz"));
-s = vcf2plink("test.08Jun17.d8b.vcf.gz", "test.08Jun17.d8b")
-@test size(s, 1) == 191
-@test size(s, 2) == 1354
-@test s[1, 1] == 0x00
-@test s[1, 5] == 0x02
-@test s[38, 5] == 0x03
-rm("test.08Jun17.d8b.vcf.gz", force=true)
-rm("test.08Jun17.d8b.bed", force=true)
-rm("test.08Jun17.d8b.bim", force=true)
-rm("test.08Jun17.d8b.fam", force=true)
-end
-
-@testset "kinship_pruning" begin
-g = grm(mouse)
-@test count(kinship_pruning(g; method=:gcta)) == 68
-@test count(kinship_pruning(g; method=:top_down)) == 125
-@test count(kinship_pruning(g; method=:bottom_up)) == 132
-@test count(kinship_pruning(g; method=:plink)) == 126
-end
-
-@testset "stackedsnparray" begin
-mouse2 = StackedSnpArray([mouse, mouse])
-@test size(mouse2) == (1940, 20300)
-@test size(mouse2, 3) == 1
-@test_throws ErrorException size(mouse2, 0)
-@test length(mouse2) == 1940 * 20300
-@test eltype(mouse2) == UInt8
-@test mouse2[777, 16384] == mouse[777, 16384 - 10150]
-end
+# const mouse = SnpArray(SnpArrays.datadir("mouse.bed")) # has missing genotypes
+
+# @testset "size" begin
+# @test size(EUR) == (379, 54051)
+# @test size(EUR, 1) == 379
+# @test size(EUR, 2) == 54051
+# @test size(mouse) == (1940, 10150)
+# @test size(mouse, 3) == 1
+# @test_throws ErrorException size(mouse, 0)
+# @test length(mouse) == 1940 * 10150
+# @test eltype(mouse) == UInt8
+# end
+
+# @testset "counts" begin
+# cc = counts(EUR, dims=1)
+# @test counts(EUR) == sum(cc, dims=2)
+# @test size(cc) == (4, size(EUR, 2))
+# @test view(cc, :, 1) == [2, 0, 70, 307]
+# @test all(sum(cc, dims = 1) .== size(EUR, 1))
+# @test all(iszero, view(cc, 2, :))
+# rc = counts(EUR, dims=2)
+# @test size(rc) == (4, size(EUR, 1))
+# @test view(rc, :, 1) == [2997, 0, 13143, 37911]
+# @test all(sum(rc, dims = 1) .== size(EUR, 2))
+# @test all(iszero, view(rc, 2, :))
+# @test_throws ArgumentError counts(EUR, dims=3)
+# end
+
+# @testset "means" begin
+# # ADDITIVE_MODEL
+# cmns = mean(EUR, dims = 1)
+# @test size(cmns) == (1, size(EUR, 2))
+# @test cmns[1] ≈ 1.804749340369393
+# @test mean(cmns) ≈ mean(EUR)
+# @test minimum(cmns) ≈ 1.0
+# @test maximum(cmns) ≈ 1.9788918205804749
+# rmns = mean(EUR, dims = 2)
+# @test size(rmns) == (size(EUR, 1), 1)
+# @test rmns[1] ≈ 1.6459454959205195
+# @test mean(rmns) ≈ mean(EUR)
+# rmnextrma = extrema(rmns)
+# @test rmnextrma[1] ≈ 1.637749532848606
+# @test rmnextrma[2] ≈ 1.6627259440158368
+# @test_throws ArgumentError mean(EUR, dims=3)
+# # DOMINANT_MODEL
+# cmns = mean(EUR, dims = 1, model = DOMINANT_MODEL)
+# @test size(cmns) == (1, size(EUR, 2))
+# @test cmns[1] ≈ 0.9947229551451188
+# @test mean(cmns) ≈ mean(EUR, model = DOMINANT_MODEL)
+# @test minimum(cmns) ≈ 0.5989445910290238
+# @test maximum(cmns) ≈ 1.0
+# rmns = mean(EUR, dims = 2, model = DOMINANT_MODEL)
+# @test size(rmns) == (size(EUR, 1), 1)
+# @test rmns[1] ≈ 0.9445523672087472
+# @test mean(rmns) ≈ mean(EUR, model = DOMINANT_MODEL)
+# rmnextrma = extrema(rmns)
+# @test rmnextrma[1] ≈ 0.9325266877578583
+# @test rmnextrma[2] ≈ 0.9530073449149877
+# @test_throws ArgumentError mean(EUR, dims=3, model=DOMINANT_MODEL)
+# # RECESSIVE_MODEL
+# cmns = mean(EUR, dims = 1, model = RECESSIVE_MODEL)
+# @test size(cmns) == (1, size(EUR, 2))
+# @test cmns[1] ≈ 0.8100263852242744
+# @test mean(cmns) ≈ mean(EUR, model = RECESSIVE_MODEL)
+# @test minimum(cmns) ≈ 0.0
+# @test maximum(cmns) ≈ 0.9868073878627969
+# rmns = mean(EUR, dims = 2, model = RECESSIVE_MODEL)
+# @test size(rmns) == (size(EUR, 1), 1)
+# @test rmns[1] ≈ 0.7013931287117722
+# @test mean(rmns) ≈ mean(EUR, model = RECESSIVE_MODEL)
+# rmnextrma = extrema(rmns)
+# @test rmnextrma[1] ≈ 0.6932341677304767
+# @test rmnextrma[2] ≈ 0.7191171301178516
+# @test_throws ArgumentError mean(EUR, dims=3, model=RECESSIVE_MODEL)
+# end
+
+# @testset "var/std" begin
+# cvars = var(EUR, dims = 1)
+# @test size(cvars) == (1, size(EUR, 2))
+# @test cvars[1] ≈ 0.16812553224162724
+# @test iszero(minimum(cvars))
+# @test maximum(cvars) ≈ 0.8705727966941688
+# @test std(EUR, dims = 1) ≈ sqrt.(cvars)
+# rvars = var(EUR, dims = 2)
+# @test size(rvars) == (size(EUR, 1), 1)
+# @test rvars[1] ≈ 0.33960146078503206
+# @test minimum(rvars) ≈ 0.31972707244268267
+# @test maximum(rvars) ≈ 0.365901927927013
+# @test std(EUR, dims = 2) ≈ sqrt.(rvars)
+# end
+
+# @testset "missingpos" begin
+# mp = missingpos(mouse)
+# cc = counts(mouse, dims=1)
+# @test isa(mp, SparseMatrixCSC)
+# @test sum(mp, dims = 1) == view(cc, 2:2, :)
+# @test sum(mp, dims = 2) == view(counts(mouse, dims=2), 2:2, :)'
+# end
+
+# @testset "create bed" begin
+# tmpbf = SnpArray("tmp.bed", 5, 3)
+# @test isfile("tmp.bed")
+# @test all(tmpbf .== 0x00)
+# fill!(tmpbf, 0x02)
+# tmpbf2 = SnpArray("tmp.bed", 5)
+# @test all(tmpbf2 .== 0x02)
+# Sys.iswindows() || rm("tmp.bed", force=true)
+# tmpbf2 = SnpArray("tmp2.bed", SnpArray(undef, 5, 3))
+# fill!(tmpbf2, 0x01)
+# @test all(tmpbf2 .== 0x01)
+# Sys.iswindows() || rm("tmp2.bed", force=true)
+# end
+
+# @testset "convert" begin
+# # check convert coding
+# tmpbf = SnpArray(undef, 4, 1)
+# tmpbf[1] = 0x00
+# tmpbf[2] = 0x01
+# tmpbf[3] = 0x02
+# tmpbf[4] = 0x03
+# # additive model
+# v = convert(Matrix{Float64}, tmpbf, model=ADDITIVE_MODEL)
+# @test all([v[1], v[3], v[4]] .== [0.0, 1.0, 2.0])
+# @test isnan(v[2])
+# # dominant model
+# v = convert(Matrix{Float64}, tmpbf, model=DOMINANT_MODEL)
+# @test all([v[1], v[3], v[4]] .== [0.0, 1.0, 1.0])
+# @test isnan(v[2])
+# # recessive model
+# v = convert(Matrix{Float64}, tmpbf, model=RECESSIVE_MODEL)
+# @test all([v[1], v[3], v[4]] .== [0.0, 0.0, 1.0])
+# @test isnan(v[2])
+# # mouse test data
+# cmns = mean(mouse, dims = 1)
+# cvrs = var(mouse, dims = 1)
+# cstd = std(mouse, dims = 1)
+# v = convert(Vector{Float64}, @view(mouse[:, 1]))
+# @test isnan(mean(v))
+# @test mean(filter(!isnan, v)) ≈ cmns[1]
+# @test isnan(var(v))
+# @test var(filter(!isnan, v)) ≈ cvrs[1]
+# @test isnan(std(v))
+# @test std(filter(!isnan, v)) ≈ cstd[1]
+# end
+
+# @testset "convert (ADMIXTURE)" begin
+# # convert A2 allele frequency to dosage
+# @test convert(Float64, 0.5, ADDITIVE_MODEL)  == 1.0
+# @test convert(Float64, 0.5, DOMINANT_MODEL)  == 0.75
+# @test convert(Float64, 0.5, RECESSIVE_MODEL) == 0.25
+# end
+
+# @testset "maf" begin
+# cc = counts(mouse, dims=1)
+# @test all(0 .≤ maf(mouse) .≤ 0.5)
+# @test all(minorallele(mouse) .== (cc[1, :] .> cc[4, :]))
+# cc = counts(EUR, dims=1)
+# @test all(0 .≤ maf(EUR) .≤ 0.5)
+# @test all(minorallele(EUR) .== (cc[1, :] .> cc[4, :]))
+# end
+
+# @testset "grm" begin
+# Φgrm = grm(EUR, method=:GRM)
+# @test size(Φgrm) == (size(EUR, 1), size(EUR, 1))
+# @test issymmetric(Φgrm)
+# @test eigmin(Φgrm) > -1e-8
+# Φmom = grm(EUR, method=:MoM)
+# @test size(Φmom) == (size(EUR, 1), size(EUR, 1))
+# @test issymmetric(Φmom)
+# Φrbs = grm(EUR, method=:Robust)
+# @test size(Φrbs) == (size(EUR, 1), size(EUR, 1))
+# @test issymmetric(Φrbs)
+# @test eigmin(Φrbs) > -1e-8
+# m, n, K = size(EUR, 1), size(EUR, 2), 3
+# P = rand(n, K)
+# Q = rand(m, K)
+# Q ./= sum(Q, dims=2)
+# Φreap = grm_admixture(EUR, transpose(P), transpose(Q))
+# @test size(Φreap) == (size(EUR, 1), size(EUR, 1))
+# @test issymmetric(Φreap)
+# @test eigmin(Φrbs) > -1e-8
+# end
+
+# @testset "filter" begin
+# rowmask, colmask =  SnpArrays.filter(mouse, min_success_rate_per_row=0.99,
+#     min_success_rate_per_col=0.99)
+# SnpArrays.filter(SnpArrays.datadir("mouse"), rowmask, colmask; des="tmp")
+# tmpbf = SnpArray("tmp.bed")
+# @test size(tmpbf) == (1912, 9997)
+# @test all(missingrate(tmpbf, 1) .≤ 0.01)
+# @test all(missingrate(tmpbf, 2) .≤ 0.01)
+# Sys.iswindows() || rm("tmp.bed", force=true)
+# rm("tmp.bim", force=true)
+# rm("tmp.fam", force=true)
+# rowmask, colmask =  SnpArrays.filter(mouse, min_success_rate_per_row=0.99,
+#     min_success_rate_per_col=0.99, min_maf=0.01, min_hwe_pval=1e-8)
+# @test (count(rowmask), count(colmask)) == (1911, 9510)
+# compress_plink(SnpArrays.datadir("mouse"), "gz")
+# SnpArrays.filter(SnpArrays.datadir("mouse.bed.gz"), SnpArrays.datadir("mouse.bim.gz"),
+# SnpArrays.datadir("mouse.fam.gz"), 1:5, 1:3; des="tmp")
+# tmpbf = SnpArray("tmp.bed")
+# @test size(tmpbf) == (5, 3)
+# @test isfile("tmp.bed")
+# @test isfile("tmp.fam")
+# @test isfile("tmp.bim")
+# rm(SnpArrays.datadir("mouse.bed.gz"), force=true)
+# rm(SnpArrays.datadir("mouse.fam.gz"), force=true)
+# rm(SnpArrays.datadir("mouse.bim.gz"), force=true)
+# Sys.iswindows() || rm("tmp.bed", force=true)
+# rm("tmp.bim", force=true)
+# rm("tmp.fam", force=true)
+# end
+
+# @testset "SnpBitMatrix-vector multiplication" begin
+# reltol = 5e-4
+# for t in [Float32, Float64]
+#     v1 = randn(t, size(EUR, 1))
+#     v2 = randn(t, size(EUR, 2))
+#     for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
+#         @test norm(SnpBitMatrix{t}(EUR, model=model) * v2 -
+#             convert(Matrix{t}, EUR, model=model) * v2) /
+#             norm(convert(Matrix{t}, EUR, model=model) * v2) < reltol
+#         @test norm(SnpBitMatrix{t}(EUR, center=true, model=model) * v2 -
+#             convert(Matrix{t}, EUR, center=true, model=model) * v2) /
+#             norm(convert(Matrix{t}, EUR, center=true, model=model) * v2) < reltol
+#         @test norm(SnpBitMatrix{t}(EUR, scale=true, model=model) * v2 -
+#             convert(Matrix{t}, EUR, scale=true, model=model) * v2) /
+#             norm(convert(Matrix{t}, EUR, scale=true, model=model) * v2) < reltol
+#         @test norm(SnpBitMatrix{t}(EUR, center=true, scale=true, model=model) * v2 -
+#             convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) /
+#             norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) < reltol
+#         @test norm(transpose(SnpBitMatrix{t}(EUR, model=model)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, model=model)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, model=model)) * v1) < reltol
+#         @test norm(transpose(SnpBitMatrix{t}(EUR, center=true, model=model)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) < reltol
+#         @test norm(transpose(SnpBitMatrix{t}(EUR, scale=true, model=model)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) < reltol
+#         @test norm(transpose(SnpBitMatrix{t}(EUR, center=true, scale=true, model=model)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) < reltol
+#     end
+# end
+# end
+
+# @testset "copyto bitmatrix" begin
+# for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
+#     mousebm = SnpBitMatrix{t}(mouse, model=ADDITIVE_MODEL, center=false, scale=false)
+#     v = copyto!(zeros(1), @view(mousebm[702]))
+#     @test v[1] == 0
+#     mousebm = SnpBitMatrix{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true)
+#     v = copyto!(zeros(1), @view(mousebm[702]))
+#     @test v[1] ≈ -1 * mousebm.μ[1] * mousebm.σinv[1]
+#     EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
+#     EURtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
+#     EURtest = copyto!(zeros(t, size(EUR)), EURbm)
+#     @test all(EURtest .≈ EURtrue)
+#     EURtest = copyto!(zeros(10, 5), @view(EURbm[11:20, 1:2:10]))
+#     @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
+#     EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
+#     EURtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
+#     EURtest = copyto!(zeros(t, size(EUR)), EURbm)
+#     @test all(EURtest .≈ EURtrue)
+#     EURtest = copyto!(zeros(10, 5), @view(EURbm[11:20, 1:2:10]))
+#     @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
+# end
+# end
+
+# @testset "convert (SnpBitMatrix)" begin
+# for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
+#     EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
+#     xtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
+#     xtest = convert(Matrix{t}, EURbm)
+#     @test all(xtrue .≈ xtest)
+#     EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
+#     xtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
+#     xtest = convert(Matrix{t}, EURbm)
+#     @test all(xtrue .≈ xtest)
+#     EURbm = SnpBitMatrix{t}(EUR, model=model, center=true, scale=true)
+#     xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=true, scale=true)
+#     xtest = convert(Matrix{t}, @view(EURbm[:, 2:2:10]))
+#     @test all(xtrue .≈ xtest)
+#     EURbm = SnpBitMatrix{t}(EUR, model=model, center=false, scale=false)
+#     xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=false, scale=false)
+#     xtest = convert(Matrix{t}, @view(EURbm[:, 2:2:10]))
+#     @test all(xtrue .≈ xtest)
+# end
+# end
+
+# @testset "copyto SnpLinAlg" begin
+# for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
+#     # imputing missing data
+#     mousela = SnpLinAlg{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true, impute=true)
+#     v = copyto!(zeros(1), @view(mousela[702])) # missing entry
+#     @test isapprox(v[1], 1.113003134727478, atol=1e-6)
+#     # not imputing missing data
+#     mousela = SnpLinAlg{t}(mouse, model=ADDITIVE_MODEL, center=true, scale=true, impute=false)
+#     v = copyto!(zeros(1), @view(mousela[702])) # missing entry
+#     @test isnan(v[1])
+#     # no missing data
+#     EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
+#     EURtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
+#     EURtest = copyto!(zeros(t, size(EUR)), EURla)
+#     @test all(EURtest .≈ EURtrue)
+#     EURtest = copyto!(zeros(10, 5), @view(EURla[11:20, 1:2:10]))
+#     @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
+#     EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
+#     EURtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
+#     EURtest = copyto!(zeros(t, size(EUR)), EURla)
+#     @test all(EURtest .≈ EURtrue)
+#     EURtest = copyto!(zeros(10, 5), @view(EURla[11:20, 1:2:10]))
+#     @test all(EURtest .≈ EURtrue[11:20, 1:2:10])
+# end
+# end
+
+# @testset "convert (SnpLinAlg)" begin
+# for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL], t in [Float32, Float64]
+#     EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
+#     xtrue = convert(Matrix{t}, EUR, model=model, center=true, scale=true)
+#     xtest = convert(Matrix{t}, EURla)
+#     @test all(xtrue .≈ xtest)
+#     EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
+#     xtrue = convert(Matrix{t}, EUR, model=model, center=false, scale=false)
+#     xtest = convert(Matrix{t}, EURla)
+#     @test all(xtrue .≈ xtest)
+#     EURla = SnpLinAlg{t}(EUR, model=model, center=true, scale=true)
+#     xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=true, scale=true)
+#     xtest = convert(Matrix{t}, @view(EURla[:, 2:2:10]))
+#     @test all(xtrue .≈ xtest)
+#     EURla = SnpLinAlg{t}(EUR, model=model, center=false, scale=false)
+#     xtrue = convert(Matrix{t}, @view(EUR[:, 2:2:10]), model=model, center=false, scale=false)
+#     xtest = convert(Matrix{t}, @view(EURla[:, 2:2:10]))
+#     @test all(xtrue .≈ xtest)
+# end
+# end
+
+# @testset "SnpLinAlg-vector multiplication (zeroimpute)" begin
+# reltol = 5e-4
+# for t in [Float32, Float64]
+#     v1 = randn(t, size(EUR, 1))
+#     v2 = randn(t, size(EUR, 2))
+#     for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
+#         @test norm(SnpLinAlg{t}(EUR, model=model, impute=false) * v2 -
+#             convert(Matrix{t}, EUR, model=model, impute=false) * v2) /
+#             norm(convert(Matrix{t}, EUR, model=model, impute=false) * v2) < reltol
+#         @test norm(SnpLinAlg{t}(EUR, center=true, model=model, impute=false) * v2 -
+#             convert(Matrix{t}, EUR, center=true, model=model, impute=false) * v2) /
+#             norm(convert(Matrix{t}, EUR, center=true, model=model, impute=false) * v2) < reltol
+#         @test norm(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false) * v2 -
+#             convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * v2) /
+#             norm(convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * v2) < reltol
+#         @test norm(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * v2 -
+#             convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * v2) /
+#             norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * v2) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, model=model, impute=false)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * v1) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=false)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * v1) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * v1) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * v1) < reltol
+#     end
+# end
+# end
+
+# @testset "SnpLinAlg-matrix multiplication zeroimpute" begin
+#     for t in [Float32, Float64] # Y = A*X
+#         X1 = randn(t, size(EUR, 1), 3)
+#         X2 = randn(t, size(EUR, 2), 3)
+#         for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
+#             @test all(isapprox(SnpLinAlg{t}(EUR, model=model, impute=false) * X2,
+#                 convert(Matrix{t}, EUR, model=model, impute=false) * X2))
+#             @test all(isapprox(SnpLinAlg{t}(EUR, center=true, model=model, impute=false) * X2,
+#                 convert(Matrix{t}, EUR, center=true, model=model, impute=false) * X2))
+#             @test all(isapprox(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false) * X2,
+#                 convert(Matrix{t}, EUR, scale=true, model=model, impute=false) * X2))
+#             @test all(isapprox(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * X2,
+#                 convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false) * X2))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, model=model, impute=false)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, model=model, impute=false)) * X1))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=false)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=false)) * X1))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=false)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=false)) * X1))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=false)) * X1))
+#         end
+#     end
+# end
+
+# @testset "SnpLinAlg-vector multiplication (Miter > 0)" begin
+#     EUR11 = [EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR;EUR]
+#     EUR11la = SnpLinAlg{Float64}(EUR11, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)
+#     v = rand(size(EUR11la, 2))
+#     vtest = EUR11la * v
+#     vtrue = convert(Matrix{Float64}, EUR11, model=ADDITIVE_MODEL, impute=true, center=true, scale=true) * v
+#     @test norm(vtest - vtrue) < 5e-4
+# end
+
+# if get(ENV,"JULIA_SNPARRAYS_TEST_CUDA","") == "true"
+#     using CUDA, Adapt
+#     @testset "lin. alg. cuda" begin
+#     reltol = 5e-4
+#     for t in [Float32, Float64]
+#         v1 = randn(t, size(EUR, 1))
+#         v2 = randn(t, size(EUR, 2))
+#         v1_d = adapt(CuVector{t}, v1)
+#         v2_d = adapt(CuVector{t}, v2)
+#         for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
+#             @test norm(collect(CuSnpArray{t}(EUR, model=model) * v2_d) -
+#                 convert(Matrix{t}, EUR, model=model) * v2) /
+#                 norm(convert(Matrix{t}, EUR, model=model) * v2) < reltol
+#             @test norm(collect(CuSnpArray{t}(EUR, center=true, model=model) * v2_d) -
+#                 convert(Matrix{t}, EUR, center=true, model=model) * v2) /
+#                 norm(convert(Matrix{t}, EUR, center=true, model=model) * v2) < reltol
+#             @test norm(collect(CuSnpArray{t}(EUR, scale=true, model=model) * v2_d) -
+#                 convert(Matrix{t}, EUR, scale=true, model=model) * v2) /
+#                 norm(convert(Matrix{t}, EUR, scale=true, model=model) * v2) < reltol
+#             @test norm(collect(CuSnpArray{t}(EUR, center=true, scale=true, model=model) * v2_d) -
+#                 convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) /
+#                 norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model) * v2) < reltol
+#             @test norm(collect(transpose(CuSnpArray{t}(EUR, model=model)) * v1_d) -
+#                 transpose(convert(Matrix{t}, EUR, model=model)) * v1) /
+#                 norm(transpose(convert(Matrix{t}, EUR, model=model)) * v1) < reltol
+#             @test norm(collect(transpose(CuSnpArray{t}(EUR, center=true, model=model)) * v1_d) -
+#                 transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) /
+#                 norm(transpose(convert(Matrix{t}, EUR, center=true, model=model)) * v1) < reltol
+#             @test norm(collect(transpose(CuSnpArray{t}(EUR, scale=true, model=model)) * v1_d) -
+#                 transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) /
+#                 norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model)) * v1) < reltol
+#             @test norm(collect(transpose(CuSnpArray{t}(EUR, center=true, scale=true, model=model)) * v1_d) -
+#                 transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) /
+#                 norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model)) * v1) < reltol
+#         end
+#     end
+#     end
+# end
+
+# @testset "SnpLinAlg-vector multiplication meanimpute" begin
+# reltol = 5e-4
+# for t in [Float32, Float64]
+#     v1 = randn(t, size(EUR, 1))
+#     v2 = randn(t, size(EUR, 2))
+#     for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
+#         @test norm(SnpLinAlg{t}(EUR, model=model, impute=true) * v2 -
+#             convert(Matrix{t}, EUR, model=model, impute=true) * v2) /
+#             norm(convert(Matrix{t}, EUR, model=model, impute=true) * v2) < reltol
+#         @test norm(SnpLinAlg{t}(EUR, center=true, model=model, impute=true) * v2 -
+#             convert(Matrix{t}, EUR, center=true, model=model, impute=true) * v2) /
+#             norm(convert(Matrix{t}, EUR, center=true, model=model, impute=true) * v2) < reltol
+#         @test norm(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true) * v2 -
+#             convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * v2) /
+#             norm(convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * v2) < reltol
+#         @test norm(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true) * v2 -
+#             convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * v2) /
+#             norm(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * v2) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, model=model, impute=true)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * v1) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=true)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * v1) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * v1) < reltol
+#         @test norm(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true)) * v1 -
+#             transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * v1) /
+#             norm(transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * v1) < reltol
+#     end
+# end
+# end
+
+# @testset "SnpLinAlg-matrix multiplication meanimpute" begin
+#     for t in [Float32, Float64] # Y = A*X
+#         X1 = randn(t, size(EUR, 1), 3)
+#         X2 = randn(t, size(EUR, 2), 3)
+#         for model in [ADDITIVE_MODEL, DOMINANT_MODEL, RECESSIVE_MODEL]
+#             @test all(isapprox(SnpLinAlg{t}(EUR, model=model, impute=true) * X2,
+#                 convert(Matrix{t}, EUR, model=model, impute=true) * X2))
+#             @test all(isapprox(SnpLinAlg{t}(EUR, center=true, model=model, impute=true) * X2,
+#                 convert(Matrix{t}, EUR, center=true, model=model, impute=true) * X2))
+#             @test all(isapprox(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true) * X2,
+#                 convert(Matrix{t}, EUR, scale=true, model=model, impute=true) * X2))
+#             @test all(isapprox(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=false) * X2,
+#                 convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true) * X2))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, model=model, impute=true)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, model=model, impute=true)) * X1))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, model=model, impute=true)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, center=true, model=model, impute=true)) * X1))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, scale=true, model=model, impute=true)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, scale=true, model=model, impute=true)) * X1))
+#             @test all(isapprox(transpose(SnpLinAlg{t}(EUR, center=true, scale=true, model=model, impute=true)) * X1,
+#                 transpose(convert(Matrix{t}, EUR, center=true, scale=true, model=model, impute=true)) * X1))
+#         end
+#     end
+# end
+
+# @testset "subarrays" begin
+# @test all(@view(EUR[1:2:10, 1:2:10]) .==
+# [[0x03 0x03 0x02 0x02 0x03];
+# [0x03 0x03 0x03 0x03 0x03];
+# [0x03 0x03 0x03 0x03 0x03];
+# [0x02 0x03 0x02 0x00 0x02];
+# [0x03 0x03 0x02 0x02 0x03]
+# ])
+# @test all(convert(Matrix{Float64}, @view(EUR[1:2:10, 1:2:10])) .≈
+# [[2.0 2.0 1.0 1.0 2.0];
+# [2.0 2.0 2.0 2.0 2.0];
+# [2.0 2.0 2.0 2.0 2.0];
+# [1.0 2.0 1.0 0.0 1.0];
+# [2.0 2.0 1.0 1.0 2.0]
+# ])
+
+# EURsub = @view EUR[1:2:100, 1:2:100]
+# EURsubbm = SnpBitMatrix{Float64}(EURsub, model=ADDITIVE_MODEL, center=true, scale=true) # BitMatrix from the SubArray
+# EURsubfm = convert(Matrix{Float64}, EURsub, model=ADDITIVE_MODEL, center=true, scale=true) # Float64 Matrix from the SubArray
+
+# v1 = randn(size(EURsub, 1))
+# v2 = randn(size(EURsub, 2))
+
+# @test isapprox(EURsubbm * v2, EURsubfm * v2)
+# @test isapprox(EURsubbm' * v1, EURsubfm' * v1)
+# end
+
+# @testset "split-merge-readwrite" begin
+# EUR_data = SnpData(SnpArrays.datadir("EUR_subset"))
+# # a small subset of EUR_data
+# SnpArrays.filter(EUR_data, collect(1:10), collect(1:20000); des="EUR_subset.small")
+# EUR_small = SnpData("EUR_subset.small")
+
+# # filter
+# chr17 = SnpArrays.filter(SnpArrays.datadir("EUR_subset"); des="tmp.filter.chr17", f_snp = x -> String(x[:chromosome]) == "17")
+# @test chr17.snps == 11041
+# @test chr17.people == 379
+# male = SnpArrays.filter(EUR_data; des="tmp.filter.male", f_person = x -> String(x[:sex]) == "1")
+# @test male.snps == 54051
+# @test male.people == 178
+# chr17_male = SnpArrays.filter(EUR_data; des="tmp.filter.chr17.male", f_snp = x -> String(x[:chromosome]) == "17", f_person = x -> String(x[:sex]) == "1")
+# @test chr17_male.snps == 11041
+# @test chr17_male.people == 178
+# @test size(chr17_male.snparray) == (178, 11041)
+# # cleanup
+# for ft in ["bim", "fam", "bed"]
+#     ft == "bed" && Sys.iswindows() && continue
+#     rm("tmp.filter.chr17." * ft, force=true)
+#     rm("tmp.filter.male." * ft, force=true)
+#     rm("tmp.filter.chr17.male." * ft, force=true)
+# end
+
+# # split
+# splitted = SnpArrays.split_plink("EUR_subset.small"; prefix="tmp.split.chr.")
+# piece = splitted["17"]
+# @test piece.people == 10
+# @test piece.snps == 11041
+# @test size(piece.person_info) == (10, 6)
+# @test size(piece.snp_info) == (11041, 6)
+# @test size(piece.snparray) == (10, 11041)
+# @test size(piece.snparray.columncounts) == (4, 11041)
+# @test piece.snparray.m == 10
+
+# splitted_bysex = SnpArrays.split_plink(EUR_data, :sex; prefix="tmp.split.sex.")
+# @test splitted_bysex["1"].people == 178
+# @test splitted_bysex["2"].people == 201
+
+# # merge
+# #@time merged = SnpArrays.merge_plink("tmp.merged", splitted) # write_plink is included here
+# #@test EUR_data.people == merged.people
+# #@test EUR_data.snps == merged.snps
+# #@test EUR_data.person_info == merged.person_info
+# #@test EUR_data.snp_info == merged.snp_info # note: the ordering of merged data might be different on other dataset b/c sorted order of chromosomes
+
+# #output = SnpData("tmp.merged")
+# #@test EUR_data.people == output.people
+# #@test EUR_data.snps == output.snps
+# #@test EUR_data.person_info == output.person_info
+# #@test EUR_data.snp_info == output.snp_info
+
+# # cleanup
+# #rm("tmp.merged.bim", force=true)
+# #rm("tmp.merged.fam", force=true)
+# #Sys.iswindows() || rm("tmp.merged.bed", force=true)
+
+# # merge from splitted files
+# @time merged_from_splitted_files = merge_plink("tmp.split.chr"; des = "tmp2.merged")
+# @test EUR_small.people == merged_from_splitted_files.people
+# @test EUR_small.snps == merged_from_splitted_files.snps
+# @test EUR_small.person_info == merged_from_splitted_files.person_info
+# @test EUR_small.snp_info == merged_from_splitted_files.snp_info
+
+# # cleanup
+# rm("tmp2.merged.bim", force=true)
+# rm("tmp2.merged.fam", force=true)
+# Sys.iswindows() || rm("tmp2.merged.bed", force=true)
+
+# rm("EUR_subset.small.bim", force=true)
+# rm("EUR_subset.small.fam", force=true)
+# Sys.iswindows() || rm("EUR_subset.small.bed", force=true)
+
+# for k in keys(splitted)
+#     for ft in ["bim", "fam", "bed"]
+#         ft == "bed" && Sys.iswindows() && continue
+#         rm("tmp.split.chr.$(k)." * ft, force=true)
+#     end
+# end
+# for k in keys(splitted_bysex)
+#     for ft in ["bim", "fam", "bed"]
+#         ft == "bed" && Sys.iswindows() && continue
+#         rm("tmp.split.sex.$(k)." * ft, force=true)
+#     end
+# end
+# end
+
+# @testset "(de)compress" begin
+# for format in SnpArrays.ALLOWED_FORMAT
+#     # compress mouse Plink files
+#     compress_plink(SnpArrays.datadir("mouse"), format)
+#     @test isfile(SnpArrays.datadir("mouse.bed." * format))
+#     @test isfile(SnpArrays.datadir("mouse.fam." * format))
+#     @test isfile(SnpArrays.datadir("mouse.bim." * format))
+#     # read in compressed Plink files
+#     mouse_zip = SnpArray(SnpArrays.datadir("mouse.bed." * format))
+#     fill!(mouse.rowcounts, 0)
+#     fill!(mouse.columncounts, 0)
+#     @time begin
+#         SnpArrays._counts(mouse, 1)
+#         SnpArrays._counts(mouse, 2)
+#     end
+#     fill!(mouse_zip.rowcounts, 0)
+#     fill!(mouse_zip.columncounts, 0)
+#     @time begin
+#         SnpArrays._counts(mouse_zip, 1)
+#         SnpArrays._counts(mouse_zip, 2)
+#     end
+#     @test mouse.m == mouse_zip.m
+#     @test norm(mouse.rowcounts - mouse_zip.rowcounts) < 1e-8
+#     @test norm(mouse.columncounts - mouse_zip.columncounts) < 1e-8
+#     @test norm(mouse.data - mouse_zip.data) < 1e-8
+#     # Decompress Plink files
+#     decompress_plink(SnpArrays.datadir("mouse"), format, SnpArrays.datadir("mouse2"))
+#     @test stat(SnpArrays.datadir("mouse") * ".bed").size ==
+#     stat(SnpArrays.datadir("mouse2") * ".bed").size
+#     @test stat(SnpArrays.datadir("mouse") * ".fam").size ==
+#     stat(SnpArrays.datadir("mouse2") * ".fam").size
+#     @test stat(SnpArrays.datadir("mouse") * ".bim").size ==
+#     stat(SnpArrays.datadir("mouse2") * ".bim").size
+#     # clean up
+#     rm(SnpArrays.datadir("mouse.bed." * format), force=true)
+#     rm(SnpArrays.datadir("mouse.fam." * format), force=true)
+#     rm(SnpArrays.datadir("mouse.bim." * format), force=true)
+#     rm(SnpArrays.datadir("mouse2.bed"), force=true)
+#     rm(SnpArrays.datadir("mouse2.fam"), force=true)
+#     rm(SnpArrays.datadir("mouse2.bim"), force=true)
+# end
+# @test_throws ArgumentError SnpArray(SnpArrays.datadir("mouse.bed.zip"))
+# end
+
+# @testset "indexin" begin
+#     b = ['1', '2', '3', '4']
+#     for a in [['4', '3', '1', '2'], ['3', '2', '6'], ['b', 'c', 'd', 'c']]
+#         aind, bmask = SnpArrays.indexin_general(a, b)
+#         @test all(a[aind] .== b[bmask])
+#     end
+# end
+
+# @testset "concat" begin
+# s = SnpArrays.filter(SnpArrays.datadir("mouse"), 1:2, 1:3; des="mouse.tmp.filtered")
+# sd = SnpData("mouse.tmp.filtered")
+# sd_hcat = hcat(sd, sd, sd; des="mouse.hcat")
+# sd_vcat = vcat(sd, sd, sd; des="mouse.vcat")
+# sd_hvcat = hvcat((2, 2), sd, sd, sd, sd; des="mouse.hvcat")
+
+# ref = [[0x02 0x02 0x02];
+# [0x02 0x02 0x03]]
+
+# @test all(sd_hcat.snparray .== [ref ref ref])
+# @test all(sd_vcat.snparray .== [ref; ref; ref])
+# @test all(sd_hvcat.snparray .== [ref ref; ref ref])
+
+# @test sd_hcat.people == 2
+# @test sd_hcat.snps == 9
+
+# @test sd_vcat.people == 6
+# @test sd_vcat.snps == 3
+
+# @test sd_hvcat.people == 4
+# @test sd_hvcat.snps == 6
+
+# for ft in ["bim", "fam", "bed"]
+#     ft == "bed" && Sys.iswindows() && continue
+#     rm("mouse.tmp.filtered." * ft, force=true)
+#     rm("mouse.hcat." * ft, force=true)
+#     rm("mouse.vcat." * ft, force=true)
+#     rm("mouse.hvcat." * ft, force=true)
+# end
+# end
+
+# @testset "reorder" begin
+# mouse_prefix = SnpArrays.datadir("mouse")
+# run(`cp $(mouse_prefix * ".bed") mouse_testreorder.bed`)
+# run(`cp $(mouse_prefix * ".bim") mouse_testreorder.bim`)
+# run(`cp $(mouse_prefix * ".fam") mouse_testreorder.fam`)
+
+# mouse_data = SnpData(mouse_prefix)
+# mouse_toreorder = SnpData("mouse_testreorder", "r+")
+# m, n = size(mouse_toreorder.snparray)
+# using Random
+# ind = randperm(m);
+# SnpArrays.reorder!(mouse_toreorder, ind)
+
+# # reread the file
+# mouse_toreorder_read = SnpData("mouse_testreorder"; famnm="mouse_testreorder.reordered.fam")
+# @test all(mouse_data.snparray[ind, :] .== mouse_toreorder_read.snparray[:, :])
+# @test all(mouse_data.person_info[ind, 1] .== mouse_toreorder_read.person_info[:, 1])
+
+# # cleanup
+# rm("mouse_testreorder.bed", force=true)
+# rm("mouse_testreorder.bim", force=true)
+# rm("mouse_testreorder.fam", force=true)
+# rm("mouse_testreorder.reordered.fam", force=true)
+# end
+
+# @testset "vcf2plink" begin
+# # Download an example VCF file
+# isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
+#     joinpath(pwd(), "test.08Jun17.d8b.vcf.gz"));
+# s = vcf2plink("test.08Jun17.d8b.vcf.gz", "test.08Jun17.d8b")
+# @test size(s, 1) == 191
+# @test size(s, 2) == 1354
+# @test s[1, 1] == 0x00
+# @test s[1, 5] == 0x02
+# @test s[38, 5] == 0x03
+# rm("test.08Jun17.d8b.vcf.gz", force=true)
+# rm("test.08Jun17.d8b.bed", force=true)
+# rm("test.08Jun17.d8b.bim", force=true)
+# rm("test.08Jun17.d8b.fam", force=true)
+# end
+
+# @testset "kinship_pruning" begin
+# g = grm(mouse)
+# @test count(kinship_pruning(g; method=:gcta)) == 68
+# @test count(kinship_pruning(g; method=:top_down)) == 125
+# @test count(kinship_pruning(g; method=:bottom_up)) == 132
+# @test count(kinship_pruning(g; method=:plink)) == 126
+# end
+
+# @testset "stackedsnparray" begin
+# mouse2 = StackedSnpArray([mouse, mouse])
+# @test size(mouse2) == (1940, 20300)
+# @test size(mouse2, 3) == 1
+# @test_throws ErrorException size(mouse2, 0)
+# @test length(mouse2) == 1940 * 20300
+# @test eltype(mouse2) == UInt8
+# @test mouse2[777, 16384] == mouse[777, 16384 - 10150]
+# end
+
+using Random, Mmap, GeneticVariantBase
+abstract type VariantIterator end
+abstract type Variant end
+#include(joinpath(@__DIR__, "..", "src", "iterator.jl"))
+
+# # need to export outside of module 
+
+# const EUR_data = SnpData(SnpArrays.datadir("EUR_subset"))
+# ind = SnpArrays.SnpArrayIndex(1)
+# snp_iter = SnpArrays.SnpArrayIterator(EUR_data)
+
+# @testset "n_variant and n_samples" begin
+#     @test GeneticVariantBase.n_samples(EUR_data) == 379
+#     @test GeneticVariantBase.n_variants(EUR_data) == 54051
+# end
+
+
+
+# @testset "SNP Array Tests" begin
+#     # Initialize necessary variables
+#     snp_iter = SnpArrays.SnpArrayIterator(EUR_data)
+#     num_snps = length(snp_iter)
+#     @test length(snp_iter) == 54051
+    
+#     # Iterate through each SNP in the SNP array
+#     for i in 1:10
+#         # Obtain the SnpArrayIndex object at the current iteration
+#         snp_index, _ = iterate(snp_iter, i)
+        
+#         # Test each function with the SNP data
+#         # @test chrom(EUR_data, snp_index) == EUR_data.snp_info[i,:chromosome]
+#         # println(chrom(EUR_data, snp_index))
+#         # @test pos(EUR_data, snp_index) == EUR_data.snp_info[i,:position]
+#         # @test rsid(EUR_data, snp_index) == EUR_data.snp_info[i,:snpid]
+#         # @test alleles(EUR_data, snp_index) == [EUR_data.snp_info[i,:allele1],EUR_data.snp_info[i,:allele2]]
+#         # @test alt_allele(EUR_data, snp_index) == EUR_data.snp_info[i,:allele2]
+#         # println(alt_allele(EUR_data, snp_index))
+#         # @test ref_allele(EUR_data, snp_index) == EUR_data.snp_info[i,:allele1]
+#         # println(ref_allele(EUR_data, snp_index))
+#     end
+# end
+
+
+# Test for maf function
+
+# Compute the MAF values using maf
+
+# @testset "MAF Values" begin
+#     snp_iterator = SnpArrays.SnpArrayIterator(EUR_data)
+#     maf_values = SnpArrays.maf(EUR)
+#     maf_data = SnpArrays.calculate_maf_data(EUR_data)
+#     test = SnpArrays.SnpArrayIndex(1)
+#     SnpArrays.maf_index(maf_data, test)
+
+#     for idx in eachindex(maf_values)
+#         s = SnpArrays.SnpArrayIndex(idx)
+#         println("index: $idx")
+#         @test SnpArrays.maf_index(maf_data, s) == maf_values[idx]
+#         @test GeneticVariantBase.maf(EUR_data,s) == maf_values[idx]
+#     end
+# end
+
+# snp_iterator = iterator(EUR_data)
+
+# @testset "HWE P Values" begin
+    
+#     genotypes = EUR_data.snparray[:,1]
+#     n00_count = 0
+#     n01_count = 0
+#     n11_count = 0
+
+#     for j in genotypes
+#         if j == 0x00
+#             n00_count += 1
+#         elseif j == 0x02
+#             n01_count += 1
+#         elseif j == 0x03
+#             n11_count += 1
+#         end
+#     end
+#     n00 = sum(genotypes .== 0x00) 
+#     n01 = sum(genotypes .== 0x02) 
+#     n11 = sum(genotypes .== 0x03) 
+#     @test n00 == n00_count
+#     @test n01 == n01_count 
+#     @test n11 == n11_count 
+# end
+
+# @testset "Alt Dosages & Alt Genotypes" begin
+
+#     v = Vector{Float64}(undef, size(EUR, 1))
+#     v2 = Vector{Float64}(undef, size(EUR, 1))
+#     s = SnpArrayIndex(1)
+#     alt_dosages!(v,EUR_data,s)
+#     copyto!(v2, @view(EUR[:, 1]))
+#     @test v == v2 
+#     alt_genotypes!(v,EUR_data,s)
+#     @test v == v2 
+
+# end
+
+# EUR_data = SnpData(joinpath(SnpArrays.datadir(), "EUR_subset"))
+# @testset "n_sample and n_variants" begin
+#     expected_n_samples = size(EUR, 1)
+#     expected_n_variants = size(EUR, 2)
+
+#     @test GeneticVariantBase.n_samples(EUR_data) == expected_n_samples
+#     @test GeneticVariantBase.n_variants(EUR_data) == expected_n_variants
+
+#     @test isa(n_samples(EUR_data), Int)
+#     @test isa(n_variants(EUR_data), Int)
+
+#     @test n_samples(EUR_data) > 0
+#     @test n_variants(EUR_data) > 0
+# end
\ No newline at end of file