Use 3-arg mul! for Companion and improve test coverage #69

JeffFessler · 2022-10-31T22:31:51Z

This is a breaking change because mul! now is the 3-arg version, just like #68.
It also adds tests to get complete coverage.

coveralls · 2022-10-31T22:33:36Z

Pull Request Test Coverage Report for Build 3377679255

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

For more information on this, see Tracking coverage changes with pull request builds.
To avoid this issue with future PRs, see these Recommended CI Configurations.
For a quick fix, rebase this PR at GitHub. Your next report should be accurate.

Details

25 of 25 (100.0%) changed or added relevant lines in 1 file are covered.
No unchanged relevant lines lost coverage.
Overall coverage increased (+0.5%) to 94.0%

Totals
Change from base Build 3364853852:	0.5%
Covered Lines:	188
Relevant Lines:	200

💛 - Coveralls

codecov · 2022-10-31T22:35:35Z

Codecov Report

Base: 93.48% // Head: 94.02% // Increases project coverage by +0.54% 🎉

Coverage data is based on head (8c44bce) compared to base (b7c98d7).
Patch coverage: 100.00% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #69      +/-   ##
==========================================
+ Coverage   93.48%   94.02%   +0.54%     
==========================================
  Files           8        8              
  Lines         215      201      -14     
==========================================
- Hits          201      189      -12     
+ Misses         14       12       -2

Impacted Files	Coverage Δ
src/companion.jl	`100.00% <100.00%> (+4.16%)`	⬆️
src/hilbert.jl	`100.00% <0.00%> (ø)`
src/frobenius.jl	`100.00% <0.00%> (ø)`

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

JeffFessler · 2022-11-01T00:56:17Z

@ivanslapnicar this applied the 3-arg mul! to Companion as well.
Does it look OK to you?

@mkitti you were the last to edit the companion.jl files. Any comments here?

mkitti · 2022-11-01T04:00:58Z

src/companion.jl

-# getindex(C::Companion, i, j) = getindex(Matrix(C), i, j)
-# isassigned(C::Companion, i, j) = isassigned(Matrix(C), i, j)
-getindex(C::Companion{T}, i::Int, j::Int) where T = (j==length(C.c)) ? -C.c[i] : (i==j+1 ? one(T) : zero(T) )
-isassigned(C::Companion, i::Int, j::Int) = (j==length(C.c)) ? isassigned(C.c,i) : true


Why is isassigned being removed here?

HI @mkitti, thank for the very thorough review. I will address (or commit) your comments one-by-one, which might take a while, and I will ping you again when it is ready for another review.
This repo goes back to early versions of Julia and many things in the original version are no longer needed in current Julia versions. Now Base has default fallbacks for isassigned for any AbstractArray types so there is no need for this function here. To double check that that fallback is working correctly, I have added tests of this form:

@test isassigned(Z, 1) @test isassigned(Z, n^2) @test !isassigned(Z, 0) @test !isassigned(Z, n^2+1)

Now, these are all standard uses of isassigned with a single index argument per the manual. I have never seen a 2-index use of isassigned before, and it was never tested in this package. If there is a use case for it, then it's fine to restore it, but then it should have a docstring and jldoctests that explains what it does.
Note that SpecialMatrices.jl does an import of Base.isassigned and this is the only one of the special matrices that extends this method, so it really looked like a bug (or historical baggage) to be importing the 1-index version and then defining a 2-index version. Again, happy to restore it if there is a reason to have this non-standard method...

Update: I accidentally pushed those tests to master, which is embarrassing. I thought I had set up the settings here to prevent direct pushes to master without a PR. I will have to look into that. I have reverted that in master and will push them here now. Sigh.

Actually the settings do require a PR before merging so now I am baffled:

I know it is off topic but if you are a Julia wizard who knows how to make sure all merges come from a PR please let me know, so that no one including myself can make that mistake...

Here is the fallback function. While technically this works, I think we may be able to do this more efficiently.

julia> @which isassigned(Companion(1:3),1,1) isassigned(a::AbstractArray, i::Integer...) in Base at abstractarray.jl:563 function isassigned(a::AbstractArray, i::Integer...) try a[i...] true catch e if isa(e, BoundsError) || isa(e, UndefRefError) return false else rethrow() end end end

