From 220f4d7d34e9de9c14fece4b43425565713f2b25 Mon Sep 17 00:00:00 2001 From: Michael Abbott <32575566+mcabbott@users.noreply.github.com> Date: Sat, 24 Feb 2024 08:26:47 -0500 Subject: [PATCH] logrange (#821) --- Project.toml | 2 +- README.md | 3 + src/Compat.jl | 240 +++++++++++++++++++++++++++++++++++++++++++++-- test/runtests.jl | 65 +++++++++++++ 4 files changed, 301 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index b10bf9e03..fd1d17c61 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Compat" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.13.0" +version = "4.14.0" [deps] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" diff --git a/README.md b/README.md index 1d645983d..b935b68f9 100644 --- a/README.md +++ b/README.md @@ -72,6 +72,8 @@ changes in `julia`. * `allequal(f, itr)` and `allunique(f, itr)` methods. ([#47679]) (since Compat 4.13.0) +* `logrange(lo, hi; length)` is like `range` but with a constant ratio, not difference. ([#39071]) (since Compat 4.14.0) Note that on Julia 1.8 and earlier, the version from Compat has slightly lower floating-point accuracy than the one in Base (Julia 1.11 and later). + * `Iterators.cycle(itr, n)` is the lazy version of `repeat(vector, n)`. ([#47354]) (since Compat 4.13.0) * `@compat public foo, bar` marks `foo` and `bar` as public in Julia 1.11+ and is a no-op in Julia 1.10 and earlier. ([#50105]) (since Compat 3.47.0, 4.10.0) @@ -163,6 +165,7 @@ Note that you should specify the correct minimum version for `Compat` in the [#36229]: https://github.com/JuliaLang/julia/issues/36229 [#37978]: https://github.com/JuliaLang/julia/issues/37978 [#39037]: https://github.com/JuliaLang/julia/issues/39037 +[#39071]: https://github.com/JuliaLang/julia/pull/39071 [#39245]: https://github.com/JuliaLang/julia/issues/39245 [#39285]: https://github.com/JuliaLang/julia/issues/39285 [#39794]: https://github.com/JuliaLang/julia/issues/39794 diff --git a/src/Compat.jl b/src/Compat.jl index 2b1eae369..8543d9160 100644 --- a/src/Compat.jl +++ b/src/Compat.jl @@ -250,7 +250,7 @@ end if VERSION < v"1.8.0-DEV.487" export eachsplit - """ + @doc """ eachsplit(str::AbstractString, dlm; limit::Integer=0) eachsplit(str::AbstractString; limit::Integer=0) @@ -275,7 +275,8 @@ if VERSION < v"1.8.0-DEV.487" "Ma" "rch" ``` - """ + """ eachsplit + function eachsplit end struct SplitIterator{S<:AbstractString,F} @@ -421,7 +422,7 @@ end return nothing end - """ + @doc """ pkgversion(m::Module) Return the version of the package that imported module `m`, @@ -433,7 +434,8 @@ end To get the version of the package that imported the current module the form `pkgversion(@__MODULE__)` can be used. - """ + """ pkgversion + function pkgversion(m::Module) path = pkgdir(m) path === nothing && return nothing @@ -457,7 +459,7 @@ if VERSION < v"1.9.0-DEV.1163" import Base: IteratorSize, HasLength, HasShape, OneTo export stack - """ + @doc """ stack(iter; [dims]) Combine a collection of arrays (or other iterable objects) of equal size @@ -547,10 +549,11 @@ if VERSION < v"1.9.0-DEV.1163" julia> hvcat(5, M...) |> size # hvcat puts matrices next to each other (14, 15) ``` - """ + """ stack + stack(iter; dims=:) = _stack(dims, iter) - """ + @doc """ stack(f, args...; [dims]) Apply a function to each element of a collection, and `stack` the result. @@ -576,7 +579,8 @@ if VERSION < v"1.9.0-DEV.1163" 1.0 2.0 3.0 10.0 20.0 30.0 0.1 0.2 0.3 4.0 5.0 6.0 400.0 500.0 600.0 0.04 0.05 0.06 ``` - """ + """ stack(f, iter) + stack(f, iter; dims=:) = _stack(dims, f(x) for x in iter) stack(f, xs, yzs...; dims=:) = _stack(dims, f(xy...) for xy in zip(xs, yzs...)) @@ -803,6 +807,226 @@ if VERSION < v"1.11.0-DEV.1579" Iterators.cycle(xs, n::Integer) = Iterators.flatten(Iterators.repeated(xs, n)) end +# https://github.com/JuliaLang/julia/pull/39071 +if !isdefined(Base, :logrange) # VERSION < v"1.12.0-DEV.2" or appropriate 1.11.x after backporting + + export logrange + + @doc """ + logrange(start, stop, length) + logrange(start, stop; length) + + Construct a specialized array whose elements are spaced logarithmically + between the given endpoints. That is, the ratio of successive elements is + a constant, calculated from the length. + + This is similar to `geomspace` in Python. Unlike `PowerRange` in Mathematica, + you specify the number of elements not the ratio. + Unlike `logspace` in Python and Matlab, the `start` and `stop` arguments are + always the first and last elements of the result, not powers applied to some base. + + # Examples + ``` + julia> logrange(10, 4000, length=3) + 3-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}: + 10.0, 200.0, 4000.0 + + julia> ans[2] ≈ sqrt(10 * 4000) # middle element is the geometric mean + true + + julia> range(10, 40, length=3)[2] ≈ (10 + 40)/2 # arithmetic mean + true + + julia> logrange(1f0, 32f0, 11) + 11-element Base.LogRange{Float32, Float64}: + 1.0, 1.41421, 2.0, 2.82843, 4.0, 5.65685, 8.0, 11.3137, 16.0, 22.6274, 32.0 + + julia> logrange(1, 1000, length=4) ≈ 10 .^ (0:3) + true + ``` + + See the [`Compat.LogRange`](@ref Compat.LogRange) type for further details. + + !!! compat "Julia 1.9" + The version of this struct in Compat.jl does not use `Base.TwicePrecision{Float64}` + before Julia 1.9, so it sometimes has larger floating-point errors on intermediate points. + + !!! compat "Julia 1.11" + The printing of Compat.jl's version of the struct is also different, + less like `LinRange` and more like `Vector`. + """ logrange + + logrange(start::Real, stop::Real, length::Integer) = LogRange(start, stop, Int(length)) + logrange(start::Real, stop::Real; length::Integer) = logrange(start, stop, length) + + @doc """ + LogRange{T}(start, stop, len) <: AbstractVector{T} + + A range whose elements are spaced logarithmically between `start` and `stop`, + with spacing controlled by `len`. Returned by [`logrange`](@ref). + + Like [`LinRange`](@ref), the first and last elements will be exactly those + provided, but intermediate values may have small floating-point errors. + These are calculated using the logs of the endpoints, which are + stored on construction, often in higher precision than `T`. + + !!! compat "Julia 1.9" + The version of this struct in Compat.jl does not use `Base.TwicePrecision{Float64}` + before Julia 1.9. Therefore it has larger floating-point errors on intermediate + points than shown below. + + !!! compat "Julia 1.11" + The printing of Compat.jl's version of the struct is also different, + less like `LinRange` and more like `Vector`. + + # Examples + ``` + julia> logrange(1, 4, length=5) + 5-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}: + 1.0, 1.41421, 2.0, 2.82843, 4.0 + + julia> Base.LogRange{Float16}(1, 4, 5) + 5-element Base.LogRange{Float16, Float64}: + 1.0, 1.414, 2.0, 2.828, 4.0 + + julia> logrange(1e-310, 1e-300, 11)[1:2:end] + 6-element Vector{Float64}: + 1.0e-310 + 9.999999999999974e-309 + 9.999999999999981e-307 + 9.999999999999988e-305 + 9.999999999999994e-303 + 1.0e-300 + + julia> prevfloat(1e-308, 5) == ans[2] + true + ``` + + Note that integer eltype `T` is not allowed. + Use for instance `round.(Int, xs)`, or explicit powers of some integer base: + + ``` + julia> xs = logrange(1, 512, 4) + 4-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}: + 1.0, 8.0, 64.0, 512.0 + + julia> 2 .^ (0:3:9) |> println + [1, 8, 64, 512] + ``` + """ LogRange + + struct LogRange{T<:Real,X} <: AbstractArray{T,1} + start::T + stop::T + len::Int + extra::Tuple{X,X} + function LogRange{T}(start::T, stop::T, len::Int) where {T<:Real} + if T <: Integer + # LogRange{Int}(1, 512, 4) produces InexactError: Int64(7.999999999999998) + throw(ArgumentError("LogRange{T} does not support integer types")) + end + if iszero(start) || iszero(stop) + throw(DomainError((start, stop), + "LogRange cannot start or stop at zero")) + elseif start < 0 || stop < 0 + # log would throw, but _log_twice64_unchecked does not + throw(DomainError((start, stop), + "LogRange does not accept negative numbers")) + elseif !isfinite(start) || !isfinite(stop) + throw(DomainError((start, stop), + "LogRange is only defined for finite start & stop")) + elseif len < 0 + throw(ArgumentError(string( # LazyString( + "LogRange(", start, ", ", stop, ", ", len, "): can't have negative length"))) + elseif len == 1 && start != stop + throw(ArgumentError(string( # LazyString( + "LogRange(", start, ", ", stop, ", ", len, "): endpoints differ, while length is 1"))) + end + ex = _logrange_extra(start, stop, len) + new{T,typeof(ex[1])}(start, stop, len, ex) + end + end + + function LogRange{T}(start::Real, stop::Real, len::Integer) where {T} + LogRange{T}(convert(T, start), convert(T, stop), convert(Int, len)) + end + function LogRange(start::Real, stop::Real, len::Integer) + T = float(promote_type(typeof(start), typeof(stop))) + LogRange{T}(convert(T, start), convert(T, stop), convert(Int, len)) + end + + Base.size(r::LogRange) = (r.len,) + Base.length(r::LogRange) = r.len + + Base.first(r::LogRange) = r.start + Base.last(r::LogRange) = r.stop + + function _logrange_extra(a::Real, b::Real, len::Int) + loga = log(1.0 * a) # widen to at least Float64 + logb = log(1.0 * b) + (loga/(len-1), logb/(len-1)) + end + + function Base.getindex(r::LogRange{T}, i::Int) where {T} + @inline + @boundscheck checkbounds(r, i) + i == 1 && return r.start + i == r.len && return r.stop + # Main path uses Math.exp_impl for TwicePrecision, but is not perfectly + # accurate, hence the special cases for endpoints above. + logx = (r.len-i) * r.extra[1] + (i-1) * r.extra[2] + x = _exp_allowing_twice64(logx) + return copysign(T(x), r.start) + end + + function Base.show(io::IO, r::LogRange{T}) where {T} + print(io, "LogRange{", T, "}(") + ioc = IOContext(io, :typeinfo => T) + show(ioc, first(r)) + print(io, ", ") + show(ioc, last(r)) + print(io, ", ") + show(io, length(r)) + print(io, ')') + end + + # Display LogRange like LinRange -- PR widened signature of print_range to allow this + # function Base.show(io::IO, ::MIME"text/plain", r::LogRange) + # isempty(r) && return show(io, r) + # summary(io, r) + # println(io, ":") + # print_range(io, r, " ", ", ", "", " \u2026 ") + # end + + _exp_allowing_twice64(x::Number) = exp(x) + + if VERSION >= v"1.9.0-DEV.318" # Julia PR #44717 allows this high-precision path: + + _exp_allowing_twice64(x::Base.TwicePrecision{Float64}) = Base.Math.exp_impl(x.hi, x.lo, Val(:ℯ)) + + function _log_twice64_unchecked(x::Float64) + xu = reinterpret(UInt64, x) + if xu < (UInt64(1)<<52) # x is subnormal + xu = reinterpret(UInt64, x * 0x1p52) # normalize x + xu &= ~Base.sign_mask(Float64) + xu -= UInt64(52) << 52 # mess with the exponent + end + Base.TwicePrecision(Base.Math._log_ext(xu)...) + end + + function _logrange_extra(a::Float64, b::Float64, len::Int) + loga = _log_twice64_unchecked(a) + logb = _log_twice64_unchecked(b) + # The reason not to do linear interpolation on log(a)..log(b) in `getindex` is + # that division of TwicePrecision is quite slow, so do it once on construction: + (loga/(len-1), logb/(len-1)) + end + end +else + # Ensure that Compat.LogRange is always this struct, not exported from Base + using Base: LogRange +end + include("deprecated.jl") end # module Compat diff --git a/test/runtests.jl b/test/runtests.jl index 35631190e..f036601fb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -780,3 +780,68 @@ end Base.haslength(cycle(0:3, 2)) == false # but not sure we should test these Base.IteratorSize(cycle(0:3, 2)) == Base.SizeUnknown() end + +# https://github.com/JuliaLang/julia/pull/39071 +@testset "logrange" begin + # basic idea + @test logrange(2, 16, 4) ≈ [2, 4, 8, 16] + @test logrange(1/8, 8.0, 7) ≈ [0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0] + @test logrange(1000, 1, 4) ≈ [1000, 100, 10, 1] + @test logrange(1, 10^9, 19)[1:2:end] ≈ 10 .^ (0:9) + + # endpoints + @test logrange(0.1f0, 100, 33)[1] === 0.1f0 + @test logrange(0.789, 123_456, 135_790)[[begin, end]] == [0.789, 123_456] + @test logrange(nextfloat(0f0), floatmax(Float32), typemax(Int))[end] === floatmax(Float32) + @test logrange(nextfloat(Float16(0)), floatmax(Float16), 66_000)[end] === floatmax(Float16) + @test first(logrange(pi, 2pi, 3000)) === logrange(pi, 2pi, 3000)[1] === Float64(pi) + if Int == Int64 + @test logrange(0.1, 1000, 2^54)[end] === 1000.0 + end + + # empty, only, constant + @test first(logrange(1, 2, 0)) === 1.0 + @test last(logrange(1, 2, 0)) === 2.0 + @test collect(logrange(1, 2, 0)) == Float64[] + @test only(logrange(2pi, 2pi, 1)) === logrange(2pi, 2pi, 1)[1] === 2pi + @test logrange(1, 1, 3) == fill(1.0, 3) + + # subnormal Float64 + x = logrange(1e-320, 1e-300, 21) .* 1e300 + @test x ≈ logrange(1e-20, 1, 21) rtol=1e-6 + + # types + @test eltype(logrange(1, 10, 3)) == Float64 + @test eltype(logrange(1, 10, Int32(3))) == Float64 + @test eltype(logrange(1, 10f0, 3)) == Float32 + @test eltype(logrange(1f0, 10, 3)) == Float32 + @test eltype(logrange(1, big(10), 3)) == BigFloat + @test logrange(big"0.3", big(pi), 50)[1] == big"0.3" + @test logrange(big"0.3", big(pi), 50)[end] == big(pi) + + # more constructors + @test logrange(1,2,length=3) === Compat.LogRange(1,2,3) == Compat.LogRange{Float64}(1,2,3) + @test logrange(1f0, 2f0, length=3) == Compat.LogRange{Float32}(1,2,3) + + # errors + @test_throws UndefKeywordError logrange(1, 10) # no default length + @test_throws ArgumentError logrange(1, 10, -1) # negative length + @test_throws ArgumentError logrange(1, 10, 1) # endpoints must not differ + @test_throws DomainError logrange(1, -1, 3) # needs complex numbers + @test_throws DomainError logrange(-1, -2, 3) # not supported, for now + @test_throws MethodError logrange(1, 2+3im, length=4) # not supported, for now + @test_throws ArgumentError logrange(1, 10, 2)[true] # bad index + @test_throws BoundsError logrange(1, 10, 2)[3] + @test_throws ArgumentError Compat.LogRange{Int}(1,4,5) # no integer ranges + @test_throws MethodError Compat.LogRange(1,4, length=5) # type does not take keyword + # (not sure if these should ideally be DomainError or ArgumentError) + @test_throws DomainError logrange(1, Inf, 3) + @test_throws DomainError logrange(0, 2, 3) + @test_throws DomainError logrange(1, NaN, 3) + @test_throws DomainError logrange(NaN, 2, 3) + + # printing + @test repr(Compat.LogRange(1,2,3)) == "LogRange{Float64}(1.0, 2.0, 3)" # like 2-arg show + @test_skip repr("text/plain", Compat.LogRange(1,2,3)) == "3-element Compat.LogRange{Float64, Base.TwicePrecision{Float64}}:\n 1.0, 1.41421, 2.0" + @test_skip repr("text/plain", Compat.LogRange(1,2,0)) == "LogRange{Float64}(1.0, 2.0, 0)" # empty case +end