#############################
# Given the order n, these functions calculate the 
# explicit fusion data and fusion table for the list 
# of association schemes of order n
#############################

#############################
# function producing Adjacency Matrices of a scheme,
# one of Hanaki and Miyamoto's GAP functions
#############################

AdjacencyMatrices := function(M)
    local A, i, j, k, n, d;

    n := Length(M);
    d := Maximum(M[1]);
    A := [];
    for i in [1..d+1] do
        A[i] := NullMat(n, n);
    od;
    for i in [1..n] do
        for j in [1..n] do
            A[M[i][j]+1][i][j] := 1;
        od;
    od;
    
    return A;
end;

###########################
# Calculates the Automorphism Group of a Scheme using 
# the implementation of nauty in GAP's Grape package,
# one of Hanaki and Miyamoto's functions
###########################
LoadPackage("grape");

AutomorphismGroupOfScheme:=function(R)
    local G, adj, gp, gr, n, x, y, i;
    
    adj := AdjacencyMatrices(R);
    n := Length(R);
    G := SymmetricGroup(n);
    for i in [2..(Length(adj) - 1)] do
        gr := Graph(Group(()), [1..n], OnPoints, function(x,y) return adj[i][x][y]=1; end);
        gp := AutomorphismGroup(gr);
        G := Intersection(G, gp);
    od;
  
    return G;
end;

###########################
# calculates the fusion of a scheme given a partition
# one of Hanaki and Miyamoto's functions
###########################

FusionScheme := function(M, L)
    local i, j, k, N, s, l, L2;

    if L[1] = [0] then
        L2 := L{[2..Length(L)]};
        L := L2;
    fi;
    s := Length(M);
    l := Length(L);
    N := NullMat(s, s);
    for i in [1..s] do
        for j in [1..s] do
            if i <> j then
                for k in [1..l] do
                    if M[i][j] in L[k] then
                        N[i][j] := k;
                    fi;
                od;
            fi;
        od;
    od;
      
    return N;
end;

#############################
# Gets Frequencies in 1st row of scheme matrix
#############################

getFrequencies:= function(M)
 local d, n, L, i, j;

  d:= Maximum( M );
  n:= Size(M);
  
  L:=[];
       
  for i in [1..d] do 
   L[i]:=0;
  od; 

for j in [1..d] do 
  for i in [1..n] do
       if M[i] = j then 
     L[j]:=L[j]+1;
    fi;
   od;
  od;

 Sort(L);
 return L;

end;

#############################
# Gets Frequencies in a partition
#############################

getFrequenciesP:= function(p)
 local k, L, i;

k:=Size(p);
L:=[];
for i in [1..k] do 
  L[i]:=Size(p[i]);
od;

Sort(L);
return L;
end;  



#######################################################
# The two ways we check for B[j] to be fusion of scheme B[i]
#######################################################
FusionCheckByGroups:=function(n,i,j)
local B, C, S, k, G, H, v, g, F, f, U, L, a, P, Q, Q1, g1, Y, B1, B2, K, m;

B:=AS[n][i]; 
C:=AS[n][j]; 
S:=SymmetricGroup(n);
k:=Maximum(C[1]);
L:=false; 
if not( TransposedMat(B)=B and not(TransposedMat(C)=C) ) then 

G:=AutomorphismGroupOfScheme(C); 
H:=AutomorphismGroupOfScheme(B);

if Size(G)/Size(H) in PositiveIntegers then 
  if IsSubgroup(G,H) then 
   g:=();
   B1:=B;  
   U:=AdjacencyMatrices(B1);
   a:=Algebra(Rationals,U);
   U:=Basis(a); 
   if C in a then 
     K:=Coefficients(U,C);
     if IsEqualSet([0..k],K) then 
       Q:=[];
       for m in [1..k] do 
         Q[m]:=Positions(K,m)-1; 
       od; 
       B2:=FusionScheme(B1,Q);
       if B2=C then 
         L:=[Q,g];
       fi;
     fi;
   fi;
  else
   F:=IsomorphicSubgroups(G,H); 
   if Size(F)>0 then  
   for f in F do 
     g:=RepresentativeAction(S,H,Image(f));   
     if not(g=fail) then 
       P:=PermutationMat(g,n);
       B1:=P^(-1)*B*P;  
       U:=AdjacencyMatrices(B1);
       a:=Algebra(Rationals,U);
       U:=Basis(a); 
       if C in a then 
          K:=Coefficients(U,C);
          if IsEqualSet([0..k],K) then 
            Q:=[];
            for m in [1..k] do 
              Q[m]:=Positions(K,m)-1; 
            od; 
            B2:=FusionScheme(B1,Q);
            if B2=C then 
               L:=[g,Q];
               break; 

            fi;
          fi;
        fi; 
      fi;
    od;       
  fi;
 fi;
fi;
fi;

return L; 
end;  
########################################
FusionCheckByPartitions:=function(n,i,j)
  local S, B, C, G, H, d, k, P, P1, L, Q, Q1, K, U, a, m, B1, B2, g;

S:=SymmetricGroup(n);
B:=AS[n][i];
C:=AS[n][j];
L:=false; 
if not( TransposedMat(B)=B and not(TransposedMat(C)=C) ) then 
G:=AutomorphismGroupOfScheme(C); 
d:=Maximum(B[1]);
k:=Maximum(C[1]);
P:=PartitionsSet([1..d],k);
P1:=[];
for Q in P do 
 if (getFrequencies(FusionScheme(B,Q)[1])=getFrequencies(C[1])) then
    AddSet(P1,Q);
 fi;
od;
for Q in P1 do
  B1:=FusionScheme(B,Q);
  H:=AutomorphismGroupOfScheme(B1); 
  if Size(H)=Size(G) then 
    if IsConjugate(S,G,H) then 
      g:=RepresentativeAction(S,G,H); 
      P:=PermutationMat(g,n);
      B2:=P*B1*P^(-1);
      if B2=C then 
        L:=[Q,g];
      else 
       if IsEqualSet(AdjacencyMatrices(B2),AdjacencyMatrices(C)) then 
          a:=Algebra(Rationals,AdjacencyMatrices(P*B*P^(-1)));
          U:=Basis(a);
          K:=Coefficients(U,C);
          Q1:=[];
          for m in [1..k] do 
            Q1[m]:=Positions(K,m)-1; 
          od; 
          L:=[g,Q1];
        fi;
       fi;
      fi;
    fi; 
od;
fi;

return L; 
end; 

################################
#  Fusion Table Refinement
################################
RefineFT:=function(F)
 local b, i, k, j1, j;

b:=Size(F);

for i in [1..b] do  
  k:=b-i+1;
for j1 in [1..(k-1)] do
  j:=k-j1;
  if j in F[k][2] then 
    SubtractSet(F[k][2],F[j][2]);    
  fi;
od;
od;
 
return F; 
end; 

##################################################
# The Explicit Fusion Function that produces 
# the fusion data and fusion table
##################################################
FusionTableN:= function(n)
  local B, F, b, i, M, j1, j, L, k;

B:=AS[n];
F:=[];

b:=Size(B);
for i in [1..b] do 
  F[i]:=[];
  F[i][1]:=i;
  F[i][2]:=[];
od;

for i in [2..b] do 
 M:=Maximum(B[i][1]);
 if M = 2 then
  AddSet(F[i][2],1); 
  continue;
 fi; 

 if M > 2 then 
 AddSet(F[i][2],1); 
 for j1 in [1..(i-2)] do 
   j:=i-j1; 
   if j in F[i][2] then 
      continue; 
   fi; 
   k:=Maximum(B[j][1]);

  if ( NrPartitionsSet([1..M],k)<100000 and M<n-1 )then 
   L:=FusionCheckByPartitions(n,i,j);
   if not(L=false) then 
     Print(i,"-->",j,":",L,"<br>\n");
     F[i][2]:=UnionSet(F[i][2],[j]);
     F[i][2]:=UnionSet(F[i][2],F[j][2]);
   fi; 
   else
   L:=FusionCheckByGroups(n,i,j);
   if not(L=false) then 
     Print(i,"-->",j,":",L,"<br>\n");
     F[i][2]:=UnionSet(F[i][2],[j]);
     F[i][2]:=UnionSet(F[i][2],F[j][2]);
   fi;
  fi;

 od; 
 fi;
Print(F[i],"\n"); 
od;

F:=RefineFT(F); 

return F; 
end; 
###############################
# Restarting functions, allowing for 
# starting at a random pair of schemes
###############################
SimpleFusionTablePartitions:=function(n,i,j)
  local B, b, d, j1, j2, k, L1, L2, L;

B:=AS[n];
b:=Size(B);
d:=Maximum(B[i][1]);
for j1 in [1..(j-1)] do 
  j2:=j-j1+1; 
  k:=Maximum(B[j2][1]);
   if NrPartitionsSet([1..d],k)<100000 then
  #  L1:=PartitionCompatibilityCheck(n,i,j);
  #  L2:=GroupCompatibilityCheck(n,i,j);
    L:=FusionCheckByPartitions(n,i,j2);
   else 
    L:=FusionCheckByGroups(n,i,j2);
   fi;
  Print(i,"-->",j2,":",L,"\n"); 
 
od;  

return b;
end; 
###############################
SimpleFusionTableGroups:=function(n,i,j)
  local B, b, d, j1, j2, k, L1, L2, L;

B:=AS[n];
b:=Size(B);
d:=Maximum(B[i][1]);
for j1 in [1..(j-1)] do 
  j2:=j-j1+1; 
  k:=Maximum(B[j2][1]);
  # if NrPartitionsSet([1..d],k)<100000 then
  #  L1:=PartitionCompatibilityCheck(n,i,j);
  #  L2:=GroupCompatibilityCheck(n,i,j);
  #  L:=FusionCheckByPartitions(n,i,j2);
  # else 
    L:=FusionCheckByGroups(n,i,j2);
  # fi;
  Print(i,"-->",j2,":",L,"\n"); 
 
od;  

return b;
end;