That said that previous definition does seem problematic.

julia> using SpecialMatrices julia> C = Companion(1:3) 3×3 Companion{Int64}: 0 0 -1 1 0 -2 0 1 -3 julia> isassigned(C, 1, 4) false julia> Base.isassigned(C::Companion, i::Int, j::Int) = (j==length(C.c)) ? isassigned(C.c,i) : true julia> isassigned(C, 1, 4) true

My question is does anything depend on this unusual behavior? Otherwise, I suggest we consider writing an optimized version of this.

Similarly, using the old definition of getindex changes the situation as well.

julia> Base.getindex(C::Companion{T}, i::Int, j::Int) where T = (j==length(C.c)) ? -C.c[i] : (i==j+1 ? one(T) : zero(T) ) julia> C[1,4] 0

This is also weird. Just to be clear what we are doing here is materially changing the functionality of Companion. At the moment this seems odd, but perhaps there is a reason it is like this.

julia> [C[ci] for ci in CartesianIndices((3,10))] 3×10 Matrix{Int64}: 0 0 -1 0 0 0 0 0 0 0 1 0 -2 0 0 0 0 0 0 0 0 1 -3 0 0 0 0 0 0 0

This seems amenable adding additional columns.

mkitti · 2022-11-01T04:03:24Z

src/companion.jl

    c :: Vector{T}
 end

-# Generate companion matrix from a polynomial
+Companion(v::AbstractVector) = Companion(collect(v))


Why is this new constructor necessary? Julia will automatically perform convert(Vector{T}, v) for assignment to the field.

Hmm, what version of Julia are you using that does that automatic conversion? I am using 1.8.2, and if I remove this constructor then it throws an error:

julia> Companion(1:3) ERROR: MethodError: no method matching Companion(::UnitRange{Int64}) Closest candidates are: Companion(::Polynomial) at ~/.julia/dev/SpecialMatrices/src/companion.jl:37 Companion(::Vector{T}) where T at ~/.julia/dev/SpecialMatrices/src/companion.jl:30

I'm pretty sure we need this.

I see, it's due to the parameter.

julia> Companion{Int}(1:3) 3×3 Companion{Int64}: 0 0 -1 1 0 -2 0 1 -3

mkitti

There is a lot going on here and usually my preference would be to separate the stylistic changes from the functional changes. However, let's try to proceed with this.

The biggest issue here is the introduction of AbstractVector and AbstractMatrix. For these abstract types, there are much looser guarantees on the interface. For example, they could be an OffsetArray from OffsetArrays.jl and thus may not use 1-based indexing. You may need to use Base.require_one_based_indexing. We should also test for the OffsetArray case.

The function signatures have changed. It may be good to do some analysis if this has any unintended consequences. Otherwise, could you please elucidate the expected call chain or why the removed function signatures are redundant?

src/companion.jl

mkitti · 2022-11-01T04:07:16Z

src/companion.jl

-size(C::Companion, r::Int) = (r==1 || r==2) ? length(C.c) :
-    throw(ArgumentError("Companion is of rank 2"))


Can you prove this is redundant? Have you considered constant propagation and invalidations?

I have added a test of the form size(C,3) == 1 to verify that the defaults in Base for an AbstractArray do what is needed here. I already had a test of the form size(C,1) == n.

The "proof" is that the manual for the AA interface specifies only that a single size function that returns a tuple is needed.
I have no idea about "constant propagation and invalidations." (These are static structs.)

src/companion.jl

mkitti · 2022-11-01T04:19:47Z

src/companion.jl

-    b[1] = y
-    b
+# 3-argument mul! mutates first argument: y <= C * x
+function mul!(y::Vector, C::Companion, x::AbstractVector)


You've dropped the element type parameter here. What are the implications?

The benefit is now one can multiply C and x with different element types, as long as the result vector y has a suitable eltype. If it does not, then the user will get an error, just like they would with a regular matrix:

