cleanup, optimization for face bentness

TODO.md

@@ -1,6 +1,11 @@
 # To Do
 * Clean up the code
 * Document
+* Test cases
+* Improve coloring
+* Add 2d mode for `view_off.py`
+* Port over the cellular automata on polyhedra faces script from the old `geogrid` module
+* Refactor `figures_*.py`
 * Replace the stitching logic in `antitile/sgs.py` so that for Class I and II subdivision, the faces do not need to be oriented. In theory this would allow SGS to operate on tilings on unorientable surfaces, although self-intersecting surfaces might have quirks to deal with.
 * Allow "mixed" triangle & quad grids, at least for Class I. (Triangle class I and quad class I work together, but triangle class II works with particular quad class III subdivisions, not every triangle class III has a compatible quad class III, and quad class II doesn't work with triangle at all. Probably need to work out the theory of what's compatible with what.)
 * Allow arbitrary face figures other than SGS (would probably need whole new program)

antitile/projection.py

@@ -7,20 +7,20 @@
 from . import xmath, breakdown
 
 #generic methods, valid on any-d, any shape
-def square_to_quad(xy, abcd):
+def square_to_quad(xy, base_pts):
     """Transforms a square in [0,1]^2 to a (possibly skew) quadrilateral
-    defined by abcd.
+    defined by base_pts.
 
     Args:
         xy: The 2d coordinates of the point. Last element of shape should be 2.
             If not, treated like xy[..., :2]: the values in xy[..., 2:] are
             ignored.
-        abcd: The coordinates of the quadrilateral. Should be in
+        base_pts: The coordinates of the quadrilateral. Should be in
             counterclockwise order to maintain orientation. First element of
             shape should be 4.
     Returns:
-        Coordinates in whatever space abcd was defined in."""
-    a, b, c, d = abcd[0], abcd[1], abcd[2], abcd[3]
+        Coordinates in whatever space base_pts was defined in."""
+    a, b, c, d = base_pts[0], base_pts[1], base_pts[2], base_pts[3]
     x, y = xy[..., 0], xy[..., 1]
     return a + (b-a)*x + (d-a)*y + (a-b+c-d)*x*y
 
@@ -117,54 +117,54 @@ def tri_naive_slerp(bary, base_pts):
     b = np.sin(angle * bary) / np.sin(angle)
     return b.dot(base_pts)
 
-def tri_areal(beta, triangle):
+def tri_areal(beta, base_pts):
     """Given a triangle and spherical areal coordinates, returns the vectors
     cooresponding to those coordinates.
 
     Args:
         beta: spherical areal coordinates
-        triangle: vertices of the spherical triangle in a 3x3 array
+        base_pts: vertices of the spherical triangle in a 3x3 array
 
     Returns: Points on the sphere
 
     >>> beta = np.array([[ 1, 0.8, 0.6, 0.4, 0.2, 0.0 ],
     ...                  [ 0, 0.2, 0.2, 0.2, 0.5, 0.5 ],
     ...                  [ 0, 0.0, 0.2, 0.4, 0.3, 0.5 ]]).T
