Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

treat intervals more like the numbers they represent #181

Closed
wants to merge 4 commits into from

Conversation

gwater
Copy link
Contributor

@gwater gwater commented Jul 17, 2018

In order to get the ball rolling on #168 and #165 I'm starting this pull request. I imagine there is still a lot of discussion ahead of us but hopefully this will help us get a better idea and something concrete to interrogate.

The most significant change here is redefining == to be in line with its definition in Base:

Generic equality operator, giving a single Bool result. Falls back to ===. Should be implemented for all types with a notion of equality, based on the abstract value that an instance represents. […]

(emphasis added by me)

Obviously that breaks a lot of test cases and probably some parts of IntervalArithmetic and IntervalRootFinding which rely on == comparing the boundaries of two intervals. I have renamed the old behavior isequal(), in line with its definition in Base:

isequal is the comparison function used by hash tables (Dict). isequal(x,y) must imply that hash(x) == hash(y).

@dpsanders
Copy link
Member

Thanks for the impetus!

The behaviour of many of these functions is specified in the IEEE 1788-2015 standard document for interval arithmetic packages, so changing things around is not obvious.

@dpsanders
Copy link
Member

What is the change that you are proposing for ==?

@gwater
Copy link
Contributor Author

gwater commented Jul 18, 2018

The change I am proposing is to implement == so that it only returns true if two intervals contain exactly one number which is equal and false only if two intervals do not share any numbers ("a notion of equality, based on the abstract value that an instance represents").

I think there is a subtle misunderstanding of the relationship between Base and the IEEE standard at play here: Base specifies what an operator called == is supposed to do in julia. The IEEE standard specifies that any implementation must have an operation like current behavior of ==. Crucially, it does not prescribe what that operation is called. This is why I propose renaming the IEEE comparison isequal(), and replacing == with an operation that fits Base specification.

@dpsanders
Copy link
Member

I agree about the IEEE standard not specifying the names of functions.

But IMO, an interval represents a set, not a number, so it is quite reasonable for == and other relational operators to deal with sets, not numbers. The point of view of set computations is found e.g. in the book by Jaulin et al.

So we could certainly add a function equalasnumbers (or whatever name) that does what you suggest.

@gwater
Copy link
Contributor Author

gwater commented Jul 18, 2018

Ok, considering Intervals as representations of sets is of course valid and leads to a consistent implementation. However, the current implementation claims that Intervals are numbers:

abstract type AbstractInterval{T} <: Real end

Under the status quo that line leads to lots of problems because it implies that Intervals should be used as drop-in replacements for numbers. However, algorithms expecting numbers are silently failing because operations like iszero() and == treat intervals like sets. (We saw this in IntervalRootFinding, when StaticArrays changed their matrix inversion algorithm.) I hope you understand why I consider the status quo to be deeply flawed.

@dpsanders
Copy link
Member

However, the current implementation claims that Intervals are numbers:
abstract type AbstractInterval{T} <: Real end
Under the status quo that line leads to lots of problems

You have indeed hit a tricky point. I interpret this subtyping as saying something like "intervals behave like Reals for many purposes". You are right that there are probably many functions that don't work correctly and fail silently. (But you could also argue that many of them may not be "generic enough" in that case.)

On the other hand, if we drop this, then we lose a lot of important functionality, for example automatic differentiation with ForwardDiff.jl.

There is quite some previous discussion in #2. I am certainly interested to hear your thoughts.

Perhaps the best way to proceed is to advertise loudly in the docs that you cannot expect any random algorithm to just work and that lots of testing is required.

@gwater
Copy link
Contributor Author

gwater commented Jul 18, 2018

Thank you pointing out the discussion in #2, I had not read that yet.

I understand that you want Intervals to be usable in 3rd party algorithms, like ForwardDiff and I fully agree that it is great that we can combine packages in julia to get awesome synergy. But that imho requires that everyone follows the specifications laid out in Base. Otherwise unpredictable behavior is unavoidable.

I find this troubling especially w.r.t. IntervalRootFinding. The package is incredibly useful; in fact I used it for a manuscript I recently submitted. When I wrote that manuscript I trusted the that

This package provides guaranteed methods for finding roots of functions […]

If my input function in some way interacts badly through IntervalArithmetic that guarantee becomes untenable. Even worse: Imagine if I had used the buggy version of StaticArrays when I was writing. Luckily I wasn't but I am horrified when I think of that. There was no way I could have noticed or checked for that obscure bug. I'm afraid it would be very difficult for me to trust any results based on IntervalArithmetic again, knowing what I know now.

Maybe the solution is a second type of Interval that truly acts as a drop-in, i.e. for set operations it promotes to Interval but for Base functions behaves like a number. It's not ideal to have two implementations but the two notions of sets/numbers need to be separable in some way.

@Kolaru
Copy link
Collaborator

Kolaru commented Jul 18, 2018

Perhaps the best way to proceed is to advertise loudly in the docs that you cannot expect any random algorithm to just work and that lots of testing is required.

This would not really solve the problem, as it would not give the ability to check if a given third party algorithm is suitable or not for the use of Interval. I fully agree with gwater, in particular on the fact that a middle ground could be to support both versions (third party friendly interval and set-like interval). Using a supplementary type parameter and some type magic it is possible to avoid duplicating most of the implementation

