more general progress

MANIFEST.in

@@ -1 +1 @@
-include antitile/*.off
+include data/*.off

antitile/__init__.py

@@ -0,0 +1,4 @@
+from __future__ import division, absolute_import, print_function
+#import warnings
+# warnings.filterwarnings("ignore")
+from . import breakdown, flat, off, projection, sgs, tiling, xmath

antitile/projection.py

@@ -226,17 +226,19 @@ def square_intersections(lindex, abcd, freq):
 
 
 FLAT = {3: lambda bkdn, abc, freq, tweak: tri_bary(bkdn.coord, abc),
-        4: lambda bkdn, abc, freq, tweak: square_to_quad(bkdn.coord, abc)}
+        4: lambda bkdn, abc, freq, tweak:
+            square_to_quad(bkdn.coord[:, np.newaxis], abc)}
 
 SLERP = {3: lambda bkdn, abc, freq, tweak: tri_naive_slerp(bkdn.coord, abc),
-         4: lambda bkdn, abc, freq, tweak: square_slerp(bkdn.coord, abc)}
+         4: lambda bkdn, abc, freq, tweak:
+             square_slerp(bkdn.coord[:, np.newaxis], abc)}
 
 AREAL = {3: lambda bkdn, abc, freq, tweak: tri_areal(bkdn.coord, abc)}
 
-GC = {3: lambda bkdn, abc, freq, tweak: tri_intersections(bkdn.lindex,
-                                                          abc, freq, tweak),
-      4: lambda bkdn, abc, freq, tweak: square_intersections(bkdn.lindex,
-                                                             abc, freq)}
+GC = {3: lambda bkdn, abc, freq, tweak:
+            tri_intersections(bkdn.lindex, abc, freq, tweak),
+      4: lambda bkdn, abc, freq, tweak:
+          square_intersections(bkdn.lindex, abc, freq)}
 
 PROJECTIONS = {'flat':  FLAT,
                'slerp': SLERP,

antitile/sgs.py

@@ -6,23 +6,46 @@
 import numpy as np
 from numpy.linalg import norm
 from scipy.optimize import minimize_scalar
+from scipy import sparse
 from . import tiling, breakdown, projection, xmath
 
+def stitch_3(edge, bf, bkdn, index_0, index_1, freq):
+    a, b = freq
+    bkdn_0 = bkdn[index_0]
+    bkdn_1 = bkdn[index_1]
+    #first figure out how to roll the shapes so they meet at the edge
+    #TODO replace this with isin whenever that gets into numpy
+    notedge = np.nonzero(~np.in1d(bf, edge).reshape(bf.shape))[1]
+    roll = 2 - notedge
+    offset = np.array([a+b, a, 2*a+b])#from 0,0 to a,b
+    lindices = np.roll(bkdn_0.lindex, roll[0], axis=-1)
+    flipped = offset - np.roll(bkdn_1.lindex, roll[-1], axis=-1)
+    matches = np.nonzero(np.all(lindices[:, np.newaxis] ==
+                                flipped[np.newaxis], axis=-1))
+    l0 = np.nonzero(index_0)[0][matches[0]]
+    l1 = np.nonzero(index_1)[0][matches[1]]
+    return l0, l1
+
+def stitch_4(edge, bf, bkdn_0, bkdn_1, freq):
+    a, b = freq
+    offset = np.array([a+2*b, b])#from 0,0 to a,b
 
 def subdiv(base, freq=(2, 0), proj='flat', tweak=False):
     projections = projection.PROJECTIONS[proj]
     faces_dict = base.faces_by_size
     if any(x > 4 for x in faces_dict.keys()):
-        raise ValueError("Tiling contains at least one face with more than"
-                         " 4 sides. Try triangulating those faces first.")
+        raise ValueError("Tiling contains at least one face with more than "
+                         "4 sides. Try triangulating those faces first.")
     elif 3 in faces_dict and 4 in faces_dict:
         msg = ("Subdivision on mixed triangle-quadrilateral polyhedra "
                "is not yet implemented")
         raise NotImplementedError(msg)
     elif 3 in faces_dict:
         n = 3
+        stitcher = stitch_3
     elif 4 in faces_dict:
         n = 4
+        stitcher = stitch_4
     else:
         raise ValueError("Polyhedron has no faces")
     #create list of basically cartesion product of bkdn with base_faces
@@ -39,9 +62,48 @@ def subdiv(base, freq=(2, 0), proj='flat', tweak=False):
     arrays = [bf, group, coord, lindex, vertices]
     rbkdn = xmath.recordify(names, arrays)
     faces = np.concatenate([bkdn.faces + i*n_bkdn for i in range(n_bf)])
-    #remove redundant vertices/faces
-    #TODO
-
+    #find redundant vertices
+    base_edges = base.edges
+    base_edge_corr, base_face_corr = base.faces_by_edge(base_edges)
+    l0 = []
+    l1 = []
+    for i in range(len(base_edges)):
+        edge = base_edges[i]
+        index = base_edge_corr == i
+        facex = base_face_corr[index]
+        fn = len(facex)
+        if fn > 2:
+            warnings.warn("More than 2 faces meet at a single edge. "
+                          "Choosing 2 faces arbitrarily...")
+            facex = facex[:2]
+        elif fn < 2:#external edge, ignore it
+            continue
+        index_0 = rbkdn.base_face == facex[0]
+        index_1 = rbkdn.base_face == facex[1]
+        lx0, lx1 = stitcher(edge, base_faces[facex], rbkdn,
+                            index_0, index_1, freq)
+        l0.extend(lx0)
+        l1.extend(lx1)
+        #TODO
+    matches = np.stack([l0, l1], axis=-1)
+    matches.sort(axis=-1)
+    #TODO replace this with np.unique when I update numpy
+    matches = np.array(list({tuple(t) for t in matches}))
+    #if first one lies outside the base face, swap
+    index = rbkdn.group[matches[..., 0]] >= 100
+    matches[index] = matches[index, ::-1]
+    vno = len(rbkdn)
+    unique_index = np.ones(vno, dtype=bool)
+    unique_index[matches[..., 1]] = False
+
+    conns = sparse.coo_matrix((np.ones(len(matches)),
+                              (matches[:, 0], matches[:, 1])),
+                              shape=(vno, vno))
+    ncp, cp = sparse.csgraph.connected_components(conns)
+    print(unique_index.sum())
+    print(ncp)
+    print(cp)
+    #
     #project vertices
     proj_fun = projections[n]
     for i in range(n_bf):
@@ -99,7 +161,7 @@ def optimize_k(poly, base, measure, exact=True, normalize=True):
     result = minimize_scalar(objective, bracket=[0, 1],
                              args=(poly, parallel_xyz, measure, normalize))
     if ~result.success:
-        warnings.warn('optimization routine did not converge')
+        warnings.warn('Optimization routine did not converge')
     return result.x
 
 def objective(k, poly, parallel_xyz, measure, normalize=True):

antitile/tiling.py

@@ -4,7 +4,7 @@
 """
 import numpy as np
 #from scipy.linalg import circulant
-#from scipy import sparse
+from scipy import sparse
 
 
 class Tiling:
@@ -32,41 +32,9 @@ def faces_by_size(self):
                 result[x].append(face)
             else:
                 result[x] = [face]
+        for i in result:
+            result[i] = np.array(result[i])
         return result
-#    @property
-#    def face_orientation(self):
-#        """Returns the orientation of the points in each face in the
-#        tiling with respect to (0,0,0). 1 is counterclockwise, -1 is
-#        clockwise, 0 means it's mixed (probably some sort of weird
-#        self-intersecting thing) or lies on a plane through the origin."""
-#        lens = self.face_size
-#        min_len = lens.min()
-#        max_len = lens.max()
-#        orientation = np.zeros(lens.shape, dtype=np.int8)
-#        for i in range(min_len, max_len+1):
-#            index = lens == i
-#            i_faces = np.array([self.faces[i] for i in np.nonzero(index)[0]])
-#            i_pts = self.vertices[i_faces]
-#            if i == 3:
-#                dets = np.linalg.det(i_pts)
-#                cclock = dets > 0
-#                clock = dets < 0
-#            else:
-#                det_index = circulant(np.arange(4))[:, :3]
-#                dets_ind = np.linalg.det(i_pts[..., det_index, :])
-#                cclock = np.all(dets_ind > 0, axis=-1)
-#                clock =  np.all(dets_ind < 0, axis=-1)
-#            ix = np.nonzero(index)[0]
-#            orientation[ix[cclock]] = 1
-#            orientation[ix[clock]] = -1
-#        return orientation
-#
-#    def orient_faces(self):
-#        """Orient faces of tiling in counterclockwise order with
-#        respect (0,0,0)"""
-#        for face, o in zip(self.faces, self.face_orientation):
-#            if o == -1:
-#                face.reverse()
 
     @property
     def vertex_adjacency(self):
@@ -76,12 +44,32 @@ def vertex_adjacency(self):
     def face_adjacency(self):
         return NotImplemented
 
-    @property
-    def vertex_face_incidence(self):
-        return NotImplemented
 
     def faces_by_edge(self, edges):
-        return NotImplemented
+        """Given edges, lists the faces adjacent to the edges
+        Arguments:
+            edges: List of edges. An array of shape (n, 2). Lowest-numbered
+            vertex comes first.
+        Returns: two arrays:
+            edges: Index of each edge
+            faces: Index of each face"""
+        facesize = self.face_size
+        facedict = self.faces_by_size
+        ex = []
+        fx = []
+        for i in facedict:
+            face_no = np.nonzero(facesize == i)[0]
+            faces = facedict[i]
+            for j in range(i):
+                faceedge = faces[..., [j-1, j]]
+                condition = faceedge[..., 0] < faceedge[..., 1]
+                a = np.where(condition[..., np.newaxis], 
+                             faceedge, faceedge[..., ::-1])
+                edgel = np.all(edges[:, np.newaxis] == a[np.newaxis], axis=-1) 
+                edge_face = np.nonzero(edgel)
+                ex.append(edge_face[0])
+                fx.append(face_no[edge_face[1]])
+        return np.concatenate(ex), np.concatenate(fx)
 
 def edges_from_facelist(faces):
     """Given a list of faces, returns a list of edges."""

bin/sgs.py

@@ -43,6 +43,7 @@ def nonnegativeint(string, lowest=0):
     if x < lowest:
         msg = "must be greater than or equal to {y}"
         raise argparse.ArgumentTypeError(msg.format(y=lowest))
+    return x
 
 def posint(string):
     return nonnegativeint(string, 1)

data/3dihedron.off

@@ -0,0 +1,7 @@
+OFF
+3 2 3
+0 1 0
+-0.86602540378443864676372317 -0.5 0
+0.86602540378443864676372317 -0.5 0
+3 0 1 2 1
+3 2 1 0 2

data/4dihedron.off

@@ -0,0 +1,8 @@
+OFF
+4 2 4
+0 1 0
+-1 0 0
+0 -1 0
+1 0 0
+4 0 1 2 3 1
+4 3 2 1 0 2

data/cube.off

@@ -0,0 +1,16 @@
+OFF
+8 6 12
+0.5773502691896258 0.5773502691896258 0.5773502691896258
+0.5773502691896258 0.5773502691896258 -0.5773502691896258
+0.5773502691896258 -0.5773502691896258 0.5773502691896258
+0.5773502691896258 -0.5773502691896258 -0.5773502691896258
+-0.5773502691896258 0.5773502691896258 0.5773502691896258
+-0.5773502691896258 0.5773502691896258 -0.5773502691896258
+-0.5773502691896258 -0.5773502691896258 0.5773502691896258
+-0.5773502691896258 -0.5773502691896258 -0.5773502691896258
+4 1 5 4 0  1
+4 2 0 4 6  2
+4 3 1 0 2  3
+4 5 1 3 7  4
+4 6 4 5 7  5
+4 7 3 2 6  6

data/icosahedron.off

@@ -0,0 +1,34 @@
+OFF
+12 20 30
+0 0.5257311121191335 0.85065080835204
+0 0.5257311121191335 -0.85065080835204
+0 -0.5257311121191335 0.85065080835204
+0 -0.5257311121191335 -0.85065080835204
+0.5257311121191335 0.85065080835204 0
+0.5257311121191335 -0.85065080835204 0
+-0.5257311121191335 0.85065080835204 0
+-0.5257311121191335 -0.85065080835204 0
+0.85065080835204 0 0.5257311121191335
+0.85065080835204 0 -0.5257311121191335
+-0.85065080835204 0 0.5257311121191335
+-0.85065080835204 0 -0.5257311121191335
+3 1 3 11 1
+3 2 0 10 2
+3 2 5 8 3
+3 3 7 11 4
+3 4 1 6 5
+3 5 2 7 6
+3 5 3 9 7
+3 6 0 4 8
+3 7 2 10 9
+3 7 3 5 10
+3 8 0 2 11
+3 8 4 0 12
+3 8 5 9 13
+3 9 1 4 14
+3 9 3 1 15
+3 9 4 8 16
+3 10 0 6 17
+3 10 6 11 18
+3 11 6 1 19
+3 11 7 10 20

data/octahedron.off

@@ -0,0 +1,16 @@
+OFF
+6 8 12
+1 0 0
+-1 0 0
+0 0 1
+0 0 -1
+0 1 0
+0 -1 0
+3 0 2 5 1
+3 0 4 2 2
+3 1 3 5 3
+3 1 4 3 4
+3 2 4 1 5
+3 3 4 0 6
+3 5 2 1 7
+3 5 3 0 8

data/tetrahedron.off

@@ -0,0 +1,10 @@
+OFF
+4 4 6
+0.5773502691896258 0.5773502691896258 0.5773502691896258
+0.5773502691896258 -0.5773502691896258 -0.5773502691896258
+-0.5773502691896258 0.5773502691896258 -0.5773502691896258
+-0.5773502691896258 -0.5773502691896258 0.5773502691896258
+3 3 1 0 1
+3 2 0 1 2
+3 3 0 2 3
+3 2 1 3 4

setup.py

@@ -4,7 +4,7 @@
 from setuptools import setup
 
 setup(name='antitile',
-      version='20170316',
+      version='20170323',
       description='Tiling subdivision scripts, for use with Antiprism',
       url='http://github.com/brsr/antitile',
       author='B R S Recht',