Flexiunits

Project.toml

@@ -1,6 +1,6 @@
 name = "Roots"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "3.0.10"
+version = "3.0.11"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -30,6 +30,7 @@ BenchmarkTools = "1.4"
 ChainRulesCore = "1"
 ChainRulesTestUtils = "1"
 CommonSolve = "0.1, 0.2"
+FlexUnits = "0.6"
 ForwardDiff = "0.10, 1"
 IntervalRootFinding = "0.5, 0.6, 0.7"
 JET = "0.9, 0.10, 0.11, 0.12"
@@ -50,6 +51,7 @@ julia = "1.10"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+FlexUnits = "76e01b6b-c995-4ce6-8559-91e72a3d4e95"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
@@ -65,4 +67,4 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "ChainRulesTestUtils", "JET", "JSON", "SpecialFunctions", "Statistics", "Test", "BenchmarkTools", "ForwardDiff", "Measurements", "Polynomials", "Unitful", "Zygote"]
+test = ["Aqua", "ChainRulesTestUtils", "JET", "JSON", "SpecialFunctions", "Statistics", "Test", "BenchmarkTools", "FlexUnits", "ForwardDiff", "Measurements", "Polynomials", "Unitful", "Zygote"]

src/Bracketing/bisection.jl

@@ -95,8 +95,8 @@ function default_tolerances(
     ::AbstractUnivariateZeroState{T′,S′},
 ) where {T′,S′}
     T, S = real(float(T′)), real(float(S′))
-    xatol = eps(T)^3 * oneunit(T)
-    xrtol = eps(T) * one(T) # unitless
+    xatol = eps(one(T))^3 * oneunit(T)
+    xrtol = eps(one(T)) * one(T) # unitless
     atol = 0 * oneunit(S)
     rtol = 0 * one(S)
     maxiters = typemax(Int)

src/Bracketing/bracketing.jl

@@ -45,8 +45,8 @@ function default_tolerances(
     ::AbstractBracketingMethod,
     ::AbstractUnivariateZeroState{T,S},
 ) where {T,S}
-    xatol = eps(real(T))^3 * oneunit(real(T))
-    xrtol = eps(real(T))  # unitless
+    xatol = eps(one(real(T)))^3 * oneunit(real(T))
+    xrtol = eps(one(real(T)))  # unitless
     atol = zero(oneunit(real(S)))
     rtol = zero(one(real(S)))
     maxevals = 60
@@ -59,8 +59,8 @@ function default_tolerances(
     ::AbstractNonStrictBracketingMethod,
     ::AbstractUnivariateZeroState{T,S},
 ) where {T,S}
-    xatol = eps(real(T))^3 * oneunit(real(T))
-    xrtol = eps(real(T))  # unitless
+    xatol = eps(one(real(T)))^3 * oneunit(real(T))
+    xrtol = eps(one(real(T)))  # unitless
     atol = zero(oneunit(real(S)))
     rtol = zero(one(real(S)))
     maxevals = 60

src/Bracketing/modAB.jl

@@ -103,8 +103,8 @@ initial_fncalls(M::ModAB) = 2
 # where X is 16,32,or 64, N is -Int(log2(eps(T))) ÷ 2 + 1
 function default_tolerances(::ModAB, ::AbstractUnivariateZeroState{T,S}) where {T,S}
     xatol = 2 * eps(zero(T)) * oneunit(real(T)) # not quite 0
-    xrtol = eps(real(T))  # unitless
-    atol = 4 * eps(real(float(S))) * oneunit(real(S))
+    xrtol = eps(one(real(T)))  # unitless
+    atol = 4 * eps(one(real(float(S)))) * oneunit(real(S))
     rtol = zero(real(S))
     maxiters = 3 * (-Int(log2(eps(T))) + 1)
     strict = true

src/convergence.jl

@@ -53,10 +53,10 @@ function default_tolerances(
     ::AbstractUnivariateZeroMethod,
     ::AbstractUnivariateZeroState{T,S},
 ) where {T,S}
-    xatol = eps(real(T)) * oneunit(real(T))
-    xrtol = eps(real(T))  # unitless
-    atol = 4 * eps(real(float(S))) * oneunit(real(S))
-    rtol = 4 * eps(real(float(S))) * one(real(S))
+    xatol = eps(one(real(T))) * oneunit(real(T))
+    xrtol = eps(one(real(T))) # unitless
+    atol = 4 * eps(one(real(float(S)))) * oneunit(real(S))
+    rtol = 4 * eps(one(real(float(S)))) * one(real(S)) # unitless
     maxiters = 40
     strict = false
     (xatol, xrtol, atol, rtol, maxiters, strict)
@@ -96,7 +96,7 @@ function is_approx_zero_f(
 }
     ab, afb = abs(state.xn1), abs(state.fxn1)
     ϵₐ, ϵᵣ = options.abstol, options.reltol
-    Δ = max(_unitless(ϵₐ), _unitless(ab) * ϵᵣ)
+    Δ = max(_unitless(ϵₐ), _unitless(ab) * _unitless(ϵᵣ))
     afb ≤ Δ * oneunit(afb)
 end
 
