Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
| Download
GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Views: 418346#############################################################################
##
#W FaceLatticeAndBoundaryBieberbachGroup.gi HAPcryst package Marc Roeder
##
##
##
#H @(#)$Id: FaceLatticeAndBoundaryBieberbachGroup.gi, v 0.1.11 2013/10/27 18:31:09 gap Exp $
##
#Y Copyright (C) 2006 Marc Roeder
#Y
#Y This program is free software; you can redistribute it and/or
#Y modify it under the terms of the GNU General Public License
#Y as published by the Free Software Foundation; either version 2
#Y of the License, or (at your option) any later version.
#Y
#Y This program is distributed in the hope that it will be useful,
#Y but WITHOUT ANY WARRANTY; without even the implied warranty of
#Y MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#Y GNU General Public License for more details.
#Y
#Y You should have received a copy of the GNU General Public License
#Y along with this program; if not, write to the Free Software
#Y Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
##
Revision.("FaceLatticeAndBoundaryBieberbachGroup_gi"):=
"@(#)$Id: FaceLatticeAndBoundaryBieberbachGroup.gi, v 0.1.11 2013/10/27 18:31:09 gap Exp $";
InstallMethod(FaceLatticeAndBoundaryBieberbachGroup,
[IsPolymakeObject,IsGroup],
function(poly,group)
local removeSomeFaces, orbitDecompositionAsIndices,
changeHasseEntries, reformatHD, reformatHD_GroupRing,
polymakeFaceLattice, calculateBoundary, dim, starttime,
vertices, hasse, codim, initialfacenumber, timetmp,
faceOrbits, tmp, k, j, groupring, elts;
removeSomeFaces:=function(upfaces,faces)
local upfaceentries, faceindex, face;
upfaceentries:=AsSet(Flat(upfaces));
for faceindex in [1..Size(faces)]
do
face:=faces[faceindex];
if not (IsSubset(upfaceentries,face) and ForAny(upfaces,f->IsSubset(f,face)))
then
Unbind(faces[faceindex]);
fi;
od;
end;
orbitDecompositionAsIndices:=function(facelist,vertexlist,group)
local sortedvertices, sortperm, vertexpositionlookup,
todolist, setfacelist, returnlist, thissize, list,
orbit, thisorbit, fpos, fpositions;
if not ForAll(facelist,IsSet)
then
Error("face list must be a list of sets");
fi;
sortedvertices:=ShallowCopy(vertexlist);
sortperm:=Sortex(sortedvertices);
vertexpositionlookup:=Permuted([1..Size(sortedvertices)],sortperm);
# convert all faces to vectors:
todolist:=Set(facelist,i->Set(vertexlist{i}));
if IsSet(facelist)
then
setfacelist:=facelist;
else
setfacelist:=AsSet(facelist);
fi;
returnlist:=[];
while todolist<>[]
do
thissize:=Size(todolist[1]);
list:=Filtered(todolist,o->Size(o)=thissize);
SubtractSet(todolist,list);
while list<>[]
do
orbit:=OrbitPartAndRepresentativesInFacesStandardSpaceGroup(group,list[1],list);
thisorbit:=[];
for fpos in [1..Size(orbit)]
do
fpositions:=Set(orbit[fpos][1],i->vertexpositionlookup[PositionSet(sortedvertices,i)]);
if fpositions in setfacelist
then
Add(thisorbit,[Position(facelist,fpositions),orbit[fpos][2]]);
fi;
od;
Add(returnlist,thisorbit);
SubtractSet(list,List(orbit,i->i[1]));
od;
od;
return returnlist;
end;
changeHasseEntries:=function(orbitsandreps,codim,vertices,hasse,group)
local genIndices, faces, nrupfaces, upfaces,
newhasseentry, idMat, orbit, gen, genIndex,
upindices, upface, genvertices, o, oface, entry;
genIndices:=[];
faces:=hasse[codim+1];
nrupfaces:=Size(hasse[codim]);
upfaces:=hasse[codim]{[1..nrupfaces]}[1];
newhasseentry:=[];
idMat:=IdentityMat(Size(hasse));
for orbit in orbitsandreps
do
gen:=faces[orbit[1][1]];
Add(newhasseentry,[gen,[]]);
genIndex:=Size(newhasseentry);
upindices:=Filtered([1..nrupfaces],
i->IsSubset(upfaces[i],gen));
#update boundary for faces containing new generator:
for upface in upindices
do
AddSet(hasse[codim][upface][2],[genIndex,idMat]);
od;
#replace all faces in the generator-orbit with representatives:
if Size(orbit)>1
then
genvertices:=Set(vertices{gen});
for o in orbit{[2..Size(orbit)]}
do
oface:=faces[o[1]];
upindices:=Filtered([1..nrupfaces],
i->IsSubset(upfaces[i],oface));
entry:=[genIndex,o[2]];
for upface in upindices
do
AddSet(hasse[codim][upface][2],entry);
od;
od;
fi;
od;
return newhasseentry;
end;
#this replaces the matrices in the partial Hasse diagram with integers and
# returns a list of matrices. This is done to comply with the structure
# of a resolution in HAP.
reformatHD:=function(hasse)
local elts, faces, face, line;
elts:=[];
for faces in hasse
do
for face in faces
do
for line in face[2]
do
Add(elts,line[2]);
od;
od;
od;
# if elts<>[]
# then
# elts:=Set(elts);
# AddSet(elts,IdentityMat(Size(elts[1])));
# for faces in hasse
# do
# for face in faces
# do
# for line in face[2]
# do
# if not IsInt(line[2])
# then
# line[2]:=PositionSet(elts,line[2]);
# fi;
# od;
# od;
# od;
# fi;
return elts;
end;
## This does not generate group ring elements. The group ring is
## generated in the resolution generator. So here, we only generate
## the coefficient-group element pairs.
reformatHD_GroupRing:=function(groupring,hasse)
local zero, family, groupelements, dim, face, term,
firstpos, firstgen, position, firspos, one,
onepos;
zero:=Zero(groupring);
family:=FamilyObj(zero);
groupelements:=[];
for dim in [2..Size(hasse)]
do
for face in hasse[dim]
do
for term in face[2]
do
if IsMatrix(term[2])
then
Add(groupelements,term[2]);
term[2]:=SignInt(term[1])*
ElementOfMagmaRing(family,0,[1],[term[2]]);
term[1]:=AbsInt(term[1]);
fi;
od;
Sort(face[2]);
firstpos:=1;
firstgen:=face[2][firstpos][1];
for position in [2..Size(face[2])]
do
if face[2][position][1]=firstgen
then
face[2][firstpos][2]:=face[2][firstpos][2]
+face[2][position][2];
Unbind(face[2][position]);
else
firstpos:=position;
firstgen:=face[2][position][1];
fi;
od;
face[2]:=Compacted(face[2]);
# face[2]:=vector;
od;
od;
groupelements:=Set(groupelements);
one:=One(UnderlyingMagma(groupring));
if groupelements[1]<>one
then
onepos:=Position(groupelements,one);
groupelements[onepos]:=groupelements[1];
groupelements[1]:=one;
fi;
return groupelements;
end;
######## patch for next genertion polymake ###
## first part
polymakeFaceLattice:=function(polygon)
local faces,dimlist, dim, lattice, i;
faces:=Polymake(polygon,"FACES");
dimlist:=Polymake(polygon,"DIMS");
dim := Length(dimlist);
lattice:=EmptyPlist(dim);
## walk through the face list from top to bottom
# the dim list contains the maximal index of a face in the
# respective dimension. The top face containing all vertices
# has no entry in the dimlist, so we add it beforehand.
##
lattice[dim] := [StructuralCopy(faces[Length(faces)])];
for i in [1..dim-1] do
lattice[i] := StructuralCopy(faces{[dimlist[i]+1..dimlist[i+1]]});
od;
## The reverse probably runs every time
## panic-coding? I don't know, let's keep it...
if Length(lattice[1][1]) = 1 then
lattice := Reversed(lattice);
fi;
return lattice;
end;
### end of patch
######################################################################
######################################################################
##################################################
## THIS IS JUST FOR BIEBERBACH GROUPS!
##
## here, hasse has to be a list contining the generators of
## dimension $k$ in the $k+1$st entry.
## each of the genertors already knows it's unoriented boundary.
##
## Signs are now assigned to the boundary to define a proper
## boundary homomorphism.
##
calculateBoundary:=function(k,j,hasse)
local dirLess, boundaryFromPair, face, facebound,
linesdone, linestodo, linesinpoint, pointset,
idMat, line, linemat, points, point, pos,
firstlinepos, firstline, linesforthisrun,
nextlines, pointsdone, linebound, lineboundDirless,
dirlessline, pointsign, otherline, otherlinepos,
otherlinebound, orientedOtherLine;
dirLess:=function(face)
return [AbsInt(face[1]),face[2]];
end;
boundaryFromPair:=function(dim,pair,hasse)
local gen, bound;
gen:=hasse[dim+1][AbsInt(pair[1])];
bound:=gen[2];
return List(bound,b->[SignInt(pair[1])*b[1],b[2]*pair[2]]);
end;
face:=hasse[k+1][j];
################################
## face is of the form
## [ <verts>, <bound>]
## where
## <verts> is a list of integers enumerating the vertices of the
## fundamental domain
##
## <bound> is a list of pairs [<downfaceindex>,<mapping>]
## with an integer <downfaceindex> denoting the position of a face
## in the list hasse[k]. <mapping> is an affine matrix taking the ##
## hasse[k][downfaceindex] to the desired face. ##
## ##
## ###########################
if k=1 # First, define the boundary of 1-faces
then
facebound:=face[2];
if facebound[1]>facebound[2]
then
linesdone:=[facebound[1],[-facebound[2][1],facebound[2][2]]];
else
linesdone:=[[-facebound[1][1],facebound[1][2]],facebound[2]];
fi;
else
#
#now the other faces. Note that in an n-face every
# n-2 face is contained in exactly 2 n-1 faces.
##############################
linestodo:=Set(face[2],i->dirLess(i));
# # the edges which are generators for this module:
# GenLinesTodo:=Filtered(linestodo,i->i[2]=idMat);
# # edges which are images under generators for this module:
# ImLinesTodo:=Filtered(linestodo,i->i[2]<>idMat);
# As the group is Bieberbach, each point has a unique name.
# generate a list to look up the incident lines for each point:
linesinpoint:=[];
pointset:=[];
idMat:=IdentityMat(Size(hasse));
for line in linestodo
# form: line=[<pos>,<mat>]
# <pos>: generator position
# <mat>: matrix taking generator to actual line.
do
linemat:=line[2];
points:=List(hasse[k][line[1]][2],p->[AbsInt(p[1]),p[2]*linemat]);
for point in points
#form of point as form of line but one dimension lower
do
if not point in pointset
then
Add(linesinpoint,[point,[line]]);
else
pos:=PositionSet(pointset,point);
AddSet(linesinpoint[pos][2],line);
fi;
od;
UniteSet(pointset,points);
Sort(linesinpoint);
od;
firstlinepos:=PositionProperty(linestodo,i->i[2]=idMat);
if firstlinepos=fail
then
firstline:=Remove(linestodo);
else
firstline:=linestodo[firstlinepos];
Remove(linestodo,firstlinepos);
fi;
linesforthisrun:=[firstline];
nextlines:=[];
pointsdone:=[];
linesdone:=[];
repeat
for line in linesforthisrun
do
linesdone:=Set(linesdone);
if not ([line[1],line[2]] in linesdone
or [-line[1],line[2]] in linesdone)
then
Add(linesdone,line);
fi;
linebound:=boundaryFromPair(k-1,line,hasse);
lineboundDirless:=List(linebound,dirLess);
SortParallel(lineboundDirless,linebound);
dirlessline:=dirLess(line);
points:=Difference(lineboundDirless,pointsdone);
for point in points
do
Add(pointsdone,point);
pointsign:=SignInt(linebound[PositionSet(lineboundDirless,point)][1]);
otherline:=linesinpoint[PositionSet(pointset,point)];
otherline:=First(otherline[2],l->l<>dirlessline);
if otherline in linestodo
then
RemoveSet(linestodo,otherline);
otherlinebound:=boundaryFromPair(k-1,otherline,hasse);
if pointsign=SignInt(First(otherlinebound,i->dirLess(i)=point)[1])
then
orientedOtherLine:=[-otherline[1],otherline[2]];
else
orientedOtherLine:=otherline;
fi;
Add(linesdone,orientedOtherLine);
Add(nextlines,orientedOtherLine);
fi;
od;
od;
linesforthisrun:=ShallowCopy(nextlines);
nextlines:=[];
until linesforthisrun=[];
fi;
## return new boundary entry:
return Set(linesdone);
end;
######################################################################
######################################################################
if not IsStandardSpaceGroup(group) and IsAffineCrystGroupOnRight(group)
then
Error("<group> must be a StandardSpaceGroup acting on right");
fi;
dim:=DimensionOfMatrixGroup(group)-1;
starttime:=Runtime();
## patch for next generation polymake:
# second part. If you have the old polymake version,
# use "FACE_LATTICE" instead of "FACES DIMS" but keep the VERTICES!
Polymake(poly,"FACES DIMS VERTICES");
vertices:=Polymake(poly,"VERTICES");
MakeImmutable(vertices);
if not dim=Size(vertices[1]-1)
then
Error("group and polyhedron do not match");
fi;
##
# Add initial (top) node:
# Note that the face lattice is generated top- down. So the
# order has to be reversed to be ascending in dimension.
##
## patch for next generation polymake:
# third part. If you have the old polymake version,
# use the line where "Polymake" is called with "FACE_LATTICE"
hasse:=Concatenation([[[[1..Size(vertices)],[]]]],
#StructuralCopy(Polymake(poly,"FACE_LATTICE"){[1..dim]}));
StructuralCopy(polymakeFaceLattice(poly){[1..dim]}));
#moduleGenerators:=List([1..dim+1],i->[]);
#moduleGenerators is a list of list. The i^{th} entry contains the
# positions of the module generators of dimension i in the Hasse diagram.
# moduleGenerators[1]:=[1];# so this means [hasse[1][1]];
#Now the next nodes:
for codim in [1..dim]
do
#At this stage, hasse[codim+1] is a list of faces.
initialfacenumber:=Size(hasse[codim+1]);##########
timetmp:=Runtime();########
removeSomeFaces(List(hasse[codim],i->i[1]),hasse[codim+1]);
hasse[codim+1]:=Set(hasse[codim+1]);
faceOrbits:=orbitDecompositionAsIndices(hasse[codim+1],vertices,group);
hasse[codim+1]:=changeHasseEntries(faceOrbits,codim,vertices,hasse,group);
Info(InfoHAPcryst,3,codim,"(",Size(hasse[codim+1]),"/",initialfacenumber,"):",StringTime(Runtime()-timetmp));
od;
for codim in [1..Int((dim+1)/2)]
do
tmp:=hasse[codim];
hasse[codim]:=hasse[dim+2-codim];
hasse[dim+2-codim]:=tmp;
od;
Info(InfoHAPcryst,3,"Face lattice done (",StringTime(Runtime()-starttime),"). Calculating boundary ");
for k in [1..dim]
do
for j in [1..Size(hasse[k+1])]
do
hasse[k+1][j][2]:=calculateBoundary(k,j,hasse);
od;
od;
if fail in Flat(hasse)
then
Error("boundary generation failed");
fi;
Info(InfoHAPcryst,3,"done (",StringTime(Runtime()-timetmp),") Reformating...\c");
# elts:=reformatHD(hasse);
groupring:=GroupRing(Integers,group);
elts:=reformatHD_GroupRing(groupring,hasse);
# elts:=[IdentityMat(dim+1)];
return rec(hasse:=hasse,elts:=elts,groupring:=groupring);
end);