Construct generic results via lohi/midrad

Project.toml

@@ -1,7 +1,7 @@
 name = "ThickNumbers"
 uuid = "b57aa878-5b76-4266-befc-f8e007760995"
 authors = ["Tim Holy <tim.holy@gmail.com> and contributors"]
-version = "1.1.1"
+version = "1.1.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

ThickNumbersInterfaceTests/test/runtests.jl

@@ -8,6 +8,7 @@ if testdir ∉ LOAD_PATH
 end
 
 using IntervalArith
+using MidRadArith
 
 @testset "Interface tests" begin
     ThickNumbersInterfaceTests.test_reserved()
@@ -16,6 +17,13 @@ using IntervalArith
     ThickNumbersInterfaceTests.test_optional(Interval{Float64})
     ThickNumbersInterfaceTests.test_optional(Interval, [Float32, Float64])
     ThickNumbersInterfaceTests.test_FPTNviolations(Interval(1.0, 2.0))
+
+    # `MidRad` is parametrized by midpoint and radius rather than by lo/hi
+    ThickNumbersInterfaceTests.test_required(MidRad{Float64})
+    ThickNumbersInterfaceTests.test_required(MidRad, [Float32, Float64])
+    ThickNumbersInterfaceTests.test_optional(MidRad{Float64})
+    ThickNumbersInterfaceTests.test_optional(MidRad, [Float32, Float64])
+    ThickNumbersInterfaceTests.test_FPTNviolations(MidRad(1.5, 0.5))
 end
 
 filter!(LOAD_PATH) do path

ThickNumbersInterfaceTests/test/testpackages/MidRadArith/Project.toml

