Add a lazy logrange function and LogRange type

[mcabbott]

I think it would be useful to have an easy way to make a range of points spaced logarithmically. This PR adds the following:

julia> logrange(0.03, 810, 7)  # after some updates
7-element Base.LogRange{Float64}:
 0.03, 0.164317, 0.9, 4.9295, 27.0, 147.885, 810.0

julia> logrange(0.03, 27, 5) |> collect
5-element Vector{Float64}:
  0.03
  0.16431676725154984
  0.9
  4.929503017546495
 27.0

Before Julia 1.0, there was a logspace function which did something different, and less useful:

julia> logspace(1, 4, 4, base=2)
4-element Array{Float64,1}:
  2.0
  4.0
  8.0
 16.0

This was removed in #25896 in favour of just writing 2 .^ (1:4). If you know what powers you want, that is fine. But if you know what start and endpoints you want, then it's more of a hassle to write something like 2 .^ range(log2(2), log2(16), length=4). Well, for powers of 2, I think we can all do that in our heads, but if you want 20 plot points between about 3.5 and 1700, that's when this would be neat.

After deciding whether or not this is a good idea at all:

• I've given this the same notation as Add range(start => stop, length) #39069 but it could equally well be logrange(lo, hi; length) instead, or in addition.
• Matlab's logspace behaves like the old one, so perhaps best to avoid that name. As @mkitti points out below, Python also has similar logspace, and a geomspace which is like this. That's another possible name. Mathematica's PowerRange[a,b,r] takes start, stop, and ratio instead. (This one could be upgraded to accept keywords ratio=1.1 instead of length=10, like range, but I'm not sure that's often useful.)
• Ideally it would use something like TwicePrecision to hit the endpoints more precisely? If it is to be an iterator, then maybe it should iterate by multiplication, not by addition & exp.

Edit:

The current implementation is a LogRange struct. Like LinRange, this should hit the endpoints exactly, but not always rational intermediate points. I think this is a reasonable trade-off, and matches what was asked succinctly here.

This thread got quite long, sorry. There are about a dozen implementations (here and in a linked gist), many trying harder to hit intermediate points, trying to figure out what was possible. Too many examples also show powers of 10, as that's an eay place to see how often you hit the exact number... these are a distraction, I think.

Perhaps I should also stress that the definition of this is mathematically very simple. Just as LinRange(lo, hi, 3) is the two ends plus their arithmetic mean, so LogRange(lo, hi, 3) is the two ends plus the geometric mean. For any quantity which is commonly plotted on a log scale, or measured in dB, this is the natural way to divide up the space between lo and hi. No mention of a base is necessary. No magic inference about whether 0.1 really means 1//10 is needed either.

[antoine-levitt]

