JuliaIntervals / IntervalArithmetic.jl

Library for validated numerics using interval arithmetic
https://juliaintervals.github.io/IntervalArithmetic.jl/
Other
299 stars 71 forks source link

abstract type #2

Closed JeffreySarnoff closed 11 months ago

JeffreySarnoff commented 7 years ago

We should follow Tim Holy and use

abstract type AbstractInterval{T} end

parameterizing the underlying type that fits well with an emerging pattern AbstractTime{T} etc.

JeffreySarnoff commented 7 years ago

I would go so far as to have a separate micro-API-ish package

module AbstractIntervals

export AbstractInterval, LoHiInterval 

abstract type AbstractInterval{T} end

struct LoHiInterval{T} <: AbstractInterval{T}
    lo::T
    hi::T
end

end # module

as that will help get other interval packages all conformant

dpsanders commented 7 years ago

In principle I certainly agree.

However, currently we need to have intervals be a subtype of Real for use with ForwardDiff.jl.

JeffreySarnoff commented 7 years ago

are you referring to dual.jl or another file?

dpsanders commented 7 years ago

Yes, to be able to create Duals of Intervals, to do e.g.

julia v0.5> using ValidatedNumerics, ForwardDiff

julia v0.5> f(x) = x^2 - 2;

julia v0.5> X = 3..4   # Interval
[3, 4]

julia v0.5> ForwardDiff.derivative(f, X)
[6, 8]

to get an enclosure of the derivative of f over the interval X.

dpsanders commented 7 years ago

And that restriction on Dual (that it takes Reals) will not be able to be removed until we have traits in Julia.

Of course, we could certainly do

abstract type AbstractInterval{T} <: Real end

but I'm not sure if Tim will be happy with that for IntervalSets.jl.

JeffreySarnoff commented 7 years ago
abstract type AbstractInterval{T, Orderedness<TotallyOrdered, e.g. linearly connected>} <: Real end

maybe you have a better fitting second param (to become trait, probably)

JeffreySarnoff commented 7 years ago

@timholy What are your thoughts? Is this a good place for your way of entraiting (only if it would easily slip into the Julia future realization ... this would propagate). Is a second (ignorable) param better/worse/ much the same?

(why can we not already pre-furgate mathematical abstractions without tangling the numerical hierarchy?)

timholy commented 7 years ago

In principle, I don't think an interval should be a subtype of Real or even Number; anything that can be thought of as a pair of numbers certainly isn't the same thing as a single number, even if you can define a lot of similar operations for them. I would say this is even clearer than the case of whether a Factorization should be an AbstractArray:

julia> Factorization <: AbstractArray
false

even though a Factorization is simply a representation of an AbstractArray F = A*B.

More practically, I think (hope) I understand the bind this puts you in. Could one relax the requirement on Real in ForwardDiff? I did some brief experiments and it seems one can get fairly far:

$ git diff
diff --git a/src/dual.jl b/src/dual.jl
index dfed492..0b56063 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -4,7 +4,7 @@ const ExternalReal = Union{subtypes(Real)...}
 # Dual #
 ########

-immutable Dual{N,T<:Real} <: Real
+immutable Dual{N,T}
     value::T
     partials::Partials{N,T}
 end
@@ -13,7 +13,7 @@ end
 # Constructors #
 ################

-Dual{N,T}(value::T, partials::Partials{N,T}) = Dual{N,T}(value, partials)
+# Dual{N,T}(value::T, partials::Partials{N,T}) = Dual{N,T}(value, partials)

 function Dual{N,A,B}(value::A, partials::Partials{N,B})
     T = promote_type(A, B)
@@ -144,11 +144,11 @@ end

 isconstant(n::Dual) = iszero(partials(n))

-@ambiguous Base.isequal{N}(a::Dual{N}, b::Dual{N}) = isequal(value(a), value(b))
-@ambiguous Base.:(==){N}(a::Dual{N}, b::Dual{N}) = value(a) == value(b)
-@ambiguous Base.isless{N}(a::Dual{N}, b::Dual{N}) = value(a) < value(b)
-@ambiguous Base.:<{N}(a::Dual{N}, b::Dual{N}) = isless(a, b)
-@ambiguous Base.:(<=){N}(a::Dual{N}, b::Dual{N}) = <=(value(a), value(b))
+Base.isequal{N}(a::Dual{N}, b::Dual{N}) = isequal(value(a), value(b))
+Base.:(==){N}(a::Dual{N}, b::Dual{N}) = value(a) == value(b)
+Base.isless{N}(a::Dual{N}, b::Dual{N}) = value(a) < value(b)
+Base.:<{N}(a::Dual{N}, b::Dual{N}) = isless(a, b)
+Base.:(<=){N}(a::Dual{N}, b::Dual{N}) = <=(value(a), value(b))

 for T in (AbstractFloat, Irrational, Real)
     Base.isequal(n::Dual, x::T) = isequal(value(n), x)
