juliaintervals / intervalarithmetic.jl Goto Github PK
View Code? Open in Web Editor NEWLibrary for validated numerics using interval arithmetic
Home Page: https://juliaintervals.github.io/IntervalArithmetic.jl/
License: Other
Library for validated numerics using interval arithmetic
Home Page: https://juliaintervals.github.io/IntervalArithmetic.jl/
License: Other
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 tight conversion as
convert(Interval, x::Float64, Val{:tight})
that does what convert
does now.
Make ^
use the fast version and pow
the slow but precise version using BigFloat
(i.e., invert what currently happens).
Add version of mid
that just does 0.5*(a.lo + a.hi)
for speed, assuming the interval is bounded.
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.
This requires fixing AdjacentFloats.jl for Float32.
cf JuliaIntervals/IntervalRootFinding.jl#8
The @intervalbox
macro is used to produce versions used by bisection
and ForwardDiff
for root finding.
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.
cf #27
Base.iszero(a::Interval) = iszero(a.lo) && iszero(a.hi)
This is much faster than a == zero(a)
.
E.g. for multiplication.
julia> @interval "1e-400"
[0, 0]
julia> I"1e-400"
[0, 0]
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
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.
Check if the ambiguities are still a problem.
Move the dependency to IntervalRootFinding?
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.
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?
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?)julia> setrounding(Interval, :none)
:none
julia> x = 0.5..0.5
[0.499999, 0.500001]
julia> sin(x)
Interval(0.47942553860420295, 0.4794255386042031)
To show the current interval rounding mode in the fast_rounding
branch.
julia> x = -∞..∞
[-∞, ∞]
julia> pow(x, 2)
[NaN, ∞]
E.g. just take prevfloat
and nextfloat
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.
To make n-dimensional cubic box.
Do rounding down and up separately, not as an interval.
See #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?
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
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 future bisection
These should give the tightest possible interval.
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
.
For checking if an interval is atomic, ie unable to be split further.
pow(x::Interval, a::Real) = exp(a * log(x))
The former names should be deprecated.
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
cf #51
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.
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, ∞]
to compare bit strings.
e.g.
Base.isempty(x::Interval{Float64}) = x.lo === Inf && x.hi === -Inf
Currently in IntervalRootFinding.
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 Interval
s 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?
That do the rationalize trick.
This should halve the time required for convert(Interval, x)
for a float x
.
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
It is the entire real line.
cc @lbenet
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
.
julia> x = I"1e300"
[∞, ∞]
julia> @interval(1e300)
[∞, ∞]
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.