From 39f0b86af2f014ba961e3d1a033c126ee2601bae Mon Sep 17 00:00:00 2001
From: kalmar <kalmar@amu.edu.pl>
Date: Mon, 19 Dec 2016 15:44:52 +0100
Subject: [PATCH] Initial (working) code

---
 GroupAlgebras.jl     | 114 +++++++++++++++++++++++++++++++++++++++++
 Matrix_Constraints.g | 117 +++++++++++++++++++++++++++++++++++++++++++
 property(T).jl       |  70 ++++++++++++++++++++++++++
 3 files changed, 301 insertions(+)
 create mode 100644 GroupAlgebras.jl
 create mode 100644 Matrix_Constraints.g
 create mode 100644 property(T).jl

diff --git a/GroupAlgebras.jl b/GroupAlgebras.jl
new file mode 100644
index 0000000..cff32e1
--- /dev/null
+++ b/GroupAlgebras.jl
@@ -0,0 +1,114 @@
+module GroupAlgebras
+
+import Base: convert, show, isequal, ==
+import Base: +, -, *, //
+import Base: size, length, norm
+
+export GroupAlgebraElement
+
+
+immutable GroupAlgebraElement{T<:Number}
+    coordinates::Vector{T}
+    product_matrix::Array{Int,2}
+    # basis::Array{Any,1}
+
+    function GroupAlgebraElement(coordinates::Vector{T},
+        product_matrix::Array{Int,2})
+
+        length(coordinates) == size(product_matrix,1) ||
+            throw(ArgumentError("Matrix has to represent products of basis
+                                elements"))
+        size(product_matrix, 1) == size(product_matrix, 2) ||
+            throw(ArgumentError("Product matrix has to be square"))
+        # length(coordinates) == length(basis) || throw(ArgumentError("Coordinates must be given in the given basis"))
+        # new(coordinates, product_matrix, basis)
+        new(coordinates, product_matrix)
+    end
+end
+
+# GroupAlgebraElement(c,pm,b) = GroupAlgebraElement(c,pm)
+GroupAlgebraElement{T}(c::Vector{T},pm) = GroupAlgebraElement{T}(c,pm)
+
+convert{T<:Number}(::Type{T}, X::GroupAlgebraElement) =
+    GroupAlgebraElement(convert(Vector{T}, X.coordinates), X.product_matrix)
+
+show{T}(io::IO, X::GroupAlgebraElement{T}) = print(io,
+    "Element of Group Algebra over ", T, "\n", X.coordinates)
+
+
+function isequal{T, S}(X::GroupAlgebraElement{T}, Y::GroupAlgebraElement{S})
+    if T != S
+        warn("Comparing elements with different coefficients Rings!")
+    end
+    X.product_matrix == Y.product_matrix || return false
+    X.coordinates == Y.coordinates || return false
+    return true
+end
+
+(==)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = isequal(X,Y)
+
+function add{T<:Number}(X::GroupAlgebraElement{T}, Y::GroupAlgebraElement{T})
+    X.product_matrix == Y.product_matrix || throw(DomainError(
+    "Elements don't seem to belong to the same Group Algebra!"))
+    return GroupAlgebraElement(X.coordinates+Y.coordinates, X.product_matrix)
+end
+
+function add{T<:Number, S<:Number}(X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{S})
+    warn("Adding elements with different base rings!")
+    return GroupAlgebraElement(+(promote(X.coordinates, Y.coordinates)...),
+    X.product_matrix)
+end
+
+(+)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = add(X,Y)
+(-)(X::GroupAlgebraElement) = GroupAlgebraElement(-X.coordinates, X.product_matrix)
+(-)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = add(X,-Y)
+
+function group_star_multiplication{T<:Number}(X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{T})
+    X.product_matrix == Y.product_matrix || DomainError(
+    "Elements don't seem to belong to the same Group Algebra!")
+
+    result = zeros(X.coordinates)
+    for (i,x) in enumerate(X.coordinates), (j,y) in enumerate(Y.coordinates)
+        index = X.product_matrix[i,j]
+        if index != 0
+            result[index]+= x*y
+        end
+    end
+    return GroupAlgebraElement(result, X.product_matrix)
+end
+
+function group_star_multiplication{T<:Number, S<:Number}(
+    X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{S})
+    S == T || warn("Multiplying elements with different base rings!")
+    return group_star_multiplication(promote(X,Y)...)
+end
+
+(*){T<:Number, S<:Number}(X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{S}) = group_star_multiplication(X,Y);
+
+(*){T<:Number}(a::T, X::GroupAlgebraElement{T}) = GroupAlgebraElement(
+    a*X.coordinates, X.product_matrix)
+
+function scalar_multiplication{T<:Number, S<:Number}(a::T,
+    X::GroupAlgebraElement{S})
+    if  T!=S
+        warn("Scalars and coefficients ring are not the same! Trying to promote...")
+    end
+    return GroupAlgebraElement(a*X.coordinates, X.product_matrix)
+end
+(*){T<:Number}(a::T,X::GroupAlgebraElement) = scalar_multiplication(a, X)
+
+//{T<:Rational, S<:Rational}(X::GroupAlgebraElement{T}, a::S) =
+    GroupAlgebraElement(X.coordinates//a, X.product_matrix)
+
+//{T<:Rational, S<:Integer}(X::GroupAlgebraElement{T}, a::S) =
+    X//convert(T,a)
+
+length(X::GroupAlgebraElement) = length(X.coordinates)
+size(X::GroupAlgebraElement) = size(X.coordinates)
+norm(X::GroupAlgebraElement, p=2) = norm(X.coordinates, p)
+
+end
diff --git a/Matrix_Constraints.g b/Matrix_Constraints.g
new file mode 100644
index 0000000..effcd5e
--- /dev/null
+++ b/Matrix_Constraints.g
@@ -0,0 +1,117 @@
+Symmetrise := function(elts)
+  return Unique(Concatenation(elts, List(elts, Inverse)));
+end;
+
+MYAllProducts := function(elts1, elts2)
+  local products, elt;
+  products := [];
+  for elt in elts1 do
+    products := Concatenation(products, elt*elts2);
+  od;
+  return products;
+end;
+
+Products := function(elts, n)
+  local products, i;
+  if n<=0 then
+      return [ ];
+  elif n = 1 then
+      return elts;
+  else
+    products := elts;
+    for i in [2..n] do
+      products := MYAllProducts(elts, products);
+    od;
+    return products;
+  fi;
+end;
+
+Laplacian := function(G, generating_set)
+  local QG, emb, result, S, g, elt;
+  QG := GroupRing(Rationals, G);;
+  emb := Embedding(G,QG);;
+
+  S := generating_set;
+
+  result := Length(S)*One(QG);
+  for g in S do
+      result := result - g^emb;
+  od;
+  return result;
+end;
+
+Vectorise := function(elt, basis)
+  local result, l, i, g, coeff, axis;
+  Assert(0, IsSupportedOn(basis, elt),
+    "AssertionError: Element of interest is not supported on the basis!");
+  result := List(0*[1..Length(basis)]);
+
+  l := CoefficientsAndMagmaElements(elt);
+  for i in [1..Length(l)/2] do
+    g := l[2*i-1];
+    coeff := l[2*i];
+    axis := Position(basis, g);
+    result[axis] := result[axis] + coeff;
+  od;
+  return result;
+end;
+
+Constraints := function(basis)
+  local result, i, j, pos;
+  result := [];
+  for i in [1..Length(basis)] do
+    Add(result,[]);
+  od;
+  for i in [1..Length(basis)] do
+     for j in [1..Length(basis)] do
+       pos := Position(basis, Inverse(basis[i])*basis[j]);
+       if not pos = fail then
+          Add(result[pos], [i,j]);
+       fi;
+     od;
+  od;
+  return result;
+end;
+
+USupport := function(x)
+    return Unique(Support(x));
+end;
+
+IsSupportedOn := function(basis, elt)
+    local elt_supp, x;
+    elt_supp := USupport(elt);
+    for x in elt_supp do
+        if not x in basis then
+            return x;
+        fi;
+    od;
+    return true;
+end;
+
+SDPGenerateAll := function(G, S, basis, name)
+  local QG, emb, delta, delta_sq, delta_vec, delta_sq_vec, product_constr;
+  QG := GroupRing(Rationals, G);;
+  emb := Embedding(G,QG);;
+
+  delta := Laplacian(G, S);;
+  delta_sq := delta^2;;
+  if not IsSupportedOn(basis, delta_sq) then
+    #   Print("delta_sq is not supported on basis\n");
+      return fail;
+  else
+      PrintTo(Concatenation("./basis.", name), basis);
+      Print("Written basis to ", Concatenation("./basis.", name), "\n");
+      delta_vec := Vectorise(delta, basis);;
+      PrintTo(Concatenation("./delta.", name), delta_vec);
+      Print("Written delta to ", Concatenation("./delta.", name), "\n");
+      delta_sq_vec := Vectorise(delta_sq, basis);;
+      PrintTo(Concatenation("./delta_sq.", name), delta_sq_vec);
+      Print("Written delta_sq to ", Concatenation("./delta_sq.", name), "\n");
+
+      product_constr := Constraints(basis);;
+      PrintTo(Concatenation("./constraints.", name), product_constr);
+      Print("Written Matrix Constraints to ", Concatenation("./Constraints.", name), "\n");
+      return "Done!";
+  fi;
+
+end;;
diff --git a/property(T).jl b/property(T).jl
new file mode 100644
index 0000000..ee57bde
--- /dev/null
+++ b/property(T).jl
@@ -0,0 +1,70 @@
+using JuMP
+
+
+function read_GAP_raw_list(filename::String)
+    return eval(parse(String(read(filename))))
+end
+
+function create_product_matrix(matrix_constraints)
+    l = length(matrix_constraints)
+    product_matrix = zeros(Int, (l, l))
+    for (index, pairs) in enumerate(matrix_constraints)
+        for (i,j) in pairs
+            product_matrix[i,j] = index
+        end
+    end
+    return product_matrix
+end
+
+function create_sparse_product_matrix(matrix_constraints)
+    row_indices = Int[]
+    column_indices = Int[]
+    values = Int[]
+    for (index, pairs) in enumerate(matrix_constraints)
+        for (i,j) in pairs
+            push!(row_indices, i)
+            push!(column_indices, j)
+            push!(values, index)
+        end
+    end
+    sparse_product_matrix = sparse(row_indices, column_indices, values)
+    return sparse_product_matrix
+end
+
+function create_SDP_problem(matrix_constraints,
+                            Δ²::GroupAlgebraElement, Δ::GroupAlgebraElement)
+    N = length(Δ)
+    @assert length(Δ) == length(Δ²)
+    @assert length(Δ) == length(matrix_constraints)
+    m = Model();
+    @variable(m, A[1:N, 1:N], SDP)
+    @SDconstraint(m, A >= zeros(size(A)))
+    @variable(m, κ >= 0.0)
+    @objective(m, Max, κ)
+
+    for (pairs, δ², δ) in zip(matrix_constraints, Δ².coordinates, Δ.coordinates)
+        @constraint(m, sum(A[i,j] for (i,j) in pairs) == δ² - κ*δ)
+    end
+    return m
+end
+
+function resulting_SOS{T<:Number, S<:Number}(sqrt_matrix::Array{T,2}, elt::GroupAlgebraElement{S})
+    zzz = zeros(T, size(sqrt_matrix)[1])
+    result::GroupAlgebraElement{T} = GroupAlgebraElement(zzz, elt.product_matrix)
+    for i in 1:length(result)
+        new_base = GroupAlgebraElement(sqrt_matrix[:,i], elt.product_matrix)
+        result += new_base*new_base
+    end
+    return result
+end
+
+function correct_to_augmentation_ideal{T<:Rational}(sqrt_matrix::Array{T,2})
+    sqrt_corrected = similar(sqrt_matrix)
+    l = size(sqrt_matrix,2)
+    for i in 1:l
+        col = view(sqrt_matrix,:,i)
+        sqrt_corrected[:,i] = col - sum(col)//l
+#         @assert sum(sqrt_corrected[:,i]) == 0
+    end
+    return sqrt_corrected
+end