OpenMendel / SnpArrays.jl

Compressed storage for SNP data
https://openmendel.github.io/SnpArrays.jl/latest
Other
44 stars 9 forks source link

SnpLinAlg broken for some data on master #87

Closed biona001 closed 3 years ago

biona001 commented 3 years ago

On master branch, SnpLinAlg works on EUR data as usual

using SnpArrays, LinearAlgebra
x = SnpArray(SnpArrays.datadir("EUR_subset.bed"));
xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, impute=true, center=true, scale=true);
n, p = size(x)

# both works
xla * rand(size(xla, 2));
LinearAlgebra.mul!(zeros(n), xla, rand(p));

but I get weird error on my test data (master branch)

using SnpArrays, LinearAlgebra
x = SnpArray("/Users/biona001/Benjamin_Folder/UCLA/research/stampeed/by_chromosome/split.chr.21.bed")
xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)
n, p = size(x)

julia> xla * rand(size(xla, 2))
ERROR: MethodError: no method matching getindex(::VectorizationBase.StridedPointer{UInt8, 2, 1, 0, (1, 2), Tuple{Static.StaticInt{1}, Int64}, Tuple{Static.StaticInt{0}, Static.StaticInt{0}}}, ::Int64, ::Int64)
Stacktrace:
 [1] _snparray_ax_additive!
   @ ~/.julia/packages/SnpArrays/8SdVC/src/linalg_direct.jl:0 [inlined]
 [2] _snparray_ax_tile!(c::Vector{Float64}, A::Matrix{UInt8}, b::Vector{Float64}, model::Val{1}, μ::Vector{Float64}, impute::Bool)
   @ SnpArrays ~/.julia/packages/SnpArrays/8SdVC/src/linalg_direct.jl:211
 [3] mul!(out::Vector{Float64}, sla::SnpLinAlg{Float64}, v::Vector{Float64})
   @ SnpArrays ~/.julia/packages/SnpArrays/8SdVC/src/linalg_direct.jl:137
 [4] *(A::SnpLinAlg{Float64}, x::Vector{Float64})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:51
 [5] top-level scope
   @ REPL[12]:1

julia> LinearAlgebra.mul!(rand(n), xla, rand(p))
ERROR: MethodError: no method matching getindex(::VectorizationBase.StridedPointer{UInt8, 2, 1, 0, (1, 2), Tuple{Static.StaticInt{1}, Int64}, Tuple{Static.StaticInt{0}, Static.StaticInt{0}}}, ::Int64, ::Int64)
Stacktrace:
 [1] _snparray_ax_additive!
   @ ~/.julia/packages/SnpArrays/8SdVC/src/linalg_direct.jl:0 [inlined]
 [2] _snparray_ax_tile!(c::Vector{Float64}, A::Matrix{UInt8}, b::Vector{Float64}, model::Val{1}, μ::Vector{Float64}, impute::Bool)
   @ SnpArrays ~/.julia/packages/SnpArrays/8SdVC/src/linalg_direct.jl:211
 [3] mul!(out::Vector{Float64}, sla::SnpLinAlg{Float64}, v::Vector{Float64})
   @ SnpArrays ~/.julia/packages/SnpArrays/8SdVC/src/linalg_direct.jl:137
 [4] top-level scope
   @ REPL[13]:1

It seems d63c016 broke the test case above. Up until 39566d1 the above code works. Also, the test data is real data so I cannot upload here. But I can send it privately.

kose-y commented 3 years ago

I'm inspecting the issue. This issue is likely because of a breaking change made in VectorizationBase.jl, and may take some time to resolve.

biona001 commented 3 years ago

I think the problem is this line.

Here typeof(s) is VectorizationBase.StridedPointer{UInt8, 2, 1, 0, (1, 2), Tuple{StaticInt{1}, Int64}, Tuple{StaticInt{0}, StaticInt{0}}} which does not support indexing. However indexing becomes possible if we add @avx in front of the for loop, evident if we compare with Chris Elrod's gemv_avx_512by512! in this issue. It seems d63c016 removed the @avx? Was that intentional?