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

Implement sum #446

Merged
merged 14 commits into from
Jul 16, 2024
Merged

Conversation

max-vassili3v
Copy link
Contributor

I struggled to find an elegant solution involving selecting relevant elements from B.data in the case where B is created by brand() with certain parameters (e.g very non square matrices, more bands than those that fit in the matrix). I decided to go with this solution that involves populating a data matrix only using relevant information accessed by B[band(i)]. Please let me know any improvements.

Copy link

codecov bot commented Jul 8, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 89.75%. Comparing base (47c15ab) to head (5ec4b3c).

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #446      +/-   ##
==========================================
+ Coverage   89.61%   89.75%   +0.14%     
==========================================
  Files          25       25              
  Lines        3571     3622      +51     
==========================================
+ Hits         3200     3251      +51     
  Misses        371      371              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@dlfivefifty
Copy link
Member

Can you make sure your unit tests are being run? I think just add include("test_sum.jl") to runtests.jl

src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
@max-vassili3v
Copy link
Contributor Author

added changes

src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
src/banded/BandedMatrix.jl Outdated Show resolved Hide resolved
@dlfivefifty
Copy link
Member

The test failure is very surprising. It's possible its calling sum and this PR has changed the order of operations.

We should make sure we are consistent with the order. Can you add some tests that sum(A) == sum(Matrix(A)) for A a banded matrix?

@max-vassili3v
Copy link
Contributor Author

I noticed earlier with == that there is some floating point error and so the test fails.

@max-vassili3v
Copy link
Contributor Author

I thought it was to be expected but this could be the problem

@dlfivefifty
Copy link
Member

Right, this floating point error is because you are computing the sums in a different order. This is unnecessary so we can change the implementation to make sure we do things in the right order. Eg:

julia> A = randn(5,5)
5×5 Matrix{Float64}:
 -0.574303  -0.909723    0.589035   0.125461  -0.85839
  2.36645   -2.01842     0.305596   0.739664   0.281112
 -0.449434   1.65078     0.293241  -0.12409    0.535829
  0.388728   1.3232     -1.61161   -0.54598   -0.237829
 -0.570773  -0.0989053  -0.515742   0.116799  -2.14109

julia> sum(A) == sum(vec(A)) # sum should traverse column-by-coumn
true

julia> sum(A)  sum(vec(A')) # it doesn't match row-by-row 
true

@dlfivefifty
Copy link
Member

I thought it was to be expected but this could be the problem

It is expected when the order of the operations change. But all things being equal we should avoid it.

Also, traversing column-by-column will be much faster than row-by-row since it accesses memory in order.

@max-vassili3v
Copy link
Contributor Author

I've changed the traversal order but I still get floating point error on the tests without dims and dims = 1. It's a different type, but I've also noticed that sum(vec(A)) == sum(Matrix(A)) for BandedMatrix A returns false

@dlfivefifty
Copy link
Member

Can you push your changes? Note I made a suggestion that fixes the order for the no-dims case

@dlfivefifty dlfivefifty self-requested a review July 11, 2024 13:35
@DanielVandH
Copy link
Contributor

It's a different type, but I've also noticed that sum(vec(A)) == sum(Matrix(A)) for BandedMatrix A returns false

Wouldn't a better check for this probably be that sum(vec(A)) == foldl(+, Matrix(A))?

This is a similar problem in e.g. SparseArrays where A = sprand(1000, 1000, 0.001); sum(vec(A)) == sum(Matrix(A)) returns false but A = sprand(1000, 1000, 0.001); sum(vec(A)) == foldl(+, Matrix(A)) is true

@dlfivefifty
Copy link
Member

Can you explain the difference?

I take it foldl forces a specific order. Do you know why sum might choose a different order?

@DanielVandH
Copy link
Contributor

DanielVandH commented Jul 14, 2024

As far as I can tell, the difference is that IndexStyle(vec(A)) = IndexCartesian() (since, for sparse arrays and banded matrices, vec(A) becomes a reshape type unlike a normal matrix where it becomes a vector) which uses mapfoldl, while Matrix(A) is an IndexLinear() which uses some sort of block-based summation.

Another way to test would be to check sum(Vector(vec(a))) == sum(Matrix(A)). I think the implementation in this PR is equivalent to doing a foldl implementation, so the tests should probably look at sum(B) == foldl(+, Matrix(B)) if I've read it correctly

@dlfivefifty
Copy link
Member

I think in this case just use ≈ since we don't care that much about preserving order


function sum!(ret::AbstractArray, A::AbstractBandedMatrix)
#Behaves similarly to Base.sum!
ret .= 0
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
ret .= 0
fill!(ret, zero(eltype(ret)))

if s[1] == 1 && (l == 1 || s[2]==1)
for j = 1:m, i = colrange(A, j)
ret .+= A[i, j]
end
Copy link
Member

Choose a reason for hiding this comment

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

Add tests fro this special case

ret[1, j] += A[i, j]
end
elseif s[1] == n && s[2] == m
ret = A
Copy link
Member

Choose a reason for hiding this comment

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

This is not changing ret!

Suggested change
ret = A
copyto!(ret, A)

elseif s[1] == n && s[2] == m
ret = A
else
throw(DimensionMismatch("reduction on matrix of size ($n, $m) with output size $s"))
Copy link
Member

Choose a reason for hiding this comment

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

Add test using @test_throws

test/test_sum.jl Outdated Show resolved Hide resolved
@max-vassili3v
Copy link
Contributor Author

seems to be the same issue as before with the floating point error

@dlfivefifty
Copy link
Member

That error is unrelated I believe, we have seen it other places.

@dlfivefifty dlfivefifty merged commit ea616cc into JuliaLinearAlgebra:master Jul 16, 2024
16 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants