Fix optimization issues

- Bugfix: sorting of cells by sensitivity
- Bugfix: quality caching
- Change output of optimization reports
- Remove caching (doesn't work and produces invalid geometry)
- Skip cells that were made illegal during optimization

src/classy_blocks/optimize/cell.py

@@ -1,4 +1,5 @@
 import abc
+import warnings
 from typing import ClassVar, Dict, List, Optional, Set, Tuple
 
 import numpy as np
@@ -128,50 +129,50 @@ def quality(self):
 
         # both 3D (cell) and 2d (face) use the same calculation but elements are different.
 
-        if self._quality is not None:
-            return self._quality
-
         quality = 0
         center = self.center
 
         def q_scale(base, exponent, factor, value):
             return factor * base ** (exponent * value) - factor
 
-        for orient, neighbour in self.neighbours.items():
-            i = self.side_names.index(orient)
+        try:
+            warnings.filterwarnings("error")
+
+            for orient, neighbour in self.neighbours.items():
+                i = self.side_names.index(orient)
 
-            ### non-orthogonality
-            # angles between sides and self.center-neighbour.center or, if there's no neighbour
-            # on this face, between face and self.center-face_center
-            if neighbour is None:
-                # Cells at the wall simply use center of the face on the wall
-                # instead of their neighbour's center
-                c2c = center - self.get_side_center(i)
-            else:
-                c2c = center - neighbour.center
+                ### non-orthogonality
+                # angles between sides and self.center-neighbour.center or, if there's no neighbour
+                # on this face, between face and self.center-face_center
+                if neighbour is None:
+                    # Cells at the wall simply use center of the face on the wall
+                    # instead of their neighbour's center
+                    c2c = center - self.get_side_center(i)
+                else:
+                    c2c = center - neighbour.center
 
-            c2cn = c2c / np.linalg.norm(c2c)
+                c2cn = c2c / np.linalg.norm(c2c)
 
-            angles = 180 * np.arccos(np.dot(self.get_side_normals(i), c2cn)) / np.pi
-            quality += np.sum(q_scale(1.25, 0.35, 0.8, angles))
+                angles = 180 * np.arccos(np.dot(self.get_side_normals(i), c2cn)) / np.pi
+                quality += np.sum(q_scale(1.25, 0.35, 0.8, angles))
+                ### cell inner angles
+                quality += np.sum(q_scale(1.5, 0.25, 0.15, abs(self.get_inner_angles(i))))
 
-            ### cell inner angles
-            quality += np.sum(q_scale(1.5, 0.25, 0.15, abs(self.get_inner_angles(i))))
+            ### aspect ratio: one number for the whole cell (not per side)
+            edge_lengths = self.get_edge_lengths()
+            side_max = max(edge_lengths)
+            side_min = min(edge_lengths) + VSMALL
+            aspect_factor = np.log10(side_max / side_min)
 
-        ### aspect ratio: one number for the whole cell (not per side)
-        edge_lengths = self.get_edge_lengths()
-        side_max = max(edge_lengths)
-        side_min = min(edge_lengths) + VSMALL
-        aspect_factor = np.log10(side_max / side_min)
+            quality += np.sum(q_scale(3, 2.5, 3, aspect_factor))
 
-        quality += np.sum(q_scale(3, 2.5, 3, aspect_factor))
+        except RuntimeWarning:
+            raise ValueError("Degenerate Cell") from RuntimeWarning
+        finally:
+            warnings.resetwarnings()
 
-        self._quality = quality
         return quality
 
-    def invalidate(self) -> None:
-        self._quality = None
-
     @property
     def min_length(self) -> float:
         return min(self.get_edge_lengths())
@@ -193,22 +194,17 @@ def get_side_normals(self, i):
         side_points = self.get_side_points(i)
         side_vector = side_points[1] - side_points[0]
 
-        if f.norm(side_vector) < VSMALL:
-            # two coincident points
-            return f.unit_vector(self.center - side_points[0])
-
         normal = np.cross(self.normal, side_vector)
 
-        if f.norm(normal) < VSMALL:
-            # a 180-degree angle
-            return f.unit_vector(self.center - side_points[0])
-
         return [f.unit_vector(normal)]
 
-    def get_inner_angles(self, _):
-        # this metric does not make sense for a line segment
-        # so only return side length
-        return np.zeros(2)
+    def get_inner_angles(self, i):
+        points = np.take(self.points, ((i - 1) % 4, i, (i + 1) % 4), axis=0)
+
+        side_1 = f.unit_vector(points[2] - points[1])
+        side_2 = f.unit_vector(points[0] - points[1])
+
+        return np.expand_dims(180 * np.arccos(np.dot(side_1, side_2)) / np.pi - 90, axis=0)
 
 
 class HexCell(CellBase):

src/classy_blocks/optimize/grid.py

@@ -114,15 +114,13 @@ def update(self, index: int, position: NPPointType) -> float:
         self.points[index] = position
 
         junction = self.junctions[index]
-        junction.invalidate()
 
         if len(junction.links) > 0:
             for indexed_link in junction.links:
                 indexed_link.link.leader = position
                 indexed_link.link.update()
 
                 self.points[indexed_link.follower_index] = indexed_link.link.follower
-                self.junctions[indexed_link.follower_index].invalidate()
 
             return self.quality
 

src/classy_blocks/optimize/iteration.py

@@ -9,32 +9,48 @@
 class ClampOptimizationData:
     """Quality tracking pre/post iteration"""
 
-    vertex_index: int
+    index: int
     grid_initial: float
     junction_initial: float
     junction_final: float = VBIG
     grid_final: float = VBIG
 
+    skipped: bool = False
+    rolled_back: bool = False
+
     def report_start(self) -> None:
-        report(f"{self.vertex_index:>6}", end="   ")
+        report(f"{self.index:>6}", end="   ")
         report(f"{self.grid_initial:.3e}", end="   ")
         report(f"{self.junction_initial:.3e}", end="   ")
 
     def report_end(self) -> None:
         report(f"{self.improvement: >11.0f}", end="   ")
         report(f"{self.grid_final:.3e}", end="   ")
 
-        rollback = "X" if self.rollback else ""
-        report(rollback)
+        comment = ""
+        if self.skipped:
+            comment = "Skip"
+        elif self.rolled_back:
+            comment = "Rollback"
+
+        report(comment)
+
+    def undo(self) -> None:
+        self.junction_final = self.junction_initial
+        self.grid_final = self.grid_initial
+
+    def rollback(self) -> None:
+        self.rolled_back = True
+        self.undo()
+
+    def skip(self) -> None:
+        self.skipped = True
+        self.undo()
 
     @property
     def improvement(self) -> float:
         return self.grid_initial - self.grid_final
 
-    @property
-    def rollback(self) -> bool:
-        return self.grid_final >= self.grid_initial
-
 
 class IterationData:
     """Data about a single iteration's progress"""
@@ -54,7 +70,7 @@ def improvement(self) -> float:
     def report_begin(self):
         report(f"Optimization iteration {self.index+1}:")
         # headers
-        report("Vertex     Initial       Local   Improvement       Final   Rollback")
+        report("Vertex     Initial       Local   Improvement       Final   Status")
 
     def report_end(self):
         report(f"Iteration {self.index+1} finished.", end=" ")

src/classy_blocks/optimize/junction.py

@@ -73,10 +73,6 @@ def add_clamp(self, clamp: ClampBase) -> None:
     def add_link(self, link: LinkBase, follower_index: int) -> None:
         self.links.append(IndexedLink(link, follower_index))
 
-    def invalidate(self) -> None:
-        for cell in self.cells:
-            cell.invalidate()
-
     @property
     def is_boundary(self) -> bool:
         """Returns True if this junction lies on boundary"""

src/classy_blocks/optimize/optimizer.py

@@ -50,29 +50,44 @@ def fquality(params):
             clamp.update_params(params)
             return self.grid.update(junction.index, clamp.position)
 
-        scipy.optimize.minimize(fquality, clamp.params, bounds=clamp.bounds, method=method)
+        try:
+            scipy.optimize.minimize(fquality, clamp.params, bounds=clamp.bounds, method=method)
 
-        reporter.junction_final = junction.quality
-        reporter.grid_final = self.grid.quality
-        reporter.report_end()
+            reporter.junction_final = junction.quality
+            reporter.grid_final = self.grid.quality
+
+            if reporter.improvement <= 0:
+                reporter.rollback()
 
-        if reporter.rollback:
+                clamp.update_params(initial_params)
+                self.grid.update(junction.index, clamp.position)
+        except ValueError:
+            # a degenerate cell (currently) cannot be untangled;
+            # try with a different junction
+            reporter.skip()
             clamp.update_params(initial_params)
             self.grid.update(junction.index, clamp.position)
 
+        reporter.report_end()
+
     def _get_sensitivity(self, clamp):
         """Returns maximum partial derivative at current params"""
         junction = self.grid.get_junction_from_clamp(clamp)
+        initial_params = copy.copy(clamp.params)
 
         def fquality(clamp, junction, params):
             clamp.update_params(params)
+            self.grid.update(junction.index, clamp.position)
             return junction.quality
 
         sensitivities = np.asarray(
             scipy.optimize.approx_fprime(clamp.params, lambda p: fquality(clamp, junction, p), epsilon=10 * TOL)
         )
+
+        clamp.update_params(initial_params)
+        self.grid.update(junction.index, clamp.position)
+
         return np.linalg.norm(sensitivities)
-        # return np.max(np.abs(sensitivities.flatten()))
 
     def optimize_iteration(self, method: MinimizationMethodType) -> None:
         clamps = sorted(self.grid.clamps, key=lambda c: self._get_sensitivity(c), reverse=True)
@@ -82,7 +97,7 @@ def optimize_iteration(self, method: MinimizationMethodType) -> None:
 
     def optimize(
         self, max_iterations: int = 20, tolerance: float = 0.1, method: MinimizationMethodType = "SLSQP"
-    ) -> None:
+    ) -> IterationDriver:
         """Move vertices, defined and restrained with Clamps
         so that better mesh quality is obtained.
 
@@ -113,6 +128,8 @@ def optimize(
 
         self.backport()
 
+        return driver
+
     @abc.abstractmethod
     def backport(self) -> None:
         """Reflect optimization results back to the original mesh/sketch"""
@@ -143,7 +160,7 @@ def backport(self):
 
     def auto_optimize(
         self, max_iterations: int = 20, tolerance: float = 0.1, method: MinimizationMethodType = "SLSQP"
-    ) -> None:
+    ) -> IterationDriver:
         """Adds a PlaneClamp to all non-boundary points and optimize the sketch.
         To include boundary points (those that can be moved along a line or a curve),
         add clamps manually before calling this method."""
@@ -154,4 +171,4 @@ def auto_optimize(
                 clamp = PlaneClamp(junction.point, junction.point, normal)
                 self.add_clamp(clamp)
 
-        super().optimize(max_iterations, tolerance, method)
+        return super().optimize(max_iterations, tolerance, method)

tests/test_optimize/optimize_fixtures.py

@@ -19,7 +19,7 @@ def positions(self):
                 [1, 0, 0],
                 [2, 0, 0],
                 [0, 1, 0],
-                [1.5, 1.5, 0],  # a moved vertex, should be [1, 1, 0]
+                [1.2, 1.6, 0],  # a moved vertex, should be [1, 1, 0]
                 [2, 1, 0],
                 [0, 2, 0],
                 [1, 2, 0],

tests/test_optimize/test_optimizer.py

@@ -1,3 +1,5 @@
+import unittest
+
 import numpy as np
 
 from classy_blocks.construct.flat.sketches.mapped import MappedSketch
@@ -6,6 +8,7 @@
 from classy_blocks.optimize.junction import ClampExistsError
 from classy_blocks.optimize.links import TranslationLink
 from classy_blocks.optimize.optimizer import MeshOptimizer, SketchOptimizer
+from classy_blocks.optimize.smoother import SketchSmoother
 from classy_blocks.util import functions as f
 from tests.test_optimize.optimize_fixtures import BoxTestsBase, SketchTestsBase
 
@@ -52,7 +55,7 @@ def test_optimize_linked(self):
 class SketchOptimizerTests(SketchTestsBase):
     def test_optimize_manual(self):
         sketch = MappedSketch(self.positions, self.quads)
-        clamp = PlaneClamp([1.5, 1.5, 0], [0, 0, 0], [0, 0, 1])
+        clamp = PlaneClamp([1.2, 1.6, 0], [0, 0, 0], [0, 0, 1])
 
         optimizer = SketchOptimizer(sketch)
         optimizer.add_clamp(clamp)
@@ -68,3 +71,71 @@ def test_optimize_auto(self):
         optimizer.auto_optimize(method="L-BFGS-B")
 
         np.testing.assert_almost_equal(sketch.positions[4], [1, 1, 0], decimal=3)
+
+
+class ComplexSketchTests(unittest.TestCase):
+    """Tests on a real-life case"""
+
+    # An degenerate starting configuration,
+    # smoothed to just barely valid
+
+    def setUp(self):
+        positions = np.array(
+            [
+                [0.01672874, 0.02687117, 0.02099406],
+                [0.01672874, 0.03602998, 0.02814971],
+                [0.01371287, 0.04465496, 0.03488828],
+                [0.00370317, 0.04878153, 0.0381123],
+                [-0.00689365, 0.04569677, 0.03570223],
+                [-0.01109499, 0.03743333, 0.02924612],
+                [-0.00613247, 0.02943629, 0.02299815],
+                [0.00472874, 0.02687117, 0.02099406],
+                [0.00872874, 0.02687117, 0.02099406],
+                [0.01272874, 0.02687117, 0.02099406],
+                [0.0110113, 0.03091146, 0.02415068],
+                [0.0110113, 0.03549087, 0.0277285],
+                [0.00950337, 0.03980335, 0.03109779],
+                [0.00449852, 0.04186664, 0.0327098],
+                [-0.00079989, 0.04032426, 0.03150476],
+                [-0.00290056, 0.03619254, 0.02827671],
+                [-0.0004193, 0.03219402, 0.02515273],
+                [0.0050113, 0.03091146, 0.02415068],
+                [0.0070113, 0.03091146, 0.02415068],
+                [0.0090113, 0.03091146, 0.02415068],
+            ],
+        )
+
+        face_map = [
+            # outer blocks
+            [8, 9, 1, 0],  # 0
+            [9, 10, 2, 1],  # 1
+            [10, 11, 3, 2],  # 2
+            [11, 12, 4, 3],  # 3
+            [12, 13, 5, 4],  # 4
+            [13, 14, 6, 5],  # 5
+            [14, 15, 7, 6],  # 6
+            # inner blocks
+            [15, 14, 9, 8],  # 7
+            [14, 13, 10, 9],  # 8
+            [13, 12, 11, 10],  # 9
+        ]
+
+        self.sketch = MappedSketch(positions, face_map)
+
+    def test_optimize(self):
+        smoother = SketchSmoother(self.sketch)
+        smoother.smooth()
+        optimizer = SketchOptimizer(self.sketch)
+        initial_quality = optimizer.grid.quality
+
+        # use a method that doesn't work well with this kind of problem
+        # (SLSQP seems to have issues with different orders of magnitude)
+        # so that a lot of rollback is required
+        iterations = optimizer.auto_optimize(method="SLSQP", tolerance=1e-3)
+
+        self.assertLess(optimizer.grid.quality, initial_quality)
+
+        last_val = 1e12
+        for iteration in iterations.iterations:
+            self.assertLess(iteration.final_quality, last_val)
+            last_val = iteration.final_quality