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 FundamentalDomainBieberbach.gi HAPcryst package Marc Roeder
##
##
##
#H @(#)$Id: FundamentalDomainBieberbach.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.("FundamentalDomainBieberbach_gi"):=
"@(#)$Id: FundamentalDomainBieberbach.gi, v 0.1.11 2013/10/27 18:31:09 gap Exp $";
#############################################################################
##
#O FundamentalDomainBieberbachGroupNC(<group>)
##
InstallMethod(FundamentalDomainBieberbachGroupNC,
[IsGroup],
function(group)
local center;
center:=0*[1..Size(Representative(group))-1];
return FundamentalDomainBieberbachGroupNC(center,group);
end);
#############################################################################
##
#O FundamentalDomainBieberbachGroup(<group>)
##
InstallMethod(FundamentalDomainBieberbachGroup,
[IsGroup],
function(group)
local center;
center:=0*[1..Size(Representative(group))-1];
return FundamentalDomainBieberbachGroup(center,group);
end);
#############################################################################
##
#O FundamentalDomainBieberbachGroupNC(<center>, <group>)
##
InstallMethod(FundamentalDomainBieberbachGroupNC,
[IsVector,IsGroup],
function(center,group)
local gram;
gram:=GramianOfAverageScalarProductFromFiniteMatrixGroup(PointGroup(group));
return FundamentalDomainBieberbachGroupNC(center,group,gram);
end);
#############################################################################
##
#O FundamentalDomainBieberbachGroup(<center>, <group>)
##
InstallMethod(FundamentalDomainBieberbachGroup,
[IsVector,IsGroup],
function(center,group)
local gram;
gram:=GramianOfAverageScalarProductFromFiniteMatrixGroup(PointGroup(group));
return FundamentalDomainBieberbachGroup(center,group,gram);
end);
#############################################################################
##
#O FundamentalDomainBieberbachGroup(<center>, <group>, <gram>)
##
InstallMethod(FundamentalDomainBieberbachGroup,
[IsVector,IsGroup,IsMatrix],
function(center,group,gram)
local needBieberbachTest, phi, dim, fd;
needBieberbachTest:=false;
if not (IsStandardSpaceGroup(group) and IsAffineCrystGroupOnRight(group))
then
Error("group must be a StandardSpaceGroup acting on right");
else
phi:=IsomorphismPcpGroup(group);
if phi=fail
then
Info(InfoHAPcryst,1,"Couldn't calculate pcp representation. Testing tortion freeness later");
needBieberbachTest:=true;
elif not IsAlmostBieberbachGroup(Image(phi))
then
Error("group must be a Bieberbach group");
fi;
fi;
dim:=Size(center);
if not (dim=Size(gram) and DimensionOfMatrixGroup(group)=dim+1)
then
Error("dimensions don't match");
fi;
if not IsTrivial(StabilizerOnSetsStandardSpaceGroup(group,[center]))
then
Error("group must be a Bieberbach group");
fi;
fd:=FundamentalDomainBieberbachGroupNC(center,group,gram);
if not needBieberbachTest or IsFundamentalDomainBieberbachGroup(fd,group)
then
return fd;
else
Error("group must be a Bieberbach group");
fi;
end);
#############################################################################
##
#O FundamentalDomainBieberbachGroupNC(<center>, <group>, <gram>)
##
InstallMethod(FundamentalDomainBieberbachGroupNC,
[IsVector,IsGroup,IsMatrix],
function(center,group,gram)
local initialInequalities, satisfiesAllInequalities,
newInequalitiesEuclidean, newInequalities, shuffledList,
dim, fdVol, ineqThreshold, isEuclidean, orbitpart,
partialFD, newinequalities, oldinequalities,
donevertices, allvertices, vertices, box, i, ilist,
inequalities, v, new;
###################
## A lot of helper functions:
###################
# This calculates the inequalities for the unit box around
initialInequalities:=function(center,gram)
local dim, directions, points, ineqs;
dim:=Size(center);
directions:=Union(IdentityMat(dim),-IdentityMat(dim));
points:=[];
points:=center+directions;
#points:=Set(ShiftedOrbitPart(center,orbitpart));
#UniteSet(points,center+directions);
# UniteSet(points,Union(List(directions,d->d+points)));
if gram<>IdentityMat(dim)
then
ineqs:=Set(points,
i->BisectorInequalityFromPointPair(center,i,gram));
else
ineqs:=Set(points,
i->BisectorInequalityFromPointPair(center,i));
fi;
RemoveSet(ineqs,fail);
return ineqs;
end;
satisfiesAllInequalities:=function(point,ineqs)
return ForAll(ineqs,i->WhichSideOfHyperplaneNC(point,i) in [0,1]);
end;
newInequalitiesEuclidean:=function(vertex,center,ineqs,vertices,box,group)
local dim, affinevertex, returnlist, g,
affimage, image, trans, t, gl, gli, gt;
dim:=Size(vertex);
affinevertex:=Concatenation(vertex,[1]);
returnlist:=[];
for g in PointGroupRepresentatives(group)
do
gl:=LinearPartOfAffineMatOnRight(g);
gli:=Inverse(gl);
gt:=g[dim+1]{[1..dim]};
affimage:=affinevertex*g;
image:=affimage{[1..dim]};
###
trans:=TranslationsToBox(image,box);
for t in trans
do
if not image+t in vertices
then
if satisfiesAllInequalities(image+t,ineqs)
then
Add(returnlist,center*gl+gt+t);
Add(returnlist,(center-gt-t)*gli);
fi;
fi;
od;
od;
if returnlist=[]
then
return [];
fi;
returnlist:=Set(returnlist);
RemoveSet(returnlist,center);
returnlist:=Set(returnlist,i->BisectorInequalityFromPointPair(center,i));
SubtractSet(returnlist,ineqs);
return Filtered(returnlist,i->ForAny(vertices,v->WhichSideOfHyperplaneNC(v,i)=-1));
end;
newInequalities:=function(vertex,center,ineqs,vertices,box,group,gram)
local dim, affinevertex, returnlist, g, gl, gli, gt,
affimage, image, trans, t;
dim:=Size(vertex);
affinevertex:=Concatenation(vertex,[1]);
returnlist:=[];
for g in PointGroupRepresentatives(group)
do
gl:=LinearPartOfAffineMatOnRight(g);
gli:=Inverse(gl);
gt:=g[dim+1]{[1..dim]};
affimage:=affinevertex*g;
image:=affimage{[1..dim]};
###
trans:=TranslationsToBox(image,box);
for t in trans
do
if not image+t in vertices
then
if satisfiesAllInequalities(image+t,ineqs)
then
Add(returnlist,center*gl+gt+t);
Add(returnlist,(center-gt-t)*gli);
fi;
fi;
od;
od;
if returnlist=[]
then
return [];
fi;
returnlist:=Set(returnlist);
RemoveSet(returnlist,center);
returnlist:=Set(returnlist,i->BisectorInequalityFromPointPair(center,i,gram));
SubtractSet(returnlist,ineqs);
return Filtered(returnlist,i->ForAny(vertices,v->WhichSideOfHyperplaneNC(v,i)=-1));
end;
shuffledList:=function(list)
local positionlist, indexlist, i, nextpos;
if Size(list)=1
then return list;
fi;
positionlist:=EmptyPlist(Size(list));
indexlist:=[1..Size(list)];
for i in [Size(list),Size(list)-1..1]
do
nextpos:=Random(1,i);
Add(positionlist,indexlist[nextpos]);
Remove(indexlist,nextpos);
od;
return list{positionlist};
end;
###################
## helper functions done. Let's go:
###################
dim:=Size(center);
## The volume of a fundamental domain is seemingly
## well-known:
#fdVol:=1/Order(PointGroup(group));
## This is only used in the first step.
## Polymake calculates a triangulation using the beneath-beyond
## algorithm for this.
## This seams to be inefficient in high dimensions, especially
## with complicated polyhedrae.
##
ineqThreshold:=ValueOption("ineqThreshold");
if ineqThreshold=fail
then
ineqThreshold:=200;
fi;
isEuclidean:=(gram=IdentityMat(dim));
orbitpart:=OrbitStabilizerInUnitCubeOnRight(group,center).orbit;;
newinequalities:=initialInequalities(center,gram);;
partialFD:=CreatePolymakeObject("partialFD",
POLYMAKE_DATA_DIR,
["polytope","2.3","RationalPolytope"]
);
AppendToPolymakeObject(partialFD,
ConvertMatrixToPolymakeString("FACETS",newinequalities)
);
# AppendInequalitiesToPolymakeObject(partialFD,newinequalities);
## slightly experimental: remove if this causes trouble:
Polymake(partialFD,"VERTICES FACETS");
newinequalities:=[];
oldinequalities:=[];
donevertices:=[];
allvertices:=[];
repeat
allvertices:=AsSet(Polymake(partialFD,"VERTICES"));
donevertices:=Set(donevertices);
vertices:=Difference(allvertices,donevertices);
##
box:=[1..dim];
for i in [1..dim]
do
ilist:=allvertices{[1..Size(allvertices)]}[i];
box[i]:=[Minimum(ilist),Maximum(ilist)];
od;
# Sort(donevertices);
# vertices:=(Filtered(Polymake(partialFD,"VERTICES"),v->not v in donevertices));
# vertices:=Difference(allvertices,donevertices);
if vertices=[]
then
Info(InfoHAPcryst,2,"out of vertices");
return partialFD;
fi;
inequalities:=Set(Polymake(partialFD,"FACETS"));
Info(InfoHAPcryst,3,"v: ",Size(allvertices),"/",Size(vertices)," f:", Size(inequalities)," \c");
newinequalities:=[];
repeat
v:=Remove(vertices,Random(1,Size(vertices)));
# v:=Remove(vertices);
Add(donevertices,v);
if satisfiesAllInequalities(v,newinequalities)
then
if isEuclidean
then
new:= newInequalitiesEuclidean(v,center,inequalities,
allvertices,box,group);
else
new:= newInequalities(v,center,inequalities,
allvertices,box,group,gram);
fi;
UniteSet(inequalities,new);
Append(newinequalities,new);
fi;
until vertices=[] or Size(newinequalities)>ineqThreshold;
Info(InfoHAPcryst,3,"new: ",Size(newinequalities));
if newinequalities<>[]
then
ClearPolymakeObject(partialFD,["polytope","2.3","RationalPolytope"]);
# AppendInequalitiesToPolymakeObject(partialFD,shuffledList(inequalities));;
AppendInequalitiesToPolymakeObject(partialFD,inequalities);;
Polymake(partialFD,"VERTICES FACETS");
fi;
until newinequalities=[];
Info(InfoHAPcryst,2,"out of inequalities");
return partialFD;
end);