/*  the symmetry group of frames

   Markus Grassl <markus.grassl@mpl.mpg.de>
   Shayne Waldron <swal089@math.auckland.ac.nz  >
   2017-03-01


   brief list of functions:
   - CanonicalGramian(Matrix)
   - MinimalCycleBasis(Graph)
   - FrameSymmetryMatrix(V)
     input: a (d x n) matrix V, the vectors are the columns of V
   - FrameSymmetry(L)
     input: a sequence of vectors in a d-dimensional vector space
   - FrameSymmetry(G)
     input: a canonical Gram matrix

   - IsIsomorphicFrame(G1,G2);
     IsIsomorphicFrameMatrix(V1,V2);
     test whether two frames are isomorphic

   - permutation_to_matrix(V,perm)
     input: matrix V of the frame and a  permutation symmetry
     output: matrices M and D such that D*V eq V*perm*D

   - permutation_to_matrix_grp(V,sym)
     input: matrix V of the frame and its permutation symmetry

*/

declare verbose FrameSymmetry,2;


/*

  for verbose level 2, output information on the computation of the cycle basis

/*
******************************************************************************
  Calculate the canonical Gramian P of the frame (v_j) given by the columns of
  the d x n matrix V=[v1,...vn]. If V is an orthogonal projection this should
  return V.

  The base ring of V must be a field on which complex conjugation is defined.
******************************************************************************
*/

// If complex conjugation is defined on the base field, then we can take
// the HermitianTranspose
HermTranspose := function(B)
  T:=Transpose(B);
  T:=Parent(T)![ComplexConjugate(z):z in Eltseq(T)];
  return T;
end function;

// Define the standard inner product on vectors (linear in the first variable)
// This might need to be refined for other fields
ip := function(u,v)
  if Type(BaseRing(u)) in {FldRat,FldCyc,FldCom,FldRe} then
    return InnerProduct(u,v);
  else
    error("The inner product is not defined over this field");
  end if;
end function;

// For F=ComplexField(m), NullspaceOfTranspose(V) fails.
// Better to do the calculations in matlab:
// [d n]=size(V), A=null(V); P=eye(n)-A*A'

myNullspaceOfTranspose:=function(V);
// a faster way to compute the nullpsace of the transpose of V
// WARNING: this is numerically unstable
  V1:=Transpose(EchelonForm(V));
  ind_pivot:=[];
  row:=1;
  for col:=1 to Ncols(V1) do
    while row le Nrows(V1) and Support(V1[row]) ne {col} do
      row+:=1;
    end while;
    if row le Nrows(V1) then
      Append(~ind_pivot,row);
    end if;
  end for;
  M0:=Matrix([V1[i]:i in [1..Nrows(V1)]|not i in ind_pivot]);
  M:=HorizontalJoin(M0,-One(MatrixRing(CoefficientRing(M0),Nrows(M0))));
  return [M[i]:i in [1..Nrows(M)]];
end function;

intrinsic CanonicalGramian(V::Mtrx)->Mtrx
{}
  F:=BaseRing(V);
  n:=Ncols(V);
// Find the kernel of V (with usual matrix multiplication)

  vprintf FrameSymmetry: "  computing the kernel ";
  vtime FrameSymmetry:
//  B:=myNullspaceOfTranspose(V);
  dep:=NullspaceOfTranspose(V);
  B:=BasisMatrix(dep);
// Apply Gram Schmidt without normalisation to the basis B
  L:=B;
  vprintf FrameSymmetry: "  Gram-Schmidt orthogonalisation ";
  vtime FrameSymmetry:
  for j:=1 to Nrows(B) do
    v:=B[j];
    prj:=0*v;
    for k in [1..j-1] do
      prj+:=ip(v,L[k])/ip(L[k],L[k])*L[k];
    end for;
    v-:=prj;
    L[j]:=v;
  end for;
// To define the projection operator we need to use the Hermitian Transpose
  P:=Zero(MatrixRing(F,n));
  vprintf FrameSymmetry: "  computing the projection ";
  vtime FrameSymmetry:
  for j:=1 to Nrows(L) do
    v:=Transpose(Matrix([L[j]]));
    p0:=v*HermTranspose(v);
    P+:=1/Trace(p0)*p0;
  end for;
  I:=Identity(Parent(P));
  P:=I-P;
  return(P);
end intrinsic;

// Calculate the frame graph of a frame given by the columns of a
// d x n matrix V=[v1,...vn].
intrinsic FrameGraph(V:Mtrx)->Grph
{Calculate the grame grpah a a frame given byt the columns of the matrix V}
  F:=BaseRing(V);
  n:=Ncols(V);
  P:=CanonicalGramian(V);
  A:=Matrix([[P[j,k] eq 0 or j eq k select 0 else 1:k in [1..n]]:j in [1..n]]);
  return Graph<n|A>;
end intrinsic;


intrinsic MinimalCycleBasis(G::Grph:UseAut:=false)-> SeqEnum,SetIndx
{Computes a minimal cycle basis for the graph G.
 The basis is returned as a sequence of binary vectors.
 The second return value is the indexed set of edges of G.
 When "UseAut" is set to true, the orbit of cycles under the autormorphism group fo the graph is used.
}