@@ -210,11 +210,11 @@ Base.float{N,T}(n::Dual{N,T}) = Dual{N,promote_type(T, Float16)}(n)
 # Addition/Subtraction #
 #----------------------#

-@ambiguous @inline Base.:+{N}(n1::Dual{N}, n2::Dual{N}) = Dual(value(n1) + value(n2), partials(n1) + partials(n2))
+@inline Base.:+{N}(n1::Dual{N}, n2::Dual{N}) = Dual(value(n1) + value(n2), partials(n1) + partials(n2))
 @ambiguous @inline Base.:+(n::Dual, x::Real) = Dual(value(n) + x, partials(n))
 @ambiguous @inline Base.:+(x::Real, n::Dual) = n + x

-@ambiguous @inline Base.:-{N}(n1::Dual{N}, n2::Dual{N}) = Dual(value(n1) - value(n2), partials(n1) - partials(n2))
+@inline Base.:-{N}(n1::Dual{N}, n2::Dual{N}) = Dual(value(n1) - value(n2), partials(n1) - partials(n2))
 @ambiguous @inline Base.:-(n::Dual, x::Real) = Dual(value(n) - x, partials(n))
 @ambiguous @inline Base.:-(x::Real, n::Dual) = Dual(x - value(n), -(partials(n)))
 @inline Base.:-(n::Dual) = Dual(-(value(n)), -(partials(n)))
@@ -225,7 +225,7 @@ Base.float{N,T}(n::Dual{N,T}) = Dual{N,promote_type(T, Float16)}(n)
 @inline Base.:*(n::Dual, x::Bool) = x ? n : (signbit(value(n))==0 ? zero(n) : -zero(n))
 @inline Base.:*(x::Bool, n::Dual) = n * x

-@ambiguous @inline function Base.:*{N}(n1::Dual{N}, n2::Dual{N})
+@inline function Base.:*{N}(n1::Dual{N}, n2::Dual{N})
     v1, v2 = value(n1), value(n2)
     return Dual(v1 * v2, _mul_partials(partials(n1), partials(n2), v2, v1))
 end
@@ -236,7 +236,7 @@ end
 # Division #
 #----------#

-@ambiguous @inline function Base.:/{N}(n1::Dual{N}, n2::Dual{N})
+@inline function Base.:/{N}(n1::Dual{N}, n2::Dual{N})
     v1, v2 = value(n1), value(n2)
     return Dual(v1 / v2, _div_partials(partials(n1), partials(n2), v1, v2))
 end
@@ -254,7 +254,7 @@ end

 for f in (:(Base.:^), :(NaNMath.pow))
     @eval begin
-        @ambiguous @inline function ($f){N}(n1::Dual{N}, n2::Dual{N})
+        @inline function ($f){N}(n1::Dual{N}, n2::Dual{N})
             v1, v2 = value(n1), value(n2)
             expv = ($f)(v1, v2)
             powval = v2 * ($f)(v1, v2 - 1)

It doesn't pass DualTests.jl but it's clear that most stuff works. I'm not going to push forward on this myself, but perhaps others could pick it up.

timholy commented 7 years ago

Though one should get in touch with @jrevels before going down this road too far.

If just removing the type restriction won't fly, one can still control dispatch even without "official" traits. One can use SimpleTraits.jl or just write them out manually:

foo(x) = _foo(DualWorthy(x), x)
_foo(::SafeForDual, x) = blah blah
_foo(::Any, x) = error(x, " is not worthy of being dual-ified")

DualWorthy can satisfy an interface somewhat like another trait that Base uses heavily, IndexStyle. Or you could even make DualWorthy(x) itself throw the error.

dpsanders commented 7 years ago

In Julia at the moment, there simply isn't any way of saying that an interval (or a dual number, for that matter) is "something number-like" (e.g. is something like a ring, in the abstract algebra sense), without making it a subtype of Real or Number.

In that sense, currently, Real in Base seems to be something approximating "real numbers", but in the wider ecosystem has come to take on rather the meaning of "something number-like".

I was discussing this offline with @jrevels, who may wish to comment about ForwardDiff.jl.

timholy commented 7 years ago

We may have had a "race condition" in our comments, but in case not: DualWorthy would allow you to opt-in for any type that you are willing to support the necessary operations for Dual. Types that don't support the needed interface could either automatically throw an error or could simply dispatch differently (that second _foo above could simply do something different).

timholy commented 7 years ago

The docs on the AbstractArray interface may be a helpful model showing how IndexStyle solves these problems.

dpsanders commented 7 years ago

Thanks @timholy, that's a very nice idea.

JeffreySarnoff commented 7 years ago

