coxeter-dice.py

import collections
import pprint
import sys

import attr
import numpy
import scipy.spatial
import graphviz
from cached_property import cached_property

class El(tuple):
    """Define our own tuple subclass, because we want elements to be ordered
    by word size first, then contents.
    """
    def __new__(cls, *args):
        return super().__new__(cls, args)

    def __repr__(self):
        return "{}{}".format(self.__class__.__name__, super().__repr__())

    def __lt__(self, other):
        if len(self) == len(other):
            return super().__lt__(other)
        return len(self) < len(other)

    def __le__(self, other):
        if len(self) == len(other):
            return super().__le__(other)
        return len(self) <= len(other)

    def __gt__(self, other):
        if len(self) == len(other):
            return super().__gt__(other)
        return len(self) > len(other)

    def __ge__(self, other):
        if len(self) == len(other):
            return super().__ge__(other)
        return len(self) >= len(other)

    def __getitem__(self, item):
        if isinstance(item, slice):
            return El(*super().__getitem__(item))
        return super().__getitem__(item)

    def __add__(self, other):
        return El(*super().__add__(other))

    def inverse(self):
        return El(*reversed(self))

@attr.s
class CoxeterGroup:
    n = attr.ib()
    orders = attr.ib()
    actives = attr.ib(default=set(), converter=set)
    snubs = attr.ib(default=set(), converter=set)

    # used to cache paths through the reduction space
    _reductions = attr.ib(init=False, repr=False, default=attr.Factory(dict))

    @cached_property
    def passives(self):
        """Generators not marked either active or snub."""
        return {a for a in range(self.n) if not (a in self.actives or a in self.snubs)}

    #
    # Linear Algebra Stuff
    # https://www.win.tue.nl/~amc/pub/CoxNotes.pdf (section 2.3)
    #
    
    @cached_property
    def order_matrix(self):
        """Return the Coxeter matrix, the matrix of generator orders."""
        m = 2 - numpy.identity(self.n, dtype=int)
        for (a, b), order in self.orders.items():
            m[a, b] = order
            m[b, a] = order
        return m

    @cached_property
    def bilinear_form(self):
        """Return the bilinear form defining inner products in the basis
        of normals to each reflection plane.
        """
        return -numpy.cos(numpy.pi / self.order_matrix)

    @cached_property
    def planes(self):
        """Return a list of normals for each reflection plane."""
        # the bilinear form matrix we have is in the basis defined
        # by the normal vectors. To get the normal vectors in a nice,
        # god-fearing orthonormal basis, we can use Cholesky decomposition
        # to find L such that L.dot(L.T) == bilinear_form, and that L
        # is the one that takes normal-basis vectors into nice-basis vectors
        # (whew!)
        L = numpy.linalg.cholesky(self.bilinear_form)
        return [L[a] for a in range(self.n)]

    @cached_property
    def seed_vertex(self):
        """Get the default, equidistant seed vertex."""
        distances = [0 if a in self.passives else (a + 1) for a in range(self.n)]
        return numpy.linalg.inv(numpy.array(self.planes)).dot(distances)

    def make_points(self, seed=None):
        """Generate (element, vertex) of the described shape."""
        if seed is None:
            seed = self.seed_vertex
        planes = self.planes
        vertices = []
        for el in self.vertices:
            point = seed.copy()
            parity = 0
            for a in el:
                if a in self.snubs:
                    parity += 1
                normal = planes[a]
                point -= 2 * normal.dot(point) * normal
            if parity % 2 == 0:
                vertices.append((el, point))
        return vertices

    #
    # Diagram stuff
    #

    @classmethod
    def from_diagram(cls, diagram):
        """https://en.wikipedia.org/wiki/Coxeter-Dynkin_diagram
        diagrams like 's4x o' or 'o-o4o'.
        can only handle lines, but still handy.
        """
        n = 0
        orders = {}
        actives = []
        snubs = []
        for c in diagram:
            if c in 'sox':
                if c in 's':
                    snubs.append(n)
                if c in 'x':
                    actives.append(n)
                n += 1
            else:
                if c in ' ':
                    order = 2
                elif c in '-':
                    order = 3
                else:
                    order = int(c)
                orders[(n-1, n)] = order
        return cls(n, orders, actives, snubs)

    def diagram(self):
        dot = graphviz.Graph(graph_attr={'rankdir': 'LR'})
        for a in range(self.n):
            if a in self.actives:
                dot.attr('node', shape='doublecircle')
                dot.attr('node', style='filled', fillcolor='black', fontcolor='white')
            elif a in self.snubs:
                dot.attr('node', shape='circle')
                dot.attr('node', style='filled', fillcolor='white', fontcolor='black')
            else:
                dot.attr('node', shape='circle')
                dot.attr('node', style='filled', fillcolor='black', fontcolor='white')
            dot.node(str(a))
        for (a, b), order in self.orders.items():
            if order <= 2:
                continue
            if order == 3:
                dot.edge(str(a), str(b))
            else:
                dot.edge(str(a), str(b), label=str(order))
        return dot

    def graph(self, stabilizer=None):
        equiv_map = {}
        nodes = self.elements
        if stabilizer:
            nodes = []
            for coset in self.cosets(stabilizer):
                m = min(coset)
                for el in coset:
                    equiv_map[el] = m
                nodes.append(m)
            nodes.sort()

        dot = graphviz.Graph()
        for el in nodes:
            if len(el) == 0:
                label = 'e'
            else:
                label = ' '.join(str(a) for a in el)
            dot.node(repr(el), label=label)
        edges = set()
        for el in nodes:
            for a in range(self.n):
                newel = self.reduce(el + El(a))
                newel = equiv_map.get(newel, newel)
                if newel != el and not (newel, el) in edges:
                    dot.edge(repr(el), repr(newel), label=str(a))
                    edges.add((el, newel))
        return dot

    def vertex_graph(self):
        nodes = self.vertices
        dot = graphviz.Graph()
        for el in nodes:
            if len(el) == 0:
                label = 'e'
            else:
                label = ' '.join(str(a) for a in el)
            dot.node(repr(el), label=label)
        for k, edges in self.edges.items():
            for v1, v2 in edges:
                dot.edge(repr(v1), repr(v2), label=str(k))
        return dot

    #
    # Word Problem Stuff
    # http://campus.murraystate.edu/academic/faculty/rdonnelly/CoxeterGroupsSeminar/TitsThm.pdf
    #

    def reduce_step(self, base):
        """yield all possible reductions that can be performed on a base"""
        if len(base) < 2:
            return
        base = list(base)
        # pair elimination (aa -> e)
        for i in range(1, len(base)):
            if base[i] == base[i - 1]:
                yield El(*(base[:i-1] + base[i+1:]))
        for a in range(self.n):
            for b in range(self.n):
                if a == b:
                    continue
                braid_size = self.order_matrix[a, b]
                for i in range(len(base) - braid_size + 1):
                    for j in range(braid_size):
                        if base[i + j] != (a if j % 2 else b):
                            break
                    else:
                        # we found a braid starting at i
                        braided = list(base)
                        for j in range(braid_size):
                            braided[i + j] = (b if j % 2 else a)
                        yield El(*braided)
                

    def reduce(self, root):
        """Reduce the given element to a canonical, minimal form."""
        if root in self._reductions:
            return self._reductions[root]

        new = collections.deque([root])
        found = {root}
        while new:
            el = new.popleft()
            for subel in self.reduce_step(el):
                if subel in self._reductions:
                    best = self._reductions[subel]
                    break
                if not subel in found:
                    found.add(subel)
                    new.append(subel)
            else:
                continue
            break
        else:
            best = min(found)

        for el in found:
            self._reductions[el] = best
        return best

    @cached_property
    def elements(self):
        """Return all minimal elements of this group."""
        generators = [El(a) for a in range(self.n)]
        els = {El()}
        els.update(generators)
        new = collections.deque(generators)
        while new:
            el = new.popleft()
            for a in range(self.n):
                if el[-1] == a:
                    continue
                subel = self.reduce(el + El(a))
                if not subel in els:
                    els.add(subel)
                    new.append(subel)
        els = list(els)
        els.sort()
        return els

    #
    # Misc. Group Stuff
    #

    def coset(self, base_set, el):
        """Return the coset of an element with respect to a base set."""
        return frozenset([self.reduce(base + el) for base in base_set])

    def cosets(self, base_set):
        """Return all cosets of a base set."""
        return {self.coset(base_set, el) for el in self.elements}

    def conjugate_subgroup(self, base_group, el):
        """Return the conjugate subgroup of a base group."""
        return {self.reduce(el.inverse() + base + el) for base in base_group}

    @cached_property
    def passive_elements(self):
        """Return the elements of the passive subgroup."""
        return [el for el in self.elements if all((a in self.passives) for a in el)]

    @cached_property
    def vertices(self):
        """Return minimal elements of each partition representing each element
        of the orbit of the seed vertex."""
        # https://en.wikipedia.org/wiki/Group_action#Orbit-stabilizer_theorem_and_Burnside's_lemma
        representatives = [min(coset) for coset in self.cosets(self.passive_elements)]
        representatives.sort()
        return representatives

    def fixed_edge_group(self, el):
        """Return the subgroup that fixes an edge, given an element that
        interchanges the vertices of that edge."""
        # must leave el(x) and x unchanged
        group = self.conjugate_subgroup(self.passive_elements, el)
        group.intersection_update(self.passive_elements)
        # or: swaps el(x) and x
        group.update(self.coset(group, el))
        return group

    @cached_property
    def passive_edge_elements(self):
        """Return the elements of the subgroups that fix each base edge."""
        # hmm. snubs...
        return {n: self.fixed_edge_group(El(n)) for n in self.actives.union(self.snubs)}

    @cached_property
    def edge_elements(self):
        """Return the minimal elements of the partition representing each
        element of the orbit of each base edge."""
        return {k: [min(coset) for coset in self.cosets(group)] for k, group in self.passive_edge_elements.items()}

    @cached_property
    def edges(self):
        """Return each edge generated by each base edge, in terms of which
        vertices they connect."""
        ret = {}
        for k, edges in self.edge_elements.items():
            ret[k] = [None] * len(edges)
            for i, edge in enumerate(edges):
                v1 = min(self.coset(self.passive_elements, edge))
                v2 = min(self.coset(self.passive_elements, El(k) + edge))
                ret[k][i] = (v1, v2)
        return ret

