Status quo: going nowhere

src/classy_blocks/modify/cell.py

@@ -50,6 +50,9 @@ def __init__(self, block: Block):
         self.side_indexes = [item[0] for item in q_map.items()]
         self.face_indexes = [item[1] for item in q_map.items()]
 
+        # caching
+        self._quality = 0
+
     def get_common_vertices(self, candidate: "Cell") -> Set[int]:
         """Returns indexes of common vertices between this and provided cell"""
         this_indexes = set(self.vertex_indexes)
@@ -121,6 +124,9 @@ def face_centers(self) -> NPPointListType:
 
     @property
     def quality(self) -> float:
+        if self._quality > 0:
+            return self._quality
+
         quality = 0
 
         center = self.center
@@ -184,6 +190,8 @@ def q_scale(base, exponent, factor, value):
 
             quality += np.sum(q_scale(3, 2.5, 3, aspect_factor))
 
+        self._quality = quality
+
         return quality
 
     @property

src/classy_blocks/modify/grid.py

@@ -67,3 +67,7 @@ def clamps(self) -> List[ClampBase]:
     def quality(self) -> float:
         """Returns summed qualities of all junctions"""
         return sum([cell.quality for cell in self.cells])
+
+    def clear_cache(self):
+        for cell in self.cells:
+            cell._quality = 0

src/classy_blocks/modify/optimizer.py

@@ -1,6 +1,6 @@
 import copy
 import time
-from typing import Literal
+from typing import List, Literal
 
 import numpy as np
 import scipy.optimize
@@ -15,6 +15,99 @@
 MinimizationMethodType = Literal["SLSQP", "L-BFGS-B", "Nelder-Mead", "Powell"]
 
 
+class ParamHandler:
+    N_ENTRIES = 6
+
+    def __init__(self, junction: Junction, clamp: ClampBase, index: int):
+        self.junction = junction
+        self.clamp = clamp
+        self.index = index
+
+        self.sign = 1
+        self.delta = TOL
+
+        self.param_log: List[float] = []
+        self.quality_log: List[float] = []
+
+    @property
+    def param(self) -> float:
+        return self.clamp.params[self.index]
+
+    @property
+    def slope(self) -> float:
+        return self.quality_log[-1] - self.quality_log[-2]
+
+    @property
+    def trend(self) -> float:
+        return np.diff(np.diff(self.quality_log))[-1]
+
+    @property
+    def last_param(self) -> float:
+        return self.param_log[-1]
+
+    @property
+    def last_delta(self) -> float:
+        return self.param_log[-1] - self.param_log[-2]
+
+    @property
+    def scale(self) -> float:
+        return abs(self.slope) / abs(self.last_delta + TOL)
+
+    @property
+    def done(self) -> bool:
+        if len(self.quality_log) < self.N_ENTRIES:
+            return False
+
+        return self.quality_log[-3] > self.quality_log[-2] and self.quality_log[-1] > self.quality_log[-2]
+
+    def update(self, param) -> None:
+        self.clamp.params[self.index] = param
+        self.clamp.update_params(self.clamp.params)
+
+        self.param_log.append(param)
+
+        if len(self.param_log) > self.N_ENTRIES:
+            self.param_log.pop(0)
+
+    def begin(self) -> None:
+        """Sets a new parameter"""
+        # make a few random tries to gather data
+        if len(self.param_log) < self.N_ENTRIES:
+            self.update(self.param + TOL)
+            return
+
+        if self.slope > 0:
+            # quality is getting worse!
+            self.sign *= -1
+            self.delta /= 2
+            return
+        else:
+            # slower improvement - larger time step
+            if self.trend > 0:
+                self.delta *= 0.8
+            else:
+                self.delta *= 1.2
+        self.update(self.param_log[-1] + self.delta * self.sign)
+
+    def end(self) -> None:
+        """Calculates new quality"""
+        new_quality = self.junction.quality
+
+        if len(self.quality_log) > self.N_ENTRIES - 1 and self.quality_log[-1] < new_quality:
+            # rollback!
+            self.param_log.pop()
+            self.quality_log.pop()
+            self.delta = TOL
+            self.sign *= -1
+            self.update(self.param_log[-1])
+            return
+
+        self.quality_log.append(new_quality)
+
+        if len(self.quality_log) > self.N_ENTRIES:
+            self.quality_log.pop(0)
+
+
 class Optimizer:
     """Provides tools for blocking optimization"""
 