-    >>> triangle = np.eye(3)
-    >>> triangle[2,1] = 1
-    >>> triangle = normalize(triangle)
-    >>> from_sph_areal_coords(beta, triangle) # doctest: +NORMALIZE_WHITESPACE
+    >>> base_pts = np.eye(3)
+    >>> base_pts[2,1] = 1
+    >>> base_pts = normalize(base_pts)
+    >>> from_sph_areal_coords(beta, base_pts) # doctest: +NORMALIZE_WHITESPACE
     array([[ 1.        ,  0.        ,  0.        ],
            [ 0.97123031,  0.23814214,  0.        ],
            [ 0.87887868,  0.44073244,  0.18255735],
            [ 0.68229421,  0.63250197,  0.3666277 ],
            [ 0.37935999,  0.88556406,  0.26807143],
            [ 0.        ,  0.92387953,  0.38268343]])
     """
-    triangle = xmath.normalize(triangle)
-    area = xmath.triangle_solid_angle(triangle[0], triangle[1], triangle[2])
+    base_pts = xmath.normalize(base_pts)
+    area = xmath.triangle_solid_angle(base_pts[0], base_pts[1], base_pts[2])
     area_i = beta * area
-    triangle_iplus1 = np.roll(triangle, -1, axis=0)
-    triangle_iplus2 = np.roll(triangle, 1, axis=0)
+    base_pts_iplus1 = np.roll(base_pts, -1, axis=0)
+    base_pts_iplus2 = np.roll(base_pts, 1, axis=0)
     #FIXME whytf is this commented statement not equivalent to below?
 #    L = ((1 + np.cos(area_i))[:, np.newaxis]*
-#         np.cross(triangle_iplus1, triangle_iplus2) -
+#         np.cross(base_pts_iplus1, base_pts_iplus2) -
 #         np.sin(area_i)[:, np.newaxis]*
-#         (triangle_iplus1 + triangle_iplus2)).transpose((0,2,1))
+#         (base_pts_iplus1 + base_pts_iplus2)).transpose((0,2,1))
     L0 = ((1 + np.cos(area_i[..., 0]))[..., np.newaxis]*
-          np.cross(triangle[1], triangle[2]) -
+          np.cross(base_pts[1], base_pts[2]) -
           np.sin(area_i[..., 0])[..., np.newaxis]*
-          (triangle[1] + triangle[2]))
+          (base_pts[1] + base_pts[2]))
     L1 = ((1 + np.cos(area_i[..., 1]))[..., np.newaxis]*
-          np.cross(triangle[2], triangle[0]) -
+          np.cross(base_pts[2], base_pts[0]) -
           np.sin(area_i[..., 1])[..., np.newaxis]*
-          (triangle[2] + triangle[0]))
+          (base_pts[2] + base_pts[0]))
     L2 = ((1 + np.cos(area_i[..., 2]))[..., np.newaxis]*
-          np.cross(triangle[0], triangle[1]) -
+          np.cross(base_pts[0], base_pts[1]) -
           np.sin(area_i[..., 2])[..., np.newaxis]*
-          (triangle[0] + triangle[1]))
+          (base_pts[0] + base_pts[1]))
     L = np.stack([L0, L1, L2], axis=-2)
-    h = np.sin(area_i)*(1 + np.sum(triangle_iplus1*triangle_iplus2, axis=-1))
+    h = np.sin(area_i)*(1 + np.sum(base_pts_iplus1*base_pts_iplus2, axis=-1))
     return np.linalg.solve(L, h)
 
 def triangles_method2(lindex, base_pts, freq):
@@ -216,7 +216,7 @@ def square_naive_slerp(xy, base_pts):
     angle = xmath.central_angle_equilateral(base_pts)
     x, y = xy[..., 0], xy[..., 1]
     sx = np.sin(x*angle)
-    sy = np.sin(y*angle)    
+    sy = np.sin(y*angle)
     scx = np.sin((1-x)*angle)
     scy = np.sin((1-y)*angle)
     a = scx * scy
@@ -226,36 +226,44 @@ def square_naive_slerp(xy, base_pts):
     mat = np.stack([a, b, c, d], axis=-1) / np.sin(angle)**2
     return mat.dot(base_pts)
 
-def square_slerp(xy, abcd):
+def square_naive_slerp_2(xy, base_pts):
+    angle = xmath.central_angle_equilateral(base_pts)
+    x, y = xy[..., 0], xy[..., 1]
+    a = (1-x)*(1-y)
+    b = x*(1-y)
+    c = x*y
+    d = (1-x)*y
+    mat = np.sin(np.stack([a, b, c, d], axis=-1)*angle) / np.sin(angle)
+    return mat.dot(base_pts)
+
+def square_slerp(xy, base_pts):
     """Transforms a square in [0,1]^2 to a spherical quadrilateral
-    defined by abcd, using spherical linear interpolation
+    defined by base_pts, using spherical linear interpolation
 
     Args:
         xy: The 2d coordinates of the point. Last element of shape should be 2.
             If not, treated like xy[..., :2]: the values in xy[..., 2:] are
             ignored.
