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

Reenable and fix example 5.5 from Smiley et al. #153

Merged
merged 7 commits into from
Jun 26, 2020
Merged
Show file tree
Hide file tree
Changes from 6 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
4 changes: 4 additions & 0 deletions .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,10 @@ env:
global:
- DOCUMENTER_DEBUG=true

script:
- if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
- travis_wait 20 julia --project -e 'using Pkg; Pkg.instantiate(); Pkg.build(; verbose = true); Pkg.test(; coverage=true)';

jobs:
allow_failures:
- julia: nightly
Expand Down
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "IntervalRootFinding"
uuid = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
version = "0.5.3"
version = "0.5.4"

[deps]
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
Expand Down
6 changes: 3 additions & 3 deletions src/contractors.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,10 @@ end
Multi-variable Newton operator.
"""
function 𝒩(f::Function, jacobian::Function, X::IntervalBox, α) # multidimensional Newton operator
m = IntervalBox(Interval.(mid(X, α)))
m = Interval.(mid(X, α))
J = jacobian(X)

return IntervalBox(m .- (J \ f(m)))
y = gauss_elimination_interval(J, f(m)) # J \ f(m)
return IntervalBox(m .- y)
end

"""
Expand Down
15 changes: 9 additions & 6 deletions src/linear_eq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ end
function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true)

x = similar(b)
x .= -1e16..1e16
x .= -Inf..Inf
x = gauss_elimination_interval!(x, A, b, precondition=precondition)

return x
Expand All @@ -103,10 +103,13 @@ function gauss_elimination_interval!(x::AbstractArray, A::AbstractMatrix, b::Abs

if precondition
(A, b) = preconditioner(A, b)
else
A = copy(A)
b = copy(b)
end
_A = A
_b = b
A = similar(A)
b = similar(b)
A .= _A
b .= _b
Copy link
Member

Choose a reason for hiding this comment

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

Is this different / better than _A = A; A = copy(_A)?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If I recall correctly, calling similar() on a StaticArray gives a MutableArray which is required for the setindex!() calls below.


n = size(A, 1)

Expand Down Expand Up @@ -162,5 +165,5 @@ function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::Ab
a \ b
end

\(A::StaticMatrix{Interval{T}}, b::StaticArray{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)
\(A::Matrix{Interval{T}}, b::Array{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)
\(A::SMatrix{S, S, Interval{T}}, b::SVector{S, Interval{T}}; kwargs...) where {S, T} = gauss_elimination_interval(A, b, kwargs...)
\(A::Matrix{Interval{T}}, b::Vector{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)
2 changes: 1 addition & 1 deletion test/roots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ end

# Bisection
rts = roots(f, X, Bisection, 1e-3)
@test length(rts) == 4
Copy link
Member

Choose a reason for hiding this comment

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

I do get 4 boxes here (even though there are only 2 roots).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, it gives an incorrect answer. Marking it as @test_broken makes transparent that Bisection is unreliable and probably should never be used in practise. In addition, it will notify us if Bisection ever starts to give the correct result here. (I admit, this is not really related to the topic of this PR–I think I just noticed this inconsistency while I was working on it.)

Copy link
Member

@dpsanders dpsanders Jun 24, 2020

Choose a reason for hiding this comment

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

I see what you mean, but I don't agree that the result is unreliable. It is reliably giving you a set of boxes that are guaranteed to contain all of the roots in the original box, and correctly telling you that it is unable to guarantee whether or not there are actually roots there.

The test is testing what the function returns, not what we would like it to ideally return ;)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ok, that makes sense. But testing for the result 4 seems overly specific if any number above 2 would acceptable. Maybe it would be better to test that the exact solutions are contained in whatever boxes are returned? I've created a separate PR #156 for that.

@test_broken length(rts) == 2

for method in newtonlike_methods
deriv = xx -> ForwardDiff.jacobian(f, xx)
Expand Down
6 changes: 4 additions & 2 deletions test/test_smiley.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ function test_all_unique(xs)
end

const tol = 1e-6
const method = Newton # NOTE: Bisection method performs badly in all examples
for method in (Newton, Krawczyk) # NOTE: Bisection method performs badly in all examples

@info("Testing method $(method)")

Expand All @@ -23,7 +23,7 @@ const method = Newton # NOTE: Bisection method performs badly in all examples
# no reference data for roots given
end

for example in (SmileyExample52, SmileyExample54) #, SmileyExample55)
for example in (SmileyExample52, SmileyExample54)#, SmileyExample55)
@testset "$(example.title)" begin
roots_found = roots(example.f, example.region, method, tol)
@test length(roots_found) == length(example.known_roots)
Expand All @@ -35,3 +35,5 @@ for example in (SmileyExample52, SmileyExample54) #, SmileyExample55)
end
end
end

end # method loop