#%Format: latex
#% This file is marked for TeX processing
#[caption]"Skeleton manipulation in \GAP"
#================================================================
#: A skeleton is stored as a record, with fields computed on demand
#: Initialized are "Count" (number of edges) and%
#:   permutations $"X"=\nx$, $"Y"=\op$
#: Other fields are accessed \via\ "get(S, name)";%
#:   after first call, "S.name" will also do
#================================================================
#: A home grown VMT and virtual method dispatcher
#: For a field "name" let "get_.name := [f, ...]", where%
#:  "f = function(S, arg)" and other members are%
#:  `field dependent' arguments:%
#:   "f" receives "get_.name" in "arg"
#: Probably, there is a more civilized way to do this in \GAP\ \dots
#----------------------------------------------------------------
get_ := rec(XY := [function(S, arg) return S.X*S.Y; end]);  #"0XY"
get := function(S, name)  #"get"
 if not IsBound(S.(name)) then S.(name) := get_.(name)[1](S, get_.(name)); fi;
 return S.(name);
end;
#================================================================
#: This function creates a new skeleton
#: Usage: "Skeleton([n, ]x, y{, name1, value1})", where%
#:   "n" is the number of edges,%
#:   "x", "y" are two permutations (possibly as lists), and%
#:   "name1, value1" \etc\. are optional user data to be stored
#----------------------------------------------------------------
Skeleton := function(arg)  #"Skeleton"
 local n, s, x, y, perm, res;
 perm := function(p)
#% Recode to "list" and back to save memory
  if IsPerm(p) then p := ListPerm(p); fi;
  p := PermList(p);
  if p <> fail then n := Maximum(n, LargestMovedPoint(p)); fi;
  return p;
 end;
 n := 0;#%
 s := 0;
 if arg[1] in Integers then n := arg[1]; s := 1; fi;
 x := perm(arg[s+1]);#%
 y := perm(arg[s+2]);
 if (x = fail) or (y = fail) then return fail; fi;
 res := rec(Count := n, X := x, Y := y);
 for n in [s+3, s+5..Length(arg) - 1] do res.(arg[n]) := arg[n+1]; od;
 return res;
