Complete rewrite using generic branch and prune with contractors (#24)

* Start refactoring using generic branch and prune

* Make generic newton_helper for Newton and Krawczyk

* Revert inadvertent modification to newton.jl

* Add branch_and_prune.jl

* Add bisection branch-and-prune example

* Refactor to include f explicitly

* Add explicit dimension to helpers

* Add 2D bisection example

* Multi-dim Newton working with some clumsiness

Due to the interaction of FowardDiff and StaticArrays.

JuliaDiff/ForwardDiff.jl#235

* Add 2D Newton example with required syntax for now

* Docstring improvements

* Use contractors

* Add fudged Newton contractor

* Add StaticArrays dependency

* Add support for bisection of complex function

* Remove Contractor from contractor object names

* Export Bisection and Newton; remove _helper function

* Fix Newton

* Remove typing for Root

* Cleanup and rewrite examples for new syntax

* Add references to test suites

* Export roots and invert order of f and X in argument list

* Use IntervalArithmetic.Interval explicitly

* Update examples (now roots.jl)

* Use SVector in multidim Newton

* Fix complex examples

* Update minimum Julia version to 0.6 in REQUIRE and .travis

* Fix examples.jl

* Bump IntervalArithmetic to 0.12 in REQUIRE

* Address problems in tests

* Allow failures for nightly in travis

* Allow root vector without :unique information. Add basic docstring for roots function

* Add test_suites.md file with a list of places with test functions

* Add tests of new root finders

* Import roots from Polynomials.jl as stopgap measure

* Fix typo

* Fix test

* Fix typo

Author: dpsanders

.travis.yml

@@ -7,9 +7,13 @@ os:
   - osx
 
 julia:
-  - release
+  - 0.6
   - nightly
 
+matrix:
+   allow_failures:
+     - julia: nightly
+
 notifications:
   email: false
 

REQUIRE

@@ -1,3 +1,5 @@
-julia 0.5
-IntervalArithmetic 0.9
+julia 0.6
+IntervalArithmetic 0.12
 ForwardDiff 0.3
+StaticArrays 0.5
+Polynomials

docs/test_suites.md

@@ -0,0 +1,19 @@
+
+# Test suites:
+
+- http://folk.uib.no/ssu029/Pdf_file/Testproblems/testprobRheinboldt03.pdf
+
+-  http://www.mat.univie.ac.at/~neum/glopt/test.html
+
+-  http://titan.princeton.edu/TestProblems/
+
+- http://www-sop.inria.fr/saga/POL/
+
+-  MINPACK benchmarks: https://github.com/JuliaNLSolvers/NLsolve.jl/blob/master/test/minpack.jl
+
+ and other files in NLsolve test suite
+
+
+# Papers
+
+- Testing unconstrained optimization software, Moré, Garbow, Hillstrom ACM Trans. Math. Soft. 7 (1), 17-41 (1981)

examples/complex_roots.jl

@@ -0,0 +1,83 @@
+
+using IntervalArithmetic, IntervalRootFinding, StaticArrays
+
+rts = roots(sin, -5..5)
+
+rts = roots(sin, -5..6, Bisection)
+rts = roots(sin, rts, Newton)
+
+# 2D:
+f(x, y) = SVector(x^2 + y^2 - 1, y - x)
+f(X) = f(X...)
+
+rts = roots(f, (-5..5) × (-5..5))
+rts = roots(f, rts, Bisection)
+
+# complex:
+x = -5..6
+Xc = Complex(x, x)
+f(z) = z^3 - 1
+
+rts = roots(f, Xc, Bisection)
+rts = roots(f, rts, Newton)
+rts = roots(f, Xc)
+
+# From R docs:
+
+# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden
+
+function g(x)
+    (x1, x2, x3) = x
+    SVector(    x1^2 + x2^2 + x3^2 - 1,
+                x1^2 + x3^2 - 0.25,
+                x1^2 + x2^2 - 4x3
+            )
+end
+
+X = (-5..5)
+@time rts = roots(g, X × X × X)
+
+
+
+
+
+h(xv) = ((x,y) = xv; SVector(2*x - y - exp(-x), -x + 2*y - exp(-y)))
+
+rts = roots(h, X × X, Bisection)
+rts = roots(h, rts, Newton)
+
+
+
+# Dennis-Schnabel:
+h(xv) = ((x, y) = xv; SVector(x^2 + y^2 - 2, exp(x - 1) + y^3 - 2))
+
+rts = roots(h, X × X, Bisection)
+rts = roots(h, rts, Newton)
+
+# Test suites:
+
+# http://folk.uib.no/ssu029/Pdf_file/Testproblems/testprobRheinboldt03.pdf
+
+# http://www.mat.univie.ac.at/~neum/glopt/test.html
+
+# http://titan.princeton.edu/TestProblems/
+
+# http://www-sop.inria.fr/saga/POL/
+
+##  MINPACK benchmarks: https://github.com/JuliaNLSolvers/NLsolve.jl/blob/master/test/minpack.jl
+
+rosenbrock(x, y) = SVector( 1 - x, 10 * (y - x^2) )
+rosenbrock(X) = rosenbrock(X...)
+
+rosenbrock(xx) = ( (x, y) = xx; SVector( 1 - x, 1000 * (y - x^2) ) )
+X = IntervalBox(-1e5..1e5, 2)
+
+rts = roots(rosenbrock, X)
+
+
+
+# and other files in NLsolve test suite
+
+
+
+# Testing unconstrained optimization software, Moré, Garbow, Hillstrom ACM Trans. Math. Soft. 7 (1), 17-41 (1981)

examples/roots.jl

@@ -6,24 +6,26 @@ module IntervalRootFinding
 
 using IntervalArithmetic
 using ForwardDiff
+using StaticArrays
+
+import Polynomials: roots
 
 ## Root finding
 export
-    newton, krawczyk,
     derivative, jacobian,  # reexport derivative from ForwardDiff
     Root, is_unique,
-    find_roots,
-    find_roots_midpoint,
-    bisection,
+    roots, find_roots,
     bisect
 
 import Base: ⊆, show
 
+const Interval = IntervalArithmetic.Interval
+
 const derivative = ForwardDiff.derivative
 const D = derivative
 
 # Root object:
-immutable Root{T<:Union{Interval,IntervalBox}}
+immutable Root{T}
     interval::T
     status::Symbol
 end
@@ -58,7 +60,8 @@ include("bisect.jl")
 include("bisection.jl")
 include("newton.jl")
 include("krawczyk.jl")
-
+include("complex.jl")
+include("branch_and_prune.jl")
 
 function find_roots{T}(f::Function, a::Interval{T}, method::Function = newton;
                     tolerance = eps(T), debug = false, maxlevel = 30)

src/IntervalRootFinding.jl

@@ -10,7 +10,8 @@ function bisection{T<:Union{Interval,IntervalBox}}(f, X::T; tolerance=1e-3, debu
 
     debug && @show X, image
 
-    if !(zero(X) ⊆ image)
+    #if !(zero(X) ⊆ image)
+    if !(contains_zero(image))
         return Root{typeof(X)}[]
     end