Generation bug fix (#163)

* updated pybind11

* small setup cleanup

* WIP face oriented edges

* bibliography

* version bump and bug fix in generation

Author: glyg

doc/bibliography/tyssue.bib

@@ -29,8 +29,8 @@
 ## Thanks to them!
 MAJOR = 0
 MINOR = 6
-MICRO = 1
-ISRELEASED = False
+MICRO = 2
+ISRELEASED = True
 VERSION = "%d.%d.%s" % (MAJOR, MINOR, MICRO)
 
 

setup.py

@@ -7,6 +7,7 @@
     LineTension,
     SurfaceTension,
     LengthElasticity,
+    PerimeterElasticity,
     FaceContractility,
     FaceAreaElasticity,
     FaceVolumeElasticity,
@@ -17,6 +18,7 @@
 
 sheet_effectors = [
     LengthElasticity,
+    PerimeterElasticity,
     FaceAreaElasticity,
     FaceVolumeElasticity,
     LineTension,

tests/dynamics/test_effectors.py

@@ -85,3 +85,4 @@ def test_spherical_mono():
     mean_basal = rho[organo.vert_df["segment"] == "basal"].mean()
     np.testing.assert_almost_equal(mean_apical, radius, decimal=1)
     np.testing.assert_almost_equal(mean_basal, radius + height, decimal=1)
+    assert np.all(organo.cell_df["vol"] > 0)

tests/generation/test_shapes.py

@@ -1,6 +1,8 @@
 import os
+import numpy as np
 import pandas as pd
 
+
 from tyssue import config, Sheet, SheetGeometry
 from tyssue.generation import three_faces_sheet
 from tyssue.io.hdf5 import load_datasets
@@ -50,3 +52,18 @@ def test_rod_update_height():
     SheetGeometry.update_all(sheet)
 
     assert (sheet.vert_df.rho.mean() - 0.96074585429756632) ** 2 < TOLERANCE
+
+
+def test_get_phis():
+    pth = os.path.join(stores_dir, "rod_sheet.hf5")
+    datasets = load_datasets(pth)
+    specs = config.geometry.rod_sheet()
+    sheet = Sheet("rod", datasets, specs)
+    SheetGeometry.update_all(sheet)
+    sheet.edge_df["sphi"] = SheetGeometry.get_phis(sheet)
+
+    assert np.all(
+        sheet.edge_df.sort_values(["face", "sphi"])
+        .groupby("face")
+        .apply(lambda df: np.roll(df["trgt"], 1) == df["srce"])
+    )

tests/geometry/test_sheet_geometry.py

@@ -118,6 +118,46 @@ def gradient(eptm):
         return -grad, grad
 
 
+# From Mapeng Bi et al. https://doi.org/10.1038/nphys3471
+class PerimeterElasticity(AbstractEffector):
+
+    dimensions = units.line_elasticity
+    magnitude = "perimeter_elasticity"
+    label = "Perimeter Elasticity"
+    element = "face"
+    specs = {
+        "face": {
+            "is_alive": 1,
+            "perimeter": 1.0,
+            "perimeter_elasticity": 0.1,
+            "prefered_perimeter": 3.81,
+        }
+    }
+
+    spatial_ref = "prefered_perimeter", units.length
+
+    @staticmethod
+    def energy(eptm):
+        return eptm.face_df.eval(
+            "0.5 * is_alive"
+            "* perimeter_elasticity"
+            "* (perimeter - prefered_perimeter)** 2"
+        )
+
+    @staticmethod
+    def gradient(eptm):
+
+        gamma_ = eptm.face_df.eval(
+            "perimeter_elasticity * is_alive" "*  (perimeter - prefered_perimeter)"
+        )
+        gamma = eptm.upcast_face(gamma_)
+
+        grad_srce = -eptm.edge_df[eptm.ucoords] * to_nd(gamma, len(eptm.coords))
+        grad_srce.columns = ["g" + u for u in eptm.coords]
+        grad_trgt = -grad_srce
+        return grad_srce, grad_trgt
+
+
 class FaceAreaElasticity(AbstractEffector):
 
     dimensionless = False

tyssue/dynamics/effectors.py

@@ -138,6 +138,8 @@ def extrude(apical_datasets, method="homotecy", scale=0.3, vector=[0, 0, -1]):
             if not col in datasets[elem]:
                 datasets[elem][col] = value
 
+    if (method == "normals") and (scale < 0):
+        datasets["edge"][["srce", "trgt"]] = datasets["edge"][["trgt", "srce"]]
     return datasets
 
 

tyssue/generation/modifiers.py

@@ -1,5 +1,4 @@
 import numpy as np
-
 from .base_geometry import BaseGeometry
 
 

tyssue/geometry/planar_geometry.py

@@ -2,7 +2,7 @@
 import pandas as pd
 
 from .planar_geometry import PlanarGeometry
-from .utils import rotation_matrix
+from .utils import rotation_matrix, rotation_matrices
 
 
 class SheetGeometry(PlanarGeometry):
@@ -186,8 +186,8 @@ def face_rotation(sheet, face, psi=0):
         else:
             return np.dot(rotation_matrix(psi, [0, 0, 1]), r1)
 
-    @classmethod
-    def face_projected_pos(cls, sheet, face, psi=0):
+    @staticmethod
+    def face_projected_pos(sheet, face, psi=0):
         """Returns the position of a face vertices projected on a plane
         perpendicular to the face normal, and translated so that the face
         center is at the center of the coordinate system
@@ -216,7 +216,7 @@ def face_projected_pos(cls, sheet, face, psi=0):
             sheet.vert_df.loc[face_orbit.values, sheet.coords].values
             - sheet.face_df.loc[face, sheet.coords].values
         )
-        u, s, rotation = np.linalg.svd(rel_pos.astype(np.float), full_matrices=False)
+        _, _, rotation = np.linalg.svd(rel_pos.astype(np.float), full_matrices=False)
         # rotation = cls.face_rotation(sheet, face, psi=psi)
         if psi != 0:
             rotation = np.dot(rotation_matrix(psi, [0, 0, 1]), rotation)
@@ -225,6 +225,49 @@ def face_projected_pos(cls, sheet, face, psi=0):
         )
         return rot_pos
 
+    @staticmethod
+    def face_rotations(sheet):
+        """Returns the (sheet.Ne, 3, 3) array of rotation matrices
+        such that each rotation aligns the coordinate system along face normals
+
+        """
+        normals = sheet.edge_df.groupby("face")[sheet.ncoords].mean()
+        normals = normals / np.linalg.norm(normals, axis=0)
+        normals = sheet.upcast_face(normals)
+
+        n_xy = np.linalg.norm(normals[["nx", "ny"]])
+        theta = -np.arctan2(n_xy, normals.nz)
+
+        direction = np.array(
+            [normals.ny.to_numpy(), -normals.nx.to_numpy(), np.zeros(sheet.Ne)]
+        ).T
+        rots = rotation_matrices(theta, direction)
+
+        return rots
+
+    @classmethod
+    def get_phis(cls, sheet):
+        """Returns the 'latitude' of the vertices in the plane perpendicular
+        to each face normal. For not-too-deformed faces, sorting vertices by this
+        gives clockwize orientation.
+
+        I think not-too-deformed means starconvex here.
+
+        """
+        rots = cls.face_rotations(sheet)
+
+        rel_srce_pos = (
+            sheet.edge_df[["sx", "sy", "sz"]]
+            - sheet.edge_df[["fx", "fy", "fz"]].to_numpy()
+        ).to_numpy()
+        rotated = np.einsum("ikj, ik -> ij", rots, rel_srce_pos)
+        return np.arctan2(rotated[:, 1], rotated[:, 0])
+
+    @classmethod
+    def sort_oriented_edges(cls, sheet):
+        phis = cls.get_phis(sheet)
+        sheet.edge_df = sheet.edge_df.sort_values(["face", phis]).reset_index()
+
 
 class ClosedSheetGeometry(SheetGeometry):
     """Geometry for a closed 2.5D sheet.

tyssue/geometry/sheet_geometry.py

@@ -74,3 +74,43 @@ def rotation_matrix(angle, direction):
         ]
     )
     return R
+
+
+def rotation_matrices(angle, direction):
+    """ Return an (N, 3, 3) array of rotation matrices
+    along N angles and N directions
+
+
+
+    Parameters
+    ----------
+    angle : np.ndarray of shape (N,)
+        array of rotation angles
+
+    directions : np.ndarray of shape (N, 3)
+        array of rotation vectors
+
+    Returns
+    -------
+    rots : np.ndarray of shape (N, 3, 3)
+        the array of rotation matrices
+    """
+    direction = direction / np.linalg.norm(direction, axis=1)[:, None]
+
+    sint, cost = np.sin(angle), np.cos(angle)
+
+    rots = np.zeros((cost.size, 3, 3))
+    for i in range(3):
+        rots[:, i, i] = cost
+    rots += np.einsum("ij, ik -> ijk", direction, direction) * (1 - cost[:, None, None])
+
+    direction *= sint[:, None]
+
+    rots[:, 0, 1] -= direction[:, 2]
+    rots[:, 0, 2] += direction[:, 1]
+    rots[:, 1, 0] += direction[:, 2]
+    rots[:, 1, 2] -= direction[:, 0]
+    rots[:, 2, 0] -= direction[:, 1]
+    rots[:, 2, 1] += direction[:, 0]
+
+    return rots