@@ -91,25 +184,47 @@ def optimize(
 
         Within each iteration, all vertices will be moved, starting with the one with the most influence on quality.
         Lower tolerance values"""
-        driver = IterationDriver(max_iterations, tolerance)
+        # max_iterations = 2
+        # driver = IterationDriver(max_iterations, tolerance)
+
+        # start_time = time.time()
+
+        # while not driver.converged:
+        #     driver.begin_iteration(self.grid.quality)
+        #     self.optimize_iteration(method)
+        #     driver.end_iteration(self.grid.quality)
+
+        # end_time = time.time()
+
+        # if self.report:
+        #     end_quality = driver.iterations[-1].final_quality
+        #     start_quality = driver.iterations[0].initial_quality
+        #     abs_improvement = start_quality - end_quality
+        #     rel_improvement = abs_improvement / start_quality
+
+        #     print(
+        #         f"Overall improvement: {start_quality:.3e} > {end_quality:.3e}"
+        #         f"({abs_improvement:.3e}, {rel_improvement*100:.0f}%)"
+        #     )
+        #     print(f"Elapsed time: {end_time - start_time:.0f}s")
+
+        handlers: List[ParamHandler] = []
+
+        for junction in self.grid.junctions:
+            if junction.clamp is not None:
+                for i in range(junction.clamp.dof):
+                    handlers.append(ParamHandler(junction, junction.clamp, i))
 
-        start_time = time.time()
+        for _ in range(10000):
+            for handler in handlers:
+                handler.begin()
 
-        while not driver.converged:
-            driver.begin_iteration(self.grid.quality)
-            self.optimize_iteration(method)
-            driver.end_iteration(self.grid.quality)
+            for handler in handlers:
+                handler.end()
 
-        end_time = time.time()
+            print(self.grid.quality)
+            self.grid.clear_cache()
 
-        if self.report:
-            end_quality = driver.iterations[-1].final_quality
-            start_quality = driver.iterations[0].initial_quality
-            abs_improvement = start_quality - end_quality
-            rel_improvement = abs_improvement / start_quality
+        import sys
 
-            print(
-                f"Overall improvement: {start_quality:.3e} > {end_quality:.3e}"
-                f"({abs_improvement:.3e}, {rel_improvement*100:.0f}%)"
-            )
-            print(f"Elapsed time: {end_time - start_time:.0f}s")
+        sys.exit()

tests/test_optimize/test_concurrent.py

@@ -0,0 +1,41 @@
+import unittest
+
+from classy_blocks.construct.flat.sketches.grid import Grid
+from classy_blocks.construct.stack import ExtrudedStack
+from classy_blocks.mesh import Mesh
+from classy_blocks.modify.clamps.curve import LineClamp
+from classy_blocks.modify.clamps.free import FreeClamp
+from classy_blocks.modify.find.geometric import GeometricFinder
+from classy_blocks.modify.optimizer import Optimizer
+
+
+class ConcurrentOptimizerTests(unittest.TestCase):
+    def setUp(self):
+        self.mesh = Mesh()
+
+        grid = Grid([0, 0, 0], [3, 3, 0], 3, 3)
+        stack = ExtrudedStack(grid, 3, 3)
+
+        self.mesh.add(stack)
+        self.mesh.assemble()
+
+        self.optimizer = Optimizer(self.mesh)
+        self.finder = GeometricFinder(self.mesh)
+
+    def get_vertex(self, position):
+        return list(self.finder.find_in_sphere(position))[0]
+
+    def test_slopes(self):
+        free_vertex = self.get_vertex([1, 1, 1])
+        free_clamp = FreeClamp(free_vertex)
+        self.optimizer.release_vertex(free_clamp)
+
+        free_vertex.translate([0.2, -0.2, 0])
+
+        line_vertex = self.get_vertex([1, 0, 0])
+        line_clamp = LineClamp(line_vertex, [0, 0, 0], [3, 0, 0])
+        self.optimizer.release_vertex(line_clamp)
+
+        line_vertex.translate([-0.5, 0, 0])
+
+        self.optimizer.grid.get_slopes()