Navigation Menu

Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow more general eltypes in sparse array multiplication #33205

Merged
merged 10 commits into from Sep 20, 2019

Conversation

pablosanjose
Copy link
Contributor

@pablosanjose pablosanjose commented Sep 9, 2019

Fixes #33169

Replaces one(T) by 1 in five-argument mul! calls with SparseArray types, where T was before the target eltype. This allows eltypes that don't define one by relying on the fact that 1 is always a multiplicative identity for AbtractArrays as far as I can tell (even for Bool eltypes) EDIT: changed 1, 0 to true, false.

With this we can in particular use StaticArrays as eltypes of SparseMatrixCSC. It does not allow Array eltypes, since they don't define zero, which is used in mul!. Solving the issue also for eltypes that don't define zero is of much broader scope than the present PR due to the pervasive assumption throughout SparseArrays that zero(eltype(A)) for A::AbstractSparseArray is defined.

As a consequence of allowing general (e.g. matrix-like) eltypes in all product combinations involving a sparse vector or matrices (see tests) a range of changes throughout SparseArrays were needed, including in copy(::AbstractSparseMatrix), halfperm!, spmatmul and _At_or_Ac_mul_B.

This fix does not address a possible version of #33169 without any sparse factors (i.e. all changes here are strictly confined to the SparseArrays lib)

(This is my first PR to Julia stdlib. It might be a good idea to show this to Nanosoldier, and to have some careful review if possible)

copy(::Adjoint/Transpose)+ simplesmatrix fixes

 added news

fix promote_op for sparsevec mul!

added + to SimpleSMatrix (required for matprod)