-        abcd: The coordinates of the quadrilateral. Should be in
+        base_pts: The coordinates of the quadrilateral. Should be in
             counterclockwise order to maintain orientation. Shape should be
             (4, ..., 3)
     Returns:
         Coordinates on the sphere."""
-    #FIXME returns coordinates off the sphere.
-    #square_slerp(xy, abcd) == square_slerp(xy[:,::-1], abcd[[0,3,2,1]]) ?
-    #if not, what can we do to make it isotropic?
-    a, b, c, d = abcd[0], abcd[1], abcd[2], abcd[3]
+    #FIXME not symmetric on irregular faces
+    a, b, c, d = base_pts[0], base_pts[1], base_pts[2], base_pts[3]
     x, y = xy[..., 0], xy[..., 1]
     ab = xmath.slerp(a, b, x)
     dc = xmath.slerp(d, c, x)
     result = xmath.slerp(ab, dc, y)
     return result
 
-def square_intersections(lindex, abcd, freq):
+def square_intersections(lindex, base_pts, freq):
     """Transforms a square to a spherical quadrilateral using the method of
     intersections"""
     n, m = freq
     preframe = breakdown.frame_square(n, m)
     frame = square_slerp(preframe[..., np.newaxis, :],
-                         abcd[:, np.newaxis, np.newaxis, np.newaxis])
+                         base_pts[:, np.newaxis, np.newaxis, np.newaxis])
     gc_normals = np.cross(frame[..., 0, :], frame[..., 1, :])
     index = np.arange(2)
     pairs = gc_normals[index, lindex[:, index]]
@@ -265,13 +273,13 @@ def square_intersections(lindex, abcd, freq):
     # test to see if the points are on the correct side of the sphere
     # take the dot product of these vectors with the center of the
     # base face. if it's positive, it's right, if not, negate it
-    center = np.sum(abcd, axis=0)#don't need to normalize this
+    center = np.sum(base_pts, axis=0)#don't need to normalize this
     sign_correct = np.sum(center*ptx, axis=-1, keepdims=True) >= 0
     result = np.where(sign_correct, ptx, -ptx)
-    result[(lindex[:, 0] == b)     & (lindex[:, 1] == 0)] = abcd[0]
-    result[(lindex[:, 0] == a + b) & (lindex[:, 1] == b)] = abcd[1]
-    result[(lindex[:, 0] == a)     & (lindex[:, 1] == a + b)] = abcd[2]
-    result[(lindex[:, 0] == 0)     & (lindex[:, 1] == a)] = abcd[3]
+    result[(lindex[:, 0] == b)     & (lindex[:, 1] == 0)] = base_pts[0]
+    result[(lindex[:, 0] == a + b) & (lindex[:, 1] == b)] = base_pts[1]
+    result[(lindex[:, 0] == a)     & (lindex[:, 1] == a + b)] = base_pts[2]
+    result[(lindex[:, 0] == 0)     & (lindex[:, 1] == a)] = base_pts[3]
 
     return result
 
@@ -346,6 +354,10 @@ def project_sphere(sphcoords, zfunc=np.arcsin, scale=180 / np.pi):
          4: lambda bkdn, abc, freq, tweak:
              square_naive_slerp(bkdn.coord, abc)}
 
+SLERP2 = {3: lambda bkdn, abc, freq, tweak: tri_naive_slerp(bkdn.coord, abc),
+         4: lambda bkdn, abc, freq, tweak:
+             square_naive_slerp_2(bkdn.coord, abc)}
+
 OTHER = {3: lambda bkdn, abc, freq, tweak: tri_areal(bkdn.coord, abc),
          4: lambda bkdn, abc, freq, tweak:
              square_slerp(bkdn.coord[:, np.newaxis], abc)}
@@ -359,28 +371,31 @@ def project_sphere(sphcoords, zfunc=np.arcsin, scale=180 / np.pi):
         4: _sq_cir}
 
 PROJECTIONS = {'flat':  FLAT,
-               'slerp': SLERP,
+               'nslerp': SLERP,
+               'nslerp2': SLERP2,
                'other': OTHER,
                'gc':    GC,
                'disk':  DISK}
 
+PARALLEL = ['nslerp', 'nslerp2', 'disk']
+
 if __name__ == "__main__":
     import matplotlib.pyplot as plt
     from matplotlib.collections import PolyCollection, LineCollection
 
     a, b = 2, 1
     bkdn = breakdown.Breakdown(a, b, 'q')
-    abcd = xmath.normalize(np.array([[1, 1, 1],
+    base_pts = xmath.normalize(np.array([[1, 1, 1],
                                      [1, -1, 1],
                                      [-1, -1, 1],
                                      [-1, 1, 1]]))
-    borders = xmath.slerp(abcd[:, np.newaxis],
-                          np.roll(abcd, -1, axis=0)[:, np.newaxis],
+    borders = xmath.slerp(base_pts[:, np.newaxis],
+                          np.roll(base_pts, -1, axis=0)[:, np.newaxis],
                           np.linspace(0, 1)[:, np.newaxis]).reshape((-1, 3))
 
-    sq1 = xmath.normalize(square_to_quad(bkdn.coord[:, np.newaxis], abcd))
-    sqslerp = xmath.normalize(square_slerp(bkdn.coord[:, np.newaxis], abcd))
-    sq2 = xmath.normalize(square_intersections(bkdn.lindex, abcd, (a, b)))
+    sq1 = xmath.normalize(square_to_quad(bkdn.coord[:, np.newaxis], base_pts))
+    sqslerp = xmath.normalize(square_slerp(bkdn.coord[:, np.newaxis], base_pts))
+    sq2 = xmath.normalize(square_intersections(bkdn.lindex, base_pts, (a, b)))
     badsq2 = np.linalg.norm(sq2, axis=-1) < 1E-6
     sq2[badsq2] = sq1[badsq2]
     #x = np.stack([sqslerp[:,0],sq2[:,0]],axis=-1)

antitile/sgs.py

@@ -29,28 +29,20 @@ def stitch_3(edge, bf, bkdn, index_0, index_1, freq):
                    #0, 1,  2, 3
 ROLLMAP = np.array([3, 6, 12, 9])
 ROLLMASK = 2**np.arange(4)
+ROLLOFFSET = [0, 1, 1+1j, 1j]
 
 def _rot_4(coord, n, freq):
     n = n%4
     a, b = freq
     coord = coord.astype(float)
     cco = xmath.float2d_to_complex(coord).flatten()
     rot_cco = cco * np.exp(1j*np.pi*n/2)
-    if n == 0:
-        offset = 0
-    elif n == 1:
-        offset = 1
-    elif n == 2:
-        offset = 1+1j
-    elif n == 3:
-        offset = 1j
-    shift_cco = rot_cco + offset
+    shift_cco = rot_cco + ROLLOFFSET[n]
     cxy = shift_cco * (a + b*1j) + b
     lindex = xmath.complex_to_float2d(cxy)
     return np.round(lindex).astype(int)
 
 def stitch_4(edge, bf, bkdn, index_0, index_1, freq):
-    #FIXME
     a, b = freq
     bkdn_0 = bkdn[index_0]
     bkdn_1 = bkdn[index_1]
@@ -96,8 +88,8 @@ def _find_dupe_verts(base, base_faces, rbkdn, freq, stitcher):
     #TODO replace this with np.unique when I update numpy
     matches = np.array(list({tuple(sorted(t)) for t in matches}))
     #if first one lies outside the base face, swap
-    index = rbkdn.group[matches[..., 0]] >= 200
-    matches[index] = matches[index, ::-1]
+    #index = rbkdn.group[matches[..., 0]] >= 200
+    #matches[index] = matches[index, ::-1]
     vno = len(rbkdn)
     conns = sparse.coo_matrix((np.ones(len(matches)),
                                (matches[:, 0], matches[:, 1])),
@@ -168,6 +160,7 @@ def subdiv(base, freq=(2, 0), proj='flat', tweak=False):
     return result
 
 def face_color_bf(poly):
+    #TODO refactor
     faces = poly.faces
     bf = np.array(poly.base_face)
     result = []
@@ -185,6 +178,7 @@ def face_color_bf(poly):
 
 
 def face_color_group(poly, fn=max):
+    #TODO refactor    
     faces = poly.faces
     group = np.array(poly.group)
     result = []
@@ -233,14 +227,16 @@ def parallel_approx(pts, normal):
     return q[..., np.newaxis] * normal
 
 def optimize_k(poly, base, measure, exact=True, normalize=True):
+    """Routine to optimize the factor k"""
     parallel_xyz = parallels(poly, base, exact)
-    result = minimize_scalar(objective, bracket=[0, 1],
+    result = minimize_scalar(objective, bracket=[0, 2],
                              args=(poly, parallel_xyz, measure, normalize))
     if not result.success:
         warnings.warn('Optimization routine did not converge')
     return result.x
 
 def objective(k, poly, parallel_xyz, measure, normalize=True):
+    """Objective function for the optimization routine"""
     test_v = parallel_sphere(poly.vertices, parallel_xyz, k)
     if normalize:
         test_v = xmath.normalize(test_v)
@@ -257,19 +253,31 @@ def energy(xyz, poly):
     Thomson's problem."""
     return tiling.energy(xyz)
 
+def bentness(xyz, poly):
+    """Measure: face bentness. Optimizes the max in the tiling, 
+    not the ratio of max to min."""
+    return tiling.bentness(xyz, poly).max()
+
 def edge_length(xyz, poly, spherical=False):
+    """Measure: edge length. Optimizes the ratio of max to min 
+    in the tiling."""
     result = tiling.edge_length(xyz, poly.edges, spherical)
     return result.max()/result.min()
 
 def face_area(xyz, poly, spherical=False):
+    """Measure: face area. Optimizes the ratio of max to min
+    in the tiling."""
     result = tiling.face_area(xyz, poly, spherical)
     return result.max()/result.min()
 
 def aspect_ratio(xyz, poly, spherical=False):
+    """Measure: aspect ratio. Optimizes the ratio of max to min
+    in the tiling."""
     result = tiling.aspect_ratio(xyz, poly, spherical)
     return result.max()/result.min()
 
 MEASURES = {'energy': energy,
+            'bent': bentness,
             'edges': edge_length,
             'aspect': aspect_ratio,
             'faces': face_area,