fix minor bug in example 5.5 from Smiley et al.

[gwater]

Just tested if this example still produces wrong results and found a minor bug that prevents it from running at all. This fixes that minor bug. (#86 is unaffected and remains unresolved.)

[dpsanders]

Thanks!

Have you tried with the new version of Krawczyk from #114?

[dpsanders]

With that branch and after compilation, I get

julia> @time roots_found = roots(SmileyExample55.f, SmileyExample55.region, Krawczyk, 1e-5)
 27.608686 seconds (76.47 M allocations: 6.279 GiB, 3.94% gc time)
41-element Array{Root{IntervalBox{4,Float64}},1}:
 Root([-3.1416, -3.14159] × [-9.06133e-08, 9.1164e-08] × [-1.33338e-06, 1.31911e-06] × [-7.83925e-07, 7.80965e-07], :unique)
 Root([-3.1416, -3.14159] × [2.36554, 2.36555] × [-1.3218e-08, 1.32441e-08] × [-0.775949, -0.775948], :unique)
 Root([-0.303088, -0.303087] × [-0.575478, -0.575476] × [-0.435553, -0.435551] × [-0.575406, -0.575405], :unique)
 Root([-0.303088, -0.303087] × [0.575476, 0.575478] × [-0.435553, -0.435551] × [0.575405, 0.575406], :unique)
 Root([-1.5708, -1.57079] × [-1.5708, -1.57079] × [-5.77248e-08, 5.75057e-08] × [-3.53968e-08, 3.54698e-08], :unique)
 Root([3.14158, 3.1416] × [-2.36555, -2.36554] × [-2.59861e-06, 3.06745e-06] × [0.775947, 0.77595], :unique)
 Root([0.303086, 0.303088] × [0.575476, 0.575478] × [0.43555, 0.435553] × [0.575405, 0.575407], :unique)
 Root([-3.17064e-07, 3.20536e-07] × [3.14159, 3.1416] × [-2.06644e-06, 2.04596e-06] × [-1.25055e-06, 1.24564e-06], :unique)
 Root([1.57079, 1.5708] × [1.57079, 1.5708] × [-3.29788e-08, 3.28537e-08] × [-2.02226e-08, 2.02643e-08], :unique)
 Root([-2.09169, -2.09168] × [-3.1416, -3.14159] × [1.50876, 1.50877] × [-2.59695e-08, 2.58817e-08], :unique)
 ⋮
 Root([1.0992, 1.09921] × [-2.02344, -2.02343] × [-0.677698, -0.677696] × [-0.452582, -0.45258], :unique)
 Root([-2.04239, -2.04238] × [1.11815, 1.11816] × [-0.677698, -0.677697] × [-0.452582, -0.452581], :unique)
 Root([-1.09922, -1.0992] × [-2.02344, -2.02343] × [0.677694, 0.677701] × [-0.452584, -0.452578], :unique)
 Root([-2.83851, -2.8385] × [2.56611, 2.56612] × [0.435551, 0.435552] × [-0.575406, -0.575405], :unique)
 Root([-1.19037e-08, 1.19175e-08] × [-0.776045, -0.776044] × [-1.03626e-08, 1.0383e-08] × [-0.775949, -0.775948], :unique)
 Root([1.0499, 1.04991] × [-6.99808e-08, 6.99574e-08] × [1.50876, 1.50877] × [-3.09882e-08, 3.08836e-08], :unique)
 Root([1.57079, 1.5708] × [-1.5708, -1.57079] × [-4.05028e-08, 4.08435e-08] × [-2.49073e-08, 2.51223e-08], :unique)
 Root([2.04238, 2.04239] × [1.11815, 1.11816] × [0.677697, 0.677698] × [-0.452582, -0.452581], :unique)
 Root([-3.1416, -3.14159] × [-3.1416, -3.14159] × [-2.16098e-10, 2.16484e-10] × [-1.81881e-10, 1.82358e-10], :unique)
 Root([-1.5708, -1.57079] × [1.57079, 1.5708] × [-3.00507e-08, 3.03036e-08] × [-1.84798e-08, 1.86393e-08], :unique)

julia> all(isunique, roots_found)
true

which seems pretty good!

[dpsanders]

Note that I was also using the fastpowers branch of IntervalArithmetic.jl; without that branch, any calculation involving interval powers will be an order of magnitude slower.

[dpsanders]

The point is that Krawczyk never inverts an interval matrix, only a Float64 matrix. I am coming to the conclusion that Krawczyk thus has a lot of benefits over interval Newton for this kind of thing.

[gwater]

The point is that Krawczyk never inverts an interval matrix, only a Float64 matrix. I am coming to the conclusion that Krawczyk thus has a lot of benefits over interval Newton for this kind of thing.

I agree that might make Krawczyk more robust in practice but not because Newton is inherently less stable but rather because there currently isn't a robust implement of interval matrix inversion (see #86).

In any case, this pull request does not address any of those issues. It really just aims to update one line so the test can be run at all.

[dpsanders]

What type is x? I don't think the splat is a good idea.

[dpsanders]

I agree that there needs to be robust matrix inversion. This should be provided by interval Gauss-Seidel, but we need a method to be able to specify this in the API.

[gwater]

What type is x? I don't think the splat is a good idea.

x could be any AbstractVector, Tuple, or IntervalBox. Really, anything with 4 elements. The example still runs in about 10-20 minutes so the splat doesn't seem to hurt performance noticably. Still, we could write SVector(x[1], x[2], x[3], x[4]) instead (although it is somewhat less elegant).

[dpsanders]

If you just write f(x) = (g ∘ g)(x) - x instead (without the broadcast), does it work?

[gwater]

No, because if x is an IntervalBox f will return an IntervalBox (rather than a vector of Intervals) and everything breaks.

[dpsanders]

How about defining

g(x::T) where {T} = T(C1 * (x[3] - α * sin(x[1]) * cos(x[2])) + x[1],
               C2 * (x[4] - α * cos(x[1]) * sin(x[2])) + x[2],
               D1 * (x[3] - α * sin(x[1]) * cos(x[2])),
               D2 * (x[4] - α * cos(x[1]) * sin(x[2])))

This is definitely a current pain point in the package. @Kolaru had suggested that everything should be defined in terms of IntervalBox, which I (think I) agree with. Or else it should at least be rewritten to use IntervalBox internally.

[dpsanders]

This seems to need the extra definition

julia> IntervalArithmetic.IntervalBox{N,T}(x...) where {N,T} = IntervalBox{N,T}(x)

to be able to use IntervalBoxes.

[gwater]

How about defining

g(x::T) where {T} = T(C1 * (x[3] - α * sin(x[1]) * cos(x[2])) + x[1],
               C2 * (x[4] - α * cos(x[1]) * sin(x[2])) + x[2],
               D1 * (x[3] - α * sin(x[1]) * cos(x[2])),
               D2 * (x[4] - α * cos(x[1]) * sin(x[2])))

This is definitely a current pain point in the package. @Kolaru had suggested that everything should be defined in terms of IntervalBox, which I (think I) agree with. Or else it should at least be rewritten to use IntervalBox internally.

This function breaks when used with a standard Vector{Float64}. And more generally, it's not clear to me why we should expect users to define their functions in this unusual way. User functions can be expected to return a scalar or some sort of AbstractVector subtype and that should just work. All I'm proposing is to make the user function in this example consistent so that it reliably returns an SVector (like all other examples). I'm a little surprised a simple one-line fix is causing so much discussion.

[dpsanders]

Ah good point. Of course I'll merge the fix but it's bringing up a more general point.

[dpsanders]

Thanks!