From 60713f4e05a37bf8e598a141a460dd2fc3ac04d3 Mon Sep 17 00:00:00 2001 From: neozhaoliang Date: Wed, 17 Jan 2024 22:01:03 +0800 Subject: [PATCH] add polytopes back --- src/aztec/run_domino_shuffling_animation.py | 4 +- src/polytopes/README.md | 5 + src/polytopes/example_curved_polychora.py | 299 ++++++++ src/polytopes/example_polytope_animation.py | 264 +++++++ .../example_run_coset_enumeration.py | 174 +++++ src/polytopes/polytopes/__init__.py | 0 src/polytopes/polytopes/coxeter_plane.py | 73 ++ src/polytopes/polytopes/helpers.py | 123 ++++ src/polytopes/polytopes/models.py | 671 ++++++++++++++++++ src/polytopes/polytopes/povray.py | 62 ++ src/polytopes/polytopes/todd_coxeter.py | 317 +++++++++ src/polytopes/povray/mathutils.inc | 90 +++ src/polytopes/povray/polychora_curved.pov | 222 ++++++ src/polytopes/povray/polyhedra_animation.pov | 130 ++++ src/polytopes/povray/polytope_animation.pov | 141 ++++ src/polytopes/tc_examples/B24.yaml | 15 + src/polytopes/tc_examples/Cube.yaml | 11 + src/polytopes/tc_examples/F27.yaml | 13 + src/polytopes/tc_examples/G8723.yaml | 14 + src/polytopes/tc_examples/HNO_1.yaml | 12 + src/polytopes/tc_examples/M12.yaml | 12 + 21 files changed, 2649 insertions(+), 3 deletions(-) create mode 100644 src/polytopes/README.md create mode 100644 src/polytopes/example_curved_polychora.py create mode 100644 src/polytopes/example_polytope_animation.py create mode 100644 src/polytopes/example_run_coset_enumeration.py create mode 100644 src/polytopes/polytopes/__init__.py create mode 100644 src/polytopes/polytopes/coxeter_plane.py create mode 100644 src/polytopes/polytopes/helpers.py create mode 100644 src/polytopes/polytopes/models.py create mode 100644 src/polytopes/polytopes/povray.py create mode 100644 src/polytopes/polytopes/todd_coxeter.py create mode 100644 src/polytopes/povray/mathutils.inc create mode 100644 src/polytopes/povray/polychora_curved.pov create mode 100644 src/polytopes/povray/polyhedra_animation.pov create mode 100644 src/polytopes/povray/polytope_animation.pov create mode 100644 src/polytopes/tc_examples/B24.yaml create mode 100644 src/polytopes/tc_examples/Cube.yaml create mode 100644 src/polytopes/tc_examples/F27.yaml create mode 100644 src/polytopes/tc_examples/G8723.yaml create mode 100644 src/polytopes/tc_examples/HNO_1.yaml create mode 100644 src/polytopes/tc_examples/M12.yaml diff --git a/src/aztec/run_domino_shuffling_animation.py b/src/aztec/run_domino_shuffling_animation.py index 28df496..87a173e 100644 --- a/src/aztec/run_domino_shuffling_animation.py +++ b/src/aztec/run_domino_shuffling_animation.py @@ -69,9 +69,7 @@ def make_animation(order, size, filename): parser.add_argument( "-order", metavar="o", type=int, default=30, help="order of aztec diamond" ) - parser.add_argument( - "-size", metavar="s", type=int, default=400, help="image size" - ) + parser.add_argument("-size", metavar="s", type=int, default=400, help="image size") parser.add_argument( "-filename", metavar="f", default="domino_shuffling.gif", help="output filename" ) diff --git a/src/polytopes/README.md b/src/polytopes/README.md new file mode 100644 index 0000000..f187faa --- /dev/null +++ b/src/polytopes/README.md @@ -0,0 +1,5 @@ ++ `example_curved_polychora.py`: Render 3d/4d polytopes in POV-Ray (projected to the 3/4-sphere). ++ `example_polyhedra_mirrors_shader.py`: A shader program draws polytopes with mirrors. ++ `example_polytope_animation.py`: Render 3d/4d polytopes in POV-Ray (projected to the 3d space). ++ `example_run_coset_enumeration.py`: Compute coset table of a finitely represented group. ++ `example_wythoff_shader_animation.py`: A shader program exported from Matt Zucker's shadertoy work. \ No newline at end of file diff --git a/src/polytopes/example_curved_polychora.py b/src/polytopes/example_curved_polychora.py new file mode 100644 index 0000000..a2a6c4d --- /dev/null +++ b/src/polytopes/example_curved_polychora.py @@ -0,0 +1,299 @@ +""" +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Render curved 4d polychoron examples in POV-Ray +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +:copyright (c) 2018 by Zhao Liang. +""" +import os +import subprocess +from fractions import Fraction + +from polytopes.models import Polychora +from polytopes.povray import pov_index_array1d, pov_vector + + +IMAGE_DIR = "polychora_frames" # directory to save the frames +POV_EXE = "povray" # POV-Ray exe binary +SCENE_FILE = "polychora_curved.pov" # the main scene file +IMAGE_SIZE = 600 # image size in pixels +FRAMES = 1 # number of frames +IMAGE_QUALITY_LEVEL = 11 # between 0-11 +SUPER_SAMPLING_LEVEL = 5 # between 1-9 +ANTIALIASING_LEVEL = 0.001 # lower for better quality +DATAFILE_NAME = "polychora-data.inc" # export data to this file + +data_file = os.path.join(os.getcwd(), "povray", DATAFILE_NAME) + +if not os.path.exists(IMAGE_DIR): + os.makedirs(IMAGE_DIR) + +POV_COMMAND = ( + " cd povray && " + + " {} +I{}".format(POV_EXE, SCENE_FILE) + + " +W{} +H{}".format(IMAGE_SIZE, IMAGE_SIZE) + + " +Q{}".format(IMAGE_QUALITY_LEVEL) + + " +A{}".format(ANTIALIASING_LEVEL) + + " +R{}".format(SUPER_SAMPLING_LEVEL) + + " +KFI0" + + " +KFF{}".format(FRAMES - 1) + + " -V" + + " +O../{}/".format(IMAGE_DIR) + + "{}" +) + +POV_TEMPLATE = """ +#declare bg_color = {}; +#declare vertex_size = {}; +#declare edge_size = {}; +#declare camera_loc = {}; +#declare obj_rotation = {}; +#declare size_func = {}; +#declare face_max= {}; +#declare face_min = {}; +#declare face_index = {}; +#declare use_area_light = {}; + +// this macro is used for adjusting the size of edges +// according to their positions in the space. +#macro get_size(q) + #local len = vlength(q); + #if (size_func = 0) + #local len = (1.0 + len * len) / 4; + #else #if (size_func = 1) + #local len = 2.0 * log(2.0 + len * len); + #else + #local len = 2.0 * log(1.13 + len * len); + #end + #end + len +#end + +#macro choose_face(i, face_size) + #local chosen = false; + #for (ind, 0, dimension_size(face_index, 1) - 1) + #if (i = face_index[ind]) + #if (face_size > face_min & face_size < face_max) + #local chosen = true; + #end + #end + #end + chosen +#end + +#declare vertices = {}; + +#declare edges = {}; + +#declare faces = {}; +""" + + +def draw( + coxeter_diagram, + trunc_type, + extra_relations=(), + description="polychora", + bg_color="SkyBlue", + camera_loc=(0, 0, 30), + rotation=(0, 0, 0), + vertex_size=0.04, + edge_size=0.02, + size_func=0, + face_index=(0,), + face_max=3.0, + face_min=0.5, + use_area_light=False, +): + """ + Export data to povray .inc file and call the rendering process. + + :param camera_loc: + location of the camera. + + :param size_func: + choose which sizing funcion to use, currently only 0, 1, 2. + + :param face_index: + a list controls which types of faces are shown. + + :param face_max: + only faces smaller than this threshold are shown. + + :param face_min: + only faces larger than this threshold are shown. + """ + P = Polychora(coxeter_diagram, trunc_type, extra_relations) + P.build_geometry() + vert_data, edge_data, face_data = P.get_povray_data() + with open(data_file, "w") as f: + f.write( + POV_TEMPLATE.format( + bg_color, + vertex_size, + edge_size, + pov_vector(camera_loc), + pov_vector(rotation), + size_func, + face_max, + face_min, + pov_index_array1d(face_index), + int(use_area_light), + vert_data, + edge_data, + face_data, + ) + ) + + print( + "rendering {}: {} vertices, {} edges, {} faces".format( + description, P.num_vertices, P.num_edges, P.num_faces + ) + ) + + process = subprocess.Popen( + POV_COMMAND.format(description), + shell=True, + stderr=subprocess.PIPE, + stdin=subprocess.PIPE, + stdout=subprocess.PIPE, + ) + + _, err = process.communicate() + if process.returncode: + print(type(err), err) + raise IOError("POVRay error: " + err.decode("ascii")) + + +def main(): + draw( + (3, 2, 2, 3, 2, 3), + (1, 0, 0, 0), + description="5-cell", + camera_loc=(0, 0, 120), + vertex_size=0.08, + rotation=(-30, 60, 0), + edge_size=0.04, + size_func=1, + ) + + draw( + (4, 2, 2, 3, 2, 3), + (1, 0, 0, 0), + description="4d-cube", + camera_loc=(0, 0, 130), + vertex_size=0.06, + rotation=(60, 0, 0), + edge_size=0.03, + size_func=1, + face_min=0.2, + face_max=0.8, + ) + + draw( + (3, 2, 2, 3, 2, 4), + (1, 0, 0, 0), + description="16-cell", + camera_loc=(0, 0, 160), + vertex_size=0.08, + edge_size=0.03, + size_func=2, + face_min=1.0, + face_max=1.2, + ) + + draw( + (3, 2, 2, 4, 2, 3), + (1, 0, 0, 0), + description="24-cell", + camera_loc=(0, 0, 200), + vertex_size=0.06, + edge_size=0.04, + size_func=2, + face_min=0.2, + face_max=0.8, + ) + + draw( + (5, 2, 2, 3, 2, 3), + (1, 0, 0, 0), + description="120-cell", + camera_loc=(0, 0, 400), + vertex_size=0.05, + edge_size=0.025, + size_func=0, + face_min=3.0, + face_max=100.0, + ) + + draw( + (3, 2, 2, 3, 2, 5), + (1, 0, 0, 0), + description="600-cell", + bg_color="White", + camera_loc=(0, 0, 500), + vertex_size=0.12, + edge_size=0.04, + size_func=2, + face_max=4.0, + face_min=3.0, + ) + + draw( + (3, 2, 2, 3, 2, 4), + (1, 0, 0, 1), + description="runcinated-16-cell", + bg_color="White", + camera_loc=(0, 0, 450), + vertex_size=0.1, + face_index=(0, 1, 2, 3), + edge_size=0.04, + size_func=1, + face_min=0, + face_max=3, + ) + + draw( + (5, 2, 2, 3, 2, 3), + (1, 0, 0, 1), + description="runcinated-120-cell", + camera_loc=(0, 0, 360), + vertex_size=0.028, + edge_size=0.014, + face_min=20, + ) + + # this is the settings I used to render the movie at + # http://pywonderland.com/images/polytopes/rectified-grand-stellated-120-cell.mp4 + # (the parameters are not exactly the same but very similar) + # take quite a while to render. + draw( + (Fraction(5, 2), 2, 2, 5, 2, Fraction(5, 2)), + (0, 1, 0, 0), + extra_relations=((0, 1, 2, 1) * 3, (1, 2, 3, 2) * 3), + description="rectified-grand-stellated-120-cell", + size_func=1, + vertex_size=0.06, + edge_size=0.03, + use_area_light=1, + face_max=0.0, + camera_loc=(0, 0, 400), + ) + + draw( + (Fraction(5, 2), 2, 2, 5, 2, Fraction(5, 2)), + (0, 1, 1, 0), + extra_relations=((0, 1, 2, 1) * 3, (1, 2, 3, 2) * 3), + description="runcinated-grand-stellated-120-cell", + size_func=1, + vertex_size=0.06, + edge_size=0.03, + use_area_light=1, + face_min=0.3, + face_max=1, + camera_loc=(0, 0, 380), + ) + + +if __name__ == "__main__": + main() diff --git a/src/polytopes/example_polytope_animation.py b/src/polytopes/example_polytope_animation.py new file mode 100644 index 0000000..4e7b15d --- /dev/null +++ b/src/polytopes/example_polytope_animation.py @@ -0,0 +1,264 @@ +""" +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Make animations of rotating polytopes +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This script computes the data of a given polytope and writes it +into a POV-Ray .inc file, then automatically calls POV-Ray +to render the frames and calls FFmpeg to convert the frames to +a mp4 movie. You need to have POV-Ray and FFmpeg installed and +set the path to their executables in `POV_EXE` and `FFMPEG_EXE`. + +:copyright (c) 2018 by Zhao Liang. +""" +import os +import subprocess +from fractions import Fraction + +import polytopes.models as models + +IMAGE_DIR = "polytope_frames" # directory to save the frames +POV_EXE = "povray" # povray command +FFMPEG_EXE = "ffmpeg" # ffmpeg command +FRAMES = 1 # number of frames (120 is quite good) +IMAGE_SIZE = 500 # image size +SUPERSAMPLING_LEVEL = 1 # supersampling level +DATAFILE_NAME = "polytope-data.inc" # export data to this file + +data_file = os.path.join(os.getcwd(), "povray", DATAFILE_NAME) + +if not os.path.exists(IMAGE_DIR): + os.makedirs(IMAGE_DIR) + +# POV-Ray command line options +POV_COMMAND = ( + " cd povray && {}".format(POV_EXE) + + " {}" + + " +W{} +H{}".format(IMAGE_SIZE, IMAGE_SIZE) + + " +Q11 +A0.003 +R{}".format(SUPERSAMPLING_LEVEL) + + " +KFI0" + + " +KFF{}".format(FRAMES - 1) + + " -V" + + " +O../{}/".format(IMAGE_DIR) + + "{}" +) + +# FFmpeg command line options +FFMPEG_COMMAND = ( + " cd {} && ".format(IMAGE_DIR) + + " {} -framerate 12".format(FFMPEG_EXE) + + " -y" + + " -i {}" + + "%0{}d.png".format(len(str(FRAMES - 1))) + + " -crf 18 -c:v libx264 -pix_fmt yuv420p" + + " ../{}.mp4" +) + +POV_TEMPLATE = """ +#declare vertices = {}; + +#declare edges = {}; + +#declare faces = {}; +""" + + +def anim( + coxeter_diagram, + trunc_type, + extra_relations=(), + snub=False, + description="polytope-animation", +): + """ + Call POV-Ray to render the frames and call FFmpeg to generate the movie. + """ + if len(coxeter_diagram) == 3: + if snub: + P = models.Snub(coxeter_diagram, extra_relations=extra_relations) + else: + P = models.Polyhedra(coxeter_diagram, trunc_type, extra_relations) + scene_file = "polyhedra_animation.pov" + + elif len(coxeter_diagram) == 6: + P = models.Polychora(coxeter_diagram, trunc_type, extra_relations) + scene_file = "polytope_animation.pov" + + elif len(coxeter_diagram) == 10: + P = models.Polytope5D(coxeter_diagram, trunc_type, extra_relations) + scene_file = "polytope_animation.pov" + + else: + raise ValueError("Invalid Coxeter diagram: {}".format(coxeter_diagram)) + + P.build_geometry() + + # POV-Ray does not support 5d vectors well, so project the vertices in python + if isinstance(P, models.Polytope5D): + P.proj4d() + vert_data, edge_data, face_data = P.get_povray_data() + with open(data_file, "w") as f: + f.write(POV_TEMPLATE.format(vert_data, edge_data, face_data)) + + subprocess.call(POV_COMMAND.format(scene_file, description), shell=True) + subprocess.call(FFMPEG_COMMAND.format(description, description), shell=True) + return P + + +def snub24cell(description="snub-24-cell"): + """Handle the special case snub 24-cell. + """ + P = models.Snub24Cell() + P.build_geometry() + vert_data, edge_data, face_data = P.get_povray_data() + with open(data_file, "w") as f: + f.write(POV_TEMPLATE.format(vert_data, edge_data, face_data)) + scene_file = "polytope_animation.pov" + subprocess.call(POV_COMMAND.format(scene_file, description), shell=True) + subprocess.call(FFMPEG_COMMAND.format(description, description), shell=True) + return P + + +def main(): + # Platonic solids + anim((3, 2, 3), (1, 0, 0), description="tetrahedron") + anim((4, 2, 3), (1, 0, 0), description="cube") + anim((3, 2, 4), (1, 0, 0), description="octahedron") + anim((5, 2, 3), (1, 0, 0), description="dodecahedron") + anim((3, 2, 5), (1, 0, 0), description="icosahedron") + # Archimedian solids + anim((3, 2, 3), (1, 1, 0), description="truncated-tetrahedron") + anim((4, 2, 3), (1, 1, 0), description="truncated-cube") + anim((3, 2, 4), (1, 1, 0), description="truncated-octahedron") + anim((5, 2, 3), (1, 1, 0), description="truncated-dodecahedron") + anim((3, 2, 5), (1, 1, 0), description="truncated-icosahedron") + anim((4, 2, 3), (0, 1, 0), description="cuboctahedron") + anim((5, 2, 3), (0, 1, 0), description="icosidodecahedron") + anim((4, 2, 3), (1, 0, 1), description="rhombicuboctahedron") + anim((5, 2, 3), (1, 0, 1), description="rhombicosidodecahedron") + anim((4, 2, 3), (1, 1, 1), description="truncated-cuboctahedron") + anim((5, 2, 3), (1, 1, 1), description="truncated-icosidodecahedron") + anim((4, 2, 3), (1, 1, 1), snub=True, description="snub-cube") + anim((5, 2, 3), (1, 1, 1), snub=True, description="snub-dodecahedron") + # prism and antiprism + anim((7, 2, 2), (1, 0, 1), description="7-prism") + anim((8, 2, 2), (1, 1, 1), snub=True, description="8-antiprism") + # Kepler-Poinsot solids + anim( + (5, 2, Fraction(5, 2)), + (1, 0, 0), + extra_relations=((0, 1, 2, 1) * 3,), + description="great-dodecahedron", + ) + anim( + (5, 2, Fraction(5, 2)), + (0, 0, 1), + extra_relations=((0, 1, 2, 1) * 3,), + description="small-stellated-dodecahedron", + ) + anim((3, 2, Fraction(5, 2)), (0, 0, 1), description="great-stellated-dodecahedron") + anim((3, 2, Fraction(5, 2)), (1, 0, 0), description="great-icosahedron") + # 5-cell family, symmetry group A_4 + anim((3, 2, 2, 3, 2, 3), (1, 0, 0, 0), description="5-cell") + anim((3, 2, 2, 3, 2, 3), (1, 1, 0, 0), description="truncated-5-cell") + anim((3, 2, 2, 3, 2, 3), (0, 1, 0, 0), description="rectified-5-cell") + anim((3, 2, 2, 3, 2, 3), (0, 1, 1, 0), description="bitruncated-5-cell") + anim((3, 2, 2, 3, 2, 3), (1, 0, 0, 1), description="runcinated-5-cell") + + # tesseract family, symmetry group B_4 + anim((4, 2, 2, 3, 2, 3), (1, 0, 0, 0), description="tesseract") + anim((4, 2, 2, 3, 2, 3), (1, 1, 0, 0), description="truncated-tesseract") + anim((4, 2, 2, 3, 2, 3), (1, 0, 1, 0), description="cantellated-tesseract") + + # 16-cell family, dual to tesseract + anim((3, 2, 2, 3, 2, 4), (1, 0, 0, 0), description="16-cell") + anim((3, 2, 2, 3, 2, 4), (1, 1, 0, 0), description="truncated-16-cell") + anim((3, 2, 2, 3, 2, 4), (1, 0, 0, 1), description="runcinated-16-cell") + + # 24-cell family, symmetry group F_4 + anim((3, 2, 2, 4, 2, 3), (1, 0, 0, 0), description="24-cell") + anim((3, 2, 2, 4, 2, 3), (1, 0, 1, 0), description="cantellated-24-cell") + snub24cell() + # 120-cell family, symmetry group H_4 + anim((5, 2, 2, 3, 2, 3), (1, 0, 0, 0), description="120-cell") + + # 600-cell family, dual to 120-cell + P = anim((3, 2, 2, 3, 2, 5), (1, 0, 0, 0), description="600-cell") + P.draw_on_coxeter_plane(svgpath="600-cell.svg") + + # 4d prism and duoprism + anim((5, 2, 2, 3, 2, 2), (1, 1, 0, 1), description="truncated-dodecahedron-prism") + anim((3, 2, 2, 2, 2, 20), (1, 0, 0, 1), description="3-20-duoprism") + + # you can also render a 3d polyhedra by embedding it into 4d and project back. + anim((5, 2, 2, 3, 2, 2), (1, 1, 0, 0), description="truncated-dodecahedron") + anim( + (3, 2, 2, Fraction(5, 2), 2, 2), (1, 0, 0, 0), description="great-dodecahedron" + ) + + # some regular star polytopes (there are 10 of them, all can be rendered in this way) + anim( + (3, 2, 2, 5, 2, Fraction(5, 2)), + (1, 0, 0, 0), + extra_relations=((1, 2, 3, 2) * 3,), + description="icosahedral-120-cell", + ) + + anim( + (5, 2, 2, Fraction(5, 2), 2, 5), + (1, 0, 0, 0), + extra_relations=((0, 1, 2, 1) * 3, (1, 2, 3, 2) * 3), + description="great-120-cell", + ) + + anim( + (5, 2, 2, 3, 2, Fraction(5, 2)), + (1, 0, 0, 0), + extra_relations=((0, 1, 2, 3, 2, 1) * 3,), + description="grand-120-cell", + ) + + P = anim( + (Fraction(5, 2), 2, 2, 5, 2, Fraction(5, 2)), + (1, 0, 0, 0), + extra_relations=((0, 1, 2, 1) * 3, (1, 2, 3, 2) * 3), + description="grand-stellated-120-cell", + ) + P.draw_on_coxeter_plane(svgpath="grand-stellated-120-cell.svg") + + # and 5d polytopes + P = anim((4, 2, 2, 2, 3, 2, 2, 3, 2, 3), (1, 0, 0, 0, 0), description="5d-cube") + P.draw_on_coxeter_plane(svgpath="5-cube.svg") + + # some star polyhedron + anim( + (Fraction(3, 2), 3, 3), + (1, 0, 1), + description="octahemioctahedron", + extra_relations=((0, 1, 2, 1) * 2,), + ) + + anim( + (Fraction(3, 2), 5, 5), + (1, 0, 1), + description="small-dodecicosidodecahedron", + extra_relations=((0, 1, 2, 1) * 2,), + ) + + anim( + (Fraction(3, 2), 4, 4), + (1, 0, 1), + description="small-cubicuboctahedron", + extra_relations=((0, 1, 2, 1) * 2,), + ) + + anim( + (Fraction(5, 3), 3, 5), + (0, 1, 0), + description="ditrigonal-dodecadodecahedron", + extra_relations=((0, 1, 2, 1) * 2,), + ) + + +if __name__ == "__main__": + main() diff --git a/src/polytopes/example_run_coset_enumeration.py b/src/polytopes/example_run_coset_enumeration.py new file mode 100644 index 0000000..d9b63cb --- /dev/null +++ b/src/polytopes/example_run_coset_enumeration.py @@ -0,0 +1,174 @@ +""" +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Run coset enumeration examples using Todd-Coxeter algorithm +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This script reads information of a finitely presented group +from a .yaml file and computes its coset table. + +Example usage: + + python example_run_coset_enumeration.py filename [-std] [-o output] + + filename: required, the .yaml file to be parsed. + -std: optional, if added then the output table is standard. + -o: optional, output filename, the default is sys.stdout. + +:copyright (c) 2018 by Zhao Liang. +""" +import argparse +import re +import sys +import yaml + +from polytopes.todd_coxeter import CosetTable + + +def get_symbols(wordslist): + """Collect the set of letters from a list of strings. + """ + symbols = [] + for word in wordslist: + for x in word: + x = x.lower() + if not x.isalpha(): + raise ValueError("Only a-z and A-Z are allowed.") + if x not in symbols: + symbols.append(x) + return sorted(symbols) + + +def char2int(symbols, c): + """ + Find the integer in the generator list that represents a symbol `c`. + """ + ind = symbols.index(c.lower()) + return 2 * ind if c.islower() else 2 * ind + 1 + + +def word2int(symbols, wordslist): + """ + Map a list of words to the list of their integer representations. + """ + return tuple(tuple(char2int(symbols, c) for c in word) for word in wordslist) + + +def parse_latex(string): + """ + Given a generator relation as a LaTeX string, convert it to a flattened string. + Examples: a^{3} -> aaa, (ab)^2 -> abab, (aB)^{-2} -> bAbA + """ + if len(string) == 0: + return "" + if "^" not in string: + return string + + pattern = "(?:[a-zA-Z]{1}|\([a-zA-Z]+\))\^(?:\d+|\{-?\d+\})|[a-zA-Z]+(?!\^)" + result = re.findall(pattern, string) + word = "" + for substring in result: + if "^" not in substring: + word += substring + else: + base, exponent = substring.split("^") + if base[0] == "(": + base = base[1:-1] + if exponent[0] == "{": + exponent = exponent[1:-1] + + exponent = int(exponent) + if exponent < 0: + exponent = -exponent + base = base[::-1].swapcase() + + word += base * exponent + + return word + +class FpGroup(object): + + """ + Finitely presented group by defining relations between its generators. + """ + + def __init__(self, relators, subgens=None, name=None): + """ + :param relators: relations between the generators as a 1D list + of strings, e.g. ["aaa", "bb", "abab"]. + + :param subgens: generators of the subgroup as a 1D list of + strings, e.g. ["ab", "Ab"] + + :param name: decriptive name of the group, e.g. "S3", "D8", ..., etc. + """ + if not name: + name = self.__class__.__name__ + self.name = name + self.relators = [parse_latex(s) for s in relators] + self.generators = get_symbols(self.relators) + if len(self.generators) < 2: + raise ValueError("Not enough generators") + if not subgens: + subgens = [] + else: + subgens = [parse_latex(s) for s in subgens] + self.subgens = subgens + + # initialize the coset table + gens = tuple(range(2 * len(self.generators))) + relators += tuple(c + c.upper() for c in self.generators) + rels = word2int(self.generators, self.relators) + subgens = word2int(self.generators, self.subgens) + self.coset_table = CosetTable(gens, rels, subgens, coxeter=False) + + def __str__(self): + s = "\nName: {}\n".format(self.name) + s += "Generators: " + ", ".join(self.generators) + s += "\nRelators: " + ", ".join(self.relators) + s += "\nSubgroup generators: " + ", ".join(self.subgens) + return s + + def compute(self, standard=True): + self.coset_table.run(standard) + + def print_table(self, outfile): + """pretty print the table. + """ + f = sys.stdout if outfile is None else open(outfile, "w") + f.write(" ") + for x in self.generators: + f.write("{:>5s}".format(x)) + f.write("{:>5s}".format(x.upper())) + f.write("\n" + (2 * len(self.generators) + 2) * 5 * "-" + "\n") + + for i, row in enumerate(self.coset_table, start=1): + f.write("{:>5d}: ".format(i)) + for x in row: + f.write("{:>5d}".format(x + 1)) + f.write("\n") + f.close() + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("filename", type=str, help="Input file name") + parser.add_argument( + "-std", type=bool, default=True, help="Standardize the coset table or not" + ) + parser.add_argument( + "-out", metavar="-o", type=str, default=None, help="output file name" + ) + args = parser.parse_args() + + with open(args.filename, "r") as f: + data = yaml.safe_load(f) + rels = data["relators"] + subg = data["subgroup-generators"] + name = data["name"] + G = FpGroup(rels, subg, name) + G.compute(args.std) + G.print_table(args.out) + + +if __name__ == "__main__": + main() diff --git a/src/polytopes/polytopes/__init__.py b/src/polytopes/polytopes/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/polytopes/polytopes/coxeter_plane.py b/src/polytopes/polytopes/coxeter_plane.py new file mode 100644 index 0000000..c729f86 --- /dev/null +++ b/src/polytopes/polytopes/coxeter_plane.py @@ -0,0 +1,73 @@ +""" +Project the vertices of a polytope to the Coxeter plane. +""" +try: + import cairocffi as cairo +except ImportError: + import cairo + +import numpy as np + +from . import helpers + + +def draw_on_coxeter_plane( + P, svgpath=None, image_size=600, linewidth=0.0012, markersize=0.015 +): + """ + Project the vertices of a polytope `P` to its Coxeter plane + and draw the pattern to a svg image. + + The most important parameters are `nodes1` and `nodes2`, they + can be of lists/tuples/sets type and must partition the Coxeter + diagram of `P` into two disjoint sets. The nodes in each set are + mutually orthogonal with each other. + """ + P.build_geometry() + M = P.mirrors + C = 2 * np.dot(M, M.T) # Cartan matrix + eigenvals, eigenvecs = np.linalg.eigh(C) + # get the eigenvector with largest (or smallest) eigenvalue + u = eigenvecs[:, 0] + v = eigenvecs[:, -1] + # a basis of the Coxeter plane + u = np.dot(u, M) + v = np.dot(v, M) + u = helpers.normalize(u) + v = helpers.normalize(v) + vertices_2d = [(np.dot(u, x), np.dot(v, x)) for x in P.vertices_coords] + + # draw on image + if svgpath is None: + svgpath = P.__class__.__name__ + ".svg" + + extent = 0.99 + surface = cairo.SVGSurface(svgpath, image_size, image_size) + ctx = cairo.Context(surface) + ctx.scale(image_size / (extent * 2.0), -image_size / (extent * 2.0)) + ctx.translate(extent, -extent) + ctx.set_source_rgb(1, 1, 1) + ctx.paint() + ctx.set_line_width(linewidth) + + # draw edges + for elist in P.edge_indices: + for i, j in elist: + x1, y1 = vertices_2d[i] + x2, y2 = vertices_2d[j] + ctx.set_source_rgb(0, 0, 0) + ctx.set_line_width(linewidth) + ctx.move_to(x1, y1) + ctx.line_to(x2, y2) + ctx.stroke() + + # draw the vertices as circles + ctx.set_line_width(linewidth * 2) + for x, y in vertices_2d: + ctx.arc(x, y, markersize, 0, 2 * np.pi) + ctx.set_source_rgb(1, 0, 0) + ctx.fill_preserve() + ctx.set_source_rgb(0, 0, 0) + ctx.stroke() + + surface.finish() diff --git a/src/polytopes/polytopes/helpers.py b/src/polytopes/polytopes/helpers.py new file mode 100644 index 0000000..b81e0ac --- /dev/null +++ b/src/polytopes/polytopes/helpers.py @@ -0,0 +1,123 @@ +""" +Some helper functions for building the geometry. +""" +import sys +import numpy as np + + +def normalize(v): + """Normalize a vector `v`. + """ + return np.array(v) / np.linalg.norm(v) + + +def reflection_matrix(v): + """ + Return the reflection transformation about a plane with normal vector `v`. + see "https://en.wikipedia.org/wiki/Householder_transformation". + """ + n = len(v) + v = np.array(v)[np.newaxis] + return np.eye(n) - 2 * np.dot(v.T, v) + + +def get_init_point(M, d): + """ + Given the normal vectors of the mirrors stored as row vectors in `M` + and a tuple of non-negative floats `d`, compute the vector `v` whose + distance vector to the mirrors is `d` and return its normalized version. + """ + return normalize(np.linalg.solve(M, d)) + + +def proj3d(v): + """Stereographic projection from 4d to 3d. + """ + v = normalize(v) + return v[:-1] / (1 + 1e-8 - v[-1]) # avoid divide by zero + + +def make_symmetry_matrix(upper_triangle): + """Make a symmetric matrix from a list of integers/rationals + as the upper half above the diagonals. + + Example: + >>> upper_trianle = [4, 2, 3] + >>> make_symmstry_matrix(upper_triangle) + >>> [[1, 4, 2], + [4, 1, 3], + [2, 3, 1]] + """ + n = int((2 * len(upper_triangle)) ** 0.5) + 1 + if len(upper_triangle) != (n * n - n) / 2: + raise ValueError("Invalid input sequence") + ind = 0 + cox_mat = np.eye(n, dtype="object") + for row in range(n - 1): + for col in range(row + 1, n): + cox_mat[row][col] = cox_mat[col][row] = upper_triangle[ind] + ind += 1 + return cox_mat + + +def get_coxeter_matrix(coxeter_diagram): + """ + Get the Coxeter matrix from a given coxeter_diagram. + The Coxeter matrix is square and entries are all integers, it describes + the relations between the generators of the symmetry group. Here is the + math: suppose two mirrors mᵢ, mⱼ form an angle p/q where p,q are coprime + integers, then the two generator reflections about these two mirrors rᵢ, rⱼ + satisfy (rᵢrⱼ)^p = 1. + + Example: + >>> coxeter_diagram = (3, 2, Fraction(5, 2)) + >>> get_coxeter_matrix(coxeter_diagram) + >>> [[1, 3, 2], + [3, 1, 5], + [2, 5, 1]] + + Note that in general one cannot recover the Coxeter diagram + from the Coxeter matrix since a star polytope may have the + same Coxeter matrix with a convex one. + """ + upper_triangle = [x.numerator for x in coxeter_diagram] + return make_symmetry_matrix(upper_triangle) + + +def get_mirrors(coxeter_diagram): + """ + Given a Coxter diagram consists of integers/rationals that represent + the angles between the mirrors (a rational p means the angle is π/p), + return a square matrix whose rows are the normal vectors of the mirrors. + This matrix is not unique, here we use a lower triangle one to simplify + the computations. + """ + # error handling function when the input coxeter matrix is invalid. + def err_handler(err_type, flag): + print( + "Invalid input Coxeter diagram. This diagram does not give a finite \ +symmetry group of an uniform polytope. See \ +https://en.wikipedia.org/wiki/Coxeter_group#Symmetry_groups_of_regular_polytopes \ +for a complete list of valid Coxeter diagrams." + ) + sys.exit(1) + + np.seterrcall(err_handler) + np.seterr(all="call") + + coxeter_matrix = np.array( + make_symmetry_matrix(coxeter_diagram) + ).astype(float) + C = -np.cos(np.pi / coxeter_matrix) + M = np.zeros_like(C) + n = len(M) + # the first normal vector is simply (1, 0, ...) + M[0, 0] = 1 + # in the i-th row, the j-th entry can be computed via the (j, j) entry. + for i in range(1, n): + for j in range(i): + M[i, j] = (C[i, j] - np.dot(M[j, :j], M[i, :j])) / M[j, j] + # the (i, i) entry is used to normalize this vector + M[i, i] = np.sqrt(1 - np.dot(M[i, :i], M[i, :i])) + + return M diff --git a/src/polytopes/polytopes/models.py b/src/polytopes/polytopes/models.py new file mode 100644 index 0000000..f39ba86 --- /dev/null +++ b/src/polytopes/polytopes/models.py @@ -0,0 +1,671 @@ +""" +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Classes for building models of 3D/4D/5D uniform polytopes +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +See the doc: "https://neozhaoliang.github.io/polytopes/" + +""" +from itertools import combinations +import numpy as np + +from . import helpers +from .coxeter_plane import draw_on_coxeter_plane +from .povray import export_polytope_data +from .todd_coxeter import CosetTable + + +class BasePolytope: + + """ + Base class for building uniform polytopes using Wythoff's construction. + """ + + def __init__(self, coxeter_diagram, init_dist, extra_relations=()): + """ + :param coxeter_diagram: a tuple of rational numbers. + Coxeter diagram for this polytope. + + :param init_dist: a tuple of non-negative floats. + distances between the initial vertex and the mirrors. + + :param extra_relations: a tuple of tuples of integers. + a presentation of a star polytope can be obtained by imposing more + relations on the generators. For example "(ρ₀ρ₁ρ₂ρ₁)^n = 1" for some + integer n, where n is the number of sides of a hole. + See Coxeter's article + + "Regular skew polyhedra in three and four dimensions, + and their topological analogues" + + """ + # Coxeter matrix of the symmetry group + self.coxeter_matrix = helpers.get_coxeter_matrix(coxeter_diagram) + + # reflection mirrors stored as row vectors in a matrix + self.mirrors = helpers.get_mirrors(coxeter_diagram) + + # reflection transformations about the mirrors + self.reflections = tuple(helpers.reflection_matrix(v) for v in self.mirrors) + + # the initial vertex + self.init_v = helpers.get_init_point(self.mirrors, init_dist) + + # a mirror is active if and only if the initial vertex has non-zero distance to it + self.active = tuple(bool(x) for x in init_dist) + + # generators of the symmetry group + self.symmetry_gens = tuple(range(len(self.coxeter_matrix))) + + # relations between the generators + self.symmetry_rels = tuple( + (i, j) * self.coxeter_matrix[i][j] + for i, j in combinations(self.symmetry_gens, 2) + ) + # add the extra relations between generators (only used for star polytopes) + self.symmetry_rels += tuple(extra_relations) + + # multiplication table bwteween the vertices + self.vtable = None + # word representations of the vertices + self.vwords = None + + # number of vertices, edges, faces + self.num_vertices = 0 + self.num_edges = 0 + self.num_faces = 0 + + # coordinates of the vertices, indices of the edges and faces + self.vertices_coords = [] + self.edge_indices = [] + self.face_indices = [] + + def build_geometry(self): + """Put the building procedure into three functions. + """ + self.get_vertices() + self.get_edges() + self.get_faces() + + def get_vertices(self): + """ + This method computes the following data that will be used later: + + 1. a coset table for the vertices. + 2. a complete list of word representations of the vertices. + 3. coordinates of the vertices. + """ + # generators of the stabilizing subgroup that fixes the initial vertex. + vgens = [(i,) for i, active in enumerate(self.active) if not active] + self.vtable = CosetTable(self.symmetry_gens, self.symmetry_rels, vgens) + self.vtable.run() + # word representations of the vertices + self.vwords = self.vtable.get_words() + self.num_vertices = len(self.vwords) + # apply words of the vertices to the initial vertex to get all vertices + self.vertices_coords = np.array([ + self.transform(self.init_v, w) for w in self.vwords + ]) + + def get_edges(self): + """ + Compute the indices of all edges. + + 1. If the initial vertex v₀ lies on the mirror mᵢ then the reflection + ρᵢ fixes v₀, hence there are no edges of type i. + + 2. Otherwise, v₀ and its virtual image v₁ about mᵢ generate a base edge + e₀, the stabilizing subgroup of e₀ is generated by ρᵢ together those + simple reflections ρⱼ such that ρⱼ fixes v₀ and commutes with ρᵢ: ρᵢρⱼ=ρⱼρᵢ. + Again we use Todd-Coxeter's procedure to get a complete list of word + representations for the edges of type i and apply them to e₀ to get all + edges of type i. + """ + for i, active in enumerate(self.active): + # if there are edges of type i + if active: + # the initial edge + e0 = (0, self.move(0, (i,))) + # generators of the edge stabilizing subgroup + egens = [(i,)] + self.get_orthogonal_stabilizing_mirrors((i,)) + # get word representations of the edges + coset_reps = self.get_coset_representatives(egens) + # apply them to the base edge to get all edges of type i + orbit = self.get_orbit(coset_reps, e0) + self.edge_indices.append(orbit) + self.num_edges += len(orbit) + + def get_faces(self): + """ + Compute the indices of all faces. + The composition of two reflections ρᵢ and ρⱼ is a rotation that + fixes a base face f₀. The stabilizing subgroup of f₀ is generated + by {ρᵢ, ρⱼ} together with those simple reflections ρₖ such that ρₖ + fixes v₀ and commutes with both ρᵢ and ρⱼ. + """ + for i, j in combinations(self.symmetry_gens, 2): + m = self.coxeter_matrix[i][j] + f0 = [] + # if both two mirrors are active then they generate a face + if self.active[i] and self.active[j]: + for k in range(m): + f0.append(self.move(0, (i, j) * k)) + f0.append(self.move(0, (j,) + (i, j) * k)) + # if exactly one of the two mirrors are active then they + # generate a face only when they are not perpendicular + elif (self.active[i] or self.active[j]) and m > 2: + for k in range(m): + f0.append(self.move(0, (i, j) * k)) + # else they do not generate a face + else: + continue + + # generators of the face stabilizing subgroup + fgens = [(i,), (j,)] + self.get_orthogonal_stabilizing_mirrors([i, j]) + coset_reps = self.get_coset_representatives(fgens) + # apply the coset representatives to the base face + orbit = self.get_orbit(coset_reps, f0) + self.face_indices.append(orbit) + self.num_faces += len(orbit) + + def transform(self, vector, word): + """ + Transform a vector by a word in the symmetry group. + Return the coordinates of the resulting vector. + + :param vector: a 1d array, e.g. (1, 0, 0). + :param word: a list of integers. + """ + for w in word: + vector = np.dot(vector, self.reflections[w]) + return vector + + def move(self, vertex, word): + """ + Transform a vertex by a word in the symmetry group. + Return the index of the resulting vertex. + + :param vertex: an integer. + :param word: a list of integers. + """ + for w in word: + vertex = self.vtable[vertex][w] + return vertex + + def get_orthogonal_stabilizing_mirrors(self, subgens): + """ + Given a list of generators in `subgens`, return the generators that + commute with all of those in `subgens` and fix the initial vertex. + + :param subgens: a list of generators, e.g. [0, 1] + """ + result = [] + for s in self.symmetry_gens: + # check commutativity + if all(self.coxeter_matrix[x][s] == 2 for x in subgens): + # check if it fixes v0 + if not self.active[s]: + result.append((s,)) + return result + + def get_coset_representatives(self, subgens, coxeter=True): + """ + Get the set of coset representatives of a subgroup generated + by `subgens`. + + :param subgens: + a list of generating words of the subgroup, e.g. [(0,), (1, 2)] + """ + table = CosetTable(self.symmetry_gens, self.symmetry_rels, subgens, coxeter) + table.run() + return table.get_words() + + def get_orbit(self, coset_reps, base): + """ + Apply the words in `coset_reps` to a base edge/face. + + :param coset_reps: a list of words. + :param base: a 1d list of integers. + """ + return [[self.move(v, word) for v in base] for word in coset_reps] + + def get_latex_format(self, symbol=r"\rho", cols=3, snub=False): + """ + Return the words of the vertices in latex format. + `cols` is the number of columns of the output latex array. + """ + + def to_latex(word): + if not word: + return "e" + if snub: + return "".join(symbol + "_{{{}}}".format(i // 2) for i in word) + return "".join(symbol + "_{{{}}}".format(i) for i in word) + + latex = "" + for i, word in enumerate(self.vwords): + if i > 0 and i % cols == 0: + latex += r"\\" + latex += to_latex(word) + if i % cols != cols - 1: + latex += "&" + + return r"\begin{{array}}{{{}}}{}\end{{array}}".format("l" * cols, latex) + + def get_povray_data(self): + return export_polytope_data(self) + + def draw_on_coxeter_plane(self, *args, **kwargs): + draw_on_coxeter_plane(self, *args, **kwargs) + + +class Polyhedra(BasePolytope): + + """Base class for 3d polyhedron. + """ + + def __init__(self, coxeter_diagram, init_dist, extra_relations=()): + if not len(coxeter_diagram) == len(init_dist) == 3: + raise ValueError("Length error: the inputs must all have length 3") + super().__init__(coxeter_diagram, init_dist, extra_relations) + + +class Snub(Polyhedra): + + """ + A snub polyhedra is generated by the subgroup that consists of only + rotations in the full symmetry group. This subgroup has presentation + + + + where r = ρ₀ρ₁, s = ρ₁ρ₂ are two rotations. + Again we solve all words in this subgroup and then use them to + transform the initial vertex to get all vertices. + """ + + def __init__( + self, + coxeter_diagram, + init_dist=(1.0, 1.0, 1.0), + extra_relations=() + ): + super().__init__(coxeter_diagram, init_dist, extra_relations) + # the representation is not in the form of a Coxeter group, + # we must overwrite the relations. + + # four generators (r, r^-1, s, s^-1) + self.symmetry_gens = (0, 1, 2, 3) + + # relations in order: + # 1. r^p = 1 + # 2. s^q = 1 + # 3. (rs)^2 = 1 + # 4. rr^-1 = 1 + # 5. ss^-1 = 1 + self.symmetry_rels = ( + (0,) * self.coxeter_matrix[0][1], + (2,) * self.coxeter_matrix[1][2], + (0, 2) * self.coxeter_matrix[0][2], + (0, 1), + (2, 3), + ) + # map the extra_relations expressed by reflections + # into relations by (r, s). + for extra_rel in extra_relations: + if len(extra_rel) % 2 == 1: + extra_rel += extra_rel + + snub_rel = [] + for x, y in zip(extra_rel[:-1:2], extra_rel[1::2]): + snub_rel.extend( + { + (0, 1): [0], + (0, 2): [0, 2], + (1, 0): [1], + (1, 2): [2], + (2, 0): [2, 0], + (2, 1): [3], + }[(x, y)] + ) + self.symmetry_rels += (tuple(snub_rel),) + + # order of the generator rotations {rotation: order} + # {r: p, s: q, rs: 2} + self.rotations = { + (0,): self.coxeter_matrix[0][1], + (2,): self.coxeter_matrix[1][2], + (0, 2): self.coxeter_matrix[0][2], + } + + def get_vertices(self): + """Get the vertices of this snub polyhedra. + """ + # the stabilizing subgroup of the initial vertex contains only 1 + self.vtable = CosetTable( + self.symmetry_gens, self.symmetry_rels, coxeter=False + ) + self.vtable.run() + self.vwords = self.vtable.get_words() + self.num_vertices = len(self.vwords) + self.vertices_coords = tuple( + self.transform(self.init_v, w) for w in self.vwords + ) + + def get_edges(self): + """ + Get the edge indices of this snub polyhedra. + Each rotation of the three "fundamental rotations" {r, s, rs} + generates a base edge e, the stabilizing subgroup of e is <1> + except that this rotation has order 2, i.e. a rotation generated + by two commuting reflections. In this case the stabilizing subgroup + is the cyclic group of order 2. + """ + for rot in self.rotations: + # if this rotation has order 2, then the edge stabilizing subgroup + # is the cyclic subgroup + if self.rotations[rot] == 2: + egens = (rot,) + coset_reps = self.get_coset_representatives(egens, coxeter=False) + # else the edge stabilizing subgroup is <1> + else: + coset_reps = self.vwords + + e0 = (0, self.move(0, rot)) # the base edge + # apply coset representatives to e0 to get all edges + orbit = self.get_orbit(coset_reps, e0) + self.edge_indices.append(orbit) + self.num_edges += len(orbit) + + def get_faces(self): + """ + Get the face indices of this snub polyhedra. + Each rotation of the three "fundamental rotations" {r, s, rs} + generates a base face f₀ if the order of this rotation is strictly + greater than two, else it only gives an edge (degenerated face). + + There's another type of face given by the relation r*s = rs: + it's generated by the three vertices {v₀, v₀s, v₀rs}. + Note (v₀, v₀s) is an edge of type s, (v₀, v₀rs) is an edge of type rs, + (v₀s, v₀rs) is an edge of type r since it's in the same orbit of the + edge (v₀, v₀r) by applying s on it. + """ + for rot, order in self.rotations.items(): + # if the order of this rotation is > 2 then it generates a face + if order > 2: + f0 = tuple(self.move(0, rot * k) for k in range(order)) + # the stabilizing group is the cyclic group + fgens = (rot,) + coset_reps = self.get_coset_representatives(fgens, coxeter=False) + orbit = self.get_orbit(coset_reps, f0) + self.face_indices.append(orbit) + self.num_faces += len(orbit) + + # handle the special triangle face (v0, v0s, v0rs) + # note its three edges are in different orbits so + # its stabilizing subgroup must be <1>. + triangle = (0, self.move(0, (2,)), self.move(0, (0, 2))) + orbit = self.get_orbit(self.vwords, triangle) + self.face_indices.append(orbit) + self.num_faces += len(orbit) + + def transform(self, vector, word): + """ + Transform a vector by a word in the group. + Return the coordinates of the resulting vector. + Note generator 0 means r = ρ₀ρ₁, generator 1 means s = ρ₁ρ₂. + """ + for g in word: + if g == 0: + vector = np.dot(vector, self.reflections[0]) + vector = np.dot(vector, self.reflections[1]) + else: + vector = np.dot(vector, self.reflections[1]) + vector = np.dot(vector, self.reflections[2]) + return vector + + +class Polychora(BasePolytope): + + """Base class for 4d polychoron. + """ + + def __init__(self, coxeter_diagram, init_dist, extra_relations=()): + if not (len(coxeter_diagram) == 6 and len(init_dist) == 4): + raise ValueError( + "Length error: the input coxeter_diagram must have length 6 and init_dist has length 4" + ) + super().__init__(coxeter_diagram, init_dist, extra_relations) + + +class Polytope5D(BasePolytope): + + """Base class for 5d uniform polytopes. + """ + + def __init__(self, coxeter_diagram, init_dist, extra_relations=()): + if len(coxeter_diagram) != 10 and len(init_dist) != 5: + raise ValueError( + "Length error: the input coxeter_diagram must have length 10 and init_dist has length 5" + ) + super().__init__(coxeter_diagram, init_dist, extra_relations) + + def proj4d(self, pole=1.3): + """Stereographic project vertices to 4d. + """ + self.vertices_coords = [v[:4] / (pole - v[-1]) for v in self.vertices_coords] + return self + + +class Snub24Cell(Polychora): + + """ + The snub 24-cell can be constructed from snub demitesseract [3^(1,1,1)]+, + the procedure is similar with snub polyhedron above. + Coxeter-Dynkin diagram: + + ρ₀ ρ₁ ρ₂ + •-----•-----• + | + • + ρ₃ + + Its symmetric group is generated by three rotations {r, s, t}, where r = ρ₀ρ₁, + s = ρ₁ρ₂, t = ρ₁ρ₃. A presentation is + + G = + + """ + + def __init__(self): + coxeter_diagram = (3, 2, 2, 3, 3, 2) + active = (1, 1, 1, 1) + super().__init__(coxeter_diagram, active, extra_relations=()) + # generators in order: {r, r^-1, s, s^-1, t, t^-1} + self.symmetry_gens = tuple(range(6)) + # relations in order: + # 1. r^3 = 1 + # 2. s^3 = 1 + # 3. t^3 = 1 + # 4. (rs)^2 = 1 + # 5. (rt)^2 = 1 + # 6. (s^-1t)^2 = 1 + # 7. rr^-1 = 1 + # 8. ss^-1 = 1 + # 9. tt^-1 = 1 + self.symmetry_rels = ( + (0,) * 3, + (2,) * 3, + (4,) * 3, + (0, 2) * 2, + (0, 4) * 2, + (3, 4) * 2, + (0, 1), + (2, 3), + (4, 5), + ) + # rotations and their order + # {r: 3, s: 3, t: 3, rs: 2, rt: 2, s^-1t: 2} + self.rotations = {(0,): 3, (2,): 3, (4,): 3, (0, 2): 2, (0, 4): 2, (3, 4): 2} + + def get_vertices(self): + """Get the coordinates of the snub 24-cell. + """ + self.vtable = CosetTable( + self.symmetry_gens, self.symmetry_rels, coxeter=False + ) + self.vtable.run() + self.vwords = self.vtable.get_words() + self.num_vertices = len(self.vwords) + self.vertices_coords = tuple( + self.transform(self.init_v, w) for w in self.vwords + ) + + def get_edges(self): + """Get the edges of the snub 24-cell. Again each fundamental rotation + in {r, s, t, rs, rt, s^-1t} generates edges of its type. + """ + for rot, order in self.rotations.items(): + # the initial edge + e0 = (0, self.move(0, rot)) + # if the rotation has order 2 then the stabilizing subgroup of + # the init edge is the cyclic subgroup , else it's <1>. + if order == 2: + egens = (rot,) + coset_reps = self.get_coset_representatives(egens, coxeter=False) + else: + coset_reps = self.vwords + + orbit = self.get_orbit(coset_reps, e0) + self.edge_indices.append(orbit) + self.num_edges += len(orbit) + + def get_faces(self): + """ + Get the faces of the snub 24-cell. + + 1. A fundamental rotation generates a face by rotating the initial vertex + k times. This face is non-degenerate if and only if the order of the + rotation is strictly greater than two. Only {r, s, t} can give such faces. + + 2. Three funtamental rotations {r, s, t} satifying rs = t also could generate + some triangle faces. + """ + for rot in ((0,), (2,), (4,)): + order = self.rotations[rot] + f0 = tuple(self.move(0, rot * k) for k in range(order)) + fgens = (rot,) + coset_reps = self.get_coset_representatives(fgens, coxeter=False) + orbit = self.get_orbit(coset_reps, f0) + self.face_indices.append(orbit) + self.num_faces += len(orbit) + + # handle the special triangle faces generated from + # 1. {v0, v0s, v0rs} + # 2. {v0, v0t, v0rt}, + # 3. {v0, v0s, v0t^-1s} + # 4. {v0, v0rs, v0t^-1s} + # the edges of these triangles are in different orbits + # hence their stabilizing subgroups are all <1>. + for triangle in [ + (0, self.move(0, (2,)), self.move(0, (0, 2))), + (0, self.move(0, (4,)), self.move(0, (0, 4))), + (0, self.move(0, (2,)), self.move(0, (5, 2))), + (0, self.move(0, (0, 2)), self.move(0, (5, 2))), + ]: + orbit = self.get_orbit(self.vwords, triangle) + self.face_indices.append(orbit) + self.num_faces += len(orbit) + + def transform(self, vector, word): + """ + The generators are 0 for r=ρ₀ρ₁, 2 for s=ρ₁ρ₂, 4 for t=ρ₁ρ₃. + """ + for g in word: + if g == 0: + vector = np.dot(vector, self.reflections[0]) + vector = np.dot(vector, self.reflections[1]) + elif g == 2: + vector = np.dot(vector, self.reflections[1]) + vector = np.dot(vector, self.reflections[2]) + else: + vector = np.dot(vector, self.reflections[1]) + vector = np.dot(vector, self.reflections[3]) + return vector + + +class Catalan3D(Polyhedra): + + """Catalan solids are duals of uniform polyhedron. + These polyhedron are face transitive but (usually) not vertex transitive. + But to keep things consistent we still put the vertices in a 1d list. + """ + + def __init__(self, P): + """Construct a Catalan solid form a given polyhedra `P`. + """ + self.P = P + self.vertices_coords = [] + self.face_indices = [] + + def build_geometry(self): + self.P.build_geometry() + self.get_vertices() + self.get_faces() + + def get_vertices(self): + """Each vertex in this dual polyhedra comes from a face f in the + original polyhedra `P`. Usually it's not the center of f but a + scaled version of it. + """ + for face_group in self.P.face_indices: + for face in face_group: + verts = [self.P.vertices_coords[ind] for ind in face] + cen = sum(verts) + normal = helpers.normalize(cen) + weights = sum(np.dot(v, normal) for v in verts) / len(face) + self.vertices_coords.append(normal / weights) + + def get_faces(self): + """Each face f in this dual polyhedra comes from a vertex v in the + original polyhedra `P`. The vertices in f are faces of `P` that meet at v. + """ + def contain_edge(f, e): + """Check if a face f contains a given edge e = (v1, v2). + """ + v1, v2 = e + for w1, w2 in zip(f, f[1:] + [f[0]]): + if (v1, v2) in [(w1, w2), (w2, w1)]: + return True + return False + + def is_adjacent(f1, f2): + """Check if two faces f1, f2 are adjacent. + """ + for e in zip(f1, f1[1:] + [f1[0]]): + if contain_edge(f2, e): + return True + return False + + result = [] + P_faces_flatten = [face for face_group in self.P.face_indices + for face in face_group] + # there is a face for each vertex v in the original polyhedra P + for k in range(len(self.P.vertices_coords)): + # firstly we gather all faces in P that meet at v + faces_unordered = [] + for ind, f in enumerate(P_faces_flatten): + if k in f: + faces_unordered.append([ind, f]) + + # then we re-align them so that they form a cycle around v + i0, f0 = faces_unordered[0] + face = [i0] + current_face = f0 + while len(face) < len(faces_unordered): + for ind, f in faces_unordered: + if ind not in face and is_adjacent(current_face, f): + face.append(ind) + current_face = f + + result.append(face) + self.face_indices.append(result) diff --git a/src/polytopes/polytopes/povray.py b/src/polytopes/polytopes/povray.py new file mode 100644 index 0000000..8018c16 --- /dev/null +++ b/src/polytopes/polytopes/povray.py @@ -0,0 +1,62 @@ +""" +Some POV-Ray stuff. +""" + + +def concat(arr, sep=",\n", border=None): + """Concatenate items in an array using seperator `sep`. + """ + if border is None: + return sep.join(str(x) for x in arr) + return border.format(sep.join(str(x) for x in arr)) + + +def pov_vector(v): + """Convert a vector to POV-Ray format. e.g. (x, y, z) --> . + """ + return concat(v, sep=", ", border="<{}>") + + +def pov_vector_array(vectors): + """Convert a list of vectors to POV-Ray format, e.g. + [(x, y, z), (a, b, c), ...] --> array[n] {, , ...}. + """ + n = len(vectors) + return "array[{}]{{\n{}}}".format(n, concat(pov_vector(v) for v in vectors)) + + +def pov_index_array1d(arr1d): + """Convert a 1d array to POV-Ray string, e.g. (1, 2, 3) --> array[3] {1, 2, 3} + """ + n = len(arr1d) + return "array[{}] {{\n{}}}".format(n, concat(arr1d, sep=", ")) + + +def pov_index_array2d(arr2d): + """ + Convert a mxn 2d array to POV-Ray format, e.g. + [(1, 2), (3, 4), (5, 6)] --> arrar[3][2] {{1, 2}, {3, 4}, {5, 6}}. + """ + dim1 = len(arr2d) + dim2 = len(arr2d[0]) + return "array[{}][{}] {{\n{}}}".format( + dim1, dim2, concat(concat(arr, sep=", ", border="{{{}}}") for arr in arr2d) + ) + + +def pov_index_array3d(arr3d): + """Convert a 3d array to an array of 2d arrays in POV-Ray format. + """ + n = len(arr3d) + return "array[{}] {{\n{}}}".format( + n, concat(pov_index_array2d(arr2d) for arr2d in arr3d) + ) + + +def export_polytope_data(P): + """Return the data of a polytope `P` in POV-Ray format. + """ + vert_data = pov_vector_array(P.vertices_coords) + edge_data = pov_index_array3d(P.edge_indices) + face_data = pov_index_array3d(P.face_indices) + return vert_data, edge_data, face_data \ No newline at end of file diff --git a/src/polytopes/polytopes/todd_coxeter.py b/src/polytopes/polytopes/todd_coxeter.py new file mode 100644 index 0000000..8a86299 --- /dev/null +++ b/src/polytopes/polytopes/todd_coxeter.py @@ -0,0 +1,317 @@ +""" +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Coset Enumeration using Todd-Coxeter Algorithm +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Reference: + [1] "Handbook of Computational Group Theory", Holt, D., Eick, B., O'Brien, E. + [2] "Computation with finitely presented groups", Charles C.Sims. + [3] GAP doc at "https://www.gap-system.org/Manuals/doc/ref/chap47.html". + [4] Ken Brown's code at "http://www.math.cornell.edu/~kbrown/toddcox/". + +""" + + +class CosetTable: + """ + Let G = < X | R > be a finitely presented group and H be a subgroup of + G with |G:H| < inf. The coset table T of G/H is a 2D array with rows + indexed by the right cosets of G/H and columns indexed by the generators + of G and their inverses. The entry of T at row k and column x (x is a + generator or the inverse of a generator) records the right action of + coset k by x: T[k][x] = kx. If T[k][x] is not defined yet we set it to + None. When the algorithm terminates all entries in the table are + assigned a non-negative integer, all rows scan correctly under all words + in R and the first row scans correctly under all generators of H. + """ + + def __init__(self, gens, rels, subgens=(), coxeter=True): + """ + :param gens: a 1D list of integers that represents the generators, + e.g. [0, 1, 2] for [a, b, c] + + :param rels: a 2D list of integers that represents the relators R, + e.g. [[0, 0], [1, 1, 1], [0, 2]*3] for a^2 = b^3 = (ac)^3 = 1 + + :param subgens: a 2D list of integers that represents the subgroup generators, + e.g. [[0, 2]] for < ac >. + + :param coxeter: whether this is a Coxeter group. If false then the inverse of + the generators are also included as generators, otherwise they are not. + + we use a list p to hold the equivalence classes of the cosets, + p[k] = l means k and l really represent the same coset. It's always + true that p[k] <= k, if p[k] = k then we call k "alive" otherwise we call + k "dead". All cosets are created "alive", but as the algorithm runs + some of them may be declared to "dead" and their rows are deleted from + the table (we do not explicitly free the memory of these row but keep + in mind that they do not belong to our table). + + A "dead" coset arise when an `coincidence` is found, and while handling + this coincidence more coincidences may also be found, so we use a queue + q to hold them. + """ + self.coxeter = coxeter + if coxeter: + self.inv = lambda x: x + else: + self.inv = lambda x: x - 1 if x % 2 else x + 1 + self.A = gens + self.R = rels # relations R between the generators + self.H = subgens # generators of H + self.p = [0] # initially we only have the 0-th coset H + self.q = [] # a queue holds all dead cosets to be processed. + self.table = [[None] * len(self.A)] + + def __getitem__(self, item): + return self.table.__getitem__(item) + + def __len__(self): + return len(self.table) + + def is_alive(self, coset): + """Check if a coset is alive. + """ + return self.p[coset] == coset + + def is_defined(self, coset, x): + """Check if an entry is defined. + """ + return self[coset][x] is not None + + def undefine(self, coset, x): + self[coset][x] = None + + def define(self, coset, x): + """ + To define a new coset we need to: + 1. Append a new empty row to the table. + 2. Fill the two entries. + 3. Add this new coset to `p`. + """ + n = len(self.table) + self.table.append([None] * len(self.A)) + self[coset][x] = n + self[n][self.inv(x)] = coset + self.p.append(n) + + def rep(self, coset): + """Find the minimal equivalent representative of a given coset + and modify `p` along the way. + """ + m = coset + while m != self.p[m]: + m = self.p[m] + + j = coset + while j != self.p[j]: + j, self.p[j] = self.p[j], m + + return m + + def merge(self, coset1, coset2): + """ + Merge two equivalent cosets. The larger one is declared to "dead" + and is added to the queue `q` to be processed later on. + """ + s = self.rep(coset1) + t = self.rep(coset2) + if s != t: + s, t = min(s, t), max(s, t) + self.p[t] = s + self.q.append(t) + + def scan_and_fill(self, coset, word): + """ + Scan the row of a coset under a given word. + 1. Firstly we start from the left and scan forward to the right as + far as possible. + (1a) If it completes correctly then it yields no + information, do nothing. + (1b) If it completes incorrectly then a coincidence + is found, process it. + (1c) If it does not complete, go to step 2. + 2. Then we start from the right and scan backward to the left as far + as possible until it "meets" the forward scan `f`. Let this scan + results in coset `b`. + (2a) if `f` and `b` overlap, then a coincidence + is found, process it. + (2b) if `f` and `b` are about to meet (with a length 1 gap between) + then a deduction is found, fill in the two new entries. + (2c) else the scan is incomplete, define a new coset + and continue scanning forward. + """ + f = coset + b = coset + i = 0 + j = len(word) - 1 + while True: + # scan forward as far as possible. + while i <= j and self.is_defined(f, word[i]): + f = self[f][word[i]] + i += 1 + # if complete + if i > j: + # if complete incorrectly a coincidence is found, process it. + if f != b: + self.coincidence(f, b) + # else the scan yields no information + return + + # if scan forward is not completed then scan backward as + # far as possible until it meets the forward scan. + while j >= i and self.is_defined(b, self.inv(word[j])): + b = self[b][self.inv(word[j])] + j -= 1 + + # if f and b overlap then a coincidence is found. + if j < i: + self.coincidence(f, b) + return + # if f and b are about to meet a deduction is found. + elif j == i: + self[f][word[i]] = b + self[b][self.inv(word[i])] = f + return + # else define a new coset and continue scanning forward. + else: + self.define(f, word[i]) + + def coincidence(self, coset1, coset2): + """ + Process a coincidence. When two cosets are found to be equivalent + the larger one is deleted: we need to copy all information in its + row to the equivalent row and modify all entries it occurs in the + table. During this process new coincidences may be found and they + are also processed. The subtle thing here is that we must keep + filling/deleting entries in pairs, i.e. we must make sure the table + always satisfy table[k][x] = l <===> table[l][y] = k (y = inv(x)). + """ + self.merge(coset1, coset2) + # at this stage the queue contains only one item, + # but new items may be added to it later. + while len(self.q) > 0: + e = self.q.pop(0) + for x in self.A: + if self.is_defined(e, x): + f = self[e][x] + y = self.inv(x) + # the entry (f, y) in the row f is also deleted since we + # must make sure all entries in the table come in pairs. + self.undefine(f, y) + # copy the information at (e, x) to (e1, x), + # but is (e1, x) already defined? + e1 = self.rep(e) + f1 = self.rep(f) + # if (e1, x) is already defined then + # a new coincidence is found. + if self.is_defined(e1, x): + self.merge(f1, self.table[e1][x]) + elif self.is_defined(f1, y): + self.merge(e1, self.table[f1][y]) + else: + # else (e1, x) is not defined, copy the information + # at (e, x) to (e1, x). + self[e1][x] = f1 + self[f1][y] = e1 + # keep in mind all changes come in pairs, + # so we also copy the information at (f, y) to (f1, y). + + def hlt(self): + """Run the HLT strategy (Haselgrove, Leech and Trotter). + """ + for word in self.H: + self.scan_and_fill(0, word) + + current = 0 + while current < len(self.table): + for rel in self.R: + if not self.is_alive(current): + break + self.scan_and_fill(current, rel) + + if self.is_alive(current): + for x in self.A: + if not self.is_defined(current, x): + self.define(current, x) + current += 1 + + def compress(self): + """ + Delete all dead cosets in the table. The live cosets are renumbered and + their entries are also updated. + """ + ind = -1 + for coset in range(len(self)): + if self.is_alive(coset): + ind += 1 + if ind != coset: + for x in self.A: + y = self[coset][x] + if y == coset: + self[ind][x] = ind + else: + self[ind][x] = y + self[y][self.inv(x)] = ind + + self.p = list(range(ind + 1)) + self.table = self.table[: len(self.p)] + + def swap(self, k, l): + """ + Swap two live cosets in the table. It's called in the `standardize()` + method after the table is compressed, but it also works for + non-compressed table. + """ + for x in self.A: + # swap the two rows k and l. + self[k][x], self[l][x] = self[l][x], self[k][x] + # modify all k/l in the table. + for coset in range(len(self)): + # this check is not required for a compressed table. + if self.is_alive(coset): + if self[coset][x] == k: + self[coset][x] = l + elif self[coset][x] == l: + self[coset][x] = k + + def standardize(self): + """Rearrange the cosets in the table to a standard form. + """ + # the next coset we want to encounter in the table. + next_coset = 1 + for coset in range(len(self)): + for x in self.A: + y = self[coset][x] + if y >= next_coset: + if y > next_coset: + self.swap(y, next_coset) + next_coset += 1 + if next_coset == len(self) - 1: + return + + def run(self, standard=False): + self.hlt() + self.compress() + if standard: + self.standardize() + + def get_words(self): + """ + Return the list of all words in this table. This method must be called + when the table is already compressed. + """ + from collections import deque + + gens = self.A if self.coxeter else self.A[::2] + result = [None] * len(self) + result[0] = () + q = deque([0]) + while q: + coset = q.popleft() + for x in gens: + new_coset = self[coset][x] + if result[new_coset] is None: + result[new_coset] = result[coset] + (x,) + q.append(new_coset) + return result diff --git a/src/polytopes/povray/mathutils.inc b/src/polytopes/povray/mathutils.inc new file mode 100644 index 0000000..06fee3b --- /dev/null +++ b/src/polytopes/povray/mathutils.inc @@ -0,0 +1,90 @@ +// Persistence of Vision Ray Tracer Scene Description File +// Vers: 3.7 +// Date: 2018/04/22 +// Auth: Zhao Liang mathzhaoliang@gmail.com + +/* +============================================================ +Some math helper functions for rendering curved 4d polytopes +============================================================ +*/ +#include "math.inc" + +// Rotate a 4d vector +#macro rotate4d(theta, phi, xi, vec) + #local a = cos(xi); + #local b = sin(theta)*cos(phi)*sin(xi); + #local c = sin(theta)*sin(phi)*sin(xi); + #local d = cos(theta)*sin(xi); + #local p = vec.x; + #local q = vec.y; + #local r = vec.z; + #local s = vec.t; + < a*p - b*q - c*r - d*s + , a*q + b*p + c*s - d*r + , a*r - b*s + c*p + d*q + , a*s + b*r - c*q + d*p > +#end + +// stereographic project a 4d vector to 3d space +#macro proj4d(p) + #local q = p / sqrt(p.x*p.x + p.y*p.y + p.z*p.z + p.t*p.t); + / (1.0 + 1e-8 - q.t) +#end + +// return the normal vector of a 3d plane passes through the +// projected points of two 4d vectors p1 and p2 +#macro get_clipping_plane(p1, p2) + #local q1 = proj4d(p1); + #local q2 = proj4d(p2); + #local q12 = proj4d((p1+p2)/2); + VPerp_To_Plane(q1-q12, q2-q12) +#end + +// compute the signed distance of a vector to a plane, +// all vectors here are in 3d. +#macro distance_point_plane(p, p0, pnormal) + vdot(p-p0, pnormal) / vlength(pnormal) +#end + +// check if a vectors p is in the halfspace defined +// by the plane passes through p0 and has orientation pNormal. +#macro on_same_side(p, p0, pnormal) + #local result = false; + #local innprod = vdot(pnormal, p-p0); + #if (innprod > 0) + #local result = true; + #end + result +#end + +#macro get_sphere_info(nsides, pts) + #local result = array[3]; + #local A = proj4d(pts[0]); + #local B = proj4d(pts[1]); + #local C = proj4d(pts[2]); + #local P = 0; + #for (i, 0, nsides-1) + #local P = P + pts[i]; + #end + #local P = proj4d(P); + #local p1 = A - P; + #local p2 = B - P; + #local p3 = C - P; + #local det = p1.x * (p2.y*p3.z - p2.z*p3.y) - + p1.y * (p2.x*p3.z - p2.z*p3.x) + + p1.z * (p2.x*p3.y - p2.y*p3.x); + #local det = det * 2; + #if (abs(det) < 0.001) + #local result = array[3] { P, P, <-1, vlength(p1), 0> }; + #else + #local v1 = vcross(p1, p2) * vdot(p3, p3); + #local v2 = vcross(p3, p1) * vdot(p2, p2); + #local v3 = vcross(p2, p3) * vdot(p1, p1); + #local cen = (v1 + v2 + v3) / det; + #local sphere_rad = vlength(cen); + #local sphere_center = P + cen; + #local result = array[3] { sphere_center, P, <1, vlength(p1), sphere_rad> }; + #end + result +#end diff --git a/src/polytopes/povray/polychora_curved.pov b/src/polytopes/povray/polychora_curved.pov new file mode 100644 index 0000000..56c5d98 --- /dev/null +++ b/src/polytopes/povray/polychora_curved.pov @@ -0,0 +1,222 @@ +// Persistence of Vision Ray Tracer Scene Description File +// Vers: 3.7 +// Date: 2018/04/22 +// Auth: Zhao Liang mathzhaoliang@gmail.com + +/* +=============================== +Draw curved polychoron examples +=============================== +*/ + +#version 3.7; + +#include "colors.inc" +#include "textures.inc" +#include "mathutils.inc" +#include "polychora-data.inc" + +global_settings { + assumed_gamma 2.2 + max_trace_level 8 +} + +background { color bg_color } + +#declare num_segments = 30; +#declare face_transmit = 0.7; +#declare nvertices = dimension_size(vertices, 1); + +#macro Vertices(theta, phi, xi) + #local out = array[nvertices]; + #for(i, 0, nvertices-1) + #local out[i] = rotate4d(theta, phi, xi, vertices[i]); + #end + out +#end + +#declare rvs = Vertices(2*clock*pi, 0, 2*clock*pi); + +// project the arc between two 4d points on sphere S^3 to 3d +#macro get_arc(p1, p2) + sphere_sweep { + cubic_spline + num_segments + 3, + proj4d(p1), edge_size*get_size(proj4d(p1)) + #local ind=0; + #while (ind < num_segments) + #local q = proj4d(p1 + ind*(p2-p1)/num_segments); + q, edge_size*get_size(q) + #local ind = ind+1; + #end + proj4d(p2), edge_size*get_size(proj4d(p2)) + proj4d(p2), edge_size*get_size(proj4d(p2)) + } +#end + +#declare edgeColors = array[4] { Silver, Brass, Firebrick, GreenCopper }; +#declare faceColors = array[6] { White*0.9, Violet, Pink, Yellow, Brown, Magenta }; + +#declare vertexFinish = finish { ambient 0.5 diffuse 0.5 reflection 0.1 roughness 0.03 specular 0.5 } +#declare edgeFinish = finish { metallic ambient 0.2 diffuse 0.7 brilliance 6 reflection 0.25 phong 0.75 phong_size 80 } +#declare faceFinish = finish { + ambient 0.5 + diffuse 0.5 + reflection 0.1 + specular 0.4 + roughness 0.003 + irid { 0.3 thickness 0.2 turbulence 0.05 } + conserve_energy +} +#declare vertexTexture = texture { pigment { color Gold } finish { vertexFinish } } + +#macro edgeTexture(i) + texture { + pigment { edgeColors[i] } + finish { edgeFinish } + } +#end + +#macro faceTexture(i) + texture { + pigment { + color faceColors[i] + transmit face_transmit + } + finish { faceFinish } + } +#end + +#macro BubbleFace(faceType, nsides, pts) + #local sphere_info = get_sphere_info(nsides, pts); + #local isPlane = (sphere_info[2].x < 0); + #local face_size = sphere_info[2].y; + #local chosen = choose_face(faceType, face_size); + #if (chosen) + #if (isPlane) + #local face_center = sphere_info[0]; + #local pdist = vlength(face_center); + #local pnormal = vnormalize(face_center); + plane { + pnormal, pdist + faceTexture(faceType) + clipped_by { + union { + #local ind = 0; + #while (ind < nsides) + #local ind2 = ind + 1; + #if (ind2 = nsides) + #local ind2 = 0; + #end + get_arc(pts[ind], pts[ind2]) + #local ind = ind + 1; + #end + } + } + } + #else + #local face_center = sphere_info[1]; + #local sphere_center = sphere_info[0]; + #local sphere_rad = sphere_info[2].z; + #local ind = 0; + #local planes = array[nsides]; + #local pts3d = array[nsides]; + #local dists = array[nsides]; + #local sides = array[nsides]; + #while (ind < nsides) + #local ind2 = ind + 1; + #if (ind2 = nsides) + #local ind2 = 0; + #end + #local planes[ind] = get_clipping_plane(pts[ind], pts[ind2]); + #local pts3d[ind] = proj4d(pts[ind]); + #local dists[ind] = distance_point_plane(0, pts3d[ind], planes[ind]); + #local sides[ind] = on_same_side(face_center, pts3d[ind], planes[ind]); + #if (sides[ind] != true) + #local planes[ind] = -planes[ind]; + #end + #local ind = ind+1; + #end + sphere { + sphere_center, sphere_rad + faceTexture(faceType) + #local ind = 0; + #while (ind < nsides) + clipped_by { plane { -planes[ind], dists[ind] } } + #local ind = ind+1; + #end + } + #end + #end +#end + +union { + #for (ind, 0, nvertices-1) + #local q = proj4d(rvs[ind]); + sphere { + q, vertex_size*get_size(q) + texture { vertexTexture } + } + #end + + #local edge_types = dimension_size(edges, 1); + #for (i, 0, edge_types-1) + #local edge_list = edges[i]; + #local nedges = dimension_size(edge_list, 1); + #for (j, 0, nedges-1) + #local v1 = edge_list[j][0]; + #local v2 = edge_list[j][1]; + object { + get_arc(rvs[v1], rvs[v2]) + edgeTexture(i) + } + #end + #end + + #local face_types = dimension_size(faces, 1); + #for (i, 0, face_types-1) + #local face_list = faces[i]; + #local nfaces = dimension_size(face_list, 1); + #local nsides = dimension_size(face_list, 2); + #for (j, 0, nfaces-1) + #local points = array[nsides]; + #for (k, 0, nsides-1) + #local points[k] = rvs[face_list[j][k]]; + #end + BubbleFace(i, nsides, points) + #end + #end + scale 20 + rotate obj_rotation +} + +camera { + location camera_loc + look_at <0, 0, 0> + angle 40 + right x*image_width/image_height + up y +} + +light_source { + <-1, 1, 2> * 200 + color rgb 1 +} + +light_source { + <1, 3, 1> * 200 + color rgb 1 +} + +#if (use_area_light) + light_source { + camera_loc * 0.25 + color rgb 1 + area_light + x*5, y*5 + 3, 3 + adaptive 1 + orient + jitter + } +#end diff --git a/src/polytopes/povray/polyhedra_animation.pov b/src/polytopes/povray/polyhedra_animation.pov new file mode 100644 index 0000000..0cb6e08 --- /dev/null +++ b/src/polytopes/povray/polyhedra_animation.pov @@ -0,0 +1,130 @@ +// Persistence of Vision Ray Tracer Scene Description File +// Vers: 3.7 +// Date: 2018/02/22 +// Auth: Zhao Liang mathzhaoliang@gmail.com + +/* +========================================= +Make Animations of Rotating 3d Polyhedron +========================================= +*/ + +#version 3.7; + +global_settings { assumed_gamma 2.2 } + +#include "colors.inc" +#include "polytope-data.inc" + +// set background color +background { White } + +// size of vertices and edges +#declare vertexRad = 0.035; +#declare edgeRad = 0.02; + +// set transparency of faces +#declare face_filter = 0.4; + +#declare vertexColor = SkyBlue; +#declare edgeColors = array[3]{ Orange, GreenYellow, Maroon }; +#declare faceColors = array[3]{ Pink, LightSteelBlue, Brass }; +#declare edgeFinish = finish { ambient 0.2 diffuse 0.5 reflection 0.1 specular 0.6 roughness 0.01 } +#declare faceFinish = finish { ambient 0.5 diffuse 0.5 specular 0.6 roughness 0.005 } + +// texture for the faces in the i-th orbit +#macro faceTexture(i) + texture { + pigment { faceColors[i] filter face_filter } + finish { faceFinish } + } +#end + +// render the polytope +union { + // get the number of vertices + #local nvertices = dimension_size(vertices, 1); + + // render the vertives using spheres + #for (ind, 0, nvertices-1) + sphere { + vertices[ind], vertexRad + texture { + pigment { color vertexColor } + finish { edgeFinish } + } + } + #end + + // get the number of orbits of edges + #local edge_types = dimension_size(edges, 1); + // for each orbit of edges + #for (i, 0, edge_types-1) + #local edge_list = edges[i]; + // get the number of edges in this orbit + #local nedges = dimension_size(edge_list, 1); + // render each edge by a cylinder connecting its two ends + #for (j, 0, nedges-1) + #local v1 = edge_list[j][0]; + #local v2 = edge_list[j][1]; + cylinder { + vertices[v1], vertices[v2], edgeRad + texture { + pigment { color edgeColors[i] } + finish { edgeFinish } + } + } + #end + #end + + // get the number of orbits of faces + #local face_types = dimension_size(faces, 1); + #for (i, 0, face_types-1) + #local face_list = faces[i]; + #local nfaces = dimension_size(face_list, 1); + #local nsides = dimension_size(face_list, 2); + #for (j, 0, nfaces-1) + // if this face is a triangle + #if (nsides = 3) + triangle { + vertices[face_list[j][0]], + vertices[face_list[j][1]], + vertices[face_list[j][2]] + faceTexture(i) + } + // else use the union of successive triangles (center, v_{i}, v_{i+1}) + // to render this face, this works for both regular and star polygons. + #else + #local center = <0, 0, 0>; + #for (k ,0, nsides-1) + #local center = center + vertices[face_list[j][k]]; + #end + #local center = center / nsides; + #for (ind, 0, nsides-1) + triangle { + center, + vertices[face_list[j][ind]], + vertices[face_list[j][mod(ind+1, nsides)]] + faceTexture(i) + } + #end + #end + #end + #end + + rotate <720*clock, 0, 360*clock> +} + +camera { + location <0, 0, 3.6> + look_at <0, 0, 0> + angle 40 + right x*image_width/image_height + up y + sky y +} + +light_source { + <1, 1, 3> * 100 + color rgb 1.3 +} diff --git a/src/polytopes/povray/polytope_animation.pov b/src/polytopes/povray/polytope_animation.pov new file mode 100644 index 0000000..9f51e98 --- /dev/null +++ b/src/polytopes/povray/polytope_animation.pov @@ -0,0 +1,141 @@ +// Persistence of Vision Ray Tracer Scene Description File +// Vers: 3.7 +// Date: 2019/01/04 +// Auth: Zhao Liang mathzhaoliang@gmail.com + +/* +======================================== +Make Animations of Rotating 4d Polytopes +======================================== +*/ + +#version 3.7; + +global_settings { + assumed_gamma 2.2 + max_trace_level 8 +} + +#include "colors.inc" +#include "textures.inc" +#include "polytope-data.inc" + +background { Black } + +#declare vertexRad = 0.02; +#declare edgeRad = 0.01; +#declare nvertices = dimension_size(vertices, 1); + +#macro rotate4d(theta, phi, xi, vec) + #local a = cos(xi); + #local b = sin(theta)*cos(phi)*sin(xi); + #local c = sin(theta)*sin(phi)*sin(xi); + #local d = cos(theta)*sin(xi); + #local p = vec.x; + #local q = vec.y; + #local r = vec.z; + #local s = vec.t; + < a*p - b*q - c*r - d*s + , a*q + b*p + c*s - d*r + , a*r - b*s + c*p + d*q + , a*s + b*r - c*q + d*p > +#end + +#macro proj3d(q) + / (2 - q.t) +#end + +#macro Vertices(theta, phi, xi) + #local out = array[nvertices]; + #for(i, 0, nvertices-1) + #local out[i] = proj3d(rotate4d(theta, phi, xi, vertices[i])); + #end + out +#end + +#declare rvs = Vertices(2*clock*pi, 0, 2*clock*pi); + +union { + #for (ind, 0, nvertices-1) + sphere { + rvs[ind], vertexRad + texture { + pigment { color SkyBlue } + finish { + ambient 0.2 + diffuse 0.5 + reflection 0.1 + specular 3 + roughness 0.003 + } + } + } + #end + + #local edge_types = dimension_size(edges, 1); + #for (i, 0, edge_types-1) + #local edge_list = edges[i]; + #local nedges = dimension_size(edge_list, 1); + #for (j, 0, nedges-1) + #local v1 = edge_list[j][0]; + #local v2 = edge_list[j][1]; + cylinder { + rvs[v1], rvs[v2], edgeRad + texture { + pigment { color Orange } + finish { + ambient .5 + diffuse .5 + reflection .0 + specular .5 + roughness 0.1 + } + } + } + #end + #end + + #local face_types = dimension_size(faces, 1); + #for (i, 0, face_types-1) + #local face_list = faces[i]; + #local nfaces = dimension_size(face_list, 1); + #local nsides = dimension_size(face_list, 2); + #for (j, 0, nfaces-1) + #if (nsides = 3) + triangle { + rvs[face_list[j][0]], + rvs[face_list[j][1]], + rvs[face_list[j][2]] + texture { NBglass finish { reflection 0 }} + } + #else + #local center = <0, 0, 0>; + #for (k ,0, nsides-1) + #local center = center + rvs[face_list[j][k]]; + #end + #local center = center / nsides; + #for (ind, 0, nsides-1) + triangle { + center, + rvs[face_list[j][ind]], + rvs[face_list[j][mod(ind+1, nsides)]] + texture { NBglass finish { reflection 0 }} + } + #end + #end + #end + #end + scale 4 +} + +camera { + location <0, 0, -8> + look_at 0 + right x*image_width/image_height + angle 40 +} + +light_source { + <500, 500, -1000> + color 1.3 +} diff --git a/src/polytopes/tc_examples/B24.yaml b/src/polytopes/tc_examples/B24.yaml new file mode 100644 index 0000000..f26f2fd --- /dev/null +++ b/src/polytopes/tc_examples/B24.yaml @@ -0,0 +1,15 @@ +# The 2-generator Burnside group of exponent 4. +# 1024 cosets +name: B24 +relators: + - a^4 + - b^4 + - (ab)^4 + - (Ab)^4 + - (aab)^4 + - (abb)^4 + - (aabb)^4 + - (Abab)^4 + - (aBab)^4 +subgroup-generators: + - a diff --git a/src/polytopes/tc_examples/Cube.yaml b/src/polytopes/tc_examples/Cube.yaml new file mode 100644 index 0000000..02878cf --- /dev/null +++ b/src/polytopes/tc_examples/Cube.yaml @@ -0,0 +1,11 @@ +# Symmetry group of the cube (S_4 x Z_2). +# 48 cosets +name: Cube +relators: + - a^2 + - b^2 + - c^2 + - (ab)^4 + - (bc)^3 + - (ac)^2 +subgroup-generators: diff --git a/src/polytopes/tc_examples/F27.yaml b/src/polytopes/tc_examples/F27.yaml new file mode 100644 index 0000000..453198e --- /dev/null +++ b/src/polytopes/tc_examples/F27.yaml @@ -0,0 +1,13 @@ +# The Fibonacci group F(2, 7), this is a cyclic group of order 29. +# 1 coset +name: F27 +relators: + - abC + - bcD + - cdE + - deF + - efG + - fgA + - gaB +subgroup-generators: + - a diff --git a/src/polytopes/tc_examples/G8723.yaml b/src/polytopes/tc_examples/G8723.yaml new file mode 100644 index 0000000..f07b5a1 --- /dev/null +++ b/src/polytopes/tc_examples/G8723.yaml @@ -0,0 +1,14 @@ +# The G8723 group, see +# J. J. Cannon, L. A. Dimino, G. Havas, and J. M. Watson. +# "Implementation and analysis of the Todd-Coxeter algorithm". +# Mathematics of Computation, 27:463--490, 1973. +# 448 cosets +name: G8723 +relators: + - a^8 + - b^7 + - (ab)^2 + - (Ab)^3 +subgroup-generators: + - aa + - Ab diff --git a/src/polytopes/tc_examples/HNO_1.yaml b/src/polytopes/tc_examples/HNO_1.yaml new file mode 100644 index 0000000..c9bd185 --- /dev/null +++ b/src/polytopes/tc_examples/HNO_1.yaml @@ -0,0 +1,12 @@ +# This is the first group depicted in +# George Havas, M.F.Newman, and E.A.O’Brien. +# Groups of deficiency zero. In G.Baumslag et al., editors, +# Geometric and Computational Perspectives on Infinite Groups, +# volume 25 of Amer. Math. Soc. DIMACS Series, pages 53–67. +# 6561 cosets +name: HNO_1 +relators: + - aaBab + - bbCbc + - cABcabc +subgroup-generators: diff --git a/src/polytopes/tc_examples/M12.yaml b/src/polytopes/tc_examples/M12.yaml new file mode 100644 index 0000000..3f16c15 --- /dev/null +++ b/src/polytopes/tc_examples/M12.yaml @@ -0,0 +1,12 @@ +# Mathieu group M12. +# 95040 cosets +name: M12 +relators: + - a^{11} + - b^2 + - c^2 + - (ab)^3 + - (ac)^3 + - (bc)^{10} + - a^2(bc)^2a(CB)^2 +subgroup-generators: