add script check_SAutF5.jl

Manifest.toml

@@ -2,13 +2,19 @@
 
 [[AbstractAlgebra]]
 deps = ["InteractiveUtils", "LinearAlgebra", "Markdown", "Random", "SparseArrays", "Test"]
-git-tree-sha1 = "56bdd9e9bb2cfdf0876d0846b337c8f395d9a4b7"
+git-tree-sha1 = "a45d75b3cad41fafe2642d0f64e17866873869c8"
 uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
-version = "0.4.6"
+version = "0.4.4"
 
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[BenchmarkTools]]
+deps = ["JSON", "Printf", "Statistics", "Test"]
+git-tree-sha1 = "5d1dd8577643ba9014574cd40d9c028cd5e4b85a"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "0.4.2"
+
 [[BinDeps]]
 deps = ["Compat", "Libdl", "SHA", "URIParser"]
 git-tree-sha1 = "12093ca6cdd0ee547c39b1870e0c9c3f154d9ca9"
@@ -68,9 +74,9 @@ version = "0.2.0"
 
 [[Compat]]
 deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f"
+git-tree-sha1 = "195a3ffcb8b0762684b6821de18f83a16455c6ea"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "2.1.0"
+version = "2.0.0"
 
 [[DataStructures]]
 deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"]
@@ -116,9 +122,9 @@ version = "0.2.0"
 
 [[FileIO]]
 deps = ["Pkg", "Random", "Test"]
