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

Refactor with_interval_precision -> setprecision etc. #102

Merged
merged 5 commits into from
Apr 8, 2016
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/root_finding.md
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ julia> roots = newton(f, @interval(-5, 5))
Root([-1.4142135623730951, -1.414213562373095], :unique)
Root([1.414213562373095, 1.4142135623730951], :unique)

julia> set_interval_precision(256)
julia> setprecision(Interval, 256)
256

julia> newton(f, roots)
Expand Down
14 changes: 7 additions & 7 deletions docs/usage.md
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,7 @@ By default, the `@interval` macro creates intervals of `Float64`s.
This may be changed using the `set_interval_precision` function:

```julia
julia> set_interval_precision(256)
julia> setprecision(Interval, 256)
256

julia> @interval 3π/2 + 1
Expand All @@ -192,7 +192,7 @@ The subscript `256` at the end denotes the precision.

To change back to `Float64`s, use
```julia
julia> set_interval_precision(Float64)
julia> setprecision(Interval, Float64)
Float64

julia> @interval(pi)
Expand All @@ -201,7 +201,7 @@ julia> @interval(pi)

To check which mode is currently set, use
```julia
julia> get_interval_precision()
julia> precision(Interval)
(Float64,-1)
```
The result is a tuple of the type (currently `Float64` or `BigFloat`) and the precision (relevant only for `BigFloat`s).
Expand All @@ -223,7 +223,7 @@ Again, the result should contain the result of applying the function to each rea
Currently, this may be correctly calculated by using `BigFloat`s with a precision of 53 bits (the same as that of `Float64`s):