@@ -128,6 +128,12 @@ function is_approx_zero_f(
 end
 
 ## --------------------------------------------------
+function __isapprox(a, b; atol=eps(), rtol=eps())
+    # bypass isapprox Real type constraint
+    a′, b′ = promote(a, b)
+    a == b && return true
+    return (isfinite(a) && isfinite(b) && abs(a-b) <= max(atol, rtol*max(abs(a′), abs(b′))))
+end
 
 # testing xₙ₊₁ - xₙ ≈ 0
 function iszero_Δx(
@@ -137,7 +143,7 @@ function iszero_Δx(
 )
     a, b, fa, fb = state.xn0, state.xn1, state.fxn0, state.fxn1
     δₐ, δᵣ = options.xabstol, options.xreltol
-    isapprox(a, b, atol=δₐ, rtol=δᵣ)
+    __isapprox(a, b, atol=δₐ, rtol=δᵣ)
 end
 
 # this is for Bisection where we can go adjacent floating point values

src/find_zeros.jl

@@ -68,15 +68,15 @@ Base.show(io::IO, alpha::Interval) = print(io, "($(alpha.a), $(alpha.b))")
 
 # check if f(a) is non zero using tolerances max(atol, eps()), rtol
 function _non_zero(fa, a::T, atol, rtol) where {T}
-    a, r = atol, abs(a) * rtol * oneunit(fa) / oneunit(a), oneunit(fa) * eps(T)
+    a, r = atol, abs(a) * rtol * oneunit(fa) / oneunit(a), oneunit(fa) * eps(one(T))
     return abs(fa) >= max(promote(a, r)...)
 end
 
 # After splitting by zeros we have intervals (zm, zn) this is used to shrink
 # to (zm+, zn-) where both are non-zeros, as defined above
 function find_non_zero(f, a::T, barrier, xatol, xrtol, atol, rtol) where {T}
     nan = (0 * a) / (0 * a) # try to get typed NaN
-    xtol = max(xatol, abs(a) * xrtol, oneunit(a) * eps(T))
+    xtol = max(xatol, abs(a) * xrtol, oneunit(a) * eps(one(T)))
     sgn = barrier > a ? 1 : -1
     ctr = 0
     x = a + 2^ctr * sgn * xtol
@@ -323,8 +323,8 @@ function find_zeros(f, a, b=nothing; no_pts=12, k=8, naive=false, kwargs...)
     T, S     = real(eltype(a0)), real(eltype(fa0))
     xatol::T = get(d, :xatol, eps(one(T))^(4 / 5) * oneunit(T))
     xrtol    = get(d, :xrtol, eps(one(T)) * one(T))
-    atol::S  = get(d, :atol, eps(float(S)) * oneunit(S))
-    rtol     = get(d, :rtol, eps(float(S)) * one(S))
+    atol::S  = get(d, :atol, eps(one(float(S))) * oneunit(S))
+    rtol     = get(d, :rtol, eps(one(float(S))) * one(S))
 
     zs = T[]  # collect zeros
 

test/test_composable.jl

@@ -5,6 +5,8 @@
 using Roots
 using Test
 using Unitful
+import FlexUnits
+import FlexUnits.UnitRegistry
 using Polynomials
 using ForwardDiff
 using Measurements
@@ -53,6 +55,39 @@ using Measurements
         @test Roots.is_small_Δx(prob.M, prob.state, prob.options) isa Bool  # does not throw
     end
 
+    # FlexUnits
+    @testset "find zero(s) with FlexUnits" begin
+        s = FlexUnits.UnitRegistry.u"s"
+        m = FlexUnits.UnitRegistry.u"m"
+        g = (9 + 8 // 10) * m / s^2
+        v0 = 10m / s
+        y0 = 16m
+        y(t) = -g * t^2 + v0 * t + y0
+
+        for order in orders
+            @test find_zero(y, 1.8s, order) ≈ 1.886053370668014s
+            @test find_zero(y, 1.8f0s, order) isa typeof(1.88f0s)
+        end
+
+        for M in [Roots.Bisection(), Roots.A42(), Roots.AlefeldPotraShi()]
+            @test find_zero(y, (1.8s, 1.9s), M) ≈ 1.886053370668014s
+            @test find_zero(y, (1.8f0s, 1.9f0s), M) isa typeof(1.88f0s)
+        end
+
+        xrts = find_zeros(y, 0s, 10s)
+        @test length(xrts) == 1
+        @test xrts[1] ≈ 1.886053370668014s
+
+        # issue #434
+        xzs1 = find_zeros(x -> cos(x / 1u"m"), -1.6u"m", 2u"m")
+        @test length(xzs1) == 2 && maximum(xzs1) ≈ 1.5707963267948966 * u"m"
+
+        FX = ZeroProblem(y, (0.0f0s, 2.0f0s))
+        prob = Roots.init(FX, Roots.AlefeldPotraShi())
+        @test Roots.is_small_Δx(prob.M, prob.state, prob.options) isa Bool  # does not throw
+    end
+
+
     # Polynomials
     @testset "find zero(s) with Polynomials" begin
         m = s = 1.0