/* based on
  Horton, J. D.,
  A polynomial-time algorithm to find the shortest cycle basis of a graph,
  SIAM Journal on Computing, vol. 16, no 2, April 1987, pp. 358-366

*/

  //optimization: perform the search for each component separately
  components:=ConnectedComponents(G);

  V:=Vertices(G);
  E:=Edges(G);
  if UseAut then
    vprintf FrameSymmetry,2: "computing the symmetry group of the graph ";
    vtime FrameSymmetry,2:
    aut,vertex_Gset,edge_Gset:=AutomorphismGroup(G);
  end if;
  V0:=KSpace(GF(2),#E);
  candidates:=[V0|];

  //for every vertex v and every edges e=(x,y),
  //  we look for the shortest paths P(v,x) and P(v,y) from v to x and y, resp.
  //  then e+P(v,x)+P(v,y) forms a cycle
  //  note that some edges might occur twice, so that they cancel in the cycle
  //  Horton calls this "degenerate cases"

  // possibly use the orbit under the automorphism group of the graph
  ty:=Cputime();
  dim:=0;
  for C in components do
    G1:=sub<G|C>;
    V1:=Vertices(G1);
    E1:=Edges(G1);
    dim+:=#E1-#V1+1;
    for v in V1 do
      L_path:=ShortestPaths(v);
      for e in E1 do
        pts:=EndVertices(e);
        if not v in pts then
          path:=[Set(L_path[Position(V1,p)]):p in pts];
          //remove the common edges
          path:=path[1] sdiff path[2];
          if not e in path then
            Include(~path,e);
            if UseAut then
              vprintf FrameSymmetry,2: "computing the orbit of the cycle ";
              vtime FrameSymmetry,2:
              orbit:=Orbit(aut,edge_Gset,path);
            else
              orbit:={path};
            end if;
            for p in orbit do
              incidence:=V0!0;
              for x in p do
                incidence[Position(E,E!x)]+:=1;
              end for;
              Append(~candidates,incidence);
            end for;
          end if;
        end if;
      end for;
      if #candidates eq dim and not exists{v:v in candidates|Weight(v) ne 3} then
        vprintf FrameSymmetry,2: "  weights: %o, 3-cycle-space: %o, target: %o (%o s)\n",
           {* Weight(x):x in candidates *},dim,dim,Cputime(ty);
	break v;
      else
        V3:=sub<V0|{v:v in candidates|Weight(v) eq 3}>;
        vprintf FrameSymmetry,2: "  weights: %o, 3-cycle-space: %o, target: %o (%o s)\n",
           {* Weight(x):x in candidates *},Dimension(V3),dim,Cputime(ty);
        if Dimension(V3) eq dim then break v;end if;
      end if;
    end for;
  end for;
  if #{Weight(x):x in candidates} ne 1 then
    Sort(~candidates,func<a,b|Weight(a)-Weight(b)>);
  end if;
  if #candidates eq dim then
    basis:=candidates;
  else
    B:=SparseMatrix(GF(2),#E,#candidates);
    for i:=1 to #candidates do
      for j in Support(candidates[i]) do
  	B[j,i]:=1;
      end for;
    end for;
    Echelonize(~B);
    basis:=[V0|candidates[Rep(Support(B[i]))]:i in [1..dim]];
  end if;
  return basis,E;
end intrinsic;



EdgeCycleToSequence:=function(L);
// give a sequence or set of edges corresponding to a path, return an ordered set of vertices along the path
  L1:=[EndVertices(x):x in L];
  if #L1 eq 0 then
    return [];
  end if;
  if Set(Multiplicities({*x:x in a,a in L1*})) ne {2} then
    error("The sequence of edges is not a cycle");
  end if;
  // the initial vertex
  v1:=Rep(L1[1]);
  vertices:=[v1];
  v:=v1;
  repeat
    next_edge:=rep{x:x in L1|v in x};
    Exclude(~L1,next_edge);
    v:=Rep(Exclude(next_edge,v));
    Append(~vertices,v);
  until v eq v1;
  return vertices;
end function;

intrinsic AutomorphismGroupMatrix(M::Mtrx)->GrpPerm
{Computes the automorphism group of m x n matrix M as permutations on n+m points}
//see also package/Group/GrpPerm/matrix_aut.m

  m:=Nrows(M);
  n:=Ncols(M);

//  ty:=Cputime();
//  G1:=AutomorphismGroup(M);
//  ty1:=Cputime(ty);

  ty:=Cputime();

  p:=251;
  F:=GF(p);
  w:=PrimitiveElement(F);
  //values0:={* x:x in Eltseq(M) *};

  values:=SetToIndexedSet(&join{Set(Eltseq(M[i])):i in [1..Nrows(M)]});
  L_values:=[values[i..Minimum(#values,i+p-2)]:i in [1..#values by p-1]];

  if #values le 2 and m eq n then
  //for square matrices, restrict the initial group to act simultaneously on rows and columsn
  //  this allows to process matrices of the form (a-b)*I+b*J faster
  //  WARNING: for matrices with all entries equal, this is NOT the full symmetry group
  //  => treat those cases differently
    G0,embed:=DirectProduct(Sym(m),Sym(n));
    G0:=sub<Sym(n+m)|[&*[f(g):f in embed]:g in Generators(Sym(n))]>;
  else
    G0:=DirectProduct(Sym(m),Sym(n));
  end if;

  if not exists{g:g in Generators(G0)|M^g ne M} then
    return G0;
  end if;

  checked:=false;
  cnt:=1;

  vprintf FrameSymmetry: "    hashing the matrix: ";
  vtime FrameSymmetry:
//  M1:=Matrix(m,n,[F!Position(values,x):x in Eltseq(M)]);
  if m*n ge 2^30 then
    M1:=KMatrixSpace(F,m,n)!0;
    for i:=1 to m do
      M1[i]:=Vector([F!Hash(x):x in Eltseq(M[i])]);
     end for;
  else
    M1:=Matrix(m,n,[F!Hash(x):x in Eltseq(M)]);
  end if;
  
  vprintf FrameSymmetry: "    computing the matrix automorphism group: ";
  vtime FrameSymmetry:
  G0:=AutomorphismGroupFF(M1);
  
  for val in L_values do
    vprintf FrameSymmetry: "    testing the %o generators: ",Ngens(G0);
    vtime FrameSymmetry:
    if not exists{g:g in Generators(G0)|M^g ne M} then
      checked:=true;
      break;
    elif #G0 le 100 then
      //exhaustive search for small groups
      vprintf FrameSymmetry: " => exhaustive search in group of order %o: ",#G0;
      G1:=sub<G0|>;
      for g in G0 do
        if not g in G1 and M^g eq M then
          G1:=sub<G0|G1,g>;
          BSGS(G1);
        end if;
      end for;
      G0:=G1;
      checked:=true;
      break;
    end if;
    vprintf FrameSymmetry: "    iteration %o/%o ",cnt,#L_values;
    if m*n ge 2^30 then
      for i:=1 to m do
        M1[i]:=Vector([p eq 0 select 0 else w^p where p:=Position(val,x):x in Eltseq(M[i])]);
      end for;
    else
      M1:=Matrix(m,n,[p eq 0 select 0 else w^p where p:=Position(val,x):x in Eltseq(M)]);
    end if;
    vprintf FrameSymmetry: "(%o s)\n",Cputime(ty);
    vprintf FrameSymmetry: "    computing the matrix automorphism group: ";
    vtime FrameSymmetry:
    G1:=AutomorphismGroupFF(M1);
    vprintf FrameSymmetry: "    computing the intersection: ";
    vtime FrameSymmetry:
    G0:=G0 meet G1;
    vprintf FrameSymmetry: "    current order: %o\n",#G0;
    cnt+:=1;
  end for;

  if not checked and exists{g:g in Generators(G0)|M^g ne M} then
  // fallback in case the modular matrices do not determine the automorphism group
    if #G0 le 10^6 then
      //exhaustive search for medium-size groups
      vprintf FrameSymmetry: "    exhaustive search in group of order %o: ",#G0;
      done:={};
      G1:=sub<G0|>;
      while exists(g){g:g in Transversal(G0,G1)|not IsId(g) and M^g eq M} do
        G1:=sub<G0|G1,g>;
        BSGS(G1);
      end while;
      G0:=G1;
    else
      vprintf FrameSymmetry: "    FALLBACK: using original algorithm ";
      vtime FrameSymmetry:
      G0:=AutomorphismGroup(M);
    end if;
  end if;

  return G0;

end intrinsic;

intrinsic IsIsomorphicMatrices(M1::Mtrx,M2::Mtrx)->BoolElt,GrpPermElt
{Checks whether ther matrices M1 and M2 are related by a permutation of the rows and columns.}

  ty:=Cputime();

  m1:=Nrows(M1);
  n1:=Ncols(M1);

  m2:=Nrows(M2);
  n2:=Ncols(M2);

  if m1 ne m2 or n2 ne n2 or Parent(M1[1,1]) cmpne Parent(M2[1,1]) then
    return false,_;
  end if;   

  p:=251;
  F:=GF(p);
  w:=PrimitiveElement(F);
  values1:={* x:x in Eltseq(M1[i]),i in [1..m1] *};
  values2:={* x:x in Eltseq(M2[i]),i in [1..m2] *};
  if values1 ne values2 then
    return false,_;
  end if;   

  vprintf FrameSymmetry: "    hashing the first matrix: ";
  vtime FrameSymmetry:
  if m1*n1 ge 2^30 then
    M1a:=KMatrixSpace(F,m1,n1)!0;
    for i:=1 to m1 do
      if IsCoercible(F,M1[1,1]) then
        M1a[i]:=ChangeRing(M1[i],F);
      else	
        M1a[i]:=Vector([F!Hash(x):x in Eltseq(M1[i])]);
      end if;
    end for;
  else
    if IsCoercible(F,M1[1,1]) then
      M1a:=ChangeRing(M1,F);
    else
      M1a:=Matrix(m1,n1,[F!Hash(x):x in Eltseq(M1)]);
    end if;
  end if;

  vprintf FrameSymmetry: "    hashing the second matrix: ";
  vtime FrameSymmetry:
  if m2*n2 ge 2^30 then
    M2a:=KMatrixSpace(F,m2,n2)!0;
    for i:=1 to m2 do
      if IsCoercible(F,M2[1,1]) then
        M2a[i]:=ChangeRing(M2[i],F);
      else	
        M2a[i]:=Vector([F!Hash(x):x in Eltseq(M2[i])]);
      end if;
    end for;
  else
    if IsCoercible(F,M2[1,1]) then
      M2a:=ChangeRing(M2,F);
    else
      M2a:=Matrix(m2,n2,[F!Hash(x):x in Eltseq(M2)]);
    end if;
  end if;

  fl,perm:=IsIsomorphicFF(M1a,M2a);

  if not fl then
     return false,_;
  elif M1^perm eq M2 then
    vprintf FrameSymmetry: "    isomorphism of order %o\n",Order(perm);
    return true,perm;
  end if;


//  vprintf FrameSymmetry: "  computing the automorphism group of the first matrix: ";
//  vtime FrameSymmetry:
//  grp1:=AutomorphismGroupMatrix(M1a);
//
//  vprintf FrameSymmetry: "  computing the automorphism group of the second matrix: ";
//  vtime FrameSymmetry:
//  grp2:=AutomorphismGroupMatrix(M2a);

  vprintf FrameSymmetry: "  computing the automorphism group of both matrices:\n";
  vtime FrameSymmetry:
  grp12:=AutomorphismGroupMatrix(DirectSum(M1,M2));

//  vprintf FrameSymmetry: "    order: %o * %o * %o\n",#grp12 div (#grp1*#grp2),#grp1,#grp2;

  fl,tau:=IsConjugate(grp12,{1..n1} join {n1+n2+1..n1+n2+m1},{n1+1..n1+n2} join {n1+n2+m1+1..n1+n2+m1+m2});
  assert fl;
  perm:=Sym(n1+m1)!([i^tau-n1:i in [1..n1]] cat [i^tau-n2-m1: i in [n1+n2+1..n1+n2+m2]]);
  assert M1^perm eq M2;

  vprintf FrameSymmetry: "    isomorphism of order %o\n",Order(perm);
  
  if M1^perm eq M2 then
    return true,perm;
  else
//    vprintf FrameSymmetry: "\nWARNING: modular matrices isoorphic, but wrong permutation\n";
    printf "\nWARNING: modular matrices isoorphic, but wrong permutation\n";
  end if;   

  values:=SetToIndexedSet(Set(values1));
  L_values:=[values[i..Minimum(#values,i+p-2)]:i in [1..#values by p-1]];

  cnt:=1;
  for val in L_values do
    if m1*n1 ge 2^30 then
      M1a:=KMatrixSpace(F,m1,n1)!0;
      for i:=1 to m1 do
  	M1a[i]:=Vector([F!Hash(x):x in Eltseq(M1[i])]);
      end for;
    else
      M1a:=Matrix(m1,n1,[F!Hash(x):x in Eltseq(M1)]);
    end if;

    if m2*n2 ge 2^30 then
      M2a:=KMatrixSpace(F,m2,n2)!0;
      for i:=1 to m2 do
  	M2a[i]:=Vector([F!Hash(x):x in Eltseq(M2[i])]);
      end for;
    else
      M2a:=Matrix(m2,n2,[F!Hash(x):x in Eltseq(M2)]);
    end if;

    fl,perm:=IsIsomorphicFF(M1a,M2a);
    if not fl then
       return false,_;
    elif M1^perm eq M2 then
      return true,perm;
    else
      vprintf FrameSymmetry: "\nWARNING: modular matrices isoorphic at try %o/%o, but wrong permutation\n",cnt,#L_values;
    end if;   
    cnt+:=1;
  end for;

  // as fallback, use the original algorithm
  return IsIsomorphic(M1,M2);
  
end intrinsic;


intrinsic PathStabilizer(G::GrpPerm,S0::GSetIndx,S::Set)->GrpPerm
{Compute the stabilizer of the S in the group G}
  action:=hom<G->Sym(#S0)|[[Position(S0,x^G.i):x in S0]:i in [1..Ngens(G)]]>;
  G_action:=sub<Sym(#S0)|[action(G.i):i in [1..Ngens(G)]]>;
  inverse_action:=hom<G_action->G|[G.i: i in [1..Ngens(G)]]>;
  S1:={Position(S0,x):x in S};
  stab:=Stabilizer(G_action,S1);
  erg:=sub<G|[inverse_action(x):x in Generators(stab)]>;
  return erg;
end intrinsic;

intrinsic FrameSymmetryMatrix(V::Mtrx:lowerbound:=0,eps:=1e-10)->GrpPerm
{}
  P:=CanonicalGramian(V);
  return FrameSymmetry(P:lowerbound:=lowerbound,eps:=eps);
end intrinsic;

intrinsic FrameSymmetry(L::[]:lowerbound:=0,eps:=1e-10)->GrpPerm
{}
  V:=Transpose(Matrix([v:v in L]));
  vprintf FrameSymmetry: "computing the canonical Gramian\n";
  vtime FrameSymmetry:
  P:=CanonicalGramian(V);
  return FrameSymmetry(P:lowerbound:=lowerbound,eps:=eps);
end intrinsic;

intrinsic FrameSymmetry(G::Mtrx:
  lowerbound:=1,
  sort_blocks:=false,
  eps:=1e-10)->GrpPerm
{Compute the symmetry group of a frame given by a canonical Gram matrix as a permutation group.
 When the optional parameter "lowerbound" is set, the algorithm terminates when the size symmetry group is not greater than the lowerbound.
 The parameter "sort_blocks" determines whether to sort the the blocks by decreasing size.
 The optional parameter "eps" defines the precision for floating point input.
}
// the input G is a Gram matrix
// CHECKING to be added

  ty:=Cputime();

  //check whether we have exact or floating point arthmetics
  F:=CoefficientRing(G);
  exact:=not (Type(F) eq FldRe or Type(F) eq FldCom);
  if exact then
    round:=func<x|x>;
  else
    M:=F!Round(1/eps);
    round:=func<x|Round(M*x)>;
    // the following MIGHT give the anti-unitary symmetries
    // round:=func<x|Rep({Round(M*x),Round(M*ComplexConjugate(x))})>;
  end if;

  n:=Nrows(G);
  G_round:=Matrix(n,n,[round(x):x in Eltseq(G)]);
  A:=Matrix([[i eq j or G_round[i,j] eq 0 select 0 else 1:j in [1..n]]:i in [1..n]]);
  frame_graph:=Graph<n|A>;

  //determine the symmetry group of the 1-cycles and 2-cycles
  G_abs:=Matrix([[round(G[i,j]*G[j,i]):j in [1..n]]:i in [1..n]]);
  vprintf FrameSymmetry: "computing the automorphism group of the 1- and 2-cycles ";
  vtime FrameSymmetry:
  sym0:=AutomorphismGroupMatrix(G_abs);
  sym:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym0)]>;
  if IsSymmetric(sym) then
    vprintf FrameSymmetry: "  initial group is the full symmetric group Sym(%o)\n",n;
  elif #sym le 10^10 then
    vprintf FrameSymmetry: "  initial order: %o\n",#sym;
  else
    vprintf FrameSymmetry: "  initial order: 2^%o\n",Ilog2(#sym);
  end if;

  if #sym eq 1 then
    return sym;
  end if;

  current_order:=#sym;
  vertices:=Vertices(frame_graph);
  vprintf FrameSymmetry,2: "computing a minimal cycle basis:\n";
  vtime FrameSymmetry,2:
  cycle_basis,edges:=MinimalCycleBasis(frame_graph);
  dim:=#cycle_basis;

  EdgeCycleProduct:=function(L);
  // compute the m-product corresponding to the edge cycle L
    L_vertices:=EdgeCycleToSequence(L);
    indices:=[Position(vertices,v):v in L_vertices];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  VertexCycleProduct:=function(L);
  // compute the m-product corresponding to the vertex cycle L
    indices:=[Position(vertices,v):v in L];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  CycleProduct:=function(L);
  // compute the m-product corresponding to the cycle L given by indices
    //initialisation by the closing of the cycle
    if L[1] eq L[#L] then
      prod:=1;
    else
      prod:=G[L[#L],L[1]];
    end if;
    for i:=2 to #L do
      prod*:=G[L[i-1],L[i]];
    end for;
    return prod;
  end function;

  //initialise the space of the cycles processed so far
  C:=sub<Universe(cycle_basis)|>;
  C_gen:=SparseMatrix(GF(2),0,#edges);

  while #cycle_basis ne 0 and #sym gt lowerbound do
    v_cycle:=cycle_basis[1];
    if v_cycle in C then
      cycle_basis:=cycle_basis[2..#cycle_basis];
    else
      // convert the incidence vector of the cycle to a sequence of indices of vertices
      cycle:=EdgeCycleToSequence({edges[i]:i in Support(v_cycle)});
      cycle:=[Position(vertices,x):x in cycle];

      m:=#cycle-1;
      vprintf FrameSymmetry: "selected the %o-cycle %o (%o s)\n",m,cycle[1..m],Cputime(ty);
      vprintf FrameSymmetry: "  computing the orbit: ";
      vtime FrameSymmetry:
      if IsSymmetric(sym) then
        orbit:=GSet(sym,{x:x in Subsequences({1..n},m)|#Set(x) eq m});
      else
        orbit:=cycle[1..m]^sym;
      end if;
      vprintf FrameSymmetry:  "  orbit size: %o\n",#orbit;

      vprintf FrameSymmetry: "  computing the cycle-products: ";
      vtime FrameSymmetry:
      cycle_products:={<c,round(CycleProduct(c))>:c in orbit};

      values0:={* x[2]:x in cycle_products *};
      //sort the values by increasing frequencies
      if sort_blocks then
        values:=Setseq(Set(values0));
        Sort(~values,func<a,b|Multiplicity(values0,a)-Multiplicity(values0,b)>);
      else
        values:=Set(values0);
      end if;
      part0:={{x[1]:x in cycle_products|x[2] eq val}:val in values};
//      vprintf FrameSymmetry: "  %o different values (%o real) with multiplicities %o\n",#values,#{x:x in values|IsReal(x)},{*#x:x in part0*};
      vprintf FrameSymmetry: "  %o different values with multiplicities %o\n",#values,{*#x:x in part0*};


      // now look at the blocks of the orbit with a fixed value of the cycle-product
      for val in values do
        part:={x[1]:x in cycle_products|x[2] eq val};

        // convert the cycles to vectors in the cycle space
        for cyc0 in part do
          cyc:=Append(cyc0,cyc0[1]);
          j:=Nrows(C_gen)+1;
          for i:=2 to #cyc do
            edge:=edges!{vertices[cyc[i-1]],vertices[cyc[i]]};
            SetEntry(~C_gen,j,Position(edges,edge),1);
          end for;
        end for;

        if #values gt 1 then
//          part:={Rotate(x,i):i in [0..#x-1],x in part};
          vprintf FrameSymmetry: "  computing the stabilizer for %o points: ",#part;
          vtime FrameSymmetry:
          sym:=Stabilizer(sym,orbit,part);
        end if;
        vprintf FrameSymmetry: "  current order: %o\n",#sym;

        Echelonize(~C_gen);
        C:=Image(Matrix(C_gen));
        vprintf FrameSymmetry: "  dimemsion of the cycle sub-space: %o (of %o)\n",Nrows(C_gen),dim;
        if Nrows(C_gen) eq dim then break;end if;
        // when the order has changed, we use the smaller group to compute the orbits
        if #sym ne current_order then
          current_order:=#sym;
          break;
        end if;
      end for;
    end if;
  end while;

  return sym;
end intrinsic;

intrinsic CycleStabilizer_v0(L::Setq,n::RngIntElt)->GrpPerm
// first idea to compute the permutation stabilizer of a set of cycles
{Compute the stabilizer of a set of cycles given by a sequence of indices, e.g. [1,2,3].}
  m:=n+#L;
  A:=MatrixRing(Integers(),n+m)!0;
  j:=n+1;
  for cyc in L do
    for i:=2 to #cyc do
      A[cyc[i-1],cyc[i]]+:=1;
    end for;
// we close the loop in order to really have a cycle, maybe the same when computing the stabilizer directly?
    A[cyc[#cyc],cyc[1]]+:=1;
    for k in cyc do
      A[k,j]:=1;
    end for;
    j+:=1;
  end for;
  grp0:=AutomorphismGroupMatrix(A);
  grp0:=Stabilizer(grp0,{1..n});
  grp:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(grp0)]>;
  return grp;
end intrinsic;

intrinsic CyclePartitionStabilizer(part::Setq,n::RngIntElt)->GrpPerm
// first idea to compute the permutation stabiliezr of a collection of sets of cycles
{Compute the stabilizer of collectio of sets of cycles given by a sequence of indices, e.g. [1,2,3].}
  L_anz:=[n] cat [#x:x in part];
  part0:=[{m+1..m+L_anz[i]} where m:=&+L_anz[1..i-1]:i in [1..#L_anz]];
  m:=&+L_anz;
  A:=MatrixRing(Integers(),n+m)!0;
  j:=n+1;
  for L in part do
    for cyc in L do
      for i:=2 to #cyc do
        A[cyc[i-1],cyc[i]]+:=1;
      end for;
      // we close the loop in order to really have a cycle, maybe the same when computing the stabilizer directly?
      A[cyc[#cyc],cyc[1]]+:=1;
      for k in cyc do
        A[k,j]:=1;
      end for;
      j+:=1;
    end for;
  end for;
  vprintf FrameSymmetry: "  computing the full automorphism group on %o points: ",Nrows(A);
  vtime FrameSymmetry:
  grp0:=AutomorphismGroupMatrix(A);
  vprintf FrameSymmetry: "  computing stabilizer of %o\n",part0;
  for x in part0 do
    vprintf FrameSymmetry: "   current order: %o  (%o)",#grp0,x;
    vtime FrameSymmetry:
    grp0:=Stabilizer(grp0,x);
  end for;
  grp:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(grp0)]>;
  vprintf FrameSymmetry: "  group order: %o\n",#grp;
  return grp;
end intrinsic;

intrinsic FrameSymmetryKnownSubgroup(G::Mtrx,H0::GrpPerm:
  lowerbound:=1,
  sort_blocks:=false,
  eps:=1e-10)->GrpPerm
{Compute the symmetry group of a frame given by a canonical Gram matrix as a permutation group.
 When the optional parameter "lowerbound" is set, the algorithm terminates when the size symmetry group is not greater than the lowerbound.
 The parameter "sort_blocks" determines whether to sort the the blocks by decreasing size.
 The optional parameter "eps" defines the precision for floating point input.
}
// the input G is a Gram matrix
// CHECKING to be added

  ty:=Cputime();

  //check whether we have exact or floating point arthmetics
  F:=CoefficientRing(G);
  exact:=not (Type(F) eq FldRe or Type(F) eq FldCom);
  if exact then
    round:=func<x|x>;
  else
    M:=F!Round(1/eps);
    round:=func<x|Round(M*x)>;
  end if;

  n:=Nrows(G);
  G_round:=Matrix(n,n,[round(x):x in Eltseq(G)]);
  A:=Matrix([[i eq j or G_round[i,j] eq 0 select 0 else 1:j in [1..n]]:i in [1..n]]);
  frame_graph:=Graph<n|A>;

  //determine the permutation symmetry of the Gramian
  sym_gram:=AutomorphismGroupMatrix(G_round);
  H0_gram:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym_gram)]>;
  vprintf FrameSymmetry: "  symmetry of the Gramian, order: %o\n",#H0_gram;

  //determine the symmetry group of the 1-cycles and 2-cycles
  G_abs:=Matrix([[round(G[i,j]*G[j,i]):j in [1..n]]:i in [1..n]]);
  vprintf FrameSymmetry: "computing the automorphism group of the 1- and 2-cycles ";
  vtime FrameSymmetry:
  sym0:=AutomorphismGroupMatrix(G_abs);
  sym:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym0)]>;
  if IsSymmetric(sym) then
    vprintf FrameSymmetry: "  initial group is the full symmetric group Sym(%o)\n",n;
  elif #sym le 10^10 then
    vprintf FrameSymmetry: "  initial order: %o\n",#sym;
  else
    vprintf FrameSymmetry: "  initial order: 2^%o\n",Ilog2(#sym);
  end if;

  if #sym eq 1 then
    return sym;
  end if;

  H0:=sub<Sym(n)|H0,H0_gram>;
  transitive:=IsTransitive(H0);
  if transitive then
    //if H0 is transitive, restrict the inital group to the stabilizer of the first point
    //  it is sufficient that H0 is transitive, but is that also necessary?
    vprintf FrameSymmetry: "known transitive symmetry group of order %o\n",#H0;
    sym:=Stabilizer(sym,1);
  end if;

  //decide whether to report the order of the whole group or only that of the subgroup fixing the first point
  current_order:=#sym;
  vertices:=Vertices(frame_graph);
  vprintf FrameSymmetry,2: "computing a minimal cycle basis:\n";
  vtime FrameSymmetry,2:
  cycle_basis,edges:=MinimalCycleBasis(frame_graph);
  dim:=#cycle_basis;

  CycleVectorToIndices:=function(v_cycle);
  // convert the incidence vector of the cycle to a sequence of indices of vertices
    cycle:=EdgeCycleToSequence({edges[i]:i in Support(v_cycle)});
    cycle:=[Position(vertices,x):x in cycle];
    return cycle;
  end function;

  EdgeCycleProduct:=function(L);
  // compute the m-product corresponding to the edge cycle L
    L_vertices:=EdgeCycleToSequence(L);
    indices:=[Position(vertices,v):v in L_vertices];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  VertexCycleProduct:=function(L);
  // compute the m-product corresponding to the vertex cycle L
    indices:=[Position(vertices,v):v in L];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  CycleProduct:=function(L);
  // compute the m-product corresponding to the cycle L given by indices
    //initialisation by the closing of the cycle
    if L[1] eq L[#L] then
      prod:=1;
    else
      prod:=G[L[#L],L[1]];
    end if;
    for i:=2 to #L do
      prod*:=G[L[i-1],L[i]];
    end for;
    return prod;
  end function;

  //initialise the space of the cycles processed so far
  C:=sub<Universe(cycle_basis)|>;
  C_gen:=SparseMatrix(GF(2),0,#edges);

  if transitive then
  // when the group is transitive,
  // look at the symmetry group of the matrix of triple products with the first index fixed
    M:=Matrix(n,n,[round(G[1,j]*G[j,k]*G[k,1]):j,k in [1..n]]);
    vprintf FrameSymmetry: "  computing the symmetry of triple-products: ";
    vtime FrameSymmetry:
    sym_M:=AutomorphismGroupMatrix(M);
    sym:=sym meet sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym_M)]>;
    if exists{v:v in cycle_basis|Weight(v) gt 3} then
      for cyc in [v:v in cycle_basis|Weight(v) eq 3] do
        j:=Nrows(C_gen)+1;
        for i in Support(cyc) do
          SetEntry(~C_gen,j,i,1);
        end for;
      end for;
      Echelonize(~C_gen);
      C:=Image(Matrix(C_gen));
      cycle_basis:=[v:v in cycle_basis|Weight(v) gt 3];
    else
      cycle_basis:=[];
    end if;
  end if;

  test_cycle_basis_complete:=false;

  while #cycle_basis ne 0 and #sym gt lowerbound do
    v_cycle:=cycle_basis[1];

    if v_cycle in C then
      cycle_basis:=cycle_basis[2..#cycle_basis];
    else
      cycle:=CycleVectorToIndices(v_cycle);
      // if the known subgroup is transitive, then look for a cycle that contains the first point
      if transitive and #cycle_basis gt 1 and not 1 in cycle and not test_cycle_basis_complete then
        j:=2;
        repeat
          if cycle_basis[j] in C then
            Exclude(~cycle_basis,cycle_basis[j]);
          else
            cycle:=CycleVectorToIndices(cycle_basis[j]);
            if not 1 in cycle then
              j+:=1;
            end if;
          end if;
        until 1 in cycle or j gt #cycle_basis;
        test_cycle_basis_complete:=j gt #cycle_basis;
      end if;
      Exclude(~cycle_basis,v_cycle);

      m:=#cycle-1;
      vprintf FrameSymmetry: "selected the %o-cycle %o (%o s)\n",m,cycle[1..m],Cputime(ty);;
      vprintf FrameSymmetry: "  computing the orbit: ";
      vtime FrameSymmetry:
      if IsSymmetric(sym) then
        orbit:=GSet(sym,{x:x in Subsequences({1..n},m)|#Set(x) eq m});
      else
        orbit:=cycle[1..m]^sym;
      end if;
      vprintf FrameSymmetry:  "  orbit size: %o\n",#orbit;

      vprintf FrameSymmetry: "  computing the cycle-products: ";
      vtime FrameSymmetry:
      cycle_products:={<c,round(CycleProduct(c))>:c in orbit};

      values0:={* x[2]:x in cycle_products *};
      if sort_blocks then
        //sort the values by decreasing frequencies
        values:=Setseq(Set(values0));
        Sort(~values,func<a,b|Multiplicity(values0,b)-Multiplicity(values0,a)>);
      else
        values:=Set(values0);
      end if;
//      vprintf FrameSymmetry: "  %o different values (%o real)\n",#values,#{x:x in values|IsReal(x)};
      vprintf FrameSymmetry: "  %o different values\n",#values;

      part0:={{Append(x[1],x[1,1]):x in cycle_products|x[2] eq val}:val in values};
      // now look at the blocks of the orbit with a fixed value of the cycle-product
      for val in values do
        part:={x[1]:x in cycle_products|x[2] eq val};

        // add the cycles that have been processed to the cycle space
        //   can we use the known symmetry here?
        for cyc0 in part do
          cyc:=Append(cyc0,cyc0[1]);
          j:=Nrows(C_gen)+1;
          for i:=2 to #cyc do
            edge:=edges!{vertices[cyc[i-1]],vertices[cyc[i]]};
            SetEntry(~C_gen,j,Position(edges,edge),1);
          end for;
        end for;

        if #values gt 1 then
          vprintf FrameSymmetry: "  computing the stabilizer for %o points: ",#part;
          vtime FrameSymmetry:
          if transitive and not exists{x:x in part|not 1 in x} then
            orbit1:=Exclude(cycle[1..m],1)^sym;
//            orbit1:=GSet(sym,{Rotate(x,i):i in [0..#x-1],x in orbit1});
            part1:={Exclude(x,1):x in part};
//            part1:={Rotate(x,i):i in [0..#x-1],x in part1};
            sym:=Stabilizer(sym,orbit1,part1);
          else
            sym:=Stabilizer(sym,orbit,part);
          end if;
        end if;
        vprintf FrameSymmetry: "  current order: %o\n",#sym;
        Echelonize(~C_gen);
        C:=Image(Matrix(C_gen));
        vprintf FrameSymmetry: "  dimemsion of the cycle sub-space: %o (of %o)\n",Nrows(C_gen),dim;
        if Nrows(C_gen) eq dim then break;end if;
        // when the order has changed, we use the smaller group to compute the orbits
        if #sym ne current_order then
          current_order:=#sym;
          break;
        end if;
      end for;
    end if;
  end while;

  // now add the known automorphism
  sym:=sub<Generic(sym)|sym,H0>;

  return sym;
end intrinsic;


intrinsic FrameSymmetryTripleProducts(G::Mtrx:
  check_cycles:=false,   // a cycle basis is only computed when this is set to true
  check_triples:=false, 
  lowerbound:=1,
  sort_blocks:=false,
  eps:=1e-10)->GrpPerm
{Compute the symmetry group of a frame given by a canonical Gram matrix as a permutation group.
 When the optional parameter "lowerbound" is set, the algorithm terminates when the size symmetry group is not greater than the lowerbound.
 The parameter "sort_blocks" determines whether to sort the the blocks by decreasing size.
 The optional parameter "eps" defines the precision for floating point input.
}
// the input G is a Gram matrix
// CHECKING to be added

  ty:=Cputime();

  //check whether we have exact or floating point arthmetics
  F:=CoefficientRing(G);
  exact:=not (Type(F) eq FldRe or Type(F) eq FldCom);
  if exact then
    round:=func<x|x>;
  else
    M:=F!Round(1/eps);
    round:=func<x|Round(M*x)>;
  end if;

  vprintf FrameSymmetry: "computing the symmetry group of the Gramian:\n";
  ty:=Cputime();
  n:=Nrows(G);
  if exact then
    G_round:=G;
  else
    G_round:=Matrix(n,n,[round(x):x in Eltseq(G)]);
  end if;
  A:=Matrix([[i eq j or G_round[i,j] eq 0 select 0 else 1:j in [1..n]]:i in [1..n]]);
  frame_graph:=Graph<n|A>;
  components:=ConnectedComponents(frame_graph);
  components:=[[i:i in [1..n]|V[i] in x]:x in components] where V:=Vertices(frame_graph);

  //determine the permutation symmetry of the Gramian
  sym_gram:=AutomorphismGroupMatrix(G_round);
  H0_gram:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym_gram)]>;
  vprintf FrameSymmetry: "  symmetry of the Gramian, order: %o (%o s)\n",#H0_gram,Cputime(ty);
  delete G_round;

//determine the symmetry group of the 1-cycles and 2-cycles
  G_abs:=Matrix([[round(G[i,j]*G[j,i]):j in [1..n]]:i in [1..n]]);
  vprintf FrameSymmetry: "computing the automorphism group of the 1- and 2-cycles ";
  vtime FrameSymmetry:
  sym0:=AutomorphismGroupMatrix(G_abs);
  sym0:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym0)]>;
  sym:=sym0;
  if IsSymmetric(sym) then
    vprintf FrameSymmetry: "  initial group is the full symmetric group Sym(%o)\n",n;
  elif #sym le 10^10 then
    vprintf FrameSymmetry: "  initial order: %o\n",#sym;
  else
    vprintf FrameSymmetry: "  initial order: 2^%o\n",Ilog2(#sym);
  end if;
  delete G_abs;
  
  if #sym eq 1 then
    return sym;
  end if;

  H0:=sub<Sym(n)|H0_gram>;
  lowerbound:=Minimum(lowerbound,#H0);

  orbits_H0:={@ x[2]:x in OrbitRepresentatives(H0) @};

  vprintf FrameSymmetry: "computing the automorphism group of the 3-cycles\n";
  sym1:=H0;
  repeat
    t:=Rep(orbits_H0);
    vprintf FrameSymmetry: "  setting up the matrix for %o: ",t;
    vtime FrameSymmetry:
    M_triple:=Matrix([[round(G[t,i]*G[i,j]*G[j,t]):j in [1..n]]:i in [1..n]]);
    ind:=rep{x:x in components|t in x};
    if #components gt 1 then
      M_triple:=Matrix([[M_triple[i,j]:j in ind]:i in ind]);
    end if;
    vprintf FrameSymmetry: "  computing the triple-product subgroup fixing %o: ",t;
    vtime FrameSymmetry:
    stab:=AutomorphismGroupMatrix(M_triple);
    stab:=sub<Sym(n)|[[i in ind select ind[Position(ind,i)^g] else i:i in [1..n]]:g in Generators(stab)]>;
    if not stab subset sym then
      stab meet:=sym;
    end if;
    sym1:=sub<Generic(sym1)|sym1,stab>;
    vprintf FrameSymmetry: "  current order: %o\n",#sym1;
    orbits_H0 diff:=t^sym1;
    orbits_H0:={@ Rep(x):x in {t^sym1:t in orbits_H0} @};
  until #orbits_H0 eq 0;

  if not IsTransitive(sym1) then
    orbits_sym1:={@ x[2]:x in OrbitRepresentatives(sym1) @};
    // a single pair suffices for transitive action, so just compute the pairs when needed
    //vprintf FrameSymmetry: "  computing all matrices: ";
    //vtime FrameSymmetry:
    //L_M_triple:=[Matrix([[round(G[t,i]*G[i,j]*G[j,t]):j in [1..n]]:i in [1..n]]):t in orbits_sym1];
    for i:=1 to #orbits_sym1 do
      M1:=Matrix([[round(G[t,i]*G[i,j]*G[j,t]):j in [1..n]]:i in [1..n]]) where t:=orbits_sym1[i];
      ind1:=rep{x:x in components|orbits_sym1[i] in x};
      if #components gt 1 then
        M1:=Matrix([[M1[i,j]:j in ind1]:i in ind1]);
      end if;
      for j:=i+1 to #orbits_sym1 do
        if not IsConjugate(sym1,orbits_sym1[i],orbits_sym1[j]) then
          M2:=Matrix([[round(G[t,i]*G[i,j]*G[j,t]):j in [1..n]]:i in [1..n]]) where t:=orbits_sym1[j];
          ind2:=rep{x:x in components|orbits_sym1[j] in x};
          if #components gt 1 then
            M2:=Matrix([[M2[i,j]:j in ind2]:i in ind2]);
          end if;
          vprintf FrameSymmetry: "    checking %o -> %o: ",orbits_sym1[i],orbits_sym1[j];
  	  vtime FrameSymmetry:
          fl,tau:=IsIsomorphicMatrices(M1,M2);
	  //ASSUMPTION: when the restricted matrices are isomorphic, then this isomorphism is a symmetry
	  //   => WARNING: we might get a much larger group, e.g., when testing isomorphism of Gramians
	  if fl then
	    tau:=Sym(#ind1)![i^tau:i in [1..#ind1]];
	    L_tau:=[1..n];
	    for i:=1 to #ind1 do
 	       L_tau[ind1[i]]:=ind2[i^tau];
	       L_tau[ind2[i]]:=ind1[i^tau];
	    end for;
            tau:=Sym(n)!L_tau;
            assert tau in sym;
            vprintf FrameSymmetry: "  current order: %o\n",#sym1;
            vprintf FrameSymmetry: "  order tau: %o\n",Order(tau);
            sym1:=sub<Generic(sym1)|sym1,tau>;
            vprintf FrameSymmetry: "  new order: %o\n",#sym1;
	  end if;
	end if;
      end for;
    end for;
  end if;

  if not (check_cycles or check_triples) then
    return sym1;
  end if;

  vertices:=Vertices(frame_graph);
  vprintf FrameSymmetry,2: "computing a minimal cycle basis:\n";
  vtime FrameSymmetry,2:
  cycle_basis,edges:=MinimalCycleBasis(frame_graph);
  dim:=#cycle_basis;

  if exists{v:v in cycle_basis|Weight(v) eq 3} then
    sym:=sym1;
  end if;

  transitive:=IsTransitive(H0);

  if transitive then
    //if H0 is transitive, restrict the inital group to the stabilizer of the first point
    //  it is sufficient that H0 is transitive, but is that also necessary?
    vprintf FrameSymmetry: "known transitive symmetry group of order %o\n",#H0;
    sym:=Stabilizer(sym,1);
    lowerbound div:=n;
  end if;

  //decide whether to report the order of the whole group or only that of the subgroup fixing the first point
  current_order:=#sym;

  CycleVectorToIndices:=function(v_cycle);
  // convert the incidence vector of the cycle to a sequence of indices of vertices
    cycle:=EdgeCycleToSequence({edges[i]:i in Support(v_cycle)});
    cycle:=[Position(vertices,x):x in cycle];
    return cycle;
  end function;

  EdgeCycleProduct:=function(L);
  // compute the m-product corresponding to the edge cycle L
    L_vertices:=EdgeCycleToSequence(L);
    indices:=[Position(vertices,v):v in L_vertices];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  VertexCycleProduct:=function(L);
  // compute the m-product corresponding to the vertex cycle L
    indices:=[Position(vertices,v):v in L];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  CycleProduct:=function(L);
  // compute the m-product corresponding to the cycle L given by indices
    //initialisation by the closing of the cycle
    if L[1] eq L[#L] then
      prod:=1;
    else
      prod:=G[L[#L],L[1]];
    end if;
    for i:=2 to #L do
      prod*:=G[L[i-1],L[i]];
    end for;
    return prod;
  end function;

  //initialise the space of the cycles processed so far
  C:=sub<Universe(cycle_basis)|>;
  C_gen:=SparseMatrix(GF(2),0,#edges);


  //we only need to check the cycles of length larger than three
  //DOUBLE-CHECK!
  if not check_triples then
    if exists{v:v in cycle_basis|Weight(v) gt 3} then
      for cyc in [v:v in cycle_basis|Weight(v) eq 3] do
        j:=Nrows(C_gen)+1;
        for i in Support(cyc) do
          SetEntry(~C_gen,j,i,1);
        end for;
      end for;
      Echelonize(~C_gen);
      C:=Image(Matrix(C_gen));
      cycle_basis:=[v:v in cycle_basis|Weight(v) gt 3];
    else
      cycle_basis:=[];
    end if;
  end if;

  test_cycle_basis_complete:=false;

  while #cycle_basis ne 0 and #sym gt lowerbound do
    v_cycle:=cycle_basis[1];

    if v_cycle in C then
      cycle_basis:=cycle_basis[2..#cycle_basis];
    else
      cycle:=CycleVectorToIndices(v_cycle);
      // if the known subgroup is transitive, then look for a cycle that contains the first point
      if transitive and #cycle_basis gt 1 and not 1 in cycle and not test_cycle_basis_complete then
        j:=2;
        repeat
          if cycle_basis[j] in C then
            Exclude(~cycle_basis,cycle_basis[j]);
          else
            cycle:=CycleVectorToIndices(cycle_basis[j]);
            if not 1 in cycle then
              j+:=1;
            end if;
          end if;
        until 1 in cycle or j gt #cycle_basis;
        test_cycle_basis_complete:=j gt #cycle_basis;
      end if;
      Exclude(~cycle_basis,v_cycle);

      m:=#cycle-1;
      vprintf FrameSymmetry: "selected the %o-cycle %o (%o s)\n",m,cycle[1..m],Cputime(ty);;
      vprintf FrameSymmetry: "  computing the orbit: ";
      vtime FrameSymmetry:
      if IsSymmetric(sym) then
        orbit:=GSet(sym,{x:x in Subsequences({1..n},m)|#Set(x) eq m});
      else
        orbit:=cycle[1..m]^sym;
      end if;
      vprintf FrameSymmetry:  "  orbit size: %o\n",#orbit;

      vprintf FrameSymmetry: "  computing the cycle-products: ";
      vtime FrameSymmetry:
      cycle_products:={<c,round(CycleProduct(c))>:c in orbit};

      values0:={* x[2]:x in cycle_products *};
      if sort_blocks then
        //sort the values by decreasing frequencies
        values:=Setseq(Set(values0));
        Sort(~values,func<a,b|Multiplicity(values0,b)-Multiplicity(values0,a)>);
      else
        values:=Set(values0);
      end if;
//      vprintf FrameSymmetry: "  %o different values (%o real)\n",#values,#{x:x in values|IsReal(x)};
      vprintf FrameSymmetry: "  %o different values\n",#values;

      part0:={{Append(x[1],x[1,1]):x in cycle_products|x[2] eq val}:val in values};
      // now look at the blocks of the orbit with a fixed value of the cycle-product
      for val in values do
        part:={x[1]:x in cycle_products|x[2] eq val};

        // add the cycles that have been processed to the cycle space
        //   can we use the known symmetry here?
        for cyc0 in part do
          cyc:=Append(cyc0,cyc0[1]);
          j:=Nrows(C_gen)+1;
          for i:=2 to #cyc do
            edge:=edges!{vertices[cyc[i-1]],vertices[cyc[i]]};
            SetEntry(~C_gen,j,Position(edges,edge),1);
          end for;
        end for;

        if #values gt 1 then
          vprintf FrameSymmetry: "  computing the stabilizer for %o points: ",#part;
          vtime FrameSymmetry:
          if transitive and not exists{x:x in part|not 1 in x} then
            orbit1:=Exclude(cycle[1..m],1)^sym;
//            orbit1:=GSet(sym,{Rotate(x,i):i in [0..#x-1],x in orbit1});
            part1:={Exclude(x,1):x in part};
//            part1:={Rotate(x,i):i in [0..#x-1],x in part1};
            sym:=Stabilizer(sym,orbit1,part1);
          else
            sym:=Stabilizer(sym,orbit,part);
          end if;
        end if;
        vprintf FrameSymmetry: "  current order: %o\n",#sym;
        Echelonize(~C_gen);
        C:=Image(Matrix(C_gen));
        vprintf FrameSymmetry: "  dimemsion of the cycle sub-space: %o (of %o)\n",Nrows(C_gen),dim;
        if Nrows(C_gen) eq dim then break;end if;
        // when the order has changed, we use the smaller group to compute the orbits
        if #sym ne current_order then
          current_order:=#sym;
          break;
        end if;
      end for;
    end if;
  end while;

  // now add the known automorphism
  sym:=sub<Generic(sym)|sym,H0>;

  return sym;
end intrinsic;


intrinsic FrameSymmetryKnownSubgroupRecursive(G::Mtrx,H0::GrpPerm:
  lowerbound:=1,
  sort_blocks:=false,
  eps:=1e-10)->GrpPerm
{Compute the symmetry group of a frame given by a canonical Gram matrix as a permutation group.
 When the optional parameter "lowerbound" is set, the algorithm terminates when the size symmetry group is not greater than the lowerbound.
 The parameter "sort_blocks" dettermines whether to sort the the blocsk by increasing size.
 The optional parameter "eps" defines the precision for floating point input.
}
// the input G is a Gram matrix
// CHECKING to be added

  //check whether we have exact or floating point arthmetics
  F:=CoefficientRing(G);
  exact:=not (Type(F) eq FldRe or Type(F) eq FldCom);
  if exact then
    round:=func<x|x>;
  else
    M:=F!Round(1/eps);
    round:=func<x|Round(M*x)>;
  end if;

  n:=Nrows(G);
  G_round:=Matrix(n,n,[round(x):x in Eltseq(G)]);
  A:=Matrix([[i eq j or G_round[i,j] eq 0 select 0 else 1:j in [1..n]]:i in [1..n]]);
  frame_graph:=Graph<n|A>;

  //determine the permutation symmetry of the Gramian
  sym_gram:=AutomorphismGroupMatrix(G_round);
  H0_gram:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym_gram)]>;
  vprintf FrameSymmetry: "  symmetry of the Gramian, order: %o\n",#H0_gram;

  H0:=sub<Sym(n)|H0,H0_gram>;

  //restrict the intial group to the stabiliser of the point p0 picked above
  //  it is sufficient that H0 is transitive, but is that also necessary?
  transitive:=IsTransitive(H0);
  if transitive then
    vprintf FrameSymmetry: "known transitive symmetry group of order %o\n",#H0;
    G1:=Submatrix(G,2,2,n-1,n-1);
    H1:=Stabilizer(H0,1);
    H1_restricted:=sub<Sym(n-1)|[[j^g-1:j in [2..n]]:g in Generators(H1)]>;
    sym1:=FrameSymmetryKnownSubgroupRecursive(G1,H1_restricted:lowerbound:=lowerbound,sort_blocks:=sort_blocks,eps:=eps);
    sym1_embedded:=sub<Sym(n)|[[1] cat [1+j^g:j in [1..n-1]]:g in Generators(sym1)]>;
    sym:=sub<Sym(n)|H0,sym1_embedded>;
    return sym;
  end if;

  //determine the symmetry group of the 1-cycles and 2-cycles
  G_abs:=Matrix([[round(G[i,j]*G[j,i]):j in [1..n]]:i in [1..n]]);
  vprintf FrameSymmetry: "computing the automorphism group of the 1- and 2-cycles ";
  vtime FrameSymmetry:
  sym0:=AutomorphismGroupMatrix(G_abs);
  sym:=sub<Sym(n)|[[i^g:i in [1..n]]:g in Generators(sym0)]>;
  if IsSymmetric(sym) then
    vprintf FrameSymmetry: "  initial group is the full symmetric group Sym(%o)\n",n;
  elif #sym le 10^10 then
    vprintf FrameSymmetry: "  initial order: %o\n",#sym;
  else
    vprintf FrameSymmetry: "  initial order: 2^%o\n",Ilog2(#sym);
  end if;

  if #sym eq 1 then
    return sym;
  end if;

  current_order:=#sym;
  vertices:=Vertices(frame_graph);
  vprintf FrameSymmetry,2: "computing a minimal cycle basis:\n";
  vtime FrameSymmetry,2:
  cycle_basis,edges:=MinimalCycleBasis(frame_graph);
  dim:=#cycle_basis;

  EdgeCycleProduct:=function(L);
  // compute the m-product corresponding to the edge cycle L
    L_vertices:=EdgeCycleToSequence(L);
    indices:=[Position(vertices,v):v in L_vertices];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  VertexCycleProduct:=function(L);
  // compute the m-product corresponding to the vertex cycle L
    indices:=[Position(vertices,v):v in L];
    prod:=1;
    for i:=2 to #indices do
      prod*:=G[indices[i-1],indices[i]];
    end for;
    return prod;
  end function;

  CycleProduct:=function(L);
  // compute the m-product corresponding to the cycle L given by indices
    //initialisation by the closing of the cycle
    if L[1] eq L[#L] then
      prod:=1;
    else
      prod:=G[L[#L],L[1]];
    end if;
    for i:=2 to #L do
      prod*:=G[L[i-1],L[i]];
    end for;
    return prod;
  end function;

  //initialise the space of the cycles processed so far
  C:=sub<Universe(cycle_basis)|>;
  C_gen:=SparseMatrix(GF(2),0,#edges);

  while #cycle_basis ne 0 and #sym gt lowerbound do
    v_cycle:=cycle_basis[1];
    if v_cycle in C then
      cycle_basis:=cycle_basis[2..#cycle_basis];
    else
      // convert the incidence vector of the cycle to a sequence of indices of vertices
      cycle:=EdgeCycleToSequence({edges[i]:i in Support(v_cycle)});
      cycle:=[Position(vertices,x):x in cycle];
      // if H0 is transitive, pick a cycle in the orbit of H0 that contains p0
      m:=#cycle-1;
      vprintf FrameSymmetry: "selected the %o-cycle %o\n",m,cycle[1..m];
      vprintf FrameSymmetry: "  computing the orbit: ";
      vtime FrameSymmetry:
      if IsSymmetric(sym) then
        orbit:=GSet(sym,{x:x in Subsequences({1..n},m)|#Set(x) eq m});
      else
        orbit:=cycle[1..m]^sym;
      end if;
      vprintf FrameSymmetry:  "  orbit size: %o\n",#orbit;

      vprintf FrameSymmetry: "  computing the cycle-products: ";
      vtime FrameSymmetry:
      cycle_products:={<c,round(CycleProduct(c))>:c in orbit};

      values0:={* x[2]:x in cycle_products *};
      if sort_blocks then
        //sort the values by decreasing frequencies
        values:=Setseq(Set(values0));
        Sort(~values,func<a,b|Multiplicity(values0,b)-Multiplicity(values0,a)>);
      else
        values:=Set(values0);
      end if;
//      vprintf FrameSymmetry: "  %o different values (%o real)\n",#values,#{x:x in values|IsReal(x)};
      vprintf FrameSymmetry: "  %o different values\n",#values;

      // now look at the blocks of the orbit with a fixed value of the cycle-product
      for val in values do
        part:={x[1]:x in cycle_products|x[2] eq val};

        // convert the cycles to vectors in the cycle space
        // can we use the known symmetry here?
        for cyc0 in part do
          cyc:=Append(cyc0,cyc0[1]);
          j:=Nrows(C_gen)+1;
          for i:=2 to #cyc do
            edge:=edges!{vertices[cyc[i-1]],vertices[cyc[i]]};
            SetEntry(~C_gen,j,Position(edges,edge),1);
          end for;
        end for;

        if #values gt 1 then
          vprintf FrameSymmetry: "  computing the stabilizer for %o points: ",#part;
          vtime FrameSymmetry:
          sym:=Stabilizer(sym,orbit,part);
        end if;
        vprintf FrameSymmetry: "  current order: %o\n",#sym;
        Echelonize(~C_gen);
        C:=Image(Matrix(C_gen));
        vprintf FrameSymmetry: "  dimemsion of the cycle sub-space: %o (of %o)\n",Nrows(C_gen),dim;
        if Nrows(C_gen) eq dim then break;end if;
        // when the order has changed, we use the smaller group to compute the orbits
        if #sym ne current_order then
          current_order:=#sym;
          break;
        end if;
      end for;
    end if;
  end while;

  // now add the known automorphism
  sym:=sub<Generic(sym)|sym,H0>;

  return sym;
end intrinsic;

intrinsic permutation_to_matrix(V::Mtrx,perm::GrpPermElt)->Mtrx,Mtrx
{}
// compute a matrix M such that M*V = V^perm*D with a diagonal matrix D
  // the very lazy approach
  d:=Nrows(V);
  n:=Ncols(V);
  V_transp:=Transpose(V);
  V_perm:=Transpose(Matrix([V_transp[i]:i in Eltseq(perm)]));
  P<[x]>:=PolynomialRing(CoefficientRing(V),d^2+n-1,"grevlex");
  M:=MatrixRing(P,d)!x[1..d^2];
  // fix the first phase to be 1
  D:=DiagonalMatrix([1] cat x[d^2+1..d^2+n-1]);
  bed:=Reduce(Eltseq(M*ChangeRing(V,P)-ChangeRing(V_perm,P)*D));
  id:=ideal<P|bed>;
  dim,vars:=Dimension(id);
  if dim gt 0 then
    printf "WARNING: %o free parameters\n",dim;
    id:=ideal<P|id,[id.i-1:i in vars]>;
  end if;
  M1:=MatrixRing(P,d)![NormalForm(x,id):x in Eltseq(M)];
  D1:=DiagonalMatrix([NormalForm(x,id):x in Eltseq(D)]);
  if #bed eq Rank(id) then
    M1:=ChangeRing(M1,CoefficientRing(V));
    D1:=ChangeRing(D1,CoefficientRing(V));
  end if;
  return M1,D1,bed;
end intrinsic;


intrinsic permutation_to_matrix_grp(V::Mtrx,sym::GrpPerm)->Grp
{}
// compute a matrix M such that M*V = V^perm*D with a diagonal matrix D
  // the very lazy approach
  F:=CoefficientRing(V);
  d:=Nrows(V);
  n:=Ncols(V);
  V_transp:=Transpose(V);
  P<[x]>:=PolynomialRing(F,d^2+n-1,"grevlex");
  M:=MatrixRing(P,d)!x[1..d^2];
  // fix the first phase to be 1
  D:=DiagonalMatrix([1] cat x[d^2+1..d^2+n-1]);
  grp:=sub<GL(d,F)|>;
  for perm in Generators(sym) do
    V_perm:=Transpose(Matrix([V_transp[i]:i in Eltseq(perm)]));
    bed:=Reduce(Eltseq(M*ChangeRing(V,P)-ChangeRing(V_perm,P)*D));
    id:=ideal<P|bed>;
    dim,vars:=Dimension(id);
    if dim gt 0 then
      printf "WARNING: %o free parameters\n",dim;
      id:=ideal<P|id,[id.i-1:i in vars]>;
    end if;
    M1:=MatrixRing(F,d)![NormalForm(x,id):x in Eltseq(M)];
    D1:=DiagonalMatrix(F,[NormalForm(x,id):x in Eltseq(D)]);
    if not HasFiniteOrder(M1) then
      "WARNING: infinite order\n";
    end if;
    grp:=sub<Generic(grp)|grp,M1>;
  end for;
  return grp;
end intrinsic;

intrinsic IsIsomorphicFrame(G1::Mtrx,G2::Mtrx)->BoolElt,GrpPermElt
{Test whether the two frames given by their canonical Gramian G1 and G2 are equivalent.
 When they are equivalent, a permutation napping the first to the second is returned, otherwise the identity.
}

  n1:=Nrows(G1);
  n2:=Nrows(G2);
  G0:=Sym(n1+n2);

  if n1 ne n2 then
    return false, Id(G0);
  end if;
  sym:=FrameSymmetry(DirectSum(G1,G2));
  fl,perm:=IsConjugate(sym,{1..n1},{n1+1..n1+n2});
  if not fl then
    perm:=Id(G0);
  end if;
  return fl,perm;
end intrinsic;

intrinsic IsIsomorphicFrameMatrix(V1::Mtrx,V2::Mtrx)->BoolElt,GrpPermElt
{Test whether the two frames given by their defining matrices V1 and V2 are equivalent.
 When they are equivalent, a permutation napping the first to the second is returned, otherwise the identity.
}

  n1:=Ncols(V1);
  n2:=Ncols(V2);
  G0:=Sym(n1+n2);

  if n1 ne n2 then
    return false, Id(G0);
  end if;
  sym:=FrameSymmetryMatrix(DirectSum(V1,V2));
  fl,perm:=IsConjugate(sym,{1..n1},{n1+1..n1+n2});
  if not fl then
    perm:=Id(G0);
  end if;
  return fl,perm;
end intrinsic;