diff --git a/src/classy_blocks/optimize/cell.py b/src/classy_blocks/optimize/cell.py index 765ad7e..3c606b9 100644 --- a/src/classy_blocks/optimize/cell.py +++ b/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): diff --git a/src/classy_blocks/optimize/grid.py b/src/classy_blocks/optimize/grid.py index ef1545a..20a7194 100644 --- a/src/classy_blocks/optimize/grid.py +++ b/src/classy_blocks/optimize/grid.py @@ -114,7 +114,6 @@ 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: @@ -122,7 +121,6 @@ def update(self, index: int, position: NPPointType) -> float: indexed_link.link.update() self.points[indexed_link.follower_index] = indexed_link.link.follower - self.junctions[indexed_link.follower_index].invalidate() return self.quality diff --git a/src/classy_blocks/optimize/iteration.py b/src/classy_blocks/optimize/iteration.py index b0ca3af..1d3118c 100644 --- a/src/classy_blocks/optimize/iteration.py +++ b/src/classy_blocks/optimize/iteration.py @@ -9,14 +9,17 @@ 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=" ") @@ -24,17 +27,30 @@ 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=" ") diff --git a/src/classy_blocks/optimize/junction.py b/src/classy_blocks/optimize/junction.py index 5e8cc76..1c941bc 100644 --- a/src/classy_blocks/optimize/junction.py +++ b/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""" diff --git a/src/classy_blocks/optimize/optimizer.py b/src/classy_blocks/optimize/optimizer.py index 8b3b404..c10111f 100644 --- a/src/classy_blocks/optimize/optimizer.py +++ b/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) diff --git a/tests/test_optimize/optimize_fixtures.py b/tests/test_optimize/optimize_fixtures.py index b501e65..b866d74 100644 --- a/tests/test_optimize/optimize_fixtures.py +++ b/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], diff --git a/tests/test_optimize/test_optimizer.py b/tests/test_optimize/test_optimizer.py index f4a913a..10f55da 100644 --- a/tests/test_optimize/test_optimizer.py +++ b/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