I agree this is useful. If we do this we should go with linrange (instead of range(start, stop, length)) for consistency. That might be the best solution here. Added bonus is that we can easily add features to these functions if we need (eg I'd love an endpoints=:none/:left/:right/:both), which is much harder to do with plain old range.

[mkitti]

There probably should be a LogRange object to help with this. As currently implemented, there is some numerical issues for some common bases.

julia> collect( logrange(1 => 1000, 4) )
4-element Array{Float64,1}:
   1.0
   9.999999999999998
  99.99999999999996
 999.9999999999998

May I suggest the following:

       function logrange(p::Pair{<:Real, <:Real}, n::Integer)
           lo, hi = promote(p.first, p.second)
           invert = abs(lo) > abs(hi)
           ratio = invert ? lo/hi : hi/lo
           if ratio % 10 == 0
               exp = exp10
               log = log10
           elseif ratio % 2 == 0
               exp = exp2
               log = log2
           else
               exp = Base.exp
               log = Base.log
           end
           if ratio > 0
               if invert
                   log_ratio = -log(ratio)
                   (x > log_ratio/2 ? lo*exp(x) : hi*exp(x-log_ratio) for x in range(0, log_ratio, length=n))
               else
                   log_ratio = log(ratio)
                   (x < log_ratio/2 ? lo*exp(x) : hi*exp(x-log_ratio) for x in range(0, log_ratio, length=n))
               end
           else
               throw(DomainError(p, "logrange requires that the first and last elements are both positive, or both negative"))
           end
       end

This increases numerical stability for base = 2 and base = 10

julia> collect( logrange(1 => 1000, 4) )
4-element Array{Float64,1}:
    1.0
   10.0
  100.0
 1000.0

julia> collect( logrange(1 => 1024, 3) )
3-element Array{Float64,1}:
    1.0
   32.0
 1024.0

julia> collect( logrange(1 => 1024, 6) )
6-element Array{Float64,1}:
    1.0
    4.0
   16.0
   64.0
  256.0
 1024.0

[mcabbott]

Yes, the implementation here is just a sketch so far. Ideally it would always hit both endpoints to full accuracy, although I'm not sure how hard that would be.

[mkitti]

Yes, the implementation here is just a sketch so far. Ideally it would always hit both endpoints to full accuracy, although I'm not sure how hard that would be.

The trick might be figuring out which endpoint we are closer to and then deciding whether we should count up or down. For start, the calculation should be start*exp(0). For stop, the calculation should be likewise stop*exp(0). We cannot rely on exp(log(x)) === x. However, exp(0.0) is reliably 1.0.

The need for selecting the correct base and for making these kind of decisions is what ultimately justifies having this method, I think.

[mkitti]

Ideally it would always hit both endpoints to full accuracy

I updated my suggested code above. I think a ternary operator takes care of the endpoints.

               log_ratio = log(ratio)
               (x < log_ratio/2 ? lo*exp(x) : hi*exp(x-log_ratio) for x in range(0, log_ratio, length=n))

[mkitti]

Hmm... there is still some asymmetry to this.

julia> collect(logrange(1 => 1000, 4))
4-element Array{Float64,1}:
    1.0
   10.0
  100.0
 1000.0

julia> collect(logrange(1000 => 1, 4))
4-element Array{Float64,1}:
 999.9999999999998
 100.00000000000004
  10.000000000000005
   1.0000000000000002

[mkitti]

I fixed it by always taking the larger ratio.

                  invert = abs(lo) > abs(hi)
                  ratio = invert ? lo/hi : hi/lo

So now everything seems to work well:

julia> collect( logrange( 1000 => 1, 4) )
4-element Array{Float64,1}:
 1000.0
  100.0
   10.0
    1.0

julia> collect( logrange( 1 => 1000, 4) )
4-element Array{Float64,1}:
    1.0
   10.0
  100.0
 1000.0

julia> collect( logrange( -1 => -1000, 4) )
4-element Array{Float64,1}:
    -1.0
   -10.0
  -100.0
 -1000.0

julia> collect( logrange( -1000 => -1, 4) )
4-element Array{Float64,1}:
 -1000.0
  -100.0
   -10.0
    -1.0

[mkitti]

The Numpy equivalent of this function is called numpy.geomspace https://numpy.org/doc/stable/reference/generated/numpy.geomspace.html

[darsnack]

Can I suggest that 0 be allowed as an endpoint? I've always found it frustrating when I write code sweeping a parameter and xs = linrange(0 => 1, 10) works but logrange(0 => 1, 10) errors (here I'm just using the proposed names for clarity). Programmatically, this isn't surprising, but if in English I say "give me 4 points between 0 and 1 logarithmically spaced apart" then I expect [0, 0.01, 0.1, 1.0] to be the result. It would be nice if I didn't have to write my own special handling to achieve this.

[antoine-levitt]

@darsnack what you are suggesting requires an explicit base, which logrange doesn't because it's base-independent (the notion of "logarithmically equispaced" does not depend on which logarithm you choose)

[darsnack]

Ah thanks for the explanation, makes sense!

[mcabbott]

If this takes keywords like range does, then something like logspace(0 => 1, length=4, ratio=10) == [0.001, 0.01, 0.1, 1.0] might not be too odd. You need to provide three of start/stop/length/ratio, and 0 is acting as some kind of nothing? The same could be true of logspace(1 => Inf, length=4, ratio=10) == [1, 10, 100, 1000].

Including zero in the results [0, 0.01, 0.1, 1.0] would be much more surprising, as the ratio is infinite for that step.

[darsnack]

Yeah I don't necessarily want to complicate the proposal with too many API options now that I understand the rationale.

I guess it depends on which arguments are chosen for specification. To me, whenever I specify start and stop, I expect them to be in the range no matter what. So even though the ratio is infinite, the inclusion of 0.0 is less surprising than a changing lower bound depending on the number of steps.

[mkitti]

The latest changes look very nice. I see that the iteration is based on multiplying the ratio by the prior element. How confident are you that you can avoid accumulated errors through long iterations?

The iteration for UnitRange and StepRangeLen uses multiplication rather than addition:
https://github.com/JuliaLang/julia/blob/master/base/range.jl#L664

StepRangeLen:
https://github.com/JuliaLang/julia/blob/master/base/range.jl#L720-L723

