jump-dev / SumOfSquares.jl

Sum of Squares Programming for Julia
Other
114 stars 24 forks source link

orthogonalize not defined #359

Open jeremypparker opened 4 weeks ago

jeremypparker commented 4 weeks ago

Line 167 of wedderburn.jl calls a function orthogonalize, which as far as I can tell is not defined anywhere

odow commented 4 weeks ago

https://github.com/jump-dev/SumOfSquares.jl/blob/a668a7cc489002f082a8de252c09f761da84836e/src/Certificate/Symmetry/wedderburn.jl#L166-L167

blegat commented 3 weeks ago

Thanks for letting us know about the issue, seems like it was removed in https://github.com/jump-dev/SumOfSquares.jl/pull/243/

blegat commented 3 weeks ago

@jeremypparker do you have a reproducible example that hits that line of code ?

jeremypparker commented 3 weeks ago

Something like the following. I define a cyclic group which acts via a rotation matrix. when this group is order 6 I hit the line in question.

using SumOfSquares
using DynamicPolynomials
using LinearAlgebra
using COSMO
using GroupsCore

model = SOSModel(COSMO.Optimizer)

@polyvar z[1:3]

struct MatGrpElement{T} <: GroupElement
    num::Int
    order::Int
end
Base.:(==)(a::MatGrpElement, b::MatGrpElement) = isapprox(a.num, b.num, atol=1e-4, rtol=1e-4)
Base.inv(el::MatGrpElement{T}) where T = MatGrpElement{T}(mod(-el.num,el.order),el.order)

function Base.:*(a::MatGrpElement{T}, b::MatGrpElement{T}) where T
    return MatGrpElement{T}(mod(a.num+b.num,a.order),a.order)
end
Base.:^(el::MatGrpElement{T}, k::Integer) where T = MatGrpElement{T}(mod(el.num*k,el.order),el.order)

Base.conj(a::MatGrpElement, b::MatGrpElement) = inv(b) * a * b
Base.:^(a::MatGrpElement, b::MatGrpElement) = conj(a, b)

import PermutationGroups
function PermutationGroups.order(el::MatGrpElement)
    for i=1:100
        if el^i == one(el)
            return i
        end
    end
    throw("Could not compute order")
end

struct MatrixGroup{T} <: Group
    order::Int
end
Base.eltype(::MatrixGroup{T}) where T = MatGrpElement{T}
Base.one(c::MatrixGroup{T}) where T = MatGrpElement{T}(0,c.order)
Base.one(c::MatGrpElement{T}) where T = MatGrpElement{T}(0,c.order)
PermutationGroups.gens(c::MatrixGroup{T}) where T = [MatGrpElement{T}(1,c.order)]
PermutationGroups.order(::Type{T}, c::MatrixGroup) where {T} = convert(T, c.order)
function Base.iterate(c::MatrixGroup{T}, prev::MatGrpElement{T}=MatGrpElement{T}(-1,c.order)) where T
    num = prev.num + 1
    if num >= c.order
        return nothing
    else
        result = MatGrpElement{T}(num,c.order)
        return result,result
    end
end

struct MatrixAction <: Symmetry.OnMonomials
    variables
end
Symmetry.SymbolicWedderburn.coeff_type(::MatrixAction) = Float64
function Symmetry.SymbolicWedderburn.action(a::MatrixAction, el::MatGrpElement, mono::AbstractMonomial)
    th = el.num*2pi/el.order
    matrix = [cos(th) -sin(th) 0;
              sin(th) cos(th) 0;
              0 0 1]

    return subs(mono, a.variables=>matrix*a.variables)
end

symmetries = MatrixGroup{Int}(6)
action = MatrixAction(z)  
pattern = Symmetry.Pattern(symmetries, action)

basisV = monomials(z,0:2)
V = sum(@variable(model,[1:length(basisV)]).*basisV)

con_ref=@constraint(model, V*(1-dot(z,z)) >= 0, symmetry=pattern)

optimize!(model)