using LinearAlgebra: mul! mul!(ones(Float64, 2,3), ones(Int32, 2, 5), rand(Float16, 5, 3)) # works fine and is convenient for user mul!(ones(Int64, 2,3), ones(Int32, 2, 5), rand(Float16, 5, 3)) # fails with InexactError, understandably

(Note that base Julia does not fail with a Method error here.)

src/companion.jl

Co-authored-by: Mark Kittisopikul <[email protected]>

…es.jl into jf/companion

JeffFessler · 2022-11-02T02:47:00Z

@mkitti this is ready for your re-review now. Thanks!

JeffFessler · 2022-11-02T12:59:23Z

@mkitti ok, I added the OffsetArray test. Any remaining concerns?

mkitti · 2022-11-02T16:04:16Z

Could you test multiplication with an OffsetMatrix ?

JeffFessler · 2022-11-02T17:43:20Z

Could you test multiplication with an OffsetMatrix ?

Ok, I added a @test_throws because such multiplication is unsupported. I hope that's what you meant.
(Even base Julia seems not to support multiplication with an OffsetMatrix, understandably.)

mkitti

Rather than forcing a collect pass to the parameterized constructor.

mkitti · 2022-11-02T22:50:13Z

src/companion.jl

    c :: Vector{T}
 end

-# Generate companion matrix from a polynomial
+Companion(v::AbstractVector) = Companion(collect(v))


I see, it's due to the parameter.

julia> Companion{Int}(1:3) 3×3 Companion{Int64}: 0 0 -1 1 0 -2 0 1 -3

src/companion.jl

mkitti

getindex and isassigned change in function here.

Let's review the dependencies to see if anyone is depending on the anomalous prior functionality.

This is also technically a breaking change, so we need to make sure to bump the version number accordingly. Could you add a version number bump to 3.0.0 in Project.toml?

mkitti · 2022-11-02T23:41:34Z

src/companion.jl

-# getindex(C::Companion, i, j) = getindex(Matrix(C), i, j)
-# isassigned(C::Companion, i, j) = isassigned(Matrix(C), i, j)
-getindex(C::Companion{T}, i::Int, j::Int) where T = (j==length(C.c)) ? -C.c[i] : (i==j+1 ? one(T) : zero(T) )
-isassigned(C::Companion, i::Int, j::Int) = (j==length(C.c)) ? isassigned(C.c,i) : true


Here is the fallback function. While technically this works, I think we may be able to do this more efficiently.

julia> @which isassigned(Companion(1:3),1,1) isassigned(a::AbstractArray, i::Integer...) in Base at abstractarray.jl:563 function isassigned(a::AbstractArray, i::Integer...) try a[i...] true catch e if isa(e, BoundsError) || isa(e, UndefRefError) return false else rethrow() end end end

That said that previous definition does seem problematic.

julia> using SpecialMatrices julia> C = Companion(1:3) 3×3 Companion{Int64}: 0 0 -1 1 0 -2 0 1 -3 julia> isassigned(C, 1, 4) false julia> Base.isassigned(C::Companion, i::Int, j::Int) = (j==length(C.c)) ? isassigned(C.c,i) : true julia> isassigned(C, 1, 4) true

My question is does anything depend on this unusual behavior? Otherwise, I suggest we consider writing an optimized version of this.

mkitti · 2022-11-02T23:48:08Z

src/companion.jl

-# getindex(C::Companion, i, j) = getindex(Matrix(C), i, j)
-# isassigned(C::Companion, i, j) = isassigned(Matrix(C), i, j)
-getindex(C::Companion{T}, i::Int, j::Int) where T = (j==length(C.c)) ? -C.c[i] : (i==j+1 ? one(T) : zero(T) )
-isassigned(C::Companion, i::Int, j::Int) = (j==length(C.c)) ? isassigned(C.c,i) : true


Similarly, using the old definition of getindex changes the situation as well.

julia> Base.getindex(C::Companion{T}, i::Int, j::Int) where T = (j==length(C.c)) ? -C.c[i] : (i==j+1 ? one(T) : zero(T) ) julia> C[1,4] 0

This is also weird. Just to be clear what we are doing here is materially changing the functionality of Companion. At the moment this seems odd, but perhaps there is a reason it is like this.

