juliaintervals / intervalarithmetic.jl Goto Github PK

julia> tan(Interval(6.638314112824137, 8.38263151220128))
ERROR: ArgumentError: Must have a ≤ b to construct Interval(a, b).
Stacktrace:
 [1] IntervalArithmetic.Interval{Float64}(::Float64, ::Float64) at /Users/dpsanders/.julia/v0.6/IntervalArithmetic/src/intervals/intervals.jl:23
 [2] tan(::IntervalArithmetic.Interval{Float64}) at /Users/dpsanders/.julia/v0.6/IntervalArithmetic/src/intervals/trigonometric.jl:135

The checks according to quadrant are not correct.

cf #51

julia> a = (3..3) + (4..4)*im
[3, 3] + [4, 4]*im

julia> log(a)
ERROR: StackOverflowError:
Stacktrace:
 [1] log(::Complex{IntervalArithmetic.Interval{Float64}}) at ./complex.jl:463 (repeats 80000 times)

This probably requires a specialised function to be defined.
Maybe we should formalise complex intervals as a new type?

Function choices including:

• Current ˆ (tight, slow; rename to tight_power) vs. pow (less tight, fast; rename to fast_power) -- choose one of them to use as ^ and the other as pow

• convert(Interval, x::Float64) (tight vs. accurate)

• Fast versions of tanh etc. (all functions not exported by CRlibm, for which MPFR is currently used)

Syntax something like

IntervalArithmetic.configure(power=^, convert=:tight)

Can be a non-exported function.

In the branch simplify_constructor, the inner constructor of Interval, that performs checks for a valid interval, is removed.

The goal of this is performance: constructing an Interval goes down from 9ns to 1.5ns, and Intervals are constructed all the time in real code.

I propose an interval (with small i) function that carries out the current checks when creating an interval, that should be used by users. We could even (I believe) not export the Interval object itself!

Any comments, @lbenet?

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.

julia 0.6> setprecision(1000000)
1000000

julia 0.6> @biginterval sin(pi/6) - 0.5
[-5.05018e-3010, 1.01004e-30103]₁₀₀₀₀₀₀

The exponent should be around 300000

Don't use complicated conversions -- leave that for @interval.
Just do Interval(prevfloat(a), nextfloat(b)).

But can treat integers exactly:

..(a::Integer, b::Integer) = Interval(a, b)
..(a::Integer, b::Float64) = Interval(a, nextfloat(b))
..(a::Float64, b::Int) = Interval(prevfloat(a), b)

..(a::Float64, b::Float64) = Interval(prevfloat(a), nextfloat(b))

This will take care of most of the useful cases, such as writing 0..\infty.

Add tight conversion as

convert(Interval, x::Float64, Val{:tight})

that does what convert does now.

Check if the ambiguities are still a problem.
Move the dependency to IntervalRootFinding?

This requires fixing AdjacentFloats.jl for Float32.

at https://github.com/JuliaIntervals/IntervalArithmetic.jl/blob/master/src/intervals/rounding.jl#L131
arithmetic operations are defined for all T<:AbstractFloats, but prev_float (in AdjacentFloats.jl) is defined only for Float{16, 32, 64}.
This gets emanated in

ERROR: LoadError: MethodError: no method matching prev_float(::BigFloat)
Closest candidates are:
  prev_float(!Matched::Float16) at /home/kalmar/.julia/v0.5/AdjacentFloats/src/AdjacentFloats.jl:35
  prev_float(!Matched::Float32) at /home/kalmar/.julia/v0.5/AdjacentFloats/src/AdjacentFloats.jl:31
  prev_float(!Matched::Float64) at /home/kalmar/.julia/v0.5/AdjacentFloats/src/AdjacentFloats.jl:21
  ...
 in *(::IntervalArithmetic.Interval{BigFloat}, ::IntervalArithmetic.Interval{BigFloat}) at /home/kalmar/.julia/v0.5/IntervalArithmetic/src/intervals/arithmetic.jl:84

(this used to work on v0.9.1)

Julia Version 0.5.2
Commit f4c6c9d4bb* (2017-05-06 16:34 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-6200U CPU @ 2.30GHz
  WORD_SIZE: 64
  BLAS: libopenblas (NO_LAPACK NO_LAPACKE NO_AFFINITY HASWELL)
  LAPACK: liblapack
  LIBM: libm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake-avx512)

IntervalArithmetic v0.10.0

We should (maybe?) be able to change the precision of a..b.
Currently it does Interval(prevfloat(a), nextfloat(b)) for performance, but we lose a modicum of precision.

The following functions are currently broken for complex intervals; see JuliaLang/julia#22095.
They give a stack overflow due to incorrect fallbacks in base Julia:

• exp, sin, cos, cosh, sinh -- due to sincos
• sqrt, log, log1p, log2, log10
• asin, acos, atan
• tanh, asinh, acosh, atanh
• sinpi, cospi, sinc, cosc

While IntervalArithmetic successfully precompiles on Julia 0.7/1.0 as of yesterday afternoon (Julia 0.7.0-DEV.3222), any attempt to print an Interval object currently leads to a segfault, at least on a recent Intel/Darwin setup.

I think I've the issue to the utility function "round_string" which is used to round for pretty printing. Since latest commit of "IntervalArithmetic.jl/src/display.jl" (96e9d80) the definition is

