#This GAP code, for a given graph gamma evaluates the psi_c(gamma)

LoadPackage("grape");




psi:=function(n,c)
local A,i,j,sz,set_of_sizes, gamma,aut,colors,distinguishing_colorings,color,au,size_dist_colrs,psi_set,b,set_of_sizes_for_color,col, col_as_set_of_sets,alpha,  set_of_images_for_color;

####Construct the Path P_n

A:=NullMat(n,n);
for i in [1..n] do
    for j in [1..n] do
        if j=i+1 or j=i-1 then A[i][j]:=1; fi;
    od;
od;


######################### Functions needed
sz:=cell->Size(cell);
set_of_sizes:= partn->Set(List(partn,sz));



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



gamma:= Graph( Group(()), [1..n], OnPoints,function(x,y) return A[x][y]=1; end,true );
aut:=AutGroupGraph(gamma);



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

colors:=PartitionsSet([1..n],c);

if c>= n/2+1 then distinguishing_colorings:=colors; 
else

     distinguishing_colorings:=[];
     for color in colors do
        au:=aut;
        for j in [1..c-1] do
           au:=Stabilizer(au,color[j],OnSets);if Size(au)=1 then break; fi;
        od;
        if Size(au)=1 then Add(distinguishing_colorings,color); fi;
     od;
fi;


if Size(distinguishing_colorings)=0 then return 0;
else



size_dist_colrs:=Size(distinguishing_colorings);

psi_set:=[distinguishing_colorings[1]];

for i in [2..size_dist_colrs] do
     color:=distinguishing_colorings[i];
     b:=true;
     set_of_sizes_for_color:= set_of_sizes(color);
     for col in psi_set do 
           if  set_of_sizes(col)=set_of_sizes_for_color  then
                 col_as_set_of_sets:=Set([]);
                 for j in [1..c] do
                     col_as_set_of_sets:=Union(col_as_set_of_sets,[Set(col[j])]);
                 od;
                 for alpha in aut do
                      set_of_images_for_color:=Set([]);
                           for j in [1..c] do
                                set_of_images_for_color:=Union(set_of_images_for_color,[Set(OnSets(color[j],alpha))]);
                           od;
                       if  set_of_images_for_color=col_as_set_of_sets   then b:=false; break; fi;
                 od;
            fi;
      od;
      if b then Add(psi_set,color); fi;
od;
return Size(psi_set);
fi;


end;