julia> [C[ci] for ci in CartesianIndices((3,10))] 3×10 Matrix{Int64}: 0 0 -1 0 0 0 0 0 0 0 1 0 -2 0 0 0 0 0 0 0 0 1 -3 0 0 0 0 0 0 0

This seems amenable adding additional columns.

mkitti · 2022-11-03T00:10:26Z

According to https://juliahub.com/ui/Packages/SpecialMatrices/F3uIk/2.0.1?page=2 there are only two dependent packages at the moment:
1.https://github.com/SciFracX/FractionalCalculus.jl
2. https://github.com/SciFracX/FractionalDiffEq.jl

Neither seems to use Companion. Breaking the old behavior is fine. What I would like to do is make sure the performance matches or exceeds that of Matrix.

Co-authored-by: Mark Kittisopikul <[email protected]>

JeffFessler · 2022-11-03T00:58:06Z

I committed your suggestion to use built-in convert (which I did not know about), but then the tests all failed because the convert does not handle an OffsetArray (that we test), whereas collect does handle it. So I am reverting to collect and will wait for further advice for you about that.

JeffFessler · 2022-11-03T01:21:55Z

What I would like to do is make sure the performance matches or exceeds that of Matrix.

In generally I would agree with this approach and if you want to write and run benchmarks for this then please do. My sense is that companion matrices are mainly used academically for quite small sizes to illustrate points about characteristic polynomials, so I'd prefer to put my own energy into writing and testing benchmarks in other places that are more performance critical.

JeffFessler · 2022-11-03T01:26:25Z

I added the version bump as suggested. I prefer to wait to release that version until the similarly breaking changes in #17 are reviewed and merged. If you would be willing to look at that one it would be great :)
Thanks for your research into the dependent packages - I didn't know how to do that.

mkitti · 2022-11-03T03:26:45Z

I would prefer to keep the convert and change the OffsetArray test.

mkitti

Keep the convert based mechanism, change the test.

test/companion.jl

src/companion.jl

mkitti · 2022-11-03T03:34:34Z

The reason for this is that the convert mechanism is extensible. OffsetArrays.jl could implement convert(::Type{Vector{T}}, oa::OffsetVector{T}) where T but chose not to. That's an OffsetArrays issue.

JuliaArrays/OffsetArrays.jl#239

Co-authored-by: Mark Kittisopikul <[email protected]>

JeffFessler · 2022-11-03T14:49:44Z

The reason for this is that the convert mechanism is extensible

The changes are fine with me. I don't see what makes convert preferable to collect - they both seem extensible to me, but there is much I do not know :)

mkitti · 2022-11-03T20:00:47Z

collect only returns an Array of some sort and must be invoked explicitly. It meant to construct that Array through the iteration interface.

convert can return the specified type. It is used in specific places within Julia implicitly:
https://docs.julialang.org/en/v1/manual/conversion-and-promotion/#When-is-convert-called?

When is convert called?
The following language constructs call convert:

Assigning to an array converts to the array's element type.
Assigning to a field of an object converts to the declared type of the field.
Constructing an object with new converts to the object's declared field types.
Assigning to a variable with a declared type (e.g. local x::T) converts to that type.
A function with a declared return type converts its return value to that type.
Passing a value to ccall converts it to the corresponding argument type.

JeffFessler · 2022-11-03T20:12:59Z

I see now - thanks! So is this PR fine now, or do you plan to explore the efficiency issue further first or propose more changes?

mkitti

Approved

JeffFessler · 2022-11-03T20:28:36Z

Thanks so much - I really appreciate the constructive dialog and improvements to the code.

JeffFessler added 2 commits October 31, 2022 18:28

Use 3-arg mul! for Companion

84ceb55

Test Companion

ddb265b

JeffFessler marked this pull request as draft October 31, 2022 22:31

Add test for AbstracVector

5bb1588

JeffFessler marked this pull request as ready for review November 1, 2022 00:54

mkitti reviewed Nov 1, 2022

View reviewed changes

mkitti suggested changes Nov 1, 2022

View reviewed changes

JeffFessler and others added 12 commits November 1, 2022 21:48

Test isassigned, just to be sure