end;
#================================================================
#: Here are more computable fields
#----------------------------------------------------------------
#: "XCounts", "YCounts", "XYCounts": numbers of vertices/regions%
#:   as lists (by sizes)
#----------------------------------------------------------------
get_.XYCounts := [function(S, arg)            #"0XYCounts"
 local l, i, s;
 l := List(CycleStructurePerm(get(S, arg[2]))); s := 0;
 for i in [1..Length(l)] do
  if not IsBound(l[i]) then l[i] := 0; fi;#%
  s := s + (i+1)*l[i];
 od;
 return Concatenation([S.Count - s], l);  #: Include length 1 orbits as well
end, "XY"];
get_.XCounts := [get_.XYCounts[1], "X"]; #"0XCounts"
get_.YCounts := [get_.XYCounts[1], "Y"]; #"0YCounts"
#----------------------------------------------------------------
#: "XCount", "YCount", "XYCount": total numbers of vertices/regions
#----------------------------------------------------------------
get_.XYCount := [function(S, arg)            #"0XYCount"
 return Sum(get(S, arg[2]));
end, "XYCounts"];
get_.XCount := [get_.XYCount[1], "XCounts"]; #"0XCount"
get_.YCount := [get_.XYCount[1], "YCounts"]; #"0YCount"
#----------------------------------------------------------------
#: "XOrbits", "YOrbits", "XYOrbits":%
#:   lists of orbits of $\x$, $\y$, $\xy$
#----------------------------------------------------------------
get_.XYOrbits := [function(S, arg)       #"0XYOrbits"
 return OrbitsDomain(Group(get(S, arg[2])), [1..S.Count]);
end, "XY"];
get_.XOrbits := [get_.XYOrbits[1], "X"]; #"0XOrbits"
get_.YOrbits := [get_.XYOrbits[1], "Y"]; #"0YOrbits"
#----------------------------------------------------------------
#: "Group": the permutation group $G$ generated by $\x$, $\y$
#: "Stabilizer": the stabilizer $H:=\stab_G1\subset G$
#: $"Aut" = \Aut\skeleton = \Norm_G(H)/H$
#----------------------------------------------------------------
get_.Group := [function(S, arg) return Group(S.X, S.Y); end]; #"0Group"
get_.Stabilizer := [function(S, arg)  #"0Stabilizer"
 return Stabilizer(get(S, "Group"), 1);
end];
get_.Aut := [function(S, arg)  #"0Aut"
 local H;
 H := get(S, "Stabilizer");
 return FactorGroup(Normalizer(get(S, "Group"), H), H);
end];
#----------------------------------------------------------------
#: "Components": the list of connected components of $\skeleton$
#: "ComponentCount": the number of components
#: "IsConnected" returns "true" if $\skeleton$ is connected%
#:   (silently assumed by other methods)
#----------------------------------------------------------------
get_.Components := [function(S, arg)  #"0Components"
 return OrbitsDomain(get(S, "Group"), [1..S.Count]);
end];
get_.ComponentCount := [function(S, arg)  #"0ComponentCount"
 return Length(get(S, "Components"));
end];
get_.IsConnected := [function(S, arg)  #"0IsConnected"
 return get(S, "ComponentCount") = 1;
end];
#----------------------------------------------------------------
#: "Genus": the genus of $\skeleton$
#: "Type": the `minimal' type~$\stype$ of~$\skeleton$
#: "IsRegular" returns "true" if $\skeleton$ is $\stype$-regular
#----------------------------------------------------------------
get_.Genus := [function(S, arg)  #"0Genus"
 return (2+S.Count-get(S, "XCount")-get(S, "YCount")-get(S, "XYCount"))/2;
end];
get_.Type := [function(S, arg)  #"0Type"
 local lcm;
 lcm := function(l) return Lcm(Filtered([1..Length(l)], i -> l[i] > 0)); end;
 return [lcm(get(S, "XCounts")), lcm(get(S, "YCounts"))];
end];
get_.IsRegular := [function(S, arg)  #"0IsRegular"
 local t;
 t := get(S, "Type");
 return (PositionProperty(get(S, "XCounts"), i -> i > 0) = t[1]) and
        (PositionProperty(get(S, "YCounts"), i -> i > 0) = t[2]);
end];
#================================================================
#: "Get" provides conventional access to the fields:%
#:   \eg, after "XY := Get('XY')"%
#:   one can use "XY(S)" instead of "get(S, 'XY')"
#----------------------------------------------------------------
Get := function(name) return function(S) return get(S, name); end; end; #"Get"
#----------------------------------------------------------------
#% Here are the shortcuts
#%----------------------------------------------------------------
Count := Get("Count"); #"Count"#%
XY := Get("XY");   #"XY"#: We cannot assign to "X"
Black := Get("XOrbits"); #"Black"#%
Blacks := Get("XCounts"); #"Blacks"#%
BlackNo := Get("XCount"); #"BlackNo"
White := Get("YOrbits"); #"White"#%
Whites := Get("YCounts"); #"Whites"#%
WhiteNo := Get("YCount"); #"WhiteNo"
Reg := Get("XYOrbits"); #"Reg"#%
Regs := Get("XYCounts"); #"Regs"#%
RegNo := Get("XYCount"); #"RegNo"
Genus := Get("Genus"); #"Genus"#%
Stab := Get("Stabilizer"); #"Stab"#%
Aut := Get("Aut"); #"Aut"
IsConnected := Get("IsConnected"); #"IsConnected"#%
Components := Get("Components"); #"Components"
IsRegularSkeleton := Get("IsRegular"); #"IsRegularSkeleton"#%
Type := Get("Type"); #"Type"
#: These functions are for printing
PartitionNotation := function(list) #"PartitionNotation" #: A list in the partition notation
 return NormalizedWhitespace(JoinStringsWithSeparator(List([1..Length(list)],
     function(i)
      if list[i] = 0 then return ""; fi;
      return Concatenation(String(i), "^", String(list[i]));
     end), " "));
end;
RegStr := function(S) return PartitionNotation(Regs(S)); end;  #"RegStr" #: Regions in PN
get_.info := [function(S, arg)  #"0info" #: General info on the skeleton
 return [S.Count, Whites(S)[1], Blacks(S)[1], RegStr(S), Size(Aut(S))];
end];
SkeletonInfo := Get("info");   #"SkeletonInfo"
#================================================================
#: Below are a few honest functions that do not store the results
#----------------------------------------------------------------
#: "stab(S, e)": the stabilizer $\stab_G e$
#: "Component(S, c)": the $c$-th component (a sublist of "[1..S.Count]")
#: "ExtractComponent(S, c)" extracts the $c$-th component%
#:   as a new skeleton ($c$ being an integer or the actual component%
#:   returned by "Component(S, c)")
#----------------------------------------------------------------
stab := function(S, e) return Stabilizer(get(S, "Group"), e); end; #"stab"
Component := function(S, c) return get(S, "Components")[c]; end; #"Component"
ExtractComponent := function(S, c)  #"ExtractComponent"
 local p;
 if IsInt(c) then c := Component(S, c); fi;
 p := MappingPermListList([1..Length(c)], c);
 return Skeleton(Length(c), p*RestrictedPerm(S.X, c)*p^-1,
     p*RestrictedPerm(S.Y, c)*p^-1, "parent", S, #: original skeleton and inclusion
     "inc", l -> List(l, i -> c[i]));
