SizeScreen([300,1000]);; 

Read("types.gap");
Read("perms2matrices.gap");
Read("matinbasis.gap");
Read("subclass.gap");
Read("classfactor.gap");
Read("classnormalizer.gap");

# n = number of points on which gg acts
# types = collection of conjugacy classes  

# collect all the subtypes

subtypes:=[];;
for S in types do
  for s in S do
    Add(subtypes,s);
  od;
od;
subtypes;;
nsubtypes:=Size(subtypes);;
ntypes:=Size(types);;

# create and empty bifmat, this will have the bif digraph info

bifmat:=NullMat(nsubtypes,nsubtypes);;
for i in [1..nsubtypes] do
  for j in [1..nsubtypes] do
    bifmat[i][j]:="";
  od;
od;
bifsubmat:=[];;

biftypemat:=NullMat(ntypes,ntypes);;
for i in [1..ntypes] do
  for j in [1..ntypes] do
    biftypemat[i][j]:="";
  od;
od;

gg:=Representative(types[ntypes]);;
# Order(gg);

# function to find the stabilizer

stabilizer:=function(sinb,fix)
  local a,b,good,stab;
  stab:=[];
  for a in sinb do
     good:=1;
     for b in fix do
        if b*a<>b then
	   good:=0;
	fi;
     od;
     if good=1 then
       Add(stab,a);
     fi;
  od;
  return(stab);
end;;

# function to get Z2 from [2,1]

groupid2str:=function(id)
if id="[ 2, 1 ]" then
  return("Z2");
fi;
if id="[ 3, 1 ]" then
  return("Z3");
fi;
if id="[ 4, 1 ]" then
  return("Z4");
fi;
if id="[ 6, 1 ]" then
  return("D3");
fi;
if id="[ 8, 3 ]" then
  return("D4");
fi;
if id="[ 12, 4 ]" then
  return("D6");
fi;
return(id);
end;;


# find the symmetry type of a conjugacy class

group2type:=function(h)
  local i;
  for i in [1..ntypes] do
    if IsConjugate(gg,Representative(types[i]),h) then
      return(i);
    fi;
  od; 
  return(-1);
end;;

# find the symmetry of a conjugacy class

group2subtype:=function(h)
  local i;
  for i in [1..nsubtypes] do
    if h=subtypes[i] then
      return(i);
    fi;
  od; 
  return(-1);
end;;

# function to convert from symmetry to symmetry type

sym2type:=function(id)
  return(group2type(subtypes[id]));
end;;

# function to find the stabilizer

stabilizergr:=function(s,fix,n,base)
  local a,b,good,stab,sg;
  stab:=[];
#  sg:=GeneratorsOfGroup(s);
  sg:=AsList(s);
  for a in sg do
     good:=1;
     for b in fix do
        if b*matinbasis(perm2matrix(a,n),base)<>b then
	   good:=0;
	fi;
     od;
     if good=1 then
       Add(stab,a);
     fi;
  od;
  return(Group(stab));
end;;


##############  main loop

# pick a symmetry type

for symnumber in [1..nsubtypes] do
s:=subtypes[symnumber];;
sl:=AsList(s);;
ns:=Length(sl);;

Print("\nSymmetry H=H",symnumber,", Symmetry type S", group2type(s),", Id(H)=",IdGroup(s),"\n");

# get the character table info

tb:=CharacterTable(s);;       
irr:=Irr(tb);;
nirr:=Length(Irr(tb));;         
cctb:=ConjugacyClasses(tb);;
ncctb:=Length(cctb);;

# find the projections for each irred reps as matrices

projs:=[];;                   
for i in [1..nirr] do           
   proj:=0*IdentityMat(n);
   for j in [1..ncctb] do     # each class
      cc:=AsList(cctb[j]);
      for k in [1..Length(cc)] do    # each element in the class
         proj:=proj +
	       perm2matrix(cc[k],n)*irr[i][j];
      od;
   od;
   proj:=irr[i][1]/ns*proj;
   Add(projs,proj);
od;

# find bases for the projected spaces

bases:=List(projs,x->BaseMat(x));;

# get the conjugacy classes of the subgroups of s

cc:=ConjugacyClassesSubgroups(s);;

# find the fixed point subspaces

jumps:=[];;                 # this will have the possible jumps 
factors:=[];;               # symmetry of the bifurcation
for rep in [1..nirr] do     # pick a representation
 base:=bases[rep];
 isots:=[];                 # this will have the isotropy subgroup indices
 if Length(base)<>0 then
#   Print(base,"basis\n"); 
   Print("  Representation ",rep,", dimension of irred subspace U_i = ",irr[rep][1],"\n");
