From 0588eedd7fec601ea9da0370417cc70b9b471417 Mon Sep 17 00:00:00 2001
From: kalmar <kalmar@amu.edu.pl>
Date: Mon, 13 Mar 2017 11:33:40 +0100
Subject: [PATCH] Using ValidatedNumerics (Interval Arithmetic)

---
 property(T).jl | 96 ++++++++++++++++++++++++++++++++++++++------------
 1 file changed, 74 insertions(+), 22 deletions(-)

diff --git a/property(T).jl b/property(T).jl
index 5b06ece..9c49150 100644
--- a/property(T).jl
+++ b/property(T).jl
@@ -4,6 +4,7 @@ import Base: rationalize
 using GroupAlgebras
 
 using ProgressMeter
+using ValidatedNumerics
 
 function create_product_matrix(basis, limit)
     product_matrix = zeros(Int, (limit,limit))
@@ -125,7 +126,14 @@ end
 function check_solution{T<:Number}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T}; verbose=true, augmented=false)
     result = compute_SOS(sqrt_matrix, Δ)
     if augmented
-        @assert GroupAlgebras.ɛ(result) == 0//1
+        epsilon = GroupAlgebras.ɛ(result)
+        if isa(epsilon, Interval)
+            @assert 0 in epsilon
+        elseif isa(epsilon, Rational)
+            @assert epsilon == 0//1
+        else
+            warn("Does checking for augmentation has meaning for $(typeof(epsilon))?")
+        end
     end
     SOS_diff = EOI(Δ, κ) - result
 
@@ -136,56 +144,100 @@ function check_solution{T<:Number}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlge
         if augmented
             println("ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) = ", GroupAlgebras.ɛ(SOS_diff))
         else
-            ɛ_dist = Float64(round(GroupAlgebras.ɛ(SOS_diff),12))
-            println("ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ≈ $ɛ_dist")
+            ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
+            if typeof(ɛ_dist) <: Interval
+                 ɛ_dist = ɛ_dist.lo
+            end
+            @printf("ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ≈ %.10f\n", ɛ_dist)
         end
-        L₁_dist = Float64(round(eoi_SOS_L₁_dist, 12))
-        println("‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $L₁_dist")
+
+        L₁_dist = eoi_SOS_L₁_dist
+        if typeof(L₁_dist) <: Interval
+            L₁_dist = L₁_dist.lo
+        end
+        @printf("‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ≈  %.10f\n", L₁_dist)
     end
 
     distance_to_cone = κ - 2^2*eoi_SOS_L₁_dist
     return distance_to_cone
 end
 
-function rationalize{T<:Integer, S<:Real}(::Type{T},
+import ValidatedNumerics.±
+function (±)(X::AbstractArray, tol::Real)
+    r{T}(x::T) = ( x==zero(T) ? @interval(x) : x ± tol)
+    return r.(X)
+end
+
+(±)(X::GroupAlgebraElement, tol::Real) = GroupAlgebraElement(X.coefficients ± tol, X.product_matrix)
+
+function Base.rationalize{T<:Integer, S<:Real}(::Type{T},
     X::AbstractArray{S}; tol::Real=eps(eltype(X)))
     r(x) = rationalize(T, x, tol=tol)
     return r.(X)
-end;
+end
 
 ℚ(x, tol::Real) = rationalize(BigInt, x, tol=tol)
 
 
 function ℚ_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A;
-    tol=10.0^-7, verbose=true)
+    tol=1e-7, verbose=true, rational=false)
 
     isapprox(eigvals(A), abs(eigvals(A)), atol=tol) ||
         warn("The solution matrix doesn't seem to be positive definite!")
     @assert A == Symmetric(A)
     A_sqrt = real(sqrtm(A))
 
-    println("")
-    println("Checking in floating-point arithmetic...")
-    @time fp_distance = check_solution(κ, A_sqrt, Δ, verbose=verbose)
-    println("Floating point distance (to positive cone) ≈ $(Float64(trunc(fp_distance,8)))")
+    # println("")
+    # println("Checking in floating-point arithmetic...")
+    # @time fp_distance = check_solution(κ, A_sqrt, Δ, verbose=verbose)
+    # println("Floating point distance (to positive cone) ≈ $(Float64(trunc(fp_distance,8)))")
+    # println("-------------------------------------------------------------")
+    # println("")
+    #
+    # if fp_distance ≤ 0
+    #     return fp_distance
+    # end
+
+    println("Checking in interval arithmetic...")
+    A_sqrtᴵ = A_sqrt ± tol
+    κᴵ = κ ± tol
+    Δᴵ = Δ ± tol
+    @time Interval_distance = check_solution(κᴵ, A_sqrtᴵ, Δᴵ, verbose=verbose)
+    # @assert isa(ℚ_distance, Rational)
+    println("The actual distance (to positive cone) is contained in $Interval_distance")
     println("-------------------------------------------------------------")
     println("")
 
-    if fp_distance ≤ 0
-        return fp_distance
+    if Interval_distance.lo ≤ 0
+        return Interval_distance.lo
     end
 
-    println("Checking in rational arithmetic...")
-    κ_ℚ = ℚ(trunc(κ,Int(abs(log10(tol)))), tol)
-    A_sqrt_ℚ, Δ_ℚ = ℚ(A_sqrt, tol), ℚ(Δ, tol)
-    @time ℚ_distance = check_solution(κ_ℚ, A_sqrt_ℚ, Δ_ℚ, verbose=verbose)
-    @assert isa(ℚ_distance, Rational)
-    println("Rational distance (to positive cone) ≈ $(Float64(trunc(ℚ_distance,8)))")
+    println("Projecting columns of A_sqrt to the augmentation ideal...")
+    A_sqrt_ℚ = ℚ(A_sqrt, tol)
+    A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ)
+    κ_ℚ = ℚ(κ, tol)
+    Δ_ℚ = ℚ(Δ, tol)
+
+    A_sqrt_ℚ_augᴵ = A_sqrt_ℚ_aug ± tol
+    κᴵ = κ_ℚ ± tol
+    Δᴵ = Δ_ℚ ± tol
+    @time Interval_dist_to_Σ² = check_solution(κᴵ, A_sqrt_ℚ_augᴵ, Δᴵ, verbose=verbose, augmented=true)
+    println("The Augmentation-projected actual distance (to positive cone) is contained in $Interval_dist_to_Σ²")
     println("-------------------------------------------------------------")
     println("")
-    if ℚ_distance ≤ 0
-        return ℚ_distance
+
+    if Interval_dist_to_Σ².lo ≤ 0 || !rational
+        return Interval_dist_to_Σ².lo
+    else
+
+        println("Checking Projected SOS decomposition in exact rational arithmetic...")
+        @time ℚ_dist_to_Σ² = check_solution(κ_ℚ, A_sqrt_ℚ_aug, Δ_ℚ, verbose=verbose, augmented=true)
+        @assert isa(ℚ_dist_to_Σ², Rational)
+        println("Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))")
+        println("-------------------------------------------------------------")
+        return ℚ_dist_to_Σ²
     end
+end
 
 function pmΔfilenames(name::String)
     if !isdir(name)