LinRange:
https://github.com/JuliaLang/julia/blob/master/base/range.jl#L730-L738

Only for OrdinalRange where the step is an Integer is the iteration done by addition:
https://github.com/JuliaLang/julia/blob/master/base/range.jl#L669-L674

In this case, might there be situations where multiplication or exponentiation might result in better stability?

[mcabbott]

All else equal it would be nicer to have random access, which is what multiplication buys you in linear ranges. But it's tricky to make that fast & accurate here, I think the answer may simply be to collect this iterator if you need that.

Notice that it iterates in twice the precision of the output, so there are quite a few spare bits for the inevitable errors. I haven't been very systematic about stress-testing it, I don't promise it's < eps everywhere. It's noticeably less accurate for complex ranges, and possibly that's a mistake like missing a conj somewhere; the fact that it accepts complex numbers at all seems like a generalisation which happens for free, so that's OK I think.

This gist compares some different approaches: https://gist.github.com/mcabbott/a40c784b4fd27656b581ad0dac09aae5 . (And, I will have to look into why this doesn't build, a bit confused about imports to Iterators module but can surely sort it out.)

[JeffBezanson]

I think we should not use Pair for this; I don't believe any other similar functions work like that?

[mcabbott]

Pair was intended to match #39069 for linear range, if that got merged.

Notation is now logrange(a, b; length=10), and also allows logrange(a, b; ratio=2). It does not allow all possible combinations, for instance logrange(a; length=10, ratio=2) is currently spelled LogRange(a, 2, 10) instead.

[mkitti]

To fix the build, you might need to split the Complex code off and rename some functions so you don't have to duplicate them.

I also suggest that you take a look StepRangeLen which has a flexible reference value to increase accuracy in a particular part of the range.

[mbauman]

Can't this just be a very simple wrapper to whatever range spits back? It feels a little unnecessary to completely re-implement the wheel here for what's effectively 10 .^ (range(log(10, start), log(10, stop), length=len))

That is, simply:

julia> struct LogRange{T, R<:AbstractRange} <: AbstractVector{T}
           r::R
       end

julia> LogRange(r) = LogRange{typeof(exp(first(r))), typeof(r)}(r)
LogRange

julia> Base.size(r::LogRange) = size(r.r)

julia> Base.getindex(r::LogRange, i::Integer) = exp(r.r[i])

julia> logrange(start, stop, len) = LogRange(range(log(start), log(stop), length=len))
logrange (generic function with 1 method)

I guess there's the floating point annoyance that we no longer hit endpoints exactly due to the subsequent floating point exponentiation... that's annoying.

julia> logrange(0.0001, 1, 5)
5-element LogRange{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}:
 0.00010000000000000009
 0.0010000000000000002
 0.010000000000000004
 0.10000000000000002
 1.0

[mcabbott]

Yes my first version was literally (exp(x) for x in range(log(lo), log(hi), length=n)), all the present complication is in order to hit pretty numbers exactly. This seems easier to do with a generator than an AbstractArray, because there's no exp(::TwicePrecision). The gist above has the things I tried, and the fastest way to get something both accurate and indexable seemed to be collecting this generator (assuming you do visit all the values).

[mbauman]

IMO for this to be considered — and for it to be called logrange — it needs to be a lazy and efficient AbstractArray and it needs to hit its endpoints exactly. I don't think it needs do nice things with the intermediate values, but I fear a very long tail of floating point issues (a la the existing ranges themselves).

[mcabbott]

I will think some more about strategies to do that, when I get a chance.

Another option to perhaps consider is for it to be Iterators.logrange, so that there is no expectation of random access.

[oscardssmith]

Do we want exp(::TwicePrecision)? If so, I could get a PR with a reasonably fast (not fully optimized) version that would be 5 or so lines and within 4x of optimal performance.

[mcabbott]

Would this probably have similar accuracy to DoubleFloats? If so that may solve the problem, but perhaps I should experiment a bit more with (& without) that to see what works well.

julia> using DoubleFloats

julia> [Float64(exp(x)) for x in range(log(Double64(0.1)), log(Double64(100.0)), length=7)]
7-element Vector{Float64}:
   0.1
   0.31622776601683794
   1.0
   3.1622776601683795
  10.0
  31.622776601683793
 100.0

[oscardssmith]

Yeah DoubleFloats and TwicePrecision are basically the same TwicePrecision is just the version in Base that's really hacky and designed only for internal use.

[mkitti]