I want to look both ways -- making it work best with the ways open through v0.6 (as this has been focused) and getting as clear and pithy as possible on what would allow us the requisite expressiveness within the sense of Julia as developing to write it right. The above @timholy is nice. To get nice and clean requires deeper support and, it feels, renewed attention now-ish. Any ideas welcome. I'd like to hear from Bill Hart on this too (email sent), he has spent a good deal of energy around computing the mathy abstract.

dpsanders commented 7 years ago

@jrevels says that the restriction on the type of thing allowed inside a Dual could be relaxed.

Nonetheless, he pointed out that an Interval in the IntervalArithmetic.jl (formerly ValidatedNumerics.jl) package will still need to be <: Real, in order to be able to be injected into a lot of functions "out there" that require Reals to be passed in.

EDIT: This is the same reason that requires that Dual <: Real, to enable Duals to also be passed around through other pre-existing code. Both Dual and Interval behave sufficiently much like a real (extending reals, in some sense) that it makes sense to do this.

Since we do not have multiple dispatch in Julia, it is thus not possible to inherit from an AbstractInterval type that is not itself <: Real and impose that Interval itself is <: Real.

JeffreySarnoff commented 7 years ago

Are packages now using Real as a scalar-indicable neighborhood amenable to computation.?

dpsanders commented 7 years ago

I guess there are many packages that restrict arguments of functions to Real who have never thought about using Duals or Intervals, but in which it makes perfect sense to do so.

timholy commented 7 years ago

I often use Real when I'm assuming the existence of a total order, which an interval doesn't support.

dpsanders commented 7 years ago

Can you give an example where you use the total order?

On 8 Apr 2017 2:03 p.m., "Tim Holy" notifications@github.com wrote:

I often use Real when I'm assuming the existence of a total order, which an interval doesn't support.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/2#issuecomment-292738738, or mute the thread https://github.com/notifications/unsubscribe-auth/AALtTojYoZrpUnHsYsgfJLAkpo1K9VeRks5rt9oPgaJpZM4M12DS .

timholy commented 7 years ago

If I want to sort something? In other words, if I'm restricting the types allowed to be in a specific container, we don't currently have a HasTotalOrder(T) trait, so I often use Real as a surrogate.

dpsanders commented 7 years ago

Good point. But I expect (though I don't currently have evidence) that there are numerical functions whose arguments are specified as Real in order to exclude eg Complex, and that would work with intervals.

On 8 Apr 2017 2:15 p.m., "Tim Holy" notifications@github.com wrote:

If I want to sort something? In other words, if I'm restricting the types allowed to be in a specific container, we don't currently have a HasTotalOrder(T) trait, so I often use Real as a surrogate.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/2#issuecomment-292739419, or mute the thread https://github.com/notifications/unsubscribe-auth/AALtTmVrW-wOIEP0o9pmv5e4Ksx68X_Dks5rt9y7gaJpZM4M12DS .

JeffreySarnoff commented 7 years ago

There are sorts of intervals over which a total ordering is available. Subsegments of a line segment are an example. You are constraining the nature of interval more than is necessary. How about preceeding_finitely_represented_value? My initial interest in coding Julia with interval types was to explore geometries and intervals with interval-valued bounds.

JeffreySarnoff commented 7 years ago

@timholy "DualWorthy would allow you to opt-in for any type that you are willing to support the necessary operations for Dual" -- I would expect the author of the type Dual to have defined enough so that what is undefined is properly handled by fallback implementations. Or, somewhat less stringently -- to have covered arithmetic and elementary functions and for duals, quaternions the first and second derivatives. Otherwise it is not that generally useful.

JeffreySarnoff commented 7 years ago

Is there something that comes with @JeffBezanson type-system-bettering and its triangular dispatch that lets a type parameter itself be parametrized (directly)?

timholy commented 7 years ago

There are sorts of intervals over which a total ordering is available.

But the answer to 3..5 < 4..5 is going to be misleading. The more reasonably-behaved orders are partial orders.

I would expect the author of the type Dual to have defined enough so that what is undefined is properly handled by fallback implementations

The whole point of traits is to define multiple behaviors. Look more carefully at how IndexStyle works: the user defines either getindex(a, i) or getindex(a, i, j...), not both (there's no harm in defining both, but typically there's no reason to do so). The whole point of the trait is to inform dispatch about which one has been defined; if you naively tried the wrong one, it would throw an error. The fallback methods are defined at a higher level, i.e., calling ind2sub or sub2ind to supply the "missing" variant automatically.

JeffreySarnoff commented 7 years ago

3/4 out of two is not bad :)

(a) should have been clearer, I use 3..5 ≺ 4..5 (precedes with durations, contains but is not contained by) (b) always good to learn

OlivierHnt commented 11 months ago

I am closing this because with PR #593 we decided to preserve Interval <: Real while still having a BareInterval that is not a subtype of Real.

Also, #585 feels like a successor to this issue.