From a5c07a1b114a5326083cdfbf4320699e8f4c6529 Mon Sep 17 00:00:00 2001
From: Eeshan Gupta <eeshan9815@gmail.com>
Date: Sun, 10 Jun 2018 18:27:06 +0530
Subject: [PATCH] Inner iteration to improve choice of point of expansion

---
 src/linear_eq.jl |  5 +++++
 src/newtonnd.jl  | 35 +++++++++++++++++++++++++++--------
 2 files changed, 32 insertions(+), 8 deletions(-)

diff --git a/src/linear_eq.jl b/src/linear_eq.jl
index c5ec54ac..70bbd5e4 100644
--- a/src/linear_eq.jl
+++ b/src/linear_eq.jl
@@ -153,12 +153,17 @@ function linear_hull(M::AbstractMatrix, r::AbstractArray)
 	n = size(M, 1)
 
     ((M, r) = preconditioner(M, r))
+    M = interval.((2*eye(n) - sup.(M)), sup.(M))
     M_lo = inf.(M)
     M_hi = sup.(M)
+    P = eye(n)
     if all(.≤(M_lo, zero(M_lo)))
         return M \ r
     end
     P = inv(M_lo)
+    if any(isinf.(P)) || any(isnan.(P))
+        return gauss_seidel_interval(M, r, precondition=false, maxiter=2)
+    end
     if all(.≤(eye(n), (P)))
         H1 = P * sup.(r)
         C = 1 ./ (2diag(P) - 1)
diff --git a/src/newtonnd.jl b/src/newtonnd.jl
index 5a153080..0fed8055 100644
--- a/src/newtonnd.jl
+++ b/src/newtonnd.jl
@@ -1,10 +1,11 @@
+using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat
 """
 Preconditions the matrix A and b with the inverse of mid(A)
 """
 function preconditioner(A::AbstractMatrix, b::AbstractArray)
 
     Aᶜ = mid.(A)
-    B = pinv(Aᶜ)     # TODO If Aᶜ is singular
+    B = inv(Aᶜ)     # TODO If Aᶜ is singular
 
     return B*A, B*b
 
@@ -23,7 +24,10 @@ function newtonnd(f::Function, deriv::Function, x::IntervalBox{NUM, T}; reltol=e
         if !all(0 .∈ f(Xᴵ))
             continue
         end
+
         Xᴵ¹ = copy(Xᴵ)
+        use_B = false
+
         debug && (print("Current interval popped: ");  println(Xᴵ))
 
         if (isempty(Xᴵ))
@@ -45,22 +49,33 @@ function newtonnd(f::Function, deriv::Function, x::IntervalBox{NUM, T}; reltol=e
 
         while true
 
-            use_B = false
             next_iter = false
 
             initial_width = diam(Xᴵ)
             debug && (print("Current interval popped: ");  println(Xᴵ))
+            X = mid(Xᴵ)
             if use_B
-                # TODO Compute X using B in Step 19
-            else
-                X = mid(Xᴵ)
+                for i in 1:3
+                    z = X .- B * f(X)
+                    if all(z .∈ Xᴵ)
+                        if max(abs.(f(z))...) ≥ max(abs.(f(X))...)
+                            break
+                        end
+                        X = z
+                    else
+                        break
+                    end
+                end
+                if any(X .∉ Xᴵ)
+                    X = mid.(Xᴵ)
+                end
             end
 
             J = SMatrix{n}{n}(deriv(Xᴵ))   # either jacobian or calculated using slopes
-
-            # Xᴵ = IntervalBox((X + linear_hull(J, -f(interval.(X)))) .∩ Xᴵ)
+            B, r = preconditioner(J, -f(interval.(X)))
+            # Xᴵ = IntervalBox((X + linear_hull(B, r, precondition=false)) .∩ Xᴵ)
             Xᴵ = IntervalBox((X + (J \ -f(interval.(X)))) .∩ Xᴵ)
-
+            use_B = true
             if (isempty(Xᴵ))
                 next_iter = true
                 break
@@ -110,3 +125,7 @@ function newtonnd(f::Function, deriv::Function, x::IntervalBox{NUM, T}; reltol=e
 end
 
 newtonnd(f::Function, x::IntervalBox{NUM, T};  args...)  where {T<:AbstractFloat} where {NUM} = newtonnd(f, x->ForwardDiff.jacobian(f,x), x; args...)
+
+f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
+f(X) = f(X...)
+X = (0..1.23) × (0..1.123)