#   Print(" dimension of projected space ",Length(base),"\n");
   if not ForAll(Flat(projs[rep]),x-> (ComplexConjugate(x)=x) ) then
      Print("  NON-REAL projection\n");
   fi;
#   Print("isotropy subgroups\n");
   sinb:=[];
   for x in sl do
     Add(sinb,matinbasis(perm2matrix(x,n),base));
   od;
   for c in AsList(cc) do           # pick a subgroup class of s
#   for h in c do
     h:=Representative(c);
     hinb:=[];
     for x in AsList(h) do
       Add(hinb,matinbasis(perm2matrix(x,n),base));
     od;
     fixes:=List(hinb,x-> NullspaceMat( x-IdentityMat(Size(x)) ) );   # fixedpoint spaces of elements
     fix:=IdentityMat(Size(fixes[1]));                                # fixedpoint space of h
     for i in [1..Length(fixes)] do
       fix:=SumIntersectionMat(fix,fixes[i])[2];
     od;
#     Print("testing ",group2type(h),"\nstab= ",stabilizergr(s,fix,n,base),"\n k= ",h,"\n");
     if stabilizergr(s,fix,n,base)=h then            # check if stab(fix(h))=h
       hid:=group2subtype(h);
       Print("    K=H",hid,", dim  Fix_P(U) K = ",Length(fix));
       Print(", dim  Fix_U_i K = ",irr[rep][1]*Length(fix)/Length(base));
       Print(", id(N_H(K)/K)=",IdGroup(FactorGroup(Normalizer(s,h),h)),"\n");
       bifsubmat[hid]:=String(Size(FactorGroup(Normalizer(s,h),h)));
#######       Print("     Base for P(Psi) ",BaseMat(psi*projs[rep]),"\n");          
       if hid=-1 then
          Print(c,"\n");
       fi;
       AddSet(isots,hid);
     fi;
#   od;  
   od; 

 # find the maximal isotropy subgroups
 Sort(isots);
 Print("    isotropy subgroups H",isots,",  maximals ");
 ma:=isots[Length(isots)];
 mi:=isots[1];
 fgid:=IdGroup(FactorGroup(subtypes[ma],subtypes[mi]));
 Unbind(isots[Length(isots)]);
 for i1 in isots do
   mal:=true;
   for i2 in isots do
     if i1<>i2 and IsSubgroup(subtypes[i2],subtypes[i1]) then
       mal:=false;
     fi;
   od;
   if mal then
     Add(jumps,i1);
     if bifmat[symnumber][i1]=0 then
       bifmat[symnumber][i1]:="";
     fi;
     if bifmat[symnumber][i1] <> "" then
        bifmat[symnumber][i1] := Concatenation(bifmat[symnumber][i1],",");
	bifmat[symnumber][i1] := "";  # coment this out to list all occurences
     fi;
     bifmat[symnumber][i1]:=Concatenation(bifmat[symnumber][i1],bifsubmat[i1]);
     Print(i1," ");
     biftypemat[sym2type(symnumber)][sym2type(i1)]:=Concatenation(bifmat[symnumber][i1]," ",
                                                                  groupid2str(String(fgid)));
   fi;
 od;
 Print(", factor group id(H/H')=",fgid);
 Print("\n");
 fi;
od;;
Sort(jumps);
typejumps:=List(jumps,x->group2type(subtypes[x]));
Print("  Jumps from H",symnumber," to H", jumps,"\n");
Print("  Typejumps from S",group2type(subtypes[symnumber])," to S",typejumps ,"\n");

od;;

Print("\nfromH : fromS : toK size(N_H(K)/K) ... \n\n");

LogTo("bifmat.txt");
for i in [1..nsubtypes] do
  Print(i," : "); 
  Print(group2type(subtypes[i]), " : ");
  for j in [1..nsubtypes] do
    if bifmat[i][j] <> "" then
      Print(j," ",bifmat[i][j],"   ");
    fi;
  od; 
  Print("\n");
od;
LogTo();

Print("\n");


Print("\nfromH   toK size(N_H(K)/K) ... one from each type\n\n");
for i in [1..nsubtypes] do
  Print(i,"   "); 
  usedtypes:=[];
  for j in [1..nsubtypes] do
    if bifmat[i][j] <> "" then
      if not sym2type(j) in usedtypes  then
        Add(usedtypes, sym2type(j));
        Print(j," ",bifmat[i][j],"   ");
      fi;
    fi;
  od; 
  Print("\n");
od;
Print("\n");


Print("\nfromS : toS size(N_H(K)/K) id(H/H')... \n\n");
for i in [1..ntypes] do
  Print(i," : "); 
  for j in [1..ntypes] do
    if biftypemat[i][j] <> "" then
      Print(j," ",biftypemat[i][j],"   ");
    fi;
  od; 
  Print("\n");
od;

Print("\n");