Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

orthogonalize not defined #359

Open
jeremypparker opened this issue Jun 11, 2024 · 4 comments
Open

orthogonalize not defined #359

jeremypparker opened this issue Jun 11, 2024 · 4 comments

Comments

@jeremypparker
Copy link

jeremypparker commented Jun 11, 2024

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

@odow
Copy link
Member

odow commented Jun 11, 2024

if !all(is_orthogonal, S)
R = orthogonalize(R)

@blegat
Copy link
Member

blegat commented Jun 12, 2024

Thanks for letting us know about the issue, seems like it was removed in #243

@blegat
Copy link
Member

blegat commented Jun 12, 2024

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

@jeremypparker
Copy link
Author

jeremypparker commented Jun 12, 2024

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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Development

No branches or pull requests

3 participants