-
-
Notifications
You must be signed in to change notification settings - Fork 35
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
Mat index bug #249
base: main
Are you sure you want to change the base?
Mat index bug #249
Conversation
Fixes a couple of the formerly-broken tests!
See discussion in this issue: #248 |
@jonniedie polite ping re: this PR. |
test/runtests.jl
Outdated
@@ -132,7 +137,7 @@ end | |||
for T in [Int64, Int32, Float64, Float32, ComplexF64, ComplexF32] | |||
@test ComponentArray(a = T[]) == ComponentVector{T}(a = T[]) | |||
@test ComponentArray(a = T[], b = T[]) == ComponentVector{T}(a = T[], b = T[]) | |||
@test ComponentArray(a = T[], b = (;)) == ComponentVector{T}(a = T[], b = T[]) | |||
@test_broken ComponentArray(a = T[], b = (;)) == ComponentVector{T}(a = T[], b = T[]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's important we don't break this one because the behavior is relied upon in some some large simulation projects. This allows ComponentArray
s that match a nested model structure to be initialized when one of the internal models doesn't have integrable state.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Gotcha, thanks. I'll take a look at it again when I get a chance.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@jonniedie Another polite ping. :-)
@jonniedie Another polite ping re: this PR. :-) It's ready to go from my perspective. |
Codecov ReportAttention: Patch coverage is
❗ Your organization needs to install the Codecov GitHub app to enable full functionality. Additional details and impacted files@@ Coverage Diff @@
## main #249 +/- ##
==========================================
+ Coverage 73.10% 73.65% +0.55%
==========================================
Files 25 25
Lines 777 801 +24
==========================================
+ Hits 568 590 +22
- Misses 209 211 +2 ☔ View full report in Codecov by Sentry. |
@jonniedie Another polite ping re: this PR. :-) |
I merged in |
|
||
## Test setup | ||
c = (a = (a = 1, b = [1.0, 4.4]), b = [0.4, 2, 1, 45]) | ||
nt = (a = 100, b = [4, 1.3], c = c) | ||
nt2 = (a = 5, b = [(a = (a = 20, b = 1), b = 0), (a = (a = 33, b = 1), b = 0)], c = (a = (a = 2, b = [1, 2]), b = [1.0 2.0; 5 6])) | ||
|
||
ax = Axis(a = 1, b = 2:3, c = ViewAxis(4:10, (a = ViewAxis(1:3, (a = 1, b = 2:3)), b = 4:7))) | ||
ax_c = (a = ViewAxis(1:3, (a = 1, b = 2:3)), b = 4:7) | ||
ax = Axis(a = 1, b = r2v(2:3), c = ViewAxis(4:10, (a = ViewAxis(1:3, (a = 1, b = r2v(2:3))), b = r2v(4:7)))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hm. I guess I don't understand why we'd need to wrap this in a new type. In the example you gave in the issue, b
is a vector element, not an nx1
matrix. I don't think we should introduce a new ShapedAxis1d
type if all vector elements are behaving incorrectly. We should just fix the issue of adjoint/transposition that's broken directly. So I guess specifically, the ShapedAxis
needs to be created from the FlatAxis
during the adjoint operation, not beforehand. And for that, we can use the normal ShapedAxis
without having to introduce a new type.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@jonniedie My bad, just saw this comment.
Hmm... not sure I follow, but fwiw I think everything is fine with the transpose—the issue comes up later on.
The problem with the test case in issue #249, aka
ni = 4
nj = 2
nk = 3
X = ComponentVector(a=0.0, b=zeros(Float64, nk), c=zeros(Float64, ni, nj))
Y = ComponentVector(d=zeros(Float64, ni, nj))
J = Y.*X'
appears to be due to the fact that the shape of a unit range (which is used for vector components like X[:b]
in the example) isn't understood by ComponentArrays.maybe_reshape
. The error I get:
julia> J[:d, :b]
ERROR: DimensionMismatch: new dimensions (4, 2) must be consistent with array size 24
Stacktrace:
[1] (::Base.var"#throw_dmrsa#328")(dims::Tuple{Int64, Int64}, len::Int64)
@ Base ./reshapedarray.jl:41
[2] reshape
@ ./reshapedarray.jl:45 [inlined]
[3] maybe_reshape
@ ~/projects/componentarrays_bugs/dev/ComponentArrays/src/componentarray.jl:237 [inlined]
[4] ComponentArray
@ ~/projects/componentarrays_bugs/dev/ComponentArrays/src/componentarray.jl:52 [inlined]
[5] macro expansion
@ ~/projects/componentarrays_bugs/dev/ComponentArrays/src/array_interface.jl:0 [inlined]
[6] _getindex(::typeof(getindex), ::ComponentMatrix{Float64, Matrix{Float64}, Tuple{Axis{…}, Axis{…}}}, ::Val{:d}, ::Val{:b})
@ ComponentArrays ~/projects/componentarrays_bugs/dev/ComponentArrays/src/array_interface.jl:119
[7] getindex
@ ~/projects/componentarrays_bugs/dev/ComponentArrays/src/array_interface.jl:103 [inlined]
[8] getindex(::ComponentMatrix{Float64, Matrix{Float64}, Tuple{Axis{…}, Axis{…}}}, ::Symbol, ::Symbol)
@ ComponentArrays ~/projects/componentarrays_bugs/dev/ComponentArrays/src/array_interface.jl:102
[9] top-level scope
@ REPL[16]:1
Some type information was truncated. Use `show(err)` to see complete types.
shell>
In maybe_reshape
:
# Reshape ComponentArrays with ShapedAxis axes
maybe_reshape(data, ::NotShapedOrPartitionedAxis...) = data
function maybe_reshape(data, axs::AbstractAxis...)
shapes = filter_by_type(ShapedAxis, axs...) .|> size
shapes = reduce((tup, s) -> (tup..., s...), shapes)
return reshape(data, shapes)
end
during J[:d, :b]
, axs
is
axs = (ShapedAxis((4, 2)), FlatAxis())
i.e. the unit range is converted to a FlatAxis
, which is then ignored when defining shapes
in maybe_reshape
.
The axes of the transpose X'
look fine to me:
julia> getaxes(X')
(FlatAxis(), Axis(a = 1, b = 2:4, c = ViewAxis(5:12, ShapedAxis((4, 2)))))
julia>
The point of the Shaped1DAxis
is to have a type that can be a part of the NotShapedOrPartitionedAxis
Union
defined in axis.jl
.
This allows it to skip the problematic and unnecessary maybe_reshape
for vector components.
If there was a way to dispatch on ShapedAxis{Shape}
when Shape
was a length-1 Tuple
then I don't think we'd need Shaped1DAxis
, but I don't know of a way to do that.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@jonniedie Another polite ping re: this PR. :-)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@jonniedie Another polite ping re: this PR.
No description provided.