Merge pull request #5 from Baharis/development

✨ Improve `Routine`, `load`ing, add logging, fix dihedral handling

README.md

@@ -122,8 +122,8 @@ The following instructions are currently supported by picometer:
 - **Evaluation instructions**
   - measure `distance` between 2 selected objects; if the selection includes
     groups of atoms, measure closes distance to the group of atoms.
-  - measure `angle` between 2–6 selected objects; if the selection includes
-    (ordered) atoms, calculate direct or dihedral angle between presumed bonds.
+  - measure `angle` between 2–3 selected objects: planes, lines, or (ordered) atoms.
+  - measure `dihedral` andle between 4 individually-selected ordered centroids/atoms.
 
 
 ## Contributing

picometer/__main__.py

@@ -1,5 +1,6 @@
 from argparse import ArgumentParser, Namespace
-from picometer.routine import Routine
+from picometer.instructions import Routine
+from picometer.logging import add_file_handler
 from picometer.process import process
 import sys
 
@@ -20,6 +21,7 @@ def parse_args() -> Namespace:
 def main() -> int:
     args = parse_args()
     if filename := args.filename:
+        add_file_handler('picometer.log')
         routine = Routine.from_yaml(filename)
         process(routine)
     return 0

picometer/atom.py

@@ -8,11 +8,12 @@
 from numpy.linalg import norm
 import pandas as pd
 
-from picometer.shapes import degrees_between, Line, Plane, Shape, Vector3, are_synparallel
+from picometer.shapes import (are_synparallel, degrees_between, Line,
+                              Plane, Shape, Vector3)
 from picometer.utility import ustr2float
 
 
-log = logging.getLogger(__name__)
+logger = logging.getLogger(__name__)
 
 
 class Locator(NamedTuple):
@@ -29,7 +30,7 @@ def from_dict(cls, d: dict) -> 'Locator':
                        at=at if at else None)
 
     def __bool__(self):
-        return self.label and self.label != 'None'
+        return self.label is not None
 
 
 group_registry: Dict[str, List[Locator]] = {}
@@ -44,25 +45,28 @@ def __init__(self,
                  bf: BaseFrame = None,
                  table: pd.DataFrame = None,
                  ) -> None:
-        log.debug(f'Init AtomSet with {bf!r} and {table!r}')
+
+        logger.debug(f'Created atom set with {bf!r} and '
+                     f'{len(table) if table is not None else 0}-element table')
         self.base = bf
         self.table = table
 
-    def __len__(self):
+    def __len__(self) -> int:
         return len(self.table) if self.table is not None else 0
 
-    def __add__(self, other):
+    def __add__(self, other) -> 'AtomSet':
         if not (self.base or self.table):
             return other
         elif not (other.base or other.table):
             return self
-        return AtomSet(self.base, pd.concat([self.table, other.table], axis=0))
+        return self.__class__(self.base, pd.concat([self.table, other.table], axis=0))
 
     def __getitem__(self, item) -> 'AtomSet':
         return self.__class__(bf=self.base, table=self.table[item])
 
     @classmethod
     def from_cif(cls, cif_path: str, block_name: str = None) -> 'AtomSet':
+        """Initialize from cif file using hikari's `BaseFrame` and `CifFrame`"""
         bf = BaseFrame()
         cf = CifFrame()
         cf.read(cif_path)
@@ -87,11 +91,11 @@ def from_cif(cls, cif_path: str, block_name: str = None) -> 'AtomSet':
         return AtomSet(bf, atoms)
 
     @property
-    def fract_xyz(self):
+    def fract_xyz(self) -> np.ndarray:
         return np.vstack([self.table['fract_' + k].to_numpy() for k in 'xyz'])
 
     @property
-    def cart_xyz(self):
+    def cart_xyz(self) -> np.ndarray:
         return self.orthogonalise(self.fract_xyz)
 
     def fractionalise(self, cart_xyz: np.ndarray) -> np.ndarray:
@@ -105,7 +109,7 @@ def orthogonalise(self, fract_xyz: np.ndarray) -> np.ndarray:
     def locate(self, locators: Sequence[Locator]) -> 'AtomSet':
         """Convenience method to select multiple fragments from locators
         while interpreting and extending groups if necessary"""