# jank quick PLY export
def ply_export(c, f):
    points = c.make_points()
    hull = scipy.spatial.ConvexHull([p[1] for p in points])
    
    f.write('ply\n')
    f.write('format ascii 1.0\n')
    f.write('comment maths made this shape, and maths will destroy it\n')
    f.write('element vertex {}\n'.format(len(points)))
    f.write('property float32 x\n')
    f.write('property float32 y\n')
    f.write('property float32 z\n')
    f.write('element face {}\n'.format(len(hull.simplices)))
    f.write('property list uint8 int32 vertex_indices\n')
    f.write('end_header\n')

    for el, pt in points:
        f.write('{} {} {}\n'.format(*pt))

    for simplex in hull.simplices:
        f.write('{}'.format(len(simplex)))
        for i in simplex:
            f.write(' {}'.format(i))
        f.write('\n')

if __name__ == "__main__":
    c = CoxeterGroup.from_diagram(' '.join(sys.argv[1:]))

    print(c)
    print('number of elements: {}'.format(len(c.elements)))
    print('number of vertices: {}'.format(len(c.vertices)))
    print('number of edges: {}'.format({k: len(s) for k, s in c.edges.items()}))

    c.diagram().render('diagram', format='png', cleanup=True)
    c.graph().render('graph', format='png', cleanup=True)
    c.graph(stabilizer=c.passive_elements).render('graph-active', format='png', cleanup=True)
    c.vertex_graph().render('graph-vertex', format='png', cleanup=True)
    with open('shape.ply', 'wt') as f:
        ply_export(c, f)