trying to make mat * vec work :-(

_At_or_Ac_mul_B! fix to accept target eltype

more simplesmatrix fixes

yet more fixes to simplesmatrix

more stringent tests
@pablosanjose
Copy link
Contributor Author

(I see some whitespace fixes are needed... Can somebody briefly point out how to address that?)

@ViralBShah ViralBShah added the sparse Sparse arrays label Sep 9, 2019
@fredrikekre
Copy link
Member

Can somebody briefly point out how to address that?

Remove trailing whitespaces.

@pablosanjose pablosanjose changed the title Allow general eltypes that implement zero in sparse mul! Allow more general eltypes in sparse array multiplication Sep 10, 2019
@andreasnoack
Copy link
Member

Seems reasonable to me. @tkf any thoughts on this? Should we use Booleans here?

@pablosanjose
Copy link
Contributor Author

I think you are right, booleans are in fact better. I somehow though I saw true * SArray{Bool}(...) being promoted to SArray{Int}, but it is not. The universal one and zero is indeed true and false. I'll push an update with booleans

@tkf
Copy link
Member

tkf commented Sep 11, 2019

I don't think it matters too much here as we are not using α=false case (which matters when the inputs produce NaNs). But maybe it'd play nice with some custom eltype implementing GF(2) or something? I think it's good to use Booleans anyway for consistency. (Thank you @pablosanjose for the change!)

@pablosanjose
Copy link
Contributor Author

pablosanjose commented Sep 11, 2019

Funny, there is a strange win32 failure that didn't get triggered before the last commit. I assume this is unrelated to this PR?
EDIT: a restart fixed it

@pablosanjose
Copy link
Contributor Author

A friendly bump looking for some volunteer to review if possible. Thanks!

Copy link
Sponsor Member

@timholy timholy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice. Completely optional to get rid of more type parameters, this could basically be merged as is.

stdlib/SparseArrays/src/linalg.jl Show resolved Hide resolved
stdlib/SparseArrays/src/sparsevector.jl Outdated Show resolved Hide resolved
@timholy
Copy link
Sponsor Member

timholy commented Sep 18, 2019

LGTM. I should allow someone else take a look, but if you don't hear soon let's merge.

stdlib/SparseArrays/src/sparsevector.jl Outdated Show resolved Hide resolved
stdlib/SparseArrays/test/simplesmatrix.jl Show resolved Hide resolved
pablosanjose and others added 2 commits September 18, 2019 16:53
Co-Authored-By: Fredrik Ekre <ekrefredrik@gmail.com>
Co-Authored-By: Fredrik Ekre <ekrefredrik@gmail.com>
@pablosanjose
Copy link
Contributor Author

Many thanks to both for the review and approval!

@pablosanjose
Copy link
Contributor Author

Could this be merged, perhaps?

@fredrikekre fredrikekre merged commit 34d2b87 into JuliaLang:master Sep 20, 2019
@maleadt maleadt mentioned this pull request Jan 15, 2020
28 tasks
@KristofferC
Copy link
Sponsor Member

This seems like it broke the tests for JWAS.jl:

ERROR: LoadError: LoadError: MethodError: no method matching zero(::Type{Any})
Closest candidates are:
  zero(::Type{Union{Missing, T}}) where T at missing.jl:105
  zero(::Type{Missing}) at missing.jl:103
  zero(::Type{LibGit2.GitHash}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LibGit2/src/oid.jl:220
  ...
Stacktrace:
 [1] zero(::Type{Any}) at ./missing.jl:105
 [2] _zeros_eltypes at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/SparseArrays/src/higherorderfns.jl:204 [inlined]
 [3] _noshapecheck_map(::typeof(+), ::SparseArrays.SparseMatrixCSC{Any,Int64}, ::SparseArrays.SparseMatrixCSC{Float32,Int64}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/SparseArrays/src/higherorderfns.jl:159
 [4] map(::typeof(+), ::SparseArrays.SparseMatrixCSC{Any,Int64}, ::SparseArrays.SparseMatrixCSC{Float32,Int64}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/SparseArrays/src/higherorderfns.jl:1153
 [5] +(::SparseArrays.SparseMatrixCSC{Any,Int64}, ::SparseArrays.SparseMatrixCSC{Float32,Int64}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/SparseArrays/src/sparsematrix.jl:1648
 [6] addVinv(::JWAS.MME) at /Users/kristoffer/.julia/packages/JWAS/4ltS3/src/1.JWAS/src/buildMME/random_effects.jl:205
 [7] getMME(::JWAS.MME, ::DataFrame) at /Users/kristoffer/.julia/packages/JWAS/4ltS3/src/1.JWAS/src/buildMME/build_MME.jl:264
 [8] init_mixed_model_equations(::JWAS.MME, ::DataFrame, ::Bool) at /Users/kristoffer/.julia/packages/JWAS/4ltS3/src/1.JWAS/src/JWAS.jl:565
 [9] runMCMC(::JWAS.MME, ::DataFrame; heterogeneous_residuals::Bool, chain_length::Int64, starting_value::Bool, burnin::Int64, output_samples_frequency::Int64, output_samples_file::String, update_priors_frequency::Int64, methods::String, estimate_variance::Bool, Pi::Float64, estimatePi::Bool, estimateScale::Bool, single_step_analysis::Bool, pedigree::Bool, categorical_trait::Bool, missing_phenotypes::Bool, constraint::Bool, causal_structure::Bool, outputEBV::Bool, output_heritability::Bool, seed::Int64, printout_model_info::Bool, printout_frequency::Int64, big_memory::Bool, double_precision::Bool) at /Users/kristoffer/.julia/packages/JWAS/4ltS3/src/1.JWAS/src/JWAS.jl:234
 [10] top-level scope at /Users/kristoffer/.julia/packages/JWAS/4ltS3/test/test_BayesianAlphabet.jl:100
 [11] include(::String) at ./client.jl:439
 [12] top-level scope at /Users/kristoffer/.julia/packages/JWAS/4ltS3/test/runtests.jl:3
 [13] include(::String) at ./client.jl:439
 [14] top-level scope at none:6

@pablosanjose
Copy link
Contributor Author

pablosanjose commented Jan 15, 2020

First of all, my apologies, this was my first merged PR and I feel bad!
I'm trying to understand the reason for the failure. It seems the stacktrace involves codepaths in files not touched by this PR, except [5], which involves + in sparsematrix.jl, but this is a function not modified here. Any further clue why that failure is related to the PR?

@KristofferC
Copy link
Sponsor Member

KristofferC commented Jan 15, 2020

First of all, my apologies, this was my first merged PR and I feel bad!

Please don't! This happens all the time, this just either means that we lack some tests in Base or that the package does something "weird". Consider it a rite of passage ;)
I'll try figure out a bit better what is happening, I just posted it directly here for reference and so someone else can also look into it if they feel like it.

The reason for suspecting this PR is that it was bisected here: #34238 (comment). As you can see in that thread, we started with over 200 new breakages for 1.4 so breaking packages is definitely not a big deal.

@KristofferC
Copy link
Sponsor Member

KristofferC commented Jan 15, 2020

So the difference is in:

julia> a = convert(SparseMatrixCSC{AbstractFloat, Int}, sparse(rand(Float32, 2,2)));

julia> b = sparse(rand(Float32, 2,2));

julia> typeof(a*b)
SparseMatrixCSC{Any,Int64}

On 1.3, this returned a SparseMatrixCSC{AbstractFloat,Int64}. I'm not sure what is "correct" here...

@KristofferC
Copy link
Sponsor Member

I made a PR to the package to fix it: reworkhow/JWAS.jl#59.

Not sure if we need to change anything in Julia.

@pablosanjose
Copy link
Contributor Author

I see! With this PR the target eltype of the product a*b is computed with

julia> Base.promote_op(LinearAlgebra.matprod, AbstractFloat, Float32)
Any

while in 1.3 we merely did

julia> promote_type(AbstractFloat, Float32)
AbstractFloat

Probably the correct solution within Julia is to make Base.promote_op smarter, returning AbstractFloat in this case. What do you think?

@pablosanjose
Copy link
Contributor Author

Uhm, that promote_op method goes into the return_type abyss where there be dragons. Such a fix is far beyond the red line for me.

@KristofferC
Copy link
Sponsor Member

Yeah, promote_op is known to be a bit of a wart (for that reason).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
sparse Sparse arrays
Projects
None yet
Development

Successfully merging this pull request may close these issues.

SparseMatrixCSC * Vector not supported with StaticArray eltypes
7 participants