-        log.debug(f'Locate {locators} in {self}')
+        logger.debug(f'Locate {locators} in {self}')
         new = AtomSet()
         assert len(locators) == 0 or isinstance(locators[0], Locator)
         for label, symm_op_code, at in locators:
@@ -121,9 +125,9 @@ def locate(self, locators: Sequence[Locator]) -> 'AtomSet':
 
     def select_atom(self, label_regex: str) -> 'AtomSet':
         mask = self.table.index == label_regex
-        if not any(mask):
+        if not any(mask):  # noqa: mask will in fact be Iterable
             mask = self.table.index.str.match(label_regex)
-        log.debug(f'Selected {sum(mask)} atoms with {label_regex=}')
+        logger.debug(f'Selected {sum(mask)} atoms with {label_regex=}')
         return self.__class__(self.base, deepcopy(self.table[mask]))
 
     def transform(self, symm_op_code: str) -> 'AtomSet':
@@ -136,7 +140,7 @@ def transform(self, symm_op_code: str) -> 'AtomSet':
         return self.__class__(self.base, data)
 
     @property
-    def centroid(self):
+    def centroid(self) -> np.ndarray:
         """A 3-vector with average atom position."""
         return self.cart_xyz.T.mean(axis=0)
 
@@ -165,7 +169,7 @@ def origin(self) -> Vector3:
         return self.centroid
 
     @origin.setter
-    def origin(self, new_origin):
+    def origin(self, new_origin) -> None:
         """Change origin to the new one provided in cartesian coordinates"""
         new_origin_fract = self.fractionalise(new_origin)
         delta = new_origin_fract - self.fractionalise(self.centroid)
@@ -178,16 +182,8 @@ def _angle(self, *others: 'Shape') -> float:
         assert all(o.kind is o.Kind.spatial for o in [self, *others])
         combined = sum(others, self)
         xyz = combined.cart_xyz.T
-        if len(combined) == 3:  # interior angle
-            return degrees_between(xyz[0] - xyz[1], xyz[2] - xyz[1])
-        elif 4 <= len(combined) <= 6:  # dihedral angle
-            plane1_dir = np.cross(xyz[0] - xyz[1], xyz[2] - xyz[1])
-            plane2_dir = np.cross(xyz[-3] - xyz[-2], xyz[-1] - xyz[-2])
-            twist_dir = np.cross(plane1_dir, plane2_dir)
-            sign = +1 if are_synparallel(twist_dir, xyz[2] - xyz[1]) else -1
-            return sign * degrees_between(plane1_dir, plane2_dir, normalize=False)
-        else:
-            return 'Input AtomSet must contain between 3 and 6 atoms'
+        assert len(combined) == 3, 'Input AtomSet must contain exactly 3 atoms'
+        return degrees_between(xyz[0] - xyz[1], xyz[2] - xyz[1])
 
     def _distance(self, other: 'Shape') -> float:
         if other.kind is self.Kind.spatial:
@@ -205,3 +201,14 @@ def _distance(self, other: 'Shape') -> float:
             norms = norm(deltas, axis=1)
             along = np.abs(np.dot(deltas, other.direction))
             return min(norms ** 2 - along ** 2)
+
+    def dihedral(self, *others: 'AtomSet') -> float:
+        assert all(o.kind is o.Kind.spatial for o in [self, *others])
+        combined = sum(others, self)
+        xyz = combined.cart_xyz.T
+        assert len(combined) == 4, 'Input AtomSet must contain exactly 4 atoms'
+        plane1_dir = np.cross(xyz[0] - xyz[1], xyz[2] - xyz[1])
+        plane2_dir = np.cross(xyz[-3] - xyz[-2], xyz[-1] - xyz[-2])
+        twist_dir = np.cross(plane1_dir, plane2_dir)
+        sign = +1 if are_synparallel(twist_dir, xyz[2] - xyz[1]) else -1
+        return sign * degrees_between(plane1_dir, plane2_dir, normalize=False)