IMO for this to be considered — and for it to be called logrange — it needs to be a lazy and efficient AbstractArray and it needs to hit its endpoints exactly. I don't think it needs do nice things with the intermediate values, but I fear a very long tail of floating point issues (a la the existing ranges themselves).

Hitting the endpoints is easy if we just save those values in the struct. The most noticeable problem is when the roots are well known.

I guess there's the floating point annoyance that we no longer hit endpoints exactly due to the subsequent floating point exponentiation... that's annoying.

If @mbauman used used exp10 and log10 for those specific numbers, the results will meet expectations:

julia> struct LogRange{T, R<:AbstractRange} <: AbstractVector{T}
                  r::R
              end

julia> LogRange(r) = LogRange{typeof(exp(first(r))), typeof(r)}(r)
LogRange

julia> Base.size(r::LogRange) = size(r.r)

julia> Base.getindex(r::LogRange, i::Integer) = exp(r.r[i])

julia> Base.getindex(r::LogRange, i::Integer) = exp10(r.r[i])

julia> logrange(start, stop, len) = LogRange(range(log10(start), log10(stop), length=len))
logrange (generic function with 1 method)

julia> logrange( 0.0001, 1, 5 )
5-element LogRange{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}:
 0.0001
 0.001
 0.01
 0.1
 1.0

julia> logrange( 0.0001, 10000, 9 )
9-element LogRange{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}:
     0.0001
     0.001
     0.01
     0.1
     1.0
    10.0
   100.0
  1000.0
 10000.0

Of course this fails for base 2.

julia> logrange( 2^-5, 2^5, 11 )
11-element LogRange{Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}:
  0.031249999999999997
  0.0625
  0.12499999999999999
  0.25
  0.5
  1.0
  2.0
  4.0
  8.000000000000002
 16.0
 32.00000000000001

There's exp2 and log2 for that purpose. A useful function might figure out if a rational base exists and take advantage of it.

[mcabbott]

That's not possible with the current logrange structure

This is false:

julia> @eval Base begin  # on 623a4f8
         function *(λ::Number, r::LogRange)
           start = λ * r.start
           stop = λ * r.stop
           len = r.len
           extra = log(λ)/(r.len-1) .+ r.extra  # should be more careful about promotion, high precision, negative
           $(Expr(:new, :(LogRange{typeof(start), typeof(extra[1])}), :start, :stop, :len, :extra))
         end
         show(io::IO, r::LogRange{T}) where {T<:Complex} = invoke(show, Tuple{IO,AbstractVector{T}}, io, r)
       end;

julia> lo * logrange(1, hi/lo, 5)
5-element Base.LogRange{ComplexF32, ComplexF64}:
 -1.0+0.01im, -1.00004+0.00500006im, -1.00005-1.43682f-9im, -1.00004-0.00500006im, -1.0-0.01im

Edit, of course @eval is a hack, in reality we are free to remove the inner constructor. So IDK what "internals" is supposed to mean.

Obviously the entire point of using * to construct something following different branch conventions is that this doesn't follow the same branch conventions as the default constructor.

[mcabbott]

Given that it is less commonly useful and involves making more choices (e.g. branch cut, whether to support multiplication by scalar) I think it would be a good idea to merge this without complex support

You do obviously need to choose a branch cut, and 4 possible conventions were calmly listed above. I still think that a natural one falls out of log, and that it's a bit of a shame not to buid generic things out of generic pieces. But I give up.

The last commit removes all complex support, aiming for the minimal thing.

One unfortunate side effect of being less generic is that the struct can no longer store a Quantity from Uniful.jl, as this is <:Number but not Real.

Whether to support some operations which return another LogRange (a bit like 1 .+ range(2,3,4) making another range) is an entirely orthogonal question. This PR does not add any such things.

[aplavin]

That's not possible with the current logrange structure

This is false: <...>

Thanks for pointing out the inaccuracy, I should've said "not possible using actual api, without resorting to julia & logrange internals". Sorry for being too lazy in stating that :)

[LilithHafner]

@oscardssmith, @mcabbott, is this ready to merge? I don't see any outstanding objections.

I still think that a natural one falls out of log, and that it's a bit of a shame not to buid generic things out of generic pieces.

We still can, I just think it will be easier to do in a separate PR focused on expanding the scope of LogRange and don't want to delay merging what we have now until we have consensus for Complex handling.

[mcabbott]

One small question is whether we are happy with this asymmetry:

julia> logrange(2, 1e300, 5)
5-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:
 2.0, 1.68179e75, 1.41421e150, 1.18921e225, 1.0e300

julia> logrange(2, Inf, 5)  # matches the limit
5-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:
 2.0, Inf, Inf, Inf, Inf