5d17d52

Use begin:end

d091c16

Co-authored-by: Mark Kittisopikul <[email protected]>

Test size(C,3) == 1 just to be sure.

60c15fa

Use begin:end-1

8af5c07

Co-authored-by: Mark Kittisopikul <[email protected]>

Add jldoctest

195544d

Co-authored-by: Mark Kittisopikul <[email protected]>

Use Base.require_one_based_indexing

1bf2783

Remove unneeded T

deb4a57

Use jldoctest

73cb15b

Co-authored-by: Mark Kittisopikul <[email protected]>

Merge branch 'jf/companion' of github.com:JuliaMatrices/SpecialMatric…

2048700

…es.jl into jf/companion

Better docstring

79b67e2

Add Polynomials for docs

9b5c0a1

Fix jldoctest output

80514d5

JeffFessler added 3 commits November 2, 2022 08:41

Add OffsetArray test

3c28845

Use test/Project.toml

9760cb9

Remove unneeded "isassigned" import from Base

fa2b6f2

JeffFessler mentioned this pull request Nov 2, 2022

Check 1-based indexing throughout #71

Closed

Add @test_throws for OffsetMatrix * Companion

a52a456

mkitti suggested changes Nov 2, 2022

View reviewed changes

mkitti reviewed Nov 2, 2022

View reviewed changes

JeffFessler and others added 2 commits November 2, 2022 20:51

Use convert

56a8565

Co-authored-by: Mark Kittisopikul <[email protected]>

Revert to using collect

fb13d25

Bump version to 3.0.0

9e51d3e

mkitti suggested changes Nov 3, 2022

View reviewed changes

test/companion.jl Outdated Show resolved Hide resolved

src/companion.jl Outdated Show resolved Hide resolved

JeffFessler and others added 3 commits November 3, 2022 10:41

Update test/companion.jl to use no_offset_view

0058198

Co-authored-by: Mark Kittisopikul <[email protected]>

Update src/companion.jl to use convert again

2836dd7

Co-authored-by: Mark Kittisopikul <[email protected]>

Add import

8c44bce

mkitti approved these changes Nov 3, 2022

View reviewed changes

JeffFessler merged commit accdbf7 into master Nov 3, 2022

JeffFessler deleted the jf/companion branch November 3, 2022 20:27

		size(C::Companion, r::Int) = (r==1 \|\| r==2) ? length(C.c) :
		throw(ArgumentError("Companion is of rank 2"))

Use 3-arg mul! for Companion and improve test coverage #69

Use 3-arg mul! for Companion and improve test coverage #69

Conversation

JeffFessler commented Oct 31, 2022

coveralls commented Oct 31, 2022 • edited Loading

Pull Request Test Coverage Report for Build 3377679255

Warning: This coverage report may be inaccurate.

Details

💛 - Coveralls

codecov bot commented Oct 31, 2022 • edited Loading

Codecov Report

JeffFessler commented Nov 1, 2022

Choose a reason for hiding this comment

JeffFessler Nov 2, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mkitti left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

JeffFessler commented Nov 2, 2022

JeffFessler commented Nov 2, 2022

mkitti commented Nov 2, 2022

JeffFessler commented Nov 2, 2022 • edited Loading

mkitti left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mkitti left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mkitti commented Nov 3, 2022

JeffFessler commented Nov 3, 2022

JeffFessler commented Nov 3, 2022

JeffFessler commented Nov 3, 2022

mkitti commented Nov 3, 2022

mkitti left a comment

Choose a reason for hiding this comment

mkitti commented Nov 3, 2022

JeffFessler commented Nov 3, 2022

mkitti commented Nov 3, 2022 • edited Loading

JeffFessler commented Nov 3, 2022

mkitti left a comment

Choose a reason for hiding this comment

JeffFessler commented Nov 3, 2022

coveralls commented Oct 31, 2022 •

edited

Loading

codecov bot commented Oct 31, 2022 •

edited

Loading

JeffFessler Nov 2, 2022 •

edited

Loading

JeffFessler commented Nov 2, 2022 •

edited

Loading

mkitti commented Nov 3, 2022 •

edited

Loading