```julia
julia> set_interval_precision(53)
julia> setprecision(Interval, 53)
53

julia> sin(@interval(1))
Expand All @@ -239,16 +239,16 @@ Note, however, that calculations with `BigFloat`s are carried out in software, a
By default, the directed rounding used corresponds to using the `RoundDown` and `RoundUp` rounding modes when performing calculations; this gives the narrowest resulting intervals, and is set by

```julia
set_interval_rounding(:narrow)
setrounding(Interval, :narrow)
```

An alternative rounding method is to perform calculations using the (standard) `RoundNearest` rounding mode, and then widen the result by one machine epsilon in each direction using `prevfloat` and `nextfloat`. This is achived by
```julia
set_interval_rounding(:wide)
setrounding(Interval, :wide)
```
It generally results in wider intervals, but seems to be significantly faster.

The current interval rounding mode may be obtained by
```julia
get_interval_rounding()
rounding(Interval)
```
22 changes: 11 additions & 11 deletions examples/Roots of Wilkinson polynomials.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,7 @@
}
],
"source": [
"set_interval_precision(256)\n",
"setprecision(Interval, 256)\n",
"@time roots2 = newton(W₃, @interval(-10, 10))"
]
},
Expand Down Expand Up @@ -384,7 +384,7 @@
}
],
"source": [
"set_interval_precision(256)\n",
"setprecision(Interval, 256)\n",
"@time roots3 = newton(W₃, roots)"
]
},
Expand Down Expand Up @@ -431,8 +431,8 @@
}
],
"source": [
"set_interval_precision(Float64)\n",
"set_interval_rounding(:narrow)\n",
"setprecision(Interval, Float64)\n",
"setrounding(Interval, :narrow)\n",
"@time roots = newton(W₁₀, @interval(-20, 20))"
]
},
Expand Down Expand Up @@ -472,8 +472,8 @@
}
],
"source": [
"set_interval_precision(Float64)\n",
"set_interval_rounding(:wide)\n",
"setprecision(Interval, Float64)\n",
"setrounding(Interval, :wide)\n",
"@time roots = newton(W₁₀, @interval(-20, 20))"
]
},
Expand Down Expand Up @@ -520,7 +520,7 @@
}
],
"source": [
"set_interval_precision(256)\n",
"setprecision(Interval, 256)\n",
"@time roots2 = newton(W₁₀, @interval(-20, 20))"
]
},
Expand Down Expand Up @@ -567,7 +567,7 @@
}
],
"source": [
"set_interval_precision(256)\n",
"setprecision(Interval, 256)\n",
"@time roots3 = newton(W₁₀, roots)"
]
},
Expand Down Expand Up @@ -631,8 +631,8 @@
}
],
"source": [
"set_interval_precision(Float64)\n",
"set_interval_rounding(:wide)\n",
"setprecision(Interval, Float64)\n",
"setrounding(Interval, :wide)\n",
"@time roots = newton(W₁₂, @interval(-20, 20))\n",
"roots = ValidatedNumerics.clean_roots(roots)"
]
Expand Down Expand Up @@ -696,7 +696,7 @@
}
],
"source": [
"set_interval_precision(256)\n",
"setprecision(Interval, 256)\n",
"@time roots2 = newton(W₁₂, roots)\n",
"roots2 = ValidatedNumerics.clean_roots(roots2)"
]
Expand Down
25 changes: 16 additions & 9 deletions src/ValidatedNumerics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,6 @@ using CRlibm
using Compat
using FixedSizeArrays


setrounding(BigFloat, RoundNearest)
setrounding(Float64, RoundNearest)

import Base:
+, -, *, /, //, fma,
<, >, ==, !=, ⊆, ^, <=,
Expand All @@ -26,26 +22,34 @@ import Base:
⊆, eps,
floor, ceil, trunc, sign, round,
expm1, log1p,
precision,
isfinite, isnan


export
Interval,
@interval, @biginterval, @floatinterval, @make_interval,
get_interval_rounding, set_interval_rounding,
diam, radius, mid, mag, mig, hull, isinside,
emptyinterval, ∅, ∞, isempty, interior, isdisjoint, ⪽,
precedes, strictprecedes, ≺,
entireinterval, isentire, nai, isnai, isthin, iscommon,
widen, infimum, supremum,
set_interval_precision, get_interval_precision,
with_interval_precision,
parameters, eps, dist, roughly,
pi_interval,
midpoint_radius, interval_from_midpoint_radius,
RoundTiesToEven, RoundTiesToAway,
cancelminus, cancelplus, isunbounded,
.., @I_str

if VERSION >= v"0.5.0-dev+1182"
import Base: rounding, setrounding, setprecision
else
import Compat:
rounding, setrounding, setprecision

export rounding, setrounding, setprecision # reexport
end


## Multidimensional
export
Expand All @@ -60,8 +64,11 @@ export
find_roots_midpoint

function __init__()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice!

set_interval_precision(256) # set up pi
set_interval_precision(Float64)
setrounding(BigFloat, RoundNearest)
setrounding(Float64, RoundNearest)

setprecision(Interval, 256) # set up pi
setprecision(Interval, Float64)
end


Expand Down
4 changes: 2 additions & 2 deletions src/intervals/macros.jl
Original file line number Diff line number Diff line change
Expand Up @@ -44,11 +44,11 @@ end

doc"""The `@round` macro creates a rounded interval according to the current
interval rounding mode. It is the main function used to create intervals in the
library (e.g. when adding two intervals, etc.). It uses the interval rounding mode (see get_interval_rounding())"""
library (e.g. when adding two intervals, etc.). It uses the interval rounding mode (see rounding(Interval))"""
macro round(T, expr1, expr2)
#@show "round", expr1, expr2
quote
mode = get_interval_rounding()
mode = rounding(Interval)

if mode == :wide #works with any rounding mode set, but the result will depend on the rounding mode
# we assume RoundNearest
Expand Down
43 changes: 22 additions & 21 deletions src/intervals/precision.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,43 +17,47 @@ const parameters = IntervalParameters()

doc"`big53` creates an equivalent `BigFloat` interval to a given `Float64` interval."
function big53(a::Interval{Float64})
x = with_interval_precision(53) do # precision of Float64
x = setprecision(Interval, 53) do # precision of Float64
convert(Interval{BigFloat}, a)
end
end


set_interval_precision(::Type{Float64}) = parameters.precision_type = Float64
setprecision(::Type{Interval}, ::Type{Float64}) = parameters.precision_type = Float64
# does not change the BigFloat precision


function set_interval_precision{T}(::Type{T}, precision::Integer=256)
function setprecision{T<:AbstractFloat}(::Type{Interval}, ::Type{T}, prec::Integer)
#println("SETTING BIGFLOAT PRECISION TO $precision")
setprecision(BigFloat, precision)
setprecision(BigFloat, prec)

parameters.precision_type = T
parameters.precision = precision
parameters.precision = prec
parameters.pi = convert(Interval{BigFloat}, pi)

precision
prec
end

function with_interval_precision(f::Function, precision::Integer=256)
old_interval_precision = get_interval_precision()
#@show old_interval_precision
set_interval_precision(precision)
setprecision{T<:AbstractFloat}(::Type{Interval{T}}, prec) = setprecision(Interval, T, prec)

setprecision(::Type{Interval}, prec::Integer) = setprecision(Interval, BigFloat, prec)

function setprecision(f::Function, ::Type{Interval}, prec::Integer)

old_precision = precision(Interval)
setprecision(Interval, prec)

try
return f()
finally
set_interval_precision(old_interval_precision)
setprecision(Interval, old_precision)
end
end

set_interval_precision(precision) = set_interval_precision(BigFloat, precision)
set_interval_precision(t::Tuple) = set_interval_precision(t...)
# setprecision(::Type{Interval}, precision) = setprecision(Interval, precision)
setprecision(::Type{Interval}, t::Tuple) = setprecision(Interval, t...)

get_interval_precision() = (parameters.precision_type, parameters.precision)
#parameters.precision_type == Float64 ? (Float64, -1) : (BigFloat, parameters.precision)
precision(::Type{Interval}) = (parameters.precision_type, parameters.precision)


const float_interval_pi = convert(Interval{Float64}, pi) # does not change
Expand Down Expand Up @@ -88,17 +92,14 @@ else
end





float(x::Interval) =
# @round(BigFloat, convert(Float64, x.lo), convert(Float64, x.hi))
convert(Interval{Float64}, x)

## Change type of interval rounding:


doc"""`get_interval_rounding()` returns the current interval rounding mode.
doc"""`rounding(Interval)` returns the current interval rounding mode.
There are two possible rounding modes:

- :narrow -- changes the floating-point rounding mode to `RoundUp` and `RoundDown`.
Expand All @@ -108,9 +109,9 @@ This gives the narrowest possible interval.
`prevfloat` and `nextfloat` to achieve directed rounding. This creates an interval of width 2`eps`.
"""

get_interval_rounding() = parameters.rounding
rounding(::Type{Interval}) = parameters.rounding

function set_interval_rounding(mode)
function setrounding(::Type{Interval}, mode)
if mode ∉ [:wide, :narrow]
throw(ArgumentError("Only possible interval rounding modes are `:wide` and `:narrow`"))
end
Expand Down
6 changes: 3 additions & 3 deletions src/intervals/special.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ that this interval is an exception to the fact that the lower bound is
larger than the upper one."""
emptyinterval{T<:Real}(::Type{T}) = Interval{T}(Inf, -Inf)
emptyinterval{T<:Real}(x::Interval{T}) = emptyinterval(T)
emptyinterval() = emptyinterval(get_interval_precision()[1])
emptyinterval() = emptyinterval(precision(Interval)[1])
const ∅ = emptyinterval(Float64)

isempty(x::Interval) = x.lo == Inf && x.hi == -Inf
Expand All @@ -19,7 +19,7 @@ const ∞ = Inf
doc"""`entireinterval`s represent the whole Real line: [-∞, ∞]."""
entireinterval{T<:Real}(::Type{T}) = Interval{T}(-Inf, Inf)
entireinterval{T<:Real}(x::Interval{T}) = entireinterval(T)
entireinterval() = entireinterval(get_interval_precision()[1])
entireinterval() = entireinterval(precision(Interval)[1])

isentire(x::Interval) = x.lo == -Inf && x.hi == Inf
isunbounded(x::Interval) = x.lo == -Inf || x.hi == Inf
Expand All @@ -29,7 +29,7 @@ isunbounded(x::Interval) = x.lo == -Inf || x.hi == Inf
doc"""`NaI` not-an-interval: [NaN, NaN]."""
nai{T<:Real}(::Type{T}) = Interval{T}(NaN, NaN)
nai{T<:Real}(x::Interval{T}) = nai(T)
nai() = nai(get_interval_precision()[1])
nai() = nai(precision(Interval)[1])

isnai(x::Interval) = isnan(x.lo) || isnan(x.hi)

Expand Down
2 changes: 1 addition & 1 deletion src/root_finding/root_finding.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ function find_roots(f::Function, a::Real, b::Real, method::Function=newton;
tolerance=eps(1.0*a), debug=false, maxlevel=30, precision::Int=-1)

if precision >= 0
with_interval_precision(precision) do
setprecision(Interval, precision) do
find_roots(f, @interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
end

Expand Down
4 changes: 2 additions & 2 deletions test/ITF1788_tests/libieeep1788_tests_bool.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@ using ValidatedNumerics

#Preamble
setprecision(53)
set_interval_precision(Float64)
set_interval_rounding(:narrow)
setprecision(Interval, Float64)
setrounding(Interval, :narrow)

facts("minimal_empty_test") do
@fact isempty(∅) --> true
Expand Down
12 changes: 6 additions & 6 deletions test/ITF1788_tests/libieeep1788_tests_cancel.jl
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
#
# Copyright 2013 - 2015 Marco Nehmeier ([email protected])
# Copyright 2015 Oliver Heimlich ([email protected])
#
#
# Original author: Marco Nehmeier (unit tests in libieeep1788)
# Converted into portable ITL format by Oliver Heimlich with minor corrections.
#
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# http://www.apache.org/licenses/LICENSE-2.0
#
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
Expand All @@ -27,8 +27,8 @@ using ValidatedNumerics

#Preamble
setprecision(53)
set_interval_precision(Float64)
set_interval_rounding(:narrow)
setprecision(Interval, Float64)
setrounding(Interval, :narrow)

facts("minimal_cancelPlus_test") do
@fact cancelplus(Interval(-Inf, -1.0), ∅) --> entireinterval(Float64)
Expand Down
Loading