update SL(3,Z) to recent interface

SL(3,Z).jl

@@ -1,3 +1,12 @@
+using JLD
+using JuMP
+import SCS: SCSSolver
+import Mosek: MosekSolver
+
+using Groups
+using ProgressMeter
+
+
 function SL₃ℤ_generatingset()
 
     function E(i::Int, j::Int, N::Int=3)
@@ -13,41 +22,43 @@ function SL₃ℤ_generatingset()
     return S
 end
 
-function generate_B₂_and_B₄(B₁)
+function prepare_Δ_sdp_constraints(identity, S)
+    @show length(S)
+
+    B₁ = vcat([identity], S)
     B₂ = products(B₁, B₁);
     B₃ = products(B₁, B₂);
     B₄ = products(B₁, B₃);
-
     @assert B₄[1:length(B₂)] == B₂
-    return B₂, B₄;
-end
 
-function prepare_Laplacian_and_constraints(S, identity)
-
-    B₂, B₄ = generate_B₂_and_B₄(vcat([identity], S))
     product_matrix = create_product_matrix(B₄,length(B₂));
     sdp_constraints = constraints_from_pm(product_matrix, length(B₄))
-    L_coeff = splaplacian_coeff(S, B₄);
+    L_coeff = splaplacian_coeff(S, B₂, length(B₄));
+    Δ = GroupAlgebraElement(L_coeff, product_matrix)
 
-    return GroupAlgebraElement(L_coeff, product_matrix), sdp_constraints
+    return Δ, sdp_constraints
 end
 
-function prepare_Δ_sdp_constraints(name::String;cached=true)
-    f₁ = isfile("$name.product_matrix")
-    f₂ = isfile("$name.delta.coefficients")
+function load_Δ_sdp_constraints(name::String;cached=true)
+    pm_filename = "$name.product_matrix.jld"
+    Δ_coeff_filename = "$name.delta.coefficients.jld"
+    f₁ = isfile(pm_filename)
+    f₂ = isfile(Δ_coeff_filename)
     if cached && f₁ && f₂
         println("Loading precomputed pm, Δ, sdp_constraints...")
-        product_matrix = readdlm("$name.product_matrix", Int)
-        L = readdlm("$name.delta.coefficients")[:, 1]
-        Δ = GroupAlgebraElement(L, product_matrix)
+        product_matrix = load(pm_filename, "pm")
+        L = load(Δ_coeff_filename, "Δ")[:, 1]
+        Δ = GroupAlgebraElement(L, Array{Int,2}(product_matrix))
         sdp_constraints = constraints_from_pm(product_matrix)
     else
         println("Computing pm, Δ, sdp_constraints...")
         ID = eye(Int, 3)
-        S₁ = SL₃ℤ_generatingset()
-        Δ, sdp_constraints = prepare_Laplacian_and_constraints(S₁, ID)
-        writedlm("$name.delta.coefficients", Δ.coefficients)
-        writedlm("$name.product_matrix", Δ.product_matrix)
+        S = SL₃ℤ_generatingset()
+        Δ, sdp_constraints = prepare_Δ_sdp_constraints(ID, S)
+
+        save(pm_filename, "pm", Δ.product_matrix)
+        save(Δ_coeff_filename, "Δ", Δ.coefficients)
+
     end
     return Δ, sdp_constraints
 end
@@ -55,10 +66,10 @@ end
 
 function compute_κ_A(name::String, Δ, sdp_constraints;
     cached = true,
-    tol = TOL,
-    verbose = VERBOSE,
-    solver = MosekSolver(INTPNT_CO_TOL_REL_GAP=tol, QUIET=!verbose))
-    # solver = SCSSolver(eps=TOL, max_iters=ITERATIONS, verbose=VERBOSE))
+    tol = 1e-7,
+    verbose = false,
+    # solver = MosekSolver(INTPNT_CO_TOL_REL_GAP=tol, QUIET=!verbose))
+    solver = SCSSolver(eps=tol, max_iters=20000, cg_rate=3, verbose=verbose))
 
     f₁ = isfile("$name.kappa")
     f₂ = isfile("$name.SDPmatrixA")
@@ -70,31 +81,42 @@ function compute_κ_A(name::String, Δ, sdp_constraints;
     else
         println("Solving SDP problem maximizing κ...")
         κ, A = solve_SDP(sdp_constraints, Δ, solver, verbose=verbose)
-        writedlm("$name.kappa", kappa)
-        writedlm("$name.SDPmatrixA", A)
+        # writedlm("$name.kappa", kappa)
+        # writedlm("$name.SDPmatrixA", A)
     end
     return κ, A
 end
 
+function main()
+    const NAME = "SL3Z"
+    const VERBOSE = true
+    const TOL=1e-7
+    const Δ, sdp_constraints = load_Δ_sdp_constraints(NAME)
+    const κ, A = compute_κ_A(NAME, Δ, sdp_constraints, cached=false, verbose=VERBOSE)
+
+    if maximum(A) < 1e-2
+        warn("Solver might not solved the problem successfully and the positive solution is due to floating-point error, proceeding anyway...")
+    end
+
+    if κ > 0
+        @assert A == Symmetric(A)
+        const A_sqrt = real(sqrtm(A))
+
+        T = ℚ_distance_to_positive_cone(Δ, κ, A, tol=TOL, verbose=VERBOSE)
+
+        if T < 0
+            println("$NAME HAS property (T)!")
+        else
+            println("$NAME may NOT HAVE property (T)!")
+        end
+
+    else
+        println("$κ < 0: $NAME may NOT HAVE property (T)!")
+    end
+end
+
 @everywhere push!(LOAD_PATH, "./")
 using GroupAlgebras
 @everywhere include("property(T).jl")
 
-const NAME = "SL3Z"
-const VERBOSE = true
-const TOL=1e-7
-const Δ, sdp_constraints = prepare_Δ_sdp_constraints(NAME)
-const κ, A = compute_κ_A(NAME, Δ, sdp_constraints)
-
-if κ > 0
-    @time T = ℚ_distance_to_positive_cone(Δ, κ, A, tol=TOL, verbose=VERBOSE)
-
-    if T < 0
-        println("$NAME HAS property (T)!")
-    else
-        println("$NAME may NOT HAVE property (T)!")
-    end
-
-else
-    println("$κ < 0: $NAME may NOT HAVE property (T)!")
-end
+main()