From db74cb2226d1f90fa3cd8103ba3db86e3b510d3c Mon Sep 17 00:00:00 2001 From: Daniel Brosch <73886037+DanielBrosch@users.noreply.github.com> Date: Wed, 7 Aug 2024 13:32:23 +0200 Subject: [PATCH 1/2] faster --- src/FlagAlgebras/AbstractFlagAlgebra.jl | 5 +- src/FlagAlgebras/BinaryTrees.jl | 4 +- src/FlagAlgebras/EdgeMarkedFlags.jl | 29 +- src/FlagAlgebras/FreeVariables.jl | 115 +++++ src/FlagAlgebras/Graphs.jl | 40 +- src/FlagAlgebras/ProductFlag.jl | 33 +- src/FlagAlgebras/QuantumFlags.jl | 1 + src/FlagAlgebras/SymmetricFunctions.jl | 10 +- src/FlagModels/FlagModel.jl | 7 +- src/FlagModels/LasserreModel.jl | 121 ++++- src/FlagModels/QuadraticModule.jl | 73 ++- src/FlagModels/RazborovModel.jl | 70 +-- src/utils/Nauty.jl | 578 ++++++++++++++---------- src/utils/NautyBackup.jl | 475 +++++++++++++++++++ src/utils/SchreierSims.jl | 80 ++-- test/src/nauty.jl | 79 ++++ 16 files changed, 1365 insertions(+), 355 deletions(-) create mode 100644 src/FlagAlgebras/FreeVariables.jl create mode 100644 src/utils/NautyBackup.jl diff --git a/src/FlagAlgebras/AbstractFlagAlgebra.jl b/src/FlagAlgebras/AbstractFlagAlgebra.jl index 022e6f5..6dfd1cc 100644 --- a/src/FlagAlgebras/AbstractFlagAlgebra.jl +++ b/src/FlagAlgebras/AbstractFlagAlgebra.jl @@ -285,7 +285,7 @@ end Returns the sub-Flag indexed by `vertices`, which is a subset of `1:size(F)`. """ function subFlag(F::T, vertices::AbstractVector{Int})::T where {T<:Flag} - error("subFlag is not defined for Flag type $T") + error("subFlag is not defined for Flag type T") return missing end @@ -336,7 +336,7 @@ function countEdges(F::T)::Vector{Int} where {T<:Flag} end """ - isolatedVertices(F::T)::Vector{Int} where{T<:Flag} + isolatedVertices(F::T)::BitVector where{T<:Flag} Returns the indicator vector of isolated vertices of `F`. """ @@ -467,4 +467,5 @@ include("EdgeMarkedFlags.jl") include("SymmetricFunctions.jl") include("HarmonicFlags.jl") include("ProductFlag.jl") +include("FreeVariables.jl") include("PartiallyColoredFlags.jl") diff --git a/src/FlagAlgebras/BinaryTrees.jl b/src/FlagAlgebras/BinaryTrees.jl index 34f3424..cf2aa92 100644 --- a/src/FlagAlgebras/BinaryTrees.jl +++ b/src/FlagAlgebras/BinaryTrees.jl @@ -829,7 +829,7 @@ function isSym(g::BinaryTreeFlag, v1::Int, v2::Int)::Bool return isSym(g.tree, v1p, v2p) end -function subFlag(F::BinaryTree, vertices::AbstractVector{Int}) +function subFlag(F::BinaryTree, vertices::AbstractVector{Int})::BinaryTree if length(vertices) == 0 return BinaryTree() end @@ -879,7 +879,7 @@ function subFlag(F::BinaryTree, vertices::AbstractVector{Int}) end -function subFlag(F::BinaryTreeFlag, vertices::AbstractVector{Int}) +function subFlag(F::BinaryTreeFlag, vertices::AbstractVector{Int})::BinaryTreeFlag # @show F # @show vertices diff --git a/src/FlagAlgebras/EdgeMarkedFlags.jl b/src/FlagAlgebras/EdgeMarkedFlags.jl index 56628e9..7e9e674 100644 --- a/src/FlagAlgebras/EdgeMarkedFlags.jl +++ b/src/FlagAlgebras/EdgeMarkedFlags.jl @@ -29,14 +29,19 @@ struct EdgeMarkedFlag{T,P} <: Flag # function EdgeMarkedFlag{T}(F::T, marked::Vector) where {T<:Flag} # return new{T,predicateType(T)}(F, marked) # end - function EdgeMarkedFlag{T}(F::T, marked::Vector{Vector{P}}) where {T<:Flag,P} - return new{T,predicateType(T)}(F, vcat(marked...)) - end - function EdgeMarkedFlag{T}(F::T, marked::Vector{Vector}) where {T<:Flag} - return new{T,predicateType(T)}(F, vcat(marked...)) - end + # function EdgeMarkedFlag{T}(F::T, marked::Vector{Vector{P}}) where {T<:Flag,P} + # allMarks::Vector{predicateType(T)} = vcat(marked...) + # return new{T,predicateType(T)}(F, allMarks) + # end + # function EdgeMarkedFlag{T}(F::T, marked::Vector{Vector}) where {T<:Flag} + # return new{T,predicateType(T)}(F, vcat(marked...)) + # end function EdgeMarkedFlag{T}(F::T, marked::Vector{P}) where {T<:Flag,P} - return new{T,predicateType(T)}(F, marked) + if P <: Vector + return new{T,predicateType(T)}(F, vcat(marked...)) + else + return new{T,predicateType(T)}(F, marked) + end end end @@ -62,7 +67,7 @@ end function findUnknownPredicates( F::EdgeMarkedFlag, fixed::Vector{U}, predLimits::Vector ) where {U<:AbstractVector{Int}} - return findUnknownPredicates(F.F, fixed, predLimits[1:end-1]), predLimits::Vector + return findUnknownPredicates(F.F, fixed, predLimits[1:(end - 1)]), predLimits::Vector end function predicateType(::Type{EdgeMarkedFlag{T,P}}) where {T<:Flag,P} @@ -83,9 +88,7 @@ function addPredicates( return EdgeMarkedFlag{T,predicateType(T)}(addPredicates(G.F, preds), G.marked) end -function permute( - F::EdgeMarkedFlag{T,P}, p::AbstractVector{Int} -) where {T<:Flag,P} +function permute(F::EdgeMarkedFlag{T,P}, p::AbstractVector{Int}) where {T<:Flag,P} return EdgeMarkedFlag{T}(permute(F.F, p), P[permute(e, p) for e in F.marked]) end @@ -107,7 +110,7 @@ function countEdges(F::EdgeMarkedFlag) return vcat(countEdges(F.F), [length(F.marked)]) end -function isolatedVertices(F::EdgeMarkedFlag) +function isolatedVertices(F::EdgeMarkedFlag)::BitVector return BitVector([false for i in 1:size(F)]) end @@ -208,7 +211,7 @@ function zeta( ) where {T<:Flag,D,P<:Predicate} # res = moebius(F; label=label) # map!(abs, values(res.coeff)) - res = sum(c*zeta(f) for (f,c) in F.coeff) + res = sum(c * zeta(f) for (f, c) in F.coeff) return res end diff --git a/src/FlagAlgebras/FreeVariables.jl b/src/FlagAlgebras/FreeVariables.jl new file mode 100644 index 0000000..29724b4 --- /dev/null +++ b/src/FlagAlgebras/FreeVariables.jl @@ -0,0 +1,115 @@ +# Sometimes one may want to add free variables to the model, for example to model sqrt(edge density). To make this possible, we model them as a "fake flag algebra", with only finite sized models and no S_infty symmetries + +export FreeVariable, FreeVariables + +struct FreeVariable <: Flag + exponent::Int + FreeVariable() = new(0) + FreeVariable(exponent::Int) = new(exponent) +end + +FreeVariables{N} = ProductFlag{NTuple{N,FreeVariable}} + +struct IncreaseDegreePredicate <: Predicate +end + +function permute(pred::IncreaseDegreePredicate, p::AbstractVector{Int}) + return pred +end + +function Base.show(io::IO, T::FreeVariable) + return print(io, "x^$(T.exponent)") +end + +function ==(A::FreeVariable, B::FreeVariable) + return A.exponent == B.exponent +end + +function hash(A::FreeVariable, h::UInt) + return hash(A.exponent, hash(:FreeVariable, h)) +end + +function distinguish(P::IncreaseDegreePredicate, v::Int, W::BitVector)::UInt + return UInt(1) +end + + +Base.one(::Type{FreeVariable})::FreeVariable = FreeVariable() + +Base.one(::FreeVariable)::FreeVariable = FreeVariable() + +Base.size(F::FreeVariable)::Int = 0 + +function labelCanonically(F::FreeVariable)::FreeVariable + return F +end + +function countEdges(F::FreeVariable)::Vector{Int} + return [F.exponent] +end + +function maxPredicateArguments(::Type{FreeVariable}) + return 1 +end + +function subFlag(F::FreeVariable, vertices::Vector{Int})::FreeVariable + return F +end + +function glue(F::FreeVariable, G::FreeVariable, p::AbstractVector{Int}) + return FreeVariable(F.exponent + G.exponent) +end + +function distinguish(F::FreeVariable, v::Int, W::BitVector)::UInt + return UInt(1) +end + +function isolatedVertices(F::FreeVariable)::BitVector + return [true;] +end + +function predicateType(::Type{FreeVariable}) + return IncreaseDegreePredicate +end + +function addPredicates(F::FreeVariable, preds::Vector{IncreaseDegreePredicate}) + return FreeVariable(F.exponent + length(preds)) +end + +function permute(F::FreeVariable, p::AbstractVector{Int}) + return F +end + +function findUnknownPredicates( + F::FreeVariable, fixed::Vector{U}, predLimits::Vector +) where {U<:AbstractVector{Int}} + if any(1 in f for f in fixed) + return [IncreaseDegreePredicate[]] + end + return [IncreaseDegreePredicate[IncreaseDegreePredicate()]] +end + +function findUnknownGenerationPredicates( + F::FreeVariable, fixed::Vector{U}, predLimits::Vector +) where {U<:AbstractVector{Int}} + return findUnknownPredicates(F, fixed, predLimits) +end + +function isSym(F::FreeVariable, v1::Int, v2::Int)::Bool + return v1 == v2 +end + +function allowMultiEdges(::Type{FreeVariable}) + return true +end + +# function generateAll( +# ::Type{FreeVariable{F}}, maxVertices::Int, maxPredicates::Vector{Int} +# ) where {F<:Flag} +# tmp = generateAll(F, maxVertices, maxPredicates) +# return [FreeVariable{F}(f) for f in tmp] +# end + +# function eliminateIsolated(F::FreeVariable) +# return FreeVariable(filter(x -> x.second != 0, F.exponent)) +# end \ No newline at end of file diff --git a/src/FlagAlgebras/Graphs.jl b/src/FlagAlgebras/Graphs.jl index b898027..1849dc7 100644 --- a/src/FlagAlgebras/Graphs.jl +++ b/src/FlagAlgebras/Graphs.jl @@ -14,6 +14,25 @@ struct Graph <: Flag Graph() = new(Symmetric(BitMatrix(undef, 0, 0))) end +function Base.show(io::IO, G::Graph) + if size(G.A,1) == 0 + print(io, "∅ ") + else + print(io, "($(size(G.A,1)),") + first = true + for c in eachindex(G.A) + if c[1] <= c[2] && G.A[c] + if !first + print(io, " ") + end + print(io, "$(c[1])-$(c[2])") + first = false + end + end + print(io, ")") + end +end + import Base.== function ==(A::Graph, B::Graph) return A.A == B.A @@ -53,19 +72,25 @@ end function findUnknownPredicates( F::Graph, fixed::Vector{U}, predLimits::Vector )::Vector{Vector{EdgePredicate}} where {U<:AbstractVector{Int}} - @assert length(predLimits) in [0,1] + @assert length(predLimits) in [0, 1] res = EdgePredicate[] if length(predLimits) == 1 && countEdges(F)[1] >= predLimits[1] - return [res] + return [res] end for e in combinations(1:size(F), 2) - if (!F.A[e...]) && !any(issubset(e, f) for f in fixed) + if (!F.A[e[1], e[2]]) && !any(issubset(e, f) for f in fixed) push!(res, EdgePredicate(e...)) end end return [res] end +# function findUnknownGenerationPredicates( +# F::Graph, fixed::Vector{U}=Vector{Int}[], predLimits::Vector=Int[] +# ) where {U<:AbstractVector{Int}} +# return findUnknownPredicates(F, fixed, predLimits) +# end + function addPredicates(G::Graph, preds::Vector{EdgePredicate}) A = Matrix(copy(G.A)) for p in preds @@ -190,7 +215,7 @@ function isBipartite(F::Graph) end function distinguish(F::Graph, v::Int, W::BitVector)::UInt - @views return hash(sum(F.A[W, v])) + @inbounds return hash(sum(F.A[W, v])) end function distinguish(F::EdgePredicate, v::Int, W::BitVector)::UInt @@ -202,5 +227,10 @@ function countEdges(F::Graph)::Vector{Int} end function isolatedVertices(F::Graph)::BitVector - return vec(.!any(F.A; dims=1)) + # res = BitVector(zeros(Bool, size(F))) + res = BitArray(undef, size(F)) + any!(res, F.A) + map!(!, res, res) + # return .!vec(any(F.A; dims=1)) + return res end diff --git a/src/FlagAlgebras/ProductFlag.jl b/src/FlagAlgebras/ProductFlag.jl index 0a67e14..90cf659 100644 --- a/src/FlagAlgebras/ProductFlag.jl +++ b/src/FlagAlgebras/ProductFlag.jl @@ -13,6 +13,19 @@ struct ProductFlag{FT} <: Flag where {FT<:Tuple{Vararg{Flag}}} ProductFlag{FT}() where {FT} = new{FT}(Tuple(F() for F in fieldtypes(FT))) end +function Base.show(io::IO, F::ProductFlag) + print(io, "(") + first = true + for f in F.Fs + if !first + print(io, ",") + end + print(io, f) + first = false + end + print(io, ")") +end + import Base.== function ==(A::ProductFlag{FT}, B::ProductFlag{FT}) where {FT} return all(A.Fs[i] == B.Fs[i] for i = 1:length(A.Fs)) @@ -60,7 +73,7 @@ function findUnknownPredicates( if length(predLimits) == length(fieldtypes(FT)) FIP = findUnknownPredicates(F.Fs[i], fixed, predLimits[i]) - else + else FIP = findUnknownPredicates(F.Fs[i], fixed, [predLimits[1]]) end @@ -87,7 +100,7 @@ function findUnknownGenerationPredicates( # end if length(predLimits) == length(fieldtypes(FT)) FIP = findUnknownGenerationPredicates(F.Fs[i], fixed, predLimits[i]) - else + else FIP = findUnknownGenerationPredicates(F.Fs[i], fixed, [predLimits[1]]) end # FIP = findUnknownGenerationPredicates(F.Fs[i], fixed, predLimits[i]) @@ -106,7 +119,7 @@ function addPredicates( G::ProductFlag{FT}, preds::Vector{Tuple{Int,P}}) where {FT,P} tmp = [] for (i, F) in enumerate(G.Fs) - newF = addPredicates(F, [p[2] for p in preds if p[1] == i]) + newF = addPredicates(F, predicateType(typeof(F))[p[2] for p in preds if p[1] == i]) push!(tmp, newF) end @@ -156,7 +169,11 @@ end function permute( F::ProductFlag{FT}, p::AbstractVector{Int} ) where {FT} - return ProductFlag{FT}(Tuple(permute(f, p) for f in F.Fs)) + pF = [permute(f, p) for f in F.Fs] + if any(x -> x === nothing, pF) + return nothing + end + return ProductFlag{FT}(Tuple(pF)) end function countEdges(F::ProductFlag) @@ -180,13 +197,13 @@ end function vertexColor(F::ProductFlag{FT}, v::Int) where {FT} # colorCombinations = sort!(unique([[vertexColor(f, i) for f in F.Fs] for i = 1:size(F)])) - + cs = hash(Int[vertexColor(f, v) for f in F.Fs]) colorCombinations = UInt[cs] # colorCombinations = Vector{Int}[cs] - + for i = 1:size(F) - if i != v + if i != v ct = hash(Int[vertexColor(f, i) for f in F.Fs]) if !(ct in colorCombinations) push!(colorCombinations, ct) @@ -196,7 +213,7 @@ function vertexColor(F::ProductFlag{FT}, v::Int) where {FT} sort!(colorCombinations) - return findfirst(x->x==cs, colorCombinations) + return findfirst(x -> x == cs, colorCombinations) # @assert length(unique(cs)) == 1 # return cs[1] diff --git a/src/FlagAlgebras/QuantumFlags.jl b/src/FlagAlgebras/QuantumFlags.jl index 19b8e4b..eda5692 100644 --- a/src/FlagAlgebras/QuantumFlags.jl +++ b/src/FlagAlgebras/QuantumFlags.jl @@ -11,6 +11,7 @@ mutable struct QuantumFlag{F<:Flag,T} QuantumFlag{F,T}(cs::Dict{F,T}) where {F<:Flag,T} = new(cs) QuantumFlag{F,T}(opts...) where {F<:Flag,T} = new(Dict{F,T}(opts...)) QuantumFlag{F}(fc::QuantumFlag{F,T}) where {F<:Flag,T} = new{F,T}(fc.coeff) + QuantumFlag{F,T}(fc::QuantumFlag{F,U}) where {F<:Flag,T, U} = new{F,T}(Dict{F,T}(fc.coeff)) end function Base.promote_rule( diff --git a/src/FlagAlgebras/SymmetricFunctions.jl b/src/FlagAlgebras/SymmetricFunctions.jl index 4cdd679..833db7b 100644 --- a/src/FlagAlgebras/SymmetricFunctions.jl +++ b/src/FlagAlgebras/SymmetricFunctions.jl @@ -108,11 +108,11 @@ function findUnknownPredicates( return [[LabelPredicate(i) for i in 1:size(F) if !any((i in part) for part in fixed)]] end -function findUnknownGenerationPredicates( - F::SymmetricFunction, fixed::Vector{U}, predLimits::Vector -) where {U<:AbstractVector{Int}} - return findUnknownPredicates(F, fixed, predLimits) -end +# function findUnknownGenerationPredicates( +# F::SymmetricFunction, fixed::Vector{U}, predLimits::Vector +# ) where {U<:AbstractVector{Int}} +# return findUnknownPredicates(F, fixed, predLimits) +# end function isSym(F::SymmetricFunction, v1::Int, v2::Int)::Bool return get(F.exponents, v1, 0) == get(F.exponents, v2, 0) diff --git a/src/FlagModels/FlagModel.jl b/src/FlagModels/FlagModel.jl index 3145efd..8dc43e1 100644 --- a/src/FlagModels/FlagModel.jl +++ b/src/FlagModels/FlagModel.jl @@ -87,7 +87,7 @@ function addLasserreBlock!( push!(m.subModels, lM) Fs = generateAll(T, Int(floor(maxVertices / 2)), Int.(floor.(maxEdges / 2))) - # @show Fs + display(Fs) for F in Fs if isAllowed(m, F) addFlag!(lM, F) @@ -370,7 +370,7 @@ function verifySOS( Base.show(io, m) println() - tmp = sum(verifySOS(m.subModels[i], sol[i]; io=io) for i in 1:length(m.subModels)) + tmp = labelCanonically(sum(verifySOS(m.subModels[i], sol[i]; io=io) for i in 1:length(m.subModels))) println(io, "SOS result:") println(io, "$(tmp) >= 0") @@ -400,6 +400,9 @@ function verifySOS( println(io, "Error:") println(io, "$(err)") + if haskey(m.objective.coeff, T()) + constTerm -= m.objective.coeff[T()] + end res = constTerm + err + m.objective println(io, "Final rigorous bound:") println(io, "$(res) >= 0") diff --git a/src/FlagModels/LasserreModel.jl b/src/FlagModels/LasserreModel.jl index 7eb4cf8..8ddaf03 100644 --- a/src/FlagModels/LasserreModel.jl +++ b/src/FlagModels/LasserreModel.jl @@ -62,6 +62,16 @@ struct FlagSymmetries{T<:Flag} end end +import Base.== +function ==(A::FlagSymmetries, B::FlagSymmetries) + return A.F == B.F +end +import Base.hash +function hash(A::FlagSymmetries, h::UInt) + return hash(A.F, hash(:FlagSymmetries, h)) +end + + struct SpechtFlag{T<:Flag} F::FlagSymmetries{T} T::AbstractAlgebra.Generic.YoungTableau{Int} @@ -98,6 +108,13 @@ mutable struct LasserreModel{T<:Flag,N,D} <: AbstractFlagModel{T,N,D} end end +function Base.show(io::IO, m::LasserreModel{T,N,D}) where {T,N,D} + return println( + io, + "Lasserre style blocks based on $(length(m.generators)) different flags with $(length(m.basis)) blocks", + ) +end + function modelSize(m::LasserreModel) return Partition([length(b) for b in values(m.basis)]) end @@ -500,7 +517,7 @@ end function computeSDP!(m::LasserreModel{T,N,D}, reservedVerts::Int) where {T,N,D}#; maxVert = -1, useGroups = true, splitByOverlaps = false) #TODO: Parallelize better - totalNum = Int64(sum(c * (c + 1) / 2 for c in length.(values(m.basis)); init = 0)) + totalNum = Int64(sum(c * (c + 1) / 2 for c in length.(values(m.basis)); init=0)) t = 1 @@ -650,3 +667,105 @@ function buildJuMPModel( return (model=jumpModel, variables=graphCoefficients, blocks=Y, constraints=constraints) end + +function roundResults(m::LasserreModel{T,N,D}, jumpModel, variables, blocks, constraints; prec=1e-5) where {T,N,D} + ex = Dict() + + den = round(BigInt, 1 / prec) + function roundDen(x) + return round(BigInt, den * x) // den + end + + for (mu, b) in blocks + if mu isa String + # ex[mu] = rationalize(BigInt, value(b); tol=prec)#; digits = digits) + ex[mu] = roundDen(value(b)) + else + ex[mu] = roundRationalPSD(value(b); baseType=BigInt, prec=prec) + end + end + + return ex +end + +function verifySOS(m::LasserreModel, sol::Dict; io::IO=stdout) + if io !== nothing + println(io, "Lasserre model") + + for mu in keys(sol) + if mu isa String + continue + end + + println(io, "\nBlock of type $mu with flags:") + println.(io, m.basis[mu]) + + println(io, "PSD matrix:") + psd = sol[mu] isa Matrix ? sol[mu] : sol[mu].psd + show(io, "text/plain", psd) + println(io) + if !(sol[mu] isa Matrix) && size(sol[mu].psd, 1) > 1 + println(io, "With Cholesky factorization:") + show(io, "text/plain", sol[mu].chol) + println(io) + end + + println(io, "Nonzero flag products:") + n = size(m.basis[mu], 1) + # @show mu + # @show sol[mu] + psd = sol[mu] isa Matrix ? sol[mu] : sol[mu].psd + for i in 1:n + for j in i:n + if !iszero(psd[i, j]) + tmp = sum(m.sdpData) do (G, B) + if haskey(B, mu) + return Rational{BigInt}(Symmetric(B[mu])[i, j]) * G + else + return BigInt(0) // 1 * G + end + end + print( + io, + "Product at $((i,j)) is $(psd[i,j]) ⋅ $(m.basis[mu][i]) ⋅ $(m.basis[mu][j]) ", + ) + # println(io, "$(psd[i,j]) ⋅ ($tmp)") + println(io, "=$(Rational{BigInt}(psd[i,j])*tmp)") + end + end + end + println(io, "Block $mu total:") + println( + io, + "$(sum(m.sdpData) do (G, B) + if haskey(B, mu) + return dot(B[mu], psd) * G + else + return BigInt(0) // 1 * G + end + end)>=0", + ) + end + + + end + + res = sum(m.sdpData) do (G, B) + sum(B) do (mu, b) + if haskey(sol, mu) + psd = sol[mu] isa Matrix ? sol[mu] : sol[mu].psd + return dot(Rational{BigInt}.(psd), Symmetric(b)) * G + else + return BigInt(0) // 1 * G + end + end + end + + if io !== nothing + println(io, "\nLasserre model result:") + println(io, "$(res)>= 0") + println(io) + end + + return res +end diff --git a/src/FlagModels/QuadraticModule.jl b/src/FlagModels/QuadraticModule.jl index 9eb32e2..2d7c3c3 100644 --- a/src/FlagModels/QuadraticModule.jl +++ b/src/FlagModels/QuadraticModule.jl @@ -33,6 +33,17 @@ mutable struct QuadraticModule{T<:Flag,U<:Flag,B,N,D} <: end end +function Base.show(io::IO, m::QuadraticModule{T,U,B,N,D}) where {T,U,B,N,D} + println(io, "Quadratic module for inequality $(m.inequality)>=0\nWith base model:") + return println(io, m.baseModel) +end + +function roundResults( + m::QuadraticModule{T,U,B,N,D}, jumpModel, variables, blocks, constraints; prec=1e-5 +) where {T,U,N,D,B} + return roundResults(m.baseModel, jumpModel, variables, blocks, constraints) +end + function computeSDP!( m::QuadraticModule{T,U,B,N,D}, reservedVerts::Int ) where {T<:Flag,U<:Flag,N,D,B<:AbstractFlagModel{U,N,D}} @@ -70,7 +81,9 @@ function modelBlockSizes(m::QuadraticModule) return modelBlockSizes(m.baseModel) end -function buildJuMPModel(m::QuadraticModule{T,U,B,N,D}, replaceBlocks=Dict(), jumpModel=Model()) where {T,U,B,N,D} +function buildJuMPModel( + m::QuadraticModule{T,U,B,N,D}, replaceBlocks=Dict(), jumpModel=Model() +) where {T,U,B,N,D} b = modelBlockSizes(m) Y = Dict() constraints = Dict() @@ -157,12 +170,14 @@ mutable struct EqualityModule{T<:Flag,U<:Flag,N,D} <: AbstractFlagModel{T,N,D} end end -function roundResults(m::EqualityModule{T,U,N,D}, jumpModel, variables, blocks, constraints; prec=1e-5) where {T,U,N,D} +function roundResults( + m::EqualityModule{T,U,N,D}, jumpModel, variables, blocks, constraints; prec=1e-5 +) where {T,U,N,D} ex = Dict() - den = round(BigInt, 1/prec) + den = round(BigInt, 1 / prec) function roundDen(x) - return round(BigInt, den*x)//den + return round(BigInt, den * x)//den end for (mu, b) in blocks @@ -173,10 +188,11 @@ function roundResults(m::EqualityModule{T,U,N,D}, jumpModel, variables, blocks, return ex end - - function Base.show(io::IO, m::EqualityModule{T,U,N,D}) where {T,U,N,D} - println(io, "Equality constraint ($(m.equality) == 0) multiplied with $(length(m.basis)) flags.") + return println( + io, + "Equality constraint ($(m.equality) == 0) multiplied with $(length(m.basis)) flags.", + ) end function computeSDP!( @@ -217,7 +233,9 @@ function modelBlockSizes(m::EqualityModule) return Dict(i => -1 for i in 1:length(m.basis)) end -function buildJuMPModel(m::EqualityModule{T,U,N,D}, replaceBlocks=Dict(), jumpModel=Model()) where {T,U,N,D} +function buildJuMPModel( + m::EqualityModule{T,U,N,D}, replaceBlocks=Dict(), jumpModel=Model() +) where {T,U,N,D} @assert length(replaceBlocks) == 0 b = modelBlockSizes(m) @@ -266,25 +284,46 @@ function verifySOS(m::EqualityModule, sol::Dict; io::IO=stdout) println(io, "Equality module coming from constraint") println(io, "$(m.equality)= 0") for i in keys(sol) - if sol[i] != 0 // 1 + if sol[i] != 0//1 println(io, "Times $(sol[i])$(m.basis[i]) :") - print(io, sum(m.sdpData) do (G, B) - sum(B) do (j, c) - j != i && return 0 * one(G) - get(sol, j, 0 // 1) * c * G - end - end + print( + io, + sum(m.sdpData) do (G, B) + sum(B) do (j, c) + j != i && return 0 * one(G) + get(sol, j, 0//1) * c * G + end + end, ) println(io, "= 0") end end res = sum(m.sdpData) do (G, B) sum(B) do (i, c) - get(sol, i, 0 // 1) * c * G + get(sol, i, 0//1) * c * G end end println(io, "Equality module result:") println(io, "$(res)= 0") println(io) - res + return res +end + +function verifySOS( + m::QuadraticModule{T,U,B,N,D}, sol::Dict; io::IO=stdout +) where {T,U,B,N,D} + println(io, "Quadratic module coming from constraint") + println(io, "$(m.inequality)>= 0") + println(io, "With base model:") + res = verifySOS(m.baseModel, sol; io=io) + + res = labelCanonically(typeof(res)(m.inequality) * res) + + if io !== nothing + println(io, "\nQuadratic module result:") + println(io, "$(res)>= 0") + println(io) + end + + return res end diff --git a/src/FlagModels/RazborovModel.jl b/src/FlagModels/RazborovModel.jl index db8f887..895c552 100644 --- a/src/FlagModels/RazborovModel.jl +++ b/src/FlagModels/RazborovModel.jl @@ -22,7 +22,7 @@ mutable struct RazborovModel{T<:Flag,N,D} <: AbstractFlagModel{T,N,D} Dict(), Dict{T,Dict{T,SparseMatrixCSC{D,Int}}}(), parent, - Dict() + Dict(), ) end @@ -32,13 +32,16 @@ mutable struct RazborovModel{T<:Flag,N,D} <: AbstractFlagModel{T,N,D} Dict(), Dict{T,Dict{T,SparseMatrixCSC{Int,Int}}}(), parent, - Dict() + Dict(), ) end end function Base.show(io::IO, m::RazborovModel{T,N,D}) where {T,N,D} - println(io, "Flagmatic style blocks with $(length(m.basis)) types and $(sum(length.(values(m.basis)))) flags") + return println( + io, + "Flagmatic style blocks with $(length(m.basis)) types and $(sum(length.(values(m.basis)))) flags", + ) end function isAllowed(m::RazborovModel{T,N,D}, F::T) where {T<:Flag,N,D} @@ -76,7 +79,7 @@ function computeUnreducedRazborovBasis( FBlock = label(F; removeIsolated=false)[1] @assert size(FBlock) == m # @assert FBlock == label(FBlock; removeIsolated=false)[1] - FExtended = permute(FBlock, 1:(m+k)) # add isolated vertices in unlabeled part + FExtended = permute(FBlock, 1:(m + k)) # add isolated vertices in unlabeled part preds = vcat(findUnknownPredicates(FExtended, [1:m])...) @@ -89,16 +92,17 @@ function computeUnreducedRazborovBasis( return razborovBasis end -function computeRazborovBasis!(M::RazborovModel{T,N,D}, n; maxLabels=n, maxBlockSize = Inf) where {T<:Flag,N,D} +function computeRazborovBasis!( + M::RazborovModel{T,N,D}, n; maxLabels=n, maxBlockSize=Inf +) where {T<:Flag,N,D} razborovBasis = computeUnreducedRazborovBasis(M, n, maxLabels) # display(razborovBasis) - reducedBasis = Dict(mu => unique(labelCanonically.(B)) for (mu, B) in razborovBasis) - if maxBlockSize < Inf - filter!(x->length(x[2]) < maxBlockSize, reducedBasis) + if maxBlockSize < Inf + filter!(x -> length(x[2]) < maxBlockSize, reducedBasis) end - + @info "basis reduced" @info "determining symmetries" for (mu, B) in reducedBasis @@ -117,7 +121,7 @@ function computeRazborovBasis!(M::RazborovModel{T,N,D}, n; maxLabels=n, maxBlock @assert length(p) == b.n pb = labelCanonically( PartiallyLabeledFlag{T}( - permute(b.F, vcat(p, (length(p)+1):size(b))), b.n + permute(b.F, vcat(p, (length(p) + 1):size(b))), b.n ), ) gen[i] = findfirst(x -> x == pb, B) @@ -226,16 +230,14 @@ function computeSDP!(m::RazborovModel{T,N,D}, reservedVerts::Int) where {T,N,D} n = length(B) M = zeros(Int, n, n) ind = 1 - for i = 1:n - for j = i:n + for i in 1:n + for j in i:n M[i, j] = ind M[j, i] = ind ind += 1 end end - P = ( - pattern=M, n=n, fullPattern=M - ) + P = (pattern=M, n=n, fullPattern=M) else P = m.blockSymmetry[mu] end @@ -258,7 +260,7 @@ function computeSDP!(m::RazborovModel{T,N,D}, reservedVerts::Int) where {T,N,D} newSize = k + (n1 - k) + (n2 - k) p = collect(1:newSize) - p[(k+1):n1] = (n2+1):newSize + p[(k + 1):n1] = (n2 + 1):newSize T1 = a.F p1 = p[1:size(a.F)] @@ -272,7 +274,7 @@ function computeSDP!(m::RazborovModel{T,N,D}, reservedVerts::Int) where {T,N,D} p1Fin = p2Inv[p1] p1Fin = vcat(p1Fin, setdiff(1:newSize, p1Fin)) - @views sort!(p1Fin[(n1+1):end]) + @views sort!(p1Fin[(n1 + 1):end]) p1Fin = Int.(p1Fin) @@ -338,7 +340,7 @@ function computeSDP!(m::RazborovModel{T,N,D}, reservedVerts::Int) where {T,N,D} # m.sdpData[F][mu] .+= (norm(A) / norm(P.reg[s])) * d*P.reg[s]#*factor m.sdpData[F][mu] .+= d * P.reg[s]#*factor else - m.sdpData[F][mu][P.pattern.==s] .= d + m.sdpData[F][mu][P.pattern .== s] .= d end end end @@ -506,7 +508,6 @@ function buildJuMPModel(m::RazborovModel, replaceBlocks=Dict(), jumpModel=Model( end for (i, F) in enumerate(m.quotient) - y = @variable(jumpModel, base_name = "quotient$i") Y["Q$i"] = y for (G, c) in F.coeff @@ -606,7 +607,7 @@ function computeUnreducedRazborovBasis( q = collect(1:m) q[c] .= 1:k - q[setdiff(1:m, c)] .= (k+1):m + q[setdiff(1:m, c)] .= (k + 1):m # @show k q[c] .= p.d[tCanLabelPermInv[1:k]] @@ -640,12 +641,14 @@ function computeUnreducedRazborovBasis( return razborovBasis end -function roundResults(m::RazborovModel{T,N,D}, jumpModel, variables, blocks, constraints; prec=1e-5) where {T,N,D} +function roundResults( + m::RazborovModel{T,N,D}, jumpModel, variables, blocks, constraints; prec=1e-5 +) where {T,N,D} ex = Dict() den = round(BigInt, 1 / prec) function roundDen(x) - return round(BigInt, den * x) // den + return round(BigInt, den * x)//den end for (mu, b) in blocks @@ -687,17 +690,20 @@ function verifySOS(m::RazborovModel, sol::Dict; io::IO=stdout) # @show mu # @show sol[mu] psd = sol[mu] isa Matrix ? sol[mu] : sol[mu].psd - for i = 1:n - for j = i:n + for i in 1:n + for j in i:n if !iszero(psd[i, j]) tmp = sum(m.sdpData) do (G, B) if haskey(B, mu) return Rational{BigInt}(B[mu][i, j]) * G else - return BigInt(0) // 1 * G + return BigInt(0)//1 * G end end - print(io, "Product at $((i,j)) is $(psd[i,j]) ⋅ $(m.basis[mu][i]) ⋅ $(m.basis[mu][j]) ") + print( + io, + "Product at $((i,j)) is $(psd[i,j]) ⋅ $(m.basis[mu][i]) ⋅ $(m.basis[mu][j]) ", + ) # println(io, "$(psd[i,j]) ⋅ ($tmp)") println(io, "=$(Rational{BigInt}(psd[i,j])*tmp)") end @@ -712,7 +718,7 @@ function verifySOS(m::RazborovModel, sol::Dict; io::IO=stdout) else return BigInt(0) // 1 * G end - end)>=0" + end)>=0", ) end @@ -728,32 +734,30 @@ function verifySOS(m::RazborovModel, sol::Dict; io::IO=stdout) io, "$(sum(enumerate(m.quotient)) do (i, F) Rational{BigInt}(get(sol, "Q$i", 0 // 1)) * F -end)=0" +end)=0", ) end - res = sum(m.sdpData) do (G, B) sum(B) do (mu, b) if haskey(sol, mu) psd = sol[mu] isa Matrix ? sol[mu] : sol[mu].psd return dot(Rational{BigInt}.(psd), b) * G else - return BigInt(0) // 1 * G + return BigInt(0)//1 * G end end end res += sum(enumerate(m.quotient)) do (i, F) - Rational{BigInt}(get(sol, "Q$i", 0 // 1)) * F + Rational{BigInt}(get(sol, "Q$i", 0//1)) * F end if io !== nothing - println(io, "\nRazborov model result:") println(io, "$(res)>= 0") println(io) end - res + return res end diff --git a/src/utils/Nauty.jl b/src/utils/Nauty.jl index 35ab23b..f49ee3f 100644 --- a/src/utils/Nauty.jl +++ b/src/utils/Nauty.jl @@ -4,286 +4,360 @@ include("SchreierSims.jl") using DataStructures -# Custom implementation of nauty, but still sloooow -# Allows us to label near-arbitrary combinatoric objects -function label(F::T; prune=true, removeIsolated=true) where {T} - # @show F - # display(F) - # @show isolatedVertices(F) - if removeIsolated - isoV = isolatedVertices(F) - if any(isoV) - F = subFlag(F, .!isoV) - end - end - - n = size(F) +# function refine!(F::T, coloring::Vector{Int}, alpha)::Vector{UInt} where {T<:Flag} +# # Q .= 0 - if n == 0 - return F, Group(), [] - end +# end - nInv1 = UInt[] - nInvStar = UInt[] - p1 = Int[] - pStar = Int[] +# function individualize!(coloring::Vector{Int}, v::Int)::Vector{Int} +# return coloring +# end - autG = Group() - - first = true - - v1 = zeros(Int, n) - vStar = zeros(Int, n) - - curBranch = Int[] - covered = Vector{Int}[] - - # To avoid allocations in refine! +function refine!(coloring::Vector{Int}, F, v::Int)::Vector{UInt} + n = size(F) alpha = BitVector(undef, n) - # newCellSizes = zeros(Int, n) - vertDistinguish = zeros(UInt, n) + alpha .= false + if v == 0 + alpha[1] = true + else + alpha[coloring[v]] = true + # individualize!(coloring, v) + for i in 1:n + if coloring[i] >= coloring[v] && i != v + coloring[i] += 1 + end + end + end - # Q = UInt[0 for i in 1:n, j in 1:n] + # nodeInv = refine!(coloring) Q = zeros(UInt, n, n) - WCellInd = BitVector(undef, n) - curCell = BitVector(undef, n) - newCells = Dict{UInt,Int}() - - function refine!(coloring::Vector{Int}) - Q .= 0 - @inbounds while any(alpha) && maximum(coloring) < n - W = findfirst(alpha) - WCellInd .= coloring .== W - alpha[W] = false - - cell = 1 - while cell <= maximum(coloring) - vertDistinguish .= 0 - curCell .= coloring .== cell - - empty!(newCells) - # newCells = Dict{UInt,Int}() - - for v in findall(curCell)#findall(x->x==cell, coloring) - d = distinguish(F, v, WCellInd) - vertDistinguish[v] = d - # union!(newCells, [d]) - newCells[d] = get(newCells, d, 0) + 1 - end + curCell::BitVector = BitVector(undef, n) + @inbounds while maximum(coloring) < n + W::Int = 0 + for i in 1:length(alpha) + if alpha[i] + W = i + break + end + end + if W == 0 + break + end + # W::Int = findfirst(alpha) + # WCellInd .= coloring .== W + WCellInd::BitVector = coloring .== W + alpha[W] = false + + cell::Int = 1 + newCells = Dict{UInt,Int}() + while cell <= maximum(coloring) + vertDistinguish .= 0 + # curCell .= coloring .== cell + for i = 1:n + curCell[i] = coloring[i] == cell + end - # To make it invariant of initial vertex labeling, sort the hashes - # newCells = Dict(d=>i+cell-1 for (i,d) in enumerate(sort!([u for u in unique(vertDistinguish) if u != 0]))) + empty!(newCells) - # newCellSizes = Dict(d=>count(x->x==d, vertDistinguish) for d in keys(newCells)) + for v in findall(curCell)#findall(x->x==cell, coloring) + d::UInt = distinguish(F, v, WCellInd) + vertDistinguish[v] = d + # union!(newCells, [d]) + newCells[d] = get(newCells, d, 0) + 1 + end - numNewCells = length(newCells) - 1 + # To make it invariant of initial vertex labeling, sort the hashes + # newCells = Dict(d=>i+cell-1 for (i,d) in enumerate(sort!([u for u in unique(vertDistinguish) if u != 0]))) - for i in eachindex(coloring) - if coloring[i] > cell - coloring[i] += numNewCells - end - end + # newCellSizes = Dict(d=>count(x->x==d, vertDistinguish) for d in keys(newCells)) - # coloring[coloring .> cell] .+= numNewCells - for i in length(alpha):-1:(cell+numNewCells+1) - alpha[i] = alpha[i-numNewCells] + numNewCells = length(newCells) - 1 + + for i in eachindex(coloring) + if coloring[i] > cell + coloring[i] += numNewCells end - # @views alpha[(cell + numNewCells + 1):end] .= alpha[(cell + 1):(end - numNewCells)] + end - alpha[(cell+1):(cell+numNewCells)] .= true + # coloring[coloring .> cell] .+= numNewCells + for i in length(alpha):-1:(cell + numNewCells + 1) + alpha[i] = alpha[i - numNewCells] + end + # @views alpha[(cell + numNewCells + 1):end] .= alpha[(cell + 1):(end - numNewCells)] - newCellsCol = collect(newCells) - sort!(newCellsCol; by=x -> x[1]) + alpha[(cell + 1):(cell + numNewCells)] .= true - # sorting newCellsCol by hashes breaks things?!? But I need to sort! + newCellsCol::Vector{Pair{UInt,Int}} = collect(newCells) + sort!(newCellsCol; by=x -> x[1]) - for v in findall(curCell)#findall(x->x==cell, coloring) - # coloring[v] = newCells[vertDistinguish[v]] + # sorting newCellsCol by hashes breaks things?!? But I need to sort! - d = vertDistinguish[v] - coloring[v] = findfirst(x -> x[1] == d, newCellsCol) + cell - 1 - Q[coloring[v], W > cell ? W + numNewCells : W] = d + for v in 1:n + !curCell[v] && continue + # findall(curCell)#findall(x->x==cell, coloring) + # coloring[v] = newCells[vertDistinguish[v]] - # Q = hash(Q, vertDistinguish[v]) + d::UInt = vertDistinguish[v] + firstPos::Int = 0 + for i in 1:length(newCellsCol) + if newCellsCol[i][1] == d + firstPos = i + break + end end + @assert firstPos > 0 - if !alpha[cell] - alpha[cell] = true + # findfirst(x -> x[1] == d, newCellsCol) + coloring[v] = firstPos + cell - 1 + Q[coloring[v], W > cell ? W + numNewCells : W] = d - maxCellSize = 0 - maxCellInd = 0 - for (i, (d, s)) in enumerate(newCellsCol) - if s > maxCellSize - maxCellSize = s - maxCellInd = i - end - end + # Q = hash(Q, vertDistinguish[v]) + end - # Remove biggest new cell - alpha[maxCellInd+cell-1] = false + if !alpha[cell] + alpha[cell] = true + + maxCellSize = 0 + maxCellInd = 0 + for (i, (d, s)) in enumerate(newCellsCol) + if s > maxCellSize + maxCellSize = s + maxCellInd = i + end end - cell += length(newCellsCol) + # Remove biggest new cell + alpha[maxCellInd + cell - 1] = false end + + cell += length(newCellsCol) end + end - cl = maximum(coloring) - nodeInvariant = hash(Q, hash([count(x -> x == i, coloring) for i in 1:cl])) + cl = maximum(coloring) + nodeInvariant = hash(Q, hash([count(x -> x == i, coloring) for i in 1:cl])) - if maximum(coloring) == n - # append permuted Flag - # return UInt[nodeInvariant, hash(permute(F, coloring))] - # return nodeInvariant, hash(permute(F, coloring)) + if maximum(coloring) == n + # append permuted Flag + # return UInt[nodeInvariant, hash(permute(F, coloring))] + # return nodeInvariant, hash(permute(F, coloring)) - return [nodeInvariant, hash(permute(F, coloring))] - # return hash(nodeInvariant, hash(permute(F, coloring))) - end - return nodeInvariant + return UInt[nodeInvariant, hash(permute(F, coloring))] + # return hash(nodeInvariant, hash(permute(F, coloring))) end + return UInt[nodeInvariant] + # return nodeInv +end - function individualize!(coloring::Vector{Int}, v::Int) - for i in 1:n - if coloring[i] >= coloring[v] && i != v - coloring[i] += 1 - end - end - return coloring - end +function investigateNode(F,coloring::Vector{Int}, nodeInv::Vector{UInt},n,nInv1, nInvStar,autG,v1,vStar,curBranch, covered,prune)::Int + # @info "investigate node at $coloring, $nodeInv" - function refine!(coloring::Vector{Int}, v::Int) - alpha .= false - alpha[coloring[v]] = true - individualize!(coloring, v) - nodeInv = refine!(coloring) - return nodeInv - end + first = length(nInv1) == 0 - function investigateNode(coloring::Vector{Int}, nodeInv=UInt[]) - # @info "investigate node at $coloring, $nodeInv" - if maximum(coloring) == n - - # check for automorphisms - if !first - p2 = coloring - autom1 = v1[p2] - changed = false - if F == permute(F, autom1) - changed = addGen!(autG, autom1) - end - autom2 = vStar[p2] - if F == permute(F, autom2) - changed |= addGen!(autG, autom2) - end + if maximum(coloring) == n - # if changed, check for automorphism pruning higher up the tree - if changed - # @info "Added onto stabilizer, new order is $(order(autG))" - stabilizer!(autG, curBranch) - for i in 1:length(curBranch) - H = stabilizer!(autG, curBranch[1:(i-1)]) - - o = covered[i] - orbit!(H, o) - if prune && curBranch[i] in o - # cut this branch multiple levels up - difLevel = length(curBranch) - i - curBranch = curBranch[1:i] - covered = covered[1:i] - # @info "Pruning $difLevel levels from level $(length(curBranch)) via automorphism." - return difLevel - end - end - end + # check for automorphisms + if !first + p2 = coloring + autom1 = v1[p2] + changed = false + if F == permute(F, autom1) + changed = addGen!(autG, autom1) end - - # improve v1 and vStar - # if first #|| nodeInv < nInv1 - if first || nodeInv < nInv1 - for i in 1:n - v1[coloring[i]] = i - end - nInv1 = nodeInv + autom2 = vStar[p2] + if F == permute(F, autom2) + changed |= addGen!(autG, autom2) end - if first || nodeInv > nInvStar - for i in 1:n - vStar[coloring[i]] = i + + # if changed, check for automorphism pruning higher up the tree + if changed + # @info "Added onto stabilizer, new order is $(order(autG))" + stabilizer!(autG, curBranch) + for i in 1:length(curBranch) + H = stabilizer!(autG, curBranch[1:(i - 1)]) + + @assert H !== nothing + + o = covered[i] + orbit!(H, o) + if prune && curBranch[i] in o + # cut this branch multiple levels up + difLevel = length(curBranch) - i + curBranch = curBranch[1:i] + covered = covered[1:i] + # @info "Pruning $difLevel levels from level $(length(curBranch)) via automorphism." + return difLevel + end end - nInvStar = nodeInv end - first = false - else - # pruning - # if !(first || - # length(nodeInv) <= length(nInvStar) || - # length(nodeInv) <= length(nInv1)) - # @show F - # display(F) - # @show typeof(F) - # global errorFlag = deepcopy(F) - # @show first - # @show nodeInv - # @show nInvStar - # @show nInv1 - # end + end - # assertion wrong! nodeInv may be longer than nInvStar! We want the lexicographically largest hash vector, not smallest. - # @assert first || - # length(nodeInv) <= length(nInvStar) || - # length(nodeInv) <= length(nInv1) - # @assert first || - # nodeInv[1:min(length(nodeInv), length(nInvStar))] >= nInvStar[1:min(length(nodeInv), length(nInvStar))] || - # length(nodeInv) <= length(nInv1) - - @views if prune && !first && - !(nodeInv[1:min(length(nodeInv), length(nInv1))] == nInv1[1:min(length(nodeInv), length(nInv1))]) && - !(nodeInv[1:min(length(nodeInv), length(nInvStar))] >= nInvStar[1:min(length(nodeInv), length(nInvStar))]) - return 0 + # improve v1 and vStar + # if first #|| nodeInv < nInv1 + if first || nodeInv < nInv1 + for i in 1:n + v1[coloring[i]] = i end - firstBigCell = Int[] + resize!(nInv1, length(nodeInv)) + nInv1 .= nodeInv + # nInv1 = nodeInv + end + if first || nodeInv > nInvStar for i in 1:n - c = findall(x -> x == i, coloring) - if length(c) > 1 - firstBigCell = c - break - end + vStar[coloring[i]] = i end + resize!(nInvStar, length(nodeInv)) + nInvStar .= nodeInv + # nInvStar = nodeInv + end + first = false + else + # pruning + # if !(first || + # length(nodeInv) <= length(nInvStar) || + # length(nodeInv) <= length(nInv1)) + # @show F + # display(F) + # @show typeof(F) + # global errorFlag = deepcopy(F) + # @show first + # @show nodeInv + # @show nInvStar + # @show nInv1 + # end + + # assertion wrong! nodeInv may be longer than nInvStar! We want the lexicographically largest hash vector, not smallest. + # @assert first || + # length(nodeInv) <= length(nInvStar) || + # length(nodeInv) <= length(nInv1) + # @assert first || + # nodeInv[1:min(length(nodeInv), length(nInvStar))] >= nInvStar[1:min(length(nodeInv), length(nInvStar))] || + # length(nodeInv) <= length(nInv1) + + @views if prune && + !first && + !( + nodeInv[1:min(length(nodeInv), length(nInv1))] == + nInv1[1:min(length(nodeInv), length(nInv1))] + ) && + !( + nodeInv[1:min(length(nodeInv), length(nInvStar))] >= + nInvStar[1:min(length(nodeInv), length(nInvStar))] + ) + return 0 + end + firstBigCell = Int[] + for i in 1:n + c = findall(x -> x == i, coloring) + if length(c) > 1 + firstBigCell = c + break + end + end - push!(covered, Int[]) - for i in firstBigCell - if i in covered[end] && prune - # @info "Prune $(vcat(curBranch,[i])) via automorphism" - continue - end - newC = copy(coloring) - newNodeInv = refine!(newC, i) - push!(curBranch, i) - t = investigateNode(newC, vcat(nodeInv, newNodeInv)) - if t > 0 - return t - 1 - end - pop!(curBranch) + push!(covered, Int[]) + for i in firstBigCell + if prune && i in covered[end] + # @info "Prune $(vcat(curBranch,[i])) via automorphism" + continue + end + newC::Vector{Int} = copy(coloring) + newNodeInv::Vector{UInt} = refine!(newC, F, i) + push!(curBranch, i) + catNodeInv::Vector{UInt} = vcat(nodeInv, newNodeInv) + t = investigateNode(F,newC, catNodeInv,n,nInv1, nInvStar,autG,v1,vStar,curBranch, covered,prune) + + if t > 0 + return t - 1 + end + pop!(curBranch) - union!(covered[end], [i]) - orbit!(stabilizer!(autG, curBranch), covered[end]) + # union!(covered[end], [i]) + if prune || !(i in covered[end]) + push!(covered[end], i) end - pop!(covered) + H = stabilizer!(autG, curBranch) + @assert H !== nothing + orbit!(H, covered[end]) end + pop!(covered) + end - return 0 + return 0 +end + +# Custom implementation of nauty, but still sloooow +# Allows us to label near-arbitrary combinatoric objects +function label(F::T; prune=true, removeIsolated=true) where {T} + # @show F + # display(F) + # @show isolatedVertices(F) + + if size(F) > 0 && removeIsolated + isoV = isolatedVertices(F) + # @show typeof(isoV) + if any(isoV) + # @show F + # @show isolatedVertices(F) + # @show isoV + # @show typeof(isoV) + # @show typeof(map(!,isoV)) + + # F = T(subFlag(F, map(!,isoV))) + return label( + subFlag(F, map(!, isoV)); prune=prune, removeIsolated=removeIsolated + ) + end + end + + n::Int = size(F) + + if n == 0 + return F, Group(), [] end - col = Int[vertexColor(F, v) for v in 1:n] - alpha[1] = true - refine!(col) + nInv1::Vector{UInt} = UInt[] + nInvStar::Vector{UInt} = UInt[] + # p1 = Int[] + # pStar = Int[] - investigateNode(col) - p = [findfirst(x -> x == i, vStar) for i in 1:n] + autG::Group = Group() + # first = true + + v1::Vector{Int} = zeros(Int, n) + vStar::Vector{Int} = zeros(Int, n) + + curBranch::Vector{Int} = Int[] + covered::Vector{Vector{Int}} = Vector{Int}[] + + # To avoid allocations in refine! + # alpha = BitVector(undef, n) + # vertDistinguish = zeros(UInt, n) + + # Q = zeros(UInt, n, n) + # WCellInd = BitVector(undef, n) + # curCell = BitVector(undef, n) + # newCells = Dict{UInt,Int}() + + + col::Vector{Int} = Int[vertexColor(F, v) for v in 1:n] + # alpha[1] = true + refine!(col, F, 0) + + investigateNode(F,col,UInt[],n,nInv1, nInvStar,autG,v1,vStar,curBranch, covered,prune) + + p = zeros(Int, n) + p[vStar] .= 1:n + # p = Int[findfirst(x -> x == i, vStar) for i in 1:n] + # @assert p == p2 return permute(F, p), permute!(autG, p), p end -function generateAll(::Type{F}, maxVertices::Int, maxEdges::Int; withProperty=(F::F) -> true) where {F} +function generateAll( + ::Type{F}, maxVertices::Int, maxEdges::Int; withProperty=(F::F) -> true +) where {F} return generateAll(F, maxVertices, [maxEdges]; withProperty=withProperty) end @@ -293,15 +367,20 @@ end Generates all Flags of type `F` with up to `maxVertices` vertices and up to `maxPredicates` non-zero predicate values. 'maxPredicates' is a vector, for the case that there are multiple predicates. If a function `withProperty:F->{true, false}` is given, keeps adding edges to flags as long as the property holds. """ -function generateAll(::Type{T}, maxVertices::Int, maxPredicates::Vector; withProperty=(F::T) -> true) where {T} +function generateAll( + ::Type{T}, maxVertices::Int, maxPredicates::Vector; withProperty=(F::T) -> true +) where {T} generatedGraphs = Vector{T}[Vector([one(T)])] for i in 1:maxVertices nextGraphs = T[] - pq = PriorityQueue{EdgeMarkedFlag{T,predicateType(T)},Vector}() + pq = PriorityQueue{EdgeMarkedFlag{T,predicateType(T)},Vector{Int}}() for f in generatedGraphs[i] newF = permute(f, 1:i) + # if newF == nothing + # continue + # end # @show newF # @show size(newF) @@ -313,13 +392,15 @@ function generateAll(::Type{T}, maxVertices::Int, maxPredicates::Vector; withPro # @show sum(currentEdges[length(maxPredicates):end]) # @show maxPredicates - if length(maxPredicates) < length(currentEdges) && (maxPredicates[end] isa Int && maxPredicates[end] <= sum(sum.(currentEdges[length(maxPredicates):end]))) + if length(maxPredicates) < length(currentEdges) && ( + maxPredicates[end] isa Int && + maxPredicates[end] <= sum(sum.(currentEdges[length(maxPredicates):end])) + ) continue end - - fixed = allowMultiEdges(T) ? Vector{Int}[] : [collect(1:(i-1))] - uP = findUnknownPredicates(newF, fixed, maxPredicates) + fixed = allowMultiEdges(T) ? Vector{Int}[] : [collect(1:(i - 1))] + uP::Vector{Vector{predicateType(T)}} = findUnknownPredicates(newF, fixed, maxPredicates) # @show uP uP2 = findUnknownGenerationPredicates(newF, fixed, maxPredicates) if uP2 !== nothing @@ -330,23 +411,27 @@ function generateAll(::Type{T}, maxVertices::Int, maxPredicates::Vector; withPro continue end - for i = 1:length(maxPredicates) + for i in 1:length(maxPredicates) if currentEdges[i] == maxPredicates[i] empty!(uP[i]) end end # @show uP - uP = Vector{Vector{Any}}(uP) + uP = Vector{Vector{predicateType(T)}}(uP) FMarked = EdgeMarkedFlag{T}(newF, uP) for (F, _) in allWaysToAddOneMarked(FMarked) if allowMultiEdges(T) F = EdgeMarkedFlag{T}(F.F, FMarked.marked) end - cP = countEdges(F)[1:(end-1)] + cP = countEdges(F)[1:(end - 1)] # @assert length(maxPredicates) == length(cP) if length(maxPredicates) == length(cP) && all(cP .<= maxPredicates) pq[F] = cP - elseif all(cP[1:length(maxPredicates)-1] .<= maxPredicates[1:end-1]) && maxPredicates[end] isa Int && sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] + elseif all( + cP[1:(length(maxPredicates) - 1)] .<= maxPredicates[1:(end - 1)] + ) && + maxPredicates[end] isa Int && + sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] pq[F] = cP end end @@ -363,25 +448,40 @@ function generateAll(::Type{T}, maxVertices::Int, maxPredicates::Vector; withPro cP = countEdges(FMarked.F) if length(maxPredicates) == length(cP) && all(cP .== maxPredicates) continue - elseif all(cP[1:length(maxPredicates)-1] .== maxPredicates[1:end-1]) && sum(cP[length(maxPredicates):end]) == maxPredicates[end] + elseif all(cP[1:(length(maxPredicates) - 1)] .== maxPredicates[1:(end - 1)]) && + sum(cP[length(maxPredicates):end]) == maxPredicates[end] continue end # if !all(countEdges(FMarked.F) .== maxPredicates) for (F, _) in allWaysToAddOneMarked(FMarked) if allowMultiEdges(T) - F = EdgeMarkedFlag{T}(F.F, FMarked.marked) + # @show FMarked.marked + # @show F.marked + if T <: ProductFlag + for (i, FT) in enumerate(fieldtypes(fieldtypes(T)[1])) + if allowMultiEdges(FT) + union!(F.marked, [e for e in FMarked.marked if e[1] == i]) + end + end + else + F = EdgeMarkedFlag{T}(F.F, FMarked.marked) + end end # @assert sum(countEdges(FMarked.F)) < sum(countEdges(F.F)) if !withProperty(F.F) continue end - cP = countEdges(F)[1:(end-1)] + cP = countEdges(F)[1:(end - 1)] # if all(cP .<= maxPredicates) # pq[F] = cP # end + # @show maxPredicates + # @show cP if length(maxPredicates) == length(cP) && all(cP .<= maxPredicates) pq[F] = cP - elseif all(cP[1:length(maxPredicates)-1] .<= maxPredicates[1:end-1]) && sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] + elseif maxPredicates[end] isa Int && + all(cP[1:(length(maxPredicates) - 1)] .<= maxPredicates[1:(end - 1)]) && + sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] pq[F] = cP end end diff --git a/src/utils/NautyBackup.jl b/src/utils/NautyBackup.jl new file mode 100644 index 0000000..7d8a3ec --- /dev/null +++ b/src/utils/NautyBackup.jl @@ -0,0 +1,475 @@ +export generateAll + +include("SchreierSims.jl") + +using DataStructures + +# function refine!(F::T, coloring::Vector{Int}, alpha)::Vector{UInt} where {T<:Flag} +# # Q .= 0 + +# end + +# function individualize!(coloring::Vector{Int}, v::Int)::Vector{Int} +# return coloring +# end + +function refine!(coloring::Vector{Int}, F, v::Int)::Vector{UInt} + n = size(F) + alpha = BitVector(undef, n) + vertDistinguish = zeros(UInt, n) + alpha .= false + if v == 0 + alpha[1] = true + else + alpha[coloring[v]] = true + # individualize!(coloring, v) + for i in 1:n + if coloring[i] >= coloring[v] && i != v + coloring[i] += 1 + end + end + end + + # nodeInv = refine!(coloring) + Q = zeros(UInt, n, n) + @inbounds while maximum(coloring) < n + W::Int = 0 + for i in 1:length(alpha) + if alpha[i] + W = i + break + end + end + if W == 0 + break + end + # W::Int = findfirst(alpha) + # WCellInd .= coloring .== W + WCellInd::BitVector = coloring .== W + alpha[W] = false + + cell::Int = 1 + while cell <= maximum(coloring) + vertDistinguish .= 0 + # curCell .= coloring .== cell + curCell::BitVector = BitVector(undef, n) + for i = 1:n + curCell[i] = coloring[i] == cell + end + + # empty!(newCells) + newCells = Dict{UInt,Int}() + + for v in findall(curCell)#findall(x->x==cell, coloring) + d::UInt = distinguish(F, v, WCellInd) + vertDistinguish[v] = d + # union!(newCells, [d]) + newCells[d] = get(newCells, d, 0) + 1 + end + + # To make it invariant of initial vertex labeling, sort the hashes + # newCells = Dict(d=>i+cell-1 for (i,d) in enumerate(sort!([u for u in unique(vertDistinguish) if u != 0]))) + + # newCellSizes = Dict(d=>count(x->x==d, vertDistinguish) for d in keys(newCells)) + + numNewCells = length(newCells) - 1 + + for i in eachindex(coloring) + if coloring[i] > cell + coloring[i] += numNewCells + end + end + + # coloring[coloring .> cell] .+= numNewCells + for i in length(alpha):-1:(cell + numNewCells + 1) + alpha[i] = alpha[i - numNewCells] + end + # @views alpha[(cell + numNewCells + 1):end] .= alpha[(cell + 1):(end - numNewCells)] + + alpha[(cell + 1):(cell + numNewCells)] .= true + + newCellsCol::Vector{Pair{UInt,Int}} = collect(newCells) + sort!(newCellsCol; by=x -> x[1]) + + # sorting newCellsCol by hashes breaks things?!? But I need to sort! + + for v in findall(curCell)#findall(x->x==cell, coloring) + # coloring[v] = newCells[vertDistinguish[v]] + + d::UInt = vertDistinguish[v] + firstPos::Int = 0 + for i in 1:length(newCellsCol) + if newCellsCol[i][1] == d + firstPos = i + break + end + end + @assert firstPos > 0 + + # findfirst(x -> x[1] == d, newCellsCol) + coloring[v] = firstPos + cell - 1 + Q[coloring[v], W > cell ? W + numNewCells : W] = d + + # Q = hash(Q, vertDistinguish[v]) + end + + if !alpha[cell] + alpha[cell] = true + + maxCellSize = 0 + maxCellInd = 0 + for (i, (d, s)) in enumerate(newCellsCol) + if s > maxCellSize + maxCellSize = s + maxCellInd = i + end + end + + # Remove biggest new cell + alpha[maxCellInd + cell - 1] = false + end + + cell += length(newCellsCol) + end + end + + cl = maximum(coloring) + nodeInvariant = hash(Q, hash([count(x -> x == i, coloring) for i in 1:cl])) + + if maximum(coloring) == n + # append permuted Flag + # return UInt[nodeInvariant, hash(permute(F, coloring))] + # return nodeInvariant, hash(permute(F, coloring)) + + return [nodeInvariant, hash(permute(F, coloring))] + # return hash(nodeInvariant, hash(permute(F, coloring))) + end + return UInt[nodeInvariant] + # return nodeInv +end + +function investigateNode(F,coloring::Vector{Int}, nodeInv::Vector{UInt},n,nInv1, nInvStar,autG,first,v1,vStar,curBranch, covered,prune)::Int + # @info "investigate node at $coloring, $nodeInv" + if maximum(coloring) == n + + # check for automorphisms + if !first + p2 = coloring + autom1 = v1[p2] + changed = false + if F == permute(F, autom1) + changed = addGen!(autG, autom1) + end + autom2 = vStar[p2] + if F == permute(F, autom2) + changed |= addGen!(autG, autom2) + end + + # if changed, check for automorphism pruning higher up the tree + if changed + # @info "Added onto stabilizer, new order is $(order(autG))" + stabilizer!(autG, curBranch) + for i in 1:length(curBranch) + H = stabilizer!(autG, curBranch[1:(i - 1)]) + + o = covered[i] + orbit!(H, o) + if prune && curBranch[i] in o + # cut this branch multiple levels up + difLevel = length(curBranch) - i + curBranch = curBranch[1:i] + covered = covered[1:i] + # @info "Pruning $difLevel levels from level $(length(curBranch)) via automorphism." + return difLevel + end + end + end + end + + # improve v1 and vStar + # if first #|| nodeInv < nInv1 + if first || nodeInv < nInv1 + for i in 1:n + v1[coloring[i]] = i + end + nInv1 = nodeInv + end + if first || nodeInv > nInvStar + for i in 1:n + vStar[coloring[i]] = i + end + nInvStar = nodeInv + end + first = false + else + # pruning + # if !(first || + # length(nodeInv) <= length(nInvStar) || + # length(nodeInv) <= length(nInv1)) + # @show F + # display(F) + # @show typeof(F) + # global errorFlag = deepcopy(F) + # @show first + # @show nodeInv + # @show nInvStar + # @show nInv1 + # end + + # assertion wrong! nodeInv may be longer than nInvStar! We want the lexicographically largest hash vector, not smallest. + # @assert first || + # length(nodeInv) <= length(nInvStar) || + # length(nodeInv) <= length(nInv1) + # @assert first || + # nodeInv[1:min(length(nodeInv), length(nInvStar))] >= nInvStar[1:min(length(nodeInv), length(nInvStar))] || + # length(nodeInv) <= length(nInv1) + + @views if prune && + !first && + !( + nodeInv[1:min(length(nodeInv), length(nInv1))] == + nInv1[1:min(length(nodeInv), length(nInv1))] + ) && + !( + nodeInv[1:min(length(nodeInv), length(nInvStar))] >= + nInvStar[1:min(length(nodeInv), length(nInvStar))] + ) + return 0 + end + firstBigCell = Int[] + for i in 1:n + c = findall(x -> x == i, coloring) + if length(c) > 1 + firstBigCell = c + break + end + end + + push!(covered, Int[]) + for i in firstBigCell + if i in covered[end] && prune + # @info "Prune $(vcat(curBranch,[i])) via automorphism" + continue + end + newC::Vector{Int} = copy(coloring) + newNodeInv::Vector{UInt} = refine!(newC, F, i) + push!(curBranch, i) + catNodeInv::Vector{UInt} = vcat(nodeInv, newNodeInv) + t = investigateNode(F,newC, catNodeInv,n,nInv1, nInvStar,autG,first,v1,vStar,curBranch, covered,prune) + + if t > 0 + return t - 1 + end + pop!(curBranch) + + union!(covered[end], [i]) + orbit!(stabilizer!(autG, curBranch), covered[end]) + end + pop!(covered) + end + + return 0 +end + +# Custom implementation of nauty, but still sloooow +# Allows us to label near-arbitrary combinatoric objects +function label(F::T; prune=true, removeIsolated=true) where {T} + # @show F + # display(F) + # @show isolatedVertices(F) + + if size(F) > 0 && removeIsolated + isoV = isolatedVertices(F) + # @show typeof(isoV) + if any(isoV) + # @show F + # @show isolatedVertices(F) + # @show isoV + # @show typeof(isoV) + # @show typeof(map(!,isoV)) + + # F = T(subFlag(F, map(!,isoV))) + return label( + subFlag(F, map(!, isoV)); prune=prune, removeIsolated=removeIsolated + ) + end + end + + n::Int = size(F) + + if n == 0 + return F, Group(), [] + end + + nInv1::Vector{UInt} = UInt[] + nInvStar::Vector{UInt} = UInt[] + # p1 = Int[] + # pStar = Int[] + + autG::Group = Group() + + first = true + + v1::Vector{Int} = zeros(Int, n) + vStar::Vector{Int} = zeros(Int, n) + + curBranch::Vector{Int} = Int[] + covered::Vector{Vector{Int}} = Vector{Int}[] + + # To avoid allocations in refine! + # alpha = BitVector(undef, n) + # vertDistinguish = zeros(UInt, n) + + # Q = zeros(UInt, n, n) + # WCellInd = BitVector(undef, n) + # curCell = BitVector(undef, n) + # newCells = Dict{UInt,Int}() + + + col::Vector{Int} = Int[vertexColor(F, v) for v in 1:n] + # alpha[1] = true + refine!(col, F, 0) + + investigateNode(F,col,UInt[],n,nInv1, nInvStar,autG,first,v1,vStar,curBranch, covered,prune) + p = Int[findfirst(x -> x == i, vStar) for i in 1:n] + + return permute(F, p), permute!(autG, p), p +end + +function generateAll( + ::Type{F}, maxVertices::Int, maxEdges::Int; withProperty=(F::F) -> true +) where {F} + return generateAll(F, maxVertices, [maxEdges]; withProperty=withProperty) +end + +#TODO: Quite a bit slower than nauty/traces, how are they doing it? +""" + generateAll(::Type{F}, maxVertices::Int, maxPredicates::Vector{Int}; withProperty = (F::T) -> true) where {F} + +Generates all Flags of type `F` with up to `maxVertices` vertices and up to `maxPredicates` non-zero predicate values. 'maxPredicates' is a vector, for the case that there are multiple predicates. If a function `withProperty:F->{true, false}` is given, keeps adding edges to flags as long as the property holds. +""" +function generateAll( + ::Type{T}, maxVertices::Int, maxPredicates::Vector; withProperty=(F::T) -> true +) where {T} + generatedGraphs = Vector{T}[Vector([one(T)])] + for i in 1:maxVertices + nextGraphs = T[] + + pq = PriorityQueue{EdgeMarkedFlag{T,predicateType(T)},Vector{Int}}() + + for f in generatedGraphs[i] + newF = permute(f, 1:i) + # if newF == nothing + # continue + # end + + # @show newF + # @show size(newF) + push!(nextGraphs, newF) + currentEdges = countEdges(newF) + # @show currentEdges + # @show length(currentEdges) + # @show currentEdges[length(maxPredicates):end] + # @show sum(currentEdges[length(maxPredicates):end]) + # @show maxPredicates + + if length(maxPredicates) < length(currentEdges) && ( + maxPredicates[end] isa Int && + maxPredicates[end] <= sum(sum.(currentEdges[length(maxPredicates):end])) + ) + continue + end + + fixed = allowMultiEdges(T) ? Vector{Int}[] : [collect(1:(i - 1))] + uP::Vector{Vector{predicateType(T)}} = findUnknownPredicates(newF, fixed, maxPredicates) + # @show uP + uP2::Vector{Vector{predicateType(T)}} = findUnknownGenerationPredicates(newF, fixed, maxPredicates) + if uP2 !== nothing + uP = vcat(uP, uP2) + end + + if length(uP) == 1 && length(uP[1]) == 0 + continue + end + + for i in 1:length(maxPredicates) + if currentEdges[i] == maxPredicates[i] + empty!(uP[i]) + end + end + # @show uP + uP = Vector{Vector{predicateType(T)}}(uP) + FMarked = EdgeMarkedFlag{T}(newF, uP) + for (F, _) in allWaysToAddOneMarked(FMarked) + if allowMultiEdges(T) + F = EdgeMarkedFlag{T}(F.F, FMarked.marked) + end + cP = countEdges(F)[1:(end - 1)] + # @assert length(maxPredicates) == length(cP) + if length(maxPredicates) == length(cP) && all(cP .<= maxPredicates) + pq[F] = cP + elseif all( + cP[1:(length(maxPredicates) - 1)] .<= maxPredicates[1:(end - 1)] + ) && + maxPredicates[end] isa Int && + sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] + pq[F] = cP + end + end + end + + while !isempty(pq) + FMarked = dequeue!(pq) + # @show (length(pq), sum(countEdges(FMarked.F))) + # @show FMarked + if withProperty(FMarked.F) + # continue + push!(nextGraphs, FMarked.F) + end + cP = countEdges(FMarked.F) + if length(maxPredicates) == length(cP) && all(cP .== maxPredicates) + continue + elseif all(cP[1:(length(maxPredicates) - 1)] .== maxPredicates[1:(end - 1)]) && + sum(cP[length(maxPredicates):end]) == maxPredicates[end] + continue + end + # if !all(countEdges(FMarked.F) .== maxPredicates) + for (F, _) in allWaysToAddOneMarked(FMarked) + if allowMultiEdges(T) + # @show FMarked.marked + # @show F.marked + if T <: ProductFlag + for (i, FT) in enumerate(fieldtypes(fieldtypes(T)[1])) + if allowMultiEdges(FT) + union!(F.marked, [e for e in FMarked.marked if e[1] == i]) + end + end + else + F = EdgeMarkedFlag{T}(F.F, FMarked.marked) + end + end + # @assert sum(countEdges(FMarked.F)) < sum(countEdges(F.F)) + if !withProperty(F.F) + continue + end + cP = countEdges(F)[1:(end - 1)] + # if all(cP .<= maxPredicates) + # pq[F] = cP + # end + # @show maxPredicates + # @show cP + if length(maxPredicates) == length(cP) && all(cP .<= maxPredicates) + pq[F] = cP + elseif maxPredicates[end] isa Int && + all(cP[1:(length(maxPredicates) - 1)] .<= maxPredicates[1:(end - 1)]) && + sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] + pq[F] = cP + end + end + # end + end + + push!(generatedGraphs, unique(labelCanonically.(nextGraphs))) + end + return unique(vcat(generatedGraphs...)) +end diff --git a/src/utils/SchreierSims.jl b/src/utils/SchreierSims.jl index 720319c..53d03b8 100644 --- a/src/utils/SchreierSims.jl +++ b/src/utils/SchreierSims.jl @@ -10,7 +10,7 @@ mutable struct Group # SGS::Vector{Vector{Int}} # Strong generating set SV::Vector{Int} # Schreier Vector O::Vector{Int} # Orbit of b - order::AbstractVector{Int} # target order of basis elements (usually 1:n, does get modified when point-wise stabilizers are searched for) + order::Vector{Int} # target order of basis elements (usually 1:n, does get modified when point-wise stabilizers are searched for) subGroup::Union{Group,Nothing} Group() = new(-1, Vector{Int}[], -1, Int[], Int[], Int[], nothing) @@ -68,7 +68,7 @@ function orbit!(G::Group, vs::Vector{Int}) end # similar to orbit, but uses inverses of G.gen, and records Schreier Vector -function orbitSchreier(G::Group) +function orbitSchreier!(G::Group) n = G.n empty!(G.O) push!(G.O, G.b) @@ -79,9 +79,11 @@ function orbitSchreier(G::Group) @inbounds w = G.O[i] for (k, p) in enumerate(G.gen) j = findfirst(x -> x == w, p) - if !(j in G.O) - push!(G.O, j) - G.SV[j] = k + if j !== nothing # should always be true, here for type stability + if !(j in G.O) + push!(G.O, j) + G.SV[j] = k + end end end i += 1 @@ -96,7 +98,7 @@ function findInvRepr(G::Group, v::Int) while v != G.b k = G.SV[v] if k == 0 - return nothing + return Int[] end p = G.gen[k] g .= p[g] @@ -105,7 +107,11 @@ function findInvRepr(G::Group, v::Int) return g end -function schreier_sims!(G) +function schreier_sims!(G0::Union{Group,Nothing}) + if G0 === nothing + return nothing + end + G::Group = G0 if length(G.gen) == 0 return nothing end @@ -127,9 +133,11 @@ function schreier_sims!(G) @assert G.b != -1 if oldB != G.b - G.subGroup = Group() - G.subGroup.n = G.n - G.subGroup.order = G.order + G2 = Group() + G2.n = G.n + G2.order = G.order + G.subGroup = G2 + # G.subGroup = Group() elseif G.subGroup.order != G.order G.subGroup.order = G.order schreier_sims!(G.subGroup) @@ -139,7 +147,7 @@ function schreier_sims!(G) # G.SGS = G.gen - orbitSchreier(G) + orbitSchreier!(G) # @show G.SV @@ -147,7 +155,7 @@ function schreier_sims!(G) for i in G.O r = findInvRepr(G, i) - if r !== nothing + if length(r) == G.n for s in G.gen # rs = inv(r)*s # rs s[r] @@ -168,13 +176,16 @@ function schreier_sims!(G) end end -function sift(G::Group, p::Vector{Int}) - if length(G.gen) == 0 || p == 1:length(p)#isone(p) +function sift(G::Union{Group, Nothing}, p::Vector{Int}) + if G === nothing + return p + end + if length(G.gen) == 0 || issorted(p)# == 1:length(p) return p end gInv = findInvRepr(G, p[G.b]) - if gInv !== nothing + if length(gInv) == G.n # p *= gInv p .= gInv[p] end @@ -186,10 +197,13 @@ function sift(G::Group, p::Vector{Int}) return p end -sift(::Nothing, p::Vector{Int}) = p +# sift(::Nothing, p::Vector{Int}) = p # adds p to G, returns true if changed -function addGen!(G::Group, p::Vector{Int}) +function addGen!(G::Union{Group, Nothing}, p::Vector{Int}) + if G === nothing + return false + end q = sift(G, p) # if !isone(q) if q != 1:(G.n) @@ -248,7 +262,10 @@ function stabilizer(G::Group, S::Vector{Int}) end # returns same group as stabilizer(G, S), but modifies stabilizer chain of G in the process -function stabilizer!(G::Group, S::Vector{Int}, keepOrder=false) +function stabilizer!(G::Union{Group, Nothing}, S::Vector{Int}, keepOrder=false) + if G === nothing + return G + end if keepOrder order = vcat(S, setdiff(1:(G.n), S)) G.order = order @@ -260,21 +277,28 @@ function stabilizer!(G::Group, S::Vector{Int}, keepOrder=false) else # May change order of elements of S in the stabilizer chain covered = Int[] - while G.b in S - push!(covered, G.b) - G = G.subGroup + G2 = G + while G2.b in S + push!(covered, G2.b) + G2 = G2.subGroup end - order = vcat(covered, setdiff(S, covered), setdiff(1:(G.n), S)) - G.order = order - schreier_sims!(G) - while G.b in S - G = G.subGroup + if G2 === nothing + return G2 end - return G + order = vcat(covered, setdiff(S, covered), setdiff(1:(G2.n), S)) + G2.order = order + schreier_sims!(G2) + while G2.b in S + G2 = G2.subGroup + end + return G2 end end -function permute!(gr::Group, per::Vector{Int}) +function permute!(gr::Union{Group,Nothing}, per::Vector{Int}) + if gr === nothing + return gr + end if length(gr.gen) == 0 return gr end diff --git a/test/src/nauty.jl b/test/src/nauty.jl index e69de29..6059a3b 100644 --- a/test/src/nauty.jl +++ b/test/src/nauty.jl @@ -0,0 +1,79 @@ +using FlagSOS +using JET +# using ProfileCanvas +using Profile +# using ProfileView + +import .VSCodeServer.view_profile +function view_profile(data = Profile.fetch(); C=false, kwargs...) + d = Dict{String, VSCodeServer.ProfileFrame}() + + if VERSION >= v"1.8.0-DEV.460" + threads = ["all", 1:Threads.nthreads()...] + else + threads = ["all"] + end + + if isempty(data) + Profile.warning_empty() + return + end + + lidict = Profile.getdict(unique(data)) + data_u64 = convert(Vector{UInt64}, data) + for thread in threads + graph = VSCodeServer.stackframetree(data_u64, lidict; thread=thread, kwargs...) + d[string(thread)] = VSCodeServer.make_tree( + VSCodeServer.ProfileFrame( + "root", "", "", 0, graph.count, missing, 0x0, missing, VSCodeServer.ProfileFrame[] + ), graph; C=C) + end + + VSCodeServer.JSONRPC.send(VSCodeServer.conn_endpoint[], VSCodeServer.repl_showprofileresult_notification_type, (; trace=d, typ="Thread")) +end + + +@time Gs = generateAll(Graph, 6, 1000) +@time Gs = generateAll(Graph, 7, 1000) # 3.2s +@profview Gs = generateAll(Graph, 7, 1000) recur=:flat +@profview Gs = generateAll(Graph, 7, 1000) format=:flat + + +@profview_allocs Gs = generateAll(Graph, 7, 1000) recur=:flat + + +@time Gs = generateAll(Graph, 8, 1000) # 57s +@time Gs = generateAll(Graph, 9, 1000) # 57s + +# Profile.clear() +# @profile Gs = generateAll(Graph, 7, 1000) +# view_profile(Profile.fetch(); recur=:flat) +# Profile.print(format=:flat) + +# in REPL +@report_opt generateAll(Graph, 3, 1000) +@report_call generateAll(Graph, 5, 1000) + +function test(A::Vector{Int}) + push!(A, 1) +end + +tmp = [1,2] + +tmp2 = FlagSOS.Group() + +function test2(G::FlagSOS.Group) + G.n = 1 +end + +test2(tmp2) + +tmp2.n + +function test3(A::Vector{Int}) + A = [1,2,3] +end + +tmp = [1,2] +test3(tmp) +tmp \ No newline at end of file From b0d64079a6375c3d20a0057b5120b67036603a57 Mon Sep 17 00:00:00 2001 From: Daniel Brosch <73886037+DanielBrosch@users.noreply.github.com> Date: Wed, 6 Nov 2024 16:55:26 +0100 Subject: [PATCH 2/2] small fixes --- Project.toml | 1 + src/FlagAlgebras/AbstractFlagAlgebra.jl | 18 ++- src/FlagAlgebras/DirectedGraphs.jl | 1 + src/FlagAlgebras/EdgeMarkedFlags.jl | 1 + src/FlagAlgebras/Graphs.jl | 31 +++- src/FlagAlgebras/InducedFlags.jl | 43 ++++-- src/FlagAlgebras/PartiallyLabeledFlags.jl | 8 +- src/FlagAlgebras/ProductFlag.jl | 7 +- src/FlagAlgebras/QuantumFlags.jl | 8 +- src/FlagModels/FlagModel.jl | 89 +++++++++-- src/FlagModels/LasserreModel.jl | 4 +- src/FlagModels/RazborovModel.jl | 170 +++++++++++++--------- src/utils/Nauty.jl | 69 ++++----- src/utils/PermutationModuleQuotients.jl | 12 +- test/src/nauty.jl | 4 +- 15 files changed, 309 insertions(+), 157 deletions(-) diff --git a/Project.toml b/Project.toml index ea6cc77..f68baad 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.1.2" [deps] AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" +ClusteredLowRankSolver = "cadeb640-247c-4e6d-a8b8-270aef5e1f95" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" diff --git a/src/FlagAlgebras/AbstractFlagAlgebra.jl b/src/FlagAlgebras/AbstractFlagAlgebra.jl index 6dfd1cc..e3470bb 100644 --- a/src/FlagAlgebras/AbstractFlagAlgebra.jl +++ b/src/FlagAlgebras/AbstractFlagAlgebra.jl @@ -181,6 +181,9 @@ end function glueFinite_internal(N, F::T, G::T, p::AbstractVector{Int}; labelFlags=true) where {T<:Flag} if N == :limit res = glue(F, G, p) + if res === nothing + return QuantumFlag{T, Rational{Int64}}() + end if labelFlags return labelCanonically(res) end @@ -315,13 +318,13 @@ function findUnknownPredicates( F::T, fixed::Vector{U}=Vector{Int}[], predLimits::Vector = Int[] ) where {T<:Flag,U<:AbstractVector{Int}} error("findUnknownPredicates is not defined for Flag type $T") - return missing + return nothing end function findUnknownGenerationPredicates( F::T, fixed::Vector{U}=Vector{Int}[], predLimits::Vector = Int[] -) where {T<:Flag,U<:AbstractVector{Int}} - return nothing + ) where {T<:Flag,U<:AbstractVector{Int}} + return predicateType(T)[] end """ @@ -457,6 +460,15 @@ function allowMultiEdges(::Type{T}) where {T<:Flag} return false end +import Base.^ +function ^(F::T, i::Int) where {T<:Flag} + @assert i >= 0 + if i == 0 + return one(T) + end + return F * F^(i-1) +end + include("InducedFlags.jl") include("Graphs.jl") include("ConstantWeightCodes.jl") diff --git a/src/FlagAlgebras/DirectedGraphs.jl b/src/FlagAlgebras/DirectedGraphs.jl index 1473a03..b9f88b8 100644 --- a/src/FlagAlgebras/DirectedGraphs.jl +++ b/src/FlagAlgebras/DirectedGraphs.jl @@ -102,6 +102,7 @@ function glue( # res = BitMatrix(zeros(Bool, n, n)) res = BitMatrix(undef, n, n) + res .= 0 @views res[1:n2, 1:n2] = g2.A @views res[p[1:n1], p[1:n1]] .|= g1.A diff --git a/src/FlagAlgebras/EdgeMarkedFlags.jl b/src/FlagAlgebras/EdgeMarkedFlags.jl index 7e9e674..86cf016 100644 --- a/src/FlagAlgebras/EdgeMarkedFlags.jl +++ b/src/FlagAlgebras/EdgeMarkedFlags.jl @@ -181,6 +181,7 @@ function moebius(F::EdgeMarkedFlag{T,P}; label=false) where {T<:Flag,P<:Predicat tmp2 = Dict{EdgeMarkedFlag{T,P},Int}() for flippedEdges in 0:k + # @show flippedEdges for (F2, c2) in tmp res += c2 * (-1)^flippedEdges * F2.F for (F3, c3) in allWaysToAddOneMarked(F2) diff --git a/src/FlagAlgebras/Graphs.jl b/src/FlagAlgebras/Graphs.jl index 1849dc7..465d062 100644 --- a/src/FlagAlgebras/Graphs.jl +++ b/src/FlagAlgebras/Graphs.jl @@ -1,4 +1,4 @@ -export Graph +export Graph, connectedComponents """ $(TYPEDEF) @@ -234,3 +234,32 @@ function isolatedVertices(F::Graph)::BitVector # return .!vec(any(F.A; dims=1)) return res end + +function connectedComponents(G::Graph)::Vector{Graph} + n = size(G) + + if n == 0 + return [G] + end + + components = zeros(Int, n) + nComp = 0 + for i = 1:n + if components[i] == 0 + nComp += 1 + components[i] = nComp + end + for j = i+1:n + if G.A[i,j] == 1 + if components[j] == 0 + components[j] = components[i] + else + c = minimum(components[[i,j]]) + d = maximum(components[[i,j]]) + components[components.==d] .= c + end + end + end + end + return [subFlag(G, findall(x->x==u, components)) for u in unique(components)] +end \ No newline at end of file diff --git a/src/FlagAlgebras/InducedFlags.jl b/src/FlagAlgebras/InducedFlags.jl index 07bd06a..d875f8e 100644 --- a/src/FlagAlgebras/InducedFlags.jl +++ b/src/FlagAlgebras/InducedFlags.jl @@ -202,28 +202,39 @@ end # E.g. adding isolated vertices results in a quantum flag equivalent to the original flag function quotient(Fs::Vector{T}, isAllowed=(f) -> true) where {T<:Flag} oneVert = permute(T(), [1]) - @show oneVert + # @show oneVert - @show Fs + # @show Fs - n = maximum(size.(Fs); init = 1) + n = maximum(size.(Fs); init=1) - res = QuantumFlag{T, Rational{Int}}[] + res = QuantumFlag{T,Rational{Int}}[] # res = Dict() for f in Fs - size(f) == n && continue - tmp = oneVert * f - if any(tmp.coeff) do (G, _) - isAllowed(G) && !in(G, Fs) + for newVerts in 1:2 + size(f) == n && continue + if newVerts == 1 + tmp = labelCanonically(oneVert * f - 1//1 * f) + else + tmp = labelCanonically(oneVert * (oneVert * f) - 1//1 * f) + end + # @show tmp + if any(tmp.coeff) do (G, _) + isAllowed(G) && !in(G, Fs) + end + @info "Missing some flags to eliminate product of $f" + continue + end + # res[f] = tmp + filter!(x -> isAllowed(x.first), tmp.coeff) + # @show length(res) + 1 + # @show f + # @show oneVert * f + if !iszero(tmp) + push!(res, labelCanonically(tmp)) + end + break end - @info "Missing some flags to eliminate product of $f" - continue - end - # res[f] = tmp - @show length(res) + 1 - @show f - @show oneVert * f - push!(res, tmp - 1//1 * f) end return res # n = length(Fs) diff --git a/src/FlagAlgebras/PartiallyLabeledFlags.jl b/src/FlagAlgebras/PartiallyLabeledFlags.jl index a686d5a..0ad215a 100644 --- a/src/FlagAlgebras/PartiallyLabeledFlags.jl +++ b/src/FlagAlgebras/PartiallyLabeledFlags.jl @@ -40,7 +40,9 @@ function permute( )::PartiallyLabeledFlag{T} where {T<:Flag} F.n > 0 && @assert 1:(F.n) == p[1:(F.n)] - return PartiallyLabeledFlag{T}(permute(F.F, p), F.n) + tmp = permute(F.F, p) + + return PartiallyLabeledFlag{T}(tmp, F.n) end function permute(pred::LabelPredicate, p::AbstractVector{Int}) @@ -101,7 +103,7 @@ function glue( FG === nothing && return nothing if FG isa QuantumFlag - return QuantumFlag(PartiallyLabeledFlag{T}(f, max(F.n, G.n)) => d for (f, d) in FG.coeff) + return QuantumFlag{PartiallyLabeledFlag{T}, typeof(FG).parameters[2]}(PartiallyLabeledFlag{T}(f, max(F.n, G.n)) => d for (f, d) in FG.coeff) end return PartiallyLabeledFlag{T}(FG, max(F.n, G.n)) end @@ -245,7 +247,7 @@ function toInduced(F::Union{PartiallyLabeledFlag{T}, QuantumFlag{PartiallyLabele tmp = zeta(F) res = QuantumFlag{PartiallyLabeledFlag{InducedFlag{T}}, Int}() for (G,c) in tmp.coeff - res += c*PartiallyLabeledFlag{InducedFlag{T}}(InducedFlag(G.F), G.n) + res += c*PartiallyLabeledFlag{InducedFlag{T}}(InducedFlag{T}(G.F), G.n) end return res end diff --git a/src/FlagAlgebras/ProductFlag.jl b/src/FlagAlgebras/ProductFlag.jl index 90cf659..49933d0 100644 --- a/src/FlagAlgebras/ProductFlag.jl +++ b/src/FlagAlgebras/ProductFlag.jl @@ -74,7 +74,7 @@ function findUnknownPredicates( if length(predLimits) == length(fieldtypes(FT)) FIP = findUnknownPredicates(F.Fs[i], fixed, predLimits[i]) else - FIP = findUnknownPredicates(F.Fs[i], fixed, [predLimits[1]]) + FIP = findUnknownPredicates(F.Fs[i], fixed, predLimits[1]) end # @assert length(FIP) == 1 @@ -101,7 +101,7 @@ function findUnknownGenerationPredicates( if length(predLimits) == length(fieldtypes(FT)) FIP = findUnknownGenerationPredicates(F.Fs[i], fixed, predLimits[i]) else - FIP = findUnknownGenerationPredicates(F.Fs[i], fixed, [predLimits[1]]) + FIP = findUnknownGenerationPredicates(F.Fs[i], fixed, predLimits[1]) end # FIP = findUnknownGenerationPredicates(F.Fs[i], fixed, predLimits[i]) # @assert length(FIP) == 1 @@ -217,4 +217,5 @@ function vertexColor(F::ProductFlag{FT}, v::Int) where {FT} # @assert length(unique(cs)) == 1 # return cs[1] -end \ No newline at end of file +end + diff --git a/src/FlagAlgebras/QuantumFlags.jl b/src/FlagAlgebras/QuantumFlags.jl index eda5692..bafdfe1 100644 --- a/src/FlagAlgebras/QuantumFlags.jl +++ b/src/FlagAlgebras/QuantumFlags.jl @@ -59,7 +59,10 @@ The maximum number of edges of the flags in 'F'. function countEdges(F::QuantumFlag) length(F.coeff) == 0 && return 0 - return maximum([countEdges(f) for f in keys(F.coeff)]) + edgeCounts = [countEdges(f) for f in keys(F.coeff)] + k = length(edgeCounts[1]) + + return [maximum([e[i] for e in edgeCounts]) for i = 1:k] end """ @@ -206,7 +209,8 @@ function Base.:*(F::T, G::QuantumFlag{T,R}) where {T<:Flag,R<:Real} for (B, d) in G.coeff AB = F * B AB === nothing && continue - res.coeff[AB] = get(res.coeff, AB, zero(R)) + d + # res.coeff[AB] = get(res.coeff, AB, zero(R)) + d + res += d*AB end filter!(p -> !iszero(p.second), res.coeff) diff --git a/src/FlagModels/FlagModel.jl b/src/FlagModels/FlagModel.jl index 8dc43e1..570ba6e 100644 --- a/src/FlagModels/FlagModel.jl +++ b/src/FlagModels/FlagModel.jl @@ -6,7 +6,10 @@ export FlagModel, addEquality!, buildStandardModel, addRazborovBlock!, - addBinomialBlock! + addBinomialBlock!, + buildClusteredLowRankModel + +using ClusteredLowRankSolver """ FlagModel{T <: Flag, N, D} <: AbstractFlagModel{T, N, D} @@ -39,7 +42,7 @@ function Base.show(io::IO, m::FlagModel{T,N,D}) where {T,N,D} if m.objective !== nothing println(io, "Objective: $(m.objective)") end - Base.show.(io, m.subModels) + return Base.show.(io, m.subModels) end function addForbiddenFlag!(m::FlagModel{T,N,D}, F::T) where {T<:Flag,N,D} @@ -97,7 +100,9 @@ function addLasserreBlock!( return lM end -function addRazborovBlock!(m::FlagModel{T,N,D}, lvl; maxLabels=lvl, maxBlockSize=Inf) where {T<:Flag,N,D} +function addRazborovBlock!( + m::FlagModel{T,N,D}, lvl; maxLabels=lvl, maxBlockSize=Inf +) where {T<:Flag,N,D} rM = RazborovModel{T,N,D}(m) push!(m.subModels, rM) computeRazborovBasis!(rM, lvl; maxLabels=maxLabels, maxBlockSize=maxBlockSize) @@ -187,7 +192,7 @@ function addEquality!( g::QuantumFlag{PartiallyLabeledFlag{T},D}, maxEdges; maxVertices=size(g) + - floor((maxEdges - countEdges(g)[2]) / 2) * maxPredicateArguments(T), + (maxEdges - countEdges(g)[2]) * maxPredicateArguments(T), ) where {T<:Flag,N,D} gl = labelCanonically(g) @@ -220,7 +225,7 @@ function addEquality!( g::QuantumFlag{T,D}, maxEdges; maxVertices=size(g) + - floor((maxEdges - countEdges(g)[1]) / 2) * maxPredicateArguments(T), + (maxEdges - countEdges(g)[1]) * maxPredicateArguments(T), ) where {T<:Flag,N,D} gl = labelCanonically(g) @@ -260,6 +265,7 @@ function buildJuMPModel( end if addBoundVars + @warn "Adding bound variables" for F in keys(variables) fl = @variable(jumpModel, base_name = "lower$F", lower_bound = 0) fu = @variable(jumpModel, base_name = "upper$F", lower_bound = 0) @@ -330,7 +336,9 @@ function modelBlockSizes(m::FlagModel) return res end + function buildStandardModel(m::FlagModel{T,N,D}) where {T<:Flag,N,D} + #TODO: Quotient when using InducedFlags obj = labelCanonically(m.objective) vars = union([collect(keys(sM.sdpData)) for sM in m.subModels]...) filter!(F -> isAllowed(m, F), vars) @@ -358,23 +366,72 @@ function buildStandardModel(m::FlagModel{T,N,D}) where {T<:Flag,N,D} return (obj=obj, vars=vars, blocks=blocks, blockSizes=blockSizes) end +function buildClusteredLowRankModel(m::FlagModel{T,N,D}) where {T,N,D} + + obj, vars, blocks, blockSizes = buildStandardModel(m) + o = T() + + + varDicts = Dict( + G => Dict{Any, Any}( + mu => Matrix(B[G]) + for (mu, B) in blocks if haskey(B, G) && blockSizes[mu] > 0 + ) + for G in vars + ) + + # bound variables + # for G in vars + # if G != o + # varDicts[G][(:lower, G)] = [1;;] + # varDicts[G][(:upper, G)] = [-1;;] + # varDicts[o][(:upper, G)] = [1;;] + # end + # end + + freeVarDicts = Dict( + G => Dict( + mu => B[G][1, 1] + for (mu, B) in blocks if haskey(B, G) && blockSizes[mu] < 0 + ) + for G in vars + ) + + + + clObj = Objective(get(obj.coeff, o, 0), get(varDicts, o, Dict()), freeVarDicts[o]) + clCons = Constraint[] + + for G in vars + G == o && continue + # @show G + # display(varDicts[G]) + # display(freeVarDicts[G]) + c = Constraint(get(obj.coeff, G, 0), varDicts[G], freeVarDicts[G]) + # @show c + push!(clCons, c) + end + + return Problem(Minimize(clObj), clCons) +end + function facialReduction(m::AbstractFlagModel) return error("TODO") end -function verifySOS( - m::FlagModel{T,N,D}, sol::Vector; io::IO=stdout -) where {T,N,D} +function verifySOS(m::FlagModel{T,N,D}, sol::Vector; io::IO=stdout) where {T,N,D} @assert length(sol) == length(m.subModels) Base.show(io, m) println() - - tmp = labelCanonically(sum(verifySOS(m.subModels[i], sol[i]; io=io) for i in 1:length(m.subModels))) - + + tmp = labelCanonically( + sum(verifySOS(m.subModels[i], sol[i]; io=io) for i in 1:length(m.subModels)) + ) + println(io, "SOS result:") println(io, "$(tmp) >= 0") - + err = Rational{BigInt}(0) constTerm = Rational{BigInt}(0) for (F, c) in tmp.coeff @@ -393,20 +450,20 @@ function verifySOS( println(io, "\nSummary:") Base.show(io, m) - + println(io, "Constant:") println(io, "$(constTerm)") - + println(io, "Error:") println(io, "$(err)") - + if haskey(m.objective.coeff, T()) constTerm -= m.objective.coeff[T()] end res = constTerm + err + m.objective println(io, "Final rigorous bound:") println(io, "$(res) >= 0") - + println(io, "Floating point bound (rounded appropriately):") first = true for (g, cR) in res.coeff diff --git a/src/FlagModels/LasserreModel.jl b/src/FlagModels/LasserreModel.jl index 8ddaf03..6a84122 100644 --- a/src/FlagModels/LasserreModel.jl +++ b/src/FlagModels/LasserreModel.jl @@ -294,7 +294,7 @@ function multiplyPolytabsAndSymmetrize( (newVariant, fact) = symPolytabloidProduct(sp1.T, sp2.T, la, limit) - combinedOverlaps = Dict([]) + combinedOverlaps = Dict{Matrix{Int}, D}() for (a, b) in newVariant # cord = [2:size(a, 1)..., 1] # shiftedMat = a[cord, cord] @@ -328,7 +328,7 @@ function multiplyPolytabsAndSymmetrize( end # reduce using automorphisms - combinedOverlapsReduced = Dict{Any,D}([]) + combinedOverlapsReduced = Dict{Matrix{Int}, D}() if useGroups @assert limit "TODO: Fix finite case with group speedup" #TODO Current solution does reduce it somewhat, but not fully? diff --git a/src/FlagModels/RazborovModel.jl b/src/FlagModels/RazborovModel.jl index 895c552..5f4bfc5 100644 --- a/src/FlagModels/RazborovModel.jl +++ b/src/FlagModels/RazborovModel.jl @@ -12,7 +12,7 @@ A (not fully symmetry reduced) Razborov-style model. If T is an InducedFlag, the mutable struct RazborovModel{T<:Flag,N,D} <: AbstractFlagModel{T,N,D} basis::Dict{T,Vector{PartiallyLabeledFlag{T}}} blockSymmetry::Dict - sdpData::Dict{T,Dict{T,SparseMatrixCSC{D,Int}}} + sdpData::Dict{T,Dict{Union{T, String},SparseMatrixCSC{D,Int}}} parentModel quotient @@ -53,11 +53,15 @@ end function modelSize(m::RazborovModel) # return Partition([m.blockSymmetry[b].n for b in keys(m.basis)]) - return Partition([length(b) for b in values(m.basis)]) + return Partition(vcat([length(b) for b in values(m.basis)], [1 for _ in m.quotient])) end function modelBlockSizes(m::RazborovModel) - return Dict(F => length(b) for (F, b) in m.basis) + res = Dict{Any,Int}(F => length(b) for (F, b) in m.basis) + for (i, F) in enumerate(m.quotient) + res["Q$i"] = -1 + end + return res end function computeUnreducedRazborovBasis( @@ -350,15 +354,33 @@ function computeSDP!(m::RazborovModel{T,N,D}, reservedVerts::Int) where {T,N,D} # Eliminate linear dependencies @info "Eliminating linear dependencies" - # Fs = sort(collect(keys(m.sdpData)), by = size) + Fs::Vector{T} = sort(collect(keys(m.sdpData)); by=size) n = maximum(size.(keys(m.sdpData))) - Fs = generateAll(T, n, [99999]) + union!(Fs, generateAll(T, n, [99999])) + + # expandedFs = T[] + # for F in Fs + # for m in size(F):n + # push!(expandedFs, labelCanonically(permute(F, 1:m))) + # end + # end + + # reduction = quotient(expandedFs, x -> isAllowed(m.parentModel, x)) reduction = quotient(Fs, x -> isAllowed(m.parentModel, x)) display(reduction) m.quotient = reduction + for (i, F) in enumerate(m.quotient) + for (G,c) in F.coeff + if !haskey(m.sdpData, G) + m.sdpData[G] = Dict() + end + m.sdpData[G]["Q$i"] = D[c;;] + end + end + # for i in size(reduction, 1):-1:1 # j = findlast(x -> !iszero(x), reduction[i, :]) # @assert reduction[i, j] == 1 @@ -437,56 +459,65 @@ function buildJuMPModel(m::RazborovModel, replaceBlocks=Dict(), jumpModel=Model( constraints = Dict() for (mu, n) in b - P = m.blockSymmetry[mu].pattern - name = "y[$mu]" - # t = maximum(P) - t = m.blockSymmetry[mu].n - if haskey(m.blockSymmetry[mu], :reg) - reg = m.blockSymmetry[mu].reg - y = @variable(jumpModel, [keys(reg)], base_name = name) - AT = typeof(1 * y[1]) - - # Y[mu] = zeros(AT, t, t) - # for s in keys(reg) - # Y[mu] .+= reg[s]*y[s] - # end - Y[mu] = sum(reg[s] * y[s] for s in keys(reg)) - if size(P, 1) > 1 - constraints[name] = @constraint(jumpModel, Y[mu] in PSDCone()) + if n > 0 + P = m.blockSymmetry[mu].pattern + name = "y[$mu]" + # t = maximum(P) + t = m.blockSymmetry[mu].n + if haskey(m.blockSymmetry[mu], :reg) + reg = m.blockSymmetry[mu].reg + y = @variable(jumpModel, [keys(reg)], base_name = name) + AT = typeof(1 * y[1]) + + # Y[mu] = zeros(AT, t, t) + # for s in keys(reg) + # Y[mu] .+= reg[s]*y[s] + # end + Y[mu] = sum(reg[s] * y[s] for s in keys(reg)) + if size(P, 1) > 1 + constraints[name] = @constraint(jumpModel, Y[mu] in PSDCone()) + else + constraints[name] = @constraint(jumpModel, Y[mu] .>= 0) + end else - constraints[name] = @constraint(jumpModel, Y[mu] .>= 0) + # numVars = maximum(P) + # y = @variable(jumpModel, [1:numVars], base_name = name) + # AT = typeof(1 * y[1]) + + # Y[mu] = zeros(AT, t, t) + # # Y[mu] = zeros(AT, size(P)) + # for s in 1:numVars + # Y[mu][P .== s] .+= y[s] + # end + + # Y[mu] = sum((P .== s) * y[s] for s in 1:numVars) + # # Y[mu] = @variable(jumpModel, [1:size(P,1),1:size(P,1)], PSD) + # if size(P, 1) > 1 + # constraints[name] = @constraint(jumpModel, Y[mu] in PSDCone()) + # else + # constraints[name] = @constraint(jumpModel, Y[mu] .>= 0) + # end + + Y[mu] = get(replaceBlocks, mu) do + name = "Y$mu" + if n > 1 + v = @variable(jumpModel, [1:n, 1:n], Symmetric, base_name = name) + constraints[name] = @constraint(jumpModel, v in PSDCone()) + return v + else + v = @variable(jumpModel, [1:1, 1:1], Symmetric, base_name = name) + constraints[name] = @constraint(jumpModel, v .>= 0) + return v + end + end end else - # numVars = maximum(P) - # y = @variable(jumpModel, [1:numVars], base_name = name) - # AT = typeof(1 * y[1]) - - # Y[mu] = zeros(AT, t, t) - # # Y[mu] = zeros(AT, size(P)) - # for s in 1:numVars - # Y[mu][P .== s] .+= y[s] + y = @variable(jumpModel, [1:1,1:1], base_name = mu) + Y[mu] = y + # i = parse(Int,mu[2:end]) + # for (G, c) in m.quotient[i].coeff + # graphCoefficients[G] = get(graphCoefficients, G, 0 * c * y) + c * y # end - - # Y[mu] = sum((P .== s) * y[s] for s in 1:numVars) - # # Y[mu] = @variable(jumpModel, [1:size(P,1),1:size(P,1)], PSD) - # if size(P, 1) > 1 - # constraints[name] = @constraint(jumpModel, Y[mu] in PSDCone()) - # else - # constraints[name] = @constraint(jumpModel, Y[mu] .>= 0) - # end - - Y[mu] = get(replaceBlocks, mu) do - name = "Y$mu" - if n > 1 - v = @variable(jumpModel, [1:n, 1:n], Symmetric, base_name = name) - constraints[name] = @constraint(jumpModel, v in PSDCone()) - return v - else - v = @variable(jumpModel, [1:1, 1:1], Symmetric, base_name = name) - constraints[name] = @constraint(jumpModel, v .>= 0) - return v - end - end end end @@ -507,14 +538,6 @@ function buildJuMPModel(m::RazborovModel, replaceBlocks=Dict(), jumpModel=Model( graphCoefficients[G] = eG end - for (i, F) in enumerate(m.quotient) - y = @variable(jumpModel, base_name = "quotient$i") - Y["Q$i"] = y - for (G, c) in F.coeff - graphCoefficients[G] = get(graphCoefficients, G, 0 * c * y) + c * y - end - end - return (model=jumpModel, variables=graphCoefficients, blocks=Y, constraints=constraints) end @@ -722,23 +745,28 @@ function verifySOS(m::RazborovModel, sol::Dict; io::IO=stdout) ) end - println(io, "Quotient:") - for (i, c) in enumerate(m.quotient) - if haskey(sol, "Q$i") && !iszero(sol["Q$i"]) - print(io, "$(sol["Q$i"]) ⋅ ($c) ") - println(io, "=$(Rational{BigInt}(sol["Q$i"])*c)= 0") + if length(m.quotient) > 0 + println(io, "Quotient:") + for (i, c) in enumerate(m.quotient) + if haskey(sol, "Q$i") && !iszero(sol["Q$i"]) + print(io, "$(sol["Q$i"]) ⋅ ($c) ") + println(io, "=$(Rational{BigInt}(sol["Q$i"])*c)= 0") + end end + println(io, "Quotient sum:") + println( + io, + "$(sum(enumerate(m.quotient)) do (i, F) + Rational{BigInt}(get(sol, "Q$i", 0 // 1)) * F + end)=0", + ) end - println(io, "Quotient sum:") - println( - io, - "$(sum(enumerate(m.quotient)) do (i, F) - Rational{BigInt}(get(sol, "Q$i", 0 // 1)) * F -end)=0", - ) end res = sum(m.sdpData) do (G, B) + if mu isa String + return BigInt(0)//1 + end sum(B) do (mu, b) if haskey(sol, mu) psd = sol[mu] isa Matrix ? sol[mu] : sol[mu].psd diff --git a/src/utils/Nauty.jl b/src/utils/Nauty.jl index f49ee3f..ad9093c 100644 --- a/src/utils/Nauty.jl +++ b/src/utils/Nauty.jl @@ -81,12 +81,12 @@ function refine!(coloring::Vector{Int}, F, v::Int)::Vector{UInt} end # coloring[coloring .> cell] .+= numNewCells - for i in length(alpha):-1:(cell + numNewCells + 1) - alpha[i] = alpha[i - numNewCells] + for i in length(alpha):-1:(cell+numNewCells+1) + alpha[i] = alpha[i-numNewCells] end # @views alpha[(cell + numNewCells + 1):end] .= alpha[(cell + 1):(end - numNewCells)] - alpha[(cell + 1):(cell + numNewCells)] .= true + alpha[(cell+1):(cell+numNewCells)] .= true newCellsCol::Vector{Pair{UInt,Int}} = collect(newCells) sort!(newCellsCol; by=x -> x[1]) @@ -128,7 +128,7 @@ function refine!(coloring::Vector{Int}, F, v::Int)::Vector{UInt} end # Remove biggest new cell - alpha[maxCellInd + cell - 1] = false + alpha[maxCellInd+cell-1] = false end cell += length(newCellsCol) @@ -142,7 +142,6 @@ function refine!(coloring::Vector{Int}, F, v::Int)::Vector{UInt} # append permuted Flag # return UInt[nodeInvariant, hash(permute(F, coloring))] # return nodeInvariant, hash(permute(F, coloring)) - return UInt[nodeInvariant, hash(permute(F, coloring))] # return hash(nodeInvariant, hash(permute(F, coloring))) end @@ -150,7 +149,7 @@ function refine!(coloring::Vector{Int}, F, v::Int)::Vector{UInt} # return nodeInv end -function investigateNode(F,coloring::Vector{Int}, nodeInv::Vector{UInt},n,nInv1, nInvStar,autG,v1,vStar,curBranch, covered,prune)::Int +function investigateNode(F, coloring::Vector{Int}, nodeInv::Vector{UInt}, n, nInv1, nInvStar, autG, v1, vStar, curBranch, covered, prune)::Int # @info "investigate node at $coloring, $nodeInv" first = length(nInv1) == 0 @@ -175,7 +174,7 @@ function investigateNode(F,coloring::Vector{Int}, nodeInv::Vector{UInt},n,nInv1, # @info "Added onto stabilizer, new order is $(order(autG))" stabilizer!(autG, curBranch) for i in 1:length(curBranch) - H = stabilizer!(autG, curBranch[1:(i - 1)]) + H = stabilizer!(autG, curBranch[1:(i-1)]) @assert H !== nothing @@ -184,8 +183,9 @@ function investigateNode(F,coloring::Vector{Int}, nodeInv::Vector{UInt},n,nInv1, if prune && curBranch[i] in o # cut this branch multiple levels up difLevel = length(curBranch) - i - curBranch = curBranch[1:i] - covered = covered[1:i] + resize!(curBranch, i) + resize!(covered, i) + # covered = covered[1:i] # @info "Pruning $difLevel levels from level $(length(curBranch)) via automorphism." return difLevel end @@ -236,15 +236,15 @@ function investigateNode(F,coloring::Vector{Int}, nodeInv::Vector{UInt},n,nInv1, # length(nodeInv) <= length(nInv1) @views if prune && - !first && - !( - nodeInv[1:min(length(nodeInv), length(nInv1))] == - nInv1[1:min(length(nodeInv), length(nInv1))] - ) && - !( - nodeInv[1:min(length(nodeInv), length(nInvStar))] >= - nInvStar[1:min(length(nodeInv), length(nInvStar))] - ) + !first && + !( + nodeInv[1:min(length(nodeInv), length(nInv1))] == + nInv1[1:min(length(nodeInv), length(nInv1))] + ) && + !( + nodeInv[1:min(length(nodeInv), length(nInvStar))] >= + nInvStar[1:min(length(nodeInv), length(nInvStar))] + ) return 0 end firstBigCell = Int[] @@ -265,8 +265,8 @@ function investigateNode(F,coloring::Vector{Int}, nodeInv::Vector{UInt},n,nInv1, newC::Vector{Int} = copy(coloring) newNodeInv::Vector{UInt} = refine!(newC, F, i) push!(curBranch, i) - catNodeInv::Vector{UInt} = vcat(nodeInv, newNodeInv) - t = investigateNode(F,newC, catNodeInv,n,nInv1, nInvStar,autG,v1,vStar,curBranch, covered,prune) + catNodeInv::Vector{UInt} = vcat(nodeInv, newNodeInv) + t = investigateNode(F, newC, catNodeInv, n, nInv1, nInvStar, autG, v1, vStar, curBranch, covered, prune) if t > 0 return t - 1 @@ -346,7 +346,7 @@ function label(F::T; prune=true, removeIsolated=true) where {T} # alpha[1] = true refine!(col, F, 0) - investigateNode(F,col,UInt[],n,nInv1, nInvStar,autG,v1,vStar,curBranch, covered,prune) + investigateNode(F, col, UInt[], n, nInv1, nInvStar, autG, v1, vStar, curBranch, covered, prune) p = zeros(Int, n) p[vStar] .= 1:n @@ -374,7 +374,8 @@ function generateAll( for i in 1:maxVertices nextGraphs = T[] - pq = PriorityQueue{EdgeMarkedFlag{T,predicateType(T)},Vector{Int}}() + # pq = PriorityQueue{EdgeMarkedFlag{T,predicateType(T)},Vector{Int}}() + pq = PriorityQueue{EdgeMarkedFlag{T,predicateType(T)},Any}() for f in generatedGraphs[i] newF = permute(f, 1:i) @@ -399,7 +400,7 @@ function generateAll( continue end - fixed = allowMultiEdges(T) ? Vector{Int}[] : [collect(1:(i - 1))] + fixed = allowMultiEdges(T) ? Vector{Int}[] : [collect(1:(i-1))] uP::Vector{Vector{predicateType(T)}} = findUnknownPredicates(newF, fixed, maxPredicates) # @show uP uP2 = findUnknownGenerationPredicates(newF, fixed, maxPredicates) @@ -423,15 +424,17 @@ function generateAll( if allowMultiEdges(T) F = EdgeMarkedFlag{T}(F.F, FMarked.marked) end - cP = countEdges(F)[1:(end - 1)] + # @show F + cP = countEdges(F)[1:(end-1)] + # @show cP # @assert length(maxPredicates) == length(cP) if length(maxPredicates) == length(cP) && all(cP .<= maxPredicates) pq[F] = cP elseif all( - cP[1:(length(maxPredicates) - 1)] .<= maxPredicates[1:(end - 1)] - ) && - maxPredicates[end] isa Int && - sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] + cP[1:(length(maxPredicates)-1)] .<= maxPredicates[1:(end-1)] + ) && + maxPredicates[end] isa Union{Int, Vector{Int}} && + sum(sum.(cP[length(maxPredicates):end])) <= sum(maxPredicates[end]) pq[F] = cP end end @@ -448,8 +451,8 @@ function generateAll( cP = countEdges(FMarked.F) if length(maxPredicates) == length(cP) && all(cP .== maxPredicates) continue - elseif all(cP[1:(length(maxPredicates) - 1)] .== maxPredicates[1:(end - 1)]) && - sum(cP[length(maxPredicates):end]) == maxPredicates[end] + elseif all(cP[1:(length(maxPredicates)-1)] .== maxPredicates[1:(end-1)]) && + sum(cP[length(maxPredicates):end]) == maxPredicates[end] continue end # if !all(countEdges(FMarked.F) .== maxPredicates) @@ -471,7 +474,7 @@ function generateAll( if !withProperty(F.F) continue end - cP = countEdges(F)[1:(end - 1)] + cP = countEdges(F)[1:(end-1)] # if all(cP .<= maxPredicates) # pq[F] = cP # end @@ -480,8 +483,8 @@ function generateAll( if length(maxPredicates) == length(cP) && all(cP .<= maxPredicates) pq[F] = cP elseif maxPredicates[end] isa Int && - all(cP[1:(length(maxPredicates) - 1)] .<= maxPredicates[1:(end - 1)]) && - sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] + all(cP[1:(length(maxPredicates)-1)] .<= maxPredicates[1:(end-1)]) && + sum(sum.(cP[length(maxPredicates):end])) <= maxPredicates[end] pq[F] = cP end end diff --git a/src/utils/PermutationModuleQuotients.jl b/src/utils/PermutationModuleQuotients.jl index b767264..10a7d0b 100644 --- a/src/utils/PermutationModuleQuotients.jl +++ b/src/utils/PermutationModuleQuotients.jl @@ -399,14 +399,16 @@ function symPolytabloidProduct( # factor out the size of the column stabilizer fact = prod([factorial(t) for t in conj(t1.part).part]) - + # fact = prod([factorial(t) for t in conj(t1.part).part])*prod([factorial(t) for t in conj(t2.part).part]) + # @show fact for j in (m - 1):-1:1 # First loop depends on order, second does not matter for s in (j + 1):m usj = u(s, j) - + fact *= factorial(usj) - + # @show fact + for k in 1:usj res2 = Dict{Matrix{T},T}() for (B, c) in res @@ -419,14 +421,14 @@ function symPolytabloidProduct( end end end - + for j in (m - 1):-1:1 # First loop depends on order, second does not matter for s in (j + 1):m rsj = r(s, j) fact *= factorial(rsj) - + for k in 1:rsj res2 = Dict{Matrix{T},T}() for (B, c) in res diff --git a/test/src/nauty.jl b/test/src/nauty.jl index 6059a3b..ccf2be5 100644 --- a/test/src/nauty.jl +++ b/test/src/nauty.jl @@ -36,14 +36,14 @@ end @time Gs = generateAll(Graph, 6, 1000) @time Gs = generateAll(Graph, 7, 1000) # 3.2s @profview Gs = generateAll(Graph, 7, 1000) recur=:flat -@profview Gs = generateAll(Graph, 7, 1000) format=:flat +# @profview Gs = generateAll(Graph, 7, 1000) format=:flat @profview_allocs Gs = generateAll(Graph, 7, 1000) recur=:flat @time Gs = generateAll(Graph, 8, 1000) # 57s -@time Gs = generateAll(Graph, 9, 1000) # 57s +# @time Gs = generateAll(Graph, 9, 1000) # 531s, after some optimization # Profile.clear() # @profile Gs = generateAll(Graph, 7, 1000)