julia> logrange(1e-300, 20, 5)
5-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:
 1.0e-300, 2.11474e-225, 4.47214e-150, 9.45742e-75, 20.0

julia> logrange(0, 20, 5)  # limit would be [0,0,0,0,20]
5-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:
 NaN, NaN, NaN, NaN, 20.0

My take is that neither of these is useful, but forgetting that 0 is infinitely far away is pretty common, and if you silently get something == [0,0,0,0,20] it's going to take much longer to debug than getting NaN. But returning [0,NaN,NaN,NaN,20] is another possibility.

Supplying Inf is obviously a problem, and getting Inf back poisons things downstream about as well as NaN does. But [2,NaN,NaN,NaN,Inf] would be a more symmetrical choice.

[adienes]

well the (current) docstring does say

the start and stop arguments are always the first and last elements of the result

so maybe [0,NaN,NaN,NaN,20] is better than [NaN, NaN, NaN, NaN, 20]

however I also think it would be reasonable to just throw if start=0, in

elseif start < 0 || stop < 0
            # log would throw, but _log_twice64_unchecked does not
            throw(DomainError((start, stop),
                "LogRange(start, stop, length) does not accept negative numbers"))

[LilithHafner]

I don't like that asymmetry.

I think the most mathematically sound definition, and my preference, would be to take the limit, follow
sqrt(0*10) === 0.0, and get

logrange(0, Inf, 5) = [0.0, NaN, NaN, NaN, Inf]
logrange(0, 10, 5) = [0.0, 0.0, 0.0, 0.0, 10.0]
logrange(10, Inf, 5) = [10.0, Inf, Inf, Inf, Inf]

OTOH, I suspect we'll get bug reports of folks saying that logrange(0,10,10) does not provide a uniformly spaced output if we do that.

Failing faster with

logrange(0, Inf, 5) = [0.0, NaN, NaN, NaN, Inf]
logrange(0, 10, 5) = [0.0, NaN, NaN, NaN, 10.0]
logrange(10, Inf, 5) = [10.0, NaN, NaN, NaN, Inf]

or even

logrange(0, Inf, 5) = [NaN, NaN, NaN, NaN, NaN]
logrange(0, 10, 5) = [NaN, NaN, NaN, NaN, NaN]
logrange(10, Inf, 5) = [NaN, NaN, NaN, NaN, NaN]

is plausible, depending on how much we want to be proactively helpful. I tend to prefer returning sensible results whenever possible (i.e. not eagerly jumping to NaN if a limiting value exists) and counting on the user to provide 0/Inf/Nan inputs only if they want the appropriate outputs because that allows slightly more general usage, but could see going either way.

Setting the halfway point sqrt(0*10) to NaN and the half way point sqrt(10*Inf) to Inf is clever, and your justification based on plausible usefulness of the result and plausible reasons for the input make sense. However, in this case I think we should strive to be sound and predictable at the expense of being less clever and less unexpectedly helpful.

Edit after reading @adienes's comment: if we are going to proactively assume the user made a mistake and return a result other than the limit, I think we should do it by throwing, not by returning NaN.

[mcabbott]

My argument for not throwing an error is that log(0), log(Inf) do not.

Can I use log(0) being infinite to argue it's OK to "poison" everything after it, and hence logrange(0, 10, 5) = [0.0, NaN, NaN, NaN, 10.0]? Obeying "always the first and last elements of the result" is nice. Perhaps it's OK that logrange(0, 10, 2) == [0, 10] isn't poisoned.

Looking for precedents:

• Note that range(1, Inf, 3) and range(1,NaN,3) are both errors. But Inf * range(1,2,3) is somewhat surprisingly a StepRangeLen of all NaN.

• LinRange(1, Inf, 3) gives NaN instead. This is also true of LinRange(1,Inf,2) = [NaN, NaN] which means it does not always return given first/last. Maybe that's a bug, although it makes no explicit claim that it will do so -- it's just an observation that its never makes floating point errors at the ends. (Edit: I don't see any explicit tests of how it handles NaN, so perhaps this isn't a design choice. I think there are no branches in the linear interpolation Base.lerpi.)

[mcabbott]

Last commit just throws errors on 0, NaN, Inf, matching what range does. (Not so sure which should be DomainError vs an ArgumentError, does not match range on this.)

[LilithHafner]

Looks mergable to me. I've reviewed everything but the src of _log_twice64_unchecked and the usage of math internals, which I'm counting on @oscardssmith's previous approval to cover.