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  

ntypes:=Size(types);;
bifmat:=NullMat(ntypes,ntypes);
for i in [1..ntypes] do
  for j in [1..ntypes] do
    bifmat[i][j]:="";
  od;
od;
bifsubmat:=[];
gg:=Representative(types[ntypes]);;
# Order(gg);

# pick a symmetry type

for symnumber in [1..ntypes] do
Print("\nSymmetry type S",symnumber," ",IdGroup(Representative(types[symnumber])),"\n");

# get the symmetry type

S:=types[symnumber];;              
s:=Representative(S);;
sl:=AsList(s);;
ns:=Length(sl);;

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

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

# get the conjugacy classes of the subgroups of s

cc:=ConjugacyClassesSubgroups(s);;
# ccid:=List(cc,x->group2type(x));;
# Print("possible symmetry type indices\n",ccid,"\n");

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

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


# 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 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:=group2type(h);
       Print("    ",hid,", dim  Fix_P(U) K = ",Length(fix));
       Print(", dim  Fix_U_i K = ",irr[rep][1]*Length(fix)/Length(base));
       Print(", N_H(K)/K ",IdGroup(FactorGroup(Normalizer(s,h),h)),"\n");
       bifsubmat[hid]:=String(IdGroup(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 ",isots,",  maximals ");
 ma:=isots[Length(isots)];
 Unbind(isots[Length(isots)]);
 for i1 in isots do
   mal:=true;
   for i2 in isots do
     if i1<>i2 and subclass(types[i1],types[i2]) then
       mal:=false;
     fi;
   od;
   if mal then
     Add(jumps,i1);
     if bifmat[symnumber][i1]=0 then
       bifmat[symnumber][i1]:="";
     fi;
     bifmat[symnumber][i1]:=Concatenation(bifmat[symnumber][i1],bifsubmat[i1]);
#     bifmat[symnumber][i1]:=1;
     Print(i1," ");
#     Print(", factor group ",IdGroup(classfactor(types[ma],types[isots[1]])));
#     Print(", N_H(K)/K ",IdGroup(classnormalizer(types[i1],types[ma])),"\n");   
   fi;
 od;
 Print("\n");
 fi;
od;;
Sort(jumps);
Print("  Possible jumps from S",symnumber," ",jumps,"\n");

od;;

Print("\n");

LogTo("bifmat.txt");
for i in [1..ntypes] do
  Print(String(bifmat[i],15),"\n");
od;
LogTo();

Print("\n");