# round to given number of signficant digits
# basic structure taken from string(x::BigFloat) in base/mpfr.jl
function round_string(x::BigFloat, digits::Int, r::RoundingMode)

    lng = digits + Int32(8)
    buf = Array{UInt8}(lng + 1)

    lng = ccall((:mpfr_snprintf,:libmpfr), Int32,
    (Ptr{UInt8}, Culong,  Ptr{UInt8}, Int32, Ptr{BigFloat}...),
    buf, lng + 1, "%.$(digits)R*g", Base.MPFR.to_mpfr(r), &x)

    repr = unsafe_string(pointer(buf))

    repr = replace(repr, "nan", "NaN")

    return repr
end

The syntax for passing references to C functions, as is used here to operate directly on a BigFloat object using the MPFR C library function snprintf, changes in 0.7/1.0. While the old syntax ought to still work, in round_string changing the syntax as follows appears to avoid the segfault:

function round_string(x::BigFloat, digits::Int, r::RoundingMode)

    lng = digits + Int32(8)
    buf = Array{UInt8}(lng + 1)

    lng = ccall((:mpfr_snprintf,:libmpfr), Int32,
    (Ptr{UInt8}, Culong,  Ptr{UInt8}, Int32, Ref{BigFloat}...),
    buf, lng + 1, "%.$(digits)R*g", Base.MPFR.to_mpfr(r), x)

    repr = unsafe_string(pointer(buf))

    repr = replace(repr, "nan", "NaN")

    return repr
end

In the docs:

Due to the way floating-point arithmetic works, the interval
`a` created directly by the constructor turns out to contain
*neither the true real number 0.1, nor 0.3*, since the floating point
number associated to 0.1 is actually rounded up, whereas the
one associated to 0.3 is rounded down.
The [`@interval`](@ref) macro, however, uses [**directed rounding**](rounding.md) to *guarantee*
that the true 0.1 and 0.3 are included in the result.

But the interval is [0.0999999, 0.300001]: doesn't that contain the values?

julia> Interval(prev_float(Inf), Inf) + Interval(prev_float(Inf), Inf)
ERROR: ArgumentError: Must satisfy `a < Inf` and `b > -Inf` to construct Interval(a, b).

It's trying to construct Interval(Inf, Inf).

This seems to be a bug with FastRounding, since

julia 0.6> setrounding(Interval, :slow)
:slow

julia 0.6> Interval(prevfloat(Inf), Inf) + Interval(prevfloat(Inf), Inf)
[1.79769e+308, ∞]

That do the rationalize trick.
This should halve the time required for convert(Interval, x) for a float x.

I suggest that we should have the following mechanisms for constructing intervals:

• interval(a, b) (with a small i): the current Interval(a, b) -- includes the checks at the moment of construction; constructs an interval from exactly the given value of a to the given value of b; fast

• a..b: Gives a guaranteed enclosure of the interval, including both a and b; see #34; fast.
  Does not try to do anything clever with floats; just uses prevfloat and nextfloat

• convert(Interval, a): Gives the smallest possible enclosure of the number a as an interval; does fancy things (rationalize or parse(string(x))) to interpret e.g. 0.1 as the real number 0.1; slow (same as currently)

• @interval(a, b): The most general method, but slowest.

• a ± b: Only for numbers a and b (not intervals)

cc @lbenet

For checking if an interval is atomic, ie unable to be split further.

Add

• parse(Interval, "1e-400")
• convert(Interval, "1e-400")

This is part of the string input and output required by the standard. We should check the details.

It is the entire real line.

cc @lbenet

The former names should be deprecated.

Currently in IntervalRootFinding.

These should give the tightest possible interval.

cf #27

to compare bit strings.

e.g.

Base.isempty(x::Interval{Float64}) = x.lo === Inf && x.hi === -Inf

Make ^ use the fast version and pow the slow but precise version using BigFloat (i.e., invert what currently happens).

For future bisection

cf JuliaIntervals/IntervalRootFinding.jl#8

The @intervalbox macro is used to produce versions used by bisection and ForwardDiff for root finding.

To show the current interval rounding mode in the fast_rounding branch.

Doing @btime on e.g. sin(0..1) has unexpected allocations (there should be no allocation).

We traced this back to convert, which in turn uses rationalize or string, both of which allocate.

We should definitely go back to having convert(Interval, x) do x..x, which does not allocate...

cc @lbenet

Add version of mid that just does 0.5*(a.lo + a.hi) for speed, assuming the interval is bounded.

Do rounding down and up separately, not as an interval.

To make n-dimensional cubic box.

See #51

http://www.mat.univie.ac.at/~dferi/publications.html

Base.iszero(a::Interval) = iszero(a.lo) && iszero(a.hi)

This is much faster than a == zero(a).
E.g. for multiplication.

I think the only sensible options for this package (the fundamental part of what was ValidatedNumerics.jl) are:

• IntervalArithmetic.jl (my preferred option)
• Intervals.jl (too generic?)
• IntervalMethods.jl
• IntervalAnalysis.jl (too broad?)

cc @lbenet, @JeffreySarnoff

E.g. just take prevfloat and nextfloat

Instead of using the BigFloat versions, use the definitions, e.g.

tanh = sinh / cosh

and https://en.wikipedia.org/wiki/Inverse_hyperbolic_function.

However, note that these will not, in general, be correctly rounded.

julia 0.6> @biginterval(10.0)^(10.0^100.0)
ERROR: InexactError()

This is due to the line

    isinteger(x) && return a^(round(Int, x))

Here, Int should be BigInt for BigFloat.