end;
#================================================================
#: "Prod(S1, S2)" returns Cartesian product $\skeleton_1\times\skeleton_2$
#: "IsIsomorphic(S1, S2)" returns "true" if $\skeleton_1\cong\skeleton_2$%
#:   (assuming both connected)
#----------------------------------------------------------------
#: Cartesian product of two permutations
prod := function(p1, n1, p2, n2)  #"prod"
 local i, j, res;
 res := [1..n1*n2];
 for i in [1..n1] do for j in [1..n2] do
  res[(i - 1)*n2 + j] := (i^p1 - 1)*n2 + j^p2;
 od; od;
 return PermList(res);
end;
Prod := function(S1, S2)  #"Prod"
 return Skeleton(S1.Count*S2.Count,
   prod(S1.X, S1.Count, S2.X, S2.Count), prod(S1.Y, S1.Count, S2.Y, S2.Count),
   "factor", [S1, S2],   #: Store the factors
   "proj", [             #: The canonical projections to the factors
       l -> List(l, i -> QuoInt(i - 1, S2.Count) + 1),
       l -> List(l, i -> RemInt(i - 1, S2.Count) + 1)]);
#%       function(l) return List(l, i -> QuoInt(i - 1, S2.Count) + 1); end,
#%       function(l) return List(l, i -> RemInt(i - 1, S2.Count) + 1); end]);
end;
IsIsomorphic := function(S1, S2)  #"IsIsomorphic"
 if S1.Count <> S2.Count then return false; fi;
 return S1.Count in List(get(Prod(S1, S2), "Components"), Length);
end;
#================================================================
#: Morphisms of skeletons are stored as lists:%
#:   "List([1..S.Count], i -> f(i))"
#: "Morphisms(S1, S2)" lists all morphisms $\skeleton_1\to\skeleton_2$
#: "IsMorphism(S1, S2, m)" returns "true" if "m" is a morphism%
#:   $\skeleton_1\to\skeleton_2$
#: "IsAutomorphism(S, a)" returns "true" if $a\in\Aut\skeleton$%
#:   ("a" may also be a permutation)
#----------------------------------------------------------------
Morphisms := function(S1, S2)  #"Morphisms"
 local S, f, p;
 f := function(l)
  p := PermList(S.proj[1](l))^-1;  l := S.proj[2](l);
  return List([1..S1.Count], i -> l[i^p]);
 end;
 if not IsInt(S1.Count/S2.Count) then return []; fi;
 S := Prod(S1, S2);
 return List(Filtered(get(S, "Components"), i -> Length(i) = S1.Count), f);
end;
IsMorphism := function(S1, S2, m)  #"IsMorphism"
 local i;
 for i in [1..S1.Count] do
  if (m[i^S1.X] <> m[i]^S2.X)or(m[i^S1.Y] <> m[i]^S2.Y) then return false; fi;
 od;
 return true;
end;
IsAutomorphism := function(S, a)   #"IsAutomorphism"
 if IsList(a) then a := PermList(a); fi;
 return (a*S.X = S.X*a) and (a*S.Y = S.Y*a);
end;
#================================================================
#: "MapOrbits(p)" converts a list of orbits to a map
#: "FactorSkeleton(S, m)": factor of $\skeleton$ by "m". Here,%
#:    "m" is either a morphism (see above),%
#:    or a permutation (possibly given by a list),%
#:    or a list of orbits
#----------------------------------------------------------------
MapOrbits := function(p)   #"MapOrbits"
 local i, j, res;
 res := List([1..Sum(p, Length)], i -> 0);
 for i in [1..Length(p)] do for j in p[i] do res[j] := i; od; od;
 return res;
end;
FactorSkeleton := function(S, p)  #"FactorSkeleton"
 local i, j, x, y, m, res;
#% assume "p" is a permutation given by a list
 if IsList(p) and (Maximum(p) = S.Count) then p := PermList(p); fi;
#% here, p IS a permutation
 if IsPerm(p) then p := OrbitsDomain(Group(p), [1..S.Count]); fi;
#% assume "p" a list of orbits
 if IsList(p) and IsList(p[1]) then p := MapOrbits(p); fi;
#% now "p" is a morphism given by a list
 m := Maximum(p); #%
 x := List([1..m], i -> 0); #%
 y := List([1..m], i -> 0);
 for i in [1..Length(p)] do
  x[p[i]] := p[i^S.X];#%
  y[p[i]] := p[i^S.Y];
 od;
 res := Skeleton(m, x, y, "parent", S,  #: Original skeleton and projection
                 "proj", l -> List(l, i -> p[i]));
 if IsMorphism(S, res, p) then return res; fi;
 return fail;
end;
#================================================================
#[endinput]
