Fix #39798 by using oneunit(start) rather than 1 for _linrange

[mkitti]

Fixes #39798 by using oneunit.

julia> range(Int32(1), Int32(1), 10)
1.0:0.0:1.0

julia> range(Int32(1), Int32(5), 10)
1.0:0.4444444444444444:5.0

[mkitti]

Should this just be one since start is an Integer?

[timholy]

oneunit seems like the right answer. one is kind of ill-defined, whereas oneunit is straightforward and clear.

[mkitti]

I added the following test set:

@testset "Non-Int64 endpoints that are identical (#39798)" begin
    for T in DataType[Float64,Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64], 
        r in [ LinRange(1, 1, 10), StepRangeLen(7, 0 , 5) ]
        let start=T(first(r)), stop=T(last(r)), step=T(step(r)), length=length(r)
            @test range(  start, stop,       length) == r
            @test range(  start, stop;       length) == r
            @test range(  start; stop,       length) == r
            @test range(; start, stop,       length) == r
            @test range(  start;       step, length) == r
            @test range(; start,       step, length) == r
            @test range(;        stop, step, length) == r
        end
    end
end

On 1.5.3, it fails:

Test Summary:                                   | Pass  Error  Total
Non-Int64 endpoints that are identical (#39798) |   10    130    140

On this branch, it succeeds:

Test Summary:                                   | Pass  Total
Non-Int64 endpoints that are identical (#39798) |  140    140

[mbauman]

I think we want one, since this is used as the denominator. Using oneunit would strip units from a Unitful range, no? Even if unitful linspaces don't work or don't hit this codepath, it seems like decent motivation to use one instead.

[mkitti]

one vs oneunit

Regarding one vs oneunit, at the moment it probably does not matter with current code since it is being applied to start which is an Integer. For any concrete Integer subtype that I can find one(start) === oneunit(start). I would love to see a counter example.

Thus, the choice has to be informed by method dispatch. Let's examine the stack trace and see where this fails.

Stack Trace

 [1] Base.TwicePrecision{Float64}(::Tuple{Int32,Int64}, ::Int64) at ./twiceprecision.jl:185
 [2] steprangelen_hp(::Type{Float64}, ::Tuple{Int32,Int64}, ::Tuple{Int32,Int64}, ::Int64, ::Int64, ::Int64) at ./twiceprecision.jl:323
 [3] _linspace(::Type{Float64}, ::Int32, ::Int32, ::Int64, ::Int64) at ./twiceprecision.jl:666
 [4] _linspace at ./twiceprecision.jl:663 [inlined]

1. The first line generating the error is actually the Base.TwicePrecision default inner constructor which is expecting two values of Float64. However, the arguments given are Tuple{Int32,Int64} and Int64. Clearly, the caller did not intend to invoke the inner constructor of Base.TwicePrecision but another method. The intended method signature to call is Base.TwicePrecision{T}(nd::Tuple{I,I}, nb::Integer) where {T,I}. Both elements of the nd Tuple must be of the same parameteric type, I.
   

   julia/base/twiceprecision.jl

   Lines 213 to 215 in 788b2c7

   	function TwicePrecision{T}(nd::Tuple{I,I}, nb::Integer) where {T,I}
   	twiceprecision(TwicePrecision{T}(nd), nb)
   	end

2. TwicePrecision{Float64}(step, nb) in steprangelen_hp is how Base.TwicePrecision is invoked. step is the third argument of steprangelen_hp
   

   julia/base/twiceprecision.jl

   Lines 320 to 325 in 788b2c7

   	function steprangelen_hp(::Type{Float64}, ref::Tuple{Integer,Integer},
   	step::Tuple{Integer,Integer}, nb::Integer,
   	len::Integer, offset::Integer)
   	StepRangeLen(TwicePrecision{Float64}(ref),
   	TwicePrecision{Float64}(step, nb), Int(len), offset)
   	end

3. The call to steprangelen_hp is steprangelen_hp(T, (start_n, den), (zero(start_n), den), 0, len, 1) in _linspace where the third argument is (zero(start_n), den)
   

   julia/base/twiceprecision.jl

   Lines 664 to 667 in 788b2c7

   	function _linspace(::Type{T}, start_n::Integer, stop_n::Integer, len::Integer, den::Integer) where T<:IEEEFloat
   	len < 2 && return _linspace1(T, start_n/den, stop_n/den, len)
   	start_n == stop_n && return steprangelen_hp(T, (start_n, den), (zero(start_n), den), 0, len, 1)
   	tmin = -start_n/(Float64(stop_n) - Float64(start_n))

The intention is for TwicePrecision{Float64}(step, nb) === TwicePrecision{Float64}(0.0, 0.0). To dispatch correctly with the current methods, zero(start_n) and den must be of the same type.

Analysis

I would think that one should be the analog to zero. However, this is not really the case if there were actually some hypothetical concrete subtype of Integer.

julia> zero(Dates.Day(1))
0 days

julia> one(Dates.Day(1))
1

julia> oneunit(Dates.Day(1))
1 day

julia> using Unitful

julia> zero(1u"g")
0 g

julia> one(1u"g")
1

julia> oneunit(1u"g")
1 g

Empirically, oneunit rather than one appears to produce the same type as zero. Therefore, given that zero(start_n) and den must be of the same type, I think that @timholy is correct in that oneunit(start_n) more consistently does that.

Regarding the matter of stripping units as brought up by @mbauman , this code path seems to clearly intended for unitless ranges of eltype determined by the parameter T and T <: IEEEFloat. All members of that Union (Float16,Float32, and Float64) are unitless. Since the values zero(start_n) and den form a ratio to determine a Base.TwicePrecision{T<: IEEEFloat}, stripping units is exactly what what we would want to do.

Unitful.jl creates unit based ranges by stripping the units, creating a range, and then adding them back in.
https://github.com/PainterQubits/Unitful.jl/blob/f673d2db96b6c86736849209d922a726be420721/src/range.jl#L19-L22

Alternative Approach

At the end of the day, this is all trying to return StepRangeLen( TwicePrecision{T}(start_n), TwicePrecision{T}(0,0), len) which is practically equivalent to StepRangeLen(start_n, 0 len). Maybe we should just return that from _linrange rather than delegating to steprangelen_hp?

This alternative approach could be pursued in a separate pull request in addition to the current fix which is still applicable even if a StepRangeLen is returned directly.

Conclusion

Based on the current flow of the code in terms of the dispatched method, oneunit(start) seems to match zero(start) in that they consistently produce the same type. The AbstractRange being produced by this code has an element type of IEEEFloat and is thus unitless.

Thus, we should proceed with merging the current approach.

[mkitti]

Can we add a bugfix label here?

[mkitti]

@mbauman Any further thoughts on this? Do you still still think one would be better than oneunit?

[mbauman]

It's not often that I disagree with Tim, but yes, I still think that one(T) is the better choice given that it becomes the den in _linspace1(T, start_n/den, stop_n/den, len).

The question isn't so much if I still think so, but rather if I managed to convince @timholy. :)

[mbauman]

If not, what you have is just fine; this really is a distinction without a (currently-observable) difference.

[mkitti]

It's not often that I disagree with Tim, but yes, I still think that one(T) is the better choice given that it becomes the den in _linspace1(T, start_n/den, stop_n/den, len).

The question isn't so much if I still think so, but rather if I managed to convince @timholy. :)

I think we may have made this too complicated. Ultimately, we are talking about oneunit(start) vs one(start) where typeof(start) <: Integer. Really then either oneunit(start) and one(start) would work in this case for all concrete subtypes of Integer of which I am aware.

[mkitti]

The dispatch analysis I did above suggests that the intended method to invoke is TwicePrecision{T}(nd::Tuple{I,I}, nb::Integer) where {T,I} and we are trying to invoke TwicePrecision{Float64}(step, nb) where step is (zero(start_n), den).

Thus typeof(zero(start_n)) must equal typeof(den). In the hypothetical scenario where a Unitful type is a subtype of Integer that only happens when using oneunit:

julia> using Unitful

julia> start_n = 1u"g"
1 g

julia> zero(start_n)
0 g

julia> one(start_n)
1

julia> oneunit(zero(start_n))
1 g

In this purely hypothetical scenario, if we used one(start) then we would also need to implement TwicePrecision{T}(nd::Tuple{I,J}, nb::Integer) where {T,I,J} to accommodate the call to TwicePrecision{Float64}( (zero(start_n), one(start_n)) , nb).

The problem with one(start) is that we don't have a TwicePrecision constructor where the numerator and denominator are of different types.

[mkitti]

An alternative fix (or in combination with one(start) is just implement TwicePrecision{T}(nd::Tuple{Integer,Integer}, nb::Integer) where {T}

julia> range(Int32(1), Int32(1), length = 10)
ERROR: MethodError: Cannot `convert` an object of type Tuple{Int32, Int64} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::T) where T<:Number at number.jl:6
  convert(::Type{T}, ::Number) where T<:Number at number.jl:7
  convert(::Type{T}, ::TwicePrecision) where T<:Number at twiceprecision.jl:250
  ...
Stacktrace:
 [1] TwicePrecision{Float64}(hi::Tuple{Int32, Int64}, lo::Int64)
   @ Base .\twiceprecision.jl:185
 [2] steprangelen_hp(#unused#::Type{Float64}, ref::Tuple{Int32, Int64}, step::Tuple{Int32, Int64}, nb::Int64, len::Int64, offset::Int64)
   @ Base .\twiceprecision.jl:323
 [3] _linspace(#unused#::Type{Float64}, start_n::Int32, stop_n::Int32, len::Int64, den::Int64)
   @ Base .\twiceprecision.jl:663
 [4] _linspace
   @ .\twiceprecision.jl:660 [inlined]
 [5] _range
   @ .\range.jl:441 [inlined]
 [6] _range2
   @ .\range.jl:100 [inlined]
 [7] #range#58
   @ .\range.jl:94 [inlined]
 [8] top-level scope
   @ REPL[2]:1

julia> import Base: TwicePrecision, twiceprecision

julia> function TwicePrecision{T}(nd::Tuple{Integer,Integer}, nb::Integer) where {T}
           twiceprecision(TwicePrecision{T}(nd), nb)
       end

julia> range(Int32(1), Int32(1), length = 10)
1.0:0.0:1.0

Anyways, I think it would be great for range(Int32(1), Int32(1), length = 10) to work.

[mkitti]

I implemented the one(start) variation in #40467 which also defines TwicePrecision{T}(nd::Tuple{Integer,Integer}, nb::Integer) where {T}

[mkitti]

@timholy Do you have any further thoughts on this fix or the alternate in #40467? Perhaps we can resolve this issue for Julia 1.7.