@@ -0,0 +1,7 @@
+name = "MidRadArith"
+uuid = "739db5bf-5d3a-4748-80ba-ed5827d60a8c"
+authors = ["Tim Holy <tim.holy@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+ThickNumbers = "b57aa878-5b76-4266-befc-f8e007760995"

ThickNumbersInterfaceTests/test/testpackages/MidRadArith/src/MidRadArith.jl

@@ -0,0 +1,67 @@
+module MidRadArith
+
+using ThickNumbers
+
+export MidRad
+
+# To avoid ambiguity with the ForwardDiff extension, it's easiest to be specific about the promotion of
+# other Numbers against `MidRad`
+const BaseReals = Union{AbstractFloat, Integer, AbstractIrrational, Rational}
+
+"""
+    MidRad(mid, rad)
+
+A `ThickNumber` spanning `[mid-rad, mid+rad]`, stored by its midpoint and radius.
+
+`MidRad` exists to exercise the generic code against a parametrization that is
+*not* `(loval, hival)`: its two-argument constructor takes a radius where
+`Interval`'s takes an upper bound. Generic code that builds a result by calling
+the two-argument constructor with a span will therefore produce a wrong answer
+for `MidRad` while still looking correct for `Interval`.
+"""
+struct MidRad{T<:Number} <: ThickNumber{T}
+    mid::T
+    rad::T
+end
+MidRad(mid, rad) = MidRad(promote(mid, rad)...)
+MidRad{T}(x::MidRad) where T = MidRad{T}(x.mid, x.rad)
+MidRad{T}(x::Number) where T = MidRad{T}(x, zero(x))
+
+ThickNumbers.loval(x::MidRad) = x.mid - x.rad
+ThickNumbers.hival(x::MidRad) = x.mid + x.rad
+
+# Read the stored parametrization directly; recovering these from the span would
+# round-trip through `mid ± rad` and can shrink the radius by an ulp, which
+# `midrad` forbids (a radius may only err outward).
+ThickNumbers.mid(x::MidRad) = x.mid
+ThickNumbers.rad(x::MidRad) = x.rad
+ThickNumbers.lohi(::Type{M}, lo, hi) where M<:MidRad = M((lo + hi)/2, (hi - lo)/2)
+ThickNumbers.basetype(::Type{MidRad{T}}) where T = MidRad
+ThickNumbers.basetype(::Type{MidRad}) = MidRad
+
+# The natural parametrization: define it directly rather than round-tripping through `lohi`
+ThickNumbers.midrad(::Type{MidRad{T}}, mid, rad) where T = MidRad{T}(mid, rad)
+ThickNumbers.midrad(::Type{MidRad}, mid::T, rad::T) where T = MidRad{T}(mid, rad)
+ThickNumbers.midrad(::Type{MidRad}, mid, rad) = midrad(MidRad, promote(mid, rad)...)
+
+# The default empty set spans [typemax, typemin], whose midpoint is `NaN`; a negative
+# radius is the representable way to make `hival < loval` in this parametrization.
+ThickNumbers.emptyset(::Type{MidRad{T}}) where T = MidRad{T}(zero(T), -oneunit(T))
+
+# Promotion of `valuetype`
+Base.promote_rule(::Type{MidRad{S}}, ::Type{MidRad{T}}) where {S<:Number, T<:Number} = MidRad{promote_type(T, S)}
+Base.promote_rule(::Type{MidRad{S}}, ::Type{T}) where {S<:Number, T<:BaseReals} = MidRad{promote_type(T, S)}
+
+# Very basic arithmetic needed for `norm` (this would be fleshed out in real applications)
+Base.:+(x::MidRad, y::MidRad) = MidRad(x.mid + y.mid, x.rad + y.rad)
+Base.:/(x::MidRad, y::Real) = MidRad(x.mid / y, x.rad / abs(y))
+function Base.:*(x::MidRad, y::MidRad)
+    T = typeof(zero(valuetype(x))*zero(valuetype(y)))
+    (isempty(x) || isempty(y)) && return emptyset(MidRad{T})
+    v1, v2, v3, v4 = loval(x)*loval(y), hival(x)*loval(y), loval(x)*hival(y), hival(x)*hival(y)
+    return lohi(MidRad, min(v1, v2, v3, v4), max(v1, v2, v3, v4))
+end
+
+Base.conj(x::MidRad{T}) where T = MidRad{T}(conj(x.mid), conj(x.rad))
+
+end # module MidRadArith

docs/src/developers.md

@@ -36,6 +36,12 @@ There are also numerous optional methods you can specialize if it makes `MyType`
 efficiently. For example, a Gaussian random variable package might want to implement [`midrad(MyType{T}, center, σ)`](@ref) to construct values directly, assuming this is the natural parametrization
 of this type.
 
+`MyType` need not be parametrized by its `lo` and `hi` values: generic code constructs its
+results only through [`lohi`](@ref) and [`midrad`](@ref), never by calling a two-argument
+`MyType(a, b)` constructor. One caveat: the default [`emptyset`](@ref) spans
+`[typemax(T), typemin(T)]`, whose midpoint is `NaN`, so a type parametrized by a midpoint
+must specialize `emptyset`; a negative radius is the natural choice.
+
 ## Ensuring compliance with the ThickNumbers interface
 
 The `ThickNumbersInterfaceTests` package can be used to determine whether your implementations comply with the requirements.  As it is possible that this test suite will evolve and add new requirements,

src/ThickNumbers.jl

@@ -334,8 +334,10 @@ isnan_tn(a::ThickNumber) = isnan(loval(a)) | isnan(hival(a))
 Base.isnan(::ThickNumber) = throw(FPTNException(isnan, isnan_tn))
 
 
-Base.typemin(::Type{TN}) where TN<:ThickNumber{T} where T<:Number = TN(typemin(T), typemin(T))
-Base.typemax(::Type{TN}) where TN<:ThickNumber{T} where T<:Number = TN(typemax(T), typemax(T))
+# Construct via `midrad` rather than `lohi`: for a type parametrized by midpoint
+# and radius, a span of [typemin, typemin] would compute a radius of `Inf - Inf`.
+Base.typemin(::Type{TN}) where TN<:ThickNumber{T} where T<:Number = midrad(TN, typemin(T), zero(T))
+Base.typemax(::Type{TN}) where TN<:ThickNumber{T} where T<:Number = midrad(TN, typemax(T), zero(T))
 Base.typemin(x::ThickNumber) = typemin(typeof(x))
 Base.typemax(x::ThickNumber) = typemax(typeof(x))
 
@@ -415,7 +417,7 @@ end
 
 function Base.intersect(a::TN, b::TN) where TN<:ThickNumber
     isdisjoint(a,b) && return emptyset(TN)
-    TN(max(loval(a), loval(b)), min(hival(a), hival(b)))
+    lohi(TN, max(loval(a), loval(b)), min(hival(a), hival(b)))
 end
 Base.intersect(a::ThickNumber, b::ThickNumber) = intersect(promote(a, b)...)
 Base.intersect(a::ThickNumber, b::ThickNumber, c::ThickNumber...) = intersect(intersect(a, b), c...)
@@ -428,7 +430,7 @@ Construct a `ThickNumber` containing `a`, `b`, and `c...`.
 hull(a::ThickNumber, b::ThickNumber, c::ThickNumber...) = hull(hull(a, b), c...)
 hull(a::ThickNumber, b::ThickNumber) = hull(promote(a, b)...)
 hull(a::TN, b::TN) where TN<:ThickNumber =
-    TN(min(loval(a), loval(b)), max(hival(a), hival(b)))
+    lohi(TN, min(loval(a), loval(b)), max(hival(a), hival(b)))
 
 ## Operators
 

test/runtests.jl

@@ -4,6 +4,7 @@ using Test
 include("setpath.jl")
 
 using IntervalArith
+using MidRadArith
 
 @testset "ThickNumbers generics" begin
     # Test all the operations defined in ThickNumbers
@@ -143,6 +144,34 @@ end
     @test clamp(Interval(-1, 4), 2, 3) === Interval(2, 3)
 end
 
+# Generic code must construct results with `lohi`/`midrad` rather than with a
+# two-argument `TN(lo, hi)` call, which assumes a lo/hi parametrization. `MidRad`
+# stores a midpoint and radius, so such a call silently reinterprets the upper
+# bound as a radius; the resulting span is wrong while `Interval` still looks fine.
+@testset "parametrization independence: $TN" for TN in (Interval, MidRad)
+    span(x) = (loval(x), hival(x))
+
+    @test span(typemin(TN{Float64})) === (-Inf, -Inf)
+    @test span(typemax(TN{Float64})) === (Inf, Inf)
+
+    a = lohi(TN, 1.0, 3.0)
+    b = lohi(TN, 2.0, 5.0)   # overlaps `a`
+    c = lohi(TN, 7.0, 8.0)   # disjoint from `a`
+
+    @test span(hull(a, b)) === (1.0, 5.0)
+    @test span(hull(a, c)) === (1.0, 8.0)
+    @test span(hull(a, a)) === (1.0, 3.0)
+    @test span(hull(a, b, c)) === (1.0, 8.0)
+
+    @test span(intersect(a, b)) === (2.0, 3.0)
+    @test span(intersect(a, a)) === (1.0, 3.0)
+    @test isempty(intersect(a, c))
+
+    # A hull contains both operands; an intersection is contained by both.
+    @test a ⫃ hull(a, b) && b ⫃ hull(a, b)
+    @test intersect(a, b) ⫃ a && intersect(a, b) ⫃ b
+end
+
 include(joinpath("extensions", "runtests.jl"))
 
 cleanup()