1############################################################################## 2## 3## grape.g (Version 4.8.3) GRAPE Library Leonard Soicher 4## 5## Copyright (C) 1992-2019 Leonard Soicher, School of Mathematical Sciences, 6## Queen Mary University of London, London E1 4NS, U.K. 7## 8# This version includes code by Jerry James (debugged by LS) 9# which allows a user to use bliss instead of nauty for computing 10# automorphism groups and canonical labellings of graphs. 11# 12# This program is free software; you can redistribute it and/or modify 13# it under the terms of the GNU General Public License as published by 14# the Free Software Foundation; either version 2 of the License, or 15# (at your option) any later version. 16# 17# This program is distributed in the hope that it will be useful, 18# but WITHOUT ANY WARRANTY; without even the implied warranty of 19# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20# GNU General Public License for more details. 21# 22# You should have received a copy of the GNU General Public License 23# along with this program; if not, see http://www.gnu.org/licenses/gpl.html 24# 25 26GRAPE_RANDOM := false; # Determines if certain random methods are to be used 27 # in GRAPE functions. 28 # The default is that these random methods are not 29 # used (GRAPE_RANDOM=false). 30 # If these random methods are used (GRAPE_RANDOM=true), 31 # this does not affect the correctness of the results, 32 # as documented, returned by GRAPE functions, but may 33 # influence the time taken and the actual (correct) 34 # values returned. Due to improvements 35 # in GRAPE and in permutation group methods in GAP4, 36 # the use of random methods is rarely necessary, 37 # and should only be employed by GRAPE experts. 38 39GRAPE_NRANGENS := 18; # The number of random generators taken for a subgroup 40 # when GRAPE_RANDOM=true. 41 42GRAPE_NAUTY := true; # Use nauty when true, else use bliss. 43 44GRAPE_DREADNAUT_EXE := 45 ExternalFilename(DirectoriesPackagePrograms("grape"),"dreadnautB"); 46 # filename of dreadnaut or dreadnautB executable 47 48GRAPE_BLISS_EXE := ExternalFilename(DirectoriesSystemPrograms(),"bliss"); 49 # filename of bliss executable 50 51# The following variant of GAP's Exec is more flexible, and does not require a 52# shell. That makes it more reliable on Windows resp. with Cygwin. Moreover, 53# it allows to redirect input and output. 54BindGlobal("GRAPE_Exec", function(cmd, args, instream, outstream) 55 local dir, status; 56 57 if not IsString(cmd) then 58 Error("<cmd> must be a file path"); 59 fi; 60 61 if not IsInputStream(instream) then 62 Error("<instream> must be an input stream"); 63 fi; 64 65 if not IsOutputStream(outstream) then 66 Error("<outstream> must be an output stream"); 67 fi; 68 69 # execute in the current directory 70 dir := DirectoryCurrent(); 71 72 # execute the command 73 status := Process(dir, cmd, instream, outstream, args); 74 75 return status; 76end); 77 78BindGlobal("GRAPE_OrbitNumbers",function(G,n) 79# 80# Returns the orbits of G on [1..n] in the form of a record 81# containing the orbit representatives and a length n list 82# orbitNumbers, such that orbitNumbers[j]=i means that point 83# j is in the orbit of the i-th representative. 84# 85local i,j,orbnum,reps,im,norb,g,orb; 86if not IsPermGroup(G) or not IsInt(n) then 87 Error("usage: GRAPE_OrbitNumbers( <PermGroup>, <Int> )"); 88fi; 89orbnum:=[]; 90for i in [1..n] do 91 orbnum[i]:=0; 92od; 93reps:=[]; 94norb:=0; 95for i in [1..n] do 96 if orbnum[i]=0 then # new orbit 97 Add(reps,i); 98 norb:=norb+1; 99 orbnum[i]:=norb; 100 orb:=[i]; 101 for j in orb do 102 for g in GeneratorsOfGroup(G) do 103 im:=j^g; 104 if orbnum[im]=0 then 105 orbnum[im]:=norb; 106 Add(orb,im); 107 fi; 108 od; 109 od; 110 fi; 111od; 112return rec(representatives:=reps,orbitNumbers:=orbnum); 113end); 114 115BindGlobal("GRAPE_NumbersToSets",function(vec) 116# 117# Returns a list of sets as described by the numbers in vec, i.e. 118# i is in the j-th set iff vec[i]=j>0. 119# 120local list,i,j; 121if not IsList(vec) then 122 Error("usage: GRAPE_NumbersToSets( <List> )"); 123fi; 124if Length(vec)=0 then 125 return []; 126fi; 127list:=[]; 128for i in [1..Maximum(vec)] do 129 list[i]:=[]; 130od; 131for i in [1..Length(vec)] do 132 j:=vec[i]; 133 if j>0 then 134 Add(list[j],i); 135 fi; 136od; 137for i in [1..Length(list)] do 138 IsSSortedList(list[i]); 139od; 140return list; 141end); 142 143BindGlobal("GRAPE_IntransitiveGroupGenerators",function(arg) 144local conjperm,i,newgens,gens1,gens2,max1,max2; 145gens1:=arg[1]; 146gens2:=arg[2]; 147if IsBound(arg[3]) then 148 max1:=arg[3]; 149else 150 max1:=LargestMovedPoint(gens1); 151fi; 152if IsBound(arg[4]) then 153 max2:=arg[4]; 154else 155 max2:=LargestMovedPoint(gens2); 156fi; 157if not (IsList(gens1) and IsList(gens2) and IsInt(max1) and IsInt(max2)) then 158 Error( 159 "usage: GRAPE_IntransitiveGroupGenerators( <List>, <List> [,<Int> [,<Int> ]] )"); 160fi; 161if Length(gens1)<>Length(gens2) then 162 Error("Length(<gens1>) <> Length(<gens2>)"); 163fi; 164conjperm:=PermList(Concatenation(List([1..max2],x->x+max1),[1..max1])); 165newgens:=[]; 166for i in [1..Length(gens1)] do 167 newgens[i]:=gens1[i]*(gens2[i]^conjperm); 168od; 169return newgens; 170end); 171 172BindGlobal("ProbablyStabilizer",function(G,pt) 173# 174# Returns a subgroup of Stabilizer(G,pt), which is very often 175# the full stabilizer. In fact, if GRAPE_RANDOM=false, then it 176# is guaranteed to be the full stabilizer. 177# 178local sch,orb,x,y,im,k,gens,g,stabgens,i; 179if not IsPermGroup(G) or not IsInt(pt) then 180 Error("usage: ProbablyStabilizer( <PermGroup>, <Int> )"); 181fi; 182if not GRAPE_RANDOM or HasSize(G) or HasStabChainMutable(G) or IsAbelian(G) then 183 return Stabilizer(G,pt); 184fi; 185# 186# At this point we know that G is non-abelian. In particular, 187# G has at least two generators. 188# 189# First, we make a Schreier vector of permutations for the orbit pt^G. 190# 191gens:=GeneratorsOfGroup(G); 192sch:=[]; 193orb:=[pt]; 194sch[pt]:=(); 195for x in orb do 196 for g in gens do 197 im:=x^g; 198 if not IsBound(sch[im]) then 199 sch[im]:=g; 200 Add(orb,im); 201 fi; 202 od; 203od; 204# Now make gens into a new randomish generating sequence for G. 205gens:=ShallowCopy(GeneratorsOfGroup(G)); 206k:=Length(gens); 207for i in [k+1..GRAPE_NRANGENS] do 208 gens[i]:=gens[Random([1..k])]; 209od; 210for i in [1..Length(gens)] do 211 gens[i]:=Product(gens); 212od; 213# Now make a list stabgens of random elements of the stabilizer of pt. 214x:=(); 215stabgens:=[]; 216for i in [1..GRAPE_NRANGENS] do 217 x:=x*Random(gens); 218 im:=pt^x; 219 while im<>pt do 220 x:=x/sch[im]; 221 im:=im/sch[im]; 222 od; 223 if x<>() then 224 Add(stabgens,x); 225 fi; 226od; 227return Group(stabgens,()); 228end); 229 230BindGlobal("ProbablyStabilizerOrbitNumbers",function(G,pt,n) 231# 232# Returns the "orbit numbers" record for a subgroup of Stabilizer(G,pt), 233# in its action on [1..n]. 234# This subgroup is very often the full stabilizer, and in fact, 235# if GRAPE_RANDOM=false, then it is guaranteed to be the full stabilizer. 236# 237if not IsPermGroup(G) or not IsInt(pt) or not IsInt(n) then 238 Error( 239 "usage: ProbablyStabilizerOrbitNumbers( <PermGroup>, <Int>, <Int> )"); 240fi; 241return GRAPE_OrbitNumbers(ProbablyStabilizer(G,pt),n); 242end); 243 244BindGlobal("GRAPE_RepWord",function(gens,sch,r) 245# 246# Given a sequence gens of group generators, and a (word type) 247# schreier vector sch made using gens, this function returns a 248# record containing the orbit representative for r (wrt sch), and 249# a word in gens taking this representative to r. 250# (We assume sch includes the orbit of r.) 251# 252local word,w; 253word:=[]; 254w:=sch[r]; 255while w > 0 do 256 Add(word,w); 257 r:=r/gens[w]; 258 w:=sch[r]; 259od; 260return rec(word:=Reversed(word),representative:=r); 261end); 262 263BindGlobal("NullGraph",function(arg) 264# 265# Returns a null graph with n vertices and group G=arg[1]. 266# If arg[2] is bound then n=arg[2], otherwise n is the maximum 267# largest moved point of the generators of G. 268# The names, autGroup, maximumClique, minimumVertexColouring, 269# and canonicalLabelling components of the 270# returned null graph are left unbound; however, the isSimple 271# component is set (to true). 272# 273local G,n,gamma,nadj,sch,orb,i,j,k,im,gens; 274G:=arg[1]; 275if not IsPermGroup(G) or (IsBound(arg[2]) and not IsInt(arg[2])) then 276 Error("usage: NullGraph( <PermGroup>, [, <Int> ] )"); 277fi; 278n:=LargestMovedPoint(GeneratorsOfGroup(G)); 279if IsBound(arg[2]) then 280 if arg[2] < n then 281 Error("<arg[2]> too small"); 282 fi; 283 n:=arg[2]; 284fi; 285gamma:=rec(isGraph:=true,order:=n,group:=G,schreierVector:=[], 286 adjacencies:=[],representatives:=[],isSimple:=true); 287# 288# Calculate gamma.representatives, gamma.schreierVector, and 289# gamma.adjacencies. 290# 291sch:=gamma.schreierVector; 292gens:=GeneratorsOfGroup(gamma.group); 293nadj:=0; 294for i in [1..n] do 295 sch[i]:=0; 296od; 297for i in [1..n] do 298 if sch[i]=0 then # new orbit 299 Add(gamma.representatives,i); 300 nadj:=nadj+1; 301 sch[i]:=-nadj; # tells where to find the adjacency set. 302 gamma.adjacencies[nadj]:=[]; 303 orb:=[i]; 304 for j in orb do 305 for k in [1..Length(gens)] do 306 im:=j^gens[k]; 307 if sch[im]=0 then 308 sch[im]:=k; 309 Add(orb,im); 310 fi; 311 od; 312 od; 313 fi; 314od; 315gamma.representatives:=Immutable(gamma.representatives); 316gamma.schreierVector:=Immutable(gamma.schreierVector); 317return gamma; 318end); 319 320BindGlobal("CompleteGraph",function(arg) 321# 322# Returns a complete graph with n vertices and group G=arg[1]. 323# If arg[2] is bound then n=arg[2], otherwise n is the maximum 324# largest moved point of the generators of the permutation group G. 325# If the boolean argument arg[3] is bound and has 326# value true then the complete graph will have all possible loops, 327# otherwise it will have no loops (the default). 328# 329# The names, autGroup, maximumClique, minimumVertexColouring, 330# and canonicalLabelling components of the 331# returned complete graph are left unbound; however, the isSimple 332# component is set (appropriately). 333# 334local G,n,gamma,i,mustloops; 335G:=arg[1]; 336if not IsPermGroup(G) or (IsBound(arg[2]) and not IsInt(arg[2])) 337 or (IsBound(arg[3]) and not IsBool(arg[3])) then 338 Error("usage: CompleteGraph( <PermGroup>, [, <Int> [, <Bool> ]] )"); 339fi; 340n:=LargestMovedPoint(GeneratorsOfGroup(G)); 341if IsBound(arg[2]) then 342 if arg[2] < n then 343 Error("<arg[2]> too small"); 344 fi; 345 n:=arg[2]; 346fi; 347if IsBound(arg[3]) then 348 mustloops:=arg[3]; 349else 350 mustloops:=false; 351fi; 352gamma:=NullGraph(G,n); 353if gamma.order=0 then 354 return gamma; 355fi; 356if mustloops then 357 gamma.isSimple:=false; 358fi; 359for i in [1..Length(gamma.adjacencies)] do 360 gamma.adjacencies[i]:=[1..n]; 361 if not mustloops then 362 RemoveSet(gamma.adjacencies[i],gamma.representatives[i]); 363 fi; 364od; 365return gamma; 366end); 367 368BindGlobal("GRAPE_Graph",function(arg) 369# 370# First suppose that arg[5] is unbound or has value false. 371# Then L=arg[2] is a list of elements of a set S on which 372# G=arg[1] acts with action act=arg[3]. Also rel=arg[4] is a boolean 373# function defining a G-invariant relation on S (so that 374# for g in G, rel(x,y) iff rel(act(x,g),act(y,g)) ). 375# Then function GRAPE_Graph returns the graph gamma with vertex-names 376# Immutable(Concatenation(Orbits(G,L,act))), and x is joined to y 377# in gamma iff rel(VertexName(gamma,x),VertexName(gamma,y)). 378# 379# If arg[5] has value true then it is assumed that L=arg[2] 380# is invariant under G=arg[1] with action act=arg[3]. Then 381# the function GRAPE_Graph behaves as above, except that gamma.names 382# becomes an immutable copy of L. 383# 384local G,L,act,rel,invt,gamma,vertexnames,i,reps,H,orb,x,y,adj; 385G:=arg[1]; 386L:=arg[2]; 387act:=arg[3]; 388rel:=arg[4]; 389if IsBound(arg[5]) then 390 invt:=arg[5]; 391else 392 invt:=false; 393fi; 394if not (IsGroup(G) and IsList(L) and IsFunction(act) and IsFunction(rel) 395 and IsBool(invt)) then 396 Error("usage: GRAPE_Graph( <Group>, <List>, <Function>, <Function> [, <Bool> ] )"); 397fi; 398if invt then 399 vertexnames:=Immutable(L); 400else 401 vertexnames:=Immutable(Concatenation(Orbits(G,L,act))); 402fi; 403gamma:=NullGraph(Action(G,vertexnames,act),Length(vertexnames)); 404Unbind(gamma.isSimple); 405gamma.names:=vertexnames; 406if not GRAPE_RANDOM then 407 if (HasSize(G) and Size(G)<>infinity) or 408 (IsPermGroup(G) and HasStabChainMutable(G)) or 409 (HasIsNaturalSymmetricGroup(G) and IsNaturalSymmetricGroup(G)) then 410 StabChainOp(gamma.group,rec(limit:=Size(G))); 411 fi; 412fi; 413reps:=gamma.representatives; 414for i in [1..Length(reps)] do 415 H:=ProbablyStabilizer(gamma.group,reps[i]); 416 x:=vertexnames[reps[i]]; 417 if IsTrivial(H) then 418 gamma.adjacencies[i]:=Filtered([1..gamma.order],j->rel(x,vertexnames[j])); 419 else 420 adj:=[]; 421 for orb in OrbitsDomain(H,[1..gamma.order]) do 422 y:=vertexnames[orb[1]]; 423 if rel(x,y) then 424 Append(adj,orb); 425 fi; 426 od; 427 Sort(adj); 428 gamma.adjacencies[i]:=adj; 429 fi; 430 IsSSortedList(gamma.adjacencies[i]); 431od; 432return gamma; 433end); 434 435DeclareOperation("Graph",[IsGroup,IsList,IsFunction,IsFunction]); 436InstallMethod(Graph,"for use in GRAPE with 4 parameters", 437 [IsGroup,IsList,IsFunction,IsFunction],0,GRAPE_Graph); 438DeclareOperation("Graph",[IsGroup,IsList,IsFunction,IsFunction,IsBool]); 439InstallMethod(Graph,"for use in GRAPE with 5 parameters", 440 [IsGroup,IsList,IsFunction,IsFunction,IsBool],0,GRAPE_Graph); 441 442BindGlobal("JohnsonGraph",function(n,e) 443# 444# Returns the Johnson graph, whose vertices are the e-subsets 445# of {1,...,n}, with x joined to y iff Intersection(x,y) 446# has size e-1. 447# 448local rel,J; 449if not IsInt(n) or not IsInt(e) then 450 Error("usage: JohnsonGraph( <Int>, <Int> )"); 451fi; 452if e<0 or n<e then 453 Error("must have 0 <= <e> <= <n>"); 454fi; 455rel := function(x,y) 456 return Length(Intersection(x,y))=e-1; 457end; 458J:=Graph(SymmetricGroup(n),Combinations([1..n],e),OnSets,rel,true); 459J.isSimple:=true; 460return J; 461end); 462 463BindGlobal("HammingGraph",function(d,q) 464# 465# Where d and q are positive integers, this function returns the 466# Hamming graph H(d,q), defined as follows. The set of vertices 467# (actually vertex-names) is the set of all d-tuples of elements 468# of [1..q], with vertices v and w joined by an edge iff their 469# Hamming distance is 1. The group associated with the returned 470# graph is S_q wr S_d in its product action on the vertices. 471# 472local W,projection,embedding,moved,act,rel; 473if not IsPosInt(d) or not IsPosInt(q) then 474 Error("usage: HammingGraph( <PosInt>, <PosInt> )"); 475fi; 476if q=1 then 477 # special trivial case 478 return Graph(Group(()),[ListWithIdenticalEntries(d,1)], 479 function(x,g) return x; end,function(x,y) return false; end,true); 480fi; 481W:=WreathProductImprimitiveAction(SymmetricGroup([1..q]),SymmetricGroup([1..d])); 482projection:=Projection(W); 483embedding:=Embedding(W,d+1); 484moved:=List([1..d],i->MovedPoints(Image(Embedding(W,i)))); 485act := function(x,g) 486# Product action of g on d-tuple x. 487local bb,b,a,y,i; 488bb:=g^projection; 489b:=bb^embedding; 490a:=g*b^(-1); # so g factorises as a*b in W 491y:=[]; 492for i in [1..d] do 493 y[i^bb]:=PositionSorted(moved[i],moved[i][x[i]]^a); 494od; 495return y; 496end; 497rel := function(x,y) 498# boolean function returning true iff the d-tuples x and y 499# are at Hamming distance 1. 500local i,count; 501count:=0; 502for i in [1..d] do 503 if x[i]<>y[i] then 504 count:=count+1; 505 if count>1 then 506 return false; 507 fi; 508 fi; 509od; 510return count=1; 511end; 512# Now contruct and return the Hamming graph. 513return Graph(W,Tuples([1..q],d),act,rel,true); 514end); 515 516BindGlobal("IsGraph",function(obj) 517# 518# Returns true iff obj is a (GRAPE) graph. 519# 520return IsRecord(obj) and IsBound(obj.isGraph) and obj.isGraph=true 521 and IsBound(obj.group) and IsBound(obj.schreierVector); 522end); 523 524BindGlobal("CopyGraph",function(gamma) 525# 526# Returns a "structural" copy delta of the graph gamma, and 527# also ensures that the appropriate components of delta are immutable. 528# 529local delta; 530if not IsGraph(gamma) then 531 Error("usage: CopyGraph( <Graph> )"); 532fi; 533delta:=ShallowCopy(gamma); 534delta.adjacencies:=StructuralCopy(delta.adjacencies); 535delta.representatives:=Immutable(delta.representatives); 536delta.schreierVector:=Immutable(delta.schreierVector); 537if IsBound(delta.names) then 538 delta.names:=Immutable(delta.names); 539fi; 540if IsBound(delta.maximumClique) then 541 delta.maximumClique:=Immutable(delta.maximumClique); 542fi; 543if IsBound(delta.minimumVertexColouring) then 544 delta.minimumVertexColouring:=Immutable(delta.minimumVertexColouring); 545fi; 546Unbind(delta.canonicalLabelling); # for safety 547return delta; 548end); 549 550BindGlobal("OrderGraph",function(gamma) 551# 552# returns the order of gamma. 553# 554if not IsGraph(gamma) then 555 Error("usage: OrderGraph( <Graph> )"); 556fi; 557return gamma.order; 558end); 559 560DeclareOperation("Vertices",[IsRecord]); 561# to avoid the clash with `Vertices' defined in the xgap package 562InstallMethod(Vertices,"for GRAPE graph",[IsRecord],0, 563function(gamma) 564# 565# Returns the vertex-set of graph gamma. 566# 567if not IsGraph(gamma) then 568 TryNextMethod(); 569fi; 570return [1..gamma.order]; 571end); 572 573DeclareOperation("IsVertex",[IsRecord,IsObject]); 574InstallMethod(IsVertex,"for GRAPE graph",[IsRecord,IsObject],0, 575function(gamma,v) 576# 577# Returns true iff v is vertex of gamma. 578# 579if not IsGraph(gamma) then 580 TryNextMethod(); 581fi; 582return IsInt(v) and v >= 1 and v <= gamma.order; 583end); 584 585BindGlobal("AssignVertexNames",function(gamma,names) 586# 587# Assign vertex names for gamma, so that the (external) name of 588# vertex i becomes names[i], by making gamma.names an immutable 589# copy of names. 590# 591if not IsGraph(gamma) or not IsList(names) then 592 Error("usage: AssignVertexNames( <Graph>, <List> )"); 593fi; 594if Length(names)<>gamma.order then 595 Error("Length(<names>) <> gamma.order"); 596fi; 597gamma.names:=Immutable(names); 598end); 599 600BindGlobal("VertexName",function(gamma,v) 601# 602# Returns the (external) name of the vertex v of gamma. 603# 604if not IsGraph(gamma) or not IsInt(v) then 605 Error("usage: VertexName( <Graph>, <Int> )"); 606fi; 607if IsBound(gamma.names) then 608 return Immutable(gamma.names[v]); 609else 610 return v; 611fi; 612end); 613 614BindGlobal("VertexNames",function(gamma) 615# 616# Returns the list of (external) names of the vertices of gamma. 617# 618if not IsGraph(gamma) then 619 Error("usage: VertexNames( <Graph> )"); 620fi; 621if IsBound(gamma.names) then 622 return Immutable(gamma.names); 623else 624 return Immutable([1..gamma.order]); 625fi; 626end); 627 628BindGlobal("VertexDegree",function(gamma,v) 629# 630# Returns the vertex (out)degree of vertex v in the graph gamma. 631# 632local rw,sch; 633if not IsGraph(gamma) or not IsInt(v) then 634 Error("usage: VertexDegree( <Graph>, <Int> )"); 635fi; 636if v<1 or v>gamma.order then 637 Error("<v> is not a vertex of <gamma>"); 638fi; 639sch:=gamma.schreierVector; 640rw:=GRAPE_RepWord(GeneratorsOfGroup(gamma.group),sch,v); 641return Length(gamma.adjacencies[-sch[rw.representative]]); 642end); 643 644BindGlobal("VertexDegrees",function(gamma) 645# 646# Returns the set of vertex (out)degrees for the graph gamma. 647# 648local adj,degs; 649if not IsGraph(gamma) then 650 Error("usage: VertexDegrees( <Graph> )"); 651fi; 652degs:=[]; 653for adj in gamma.adjacencies do 654 AddSet(degs,Length(adj)); 655od; 656return degs; 657end); 658 659BindGlobal("IsVertexPairEdge",function(gamma,x,y) 660# 661# Assuming that x,y are vertices of gamma, returns true 662# iff [x,y] is an edge of gamma. 663# 664local w,sch,gens; 665sch:=gamma.schreierVector; 666gens:=GeneratorsOfGroup(gamma.group); 667w:=sch[x]; 668while w > 0 do 669 x:=x/gens[w]; 670 y:=y/gens[w]; 671 w:=sch[x]; 672od; 673return y in gamma.adjacencies[-w]; 674end); 675 676DeclareOperation("IsEdge",[IsRecord,IsObject]); 677InstallMethod(IsEdge,"for GRAPE graph",[IsRecord,IsObject],0, 678function(gamma,e) 679# 680# Returns true iff e is an edge of gamma. 681# 682if not IsGraph(gamma) then 683 TryNextMethod(); 684fi; 685if not IsList(e) or Length(e)<>2 or not IsVertex(gamma,e[1]) 686 or not IsVertex(gamma,e[2]) then 687 return false; 688fi; 689return IsVertexPairEdge(gamma,e[1],e[2]); 690end); 691 692BindGlobal("Adjacency",function(gamma,v) 693# 694# Returns (a copy of) the set of vertices of gamma adjacent to vertex v. 695# 696local w,adj,rw,gens,sch; 697sch:=gamma.schreierVector; 698if sch[v] < 0 then 699 return ShallowCopy(gamma.adjacencies[-sch[v]]); 700fi; 701gens:=GeneratorsOfGroup(gamma.group); 702rw:=GRAPE_RepWord(gens,sch,v); 703adj:=gamma.adjacencies[-sch[rw.representative]]; 704for w in rw.word do 705 adj:=OnTuples(adj,gens[w]); 706od; 707return SSortedList(adj); 708end); 709 710BindGlobal("IsSimpleGraph",function(gamma) 711# 712# Returns true iff graph gamma is simple (i.e. has no loops and 713# if [x,y] is an edge then so is [y,x]). Also sets the isSimple 714# field of gamma if this field was not already bound. 715# 716local adj,i,x,H,orb; 717if not IsGraph(gamma) then 718 Error("usage: IsSimpleGraph( <Graph> )"); 719fi; 720if IsBound(gamma.isSimple) then 721 return gamma.isSimple; 722fi; 723for i in [1..Length(gamma.adjacencies)] do 724 adj:=gamma.adjacencies[i]; 725 x:=gamma.representatives[i]; 726 if x in adj then # a loop exists 727 gamma.isSimple:=false; 728 return false; 729 fi; 730 H:=ProbablyStabilizer(gamma.group,x); 731 for orb in OrbitsDomain(H,adj) do 732 if not IsVertexPairEdge(gamma,orb[1],x) then 733 gamma.isSimple:=false; 734 return false; 735 fi; 736 od; 737od; 738gamma.isSimple:=true; 739return true; 740end); 741 742BindGlobal("DirectedEdges",function(gamma) 743# 744# Returns the set of directed (ordered) edges of gamma. 745# 746local i,j,edges; 747if not IsGraph(gamma) then 748 Error("usage: DirectedEdges( <Graph> )"); 749fi; 750edges:=[]; 751for i in [1..gamma.order] do 752 for j in Adjacency(gamma,i) do 753 Add(edges,[i,j]); 754 od; 755od; 756IsSSortedList(edges); # edges is a set. 757return edges; 758end); 759 760BindGlobal("UndirectedEdges",function(gamma) 761# 762# Returns the set of undirected edges of gamma, which must be 763# a simple graph. 764# 765local i,j,edges; 766if not IsGraph(gamma) then 767 Error("usage: UndirectedEdges( <Graph> )"); 768fi; 769if not IsSimpleGraph(gamma) then 770 Error("<gamma> must be a simple graph"); 771fi; 772edges:=[]; 773for i in [1..gamma.order-1] do 774 for j in Adjacency(gamma,i) do 775 if i<j then 776 Add(edges,[i,j]); 777 fi; 778 od; 779od; 780IsSSortedList(edges); # edges is a set. 781return edges; 782end); 783 784BindGlobal("AddEdgeOrbit",function(arg) 785# 786# Let gamma=arg[1] and e=arg[2]. 787# If arg[3] is bound then it is assumed to be Stabilizer(gamma.group,e[1]). 788# This procedure adds edge orbit e^gamma.group to the edge-set of gamma. 789# 790local w,word,sch,gens,gamma,e,x,y,orb,u,v; 791gamma:=arg[1]; 792e:=arg[2]; 793if not IsGraph(gamma) or not IsList(e) 794 or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then 795 Error("usage: AddEdgeOrbit( <Graph>, <List>, [, <PermGroup> ] )"); 796fi; 797if Length(e)<>2 or not IsVertex(gamma,e[1]) or not IsVertex(gamma,e[2]) then 798 Error("invalid <e>"); 799fi; 800sch:=gamma.schreierVector; 801gens:=GeneratorsOfGroup(gamma.group); 802x:=e[1]; 803y:=e[2]; 804w:=sch[x]; 805word:=[]; 806while w > 0 do 807 Add(word,w); 808 x:=x/gens[w]; 809 y:=y/gens[w]; 810 w:=sch[x]; 811od; 812if not(y in gamma.adjacencies[-sch[x]]) then 813 # e is not an edge of gamma 814 if not IsBound(arg[3]) then 815 orb:=Orbit(Stabilizer(gamma.group,x),y); 816 else 817 if ForAny(GeneratorsOfGroup(arg[3]),x->e[1]^x<>e[1]) then 818 Error("<arg[3]> not equal to Stabilizer(<gamma.group>,<e[1]>)"); 819 fi; 820 orb:=[]; 821 for u in Orbit(arg[3],e[2]) do 822 v:=u; 823 for w in word do 824 v:=v/gens[w]; 825 od; 826 Add(orb,v); 827 od; 828 fi; 829 UniteSet(gamma.adjacencies[-sch[x]],orb); 830 if e[1]=e[2] then 831 gamma.isSimple:=false; 832 elif IsBound(gamma.isSimple) and gamma.isSimple then 833 if not IsVertexPairEdge(gamma,e[2],e[1]) then 834 gamma.isSimple:=false; 835 fi; 836 else 837 Unbind(gamma.isSimple); 838 fi; 839 Unbind(gamma.autGroup); 840 Unbind(gamma.canonicalLabelling); 841 Unbind(gamma.maximumClique); 842 Unbind(gamma.minimumVertexColouring); 843fi; 844end); 845 846BindGlobal("RemoveEdgeOrbit",function(arg) 847# 848# Let gamma=arg[1] and e=arg[2]. 849# If arg[3] is bound then it is assumed to be Stabilizer(gamma.group,e[1]). 850# This procedure removes the edge orbit e^gamma.group from the edge-set 851# of gamma, if this orbit exists, and otherwise does nothing. 852# 853local w,word,sch,gens,gamma,e,x,y,orb,u,v; 854gamma:=arg[1]; 855e:=arg[2]; 856if not IsGraph(gamma) or not IsList(e) 857 or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then 858 Error("usage: RemoveEdgeOrbit( <Graph>, <List>, [, <PermGroup> ] )"); 859fi; 860if Length(e)<>2 or not IsVertex(gamma,e[1]) or not IsVertex(gamma,e[2]) then 861 Error("invalid <e>"); 862fi; 863sch:=gamma.schreierVector; 864gens:=GeneratorsOfGroup(gamma.group); 865x:=e[1]; 866y:=e[2]; 867w:=sch[x]; 868word:=[]; 869while w > 0 do 870 Add(word,w); 871 x:=x/gens[w]; 872 y:=y/gens[w]; 873 w:=sch[x]; 874od; 875if y in gamma.adjacencies[-sch[x]] then 876 # e is an edge of gamma 877 if not IsBound(arg[3]) then 878 orb:=Orbit(Stabilizer(gamma.group,x),y); 879 else 880 if ForAny(GeneratorsOfGroup(arg[3]),x->e[1]^x<>e[1]) then 881 Error("<arg[3]> not equal to Stabilizer(<gamma.group>,<e[1]>)"); 882 fi; 883 orb:=[]; 884 for u in Orbit(arg[3],e[2]) do 885 v:=u; 886 for w in word do 887 v:=v/gens[w]; 888 od; 889 Add(orb,v); 890 od; 891 fi; 892 SubtractSet(gamma.adjacencies[-sch[x]],orb); 893 if IsBound(gamma.isSimple) and gamma.isSimple then 894 if IsVertexPairEdge(gamma,e[2],e[1]) then 895 gamma.isSimple:=false; 896 fi; 897 else 898 Unbind(gamma.isSimple); 899 fi; 900 Unbind(gamma.autGroup); 901 Unbind(gamma.canonicalLabelling); 902 Unbind(gamma.maximumClique); 903 Unbind(gamma.minimumVertexColouring); 904fi; 905end); 906 907BindGlobal("EdgeOrbitsGraph",function(arg) 908# 909# Let G=arg[1], E=arg[2]. 910# Returns the (directed) graph with vertex-set {1,...,n} and edge-set 911# the union over e in E of e^G, where n=arg[3] if arg[3] is bound, 912# and n=LargestMovedPoint(GeneratorsOfGroup(G)) otherwise. 913# (E can consist of just a singleton edge.) 914# 915local G,E,n,gamma,e; 916G:=arg[1]; 917E:=arg[2]; 918if IsInt(E[1]) then # assume E consists of a single edge. 919 E:=[E]; 920fi; 921if IsBound(arg[3]) then 922 n:=arg[3]; 923else 924 n:=LargestMovedPoint(GeneratorsOfGroup(G)); 925fi; 926if not IsPermGroup(G) or not IsList(E) or not IsInt(n) then 927 Error("usage: EdgeOrbitsGraph( <PermGroup>, <List> [, <Int> ] )"); 928fi; 929gamma:=NullGraph(G,n); 930for e in E do 931 AddEdgeOrbit(gamma,e); 932od; 933return gamma; 934end); 935 936BindGlobal("GeneralizedOrbitalGraphs",function(G) 937# 938# The input to this function is a non-trivial permutation group G, 939# acting transitively on [1..n] (where n:=LargestMovedPoint(G)). 940# 941# Then this function returns a list of all the generalized orbital 942# graphs for G (that is, the simple graphs with vertex set [1..n] 943# and edge-set a union of orbitals of G). 944# 945local comb,n,H,result,reps,i,L,M,mm; 946if not IsPermGroup(G) then 947 Error("usage: GeneralizedOrbitalGraphs( <PermGroup> )"); 948fi; 949n:=LargestMovedPoint(G); 950if n=0 or not IsTransitive(G,[1..n]) then 951 Error("<G> must be a non-trivial transitive group on [1..LargestMovedPoint( <G> )"); 952fi; 953H:=Stabilizer(G,1); 954reps:=Set(List(OrbitsDomain(H,[2..n]),Minimum)); 955# 956# Now make a duplicate-free list L of the graphs 957# with vertex-set [1..n] and edge-set the union 958# of a nondiagonal G-orbital and its paired orbital. 959# At the same time, make a list M whose i-th element is 960# a list of edges of L[i], such that the union of the 961# G-orbits of the edges in M[i] is the edge-set of L[i]. 962# 963L:=[]; 964M:=[]; 965for i in [1..Length(reps)] do 966 if ForAll(L,x->not IsEdge(x,[1,reps[i]])) then 967 mm:=[[1,reps[i]],[reps[i],1]]; 968 Add(L,EdgeOrbitsGraph(G,mm)); 969 Add(M,mm); 970 fi; 971od; 972result:=[]; 973for comb in Combinations(M) do 974 if comb<>[] then 975 Add(result,EdgeOrbitsGraph(G,Concatenation(comb))); 976 fi; 977od; 978return result; 979end); 980 981BindGlobal("CollapsedAdjacencyMat",function(arg) 982# 983# Returns the collapsed adjacency matrix A for gamma=arg[2] wrt 984# group G=arg[1], assuming G <= Aut(gamma). 985# The rows and columns of A are indexed by the orbits 986# orbs[1],...,orbs[n], say, of G on the vertices of 987# gamma, and the entry A[i][j] of A is defined as follows: 988# Let reps[i] be a representative of the i-th G-orbit orbs[i]. 989# Then A[i][j] equals the number of neighbours (in gamma) 990# of reps[i] in orbs[j]. 991# Note that this definition does not depend on the choice of 992# representative reps[i]. 993# 994# *** New for Grape 2.3: In the special case where this function 995# is given just one argument, then we must have gamma=arg[1], 996# we must have gamma.group transitive on the vertices of gamma, 997# and then the returned collapsed adjacency matrix for gamma is 998# w.r.t. the stabilizer in gamma.group of 1. Additionally 999# [1]=orbs[1]. This feature is to 1000# conform with the definition of collapsed adjacency matrix in 1001# Praeger and Soicher, "Low Rank Representations and Graphs for 1002# Sporadic Groups", CUP, Cambridge, 1997. (In GRAPE we allow a collapsed 1003# adjacency matrix to be more general, as we can collapse w.r.t. to 1004# an arbitrary subgroup of Aut(gamma), and gamma need not 1005# even be vertex-transitive.) 1006# 1007local G,gamma,orbs,i,j,n,A,orbnum,reps; 1008if Length(arg)=1 then 1009 gamma:=arg[1]; 1010 if not IsGraph(gamma) then 1011 Error("usage: CollapsedAdjacencyMat( [<PermGroup>,] <Graph>)"); 1012 fi; 1013 if gamma.order=0 then 1014 return []; 1015 fi; 1016 if not IsTransitive(gamma.group,[1..gamma.order]) then 1017 Error( 1018 "<gamma.group> not transitive on vertices of single argument <gamma>"); 1019 fi; 1020 G := Stabilizer(gamma.group,1); 1021else 1022 G := arg[1]; 1023 gamma := arg[2]; 1024 if not IsPermGroup(G) or not IsGraph(gamma) then 1025 Error("usage: CollapsedAdjacencyMat( [<PermGroup>,] <Graph> )"); 1026 fi; 1027 if gamma.order=0 then 1028 return []; 1029 fi; 1030fi; 1031orbs:=GRAPE_OrbitNumbers(G,gamma.order); 1032orbnum:=orbs.orbitNumbers; 1033reps:=orbs.representatives; 1034n:=Length(reps); 1035A:=NullMat(n,n); 1036for i in [1..n] do 1037 for j in Adjacency(gamma,reps[i]) do 1038 A[i][orbnum[j]]:=A[i][orbnum[j]]+1; 1039 od; 1040od; 1041return A; 1042end); 1043 1044BindGlobal("OrbitalGraphColadjMats",function(arg) 1045# 1046# This function returns a sequence of collapsed adjacency 1047# matrices for the the orbital graphs of the (assumed) transitive 1048# G=arg[1]. The matrices are collapsed w.r.t. Stabilizer(G,1), so 1049# that these are collapsed adjacency matrices in the sense of 1050# Praeger and Soicher, "Low Rank Representations and Graphs for 1051# Sporadic Groups", CUP, Cambridge, 1997. 1052# The matrices are collapsed w.r.t. a fixed ordering of the G-suborbits, 1053# with the trivial suborbit [1] coming first. 1054# 1055# If arg[2] is bound, then it is assumed to be Stabilizer(G,1). 1056# 1057local G,H,orbs,deg,i,j,k,n,intmats,A,orbnum,reps,gamma; 1058G:=arg[1]; 1059if not IsPermGroup(G) or (IsBound(arg[2]) and not IsPermGroup(arg[2])) then 1060 Error("usage: OrbitalGraphColadjMats( <PermGroup> [, <PermGroup> ] )"); 1061fi; 1062if IsBound(arg[2]) then 1063 H:=arg[2]; 1064 if ForAny(GeneratorsOfGroup(H),x->1^x<>1) then 1065 Error("<H> does not fix the point 1"); 1066 fi; 1067else 1068 H:=Stabilizer(G,1); 1069fi; 1070deg:=Maximum(LargestMovedPoint(GeneratorsOfGroup(G)),1); 1071if not IsTransitive(G,[1..deg]) then 1072 Error("<G> not transitive"); 1073fi; 1074gamma:=NullGraph(G,deg); 1075orbs:=GRAPE_OrbitNumbers(H,gamma.order); 1076orbnum:=orbs.orbitNumbers; 1077reps:=orbs.representatives; 1078if reps[1]<>1 then # this cannot happen! 1079 Error("internal error"); 1080fi; 1081n:=Length(reps); 1082intmats:=[]; 1083for i in [1..n] do 1084 AddEdgeOrbit(gamma,[1,reps[i]],H); 1085 A:=NullMat(n,n); 1086 for j in [1..n] do 1087 for k in Adjacency(gamma,reps[j]) do 1088 A[j][orbnum[k]]:=A[j][orbnum[k]]+1; 1089 od; 1090 od; 1091 intmats[i]:=A; 1092 if i < n then 1093 RemoveEdgeOrbit(gamma,[1,reps[i]],H); 1094 fi; 1095od; 1096return intmats; 1097end); 1098 1099BindGlobal("LocalInfo",function(arg) 1100# 1101# Calculates "local info" for gamma=arg[1] from point of view of vertex 1102# set (or list or singleton vertex) V=arg[2]. 1103# 1104# Returns record containing the "layer numbers" for gamma w.r.t. 1105# V, as well as the the local diameter and local girth of gamma w.r.t. V. 1106# ( layerNumbers[i]=j>0 if vertex i is in layer[j] (i.e. at distance 1107# j-1 from V), layerNumbers[i]=0 if vertex i 1108# is not joined by a (directed) path from some element of V.) 1109# Also, if a local parameter ci[V], ai[V], or bi[V] exists then 1110# this information is recorded in localParameters[i+1][1], 1111# localParameters[i+1][2], or localParameters[i+1][3], 1112# respectively (otherwise a -1 is recorded). 1113# 1114# *** If gamma is not simple then local girth and the local parameters 1115# may not be what you think. The local girth has no real meaning if 1116# |V| > 1. 1117# 1118# *** But note: If arg[3] is bound and arg[3] > 0 1119# then the procedure stops after layers[arg[3]] has been determined. 1120# 1121# If arg[4] is bound (a set or list or singleton vertex), then the 1122# procedure stops when the first layer containing a vertex in arg[4] is 1123# complete. Moreover, if arg[4] is bound then the local info record 1124# contains a distance field, whose value (if arg[3]=0) is 1125# min{ d(v,w) | v in V, w in arg[4] }. 1126# 1127# If arg[5] is bound then it is assumed to be a subgroup of Aut(gamma) 1128# stabilising V setwise. 1129# 1130local gamma,V,layers,localDiameter,localGirth,localParameters,i,j,x,y,next, 1131 nprev,nhere,nnext,sum,orbs,orbnum,laynum,lnum, 1132 stoplayer,stopvertices,distance,loc,reps,layerNumbers; 1133gamma:=arg[1]; 1134V:=arg[2]; 1135if IsInt(V) then 1136 V:=[V]; 1137fi; 1138if not IsGraph(gamma) or not IsList(V) then 1139 Error("usage: LocalInfo( <Graph>, <Int> or <List>, ... )"); 1140fi; 1141if not IsSSortedList(V) then 1142 V:=SSortedList(V); 1143fi; 1144if V=[] or not IsSubset([1..gamma.order],V) then 1145 Error("<V> must be non-empty set of vertices of <gamma>"); 1146fi; 1147if IsBound(arg[3]) then 1148 stoplayer:=arg[3]; 1149 if not IsInt(stoplayer) or stoplayer < 0 then 1150 Error("<stoplayer> must be integer >= 0"); 1151 fi; 1152else 1153 stoplayer:=0; 1154fi; 1155if IsBound(arg[4]) then 1156 stopvertices:=arg[4]; 1157 if IsInt(stopvertices) then 1158 stopvertices:=[stopvertices]; 1159 fi; 1160 if not IsSSortedList(stopvertices) then 1161 stopvertices:=SSortedList(stopvertices); 1162 fi; 1163 if not IsSubset([1..gamma.order],stopvertices) 1164 then 1165 Error("<stopvertices> must be a set of vertices of <gamma>"); 1166 fi; 1167else 1168 stopvertices:=[]; 1169fi; 1170if IsBound(arg[5]) then 1171 if not IsPermGroup(arg[5]) then 1172 Error("<arg[5]> must be a permutation group (<= Stab(<V>)"); 1173 fi; 1174 orbs:=GRAPE_OrbitNumbers(arg[5],gamma.order); 1175else 1176 if Length(V)=1 then 1177 if IsBound(gamma.autGroup) then 1178 orbs:=ProbablyStabilizerOrbitNumbers(gamma.autGroup,V[1],gamma.order); 1179 else 1180 orbs:=ProbablyStabilizerOrbitNumbers(gamma.group,V[1],gamma.order); 1181 fi; 1182 else 1183 orbs:=rec(representatives:=[1..gamma.order], 1184 orbitNumbers:=[1..gamma.order]); 1185 fi; 1186fi; 1187orbnum:=orbs.orbitNumbers; 1188reps:=orbs.representatives; 1189laynum:=[]; 1190for i in [1..Length(reps)] do 1191 laynum[i]:=0; 1192od; 1193localGirth:=-1; 1194distance:=-1; 1195localParameters:=[]; 1196next:=[]; 1197for i in V do 1198 AddSet(next,orbnum[i]); 1199od; 1200sum:=Length(next); 1201for i in next do 1202 laynum[i]:=1; 1203od; 1204layers:=[]; 1205layers[1]:=next; 1206i:=1; 1207if Length(Intersection(V,stopvertices)) > 0 then 1208 stoplayer:=1; 1209 distance:=0; 1210fi; 1211while stoplayer<>i and Length(next)>0 do 1212 next:=[]; 1213 for x in layers[i] do 1214 nprev:=0; 1215 nhere:=0; 1216 nnext:=0; 1217 for y in Adjacency(gamma,reps[x]) do 1218 lnum:=laynum[orbnum[y]]; 1219 if i>1 and lnum=i-1 then 1220 nprev:=nprev+1; 1221 elif lnum=i then 1222 nhere:=nhere+1; 1223 elif lnum=i+1 then 1224 nnext:=nnext+1; 1225 elif lnum=0 then 1226 AddSet(next,orbnum[y]); 1227 nnext:=nnext+1; 1228 laynum[orbnum[y]]:=i+1; 1229 fi; 1230 od; 1231 if (localGirth=-1 or localGirth=2*i-1) and nprev>1 then 1232 localGirth:=2*(i-1); 1233 fi; 1234 if localGirth=-1 and nhere>0 then 1235 localGirth:=2*i-1; 1236 fi; 1237 if not IsBound(localParameters[i]) then 1238 localParameters[i]:=[nprev,nhere,nnext]; 1239 else 1240 if nprev<>localParameters[i][1] then 1241 localParameters[i][1]:=-1; 1242 fi; 1243 if nhere<>localParameters[i][2] then 1244 localParameters[i][2]:=-1; 1245 fi; 1246 if nnext<>localParameters[i][3] then 1247 localParameters[i][3]:=-1; 1248 fi; 1249 fi; 1250 od; 1251 if Length(next)>0 then 1252 i:=i+1; 1253 layers[i]:=next; 1254 for j in stopvertices do 1255 if laynum[orbnum[j]]=i then 1256 stoplayer:=i; 1257 distance:=i-1; 1258 break; 1259 fi; 1260 od; 1261 sum:=sum+Length(next); 1262 fi; 1263od; 1264if sum=Length(reps) then 1265 localDiameter:=Length(layers)-1; 1266else 1267 localDiameter:=-1; 1268fi; 1269# now change orbnum to give the layer numbers instead of orbit numbers. 1270layerNumbers:=orbnum; 1271for i in [1..gamma.order] do 1272 layerNumbers[i]:=laynum[orbnum[i]]; 1273od; 1274loc:=rec(layerNumbers:=layerNumbers,localDiameter:=localDiameter, 1275 localGirth:=localGirth,localParameters:=localParameters); 1276if Length(stopvertices) > 0 then 1277 loc.distance:=distance; 1278fi; 1279return loc; 1280end); 1281 1282BindGlobal("LocalInfoMat",function(A,rows) 1283# 1284# Calculates local info on a graph using a collapsed adjacency matrix 1285# A for that graph. 1286# This local info is from the point of view of the set of vertices 1287# represented by the set rows of row indices of A. 1288# The elements of layers[i] will be the row indices representing 1289# the vertices of the i-th layer. 1290# No distance field will be calculated. 1291# 1292# *** If A is not the collapsed adjacency matrix for a simple graph 1293# then localGirth and localParameters may not be what you think. 1294# If rows does not represent a single vertex then localGirth has 1295# no real meaning. 1296# 1297local layers,localDiameter,localGirth,localParameters,i,j,x,y,next, 1298 nprev,nhere,nnext,sum,laynum,lnum,n; 1299if IsInt(rows) then 1300 rows:=[rows]; 1301fi; 1302if not IsMatrix(A) or not IsList(rows) then 1303 Error("usage: LocalInfoMat( <Matrix>, <Int> or <List> )"); 1304fi; 1305if not IsSSortedList(rows) then 1306 rows:=SSortedList(rows); 1307fi; 1308n:=Length(A); 1309if rows=[] or not IsSubset([1..n],rows) then 1310 Error("<rows> must be non-empty set of row indices"); 1311fi; 1312laynum:=ListWithIdenticalEntries(n,0); 1313localGirth:=-1; 1314localParameters:=[]; 1315next:=ShallowCopy(rows); 1316for i in next do 1317 laynum[i]:=1; 1318od; 1319layers:=[]; 1320layers[1]:=next; 1321i:=1; 1322sum:=Length(rows); 1323while Length(next)>0 do 1324 next:=[]; 1325 for x in layers[i] do 1326 nprev:=0; 1327 nhere:=0; 1328 nnext:=0; 1329 for y in [1..n] do 1330 j:=A[x][y]; 1331 if j>0 then 1332 lnum:=laynum[y]; 1333 if i>1 and lnum=i-1 then 1334 nprev:=nprev+j; 1335 elif lnum=i then 1336 nhere:=nhere+j; 1337 elif lnum=i+1 then 1338 nnext:=nnext+j; 1339 elif lnum=0 then 1340 AddSet(next,y); 1341 nnext:=nnext+j; 1342 laynum[y]:=i+1; 1343 fi; 1344 fi; 1345 od; 1346 if (localGirth=-1 or localGirth=2*i-1) and nprev>1 then 1347 localGirth:=2*(i-1); 1348 fi; 1349 if localGirth=-1 and nhere>0 then 1350 localGirth:=2*i-1; 1351 fi; 1352 if not IsBound(localParameters[i]) then 1353 localParameters[i]:=[nprev,nhere,nnext]; 1354 else 1355 if nprev<>localParameters[i][1] then 1356 localParameters[i][1]:=-1; 1357 fi; 1358 if nhere<>localParameters[i][2] then 1359 localParameters[i][2]:=-1; 1360 fi; 1361 if nnext<>localParameters[i][3] then 1362 localParameters[i][3]:=-1; 1363 fi; 1364 fi; 1365 od; 1366 if Length(next)>0 then 1367 i:=i+1; 1368 layers[i]:=next; 1369 sum:=sum+Length(next); 1370 fi; 1371od; 1372if sum=n then 1373 localDiameter:=Length(layers)-1; 1374else 1375 localDiameter:=-1; 1376fi; 1377return rec(layerNumbers:=laynum,localDiameter:=localDiameter, 1378 localGirth:=localGirth,localParameters:=localParameters); 1379end); 1380 1381BindGlobal("InducedSubgraph",function(arg) 1382# 1383# Returns the subgraph of gamma=arg[1] induced on the list V=arg[2] of 1384# distinct vertices of gamma. 1385# If arg[3] is unbound, then the trivial group is the group associated 1386# with the returned induced subgraph. 1387# If arg[3] is bound, this function assumes that G=arg[3] fixes 1388# V setwise, and is a group of automorphisms of the induced subgraph 1389# when restriced to V. In this case, the image of G acting on V is 1390# the group associated with the returned induced subgraph. 1391# 1392# The i-th vertex of the induced subgraph corresponds to vertex V[i] of 1393# gamma, with the i-th vertex-name of the induced subgraph being the 1394# vertex-name in gamma of V[i]. 1395# 1396local gamma,V,G,Ggens,gens,indu,i,j,W,VV,X; 1397gamma:=arg[1]; 1398V:=arg[2]; 1399if not IsGraph(gamma) or not IsList(V) 1400 or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then 1401 Error("usage: InducedSubgraph( <Graph>, <List>, [, <PermGroup> ] )"); 1402fi; 1403VV:=SSortedList(V); 1404if Length(V)<>Length(VV) then 1405 Error("<V> must not contain repeated elements"); 1406fi; 1407if not IsSubset([1..gamma.order],VV) then 1408 Error("<V> must be a list of vertices of <gamma>"); 1409fi; 1410if IsBound(arg[3]) then 1411 G:=arg[3]; 1412else 1413 G:=Group([],()); 1414fi; 1415W:=[]; 1416for i in [1..Length(V)] do 1417 W[V[i]]:=i; 1418od; 1419Ggens:=GeneratorsOfGroup(G); 1420gens:=[]; 1421for i in [1..Length(Ggens)] do 1422 gens[i]:=[]; 1423 for j in V do 1424 gens[i][W[j]]:=W[j^Ggens[i]]; 1425 od; 1426 gens[i]:=PermList(gens[i]); 1427od; 1428indu:=NullGraph(Group(gens,()),Length(V)); 1429if IsBound(gamma.isSimple) and gamma.isSimple then 1430 indu.isSimple:=true; 1431else 1432 Unbind(indu.isSimple); 1433fi; 1434for i in [1..Length(indu.representatives)] do 1435 X:=W{Intersection(VV,Adjacency(gamma,V[indu.representatives[i]]))}; 1436 Sort(X); 1437 indu.adjacencies[i]:=X; 1438od; 1439if not IsBound(gamma.names) then 1440 indu.names:=Immutable(V); 1441else 1442 indu.names:=Immutable(gamma.names{V}); 1443fi; 1444return indu; 1445end); 1446 1447BindGlobal("Distance",function(arg) 1448# 1449# Let gamma=arg[1], X=arg[2], Y=arg[3]. 1450# Returns the distance d(X,Y) in the graph gamma, where X,Y 1451# are singleton vertices or lists of vertices. 1452# (Returns -1 if no (directed) path joins X to Y in gamma.) 1453# If arg[4] is bound, then it is assumed to be a subgroup 1454# of Aut(gamma) stabilizing X setwise. 1455# 1456local gamma,X,Y; 1457gamma:=arg[1]; 1458X:=arg[2]; 1459if IsInt(X) then 1460 X:=[X]; 1461fi; 1462Y:=arg[3]; 1463if IsInt(Y) then 1464 Y:=[Y]; 1465fi; 1466if not (IsGraph(gamma) and IsList(X) and IsList(Y)) then 1467 Error("usage: Distance( <Graph>, <Int> or <List>, ", 1468 "<Int> or <List> [, <PermGroup> ] )"); 1469fi; 1470if IsBound(arg[4]) then 1471 return LocalInfo(gamma,X,0,Y,arg[4]).distance; 1472else 1473 return LocalInfo(gamma,X,0,Y).distance; 1474fi; 1475end); 1476 1477DeclareOperation("Diameter",[IsRecord]); 1478# to avoid the clash with `Diameter' defined in gap4r5 1479InstallMethod(Diameter,"for GRAPE graph",[IsRecord],0, 1480function(gamma) 1481# 1482# Returns the diameter of gamma. 1483# A diameter of -1 means that gamma is not (strongly) connected. 1484# 1485local r,d,loc,reps; 1486if not IsGraph(gamma) then 1487 TryNextMethod(); 1488fi; 1489if gamma.order=0 then 1490 Error("<gamma> has no vertices"); 1491fi; 1492d:=-1; 1493if IsBound(gamma.autGroup) then 1494 reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives; 1495else 1496 reps:=gamma.representatives; 1497fi; 1498for r in reps do 1499 loc:=LocalInfo(gamma,r); 1500 if loc.localDiameter=-1 then 1501 return -1; 1502 fi; 1503 if loc.localDiameter > d then 1504 d:=loc.localDiameter; 1505 fi; 1506od; 1507return d; 1508end); 1509 1510DeclareOperation("Girth",[IsRecord]); 1511InstallMethod(Girth,"for GRAPE graph",[IsRecord],0, 1512function(gamma) 1513# 1514# Returns the girth of gamma, which must be a simple graph. 1515# A girth of -1 means that gamma is a forest. 1516# 1517local r,g,locgirth,stoplayer,adj,reps; 1518if not IsGraph(gamma) then 1519 TryNextMethod(); 1520fi; 1521if gamma.order=0 then 1522 return -1; 1523fi; 1524if not IsSimpleGraph(gamma) then 1525 Error("<gamma> not a simple graph"); 1526fi; 1527adj:=gamma.adjacencies[1]; 1528if adj<>[] and Intersection(adj,Adjacency(gamma,adj[1]))<>[] then 1529 return 3; 1530fi; 1531g:=-1; 1532stoplayer:=0; 1533if IsBound(gamma.autGroup) then 1534 reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives; 1535else 1536 reps:=gamma.representatives; 1537fi; 1538for r in reps do 1539 locgirth:=LocalInfo(gamma,r,stoplayer).localGirth; 1540 if locgirth=3 then 1541 return 3; 1542 fi; 1543 if locgirth<>-1 then 1544 if g=-1 or locgirth<g then 1545 g:=locgirth; 1546 stoplayer:=Int((g+1)/2)+1; # now no need for LocalInfo to create a 1547 # layer beyond the stoplayer-th one. 1548 fi; 1549 fi; 1550od; 1551return g; 1552end); 1553 1554BindGlobal("IsRegularGraph",function(gamma) 1555# 1556# Returns true iff the graph gamma is (out)regular. 1557# 1558local deg,i; 1559if not IsGraph(gamma) then 1560 Error("usage: IsRegularGraph( <Graph> )"); 1561fi; 1562if gamma.order=0 then 1563 return true; 1564fi; 1565deg:=Length(gamma.adjacencies[1]); 1566for i in [2..Length(gamma.adjacencies)] do 1567 if deg <> Length(gamma.adjacencies[i]) then 1568 return false; 1569 fi; 1570od; 1571return true; 1572end); 1573 1574BindGlobal("IsNullGraph",function(gamma) 1575# 1576# Returns true iff the graph gamma has no edges. 1577# 1578local i; 1579if not IsGraph(gamma) then 1580 Error("usage: IsNullGraph( <Graph> )"); 1581fi; 1582for i in [1..Length(gamma.adjacencies)] do 1583 if Length(gamma.adjacencies[i])<>0 then 1584 return false; 1585 fi; 1586od; 1587return true; 1588end); 1589 1590BindGlobal("IsCompleteGraph",function(arg) 1591# 1592# Returns true iff the graph gamma=arg[1] is a complete graph. 1593# The optional boolean parameter arg[2] 1594# is true iff all loops must exist for gamma to be considered 1595# a complete graph (default: false); otherwise loops are ignored 1596# (except to possibly set gamma.isSimple). 1597# 1598local deg,i,notnecsimple,gamma,mustloops; 1599gamma:=arg[1]; 1600if IsBound(arg[2]) then 1601 mustloops := arg[2]; 1602else 1603 mustloops := false; 1604fi; 1605if not IsGraph(gamma) or not IsBool(mustloops) then 1606 Error("usage: IsCompleteGraph( <Graph> [, <Bool> ] )"); 1607fi; 1608notnecsimple := not IsBound(gamma.isSimple) or not gamma.isSimple; 1609for i in [1..Length(gamma.adjacencies)] do 1610 deg := Length(gamma.adjacencies[i]); 1611 if deg < gamma.order-1 then 1612 return false; 1613 fi; 1614 if deg=gamma.order-1 then 1615 if mustloops then 1616 return false; 1617 fi; 1618 if notnecsimple and 1619 (gamma.representatives[i] in gamma.adjacencies[i]) then 1620 gamma.isSimple := false; 1621 return false; 1622 fi; 1623 fi; 1624od; 1625return true; 1626end); 1627 1628BindGlobal("IsLoopy",function(gamma) 1629# 1630# Returns true iff graph gamma has a loop. 1631# 1632local i; 1633if not IsGraph(gamma) then 1634 Error("usage: IsLoopy( <Graph> )"); 1635fi; 1636for i in [1..Length(gamma.adjacencies)] do 1637 if gamma.representatives[i] in gamma.adjacencies[i] then 1638 gamma.isSimple := false; 1639 return true; 1640 fi; 1641od; 1642return false; 1643end); 1644 1645DeclareOperation("IsConnectedGraph",[IsRecord]); 1646InstallMethod(IsConnectedGraph,"for GRAPE graph",[IsRecord],0, 1647function(gamma) 1648# 1649# Returns true iff gamma is (strongly) connected. 1650# 1651if not IsGraph(gamma) then 1652 TryNextMethod(); 1653fi; 1654if gamma.order=0 then 1655 return true; 1656fi; 1657if IsSimpleGraph(gamma) then 1658 return LocalInfo(gamma,1).localDiameter > -1; 1659else 1660 return Diameter(gamma) > -1; 1661fi; 1662end); 1663 1664 1665DeclareOperation("ConnectedComponent",[IsRecord,IsPosInt]); 1666InstallMethod(ConnectedComponent,"for GRAPE graph",[IsRecord,IsPosInt],0, 1667function(gamma,v) 1668# 1669# Returns the set of all vertices in gamma which can be reached by 1670# a path starting at vertex v. The graph gamma must be simple. 1671# 1672local comp,laynum; 1673if not IsGraph(gamma) then 1674 TryNextMethod(); 1675fi; 1676if not IsSimpleGraph(gamma) then 1677 Error("<gamma> not a simple graph"); 1678fi; 1679if not IsVertex(gamma,v) then 1680 Error("<v> is not a vertex of <gamma>"); 1681fi; 1682laynum:=LocalInfo(gamma,v).layerNumbers; 1683comp:=Filtered([1..gamma.order],j->laynum[j]>0); 1684IsSSortedList(comp); 1685return comp; 1686end); 1687 1688DeclareOperation("ConnectedComponents",[IsRecord]); 1689InstallMethod(ConnectedComponents,"for GRAPE graph",[IsRecord],0, 1690function(gamma) 1691# 1692# Returns the set of the vertex-sets of the connected components 1693# of gamma, which must be a simple graph. 1694# 1695local comp,used,i,j,x,cmp,laynum; 1696if not IsGraph(gamma) then 1697 TryNextMethod(); 1698fi; 1699if not IsSimpleGraph(gamma) then 1700 Error("<gamma> not a simple graph"); 1701fi; 1702comp:=[]; 1703used:=BlistList([1..gamma.order],[]); 1704for i in [1..gamma.order] do 1705 # Loop invariant: used[j]=true for all j<i. 1706 if not used[i] then # new component 1707 laynum:=LocalInfo(gamma,i).layerNumbers; 1708 cmp:=Filtered([i..gamma.order],j->laynum[j]>0); 1709 IsSSortedList(cmp); 1710 for x in Orbit(gamma.group,cmp,OnSets) do 1711 Add(comp,x); 1712 for j in x do 1713 used[j]:=true; 1714 od; 1715 od; 1716 fi; 1717od; 1718return SSortedList(comp); 1719end); 1720 1721BindGlobal("ComponentLocalInfos",function(gamma) 1722# 1723# Returns a sequence of localinfos for the connected components of 1724# gamma (w.r.t. some vertex in each component). 1725# The graph gamma must be simple. 1726# 1727local comp,used,i,j,k,laynum; 1728if not IsGraph(gamma) then 1729 Error("usage: ComponentLocalInfos( <Graph> )"); 1730fi; 1731if not IsSimpleGraph(gamma) then 1732 Error("<gamma> not a simple graph"); 1733fi; 1734comp:=[]; 1735used:=BlistList([1..gamma.order],[]); 1736k:=0; 1737for i in [1..gamma.order] do 1738 if not used[i] then # new component 1739 k:=k+1; 1740 comp[k]:=LocalInfo(gamma,i); 1741 laynum:=comp[k].layerNumbers; 1742 for j in [1..gamma.order] do 1743 if laynum[j] > 0 then 1744 used[j]:=true; 1745 fi; 1746 od; 1747 fi; 1748od; 1749return comp; 1750end); 1751 1752BindGlobal("Bicomponents",function(gamma) 1753# 1754# If gamma is bipartite, returns a length 2 list of 1755# bicomponents, or parts, of gamma, else returns the empty list. 1756# *** This function is for simple gamma only. 1757# 1758# Note: if gamma.order=0 this function returns [[],[]], and if 1759# gamma.order=1 this function returns [[],[1]] (unlike GRAPE 2.2 1760# which returned [], which was inconsistent with considering 1761# a zero vertex graph to be bipartite). 1762# 1763local bicomps,i,lnum,loc,locs; 1764if not IsGraph(gamma) then 1765 Error("usage: Bicomponents( <Graph> )"); 1766fi; 1767if not IsSimpleGraph(gamma) then 1768 Error("<gamma> not a simple graph"); 1769fi; 1770bicomps:=[[],[]]; 1771if gamma.order=0 then 1772 return bicomps; 1773fi; 1774if IsNullGraph(gamma) then 1775 return [[1..gamma.order-1],[gamma.order]]; 1776fi; 1777locs:=ComponentLocalInfos(gamma); 1778for loc in locs do 1779 for i in [2..Length(loc.localParameters)] do 1780 if loc.localParameters[i][2]<>0 then 1781 # gamma not bipartite. 1782 return []; 1783 fi; 1784 od; 1785 for i in [1..Length(loc.layerNumbers)] do 1786 lnum:=loc.layerNumbers[i]; 1787 if lnum>0 then 1788 if lnum mod 2 = 1 then 1789 AddSet(bicomps[1],i); 1790 else 1791 AddSet(bicomps[2],i); 1792 fi; 1793 fi; 1794 od; 1795od; 1796return bicomps; 1797end); 1798 1799DeclareOperation("IsBipartite",[IsRecord]); 1800InstallMethod(IsBipartite,"for GRAPE graph",[IsRecord],0, 1801function(gamma) 1802# 1803# Returns true iff gamma is bipartite. 1804# *** This function is only for simple gamma. 1805# 1806# Note: Now the one vertex graph is considered to be bipartite 1807# (as well as the zero vertex graph). This is a change from the inconsistent 1808# GRAPE 2.2 view that a zero vertex graph is bipartite, but not a one 1809# vertex graph. 1810# 1811if not IsGraph(gamma) then 1812 TryNextMethod(); 1813fi; 1814if not IsSimpleGraph(gamma) then 1815 Error("<gamma> not a simple graph"); 1816fi; 1817return Length(Bicomponents(gamma))=2; 1818end); 1819 1820BindGlobal("Layers",function(arg) 1821# 1822# Returns the list of vertex layers of gamma=arg[1], 1823# starting from V=arg[2], which may be a vertex list or singleton vertex. 1824# Layers[i] is the set of vertices at distance i-1 from V. 1825# If arg[3] is bound then it is assumed to be a subgroup 1826# of Aut(gamma) stabilizing V setwise. 1827# 1828local gamma,V; 1829gamma:=arg[1]; 1830V:=arg[2]; 1831if IsInt(V) then 1832 V:=[V]; 1833fi; 1834if not IsGraph(gamma) or not IsList(V) 1835 or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then 1836 Error("usage: Layers( <Graph>, <Int> or <List>, [, <PermGroup>] )"); 1837fi; 1838if IsBound(arg[3]) then 1839 return GRAPE_NumbersToSets(LocalInfo(gamma,V,0,[],arg[3]).layerNumbers); 1840else 1841 return GRAPE_NumbersToSets(LocalInfo(gamma,V).layerNumbers); 1842fi; 1843end); 1844 1845BindGlobal("LocalParameters",function(arg) 1846# 1847# Returns the local parameters of simple, connected gamma=arg[1], 1848# w.r.t to vertex list (or singleton vertex) V=arg[2]. 1849# The nonexistence of a local parameter is denoted by -1. 1850# If arg[3] is bound then it is assumed to be a subgroup 1851# of Aut(gamma) stabilizing V setwise. 1852# 1853local gamma,V,loc; 1854gamma:=arg[1]; 1855V:=arg[2]; 1856if IsInt(V) then 1857 V:=[V]; 1858fi; 1859if not IsGraph(gamma) or not IsList(V) 1860 or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then 1861 Error("usage: LocalParameters( <Graph>, <Int> or <List>, [, <PermGroup>] )"); 1862fi; 1863if not IsSimpleGraph(gamma) then 1864 Error("<gamma> not a simple graph"); 1865fi; 1866if Length(V)>1 and not IsConnectedGraph(gamma) then 1867 Error("<gamma> not a connected graph"); 1868fi; 1869if IsBound(arg[3]) then 1870 loc:=LocalInfo(gamma,V,0,[],arg[3]); 1871else 1872 loc:=LocalInfo(gamma,V); 1873fi; 1874if loc.localDiameter=-1 then 1875 Error("<gamma> not a connected graph"); 1876fi; 1877return loc.localParameters; 1878end); 1879 1880BindGlobal("GlobalParameters",function(gamma) 1881# 1882# Determines the global parameters of connected, simple graph gamma. 1883# The nonexistence of a global parameter is denoted by -1. 1884# 1885local i,j,k,reps,pars,lp,loc; 1886if not IsGraph(gamma) then 1887 Error("usage: GlobalParameters( <Graph> )"); 1888fi; 1889if gamma.order=0 then 1890 return []; 1891fi; 1892if not IsSimpleGraph(gamma) then 1893 Error("<gamma> not a simple graph"); 1894fi; 1895if IsBound(gamma.autGroup) then 1896 reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives; 1897else 1898 reps:=gamma.representatives; 1899fi; 1900loc:=LocalInfo(gamma,reps[1]); 1901if loc.localDiameter=-1 then 1902 Error("<gamma> not a connected graph"); 1903fi; 1904pars:=loc.localParameters; 1905for i in [2..Length(reps)] do 1906 lp:=LocalInfo(gamma,reps[i]).localParameters; 1907 for j in [1..Maximum(Length(lp),Length(pars))] do 1908 if not IsBound(lp[j]) or not IsBound(pars[j]) then 1909 pars[j]:=[-1,-1,-1]; 1910 else 1911 for k in [1..3] do 1912 if pars[j][k]<>lp[j][k] then 1913 pars[j][k]:=-1; 1914 fi; 1915 od; 1916 fi; 1917 od; 1918od; 1919return pars; 1920end); 1921 1922DeclareOperation("IsDistanceRegular",[IsRecord]); 1923InstallMethod(IsDistanceRegular,"for GRAPE graph",[IsRecord],0, 1924function(gamma) 1925# 1926# Returns true iff gamma is distance-regular 1927# (a graph must be simple to be distance-regular). 1928# 1929local i,reps,pars,lp,loc,d; 1930if not IsGraph(gamma) then 1931 TryNextMethod(); 1932fi; 1933if gamma.order=0 then 1934 return true; 1935fi; 1936if not IsSimpleGraph(gamma) then 1937 return false; 1938fi; 1939if IsBound(gamma.autGroup) then 1940 reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives; 1941else 1942 reps:=gamma.representatives; 1943fi; 1944loc:=LocalInfo(gamma,reps[1]); 1945pars:=loc.localParameters; 1946d:=loc.localDiameter; 1947if d=-1 then # gamma not connected 1948 return false; 1949fi; 1950if -1 in Flat(pars) then 1951 return false; 1952fi; 1953for i in [2..Length(reps)] do 1954 loc:=LocalInfo(gamma,reps[i]); 1955 if loc.localDiameter<>d then 1956 return false; 1957 fi; 1958 if pars <> loc.localParameters then 1959 return false; 1960 fi; 1961od; 1962return true; 1963end); 1964 1965BindGlobal("DistanceSet",function(arg) 1966# 1967# Let gamma=arg[1], distances=arg[2], V=arg[3]. 1968# Returns the set of vertices w of gamma, such that d(V,w) is in 1969# distances (a list or singleton distance). 1970# If arg[4] is bound, then it is assumed to be a subgroup 1971# of Aut(gamma) stabilizing V setwise. 1972# 1973local gamma,distances,V,maxlayer,distset,laynum,x,i; 1974gamma:=arg[1]; 1975distances:=arg[2]; 1976V:=arg[3]; 1977if IsInt(distances) then # assume distances consists of a single distance. 1978 distances:=[distances]; 1979fi; 1980if not (IsGraph(gamma) and IsList(distances) and (IsList(V) or IsInt(V))) then 1981 Error("usage: DistanceSet( <Graph>, <Int> or <List>, ", 1982 "<Int> or <List> [, <PermGroup> ] )"); 1983fi; 1984if not IsSSortedList(distances) then 1985 distances:=SSortedList(distances); 1986fi; 1987distset:=[]; 1988if Length(distances)=0 then 1989 return distset; 1990fi; 1991maxlayer:=Maximum(distances)+1; 1992if IsBound(arg[4]) then 1993 laynum:=LocalInfo(gamma,V,maxlayer,[],arg[4]).layerNumbers; 1994else 1995 laynum:=LocalInfo(gamma,V,maxlayer).layerNumbers; 1996fi; 1997for i in [1..gamma.order] do 1998 if laynum[i]-1 in distances then 1999 Add(distset,i); 2000 fi; 2001od; 2002IsSSortedList(distset); 2003return distset; 2004end); 2005 2006BindGlobal("DistanceSetInduced",function(arg) 2007# 2008# Let gamma=arg[1], distances=arg[2], V=arg[3]. 2009# Returns the graph induced on the set of vertices w of gamma, 2010# such that d(V,w) is in distances (a list or singleton distance). 2011# If arg[4] is bound, then it is assumed to be a subgroup 2012# of Aut(gamma) stabilizing V setwise. 2013# 2014local gamma,distances,V,distset,H; 2015gamma:=arg[1]; 2016distances:=arg[2]; 2017V:=arg[3]; 2018if IsInt(distances) then # assume distances consists of a single distance. 2019 distances:=[distances]; 2020fi; 2021if IsInt(V) then 2022 V:=[V]; 2023fi; 2024if IsBound(arg[4]) then 2025 H:=arg[4]; 2026elif Length(V)=1 then 2027 H:=ProbablyStabilizer(gamma.group,V[1]); 2028else 2029 H:=Group([],()); 2030fi; 2031if not (IsGraph(gamma) and IsList(distances) and IsList(V) and IsPermGroup(H)) 2032 then 2033 Error("usage: DistanceSetInduced( <Graph>, ", 2034 "<Int> or <List>, <Int> or <List> [, <PermGroup> ] )"); 2035fi; 2036distset:=DistanceSet(gamma,distances,V,H); 2037return InducedSubgraph(gamma,distset,H); 2038end); 2039 2040BindGlobal("DistanceGraph",function(gamma,distances) 2041# 2042# Returns graph delta with the same vertex-set, names, and group as 2043# gamma, and [x,y] is an edge of delta iff d(x,y) (in gamma) 2044# is in distances. 2045# 2046local r,delta,d,i; 2047if IsInt(distances) then 2048 distances:=[distances]; 2049fi; 2050if not IsGraph(gamma) or not IsList(distances) then 2051 Error("usage: DistanceGraph( <Graph>, <Int> or <List> )"); 2052fi; 2053delta:=rec(isGraph:=true,order:=gamma.order,group:=gamma.group, 2054 schreierVector:=Immutable(gamma.schreierVector),adjacencies:=[], 2055 representatives:=Immutable(gamma.representatives)); 2056if IsBound(gamma.names) then 2057 delta.names:=Immutable(gamma.names); 2058fi; 2059for i in [1..Length(delta.representatives)] do 2060 delta.adjacencies[i]:=DistanceSet(gamma,distances,delta.representatives[i]); 2061od; 2062if not (0 in distances) and IsBound(gamma.isSimple) and gamma.isSimple then 2063 delta.isSimple:=true; 2064fi; 2065return delta; 2066end); 2067 2068BindGlobal("ComplementGraph",function(arg) 2069# 2070# Returns the complement of the graph gamma=arg[1]. 2071# arg[2] is true iff loops/nonloops are to be complemented (default:false). 2072# 2073local gamma,comploops,i,delta,notnecsimple; 2074gamma:=arg[1]; 2075if IsBound(arg[2]) then 2076 comploops:=arg[2]; 2077else 2078 comploops:=false; 2079fi; 2080if not IsGraph(gamma) or not IsBool(comploops) then 2081 Error("usage: ComplementGraph( <Graph> [, <Bool> ] )"); 2082fi; 2083notnecsimple:=not IsBound(gamma.isSimple) or not gamma.isSimple; 2084delta:=rec(isGraph:=true,order:=gamma.order,group:=gamma.group, 2085 schreierVector:=Immutable(gamma.schreierVector),adjacencies:=[], 2086 representatives:=Immutable(gamma.representatives)); 2087if IsBound(gamma.names) then 2088 delta.names:=Immutable(gamma.names); 2089fi; 2090if IsBound(gamma.autGroup) then 2091 delta.autGroup:=gamma.autGroup; 2092fi; 2093if IsBound(gamma.isSimple) then 2094 if gamma.isSimple and not comploops then 2095 delta.isSimple:=true; 2096 fi; 2097fi; 2098for i in [1..Length(delta.representatives)] do 2099 delta.adjacencies[i]:=Difference([1..gamma.order],gamma.adjacencies[i]); 2100 if not comploops then 2101 RemoveSet(delta.adjacencies[i],delta.representatives[i]); 2102 if notnecsimple and (gamma.representatives[i] in gamma.adjacencies[i]) 2103 then 2104 AddSet(delta.adjacencies[i],delta.representatives[i]); 2105 fi; 2106 fi; 2107od; 2108return delta; 2109end); 2110 2111BindGlobal("PointGraph",function(arg) 2112# 2113# Assuming that gamma=arg[1] is simple, connected, and bipartite, 2114# this function returns the connected component containing 2115# v=arg[2] of the distance-2 graph of gamma=arg[1] 2116# (default: arg[2]=1, unless gamma has zero 2117# vertices, in which case a zero vertex graph is returned). 2118# Thus, if gamma is the incidence graph of a (connected) geometry, and 2119# v represents a point, then the point graph of the geometry is returned. 2120# 2121local gamma,delta,bicomps,comp,v,gens,hgens,i,g,j,outer; 2122gamma:=arg[1]; 2123if IsBound(arg[2]) then 2124 v:=arg[2]; 2125else 2126 v:=1; 2127fi; 2128if not IsGraph(gamma) or not IsInt(v) then 2129 Error("usage: PointGraph( <Graph> [, <Int> ])"); 2130fi; 2131if gamma.order=0 then 2132 return CopyGraph(gamma); 2133fi; 2134bicomps:=Bicomponents(gamma); 2135if Length(bicomps)=0 or not IsSimpleGraph(gamma) 2136 or not IsConnectedGraph(gamma) then 2137 Error("<gamma> not simple,connected,bipartite"); 2138fi; 2139if v in bicomps[1] then 2140 comp:=bicomps[1]; 2141else 2142 comp:=bicomps[2]; 2143fi; 2144delta:=DistanceGraph(gamma,2); 2145# construct Schreier generators for the subgroup of gamma.group 2146# fixing comp. 2147gens:=GeneratorsOfGroup(gamma.group); 2148hgens:=[]; 2149for i in [1..Length(gens)] do 2150 g:=gens[i]; 2151 if v^g in comp then 2152 AddSet(hgens,g); 2153 if IsBound(outer) then 2154 AddSet(hgens,outer*g/outer); 2155 fi; 2156 else # g is an "outer" element 2157 if IsBound(outer) then 2158 AddSet(hgens,g/outer); 2159 AddSet(hgens,outer*g); 2160 else 2161 outer:=g; 2162 for j in [1..i-1] do 2163 AddSet(hgens,outer*gens[j]/outer); 2164 od; 2165 g:=g^2; 2166 if g <> () then 2167 AddSet(hgens,g); 2168 fi; 2169 fi; 2170 fi; 2171od; 2172return InducedSubgraph(delta,comp,Group(hgens,())); 2173end); 2174 2175BindGlobal("EdgeGraph",function(gamma) 2176# 2177# Returns the edge graph, also called the line graph, of the 2178# (assumed) simple graph gamma. 2179# This edge graph delta has the unordered edges of gamma as 2180# vertices, and e is joined to f in delta precisely when 2181# e<>f, and e,f have a common vertex in gamma. 2182# 2183local delta,i,j,k,edgeset,adj,r,e,f; 2184if not IsGraph(gamma) then 2185 Error("usage: EdgeGraph( <Graph> )"); 2186fi; 2187if not IsSimpleGraph(gamma) then 2188 Error("<gamma> not a simple graph"); 2189fi; 2190edgeset:=UndirectedEdges(gamma); 2191delta:=NullGraph(Action(gamma.group,edgeset,OnSets),Length(edgeset)); 2192delta.names:=Immutable(List(edgeset,e->List(e,i->VertexName(gamma,i)))); 2193for i in [1..Length(delta.representatives)] do 2194 r:=delta.representatives[i]; 2195 e:=edgeset[r]; 2196 adj:=delta.adjacencies[i]; 2197 for k in [1,2] do 2198 for j in Adjacency(gamma,e[k]) do 2199 f:=SSortedList([e[k],j]); 2200 if e<>f then 2201 AddSet(adj,PositionSet(edgeset,f)); 2202 fi; 2203 od; 2204 od; 2205od; 2206return delta; 2207end); 2208 2209BindGlobal("QuotientGraph",function(gamma,R) 2210# 2211# Returns the quotient graph delta of gamma defined by the 2212# smallest gamma.group invariant equivalence relation S 2213# (on the vertices of gamma) containing the relation R (given 2214# as a list of ordered pairs of vertices of gamma). 2215# The vertices of this quotient delta are the equivalence 2216# classes of S, and [X,Y] is an edge of delta iff 2217# [x,y] is an edge of gamma for some x in X, y in Y. 2218# 2219local root,Q,F,V,W,i,j,r,q,x,y,names,gens,delta,g,h,m,pos; 2220 2221root := function(x) 2222# 2223# Returns the root of the tree containing x in the forest represented 2224# by F, and compresses the path travelled in this tree to find the root. 2225# F[x]=-x if x is a root, else F[x] is the parent of x. 2226# 2227local t,u; 2228t:=F[x]; 2229while t>0 do 2230 t:=F[t]; 2231od; 2232# compress path 2233u:=F[x]; 2234while u<>t do 2235 F[x]:=-t; 2236 x:=u; 2237 u:=F[x]; 2238od; 2239return -t; 2240end; 2241 2242if not IsGraph(gamma) or not IsList(R) then 2243 Error("usage: QuotientGraph( <Graph>, <List> )"); 2244fi; 2245if gamma.order<=1 or Length(R)=0 then 2246 delta:=CopyGraph(gamma); 2247 delta.names:=Immutable(List([1..gamma.order],i->[VertexName(gamma,i)])); 2248 return delta; 2249fi; 2250if IsInt(R[1]) then # assume R consists of a single pair. 2251 R:=[R]; 2252fi; 2253F:=[]; 2254for i in [1..gamma.order] do 2255 F[i]:=-i; 2256od; 2257Q:=[]; 2258for r in R do 2259 x:=root(r[1]); 2260 y:=root(r[2]); 2261 if x<>y then 2262 if x>y then 2263 Add(Q,x); 2264 F[x]:=y; 2265 else 2266 Add(Q,y); 2267 F[y]:=x; 2268 fi; 2269 fi; 2270od; 2271for q in Q do 2272 for g in GeneratorsOfGroup(gamma.group) do 2273 x:=root(F[q]^g); 2274 y:=root(q^g); 2275 if x<>y then 2276 if x>y then 2277 Add(Q,x); 2278 F[x]:=y; 2279 else 2280 Add(Q,y); 2281 F[y]:=x; 2282 fi; 2283 fi; 2284 od; 2285od; 2286for i in Reversed([1..gamma.order]) do 2287 if F[i] < 0 then 2288 F[i]:=-F[i]; 2289 else 2290 F[i]:=root(F[i]); 2291 fi; 2292od; 2293V:=SSortedList(F); 2294W:=[]; 2295names:=[]; 2296m:=Length(V); 2297for i in [1..m] do 2298 W[V[i]]:=i; 2299 names[i]:=[]; 2300od; 2301for i in [1..gamma.order] do 2302 Add(names[W[F[i]]],VertexName(gamma,i)); 2303od; 2304gens:=[]; 2305for g in GeneratorsOfGroup(gamma.group) do 2306 h:=[]; 2307 for i in [1..m] do 2308 h[i]:=W[F[V[i]^g]]; 2309 od; 2310 Add(gens,PermList(h)); 2311od; 2312delta:=NullGraph(Group(gens,()),m); 2313delta.names:=Immutable(names); 2314IsSSortedList(delta.representatives); 2315for i in [1..gamma.order] do 2316 pos:=PositionSet(delta.representatives,W[F[i]]); 2317 if pos<>fail then 2318 for j in Adjacency(gamma,i) do 2319 AddSet(delta.adjacencies[pos],W[F[j]]); 2320 od; 2321 fi; 2322od; 2323if IsLoopy(delta) then 2324 delta.isSimple:=false; 2325elif IsBound(gamma.isSimple) and gamma.isSimple then 2326 delta.isSimple:=true; 2327else 2328 Unbind(delta.isSimple); 2329fi; 2330return delta; 2331end); 2332 2333BindGlobal("BipartiteDouble",function(gamma) 2334# 2335# Returns the bipartite double of gamma, as defined in BCN. 2336# 2337local gens,g,delta,n,i,adj; 2338if not IsGraph(gamma) then 2339 Error("usage: BipartiteDouble( <Graph> )"); 2340fi; 2341if gamma.order=0 then 2342 return CopyGraph(gamma); 2343fi; 2344n:=gamma.order; 2345gens:=GRAPE_IntransitiveGroupGenerators 2346 (GeneratorsOfGroup(gamma.group),GeneratorsOfGroup(gamma.group),n,n); 2347g:=[]; 2348for i in [1..n] do 2349 g[i]:=i+n; 2350 g[i+n]:=i; 2351od; 2352Add(gens,PermList(g)); 2353delta:=NullGraph(Group(gens,()),2*n); 2354for i in [1..Length(delta.adjacencies)] do 2355 adj:=Adjacency(gamma,delta.representatives[i]); 2356 if Length(adj)=0 then 2357 delta.adjacencies[i]:=[]; 2358 else 2359 delta.adjacencies[i]:=adj+n; 2360 fi; 2361od; 2362if not IsBound(gamma.isSimple) or not gamma.isSimple then 2363 Unbind(delta.isSimple); 2364fi; 2365delta.names:=[]; 2366for i in [1..n] do 2367 delta.names[i]:=[VertexName(gamma,i),"+"]; 2368od; 2369for i in [n+1..2*n] do 2370 delta.names[i]:=[VertexName(gamma,i-n),"-"]; 2371od; 2372delta.names:=Immutable(delta.names); 2373return delta; 2374end); 2375 2376BindGlobal("UnderlyingGraph",function(gamma) 2377# 2378# Returns the underlying graph delta of gamma. 2379# This graph has the same vertex-set as gamma, and 2380# has an edge [x,y] precisely when gamma has an 2381# edge [x,y] or [y,x]. 2382# This function also sets the isSimple fields of 2383# gamma (via IsSimpleGraph) and delta. 2384# 2385local delta,adj,i,x,orb,H; 2386if not IsGraph(gamma) then 2387 Error("usage: UnderlyingGraph( <Graph> )"); 2388fi; 2389delta:=CopyGraph(gamma); 2390if IsSimpleGraph(gamma) then 2391 delta.isSimple:=true; 2392 return delta; 2393fi; 2394for i in [1..Length(delta.adjacencies)] do 2395 adj:=delta.adjacencies[i]; 2396 x:=delta.representatives[i]; 2397 H:=ProbablyStabilizer(delta.group,x); 2398 for orb in OrbitsDomain(H,adj) do 2399 if not IsVertexPairEdge(delta,orb[1],x) then 2400 AddEdgeOrbit(delta,[orb[1],x]); 2401 fi; 2402 od; 2403od; 2404delta.isSimple := not IsLoopy(delta); 2405return delta; 2406end); 2407 2408BindGlobal("NewGroupGraph",function(G,gamma) 2409# 2410# Returns a copy of delta of gamma, except that delta.group=G. 2411# 2412local delta,i; 2413if not IsPermGroup(G) or not IsGraph(gamma) then 2414 Error("usage: NewGroupGraph( <PermGroup>, <Graph> )"); 2415fi; 2416delta:=NullGraph(G,gamma.order); 2417if IsBound(gamma.isSimple) then 2418 delta.isSimple:=gamma.isSimple; 2419else 2420 Unbind(delta.isSimple); 2421fi; 2422if IsBound(gamma.autGroup) then 2423 delta.autGroup:=gamma.autGroup; 2424fi; 2425# if IsBound(gamma.canonicalLabelling) then 2426# delta.canonicalLabelling:=gamma.canonicalLabelling; 2427# fi; 2428if IsBound(gamma.names) then 2429 delta.names:=Immutable(gamma.names); 2430fi; 2431if IsBound(gamma.maximumClique) then 2432 delta.maximumClique:=Immutable(gamma.maximumClique); 2433fi; 2434if IsBound(gamma.minimumVertexColouring) then 2435 delta.minimumVertexColouring:=Immutable(gamma.minimumVertexColouring); 2436fi; 2437for i in [1..Length(delta.representatives)] do 2438 delta.adjacencies[i]:=Adjacency(gamma,delta.representatives[i]); 2439od; 2440return delta; 2441end); 2442 2443BindGlobal("GeodesicsGraph",function(gamma,x,y) 2444# 2445# Returns the graph induced on the set of geodesics between 2446# vertices x and y, but not including x or y. 2447# *** This function is only for simple gamma. 2448# 2449local i,n,locx,geoset,H,laynumx,laynumy,w,g,h,rwx,rwy,gens,sch,orb,pt,im; 2450if not IsGraph(gamma) or not IsInt(x) or not IsInt(y) then 2451 Error("usage: GeodesicsGraph( <Graph>, <Int>, <Int> )"); 2452fi; 2453if not IsSimpleGraph(gamma) then 2454 Error("<gamma> not a simple graph"); 2455fi; 2456locx:=LocalInfo(gamma,x,0,y); 2457if locx.distance=-1 then 2458 Error("<x> not joined to <y>"); 2459fi; 2460laynumx:=locx.layerNumbers; 2461laynumy:=LocalInfo(gamma,y,0,x).layerNumbers; 2462geoset:=[]; 2463n:=locx.distance; 2464for i in [1..gamma.order] do 2465 if laynumx[i]>1 and laynumy[i]>1 and laynumx[i]+laynumy[i]=n+2 then 2466 Add(geoset,i); 2467 fi; 2468od; 2469H:=ProbablyStabilizer(gamma.group,y); 2470gens:=GeneratorsOfGroup(gamma.group); 2471rwx:=GRAPE_RepWord(gens,gamma.schreierVector,x); 2472rwy:=GRAPE_RepWord(gens,gamma.schreierVector,y); 2473g:=(); 2474if rwx.representative=rwy.representative then 2475 for w in Reversed(rwx.word) do 2476 g:=g/gens[w]; 2477 od; 2478 for w in rwy.word do 2479 g:=g*gens[w]; 2480 od; 2481 pt:=y^g; 2482 sch:=[]; 2483 orb:=[pt]; 2484 sch[pt]:=(); 2485 for i in orb do 2486 for h in GeneratorsOfGroup(H) do 2487 im:=i^h; 2488 if not IsBound(sch[im]) then 2489 sch[im]:=h; 2490 Add(orb,im); 2491 fi; 2492 od; 2493 od; 2494 if IsBound(sch[x]) then 2495 i:=x; 2496 h:=(); 2497 while i<>pt do 2498 h:=h/sch[i]; 2499 i:=x^h; 2500 od; 2501 g:=g*h^-1; 2502 fi; 2503fi; 2504H:=ProbablyStabilizer(H,x); 2505if x^g=y and y^g=x then 2506 H:=ShallowCopy(GeneratorsOfGroup(H)); 2507 Add(H,g); 2508 H:=Group(H,()); 2509fi; 2510return InducedSubgraph(gamma,geoset,H); 2511end); 2512 2513BindGlobal("IndependentSet",function(arg) 2514# 2515# Returns a (hopefully large) independent set (coclique) of gamma=arg[1]. 2516# The returned independent set will contain the (assumed) independent set 2517# arg[2] (default []) and not contain any element of arg[3] 2518# (default [], in which case the returned independent set is maximal). 2519# An error is signalled if arg[2] and arg[3] have non-trivial intersection. 2520# A "greedy" algorithm is used, and the graph gamma must be simple. 2521# 2522local gamma,is,forbidden,i,j,k,poss,adj,mindeg,minvert,degs; 2523gamma:=arg[1]; 2524if not IsBound(arg[2]) then 2525 is:=[]; 2526else 2527 is:=SSortedList(arg[2]); 2528fi; 2529if not IsBound(arg[3]) then 2530 forbidden:=[]; 2531else 2532 forbidden:=SSortedList(arg[3]); 2533fi; 2534if not (IsGraph(gamma) and IsSSortedList(is) and IsSSortedList(forbidden)) then 2535 Error("usage: IndependentSet( <Graph> [, <List> [, <List> ]] )"); 2536fi; 2537if not IsSimpleGraph(gamma) then 2538 Error("<gamma> not a simple graph"); 2539fi; 2540if Length(Intersection(is,forbidden))>0 then 2541 Error("<is> and <forbidden> have non-trivial intersection"); 2542fi; 2543if gamma.order=0 then 2544 return []; 2545fi; 2546poss:=Difference([1..gamma.order],forbidden); 2547SubtractSet(poss,is); 2548for i in is do 2549 SubtractSet(poss,Adjacency(gamma,i)); 2550od; 2551# is contains the independent set so far. 2552# poss contains the possible new elements of the independent set. 2553degs:=[]; 2554for i in poss do 2555 degs[i]:=Length(Intersection(Adjacency(gamma,i),poss)); 2556od; 2557while poss <> [] do 2558 minvert:=poss[1]; 2559 mindeg:=degs[poss[1]]; 2560 for i in [2..Length(poss)] do 2561 k:=degs[poss[i]]; 2562 if k < mindeg then 2563 mindeg:=k; 2564 minvert:=poss[i]; 2565 fi; 2566 od; 2567 AddSet(is,minvert); 2568 RemoveSet(poss,minvert); 2569 adj:=Intersection(Adjacency(gamma,minvert),poss); 2570 SubtractSet(poss,adj); 2571 for i in adj do 2572 for j in Intersection(poss,Adjacency(gamma,i)) do 2573 degs[j]:=degs[j]-1; 2574 od; 2575 od; 2576od; 2577return is; 2578end); 2579 2580BindGlobal("CollapsedIndependentOrbitsGraph",function(arg) 2581# 2582# Given a subgroup G=arg[1] of the automorphism group of 2583# the graph gamma=arg[2] (assumed simple), this function returns 2584# a graph delta defined as follows. The vertices of delta are 2585# those G-orbits of V(gamma) that are independent sets, 2586# and x is joined to y in delta iff x union y is *not* an 2587# independent set in gamma. 2588# If arg[3] is bound then it is assumed to be a subgroup of 2589# Aut(gamma) preserving the set of orbits of G on the vertices 2590# of gamma (for example, the normalizer of G in gamma.group). 2591# 2592local G,gamma,N,orb,orbs,i,j,L,rel; 2593G:=arg[1]; 2594gamma:=arg[2]; 2595if IsBound(arg[3]) then 2596 N:=arg[3]; 2597else 2598 N:=G; 2599fi; 2600if not IsPermGroup(G) or not IsGraph(gamma) or not IsPermGroup(N) then 2601 Error("usage: CollapsedIndependentOrbitsGraph( ", 2602 "<PermGroup>, <Graph> [, <PermGroup>] )"); 2603fi; 2604if not IsSimpleGraph(gamma) then 2605 Error("<gamma> must be a simple graph"); 2606fi; 2607orbs:=List(OrbitsDomain(G,[1..gamma.order]),Set); 2608i:=1; 2609L:=[]; 2610rel:=[]; 2611for orb in orbs do 2612 if Length(Intersection(orb,Adjacency(gamma,orb[1])))=0 then 2613 # an independent set is induced on orb 2614 Add(L,orb); 2615 for j in [i+1..i+Length(orb)-1] do 2616 Add(rel,[i,j]); 2617 od; 2618 i:=i+Length(orb); 2619 fi; 2620od; 2621return QuotientGraph(InducedSubgraph(gamma,Concatenation(L),N),rel); 2622end); 2623 2624BindGlobal("CollapsedCompleteOrbitsGraph",function(arg) 2625# 2626# Given a subgroup G=arg[1] of the automorphism group of 2627# the graph gamma=arg[2] (assumed simple), this function returns 2628# a graph delta defined as follows. The vertices of delta are 2629# those G-orbits of V(gamma) that are complete subgraphs, 2630# and x is joined to y in delta iff x<>y and x union y is a 2631# complete subgraph of gamma. 2632# If arg[3] is bound then it is assumed to be a subgroup of 2633# Aut(gamma) preserving the set of orbits of G on the vertices 2634# of gamma (for example, the normalizer of G in gamma.group). 2635# 2636local G,gamma,N; 2637G:=arg[1]; 2638gamma:=arg[2]; 2639if IsBound(arg[3]) then 2640 N:=arg[3]; 2641else 2642 N:=G; 2643fi; 2644if not IsPermGroup(G) or not IsGraph(gamma) or not IsPermGroup(N) then 2645 Error("usage: CollapsedCompleteOrbitsGraph( ", 2646 "<PermGroup>, <Graph> [, <PermGroup>] )"); 2647fi; 2648if not IsSimpleGraph(gamma) then 2649 Error("<gamma> not a simple graph"); 2650fi; 2651return 2652 ComplementGraph(CollapsedIndependentOrbitsGraph(G,ComplementGraph(gamma),N)); 2653end); 2654 2655BindGlobal("GraphImage",function(gamma,perm) 2656# 2657# Returns the image delta of the graph gamma, under the permutation perm, 2658# which must be a permutation of [1..gamma.order]. 2659# 2660local delta,i,perminv; 2661if not IsGraph(gamma) or not IsPerm(perm) then 2662 Error("usage: GraphImage( <Graph>, <Perm> )"); 2663fi; 2664if LargestMovedPoint(perm) > gamma.order then 2665 Error("<perm> must be a permutation of [1..<gamma>.order]"); 2666fi; 2667delta:=NullGraph(Group(List(GeneratorsOfGroup(gamma.group),x->x^perm),()), 2668 gamma.order); 2669if HasSize(gamma.group) or HasStabChainMutable(gamma.group) then 2670 SetSize(delta.group,Size(gamma.group)); 2671fi; 2672if IsBound(gamma.isSimple) then 2673 delta.isSimple:=gamma.isSimple; 2674else 2675 Unbind(delta.isSimple); 2676fi; 2677if IsBound(gamma.autGroup) then 2678 delta.autGroup:=Group(List(GeneratorsOfGroup(gamma.autGroup),x->x^perm),()); 2679 if HasSize(gamma.autGroup) or HasStabChainMutable(gamma.autGroup) then 2680 SetSize(delta.autGroup,Size(gamma.autGroup)); 2681 fi; 2682fi; 2683# if IsBound(gamma.canonicalLabelling) then 2684# delta.canonicalLabelling:=gamma.canonicalLabelling*perm; 2685# fi; 2686perminv:=perm^-1; 2687delta.names:=Immutable(List([1..delta.order],i->VertexName(gamma,i^perminv))); 2688if IsBound(gamma.maximumClique) then 2689 delta.maximumClique:=Immutable(OnSets(gamma.maximumClique,perm)); 2690fi; 2691if IsBound(gamma.minimumVertexColouring) then 2692 delta.minimumVertexColouring:=Immutable(List([1..delta.order], 2693 i->gamma.minimumVertexColouring[i^perminv])); 2694fi; 2695for i in [1..Length(delta.representatives)] do 2696 delta.adjacencies[i]:= 2697 OnSets(Adjacency(gamma,delta.representatives[i]^perminv),perm); 2698od; 2699return delta; 2700end); 2701 2702BindGlobal("CompleteSubgraphsMain",function(gamma,kvector,allsubs,allmaxes, 2703 partialcolour,weightvectors,dovector) 2704# 2705# This function, not for the user, subsumes the tasks formerly 2706# done by CompleteSubgraphs and CompleteSubgraphsOfGivenSize. 2707# These latter functions are now shells to check parameters and 2708# to call this main function, which can also compute complete 2709# subgraphs with given vertex-weightvector sum. More precise details 2710# are given below. 2711# 2712# Let gamma be a simple graph, and kvector an integer vector 2713# of dimension (i.e. a dense integer list of length) d>=1, with all 2714# entries non-negative if d>1. 2715# 2716# The parameter weightvectors is then a list of length 2717# gamma.order of non-zero d-vectors of non-negative integers, 2718# with the *weight-vector* (or *weightvector*) of a 2719# vertex v of gamma being weightvectors[v] 2720# and the *weight* of v being the sum of the entries of its 2721# weight-vector. Moreover, we assume that gamma.group acts on the 2722# set of all integer d-vectors by permuting vector positions, such that, 2723# for all v in [1..gamma.order] and g in gamma.group, we have 2724# 2725# weightvectors[v^g] = weightvectors[v]^g 2726# 2727# (where the first action is OnPoints, and the second action is the 2728# assumed one on integer d-vectors). Note that, in particular, this implies 2729# that Set(weightvectors) is invariant under gamma.group, 2730# and that the weight of a vertex is constant over a gamma.group 2731# orbit of vertices. 2732# 2733# We assume that kvector is fixed by gamma.group, and let k=Sum(kvector). 2734# 2735# First, suppose that kvector=[k], with k<0. 2736# 2737# Then this function returns a set K of complete subgraphs of gamma, 2738# with a complete subgraph being represented by the set of its vertices. 2739# If allmaxes=true then only maximal complete subgraphs are returned, 2740# and if allmaxes=false then arbitrary complete subgraphs are returned. 2741# 2742# The parameter allsubs is used to control how many 2743# complete subgraphs are returned. 2744# If allsubs=1, then K will contain a set of gamma.group 2745# orbit-representatives of the maximal (if allmaxes=true) 2746# or of all (if allmaxes=false) complete subgraphs of gamma. 2747# If allsubs=2 then K will be (exactly) a set of gamma.group 2748# orbit-representatives of the maximal (if allmaxes=true) or all 2749# (if allmaxes=false) complete subgraphs of gamma (this is usually 2750# more expensive than when allsubs=1). 2751# If allsubs=0, then K will contain exactly one complete subgraph, 2752# which is guaranteed to be maximal if allmaxes=true. 2753# 2754# The parameters partialcolour, weightvectors and dovector 2755# are ignored if k<0. 2756# 2757# Now suppose that k>=0. 2758# 2759# Then this function returns a set K of complete subgraphs of gamma, 2760# each of which having vertex-weightvectors summing to kvector. 2761# Such a complete subgraph is called a *solution* here. 2762# A complete subgraph is represented by the set of its vertices. 2763# Note that the set of all solutions is gamma.group-invariant. 2764# 2765# If allmaxes=true then only maximal complete subgraphs 2766# forming solutions are returned, and if allmaxes=false then 2767# the returned solutions need not be maximal complete subgraphs. 2768# 2769# The parameter allsubs is used to control how many solutions 2770# are returned. If allsubs=1, 2771# then K will contain (perhaps properly) a set of gamma.group 2772# orbit-representatives of all the solutions (if allmaxes=false) or 2773# of the solutions that form maximal complete subgraphs (if allmaxes=true). 2774# If allsubs=2 then K will be (exactly) a set of gamma.group 2775# orbit-representatives of all the solutions (if allmaxes=false) or 2776# of the solutions that form maximal complete subgraphs (if allmaxes=true) 2777# (this is usually more expensive than when allsubs=1). 2778# If allsubs=0, then K will contain at most one element and will 2779# contain one element iff gamma contains a solution, unless 2780# allmaxes=true, in which case K will contain one element iff gamma 2781# contains a solution which forms a maximal complete subgraph 2782# (in which case K will contain such a solution). 2783# 2784# The boolean parameter partialcolour determines whether 2785# or not partial proper vertex-colouring is used to try to cut 2786# down the search tree. The default is true (employ this vertex-colouring). 2787# 2788# The parameter dovector must be a d-vector containing an ordering 2789# of [1..d]. There is no harm (except perhaps for efficiency) in 2790# giving dovector the value [1..d]. 2791# 2792local IsFixedPoint,HasLargerEntry,k,smallorder,weights,weighted, 2793 originalG,originalgamma,includingallmaximalreps, 2794 CompleteSubgraphsSearch,K,clique,cliquenumber,chromaticnumber; 2795 2796IsFixedPoint := function(G,point) 2797# 2798# This boolean function returns true iff point is a fixed-point of the 2799# group G in its action OnPoints. 2800# 2801return ForAll(GeneratorsOfGroup(G),x->point^x=point); 2802end; 2803 2804HasLargerEntry := function(v,w) 2805# 2806# This boolean function assumes that v and w are equal-length integer vectors, 2807# and returns true iff v[i]>w[i] for some 1<=i<=Length(v). 2808# 2809local i; 2810for i in [1..Length(v)] do 2811 if v[i]>w[i] then 2812 return true; 2813 fi; 2814od; 2815return false; 2816end; 2817 2818CompleteSubgraphsSearch := function(gamma,kvector,sofar,forbidden) 2819# 2820# This recursive function is called by CompleteSubgraphsMain to do all 2821# the real work. It is assumed that on the call from CompleteSubgraphsMain, 2822# gammma.names is bound, and is equal to [1..gamma.order]. 2823# 2824# This function returns a dense list of distinct complete subgraphs of 2825# gamma, each of which is given as a dense list of distinct vertex-names. 2826# 2827# The variables smallorder, originalG, allsubs, allmaxes, weights, 2828# weightvectors, weighted, partialcolour, dovector, IsFixedPoint, and 2829# HasLargerEntry are global. (originalG is the group of automorphisms 2830# associated with the original graph.) 2831# 2832# If allsubs=2 then the returned complete subgraphs will be 2833# (pairwise) inequivalent under gamma.group. 2834# 2835# The parameter sofar is the set of vertices of the original graph 2836# chosen "so far". 2837# 2838# The parameter forbidden is a gamma.group-invariant set of 2839# vertex-names, none of which is allowed to be in any returned 2840# complete subgraph. The value of forbidden may be changed by 2841# this function. 2842# 2843# If allsubs>0 then this function returns a list of complete 2844# subgraphs which, when each of its elements is augmented by 2845# the elements in sofar, is a list of solutions containing 2846# a transversal of the distinct originalG-orbits of all the originally 2847# required solutions (maximal complete or otherwise) which contain sofar, 2848# except possibly for those orbits containing a complete subgraph 2849# which contains sofar and an element in forbidden. 2850# If allsubs=0 then this function returns (a list of) 2851# at most one complete subgraph, and returns (a list of) one 2852# complete subgraph c if and only if c union sofar is a solution 2853# (and also a maximal complete subgraph if allmaxes=true). 2854# 2855# If allsubs=2 or allmaxes: 2856# It is assumed that the set of vertex-names of gamma is the set 2857# of all vertices in the original graph adjacent to each element of sofar. 2858# It is also assumed that gamma.group (in its action on gamma.names) 2859# is contained in the image of the stabilizer in originalG of the set 2860# sofar (in that group's action on gamma.names), with equality if 2861# allsubs=2. 2862# 2863# At each call, we will attempt possible ways of adding a 2864# vertex(-name) v to sofar, such that 2865# weightvectors[names[v]][doposition]<>0, 2866# where doposition:=First(dovector,x->kvector[x])<>0); 2867# 2868# We "dynamically" order the search tree, as described in: 2869# W. Myrvold, T. Prsa and N. Walker, A Dynamic programming approach 2870# for timing and designing clique algorithms, Algorithms and Experiments 2871# (ALEX '98): Building Bridges Between Theory and Applications, 1998, 2872# pp. 88-95, 1998. 2873# 2874local k,n,i,j,delta,adj,rep,a,b,ans,ans1,ans2,names,W,H,HH,newsofar, 2875 G,orb,kk,ll,mm,active,nadj,verticesremoved,J,doposition, 2876 A,nactive,nactivevector,wt,indorbwtsum,CompleteSubgraphsSearch1; 2877 2878CompleteSubgraphsSearch1 := function(mask,kvector,forbidmask) 2879# 2880# This function does the work of CompleteSubgraphsSearch 2881# when the group associated with the graph is trivial. 2882# The parameters mask and forbidmask are boolean lists of length n, 2883# and the global variable A has value an n x n adjacency matrix 2884# for a graph gamma (given as a length n list of boolean lists). 2885# The actual graph in which we are determining complete subgraphs 2886# is the induced subgraph delta of gamma on the 2887# vertices i for which mask[i]=true, and we are determining 2888# complete subgraphs of delta containing no vertex j 2889# for which forbidmask[j]=true. We assume that (as sets), 2890# forbidmask is a subset of mask, and that if allmaxes=false then 2891# forbidmask is the empty set. 2892# 2893# The variables n, A, names, allmaxes, allsubs, partialcolour, 2894# weights, weightvectors, weighted, dovector, and HasLargerEntry 2895# are global. 2896# The parameter mask may be changed by this function, and if 2897# allmaxes=true then forbidmask may be changed by this function. 2898# 2899local k,active,activemask,a,b,c,col,verticesremoved,i,j,ans,ans1,kk,ll,mm, 2900 vertices,nactive,nactivevector,wt,wtvector,cw,cwsum,endconsider,nadj, 2901 doposition,minptr; 2902 2903activemask:=DifferenceBlist(mask,forbidmask); 2904active:=ListBlist([1..n],activemask); 2905IsSSortedList(active); 2906k:=Sum(kvector); 2907if k=0 or (k<0 and active=[]) then 2908 if allmaxes and SizeBlist(mask)>0 then 2909 # We are only looking for maximal complete subgraphs, 2910 # but here the complete subgraph of size 0 is *not* maximal. 2911 return []; 2912 else 2913 # allmaxes=false or there are no vertices 2914 return [[]]; 2915 fi; 2916fi; 2917nactive:=Sum(weights{names{active}}); 2918if nactive<k then 2919 # in particular, if nactive=0 then 2920 return []; 2921fi; 2922nactivevector:=ShallowCopy(Sum(weightvectors{names{active}})); 2923# (ShallowCopy since the value of nactivevecter may be changed) 2924if HasLargerEntry(kvector,nactivevector) then 2925 return []; 2926fi; 2927# now we have nactive>0, and either k<0 or nactive >= k > 0. 2928vertices:=ListBlist([1..n],mask); 2929IsSSortedList(vertices); 2930repeat 2931 nadj := []; # nadj[j] will record the number of active vertices adjacent 2932 # to active[j]. 2933 verticesremoved := false; 2934 for j in [1..Length(active)] do 2935 i:=active[j]; 2936 b:=IntersectionBlist(activemask,A[i]); 2937 nadj[j]:=SizeBlist(b); 2938 if k>=0 then 2939 if weighted then 2940 ll:=Sum(weights{ListBlist(names,b)}); 2941 else 2942 ll:=nadj[j]; 2943 fi; 2944 wt:=weights[names[i]]; 2945 wtvector:=weightvectors[names[i]]; 2946 mm:=HasLargerEntry(wtvector,kvector); 2947 if ll+wt<k or (mm and (not allmaxes)) then 2948 # eliminate vertex i 2949 verticesremoved:=true; 2950 mask[i]:=false; 2951 forbidmask[i]:=false; 2952 activemask[i]:=false; 2953 nactive:=nactive-wt; 2954 AddRowVector(nactivevector,wtvector,-1); 2955 if HasLargerEntry(kvector,nactivevector) then 2956 return []; 2957 fi; 2958 elif mm then 2959 # allmaxes=true so we forbid vertex i, but don't eliminate it 2960 verticesremoved:=true; 2961 forbidmask[i]:=true; 2962 activemask[i]:=false; 2963 nactive:=nactive-wt; 2964 AddRowVector(nactivevector,wtvector,-1); 2965 if HasLargerEntry(kvector,nactivevector) then 2966 return []; 2967 fi; 2968 fi; 2969 fi; 2970 od; 2971 if verticesremoved then 2972 # Update active and vertices. 2973 # At this point we know that nactive>=k>0, and that no entry 2974 # of kvector is greater than the corresponding entry of nactivevector. 2975 active:=ListBlist([1..n],activemask); 2976 IsSSortedList(active); 2977 vertices:=ListBlist([1..n],mask); 2978 IsSSortedList(vertices); 2979 fi; 2980until not verticesremoved; 2981# At this point we know that k<0 or nactive>=k>0, and if k>0, that no 2982# entry of kvector is greater than the corresponding entry of nactivevector. 2983if nactive=k and Length(active)=Length(vertices) then 2984 # k>0, no forbidden vertices, and we are down to a required solution 2985 return [names{active}]; 2986fi; 2987kk:=Length(active)-1; 2988if ForAll(nadj,x->x=kk) then 2989 # The induced subgraph on the active vertices is a complete subgraph 2990 # with vertex-weight sum >= k (we may have k<0). 2991 if Length(vertices)>Length(active) and (k<0 or nactive=k) then 2992 # Possible solution, but at least one 2993 # vertex is forbidden (so also allmaxes=true). 2994 ll:=IntersectionBlist(IntersectionBlist(List(active,x->A[x])),mask); 2995 if SizeBlist(ll)=0 then 2996 # the complete subgraph induced on the active vertices is maximal 2997 return [names{active}]; 2998 else 2999 # the complete subgraph induced on the active vertices is not maximal 3000 return []; 3001 fi; 3002 elif k<0 then 3003 # no forbidden vertices here 3004 if allmaxes or allsubs=0 then 3005 return [names{active}]; 3006 else 3007 return Combinations(names{active}); 3008 fi; 3009 else 3010 # k>0 and nactive>k here, and so each maximal complete 3011 # subgraph of vertex-weight-sum k contains a forbidden vertex 3012 if allmaxes then 3013 return []; 3014 else 3015 if not weighted then 3016 if allsubs=0 then 3017 return [names{active{[1..k]}}]; 3018 else 3019 return Combinations(names{active},k); 3020 fi; 3021 fi; 3022 fi; 3023 fi; 3024fi; 3025endconsider:=Length(active); 3026# 3027# Now order the vertices in active for processing. 3028# 3029# Begin by pushing ignorable active indices beyond endconsider. 3030# 3031doposition:=First(dovector,x->kvector[x]<>0); 3032if (allmaxes and Length(kvector)=1) or Length(kvector)>1 then 3033 if Length(kvector)=1 then 3034 # allmaxes=true here 3035 mm:=Difference(vertices,active); 3036 if mm=[] then 3037 mm:=A[active[1]]; 3038 else 3039 mm:=A[mm[1]]; 3040 fi; 3041 fi; 3042 i:=1; 3043 while i<=endconsider do 3044 if (Length(kvector)=1 and mm[active[i]]) or (Length(kvector)>1 and 3045 weightvectors[names[active[i]]][doposition]=0) then 3046 if i<endconsider then 3047 a:=active[endconsider]; 3048 active[endconsider]:=active[i]; 3049 active[i]:=a; 3050 nadj[i]:=nadj[endconsider]; 3051 fi; 3052 endconsider:=endconsider-1; 3053 else 3054 i:=i+1; 3055 fi; 3056 od; 3057fi; 3058# 3059# Now order the elements in active{[1..endconsider]}, and the corresponding 3060# elements of nadj. 3061# 3062for i in [1..endconsider] do 3063 minptr:=i; 3064 for j in [i+1..endconsider] do 3065 if nadj[j]<nadj[minptr] then 3066 minptr:=j; 3067 fi; 3068 od; 3069 a:=active[i]; 3070 active[i]:=active[minptr]; 3071 active[minptr]:=a; 3072 a:=nadj[i]; 3073 nadj[i]:=nadj[minptr]; 3074 nadj[minptr]:=a; 3075 mm:=A[active[i]]; 3076 for j in [i+1..endconsider] do 3077 if mm[active[j]] then 3078 nadj[j]:=nadj[j]-1; 3079 fi; 3080 od; 3081od; 3082if k>=0 and partialcolour then 3083 # We do (perhaps partial) proper vertex-colouring. 3084 col:=[]; 3085 cw:=[]; # cw[j] will record the largest weight (entry) in the doposition of 3086 # (a weightvector of) a vertex having colour j 3087 cwsum:=0; 3088 mm:=0; # max. colour used so far 3089 if Length(kvector)>1 then 3090 ll:=endconsider; 3091 else 3092 ll:=Length(active); 3093 fi; 3094 for i in Reversed([1..ll]) do # col[i] := colour of active[i] 3095 c:=BlistList([1..mm+1],[]); 3096 b:=A[active[i]]; 3097 for j in [i+1..ll] do 3098 if b[active[j]] then 3099 c[col[j]]:=true; 3100 fi; 3101 od; 3102 j:=1; 3103 while c[j] do 3104 j:=j+1; 3105 od; 3106 col[i]:=j; 3107 wt:=weightvectors[names[active[i]]][doposition]; 3108 if j>mm then 3109 mm:=j; 3110 cwsum:=cwsum+wt; 3111 cw[mm]:=wt; 3112 elif cw[j]<wt then 3113 cwsum:=cwsum+(wt-cw[j]); 3114 cw[j]:=wt; 3115 fi; 3116 if cwsum>=kvector[doposition] then 3117 # stop colouring 3118 if endconsider>i then 3119 endconsider:=i; 3120 fi; 3121 break; 3122 fi; 3123 od; 3124 if cwsum < kvector[doposition] then 3125 # there is no solution 3126 return []; 3127 fi; 3128fi; 3129if k<0 and (not allmaxes) then 3130 ans:=[[]]; 3131 if allsubs=0 then 3132 return ans; 3133 fi; 3134else 3135 ans:=[]; 3136fi; 3137for i in active{[1..endconsider]} do 3138 wtvector:=weightvectors[names[i]]; 3139 ans1:=CompleteSubgraphsSearch1(IntersectionBlist(mask,A[i]), 3140 kvector-wtvector, 3141 IntersectionBlist(forbidmask,A[i])); 3142 if Length(ans1)>0 then 3143 for a in ans1 do 3144 Add(a,names[i]); 3145 Add(ans,a); 3146 od; 3147 if allsubs=0 then 3148 return ans; 3149 fi; 3150 fi; 3151 AddRowVector(nactivevector,wtvector,-1); 3152 if HasLargerEntry(kvector,nactivevector) then 3153 break; 3154 fi; 3155 if allmaxes then 3156 forbidmask[i]:=true; 3157 else 3158 mask[i]:=false; 3159 fi; 3160od; 3161return ans; 3162end; 3163 3164# 3165# begin CompleteSubgraphsSearch 3166# 3167k:=Sum(kvector); 3168n:=gamma.order; 3169if k=0 or (k<0 and n=Length(forbidden)) then 3170 if allmaxes and n>0 then 3171 # We are only looking for maximal complete subgraphs, 3172 # but here the complete subgraph of size 0 is *not* maximal. 3173 return []; 3174 else 3175 # allmaxes=false or there are no vertices 3176 return [[]]; 3177 fi; 3178fi; 3179names:=gamma.names; 3180active:=Filtered([1..n],x->not (names[x] in forbidden)); 3181IsSSortedList(active); 3182nactive:=Sum(weights{names{active}}); 3183if nactive<k then 3184 # in particular, if nactive=0 then 3185 return []; 3186fi; 3187nactivevector:=ShallowCopy(Sum(weightvectors{names{active}})); 3188# (ShallowCopy since the value of nactivevecter may be changed) 3189if HasLargerEntry(kvector,nactivevector) then 3190 return []; 3191fi; 3192# now k<0 or nactive >= k > 0. 3193G:=gamma.group; 3194if IsTrivial(G) then 3195 # The group will be trivial from here on, and we will use the 3196 # specialized function CompleteSubgraphsSearch1. 3197 if (not allmaxes) and forbidden<>[] then 3198 # strip out the forbidden vertices. 3199 gamma:=InducedSubgraph(gamma,active,G); 3200 n:=gamma.order; 3201 names:=gamma.names; 3202 active:=[1..n]; 3203 fi; 3204 # now A := adjacency matrix of gamma. 3205 A:=List([1..n],i->BlistList([1..n],gamma.adjacencies[i])); 3206 return CompleteSubgraphsSearch1(BlistList([1..n],[1..n]), kvector, 3207 BlistList([1..n],Difference([1..n],active))); 3208fi; 3209J:=Filtered([1..Length(gamma.representatives)], 3210 x->gamma.representatives[x] in active); 3211IsSSortedList(J); 3212nadj:=[]; # nadj[i] will store the number of active vertices adjacent 3213 # to gamma.representatives[J[i]] 3214verticesremoved:=false; 3215for i in [1..Length(J)] do 3216 rep:=gamma.representatives[J[i]]; 3217 a:=gamma.adjacencies[J[i]]; 3218 if forbidden<>[] then 3219 a:=Filtered(a,x->not (names[x] in forbidden)); 3220 fi; 3221 nadj[i]:=Length(a); 3222 if k>=0 then 3223 if weighted then 3224 ll:=Sum(weights{names{a}}); 3225 else 3226 ll:=nadj[i]; 3227 fi; 3228 if ll+weights[names[rep]] < k 3229 or HasLargerEntry(weightvectors[names[rep]],kvector) then 3230 # forbid the vertex-orbit containing rep 3231 verticesremoved:=true; 3232 UniteSet(forbidden,names{Orbit(G,rep)}); 3233 fi; 3234 fi; 3235od; 3236if verticesremoved then 3237 return CompleteSubgraphsSearch(gamma,kvector,sofar,forbidden); 3238fi; 3239# At this point, active is the set of non-forbidden vertices, 3240# and k<0 or nactive>=k>0. Moreover, if k>0, then no 3241# entry of kvector is greater than the corresponding entry of nactivevector. 3242if nactive=k and (not allmaxes or forbidden=[]) then 3243 # k>0 and we are down to a complete graph on active 3244 # vertices which forms a required solution. 3245 return [names{active}]; 3246fi; 3247kk:=Length(active)-1; 3248if ForAll(nadj,x->x=kk) then 3249 # The induced subgraph on the active vertices is a complete subgraph 3250 # with vertex-weight sum >= k (we may have k<0). 3251 if allmaxes then 3252 if k>=0 and nactive>k then 3253 # any maximal complete solution subgraph must contain 3254 # a forbidden vertex 3255 return []; 3256 else 3257 # k<0 or nactivevector=kvector. 3258 # We now check whether any forbidden vertex in 3259 # gamma.representatives is joined to all active vertices. 3260 for i in Difference([1..Length(gamma.representatives)],J) do 3261 if IsSubset(gamma.adjacencies[i],active) then 3262 # Each required maximal complete subgraph contains a 3263 # forbidden vertex. 3264 return []; 3265 fi; 3266 od; 3267 # At this point we know that the complete subgraph induced on the 3268 # active vertices is maximal. 3269 return [names{active}]; 3270 fi; 3271 fi; 3272 # at this point, allmaxes=false and (nactive>k or k<0). 3273 if allsubs=0 then 3274 if k<0 then 3275 return [names{active}]; 3276 elif not weighted then 3277 return [names{active{[1..k]}}]; 3278 fi; 3279 fi; 3280fi; 3281if allsubs=2 then 3282 # set up translation vector W from vertex-names to vertices (of gamma). 3283 W:=[]; 3284 for i in [1..Length(names)] do 3285 W[names[i]]:=i; 3286 od; 3287fi; 3288# next initialize ans 3289if k<0 and (not allmaxes) then 3290 ans:=[[]]; 3291 if allsubs=0 then 3292 return ans; 3293 fi; 3294else 3295 ans:=[]; 3296fi; 3297if allmaxes and Length(kvector)=1 then 3298 mm:=First(Difference([1..Length(gamma.adjacencies)],J), 3299 x->IsFixedPoint(gamma.group,gamma.representatives[x])); 3300 if mm<>fail then 3301 mm:=gamma.adjacencies[mm]; 3302 else 3303 mm:=First(J,x->IsFixedPoint(gamma.group,gamma.representatives[x])); 3304 if mm<>fail then 3305 mm:=gamma.adjacencies[mm]; 3306 else 3307 mm:=[]; 3308 fi; 3309 fi; 3310 if mm<>[] then 3311 # 3312 # We can ignore the elements of the gamma.group-invariant set mm, 3313 # since Length(kvector)=1, allmaxes=true and mm is the 3314 # adjacency of a single (heuristically) chosen vertex. 3315 # 3316 ll:=Filtered([1..Length(J)],x->not (gamma.representatives[J[x]] in mm)); 3317 J:=J{ll}; 3318 nadj:=nadj{ll}; 3319 fi; 3320else 3321 mm:=[]; 3322fi; 3323indorbwtsum:=0; 3324doposition:=First(dovector,x->kvector[x]<>0); 3325for j in [1..Length(J)] do 3326 i:=1; 3327 for kk in [2..Length(J)] do 3328 if nadj[kk]<nadj[i] then 3329 i:=kk; 3330 fi; 3331 od; 3332 nadj[i]:=n; 3333 rep:=gamma.representatives[J[i]]; 3334 adj:=gamma.adjacencies[J[i]]; 3335 orb:=SSortedList(Orbit(G,rep)); 3336 # wt will record the maximum entry in doposition of a vector in orb. 3337 if Length(kvector)=1 then 3338 wt:=weightvectors[names[rep]][1]; # does not depend on orbit rep. 3339 else 3340 wt:=0; 3341 for a in orb do 3342 kk:=weightvectors[names[a]][doposition]; 3343 if kk>wt or (a=rep and kk=wt) then 3344 # these statements are executed at least once 3345 wt:=kk; 3346 b:=a; 3347 fi; 3348 od; 3349 if b<>rep then 3350 rep:=b; 3351 adj:=Adjacency(gamma,rep); 3352 fi; 3353 fi; 3354 if wt<>0 then 3355 # 3356 # We consider searching for solutions containing rep. 3357 # 3358 # However, if k>=0 and mm=[] we shall ignore a 3359 # set of independent orbits, such that the sum 3360 # of the wt's for these orbits is less than kvector[doposition]. 3361 # This is because no solution can be made from vertices coming 3362 # only from these independent orbits. 3363 # 3364 if k>=0 and Length(mm)=0 and indorbwtsum+wt < kvector[doposition] 3365 and Length(Intersection(orb,adj))=0 then 3366 # 3367 # Ignore the independent (active) orbit orb, and add wt to the 3368 # running total indorbwtsum. 3369 # 3370 indorbwtsum:=indorbwtsum+wt; 3371 else 3372 newsofar:=Union(sofar,[names[rep]]); 3373 if allsubs<>2 and (not allmaxes) then 3374 # We can strip out all forbidden vertices since in this case: 3375 # (1) we are stabilizing each successive sofar with 3376 # (a constituent image of) a subgroup of 3377 # the previous stabilizer of sofar, and 3378 # so the set of non-forbidden vertices will *always* be 3379 # invariant under further gamma.groups; 3380 # (2) we need not check whether our complete subgraphs are maximal 3381 delta:=InducedSubgraph(gamma, 3382 Filtered(adj,x->not (names[x] in forbidden)), 3383 ProbablyStabilizer(gamma.group,rep)); 3384 ans1:=CompleteSubgraphsSearch(delta, 3385 kvector-weightvectors[names[rep]],newsofar,[]); 3386 elif allsubs<>2 then 3387 delta:=InducedSubgraph(gamma,adj, 3388 ProbablyStabilizer(gamma.group,rep)); 3389 ans1:=CompleteSubgraphsSearch(delta, 3390 kvector-weightvectors[names[rep]],newsofar, 3391 Intersection(delta.names,forbidden)); 3392 else 3393 # allsubs=2 3394 delta:=InducedSubgraph(gamma,adj,Stabilizer(gamma.group,rep)); 3395 HH:=Stabilizer(originalG,newsofar,OnSets); 3396 if not IsFixedPoint(HH,names[rep]) then 3397 H:=Action(HH,names{adj},OnPoints); 3398 delta:=NewGroupGraph(H,delta); 3399 ans1:=CompleteSubgraphsSearch(delta, 3400 kvector-weightvectors[names[rep]],newsofar, 3401 Intersection(delta.names,Union(Orbits(HH,forbidden)))); 3402 else 3403 ans1:=CompleteSubgraphsSearch(delta, 3404 kvector-weightvectors[names[rep]],newsofar, 3405 Intersection(delta.names,forbidden)); 3406 fi; 3407 fi; 3408 if Length(ans1)>0 then 3409 for a in ans1 do 3410 Add(a,names[rep]); 3411 od; 3412 if allsubs=0 then 3413 return ans1; 3414 fi; 3415 if allsubs<>2 or Length(ans1)=1 or IsFixedPoint(gamma.group,rep) then 3416 # isomorph rejection is unnecessary 3417 for a in ans1 do 3418 Add(ans,a); 3419 od; 3420 elif Size(gamma.group)<=smallorder then 3421 # perform isomorph rejection using explicit orbits 3422 ans1:=List(ans1,x->Set(W{x})); 3423 ans1:=List(Orbits(gamma.group,ans1,OnSets),x->names{x[1]}); 3424 for a in ans1 do 3425 Add(ans,a); 3426 od; 3427 else 3428 # perform isomorph rejection using SmallestImageSet 3429 ans2:=List(ans1,x-> 3430 SmallestImageSet(gamma.group,Set(W{x}))); 3431 SortParallel(ans2,ans1); 3432 Add(ans,ans1[1]); 3433 for a in [2..Length(ans1)] do 3434 if ans2[a]<>ans2[a-1] then 3435 # new gamma.group orbit representative 3436 Add(ans,ans1[a]); 3437 fi; 3438 od; 3439 fi; 3440 fi; 3441 if j < Length(J) then 3442 AddRowVector(nactivevector,Sum(weightvectors{names{orb}}),-1); 3443 if HasLargerEntry(kvector,nactivevector) then 3444 break; 3445 fi; 3446 for kk in [1..Length(J)] do 3447 if nadj[kk]<>n then 3448 adj:=gamma.adjacencies[J[kk]]; 3449 for a in orb do 3450 if a in adj then 3451 nadj[kk]:=nadj[kk]-1; 3452 fi; 3453 od; 3454 fi; 3455 od; 3456 UniteSet(forbidden,names{orb}); 3457 fi; 3458 fi; 3459 fi; 3460od; 3461return ans; 3462end; 3463 3464# 3465# begin CompleteSubgraphsMain 3466# 3467# Minimal checking of parameters since this function should only 3468# be called internally or by experts. 3469# 3470if not (IsGraph(gamma) and IsList(kvector) and IsInt(allsubs) and 3471 IsBool(allmaxes) and IsBool(partialcolour) and 3472 IsList(weightvectors) and IsList(dovector)) then 3473 Error("usage: CompleteSubgraphsMain( <Graph>, <List>, <Int>, <Bool>, <Bool>, <List>, <List>)"); 3474fi; 3475if not IsSimpleGraph(gamma) then 3476 Error("<gamma> must be a simple graph"); 3477fi; 3478smallorder:=8; # Any group with order <= smallorder is considered small, 3479 # and we calculate orbits on cliques explicitly for these 3480 # groups when allsubs=2. 3481originalgamma:=gamma; 3482originalG:=gamma.group; 3483gamma:=ShallowCopy(gamma); 3484gamma.names:=Immutable([1..gamma.order]); 3485k:=Sum(kvector); 3486if k<0 then 3487 # We are computing maximal complete subgraphs (not of given size). 3488 if Length(kvector)<>1 then 3489 Error("cannot have Sum(<kvector>)<0 if Length(<kvector>)<>1"); 3490 fi; 3491 includingallmaximalreps:=(allsubs in [1,2]); 3492 partialcolour:=false; 3493 weightvectors:=List([1..gamma.order],x->[1]); 3494 weights:=ListWithIdenticalEntries(gamma.order,1); 3495 weighted:=false; 3496 dovector:=[1]; 3497else 3498 includingallmaximalreps:=false; 3499 weights:=List(weightvectors,x->Sum(x)); 3500 weighted:=not ForAll(weightvectors,x->x=[1]); 3501 if weighted and ForAll(weights,x->x=weights[1]) 3502 and k mod weights[1] <> 0 then 3503 # there is no solution 3504 return []; 3505 fi; 3506fi; 3507if not weighted and k>=0 then 3508 if IsBound(gamma.maximumClique) then 3509 cliquenumber:=Length(gamma.maximumClique); 3510 if k>cliquenumber then 3511 # gamma has no clique of size k. 3512 return []; 3513 fi; 3514 if allsubs=0 and (not allmaxes or k=cliquenumber) then 3515 return [gamma.maximumClique{[1..k]}]; 3516 fi; 3517 elif IsBound(gamma.minimumVertexColouring) then 3518 chromaticnumber:=Length(Set(gamma.minimumVertexColouring)); 3519 if k>chromaticnumber then 3520 # gamma has no clique of size k. 3521 return []; 3522 fi; 3523 fi; 3524fi; 3525K:=CompleteSubgraphsSearch(gamma,kvector,[],[]); 3526for clique in K do 3527 Sort(clique); 3528od; 3529Sort(K); 3530if not weighted and not IsBound(originalgamma.maximumClique) then 3531 if includingallmaximalreps then 3532 # K contains a maximum clique of originalgamma. 3533 cliquenumber:=Maximum(List(K,Length)); 3534 originalgamma.maximumClique:=Immutable(First(K,x->Length(x)=cliquenumber)); 3535 elif IsBound(originalgamma.minimumVertexColouring) then 3536 chromaticnumber:=Length(Set(originalgamma.minimumVertexColouring)); 3537 if ForAny(K,x->Length(x)=chromaticnumber) then 3538 cliquenumber:=chromaticnumber; 3539 originalgamma.maximumClique:=Immutable(First(K,x->Length(x)=cliquenumber)); 3540 fi; 3541 fi; 3542fi; 3543return K; 3544end); 3545 3546BindGlobal("CompleteSubgraphsOfGivenSize",function(arg) 3547# 3548# Interface to CompleteSubgraphsMain. 3549# 3550local gamma,k,kvector,allsubs,allmaxes,partialcolour,weights,weightvectors; 3551if not (Length(arg) in [2..6]) then 3552 Error("must have 2, 3, 4, 5 or 6 parameters"); 3553fi; 3554gamma:=arg[1]; 3555k:=arg[2]; 3556if IsBound(arg[3]) then 3557 allsubs:=arg[3]; 3558else 3559 allsubs:=1; 3560fi; 3561if allsubs=false then 3562 allsubs:=0; 3563elif allsubs=true then 3564 allsubs:=1; 3565elif not (allsubs in [0,1,2]) then 3566 Error("<arg[3]> must be boolean or in [0,1,2]"); 3567fi; 3568if IsBound(arg[4]) then 3569 allmaxes:=arg[4]; 3570else 3571 allmaxes:=false; 3572fi; 3573if IsBound(arg[5]) then 3574 partialcolour:=arg[5]; 3575else 3576 partialcolour:=true; 3577fi; 3578if IsRat(partialcolour) then 3579 partialcolour:=true; # for backward compatibility 3580fi; 3581if not (IsGraph(gamma) and (IsInt(k) or IsList(k)) and IsBool(allmaxes) 3582 and IsBool(partialcolour) 3583 and (not IsBound(arg[6]) or IsList(arg[6])) ) then 3584 Error("usage: CompleteSubgraphsOfGivenSize( <Graph>, <Int> or <List>", 3585 "[, <Int> or <Bool> [, <Bool> [, <Bool> or <Rat> [, <List> ]]]] )"); 3586fi; 3587if IsInt(k) then 3588 if k<0 then 3589 Error("<k> must be non-negative"); 3590 fi; 3591 kvector:=[k]; 3592else 3593 kvector:=k; 3594 if Length(kvector)=0 or ForAny(kvector,x->x<0) then 3595 Error("<kvector> must be a non-empty list of non-negative integers"); 3596 fi; 3597fi; 3598if not IsSimpleGraph(gamma) then 3599 Error("<gamma> not a simple graph"); 3600fi; 3601if IsBound(arg[6]) then 3602 weights:=arg[6]; 3603 if Length(weights)<>gamma.order then 3604 Error("<weights> not of length <gamma>.order"); 3605 fi; 3606 if Length(weights)>0 and IsInt(weights[1]) then 3607 if ForAny(weights,w->not IsInt(w) or w<=0) then 3608 Error("all weights must be positive integers (or all lists)"); 3609 fi; 3610 if ForAny(GeneratorsOfGroup(gamma.group),g-> 3611 ForAny([1..gamma.order],i->weights[i^g]<>weights[i])) then 3612 Error("integer vertex-weights not <gamma>.group-invariant"); 3613 fi; 3614 weightvectors:=List(weights,x->[x]); 3615 else 3616 weightvectors:=weights; 3617 fi; 3618else 3619 weightvectors:=List([1..gamma.order],x->[1]); 3620fi; 3621if ForAny(weightvectors,x->Length(x)<>Length(kvector)) then 3622 Error("All weight-vectors must be the same length as <kvector>"); 3623fi; 3624return CompleteSubgraphsMain(gamma,kvector,allsubs,allmaxes,partialcolour, 3625 weightvectors,[1..Length(kvector)]); 3626end); 3627 3628BindGlobal("CliquesOfGivenSize",CompleteSubgraphsOfGivenSize); 3629 3630BindGlobal("CompleteSubgraphs",function(arg) 3631# 3632# Interface to CompleteSubgraphsMain. 3633# 3634local gamma,k,allsubs,allmaxes; 3635if not (Length(arg) in [1..3]) then 3636 Error("must have 1, 2 or 3 parameters"); 3637fi; 3638gamma:=arg[1]; 3639if IsBound(arg[2]) then 3640 k:=arg[2]; 3641else 3642 k:=-1; 3643fi; 3644if IsBound(arg[3]) then 3645 allsubs:=arg[3]; 3646else 3647 allsubs:=1; 3648fi; 3649if allsubs=false then 3650 allsubs:=0; 3651elif allsubs=true then 3652 allsubs:=1; 3653elif not (allsubs in [0,1,2]) then 3654 Error("<arg[3]> must be boolean or in [0,1,2]"); 3655fi; 3656if k<0 then 3657 allmaxes:=true; 3658else 3659 allmaxes:=false; 3660fi; 3661if not (IsGraph(gamma) and IsInt(k)) then 3662 Error("usage: CompleteSubgraphs( <Graph> [,<Int> [,<Int> or <Bool> ]] )"); 3663fi; 3664if not IsSimpleGraph(gamma) then 3665 Error("<gamma> not a simple graph"); 3666fi; 3667return CompleteSubgraphsMain(gamma,[k],allsubs,allmaxes, 3668 true,List([1..gamma.order],x->[1]),[1]); 3669end); 3670 3671BindGlobal("Cliques",CompleteSubgraphs); 3672 3673BindGlobal("CayleyGraph",function(arg) 3674# 3675# Given a group G=arg[1] and a list gens=arg[2] of 3676# generators for G, this function constructs a Cayley graph 3677# for G w.r.t. the generators gens. The generating list 3678# arg[2] is optional, and if omitted, then we take 3679# gens:=GeneratorsOfGroup(G). 3680# The boolean argument arg[3] is also optional, and if true (the default) 3681# then the returned graph is undirected (as if gens was closed 3682# under inversion whether or not it is). 3683# 3684# The Cayley graph caygraph which is returned is defined as follows: 3685# the vertices (actually the vertex-names) of caygraph are the elements 3686# of G; if arg[3]=true (the default) then vertices x,y are 3687# joined by an edge iff there is a g in gens with y=g*x 3688# or y=g^-1*x; if arg[3]=false then vertices x,y are 3689# joined by an edge iff there is a g in gens with y=g*x. 3690# 3691# *Note* It is not checked whether G = <gens>. However, even if G 3692# is not generated by gens, the function still works as described 3693# above (as long as gens is contained in G), but returns a 3694# "Cayley graph" which is not connected. 3695# 3696local G,gens,elms,undirected,caygraph; 3697G:=arg[1]; 3698if IsBound(arg[2]) then 3699 gens:=arg[2]; 3700else 3701 gens:=GeneratorsOfGroup(G); 3702fi; 3703if IsBound(arg[3]) then 3704 undirected:=arg[3]; 3705else 3706 undirected:=true; 3707fi; 3708if not(IsGroup(G) and IsList(gens) and IsBool(undirected)) then 3709 Error("usage: CayleyGraph( <Group> [, <List> [, <Bool> ]] )"); 3710fi; 3711elms:=AsList(G); 3712SetSize(G,Length(elms)); 3713if not IsSSortedList(gens) then 3714 gens:=SSortedList(gens); 3715fi; 3716caygraph := Graph(G,elms,OnRight, 3717 function(x,y) return y*x^-1 in gens; end,true); 3718# 3719# Note that caygraph.group comes from the right regular action of 3720# G as a group of automorphisms of the Cayley graph constructed. 3721# 3722SetSize(caygraph.group,caygraph.order); 3723if undirected then 3724 caygraph:=UnderlyingGraph(caygraph); 3725fi; 3726return caygraph; 3727end); 3728 3729BindGlobal("SwitchedGraph",function(arg) 3730# 3731# Returns the switched graph delta of graph gamma=arg[1], 3732# w.r.t to vertex list (or singleton vertex) V=arg[2]. 3733# 3734# The returned graph delta has vertex-set the same as gamma. 3735# If vertices x,y of delta are both in V or both not in 3736# V, then [x,y] is an edge of delta iff [x,y] is an edge 3737# of gamma; otherwise [x,y] is an edge of delta iff [x,y] 3738# is not an edge of gamma. 3739# 3740# If arg[3] is bound then it is assumed to be a subgroup 3741# of Aut(gamma) stabilizing V setwise. 3742# 3743local gamma,delta,n,V,W,H,A,i; 3744gamma:=arg[1]; 3745V:=arg[2]; 3746if IsInt(V) then 3747 V:=[V]; 3748fi; 3749if not IsGraph(gamma) or not IsList(V) then 3750 Error("usage: SwitchedGraph( <Graph>, <Int> or <List>, [,<PermGroup>] )"); 3751fi; 3752n:=gamma.order; 3753V:=SSortedList(V); 3754if not IsSubset([1..n],V) then 3755 Error("<V> must be a subset of [1..<n>]"); 3756fi; 3757if Length(V) > n/2 then 3758 V:=Difference([1..n],V); 3759fi; 3760if V=[] then 3761 return CopyGraph(gamma); 3762fi; 3763if IsBound(arg[3]) then 3764 H:=arg[3]; 3765elif Length(V)=1 then 3766 H:=ProbablyStabilizer(gamma.group,V[1]); 3767else 3768 H:=Group([],()); 3769fi; 3770if not IsPermGroup(H) then 3771 Error("usage: SwitchedGraph( <Graph>, <Int> or <List>, [, <PermGroup>] )"); 3772fi; 3773delta:=NullGraph(H,n); 3774if IsBound(gamma.isSimple) then 3775 delta.isSimple:=gamma.isSimple; 3776else 3777 Unbind(delta.isSimple); 3778fi; 3779if IsBound(gamma.names) then 3780 delta.names:=Immutable(gamma.names); 3781fi; 3782W:=Difference([1..n],V); 3783for i in [1..Length(delta.representatives)] do 3784 A:=Adjacency(gamma,delta.representatives[i]); 3785 if delta.representatives[i] in V then 3786 delta.adjacencies[i]:=Union(Intersection(A,V),Difference(W,A)); 3787 else 3788 delta.adjacencies[i]:=Union(Intersection(A,W),Difference(V,A)); 3789 fi; 3790od; 3791return delta; 3792end); 3793 3794BindGlobal("VertexTransitiveDRGs",function(gpin) 3795# 3796# If gpin is a permutation group G, then it must be transitive 3797# and non-trivial, and we set coladjmats:=OrbitalGraphColadjMats(gpin). 3798# 3799# Otherwise, we take coladjmats:=gpin, which must be a list of collapsed 3800# adjacency matrices for the orbital digraphs of a non-trivial 3801# transitive permutation group G (on a set V say), collapsed 3802# w.r.t. a point stabilizer (such as the list of matrices produced by 3803# OrbitalGraphColadjMats ). 3804# 3805# In either case, this function returns a record (called result), 3806# which gives information on G. 3807# The most important component of this record is the list 3808# orbitalCombinations, whose elements give the combinations of 3809# the (indices of) the G-orbitals whose union gives the edge-set 3810# of a distance-regular graph with vertex-set V. 3811# The component intersectionArrays gives the corresponding 3812# intersection arrays. The component degree is the degree of 3813# the permutation group G, rank is its (permutation) rank, and 3814# isPrimitive is true if G is primitive, and false otherwise. 3815# It is assumed that the orbital/suborbit indexing used is the same 3816# as that for the rows (and columns) of each of the matrices and 3817# also for the indexing of the matrices themselves, with the trivial 3818# suborbit first, so that, in particular, coladjmats[1] must be an 3819# identity matrix. 3820# 3821# The techniques used in this function are described in: 3822# Praeger and Soicher, "Low Rank Representations and Graphs for 3823# Sporadic Groups", CUP, Cambridge, 1997. 3824# 3825# May 2018: The efficiency of this function has been improved for 3826# the case when not all G-orbitals are self-paired. 3827# 3828local coladjmats,include,i,j,M,C,rank,comb,loc,sum,degree,prim,result; 3829if not IsList(gpin) and not IsPermGroup(gpin) then 3830 Error("usage: VertexTransitiveDRGs( <List> or <PermGroup> )"); 3831fi; 3832if IsPermGroup(gpin) then 3833 # Remark: OrbitalGraphColadjMats will check if gpin is transitive, 3834 # so we do not do this here. 3835 if IsTrivial(gpin) then 3836 Error("Input group must must be non-trivial,"); 3837 fi; 3838 coladjmats := OrbitalGraphColadjMats(gpin); 3839else 3840 coladjmats := gpin; 3841fi; 3842if Length(coladjmats)<2 or not IsMatrix(coladjmats[1]) 3843 or not IsInt(coladjmats[1][1][1]) then 3844 Error("<coladjmats> must be a list of integer matrices of length > 1,"); 3845fi; 3846rank:=Length(coladjmats); 3847prim:=true; 3848for i in [2..rank] do 3849 if LocalInfoMat(coladjmats[i],1).localDiameter=(-1) then 3850 # The i-th orbital graph is not (strongly) connected. 3851 prim:=false; 3852 break; 3853 fi; 3854od; 3855degree:=Sum(Sum(List(coladjmats,a->a[1]))); 3856result:=rec(degree:=degree, rank:=rank, isPrimitive:=prim, 3857 orbitalCombinations:=[], intersectionArrays:=[]); 3858include:=ListWithIdenticalEntries(rank,true); 3859include[1]:=false; # corresponding to the trivial orbital 3860M:=[]; 3861for i in [1..rank] do 3862 if coladjmats[i][1][i]=0 then 3863 Error("Error in <coladjmats>[",i,"]"); 3864 fi; 3865 if not include[i] then 3866 continue; 3867 fi; 3868 if coladjmats[i][i][1]=1 then 3869 # The orbital corresponding to coladjmats[i] is self-paired. 3870 Add(M,[i]); 3871 else 3872 j:=First([i+1..rank],x->include[x] and coladjmats[x][i][1]=1); 3873 # The orbital corresponding to coladjmats[i] is paired with 3874 # the orbital corresponding to coladjmats[j]. 3875 Add(M,[i,j]); 3876 include[j]:=false; 3877 fi; 3878od; 3879for comb in Combinations(M) do 3880 if comb<>[] then 3881 C:=Union(comb); 3882 sum:=Sum(coladjmats{C}); 3883 loc:=LocalInfoMat(sum,1); 3884 if loc.localDiameter <> -1 and not (-1 in Flat(loc.localParameters)) then 3885 # We've found a DRG. 3886 Add(result.orbitalCombinations,C); 3887 Add(result.intersectionArrays,loc.localParameters); 3888 fi; 3889 fi; 3890od; 3891return result; 3892end); 3893 3894DeclareOperation("MaximumClique",[IsRecord]); 3895InstallMethod(MaximumClique,"for GRAPE graph",[IsRecord],0, 3896function(gamma) 3897# 3898# Returns a clique C of maximum size of the simple graph gamma, 3899# and if gamma.maximumClique is unbound, sets gamma.maximumClique 3900# to be an immutable copy of C. 3901# 3902local G,delta,C,CC,lower,upper,mid; 3903if not IsGraph(gamma) then 3904 TryNextMethod(); 3905fi; 3906if not IsSimpleGraph(gamma) then 3907 Error("<gamma> must be a simple graph"); 3908fi; 3909if IsBound(gamma.maximumClique) then 3910 return ShallowCopy(gamma.maximumClique); 3911fi; 3912if gamma.order=0 then 3913 C:=[]; 3914 gamma.maximumClique:=Immutable(C); 3915 return C; 3916elif IsNullGraph(gamma) then 3917 C:=[1]; 3918 gamma.maximumClique:=Immutable(C); 3919 return C; 3920elif IsCompleteGraph(gamma) then 3921 C:=[1..gamma.order]; 3922 gamma.maximumClique:=Immutable(C); 3923 return C; 3924fi; 3925G:=AutomorphismGroup(gamma); 3926if G=gamma.group then 3927 delta:=gamma; 3928else 3929 delta:=NewGroupGraph(G,gamma); # to take full advantage of Aut(gamma) 3930fi; 3931lower:=1; 3932C:=[1]; 3933if IsBound(gamma.minimumVertexColouring) then 3934 upper:=Length(Set(gamma.minimumVertexColouring))+1; 3935else 3936 upper:=Maximum(VertexDegrees(delta))+2; 3937fi; 3938while upper-lower>1 do 3939 # Loop invariant: lower and upper are integers, 3940 # max clique size is in [lower,upper), 3941 # and C is a clique of size lower. 3942 mid:=Int((lower+upper)/2); 3943 CC:=CompleteSubgraphsOfGivenSize(delta,mid,0); 3944 if CC=[] then 3945 upper:=mid; 3946 else 3947 lower:=mid; 3948 C:=CC[1]; 3949 fi; 3950od; 3951gamma.maximumClique:=Immutable(C); 3952return C; 3953end); 3954 3955DeclareOperation("MaximumCompleteSubgraph",[IsRecord]); 3956InstallMethod(MaximumCompleteSubgraph,"for GRAPE graph",[IsRecord],0, 3957function(gamma) 3958# 3959# Implements alternative name for MaximumClique. 3960# 3961if not IsGraph(gamma) then 3962 TryNextMethod(); 3963fi; 3964return MaximumClique(gamma); 3965end); 3966 3967DeclareAttribute("CliqueNumber",IsRecord); 3968InstallMethod(CliqueNumber,"for GRAPE graph",[IsRecord],0, 3969function(gamma) 3970# 3971# Returns the size of a largest clique of the simple graph gamma. 3972# 3973if not IsGraph(gamma) then 3974 TryNextMethod(); 3975fi; 3976if not IsSimpleGraph(gamma) then 3977 Error("<gamma> must be a simple graph"); 3978fi; 3979return Length(MaximumClique(gamma)); 3980end); 3981 3982BindGlobal("GRAPE_CliqueCovering",function(arg) 3983# 3984# Let gamma:=arg[1] be a simple graph and let k:=arg[2] be 3985# a non-negative integer. 3986# 3987# This function returns a covering of gamma by at most k pairwise disjoint 3988# non-empty cliques if such a covering exists, and otherwise returns fail. 3989# 3990# A returned covering is given as a set of sets, forming a partition 3991# of the vertex set of gamma into at most k non-empty cliques. 3992# 3993# If arg[3] is bound then it must be a non-negative integer, such that 3994# no clique in any clique k-covering of gamma has size > arg[3]. 3995# 3996local gamma,k,m,cliquecovering,delta,cov,C,c,exhaustive_search,smallorder; 3997 3998cliquecovering := function(delta,k,start,olddelta) 3999# 4000# Let delta be a simple graph, and let k be a non-negative integer. 4001# 4002# Suppose that the boolean variable exhaustive_search 4003# (global to this function) has value true. Then this function 4004# returns a covering of delta by at most k pairwise disjoint 4005# non-empty cliques, if such a covering exists; otherwise fail is returned. 4006# 4007# Now suppose that exhaustive_search = false. Then start must 4008# be an integer, and this function tries to find a covering of 4009# delta by at most k pairwise disjoint non-empty cliques of size <= start, 4010# and returns such a covering if found. If no such covering is found 4011# (although one may still exist), then fail is returned. 4012# 4013# If start is an integer, then it must be non-negative, we let m=start, 4014# ignore olddelta, and it is assumed that there is *no* partition of 4015# the vertices of delta into <=k cliques such that some part has size > m. 4016# 4017# Now suppose start is not an integer. Then oldelta must be a simple graph, 4018# start must be a clique of olddelta, delta is the subgraph induced on the 4019# set of vertices of olddelta not in start, and we let m=Length(start). 4020# Furthermore, it is assumed that there is *no* partition of the vertices 4021# of olddelta into <=k+1 cliques such that some part has size > m or 4022# some part P has size m and P<start (in the usual lex-order on GAP sets). 4023# 4024# For this function (regardless of the value of exhaustive_search), 4025# a returned covering is given as a list of lists of vertex-names of delta. 4026# (These lists are not necessarily sorted, but contain no repeated elements.) 4027# In addition, we assume that on the initial call to this recursive function 4028# that m is an integer and delta.names=[1..delta.order]. 4029# 4030local m,C,CC,c,d,t,s,cov,newdelta,D,K,A,translation,wts,i,j; 4031if IsInt(start) then 4032 m:=start; 4033else 4034 m:=Length(start); 4035fi; 4036if delta.order=0 then 4037 return []; 4038elif m*k<delta.order then 4039 # in particular, if m=0 or k=0 4040 return fail; 4041elif k=1 then 4042 if IsCompleteGraph(delta) then 4043 return [ShallowCopy(delta.names)]; 4044 else 4045 return fail; 4046 fi; 4047fi; 4048if not IsInt(start) then 4049 translation:=Difference(Vertices(olddelta),start); 4050 # translation[i] is the vertex in olddelta corresponding to 4051 # the i-th vertex in delta. 4052fi; 4053s:=m; 4054while s*k>=delta.order do 4055 if exhaustive_search then 4056 if IsTrivial(delta.group) then 4057 C:=CompleteSubgraphsOfGivenSize(delta,s,1,true); 4058 else 4059 C:=CompleteSubgraphsOfGivenSize(delta,s,2,true); 4060 fi; 4061 if not IsInt(start) then 4062 CC:=[]; 4063 for c in C do 4064 t:=translation{c}; 4065 d:=Union(t,Filtered(start,x->IsSubset(Adjacency(olddelta,x),t))); 4066 if Length(d)>m then 4067 continue; 4068 fi; 4069 if Length(d)=m and 4070 (d<start or SmallestImageSet(olddelta.group,d)<start) then 4071 continue; 4072 fi; 4073 Add(CC,c); 4074 od; 4075 C:=CC; 4076 fi; 4077 if s*k=delta.order and Size(delta.group)<=smallorder then 4078 # Use exact cover. 4079 D:=Graph(delta.group,C,OnSets, 4080 function(x,y) return Intersection(x,y)=[]; end); 4081 wts:=List([1..D.order],x->ListWithIdenticalEntries(delta.order,0)); 4082 for i in [1..D.order] do 4083 for j in D.names[i] do 4084 wts[i][j]:=1; 4085 od; 4086 od; 4087 K:=CompleteSubgraphsOfGivenSize(D, 4088 ListWithIdenticalEntries(delta.order,1),0,true,true,wts); 4089 if K=[] then 4090 return fail; 4091 else 4092 return List(K[1],x->delta.names{D.names[x]}); 4093 fi; 4094 elif not IsTrivial(delta.group) then 4095 C:=Set(List(C,x->SmallestImageSet(delta.group,x))); 4096 fi; 4097 else 4098 C:=CompleteSubgraphsOfGivenSize(delta,s,0,true); 4099 fi; 4100 for c in C do 4101 A:=Difference(Vertices(delta),c); 4102 newdelta:=InducedSubgraph(delta,A,Stabilizer(delta.group,c,OnSets)); 4103 if exhaustive_search then 4104 cov:=cliquecovering(newdelta,k-1,c,delta); 4105 else 4106 cov:=cliquecovering(newdelta,k-1,s,0); 4107 fi; 4108 if cov<>fail then 4109 return Concatenation(cov,[delta.names{c}]); 4110 elif not exhaustive_search then 4111 # We give up. 4112 return fail; 4113 fi; 4114 od; 4115 s:=s-1; 4116od; 4117return fail; 4118end; 4119 4120if not (Length(arg) in [2,3]) then 4121 Error("must have 2 or 3 arguments"); 4122fi; 4123gamma:=arg[1]; 4124k:=arg[2]; 4125if not IsGraph(gamma) or not IsInt(k) then 4126 Error("usage: GRAPE_CliqueCovering( <Graph>, <Int> [, <Int> ] )"); 4127elif not IsSimpleGraph(gamma) then 4128 Error("<arg[1]> must be a simple graph"); 4129elif k<0 then 4130 Error("<arg[2]> must be non-negative"); 4131fi; 4132if IsBound(arg[3]) then 4133 m:=arg[3]; 4134 if not IsInt(m) then 4135 Error("usage: GRAPE_CliqueCovering( <Graph>, <Int> [, <Int> ] )"); 4136 elif m<0 then 4137 Error("<arg[3]> must be non-negative"); 4138 fi; 4139 if k*m<gamma.order then 4140 return fail; 4141 fi; 4142fi; 4143if gamma.order=0 then 4144 return []; 4145elif k=0 then 4146 return fail; 4147fi; 4148if IsCompleteGraph(gamma) then 4149 return [[1..gamma.order]]; 4150elif k=1 then 4151 return fail; 4152fi; 4153C:=Bicomponents(ComplementGraph(gamma)); 4154if C<>[] then 4155 # The complement of the non-complete graph gamma is bipartite. 4156 return Set(List(C,Set)); 4157elif k=2 then 4158 return fail; 4159fi; 4160delta:=NewGroupGraph(AutomorphismGroup(gamma),gamma); 4161delta.names:=Immutable([1..delta.order]); 4162if not IsBound(m) then 4163 m:=CliqueNumber(delta); 4164fi; 4165# 4166smallorder:=24; # To optimise when exact cover is used. 4167# smallorder can be given any positive integer value, but a value 4168# in the range 8 to 120 seems to work well. 4169# 4170exhaustive_search:=false; 4171cov:=cliquecovering(delta,k,m,0); 4172if cov=fail then 4173 exhaustive_search:=true; 4174 cov:=cliquecovering(delta,k,m,0); 4175 if cov=fail then 4176 return fail; 4177 fi; 4178fi; 4179for c in cov do 4180 Sort(c); 4181od; 4182Sort(cov); 4183# 4184# Check. 4185# 4186if not IsSet(cov) 4187 or Length(cov)>k 4188 or not ForAll(cov,x->x<>[] and IsSet(x)) 4189 or Union(cov)<>Vertices(gamma) 4190 or Sum(List(cov,Length))<>gamma.order 4191 or not ForAll(cov,x->IsCompleteGraph(InducedSubgraph(gamma,x))) then 4192 # This should never happen. 4193 Error("BUG: returned cov is not a clique k-covering given as a set of pairwise disjoint non-empty cliques"); 4194fi; 4195# 4196# End of check. 4197# 4198return cov; 4199end); 4200 4201BindGlobal("VertexColouring",function(arg) 4202# 4203# Let gamma:=arg[1] be a simple graph. Then this function returns 4204# a proper vertex-colouring of gamma. A proper vertex-colouring of 4205# gamma is given as a dense list C of length gamma.order, 4206# such that Set(C)=[1..Maximum(C)], where C[i] is the 4207# "colour" of the i-th vertex, and C[i]<>C[j] if [i,j] is an 4208# edge of gamma. 4209# 4210# If k:=arg[2] is bound, then it must be a non-negative integer, 4211# and a colouring using at most k colours is returned, or `fail' 4212# iff no such colouring exists. 4213# 4214# If arg[2] is unbound then a greedy algorithm only is used. 4215# 4216# If arg[3] is bound then it must be a non-negative integer, such that 4217# there is no monochromatic set of vertices of size > arg[3] in 4218# any vertex k-colouring of gamma. 4219# 4220local gamma,k,m,i,j,g,c,C,orb,a,adj,adjs,adjcolours,maxcolour,im,gens,cov; 4221if not (Length(arg) in [1,2,3]) then 4222 Error("must have 1, 2 or 3 arguments"); 4223fi; 4224gamma:=arg[1]; 4225if not IsGraph(gamma) then 4226 Error("usage: VertexColouring( <Graph> [, <Int> [, <Int> ]] )"); 4227elif not IsSimpleGraph(gamma) then 4228 Error("<arg[1]> not a simple graph"); 4229fi; 4230if IsBound(arg[2]) then 4231 k:=arg[2]; 4232 if not IsInt(k) then 4233 Error("usage: VertexColouring( <Graph> [, <Int> [, <Int> ]] )"); 4234 elif k<0 then 4235 Error("<arg[2]> must be non-negative"); 4236 fi; 4237else 4238 k:=gamma.order; 4239fi; 4240if IsBound(arg[3]) then 4241 m:=arg[3]; 4242 if not IsInt(m) then 4243 Error("usage: VertexColouring( <Graph> [, <Int> [, <Int> ]] )"); 4244 elif m<0 then 4245 Error("<arg[3]> must be non-negative"); 4246 fi; 4247 if k*m<gamma.order then 4248 return fail; 4249 fi; 4250fi; 4251if IsBound(gamma.minimumVertexColouring) then 4252 C:=gamma.minimumVertexColouring; 4253 if k<Length(Set(C)) then # k < chromatic number of gamma 4254 return fail; 4255 else 4256 return ShallowCopy(C); 4257 fi; 4258elif IsBound(gamma.maximumClique) then 4259 if k<Length(gamma.maximumClique) then 4260 return fail; 4261 fi; 4262fi; 4263if gamma.order=0 then 4264 return []; 4265fi; 4266# 4267# First try a greedy algorithm. 4268# 4269C:=ListWithIdenticalEntries(gamma.order,0); 4270maxcolour:=0; 4271gens:=GeneratorsOfGroup(gamma.group); 4272for i in [1..Length(gamma.representatives)] do 4273 orb:=[gamma.representatives[i]]; 4274 adjs:=[]; 4275 adj:=gamma.adjacencies[i]; 4276 adjs[orb[1]]:=adj; 4277 # colour vertex orb[1] 4278 adjcolours:=BlistList([1..maxcolour+1],[]); 4279 for a in adj do 4280 if C[a]>0 then 4281 adjcolours[C[a]]:=true; 4282 fi; 4283 od; 4284 c:=1; 4285 while adjcolours[c] do 4286 c:=c+1; 4287 od; 4288 if c>maxcolour then 4289 maxcolour:=c; 4290 if maxcolour>k then 4291 break; 4292 fi; 4293 fi; 4294 C[orb[1]]:=c; 4295 for j in orb do 4296 for g in gens do 4297 im:=j^g; 4298 if C[im]=0 then 4299 Add(orb,im); 4300 adj:=OnTuples(adjs[j],g); 4301 adjs[im]:=adj; 4302 # colour vertex im 4303 adjcolours:=BlistList([1..maxcolour+1],[]); 4304 for a in adj do 4305 if C[a]>0 then 4306 adjcolours[C[a]]:=true; 4307 fi; 4308 od; 4309 c:=1; 4310 while adjcolours[c] do 4311 c:=c+1; 4312 od; 4313 if c>maxcolour then 4314 maxcolour:=c; 4315 if maxcolour>k then 4316 break; 4317 fi; 4318 fi; 4319 C[im]:=c; 4320 fi; 4321 od; 4322 Unbind(adjs[j]); 4323 if maxcolour>k then 4324 break; 4325 fi; 4326 od; 4327 if maxcolour>k then 4328 break; 4329 fi; 4330od; 4331if maxcolour<=k then 4332 # C is a k-colouring. 4333 if IsBound(gamma.maximumClique) and Length(gamma.maximumClique)=Length(Set(C)) then 4334 # C is a minimum vertex-colouring of gamma. 4335 gamma.minimumVertexColouring:=Immutable(C); 4336 fi; 4337 return C; 4338fi; 4339# Otherwise, we need to work harder. 4340if IsBound(arg[3]) then 4341 cov:=GRAPE_CliqueCovering(ComplementGraph(gamma),k,arg[3]); 4342else 4343 cov:=GRAPE_CliqueCovering(ComplementGraph(gamma),k); 4344fi; 4345if cov=fail then 4346 return fail; 4347fi; 4348# Otherwise, we make C into a k-colouring from cov. 4349C:=[]; 4350for i in [1..Length(cov)] do 4351 for j in cov[i] do 4352 C[j]:=i; 4353 od; 4354od; 4355if IsBound(gamma.maximumClique) and Length(gamma.maximumClique)=Length(Set(C)) then 4356 # C is a minimum vertex-colouring of gamma. 4357 gamma.minimumVertexColouring:=Immutable(C); 4358fi; 4359return C; 4360end); 4361 4362BindGlobal("GRAPE_MinimumCliqueCovering",function(gamma) 4363# 4364# Let gamma be a simple graph. Then this function returns a clique covering 4365# of gamma of minimum size. The returned covering is given as a set of sets, 4366# forming a partition of the vertex set of gamma into non-empty cliques. 4367# 4368local C,CC,lwr,lower,upper,mid,i,delta; 4369if not IsGraph(gamma) then 4370 Error("usage: GRAPE_MinimumCliqueCovering( <Graph> )"); 4371elif not IsSimpleGraph(gamma) then 4372 Error("<gamma> must be a simple graph"); 4373fi; 4374if gamma.order=0 then 4375 return []; 4376elif IsCompleteGraph(gamma) then 4377 return [[1..gamma.order]]; 4378fi; 4379delta:=ComplementGraph(gamma); 4380C:=GRAPE_NumbersToSets(VertexColouring(delta)); 4381Sort(C); 4382# Now C is a partition of the vertex set of gamma into 4383# Length(C) cliques. 4384if Length(C)=2 then 4385 return C; 4386fi; 4387upper:=Length(C); 4388if Maximum(VertexDegrees(gamma))<(gamma.order-1)/2 then 4389 lwr:=gamma.order/CliqueNumber(gamma); 4390 if IsInt(lwr) then 4391 lower:=lwr-1; 4392 else 4393 lower:=Int(lwr); 4394 fi; 4395else 4396 lower:=CliqueNumber(delta)-1; 4397fi; 4398while upper-lower>1 do 4399 # Loop invariant: lower and upper are integers, 4400 # The clique covering number of gamma is in (lower,upper] 4401 # and C is a clique covering of size upper. 4402 mid:=Int((lower+upper)/2); 4403 CC:=GRAPE_CliqueCovering(gamma,mid); 4404 if CC=fail then 4405 lower:=mid; 4406 else 4407 upper:=Length(CC); # which is <= mid and > lower. 4408 C:=CC; 4409 fi; 4410od; 4411return C; 4412end); 4413 4414DeclareOperation("MinimumVertexColouring",[IsRecord]); 4415InstallMethod(MinimumVertexColouring,"for GRAPE graph",[IsRecord],0, 4416function(gamma) 4417# 4418# Let gamma be a simple graph. Then this function returns a proper vertex- 4419# colouring C of gamma, using as few colours as possible, and if 4420# gamma.minimumVertexColouring is unbound, sets gamma.minimumVertexColouring 4421# to be an immutable copy of C. 4422# 4423# A proper vertex-colouring of gamma is given as a list C of 4424# length gamma.order of positive integers, such that C[i] is the 4425# "colour" of the i-th vertex, and C[i]<>C[j] if [i,j] is an edge 4426# of gamma. 4427# 4428local cov,C,i,j; 4429if not IsGraph(gamma) then 4430 TryNextMethod(); 4431fi; 4432if not IsSimpleGraph(gamma) then 4433 Error("<gamma> must be a simple graph"); 4434fi; 4435if IsBound(gamma.minimumVertexColouring) then 4436 return ShallowCopy(gamma.minimumVertexColouring); 4437fi; 4438cov:=GRAPE_MinimumCliqueCovering(ComplementGraph(gamma)); 4439C:=[]; 4440for i in [1..Length(cov)] do 4441 for j in cov[i] do 4442 C[j]:=i; 4443 od; 4444od; 4445gamma.minimumVertexColouring:=Immutable(C); 4446return C; 4447end); 4448 4449DeclareAttribute("ChromaticNumber",IsRecord); 4450InstallMethod(ChromaticNumber,"for GRAPE graph",[IsRecord],0, 4451function(gamma) 4452# 4453# Let gamma be a simple graph. Then this function returns the 4454# chromatic number of gamma, that is, the minimum number of 4455# colours needed to properly vertex-colour gamma. 4456# 4457if not IsGraph(gamma) then 4458 TryNextMethod(); 4459fi; 4460if not IsSimpleGraph(gamma) then 4461 Error("<gamma> must be a simple graph"); 4462fi; 4463return Length(Set(MinimumVertexColouring(gamma))); 4464end); 4465 4466BindGlobal("IsGraphWithColourClasses",function(obj) 4467return IsRecord(obj) and IsBound(obj.graph) and IsGraph(obj.graph) and IsBound(obj.colourClasses); 4468end); 4469 4470BindGlobal("MonochromaticColourClasses",function(gamma) 4471# Returns colour-classes list with all vertices having the same 4472# colour for the vertices of the graph gamma. 4473if not IsGraph(gamma) then 4474 Error("usage: MonochromaticColourClasses( <Graph> )"); 4475fi; 4476if gamma.order=0 then 4477 return []; 4478else 4479 return [[1..gamma.order]]; 4480fi; 4481end); 4482 4483BindGlobal("CheckColourClasses",function(gamma,col) 4484# 4485# Checks whether col is a valid list of colour-classes for the 4486# graph gamma. 4487# 4488if not (IsGraph(gamma) and IsList(col)) then 4489 Error("usage: CheckColourClasses( <Graph>, <List> )"); 4490fi; 4491if not ForAll(col,x->IsSet(x) and x<>[]) then 4492 Error("each colour-class must be a non-empty set"); 4493fi; 4494if Union(col)<>[1..gamma.order] then 4495 Error("the union of the colour-classes is not equal to the vertex set of <gamma>"); 4496fi; 4497if Sum(List(col,Length))>gamma.order then 4498 Error("the colour-classes must be pairwise disjoint"); 4499fi; 4500return; 4501end); 4502 4503# Set up temporary directory for use with nauty/dreadnaut or bliss. 4504BindGlobal("GRAPE_nautytmpdir",DirectoryTemporary()); 4505Add(GAPInfo.PostRestoreFuncs,function() 4506 MakeReadWriteGlobal("GRAPE_nautytmpdir"); 4507 Unbind(GRAPE_nautytmpdir); 4508 BindGlobal("GRAPE_nautytmpdir",DirectoryTemporary()); 4509end); 4510 4511BindGlobal("PrintStreamNautyGraph",function(stream,gamma,col) 4512 local adj, i, j; 4513 PrintTo(stream,"d\n$1n",gamma.order,"g\n"); 4514 for i in [1..gamma.order] do 4515 adj:=Adjacency(gamma,i); 4516 if adj=[] then 4517 if i<gamma.order then 4518 AppendTo(stream,";\n"); 4519 else 4520 AppendTo(stream,".\n"); 4521 fi; 4522 else 4523 for j in [1..Length(adj)] do 4524 if j<Length(adj) then 4525 AppendTo(stream,adj[j],"\n"); 4526 elif i<gamma.order then 4527 AppendTo(stream,adj[j],";\n"); 4528 else 4529 AppendTo(stream,adj[j],".\n"); 4530 fi; 4531 od; 4532 fi; 4533 od; 4534 AppendTo(stream,"f[\n"); 4535 if col<>MonochromaticColourClasses(gamma) then 4536 for i in [1..Length(col)] do 4537 for j in [1..Length(col[i])] do 4538 AppendTo(stream,col[i][j]); 4539 if j<Length(col[i]) then 4540 AppendTo(stream,","); 4541 elif i<Length(col) then 4542 AppendTo(stream,"|"); 4543 fi; 4544 AppendTo(stream,"\n"); 4545 od; 4546 od; 4547 fi; 4548 AppendTo(stream,"]\n"); 4549end); 4550 4551BindGlobal("ReadOutputNauty",function(file) 4552# by Alexander Hulpke 4553 local f, bas, sgens, l, s, p, i, deg, processperm, pi; 4554 4555 processperm:=function() 4556 if Length(pi)=0 then return; fi; 4557 if deg=fail then 4558 deg:=Length(pi); 4559 else 4560 if Length(pi)<>deg then 4561 Info(InfoWarning,1,"degree discrepancy in nauty output!", 4562 Length(pi),"vs",deg); 4563 fi; 4564 fi; 4565 Add(sgens,PermList(pi)); 4566 pi:=[]; 4567 end; 4568 4569 deg:=fail; 4570 f:=InputTextFile(file); 4571 if f=fail then 4572 Error("cannot find output produced by `dreadnaut'"); 4573 fi; 4574 bas:=[]; 4575 sgens:=[]; 4576 pi:=[]; 4577 while not IsEndOfStream(f) do 4578 l:=ReadLine(f); 4579 if l<>fail then 4580 l:=Chomp(l); 4581# Print(l,"\n"); 4582 if Length(l)>4 and l{[1..5]}="level" then 4583 processperm(); 4584 s:=SplitString(l,";"); 4585 s:=s[Length(s)-1]; # should be " x...x fixed" 4586 if Length(s)<4 or s{[Length(s)-4..Length(s)]}<>"fixed" then 4587 Error("unparsable line ",l); 4588 fi; 4589 s:=s{[1..Length(s)-6]}; 4590 while s[1]=' ' do 4591 s:=s{[2..Length(s)]}; 4592 od; 4593 Add(bas,Int(s)); 4594 elif ForAll(l,x->x in CHARS_DIGITS or x=' ') then 4595 if Length(pi)>0 and (Length(l)<5 or l{[1..4]}<>" ") then 4596 processperm(); # permutation starts -- clean out old 4597 fi; 4598 s:=SplitString(l,[]," "); 4599 Append(pi,List(s,Int)); 4600 elif deg<>fail then 4601 processperm(); 4602 fi; 4603 4604 fi; 4605 od; 4606 bas:=Reversed(bas); 4607 CloseStream(f); 4608 return [sgens,bas]; 4609end); 4610 4611BindGlobal("ReadCanonNauty",function(file) 4612# by Alexander Hulpke 4613 local f, can, l, deg, s, i; 4614 f:=InputTextFile(file); 4615 if f=fail then 4616 Error("cannot find canonization produced by `dreadnaut'"); 4617 fi; 4618 can:=[]; 4619 # first line: degree 4620 l:=ReadLine(f);l:=Chomp(l); 4621# Print(l,"\n"); 4622 deg:=Int(l); 4623 # now read in until you have enough integers for the permutation -- the 4624 # rest is the relabelled graph and can be discarded 4625 while Length(can)<deg do 4626 l:=ReadLine(f);l:=Chomp(l); 4627# Print(l,"\n"); 4628 s:=SplitString(l,' '); 4629 for i in s do 4630 if Length(i)>0 and Length(can)<deg then 4631 Add(can,Int(i)); 4632 fi; 4633 od; 4634 od; 4635 CloseStream(f); 4636 return PermList(can); 4637end); 4638 4639BindGlobal("SetAutGroupCanonicalLabellingNauty",function(gr,setcanon) 4640# 4641# Sets the autGroup component (if not already bound) and the 4642# canonicalLabelling component (if not already bound and setcanon=true) 4643# of the graph or graph with colour-classes gr. 4644# Uses the nauty system. 4645# 4646 local gamma,col,ftmp1,ftmp2,fdre,fg,ftmp1_stream,ftmp2_stream,fdre_stream,gp; 4647 if IsBound(gr.canonicalLabelling) then 4648 setcanon:=false; 4649 fi; 4650 if IsBound(gr.autGroup) and not setcanon then 4651 return; 4652 fi; 4653 if IsGraph(gr) then 4654 gamma:=gr; 4655 col:=MonochromaticColourClasses(gamma); 4656 else 4657 gamma:=gr.graph; 4658 col:=gr.colourClasses; 4659 CheckColourClasses(gamma,col); 4660 fi; 4661 if gamma.order<=1 then 4662 if not IsBound(gr.autGroup) then 4663 gr.autGroup:=Group([],()); 4664 fi; 4665 if setcanon then 4666 gr.canonicalLabelling:=(); 4667 fi; 4668 return; 4669 fi; 4670 4671 ftmp1:=Filename(GRAPE_nautytmpdir,"ftmp1"); 4672 ftmp2:=Filename(GRAPE_nautytmpdir,"ftmp2"); 4673 4674 # In principle redundant, but a failed call might have left files sitting 4675 # -- just throw out what will be overwritten anyhow. 4676 RemoveFile(ftmp1); 4677 RemoveFile(ftmp2); 4678 4679 fdre:=""; 4680 fdre_stream:=OutputTextString(fdre,false); 4681 SetPrintFormattingStatus(fdre_stream,false); 4682 PrintStreamNautyGraph(fdre_stream,gamma,col); 4683 4684 if not setcanon then 4685 if IsSimpleGraph(gamma) then 4686 AppendTo( fdre_stream, "> ", ftmp1, " p,xq\n" ); 4687 else 4688 AppendTo( fdre_stream, "> ", ftmp1, " p,*=13,k=1 10,xq\n" ); 4689 fi; 4690 else 4691 if IsSimpleGraph(gamma) then 4692 AppendTo( fdre_stream, "> ", ftmp1, " p,cx\n>> ", ftmp2, " bq\n" ); 4693 else 4694 AppendTo( fdre_stream, "> ", ftmp1, " p,*=13,k=1 10,cx\n>> ", ftmp2, 4695 " bq\n" ); 4696 fi; 4697 fi; 4698 4699 CloseStream(fdre_stream); 4700 4701 # initialize tmp2 file with degree 4702 ftmp2_stream:=OutputTextFile(ftmp2,false); 4703 SetPrintFormattingStatus(ftmp2_stream,false); 4704 PrintTo(ftmp2_stream,gamma.order,"\n"); 4705 CloseStream(ftmp2_stream); 4706 4707 GRAPE_Exec(GRAPE_DREADNAUT_EXE, [], InputTextString(fdre), OutputTextUser()); 4708 4709 if not IsBound(gr.autGroup) then 4710 fg:=ReadOutputNauty(ftmp1); 4711 # fg[1]=stronggens, fg[2]=base 4712 gp:=GroupWithGenerators(fg[1],()); 4713 SetStabChainMutable(gp,StabChainBaseStrongGenerators(fg[2],fg[1],())); 4714 gr.autGroup:=gp; 4715 fi; 4716 if setcanon then 4717 gr.canonicalLabelling:=ReadCanonNauty(ftmp2); 4718 fi; 4719 4720 RemoveFile(ftmp1); 4721 RemoveFile(ftmp2); 4722 4723end); 4724 4725BindGlobal("PrintStreamBlissGraph",function(stream,gamma,col) 4726# Original procedure by Jerry James. Updated by Leonard Soicher. 4727 local adj, i, j, nedges; 4728 if IsRegularGraph(gamma) and gamma.order > 0 then 4729 nedges := gamma.order*Length(Adjacency(gamma,1)); 4730 else 4731 nedges:=0; 4732 for i in [1..gamma.order] do 4733 nedges := nedges + Length(Adjacency(gamma,i)); 4734 od; 4735 fi; 4736 # nedges = no. of directed edges of gamma. 4737 PrintTo(stream,"p edge ",gamma.order," ",nedges,"\n"); 4738 if col<>MonochromaticColourClasses(gamma) then 4739 for i in [1..Length(col)] do 4740 for j in [1..Length(col[i])] do 4741 AppendTo(stream, "n ", col[i][j], " ", i, "\n"); 4742 od; 4743 od; 4744 fi; 4745 for i in [1..gamma.order] do 4746 adj:=Adjacency(gamma,i); 4747 if adj<>[] then 4748 for j in [1..Length(adj)] do 4749 AppendTo(stream, "e ", i, " ", adj[j], "\n"); 4750 od; 4751 fi; 4752 od; 4753end); 4754 4755BindGlobal("ReadOutputBliss",function(f,canon) 4756# Original procedure by Jerry James. Updated by Leonard Soicher. 4757 local gens, l, i, pi, can; 4758 4759 if f=fail then 4760 Error("cannot find output produced by `bliss'"); 4761 fi; 4762 gens:=[]; 4763 can:=[]; 4764 while not IsEndOfStream(f) do 4765 l:=ReadLine(f); 4766 if l<>fail then 4767 l:=Chomp(l); 4768 if Length(l)>11 and l{[1..11]}="Generator: " then 4769 pi:=EvalString(l{[12..Length(l)]}); 4770 Add(gens,pi); 4771 elif canon and Length(l)>20 and l{[1..20]}="Canonical labeling: " then 4772 can:=InverseSameMutability(EvalString(l{[21..Length(l)]})); 4773 fi; 4774 fi; 4775 od; 4776 CloseStream(f); 4777 return [gens,can]; 4778end); 4779 4780BindGlobal("SetAutGroupCanonicalLabellingBliss",function(gr, setcanon) 4781# 4782# Sets the autGroup component (if not already bound) and the 4783# canonicalLabelling component (if not already bound and setcanon=true) 4784# of the graph or graph with colour-classes gr. 4785# Uses the bliss system. 4786# 4787 local gamma,col,ftmp,ftmp_stream,fdre,fg,fdre_stream,gp; 4788 if IsBound(gr.canonicalLabelling) then 4789 setcanon:=false; 4790 fi; 4791 if IsBound(gr.autGroup) and not setcanon then 4792 return; 4793 fi; 4794 if IsGraph(gr) then 4795 gamma:=gr; 4796 col:=MonochromaticColourClasses(gamma); 4797 else 4798 gamma:=gr.graph; 4799 col:=gr.colourClasses; 4800 CheckColourClasses(gamma,col); 4801 fi; 4802 if gamma.order<=1 then 4803 if not IsBound(gr.autGroup) then 4804 gr.autGroup:=Group([],()); 4805 fi; 4806 if setcanon then 4807 gr.canonicalLabelling:=(); 4808 fi; 4809 return; 4810 fi; 4811 4812 fdre:=Filename(GRAPE_nautytmpdir,"fdre"); 4813 4814 # In principle redundant, but a failed call might have left files sitting 4815 # -- just throw out what will be overwritten anyhow. 4816 RemoveFile(fdre); 4817 4818 fdre_stream:=OutputTextFile(fdre,false); 4819 SetPrintFormattingStatus(fdre_stream,false); 4820 PrintStreamBlissGraph(fdre_stream,gamma,col); 4821 CloseStream(fdre_stream); 4822 4823 ftmp:=""; 4824 ftmp_stream:=OutputTextString(ftmp,false); 4825 SetPrintFormattingStatus(ftmp_stream,false); 4826 4827 if setcanon then 4828 GRAPE_Exec(GRAPE_BLISS_EXE, [ "-directed", "-can", fdre ], InputTextUser(), ftmp_stream); 4829 else 4830 GRAPE_Exec(GRAPE_BLISS_EXE, [ "-directed", fdre ], InputTextUser(), ftmp_stream); 4831 fi; 4832 4833 fg:=ReadOutputBliss(InputTextString(ftmp),setcanon); 4834 # fg[1]=gens for the aut group, 4835 # fg[2]=canonical labelling if setcanon=true, else the empty list 4836 if not IsBound(gr.autGroup) then 4837 gp:=GroupWithGenerators(fg[1],()); 4838 gr.autGroup:=gp; 4839 fi; 4840 if setcanon then 4841 gr.canonicalLabelling:=fg[2]; 4842 fi; 4843 4844 RemoveFile(fdre); 4845 4846end); 4847 4848BindGlobal("SetAutGroupCanonicalLabelling",function(arg) 4849# 4850# Let gr:=arg[1] and setcanon:=arg[2] (default: true). 4851# Sets the autGroup component (if not already bound) and the 4852# canonicalLabelling component (if not already bound and setcanon=true) 4853# of the graph or graph with colour-classes gr. 4854# 4855 local gr,setcanon; 4856 gr:=arg[1]; 4857 if IsBound(arg[2]) then 4858 setcanon:=arg[2]; 4859 else 4860 setcanon:=true; 4861 fi; 4862 if not (IsGraph(gr) or IsGraphWithColourClasses(gr)) or not IsBool(setcanon) then 4863 Error("usage: SetAutGroupCanonicalLabelling( <Graph> or <GraphWithColourClasses> [, <Bool> ] )"); 4864 fi; 4865 if GRAPE_NAUTY then 4866 SetAutGroupCanonicalLabellingNauty(gr, setcanon); 4867 else 4868 SetAutGroupCanonicalLabellingBliss(gr, setcanon); 4869 fi; 4870end); 4871 4872BindGlobal("AutGroupGraph",function(arg) 4873# 4874# Let gr:=arg[1] be a graph or a graph with colour-classes. 4875# 4876# If arg[2] is unbound (the ususal case) then this function returns 4877# the automorphism group of gr (making use of B.McKay's 4878# dreadnaut, nauty programs). 4879# 4880# If arg[2] is bound then gr must be a graph and arg[2] is 4881# a vertex-colouring (not necessarily proper) for gr 4882# (i.e. a list of colour-classes for the vertices of gr), 4883# in which case the subgroup of Aut(gr) preserving this colouring 4884# is returned instead of the full automorphism group. 4885# (Here a vertex-colouring is a list of sets, forming an ordered 4886# partition of the vertices. The set for the last colour may be omitted.) 4887# 4888local gr,gamma,col; 4889if IsBound(arg[2]) then 4890 if not IsGraph(arg[1]) or not IsList(arg[2]) then 4891 Error("usage: AutGroupGraph( <Graph> [, <List> ] ) or AutGroupGraph( <GraphWithColourClasses> )"); 4892 fi; 4893 gamma:=arg[1]; 4894 col:=arg[2]; 4895 if Union(col)<>[1..gamma.order] then 4896 # for backward compatibility 4897 Add(col,Difference([1..gamma.order],Union(col))); 4898 fi; 4899 CheckColourClasses(gamma,col); 4900 if col<>MonochromaticColourClasses(gamma) then 4901 gr:=rec(graph:=gamma,colourClasses:=col); 4902 else 4903 gr:=gamma; 4904 fi; 4905else 4906 gr:=arg[1]; 4907fi; 4908# Now deal with <gr>. 4909if not (IsGraph(gr) or IsGraphWithColourClasses(gr)) then 4910 Error("usage: AutGroupGraph( <Graph> [, <List> ] ) or AutGroupGraph( <GraphWithColourClasses> )"); 4911fi; 4912SetAutGroupCanonicalLabelling(gr,false); 4913return gr.autGroup; 4914end); 4915 4916InstallOtherMethod(AutomorphismGroup,"for graph or graph with colour-classes", 4917 [IsRecord],100, 4918function(gamma) 4919if not IsGraph(gamma) and not IsGraphWithColourClasses(gamma) then 4920 TryNextMethod(); 4921fi; 4922return AutGroupGraph(gamma); 4923end); 4924 4925BindGlobal("IsGraphIsomorphism",function(gr1,gr2,perm) 4926# 4927# Let gr1 and gr2 both be graphs or both be graphs with colour-classes. 4928# Then this function returns true if perm is an 4929# isomorphism from gr1 to gr2 (and false if not). 4930# 4931local gamma1,gamma2,col1,col2,u,g,i,j,x,aut1,aut2,adj1,adj2,reps1; 4932if not ((IsGraph(gr1) and IsGraph(gr2)) or (IsGraphWithColourClasses(gr1) and IsGraphWithColourClasses(gr2))) or not IsPerm(perm) then 4933 Error("usage: IsGraphIsomorphism( <Graph>, <Graph>, <Perm> ) or IsGraphIsomorphism( <GraphWithColourClasses>, <GraphWithColourClasses>, <Perm> )"); 4934fi; 4935if IsGraphWithColourClasses(gr1) then 4936 # both gr1 and gr2 are graphs with colour-classes 4937 gamma1:=gr1.graph; 4938 col1:=gr1.colourClasses; 4939 CheckColourClasses(gamma1,col1); 4940 gamma2:=gr2.graph; 4941 col2:=gr2.colourClasses; 4942 CheckColourClasses(gamma2,col2); 4943else 4944 # both gr1 and gr2 are graphs 4945 gamma1:=gr1; 4946 col1:=MonochromaticColourClasses(gamma1); 4947 gamma2:=gr2; 4948 col2:=MonochromaticColourClasses(gamma2); 4949fi; 4950if LargestMovedPoint(perm)>gamma1.order then 4951 return false; 4952fi; 4953if gamma1.order<>gamma2.order or 4954 VertexDegrees(gamma1) <> VertexDegrees(gamma2) then 4955 # the graphs are not isomorphic 4956 return false; 4957elif gamma1.order<=1 then 4958 return true; 4959fi; 4960if List(col1,c->OnSets(c,perm))<>col2 then 4961 return false; 4962fi; 4963# So now we know that perm takes col1 to col2. 4964if IsBound(gamma1.autGroup) and IsBound(gamma2.autGroup) then 4965 aut1:=gamma1.autGroup; 4966 aut2:=gamma2.autGroup; 4967 if aut1^perm<>aut2 then 4968 return false; 4969 fi; 4970else 4971 aut1:=Group(()); 4972 aut2:=Group(()); 4973fi; 4974# So now, either aut1 and aut2 are both trivial, or they 4975# are the full aut groups of gamma1 and gamma2, respectively, 4976# and aut1^perm=aut2. 4977reps1:=GRAPE_OrbitNumbers(aut1,gamma1.order).representatives; 4978for i in reps1 do 4979 adj1:=Adjacency(gamma1,i); 4980 adj2:=Adjacency(gamma2,i^perm); 4981 if OnSets(adj1,perm)<>adj2 then 4982 return false; 4983 fi; 4984od; 4985return true; 4986end); 4987 4988BindGlobal("GraphIsomorphism",function(arg) 4989# 4990# Let gr1:=arg[1] and gr2:=arg[2] both be graphs or both be 4991# graphs with colour-classes. 4992# Then this function returns an isomorphism from gr1 to gr2, if 4993# gr1 and gr2 are isomorphic, else returns fail. 4994# 4995# The optional boolean parameter firstunbindcanon=arg[3] determines 4996# whether or not the canonicalLabelling components of both gr1 and 4997# gr2 are first made unbound before proceeding. If 4998# firstunbindcanon=true (the default, safe and possibly slower option) 4999# then these components are first unbound. 5000# If firstunbindcanon=false, then an old canonical labelling 5001# is used when it exists. However, canonical labellings can depend on 5002# the version of nauty, the version of GRAPE, certain settings 5003# of nauty, and the compiler and computer used. 5004# Thus, if firstunbindcanon=false, the user must be 5005# sure that any canonicalLabelling component(s) which may already 5006# exist for gr1 or gr2 were created in exactly the same 5007# environment in which the user is presently computing. 5008# 5009local gr1,gr2,gamma1,gamma2,col1,col2,firstunbindcanon,g,i,j,x; 5010gr1:=arg[1]; 5011gr2:=arg[2]; 5012if IsBound(arg[3]) then 5013 firstunbindcanon:=arg[3]; 5014else 5015 firstunbindcanon:=true; 5016fi; 5017if not ((IsGraph(gr1) and IsGraph(gr2)) or (IsGraphWithColourClasses(gr1) and IsGraphWithColourClasses(gr2))) or not IsBool(firstunbindcanon) then 5018 Error("usage: GraphIsomorphism( <Graph>, <Graph> [, <Bool>] ) or GraphIsomorphism( <GraphWithColourClasses>, <GraphWithColourClasses> [, <Bool>] )"); 5019fi; 5020if IsGraphWithColourClasses(gr1) then 5021 # both gr1 and gr2 are graphs with colour-classes 5022 gamma1:=gr1.graph; 5023 col1:=gr1.colourClasses; 5024 CheckColourClasses(gamma1,col1); 5025 gamma2:=gr2.graph; 5026 col2:=gr2.colourClasses; 5027 CheckColourClasses(gamma2,col2); 5028else 5029 # both gr1 and gr2 are graphs 5030 gamma1:=gr1; 5031 col1:=MonochromaticColourClasses(gamma1); 5032 gamma2:=gr2; 5033 col2:=MonochromaticColourClasses(gamma2); 5034fi; 5035if firstunbindcanon then 5036 Unbind(gr1.canonicalLabelling); 5037 Unbind(gr2.canonicalLabelling); 5038fi; 5039if gamma1.order<>gamma2.order or 5040 VertexDegrees(gamma1) <> VertexDegrees(gamma2) then 5041 # the graphs are not isomorphic 5042 return fail; 5043elif List(col1,Length)<>List(col2,Length) then 5044 # incompatible colourings 5045 return fail; 5046elif gamma1.order<=1 then 5047 return (); 5048fi; 5049SetAutGroupCanonicalLabelling(gr1,true); 5050SetAutGroupCanonicalLabelling(gr2,true); 5051x:=LeftQuotient(gr1.canonicalLabelling,gr2.canonicalLabelling); 5052if IsGraphIsomorphism(gr1,gr2,x) then 5053 return x; 5054else 5055 return fail; 5056fi; 5057end); 5058 5059BindGlobal("IsIsomorphicGraph",function(arg) 5060# 5061# Let gr1:=arg[1] and gr2:=arg[2] both be graphs or both be 5062# graphs with colour-classes. 5063# Then this function returns true if gr1 and gr2 are isomorphic, 5064# else returns false. 5065# 5066# The optional boolean parameter firstunbindcanon=arg[3] determines 5067# whether or not the canonicalLabelling components of both gr1 and 5068# gr2 are first made unbound before proceeding. If 5069# firstunbindcanon=true (the default, safe and possibly slower option) 5070# then these components are first unbound. 5071# If firstunbindcanon=false, then an old canonical labelling 5072# is used when it exists. However, canonical labellings can depend on 5073# the version of nauty, the version of GRAPE, certain settings 5074# of nauty, and the compiler and computer used. 5075# Thus, if firstunbindcanon=false, the user must be 5076# sure that any canonicalLabelling component(s) which may already 5077# exist for gr1 or gr2 were created in exactly the same 5078# environment in which the user is presently computing. 5079# 5080if Length(arg)=2 then 5081 return IsPerm(GraphIsomorphism(arg[1],arg[2])); 5082elif Length(arg)=3 then 5083 return IsPerm(GraphIsomorphism(arg[1],arg[2],arg[3])); 5084else 5085 Error("number of arguments must be 2 or 3"); 5086fi; 5087end); 5088 5089BindGlobal("GraphIsomorphismClassRepresentatives",function(arg) 5090# 5091# Given a list L:=arg[1] of graphs, or of graphs with colour-classes, 5092# this function returns a list 5093# containing pairwise non-isomorphic elements of L, representing 5094# all the isomorphism classes of elements of L. 5095# 5096# The optional boolean parameter firstunbindcanon=arg[2] determines 5097# whether or not the canonicalLabelling components of all 5098# the graphs in L are first made unbound before proceeding. 5099# If firstunbindcanon=true (the default, safe and possibly slower option) 5100# then these components are first unbound. 5101# If firstunbindcanon=false, then an old canonical labelling 5102# is used when it exists. However, canonical labellings can depend on 5103# the version of nauty, the version of GRAPE, certain settings 5104# of nauty, and the compiler and computer used. 5105# Thus, if firstunbindcanon=false, the user must be 5106# sure that any canonicalLabelling component(s) which may already 5107# exist for graphs in L were created in exactly the same 5108# environment in which the user is presently computing. 5109# 5110local L,firstunbindcanon,reps,i,x,found; 5111L:=arg[1]; 5112if IsBound(arg[2]) then 5113 firstunbindcanon:=arg[2]; 5114else 5115 firstunbindcanon:=true; 5116fi; 5117if not (IsList(L) and IsBool(firstunbindcanon)) then 5118 Error("usage: GraphIsomorphismClassRepresentatives( <List> [, <Bool> ] )"); 5119fi; 5120if not (ForAll(L,IsGraph) or ForAll(L,IsGraphWithColourClasses)) then 5121 Error("<L> must be a list of graphs or a list of graphs with colour-classes"); 5122fi; 5123if firstunbindcanon then 5124 for x in L do 5125 Unbind(x.canonicalLabelling); 5126 od; 5127fi; 5128if Length(L)<=1 then 5129 return ShallowCopy(L); 5130fi; 5131reps:=[L[1]]; 5132for i in [2..Length(L)] do 5133 found:=false; 5134 for x in reps do 5135 if IsIsomorphicGraph(x,L[i],false) then 5136 found:=true; 5137 break; 5138 fi; 5139 od; 5140 if not found then 5141 Add(reps,L[i]); 5142 fi; 5143od; 5144return reps; 5145end); 5146 5147BindGlobal("PartialLinearSpaces",function(arg) 5148# 5149# Let s and t be positive integers. Then a *partial linear space* 5150# (P,L), with *parameters* s,t, consists of a set P of *points*, 5151# together with a set L of (s+1)-subsets of P called *lines*, 5152# such that every point is in exactly t+1 lines, and 5153# every pair of (distinct) points is contained in at most one line. 5154# The *point graph* of a partial linear space S having point-set 5155# P is the graph with vertex-set P and having (p,q) an edge iff 5156# p<>q and p,q lie on a common line of S. Two partial linear 5157# spaces (P,L) and (P',L') (with parameters s,t) are said 5158# to be *isomorphic* if there is a bijection P-->P' which induces 5159# a bijection L-->L'. 5160# 5161# This function returns a list of representatives of distinct isomorphism 5162# classes of partial linear spaces with (simple) point graph ptgraph=arg[1], 5163# and parameters s=arg[2],t=arg[3]. The default is that representatives 5164# for all isomorphism classes are returned. 5165# 5166# The integer argument nspaces=arg[4] is optional, and has 5167# default value -1, which means that representatives for all 5168# isomorphism classes are returned. If nspaces>=0 then exactly nspaces 5169# representatives are returned if there are at least nspaces isomorphism 5170# classes, otherwise representatives for all isomorphism classes are returned. 5171# 5172# In the output of this function, a partial linear space S is given 5173# by its incidence graph delta. The point-vertices of delta are 5174# 1,...,ptgraph.order, with the name of point-vertex i being the 5175# name of vertex i of ptgraph. A line-vertex of delta is named by a 5176# list (not necessarily ordered) of the point-vertex names for the points 5177# on that line. We warn that this is a *different* naming convention to 5178# versions of GRAPE before 4.1. The group delta.group associated 5179# with the incidence graph delta is the automorphism group of S 5180# acting on point-vertices and line-vertices, and preserving both sets. 5181# 5182# If arg[5] is bound then it controls the printlevel (default 0). 5183# Permitted values for arg[5] are 0,1,2. 5184# 5185# If arg[6] is bound then it is assumed to be a list (without repeats) 5186# of the (s+1)-cliques of ptgraph. If known, this can help the function 5187# to run faster. 5188# 5189local ptgraph,aut,X,printlevel,I,K,s,t,deg,search,cliques,nlines, 5190 ans,lines,pts,i,j,k,adj,nspaces,names; 5191ptgraph:=arg[1]; 5192s:=arg[2]; 5193t:=arg[3]; 5194if IsBound(arg[4]) then 5195 nspaces:=arg[4]; 5196else 5197 nspaces:=-1; 5198fi; 5199if IsBound(arg[5]) then 5200 printlevel:=arg[5]; 5201else 5202 printlevel:=0; 5203fi; 5204if not (IsGraph(ptgraph) and IsInt(s) and IsInt(t) and 5205 IsInt(nspaces) and IsInt(printlevel)) 5206 or (IsBound(arg[6]) and not IsList(arg[6])) then 5207 Error("usage: PartialLinearSpaces(", 5208 "<Graph>, <Int>, <Int>, [<Int>, [<Int>, [<List>]]])"); 5209fi; 5210if not printlevel in [0,1,2] then 5211 Error("<printlevel> must be 0, 1, or 2"); 5212fi; 5213if s<1 or t<1 then 5214 Error("<s> and <t> must be positive integers"); 5215fi; 5216if not IsSimpleGraph(ptgraph) then 5217 Error("<ptgraph> must be a simple graph"); 5218fi; 5219nlines:=(t+1)*ptgraph.order/(s+1); # number of lines 5220if not IsInt(nlines) or nspaces=0 then # no partial linear spaces 5221 return []; 5222fi; 5223if IsBound(arg[6]) then 5224 cliques:=arg[6]; 5225 if ForAny(cliques,x->not IsSSortedList(x) or Length(x)<>s+1) then 5226 Error("<arg[6]> incorrect"); 5227 fi; 5228 if Length(cliques) < nlines then 5229 return []; 5230 fi; 5231fi; 5232deg:=s*(t+1); 5233if ForAny(ptgraph.adjacencies,x->Length(x)<>deg) then 5234 return []; 5235fi; 5236aut:=AutGroupGraph(ptgraph); 5237ptgraph:=NewGroupGraph(aut,ptgraph); 5238if not IsBound(cliques) then 5239 if IsCompleteGraph(ptgraph) then 5240 cliques:=Combinations([1..ptgraph.order],s+1); 5241 else 5242 K:=CompleteSubgraphsOfGivenSize(ptgraph,s+1,true,false,true); 5243 cliques:=Concatenation(Orbits(ptgraph.group,K,OnSets)); 5244 fi; 5245 if Length(cliques) < nlines then 5246 return []; 5247 fi; 5248fi; 5249if nlines=t*(s+1)+1 then # line graph is complete graph 5250 adj:=function(x,y) return x<>y and Length(Intersection(x,y))<>1; end; 5251else 5252 adj:=function(x,y) return x<>y and Length(Intersection(x,y))>1; end; 5253fi; 5254X:=Graph(aut,cliques,OnSets,adj,true); 5255X.isSimple:=true; 5256Unbind(X.names); 5257StabChainOp(X.group,rec(limit:=Size(aut))); 5258# 5259# The set S of independent sets of size nlines in X is in 1-to-1 5260# correspondence with (the line sets of) the partial linear spaces 5261# with point graph ptgraph, and parameters s,t. 5262# 5263# Moreover, the orbits of X.group (induced by Aut(ptgraph)) on S 5264# are in 1-to-1 correspondence with the isomorphism classes 5265# of the partial linear spaces with point graph ptgraph, 5266# and parameters s,t. 5267# 5268# We shall classify S modulo X.group, in order to classify the 5269# required partial linear spaces up to isomorphism 5270# (except that we stop if nspaces>0 isomorphism classes are required 5271# and we find that number of them). 5272# 5273if Size(X.group) < Size(aut) then 5274 # 5275 # non-faithful action of aut on (s+1)-cliques, which is impossible if 5276 # a subset of these cliques form the line set of a partial linear space 5277 # with parameters s,t > 1. 5278 # 5279 return []; 5280fi; 5281if printlevel > 0 then 5282 Print("X.order=",X.order," VertexDegrees(X)=",VertexDegrees(X),"\n"); 5283fi; 5284I:=IndependentSet(ptgraph); 5285if printlevel > 0 then 5286 Print("Length(I)=",Length(I),"\n"); 5287fi; 5288I:=Concatenation(I,Difference([1..ptgraph.order],I)); 5289# 5290# We shall build the possible partial linear spaces by determining 5291# the possible line-sets through I[1],I[2],I[3],... (in order) 5292# and backtracking when necessary. 5293# 5294# It appears to be a good strategy to start I with the vertices 5295# of a maximal independent set of ptgraph. 5296# 5297ans:=[]; 5298# 5299# The "smallest" representatives of new isomorphism classes of the required 5300# partial linear spaces are put in ans as and when they are found. 5301# 5302 5303search := function ( i, sofar, live, H ) 5304# 5305# This is the function for the backtrack search. 5306# 5307# Given in sofar the vertices of X indexing the (s+1)-cliques 5308# forming all the lines through the points I[1],...,I[i-1], this function 5309# determines representatives for the new isomorphism classes 5310# (not in ans) of the required partial linear spaces S, 5311# such that the line-set of S contains all the cliques indexed by elements 5312# of sofar, but no clique not indexed by an element in the union of 5313# sofar and live. (The clique indexed by v is simply cliques[v].) 5314# 5315# This function assumes that sofar and live are disjoint sets, and 5316# that H <= X.group stabilizes each of sofar and live (setwise). 5317# It is also assumed that X.names is unbound, or X.names=[1..X.order]. 5318# 5319# Note: On entry to this function it is assumed that nspaces=-1 5320# or nspaces > 0. If nspaces > 0 then it is assumed that (on entry) 5321# nspaces isomorphism classes have not yet been found. 5322# If nspaces > 0 then this function terminates if nspaces 5323# isomorphism classes are found. 5324# It is also assumed that, on entry, the elements of ans are distinct 5325# and are the least lexicographically in their respective X.group-orbits. 5326# 5327local L, K, ind, k, forbid, nlinesreq, F, pointstocover, wts, ii, jj; 5328if printlevel > 1 then 5329 Print("\ni=",i," Size(H)=",Size(H)); 5330fi; 5331if Length(sofar)=nlines then 5332 # 5333 # partial linear space found. 5334 # check if its isomorphism class is new. 5335 # 5336 sofar:=SmallestImageSet(X.group,sofar); 5337 if not (sofar in ans) then 5338 # process new isomorphism class 5339 if printlevel > 1 then 5340 Print("\n",cliques{sofar},"\n"); 5341 fi; 5342 Add(ans,sofar); 5343 fi; 5344 return; 5345fi; 5346F:=Filtered(sofar,x->I[i] in cliques[x]); 5347nlinesreq:=(t+1)-Length(F); 5348# 5349# nlinesreq is the number of new lines through I[i] that 5350# must be found. 5351# 5352if nlinesreq=0 then 5353 search(i+1,sofar,live,H); 5354 return; 5355fi; 5356L:=Filtered( live, x->I[i] in cliques[x]); 5357if printlevel > 1 then 5358 Print(" nlinesreq=",nlinesreq," Length(L)=",Length(L)); 5359fi; 5360if Length(L) < nlinesreq then 5361 return; 5362fi; 5363H := Stabilizer( H, L, OnSets ); 5364ind := ComplementGraph(InducedSubgraph( X, L, H )); 5365pointstocover := 5366 Difference(Adjacency(ptgraph,I[i]),Union(List(F,x->cliques[x]))); 5367wts := List(L,x->ListWithIdenticalEntries(Length(pointstocover),0)); 5368for ii in [1..Length(L)] do 5369 for jj in cliques[L[ii]] do 5370 if jj<>I[i] then 5371 wts[ii][PositionSorted(pointstocover,jj)] := 1; 5372 fi; 5373 od; 5374od; 5375K := CompleteSubgraphsOfGivenSize( ind, 5376 ListWithIdenticalEntries(Length(pointstocover),1), 2, true, true, wts); 5377# 5378# K contains the sets of possible additional lines through I[i]. 5379# 5380if K = [] then 5381 return; 5382fi; 5383if printlevel > 1 then 5384 Print(" Length(K)=",Length(K)); 5385fi; 5386for k in K do 5387 L := ind.names{k}; 5388 forbid := Union( L, Union(List(L,x->Adjacency(X,x))) ); 5389 search( i+1, Union( sofar, L), Difference(live,forbid), 5390 Stabilizer(H,L,OnSets) ); 5391 if nspaces>=0 and Length(ans)=nspaces then 5392 return; 5393 fi; 5394od; 5395end; 5396 5397search(1,[],[1..X.order],X.group); 5398for i in [1..Length(ans)] do 5399 # 5400 # Determine the incidence graph of the partial linear space 5401 # whose lines are cliques{ans[i]}, and store the result in ans[i]. 5402 # 5403 pts:=List([1..ptgraph.order],x->[]); 5404 for j in ans[i] do 5405 for k in cliques[j] do 5406 AddSet(pts[k],j); 5407 od; 5408 od; 5409 aut:=Action(Stabilizer(X.group,ans[i],OnSets),pts,OnSets); 5410 lines:=SSortedList( cliques{ans[i]} ); 5411 ans[i]:=Graph(aut, 5412 Concatenation([1..ptgraph.order],lines), 5413 function(x,g) 5414 if IsInt(x) then 5415 return x^g; 5416 else 5417 return OnSets(x,g); 5418 fi; 5419 end, 5420 function(x,y) 5421 if IsInt(x) then 5422 return IsSSortedList(y) and x in y; 5423 else 5424 return IsInt(y) and y in x; 5425 fi; 5426 end, 5427 true); 5428 ans[i].isSimple:=true; 5429 # now rename the vertices of ans[i] 5430 names:=[]; 5431 for j in [1..ptgraph.order] do 5432 names[j]:=VertexName(ptgraph,j); 5433 od; 5434 for j in [ptgraph.order+1..ans[i].order] do 5435 names[j]:=List(ans[i].names[j],x->VertexName(ptgraph,x)); 5436 od; 5437 ans[i].names:=Immutable(names); 5438od; 5439return ans; 5440end); 5441