-git-tree-sha1 = "da32159d4a2e526338506685e280e39ed2f18961"
+git-tree-sha1 = "c94b0787956629036fb2b20fccde9e52b89d079a"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.0.6"
+version = "1.0.5"
 
 [[ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"]
@@ -195,10 +201,10 @@ uuid = "1b4a561d-cfcb-5daf-8433-43fcf8b4bea3"
 version = "0.4.1"
 
 [[LibCURL]]
-deps = ["BinaryProvider", "Libdl", "Printf", "Test"]
+deps = ["BinaryProvider", "Compat", "Libdl", "Printf"]
-git-tree-sha1 = "d051c8057512ca38a273aaa514145a0b25f24d46"
+git-tree-sha1 = "6339c87cb76923a3cf947fcd213cbc364355c9c9"
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.5.0"
+version = "0.4.1"
 
 [[LibExpat]]
 deps = ["Compat"]
@@ -251,10 +257,10 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 version = "0.3.2"
 
 [[Nemo]]
-deps = ["AbstractAlgebra", "BinaryProvider", "InteractiveUtils", "Libdl", "LinearAlgebra", "Markdown", "Test"]
+deps = ["AbstractAlgebra", "InteractiveUtils", "Libdl", "LinearAlgebra", "Markdown", "Test"]
-git-tree-sha1 = "24b361effd69b4ca78014a7bab1b81d1de7f35c0"
+git-tree-sha1 = "a41d652f0eefc6e08956bde976a10fed3535db57"
 uuid = "2edaba10-b0f1-5616-af89-8c11ac63239a"
-version = "0.13.2"
+version = "0.12.3"
 
 [[OrderedCollections]]
 deps = ["Random", "Serialization", "Test"]
@@ -282,7 +288,7 @@ uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 
 [[PropertyT]]
 deps = ["AbstractAlgebra", "Dates", "GroupRings", "Groups", "IntervalArithmetic", "JLD", "JuMP", "LinearAlgebra", "Markdown", "MathProgBase", "Nemo", "Printf", "SparseArrays"]
-git-tree-sha1 = "1307df01fc484b44fd24328f30304912ebd9bf2a"
+git-tree-sha1 = "4c4a5086b0c4dee30bba27ddf42e7ebb403765c2"
 repo-rev = "master"
 repo-url = "https://github.com/kalmarek/PropertyT.jl"
 uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d"

Project.toml

@@ -1,23 +1,9 @@
-name = "1712.07167"
-uuid = "c0ba6078-cfcc-57dc-8f34-6df289db1335"
-authors = ["Marek Kaluba <kalmar@amu.edu.pl>"]
-version = "0.1.0"
-
 [deps]
-AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 GroupRings = "0befed6a-bd73-11e8-1e41-a1190947c2f5"
 Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
-Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 PropertyT = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
-[extras]
-SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Test", "SCS"]

bin/check_SpNs.jl

@@ -1,91 +0,0 @@
-using Pkg
-Pkg.activate(".")
-using Test
-
-using PropertyT
-using Nemo
-
-using PropertyT.LinearAlgebra
-using PropertyT.SparseArrays
-using PropertyT.JuMP
-
-using PropertyT.AbstractAlgebra
-using PropertyT.Groups
-using PropertyT.GroupRings
-
-using PropertyT.JLD
-
-
-# include("../1712.07167/Roots.jl")
-include(joinpath("..", "src", "SpNs.jl"))
-
-using .SpNs
-using .SpNs.Roots
-
-BLAS.set_num_threads(Threads.nthreads())
-ENV["OMP_NUM_THREADS"] = Threads.nthreads()
-
-const N = 2
-const RADIUS = 3
-const UPPER_BOUND = 0.801
-
-root_system = SpNs.gens_roots(N)
-S = [g.matrix for g in root_system]
-G = parent(S[1])
-autS = WreathProduct(PermGroup(2), PermGroup(N))
-
-using SCS
-with_SCS = with_optimizer(SCS.Optimizer, linear_solver=SCS.Direct, max_iters=400000, eps=1e-11, alpha=1.95, acceleration_lookback=1, warm_start=true)
-
-sett = PropertyT.Settings("Sp($(2N),Z)_r$RADIUS", G, S, autS, with_SCS;
-    radius=RADIUS, upper_bound=UPPER_BOUND, warmstart=true)
-
-PropertyT.print_summary(sett)
-
-fp = PropertyT.fullpath(sett)
-isdir(fp) || mkpath(fp)
-
-Δ = PropertyT.Laplacian(S, RADIUS)
-RG = parent(Δ)
-Δs = SpNs.laplacians(RG, root_system)
-
-Sq  = sum( Δα^2 for (α, Δα) in Δs );
-# Op  = sum( Δα * sum(Δβ for (β, Δβ) in Δs if isorthogonal(α, β)) for (α, Δα) in Δs );
-Adj = sum( Δα * sum(Δβ for (β, Δβ) in Δs if !isproportional(α, β)) for (α, Δα) in Δs );
-
-const elt = Adj
-
-solver_logfile(sett) = joinpath(PropertyT.fullpath(sett), "Adj_solver_$(PropertyT.Dates.now()).log")
-
-λ, P = PropertyT.approximate_by_SOS(sett, elt, Δ, 
-    solverlog=solver_logfile(sett))
-
-save(PropertyT.filename(sett, :solution), "λ", λ, "P", P)
-λ < 0 && @warn "Solver did not produce a valid solution!"
-
-Q = real(sqrt(P));
-certified_λ = PropertyT.certify_SOS_decomposition(elt, Δ, λ, Q, R=RADIUS)
-
-@info "The obtained SOS can be use to certify λ ≥" certified_λ
-
-
-
-
-# 
-# od = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
-# orbit_data = PropertyT.decimate(od);
-# 
-# using SCS
-# using JuMP
-# warmstart = nothing
-# 
-# UB = 0.88 # for elt = Δ²;
-# SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data, upper_bound=UB);
-# # SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data);
-# 
-# sett = PropertyT.Settings("Sp($(2N),Z)", G, S, autS, solver(2000, accel=20);
-# upper_bound=1.3, warmstart=false)
-# 
-# with_SCS = with_optimizer(SCS.Optimizer, linear_solver=SCS.Direct, max_iters=100000, eps=1e-12, alpha=1.95, acceleration_lookback=1, warm_start=true)
-# 
-# status, warmstart = PropertyT.solve(SDP_problem, with_SCS, warmstart)

check_SAutF5.jl

@@ -0,0 +1,41 @@
+using Pkg
+Pkg.activate(".")
+using Groups
+using GroupRings
+using PropertyT
+using SparseArrays
+using LinearAlgebra
+using IntervalArithmetic
+using JLD
+
+@show Threads.nthreads()
+BLAS.set_num_threads(Threads.nthreads());
+
+G = SAut(FreeGroup(5))
+pm = load("oSAutF5_r2/pm.jld", "pm");
+RG = GroupRing(G, pm)
+@info RG
+
+S_size = 80
+Δ_coeff = SparseVector(maximum(pm), collect(1:(1+S_size)), [S_size; -ones(S_size)])
+Δ = GroupRingElem(Δ_coeff, RG);
+Δ² = Δ^2;
+
+@info "Loading solution"
+λ₀ = load("oSAutF5_r2/1.3/lambda.jld", "λ")
+P₀ = load("oSAutF5_r2/1.3/SDPmatrix.jld", "P");
+
+@info "Taking square root of P"
+@time Q = real(sqrt(P₀));
+
+Q_aug, check_columns_augmentation = PropertyT.augIdproj(Interval, Q);
+@show check_columns_augmentation
+if !check_columns_augmentation
+  @warn "Columns of Q are not guaranteed to represent elements of the augmentation ideal!"
+end
+
+@info "Computing SOS decomposition"
+@time sos = PropertyT.compute_SOS(RG, Q_aug);
+
+residual = Δ² - @interval(λ₀)*Δ - sos;
+@show norm(residual, 1)

src/Roots.jl

@@ -1,59 +0,0 @@
-module Roots
-
-using StaticArrays
-using LinearAlgebra
-
-export Root, isproportional, isorthogonal, ~, ⟂
-
-abstract type AbstractRoot end
-
-struct Root{N} <: AbstractRoot
-    coord::SVector{N, Float64}
-    
-    Root(a::SVector{N}) where N = new{N}(a)
-end
-
-Root(a) = Root(SVector(a...))
-
-function Base.:(==)(r::Root{N}, s::Root{M}) where {M, N}
-    M == N || return false
-    r.coord == s.coord || return false
-    return true
-end
-
-Base.hash(r::Root, h::UInt) = hash(Root, hash(r.coord, h))
-
-Base.:+(r::Root{N}, s::Root{N}) where N = Root(r.coord+s.coord)
-Base.:-(r::Root{N}, s::Root{N}) where N = Root(r.coord-s.coord)
-Base.:-(r::Root{N}) where N = Root(-r.coord)
-
-Base.:*(a::Number, r::Root) = Root(a*r.coord)
-Base.:*(r::Root, a::Number) = a*r
-
-Base.length(r::Root) = norm(r, 2)
-
-LinearAlgebra.norm(r::Root, p::Real=2) = norm(r.coord, p)
-
-LinearAlgebra.dot(r::Root{N}, s::Root{N}) where N = dot(r.coord, s.coord)
-
-cos_angle(a,b) = dot(a, b) / (norm(a) * norm(b))
-
-function isproportional(α::Root{N}, β::Root{N}) where N 
-    return isapprox(abs(cos_angle(α, β)), 1.0, atol=eps(1.0))
-end
-
-Base.:~(α::Root{N}, β::Root{N}) where N = isproportional(α, β)
-
-function isorthogonal(α::Root{N}, β::Root{N}) where N
-    return isapprox(cos_angle(α, β), 0.0, atol=eps(1.0))
-end
-
-⟂(α::Root{N}, β::Root{N}) where N = isorthogonal(α, β)
-
-function Base.show(io::IO, r::Root{N}) where N
-    print(io, "$(r.coord): root of length $(norm(r))")
-end
-
-E(N, i::Integer) = Root(ntuple(k -> k == i ? 1 : 0, N))
-
-end # of module RootSystems

src/SpN_actions.jl

@@ -1,39 +0,0 @@
-function blkdiag(A, B)
-    new_size = (size(A,1) + size(B,1), size(A,2) + size(B,2))
-    res = zeros(eltype(A), new_size)
-    res[1:size(A,1), 1:size(A,2)] .= A
-    res[size(A,1)+1:end, size(A,2)+1:end] .= B
-    return res
-end
-
-function matrix_emb(n::DirectPowerGroupElem, p::perm)
-    Id = parent(n.elts[1])()
-    elt = Diagonal([(-1)^(el == Id ? 0 : 1) for el in n.elts])
-    elt = elt[:, p.d]
-    return blkdiag(elt, elt)
-end
-
-function (g::WreathProductElem)(A::MatElem)
-    g_inv = inv(g)
-    G = matrix_emb(g.n, g_inv.p)
-    G_inv = matrix_emb(g_inv.n, g.p)
-    M = parent(A)
-    return M(G)*A*M(G_inv)
-end
-
-function Base.:*(x::AbstractAlgebra.MatElem, P::Generic.perm)
-   z = similar(x)
-   m = nrows(x)
-   n = ncols(x)
-   for i = 1:m
-      for j = 1:n
-         z[i, j] = x[i,P[j]]
-      end
-   end
-   return z
-end
-
-function (p::perm)(A::MatElem)
-    length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
-    return p*A*inv(p)
-end

src/SpNs.jl

@@ -1,119 +0,0 @@
-module SpNs
-
-using LinearAlgebra
-using AbstractAlgebra
-using Nemo
-using Groups
-using GroupRings
-
-include("Roots.jl")
-using .Roots
-
-export Sp, isproportional, ~, isorthogonal, ⟂, positive
-
-# two families of generators
-
-function a(i::Int, j::Int, M::MatSpace, val=one(M.base_ring)) # short roots
-    @assert i ≠ j
-    m = one(M)
-    N = size(m,2)÷2
-    
-    m[i,j] = val
-    m[N+j, N+i] = -val
-    @assert issymplectic(m)
-    return m
-end
-
-function b(i::Int, j::Int, M::MatSpace, val=one(M.base_ring)) # long roots
-    m = one(M)
-    N = size(m,2)÷2
-    m[i, N+j] = val
-    m[j, N+i] = val
-    @assert issymplectic(m)
-    return m
-end
-
-# symplectic generator with associated root datum:
-
-struct Sp{N}
-    root::Root{N}
-    matrix::Nemo.fmpz_mat
-end
-
-function Sp{N}(typ::Val{:a}, i::Int, j::Int) where N
-    @assert i≠j
-    r = Roots.E(N, i) - Roots.E(N, j)
-    M = Nemo.MatrixSpace(Nemo.ZZ, 2N, 2N)
-    m = a(i,j, M)
-    return Sp{N}(r, m)
-end
-
-function Sp{N}(typ::Val{:b}, i::Int, j::Int=i) where N
-    r = Roots.E(N, i) + Roots.E(N, j)
-    M = Nemo.MatrixSpace(Nemo.ZZ, 2N, 2N)
-    m = b(i,j, M)
-    return Sp{N}(r, m)
-end
-
-Sp{N}(s::Symbol, i::Int, j::Int=i) where N = Sp{N}(Val(s), i, j)
-
-Base.inv(s::Sp) = Sp(s.root, inv(s.matrix))
-Base.:-(s::Sp) = Sp(-s.root, transpose(s.matrix))
-
-Roots.isproportional(x::Sp, y::Sp) = isproportional(x.root, y.root)
-Roots.:~(x::Sp,y::Sp) = x.root ~ y.root
-
-Roots.isorthogonal(x::Sp, y::Sp) = isorthogonal(x.root, y.root)
-Roots.:⟂(x::Sp,y::Sp) = x.root ⟂ y.root
-
-function issymplectic(M)
-    r = nrows(M)
-    c = ncols(M)
-    @assert r == c
-    @assert iseven(r)
-    n = r÷2
-    I = Diagonal(ones(Int, n))
-    
-    O = zeros(Int,n,n);
-    Ω = parent(M)([O  I; -I O])
-    return Ω == transpose(M)*Ω*M
-end
-
-indexing(n) = [(i,j) for i in 1:n for j in i+1:n]
-
-function gens_roots(N)
-    a_ijs = [Sp{N}(:a, i,j) for (i,j) in indexing(N)]
-    b_is = [Sp{N}(:b, i) for i in 1:N]
-    c_ijs = [Sp{N}(:b, i,j) for (i,j) in indexing(N)]
-    root_system = [a_ijs; b_is; c_ijs]
-    root_system = [root_system; [-r for r in root_system]];
-    root_system = [root_system; inv.(root_system)]
-    return root_system
-end
-
-function positive(root_system::Vector{SP}) where SP
-    roots = [r.root for r in root_system]
-    
-    pos_roots = eltype(roots)[]
-    for (idx,γ) in enumerate(roots)
-        if any(isproportional(s, γ) for s in roots[1:idx-1])
-            continue
-        end
-        push!(pos_roots, γ)
-        @info "" positive_root = γ
-    end
-    # TODO: what about their sums??
-    
-    return pos_roots
-end
-
-function laplacians(RG::GroupRing, root_system::Vector{SpNs.Sp{N}}) where N
-    positive_roots = positive(root_system)
-    Sαs = Dict(α => [w.matrix for w in root_system if isproportional(w.root, α)] for α in positive_roots)
-    Δs = Dict(α => RG(length(S)) - RG(S) for (α, S) in Sαs)
-    return Δs
-end
-
-include("SpN_actions.jl")
-
-end # of module SPNs

test/runtests.jl

@@ -1,108 +0,0 @@
-using Test
-
-using PropertyT
-include("../src/SpNs.jl")
-using .SpNs
-using .SpNs.Roots
-
-@testset "Roots" begin
-    @test Root((2,1)) isa Root
-    r = Root((2,1))
-    @test -r isa Root
-    @test 2r isa Root
-    @test length(2r) == 2length(r)
-
-    @test r+r isa Root
-    @test length(r+r) == 2length(r)
-    @test r-r isa Root;
-    @test length(r-r) == 0.0
-
-    r = Root((1,1));
-    s = Root((-1,1));
-    @test string(Root((1,1))) == "[1.0, 1.0]: root of length 1.4142135623730951"
-
-    @test isproportional(r, r)
-    @test isproportional(r, -r)
-    @test isproportional(r, 2r)
-    @test isproportional(-r, 2r)
-    @test isproportional(r+s, -s-r)
-
-    @test isorthogonal(r, s)
-    @test isorthogonal(r, -s)
-    @test isorthogonal(r, 2s)
-    @test isorthogonal(-r, 2s)
-
-    @test !isorthogonal(r, r+s)
-    @test !isorthogonal(r, -r)
-    @test !isorthogonal(r+s, 2s)
-    @test !isorthogonal(-r, 2s+r)
-
-    @test isorthogonal(r+s, -s+r)
-
-    N = 3
-    @test length(Roots.E(N, 1) - Roots.E(N,2)) == sqrt(2)
-end
-
-
-@testset "Symplectic Generators" begin
-
-    s = Sp{2}(:a, 1, 2) 
-    @test -s isa Sp 
-    @test inv(s) isa Sp
-    t = Sp{2}(:b, 1, 2)
-    @test -t isa Sp
-    @test inv(t) isa Sp
-
-    @test isproportional(s, -s)
-    @test isproportional(s, -s)
-    @test isproportional(s, inv(s))
-    
-    @test isproportional(t, -t)
-    @test isproportional(t, inv(t))
-    @test !isproportional(s, t)
-    @test !isproportional(inv(s), -t)
-    
-end
-
-@testset "Symplectic Group" begin
-
-    @testset "Symplectic Root System" begin
-        rs = SpNs.gens_roots(2)
-        S = [r.matrix for r in rs];
-        @test all(SpNs.issymplectic.(S))
-        @test length(S) == 16
-        
-        rs = SpNs.gens_roots(3)
-        S = [r.matrix for r in rs];
-        @test all(SpNs.issymplectic.(S))
-        @test length(S) == 36
-    end
-
-    @testset "Sp(4,Z)" begin
-
-        N = 2
-        root_system = SpNs.gens_roots(N)
-        S = [g.matrix for g in root_system]
-        
-        Δ = PropertyT.Laplacian(S, 2)
-        RG = parent(Δ)
-        
-        Δs = SpNs.laplacians(RG, root_system)
-        
-        @test sum(Δα for (α, Δα) in Δs) == Δ
-        
-        @testset "Orthogonality & commutation" begin
-            r = Root([2, 0])
-            s = Root([0, 2])
-
-            a = Δs[r]
-            b = Δs[s]
-            @test a*b == b*a
-        end
-        
-        Sq  = sum( Δα^2 for (α, Δα) in Δs );
-        Adj = sum( Δα * sum(Δβ for (β, Δβ) in Δs if !isproportional(α, β)) for (α, Δα) in Δs );
-        
-        @test Sq + Adj == Δ^2
-    end
-end