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

Fix ustrip broadcasting when range step has different unit than eltype #715

Merged
merged 5 commits into from
May 13, 2024

Conversation

sostock
Copy link
Collaborator

@sostock sostock commented Apr 20, 2024

Fixes #712.

This PR also changes the return type when broadcasting over a StepRange and the conversion factor is a floating-point number: in this case, a StepRangeLen is returned (instead of a StepRange) because it is the better choice for floating-point numbers.

src/range.jl Outdated
broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::AbstractRange, x::Ref{<:Units}) = r * x[]
broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Ref{<:Units}, r::AbstractRange) = x[] * r
broadcasted(::DefaultArrayStyle{1}, x::BCAST_PROPAGATE_CALLS, r::StepRangeLen) = StepRangeLen{typeof(x(zero(eltype(r))))}(x(r.ref), x(r.step), r.len, r.offset)
broadcasted(::DefaultArrayStyle{1}, x::BCAST_PROPAGATE_CALLS, r::StepRange) = StepRange(x(r.start), x(r.step), x(r.stop))
broadcasted(::DefaultArrayStyle{1}, x::BCAST_PROPAGATE_CALLS, r::StepRange) = x(r.start):x(r.step):x(r.stop)
Copy link
Collaborator Author

@sostock sostock Apr 29, 2024

Choose a reason for hiding this comment

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

This is bad for performance:

julia> using BenchmarkTools

julia> using Unitful: mJ, eV

julia> r = 1eV:1eV:5eV
(1:5) eV

julia> @btime mJ.($r);
  25.844 ns (0 allocations: 0 bytes) # v1.19.0 (returns a StepRange)
  55.767 ns (0 allocations: 0 bytes) # this PR (returns a StepRangeLen)

We could add methods that call the StepRangeLen constructor directly instead of going through start:step:stop:

julia> simple_mJ(r::StepRange) = StepRangeLen{typeof(mJ(r.start))}(Base.TwicePrecision(mJ(r.start)), Base.TwicePrecision(mJ(r.step)), length(r));

julia> @btime _mJ($r)
  14.013 ns (0 allocations: 0 bytes)

julia> mJ.(r) === simple_mJ(r)
true

This is a simple implementation (maybe the accuracy can be improved in certain cases by setting the appropriate number of trailing zero bits of step.hi and setting the offset so that the reference value is the one with the smallest magnitude), but, at least in the following case, it is actually more accurate than going through start:step:stop:

julia> r2 = 1eV:1eV:500_000eV
(1:500000) eV

julia> last(mJ.(r2))
8.01086714823366e-11 mJ

julia> last(simple_mJ(r2))
8.01088317e-11 mJ

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I have now rewritten the function. It is slightly slower than the old version that returned a StepRange, but should be more accurate than the simple version in the previous comment:

julia> @btime mJ.($r)
  25.844 ns (0 allocations: 0 bytes) # v1.19.0 (returns a StepRange)
  33.379 ns (0 allocations: 0 bytes) # this PR (returns a StepRangeLen)

@sostock sostock merged commit 9aa8703 into PainterQubits:master May 13, 2024
13 checks passed
@sostock sostock deleted the fix_ustrip_broadcasting branch May 13, 2024 07:40
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.

Wrong result when broadcasting ustrip over ranges whose eltype and step have different units
1 participant