abstract type Behavior end
struct SetLike <: Behavior end
struct NumberLike <: Behavior end
const DefaultBehavior = NumberLike

abstract type AbstractInterval{T} <: Real end

struct Interval{T<:Real, B<: Behavior} <: AbstractInterval{T}
    lo :: T
    hi :: T
end

# Two very naive constructor for the example.
Interval(a::T, b::T) where {T<:Real} = Interval{T, DefaultBehavior}(a, b)
Interval(::Type{B}, a::T, b::T) where {B <: Behavior, T <: Real} = Interval{T, B}(a, b)

function ==(a::Interval{T1, SetLike}, b::Interval{T2, SetLike}) where {T1 <: Real, T2 <: Real}
     isempty(a) && isempty(b) && return true
     a.lo == b.lo && a.hi == b.hi
end

function ==(a::Interval{T1, NumberLike}, b::Interval{T2, NumberLike}) where {T1 <: Real, T2 <: Real}
    isthin(a) && isequal(a, b) && return true
    isempty(a  b) && return false
    throw(UndecidableError("may represent equal numbers (or not)"))
end

The current type system allows to omit the trailing type parameters, which implies that here

Interval{T} == Interval{T, B} where {B <: Behavior}

And so the current syntax could still be used in methods definitions where it does not matter. Ideally the user should be able to choose the default behavior (it is perhaps possible to do similar to how the path to the python interpreter is passed to PyCall).

@dpsanders
Copy link
Member

@gwater I the interactions you refer to are very nasty and unexpected, and currently unfindable. So I guess throwing an error is a solution.

@Kolaru That's a very interesting idea / solution.

@gwater
Copy link
Contributor Author

gwater commented Jul 19, 2018

@Kolaru Thank you for providing that incredibly elegant solution. It really would allow most code to be shared, leaving only the handful of functions I modified to be re-implemented (plus figuring out the implications for various constructors).

@gwater
Copy link
Contributor Author

gwater commented Jul 25, 2018

I'm trying to implement your idea @Kolaru and I'm running into errors like this

MethodError: no method matching IntervalArithmetic.Interval{BigFloat,B} where B<:IntervalArithmetic.Behavior(::BigFloat, ::BigFloat)

basically anytime the constructor is called like Interval{T}(a, b). Any ideas why it's not picking up the correct constructor?

@Kolaru
Copy link
Collaborator

Kolaru commented Jul 26, 2018

@gwater I think that you need to explicitly implement the general method Interval{T}(a, b), as

Interval{T}(a, b) where {T <: Real} = Interval{T, DefaultBehavior}(a, b)

@Kolaru
Copy link
Collaborator

Kolaru commented Oct 22, 2018

@gwater Could you push your changes so we can see where the problem comes from ? I tried my previous suggestion and I can't reproduce your error.

I use the following definition of the Interval type:

struct Interval{T<:Real, B<:Behavior} <: AbstractInterval{T}
    lo :: T
    hi :: T

    function Interval{T, B}(a::Real, b::Real) where {T<:Real, B<:Behavior}
        if validity_check && !is_valid_interval(a, b)
            throw(ArgumentError("Interval of form [$a, $b] not allowed. Must have a ≤ b to construct Interval(a, b)."))
        end

        new(a, b)
    end
end

@gwater
Copy link
Contributor Author

gwater commented Oct 22, 2018

Sorry, I cannot work out how to get the current master to work with 0.6 or 1.0 (because there are no .toml files even though 0.7 is required). Basically, julia refuses constructors like

Interval{T}(a::T, b::T) = Interval{T, DefaultBehavior}(a, b)

Pre-compilation crashes with a stack overflow.

@dpsanders
Copy link
Member

Support for Julia 0.6 has been dropped. On Julia 1.0, just do

] dev IntervalArithmetic

which will check out git master into a local directory that you can edit.

Interval{T}(a::T, b::T) = Interval{T, DefaultBehavior}(a, b)

If that is an outer constructor, then (I believe you need) the new syntax:

Interval(a::T, b::T) where {T} = Interval{T, DefaultBehavior}(a, b)

@gwater
Copy link
Contributor Author

gwater commented Oct 22, 2018

Thanks for the help. So in Julia 1.0 try the following (on this branch):

Interval{Float64}(1., 2.)

and you get

ERROR: MethodError: no method matching Interval{Float64,B} where B<:IntervalArithmetic.Behavior(::Float64, ::Float64)
Closest candidates are:
  Interval{Float64,B} where B<:IntervalArithmetic.Behavior(::Any) where T at /home/jgraw/.julia/dev/IntervalArithmetic/src/intervals/intervals.jl:68
Stacktrace:
 [1] top-level scope at none:0

The branch currently has the constructor which was suggested.

@gwater
Copy link
Contributor Author

gwater commented Oct 22, 2018

I think before we put more time into this solution we should decide #237 though. Because if we agree that Intervals are not Numbers this pull request is obsolete.

@lbenet
Copy link
Member

lbenet commented Oct 22, 2018

I think the error is related to the need of an outer constructor for Interval{T,B}.

@dpsanders
Copy link
Member

dpsanders commented Oct 22, 2018 via email

@dpsanders
Copy link
Member

dpsanders commented Oct 22, 2018 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants