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


#############################
# function defining action of S_n on M(Q,n),
# ineffective for large n
#############################

TwoOrbitAction:=function(M,g)
local n,P,N;

n:=Size(M[1]);
P:=PermutationMat(g,n);
N:=P^(-1)*M*P; 

return N;
end;
#######################################
# One of Hanaki and Miyamoto's
#######################################
TwoOrbitConfiguration := function(G)
    local f, g, n, L, d, M, x, y, i;

    f:=IsomorphismPermGroup(G); 
    g:=Image(f);
    n:=LargestMovedPoint(g);
    L:=Orbits(g,Tuples([1..n],2),OnTuples);
    d:=Length(L); 
    M := NullMat(n, n);
    for x in [1..n] do
        for y in [1..n] do
            for i in [1..d] do
                if [x,y] in L[i] then
                    M[x][y] := i - 1;
                fi;
            od;
        od;
    od;

    return M;
end;






#############################
# 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;
###############################

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

n:=Length(M); 
d:=Maximum(M); 
A:=[];
for i in [1..(d+1)] do 
  A[i]:=[];
  for j in [1..n] do 
    A[i][j]:=0;
  od; 
od; 
for k in [1..n] do 
  for j in [1..n] do 
    
    A[M[k]+1][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;


#############################
# function replacing AutomorphismGroupOfScheme
# ineffective for large n, however
#############################

AutGp:=function(M)
local n,S,G;

n:=Size(M[1]);
S:=SymmetricGroup(n);
G:=Stabilizer(S,M,TwoOrbitAction);

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;

#############################
# function producing Basic Matrix of a Scheme 
# from the list of its Adjacency Matrices 
#############################

BasicMatrixOfScheme := function(A)
    local M, d, i;

    d:=Size(A);
    M:=0*A[1]; 
    for i in [1..(d-1)] do 
      M:=M+i*A[i+1]; 
    od; 
    return M; 
end;


#############################
# Get Frequencies from 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;

############################
# List [1, 2, 3, ... , n]
############################

getList:= function ( n )
  local L,i;
  L:=[];

  for i in [1..n] do
    Add(L,i);
  od;
return L;

end;
#############################
# Get Frequencies from 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;  


###############################
# Create copy of a List  "l"
###############################
createCopy:= function ( l )
  local L, i;

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

  return L;

end;

###############################
# Checks if nxn matrix is basic matrix of a scheme
###############################

IsScheme:= function(D)
  local d, B, A; 

d:=Maximum(D[1]);
B:=AdjacencyMatrices(D); 
A:=AlgebraWithOne(Rationals,B); 
if Dimension(A)=(d+1) then 
   return true; 
  else
   return false;
fi;
end;
#######################################
IsAssociationScheme:= function(D)
  local d, B, A, L, i, b, v; 

d:=Maximum(D[1]);
B:=AdjacencyMatrices(D); 
L:=[];
b:=Size(B);
for i in [1..b] do 
  L[i]:=TransposedMat(B[i]);
od;
v:=false;
if IsSubset(L,B) then 
A:=AlgebraWithOne(Rationals,B); 
if Dimension(A)=(d+1) then 
   v:=true; 
fi;
fi;

return v;
end;

#######################################################
# Checks if scheme B[j] can be a fusion of B[i]
#######################################################
GroupCompatibilityCheck := function(n,i,j)
 local B, C, G, H, j1, G1, i1, H1, T, f; 

B:=AS[n][i];
C:=AS[n][j];
#H:=AutGp(B);
H:=AutomorphismGroupOfScheme(B);
i1:=IsomorphismSpecialPcGroup(H);
H1:=Image(i1); 
#G:=AutGp(C);
G:=AutomorphismGroupOfScheme(C);
j1:=IsomorphismSpecialPcGroup(G);
G1:=Image(j1);
T:=false; 
f:=IsomorphicSubgroups(G1,H1);
if Size(f)>0 then 
#if Size(G)/Size(H) in PositiveIntegers then 
  T:=true; 
fi; 

return T;
end; 


#######################################################
# Checks if scheme B[j] can be a fusion of B[i]
#######################################################
BasicMatricesFusionCheck := function(B,C)
 local G, H, T, f; 

H:=AutGp(B);
#H:=AutomorphismGroupOfScheme(B); 
G:=AutGp(C);
#G:=AutomorphismGroupOfScheme(C);
T:=false; 
f:=IsomorphicSubgroups(G,H);
if Size(f)>0 then 
#if Size(G)/Size(H) in PositiveIntegers then 
  T:=true; 
fi; 

return T;
end; 



##################################
#  Partitions Compatibility Check
##################################

PartitionCompatibilityCheck := function(n,i,j)
  local B, C, B1, P, T, d, k, v;

B:=AS[n][i];
C:=AS[n][j];
d:=Maximum(B[1]);
k:=Maximum(C[1]);
T:=false; 
for P in PartitionsSet([1..d],k) do 
   B1:=FusionScheme(B,P);
   if getFrequencies(B1[1])=getFrequencies(C[1]) then 
     T:=true; 
     break; 
   fi; 
  od; 

return T; 
end; 

#######################################################
# Direct checks for B[j] to be fusion of scheme B[i]
#######################################################
FusionCheckByGroups:=function(n,i,j)
local B,C,S,k,G,H,v,g,f1,g1,G1,H1,H2,b,F,f,U,L,a,P,Q,Q1,Y,B1,B2,K,m;

B:=AS[n][i]; 
C:=AS[n][j]; 
S:=SymmetricGroup(n);
k:=Maximum(C[1]);
L:=false; 
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
   if IsSolvable(G) then 
    f1:=IsomorphismPcGroup(G);
    else
    f1:=IsomorphismPermGroup(G);
   fi;
   g1:=InverseGeneralMapping(f1);
   G1:=Image(f1);
   if IsSolvable(H) then 
     H1:=Image(IsomorphismPcGroup(H));
     else
     H1:=Image(IsomorphismPermGroup(H));
   fi;
   F:=IsomorphicSubgroups(G1,H1);
   if Size(F)>0 then 
     for f in F do 
       H2:=Image(g1,Image(f));
       b:=RepresentativeAction(S,H2,H);
       if not(b=fail) then 
       P:=PermutationMat(b^(-1),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:=[b,Q];
               break; 

            fi;
          fi;
        fi; 
      fi;
    od;       
  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];
G:=AutomorphismGroupOfScheme(C); 
d:=Maximum(B[1]);
k:=Maximum(C[1]);
#S1:=SymmetricGroup(k);
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;
L:=false; 
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;

#       T:=Orbit(SymmetricGroup(k),Q,Permuted);
#       for Q1 in T do 
#        B1:=FusionScheme(B,Q1); 
#        H:=AutomorphismGroupOfScheme(B1); 
#        g:=RepresentativeAction(S,G,H); 
#        P:=PermutationMat(g,n);
#        B2:=P*B1*P^(-1);
#        if B2=C then 
#          L:=[Q,g];
#          break;
#        fi;
#       od;  
#      fi;
#    fi; 
#  fi; 
#  if not(L=false) then 
#    break; 
#  fi;
#od;

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; 

################################
#  Fusion Table Reversal
################################
ReverseFT:=function(F)
 local b, i, j, G;

b:=Size(F);
G:=[];

for i in [1..b] do
  G[i]:=[];
  G[i][1]:=i;
  G[i][2]:=[];
od; 

for i in [1..b] do 
  for j in [1..b] do 
    if i in F[j][2] then 
       AddSet(G[i][2],j); 
    fi; 
  od; 
od; 

return G; 

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 k=M then 
      continue; 
   fi; 
#if CompatibilityCheck(n,i,j)=false then 
#   continue; 
#else
# if GroupCompatibilityCheck(n,i,j)=false then 
#    continue; 
# else
  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;
# fi;
#fi;  
 od; 
 fi;
Print(RefineFT(F),"\n"); 
od;

F:=RefineFT(F); 

return F; 
end; 
###############################
# Restarting functions, allowing for 
# checks with different method
###############################
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; 

########################
IsEqualSchemes:=function(B,C)
local t,n,S,H,G,g,P,B1;

t:=false;
n:=Length(B[1]);
S:=SymmetricGroup(n);
if getFrequencies(B[1])=getFrequencies(C[1]) then 
H:=AutomorphismGroupOfScheme(B);
G:=AutomorphismGroupOfScheme(C);
g:=RepresentativeAction(S,G,H);
if not(g=fail) then 
P:=PermutationMat(g,n);
B1:=P*B*P^(-1);
if IsEqualSet(AdjacencyMatrices(B1),AdjacencyMatrices(C)) then 
return true;
fi;
fi;
fi;

return t;
end;
##########################################
AsymmetricFusions:=function(B)
local L,n,A,d1,D,L1,L2,L3,i,B1,C,j,l1,l2,k1,k2,p,M,P,Q1,Q2,Q3,l,P1,L0;

L:=[];
n:=Length(B[1]);
A:=AdjacencyMatrices(B);
d1:=Size(A);
D:=[1..(d1-1)];
L1:=[];
L2:=[];
L3:=[];
for i in D do 
if TransposedMat(A[i+1])=A[i+1] then 
Add(L1,i);
else
Add(L2,i);
fi;
od;

#for i in D do 
#if not(i in L1) and not(i in L3) then 
#Add(L2,i);
#for j in D do  
#if not(j=i) and TransposedMat(A[i+1])=A[j+1] then
# Add(L3,j);
#fi;
#od;
#fi;
#od;
#Print(L1,L2,L3,"\n");

l1:=Size(L1);
l2:=Size(L2);
for k1 in [1..l1] do 
for Q1 in PartitionsSet(L1,k1) do
for k2 in [1..l2] do 
for P in PartitionsSet(L2,k2) do
#Q2:=[];
#Q3:=[];
#l:=Length(P);
#for i in [1..l] do 
#Q2[i]:=[];
#Q3[i]:=[];
#for p in P[i] do 
#Add(Q2[i],L2[p]);
#Add(Q3[i],L3[p]);
#od;
#od;
P1:=UnionSet(Q1,P);
B1:=FusionScheme(B,P1);
if IsAssociationScheme(B1) then 
Print(P1);
M:=getFrequencies(B1[1]);
for i in [2..Size(AS[n])] do 
if getFrequencies(AS[n][i][1])=M then 
if IsEqualSchemes(B1,AS[n][i]) then 
Print(i,"\n");
Add(L,i);
fi;  
fi;
od; 
fi;
od; 
od;
od;
od;

L0:=AsSet(L);
return L0;
end;

##########################################
SymmetricFusions:=function(B)
local L,n,A,d1,D,L1,L2,L3,i,B1,C,j,l1,l2,k1,k2,p,M,P,Q1,Q2,Q3,l,P1,L0;

L:=[];
n:=Length(B[1]);
A:=AdjacencyMatrices(B);
d1:=Size(A);
D:=[1..(d1-1)];
L1:=[];
L2:=[];
L3:=[];
for i in D do 
if TransposedMat(A[i+1])=A[i+1] then 
Add(L1,i);
fi;
od;

for i in D do 
if not(i in L1) and not(i in L3) then 
Add(L2,i);
for j in D do  
if not(j=i) and TransposedMat(A[i+1])=A[j+1] then
 Add(L3,j);
fi;
od;
fi;
od;

l1:=Size(L1);
l2:=Size(L2);
for k1 in [1..l1] do 
for Q1 in PartitionsSet(L1,k1) do
for k2 in [1..l2] do 
for P in PartitionsSet([1..l2],k2) do
Q2:=[];
Q3:=[];
l:=Length(P);
for i in [1..l] do 
Q2[i]:=[];
for p in P[i] do 
Add(Q2[i],L2[p]);
Add(Q2[i],L3[p]);
od;
od;
P1:=UnionSet(Q1,Q2);
B1:=FusionScheme(B,P1);
if IsAssociationScheme(B1) then 
Print(P1);
M:=getFrequencies(B1[1]);
for i in [2..Size(AS[n])] do 
if getFrequencies(AS[n][i][1])=M then 
if IsEqualSchemes(B1,AS[n][i]) then 
Print(i,"\n");
Add(L,i);
fi;  
fi;
od; 
fi;
od; 
od;
od;
od;

L0:=AsSet(L);

return L0;
end;

###################################

AllFusions:=function(B)
local L,n,A,d1,D,k,P,P1,n1,i,B1,C,j,M,L0;

L:=[];
n:=Length(B[1]);
A:=AdjacencyMatrices(B);
d1:=Length(A);
D:=[1..(d1-1)];

for k in D do 
P1:=PartitionsSet(D,d1-k);
n1:=Size(P1);
for j in [1..n1] do
P:=P1[j]; 
B1:=FusionScheme(B,P);
if IsAssociationScheme(B1) then 
Print(P);
M:=getFrequencies(B1[1]);
for i in [2..Size(AS[n])] do 
if getFrequencies(AS[n][i][1])=M then 
if IsEqualSchemes(B1,AS[n][i]) then 
Print(i,"\n");
Add(L,i);
fi;  
fi;
od; 
fi;
od; 
od;

L0:=AsSet(L);
return L0;
end;