From cfa906e62b16d00ba4c4709cba877796e91c858a Mon Sep 17 00:00:00 2001 From: Matthew Keeler Date: Thu, 3 Sep 2026 19:23:24 +0200 Subject: [PATCH] feat(beams2d): record the optimization trajectory on each OptiStep beams2d recorded only the objective value and the design at each step. The sensitivity there, the move the optimizer made, and the objective that move bought were all computed and discarded, so recovering any of them afterwards meant re-running the solve that had already produced them. Each step now carries what the optimizer already knew: x the design the step was evaluated at x_sensitivities the filtered objective sensitivity there x_update the move taken from that design obj_values_update the objective change that move produced Nothing new is computed; every value already existed in the loop. The step is recorded after `inner_opt` rather than before the sensitivity block, since the move is not known until then, and the design is taken beforehand because the overhang filter rebinds `xPrint`. Three choices worth stating. `x_sensitivities` holds the objective sensitivity by itself, shaped like the design it belongs to, so there is one value per design variable. The volume sensitivity `dv` is not stacked alongside it. thermoelastic2d and beams3d do stack theirs, but consumers flatten this field, so an extra channel doubles its length with nothing recording that it did; photonics2d, the other 2D problem reporting sensitivities, reports the objective gradient alone. `x_update` is measured in printed density, the space the recorded design is in, so `design + x_update` is the next step's design. `inner_opt` also returns the raw density field, and differencing that would mix the two spaces. The last step's `obj_values_update` stays None. thermoelastic2d fills its own in by running an extra iteration after convergence and reverting the design, which buys one scalar for the price of a full solve; the same value is derivable from the next step for every step that has one, and the revert would put beams2d's returned design at risk. `ExtendedOptiStep` and its `design` field are untouched; `x` holds the same array, so callers reading either keep working. Optimizer behavior is unchanged: step count, objective values, per-step designs and the returned design are identical to before, with and without the overhang constraint. --- engibench/problems/beams2d/v0.py | 43 ++++++++++++++--- engibench/problems/beams2d/v1.py | 43 ++++++++++++++--- tests/test_beams2d.py | 79 ++++++++++++++++++++++++++++++++ 3 files changed, 153 insertions(+), 12 deletions(-) create mode 100644 tests/test_beams2d.py diff --git a/engibench/problems/beams2d/v0.py b/engibench/problems/beams2d/v0.py index aeb528c7..1c9d3a08 100644 --- a/engibench/problems/beams2d/v0.py +++ b/engibench/problems/beams2d/v0.py @@ -226,12 +226,9 @@ def optimize( self.reset_called = True # override for multiple reset calls in optimize c = self.simulate(xPrint, ce=ce, config=dataclasses.asdict(simulate_config)) - # Record the current state in optisteps_history - current_step = ExtendedOptiStep(obj_values=np.array(c), step=loop) - current_step.design = np.array(xPrint) - optisteps_history.append(current_step) - - loop += 1 + # The design this step was evaluated at, taken before the overhang + # filter below rebinds xPrint to the next one. + design = np.array(xPrint) dc = (-base_config.penal * xPrint ** (base_config.penal - 1) * (self.__st.Emax - self.__st.Emin)) * ce dv = np.ones(base_config.nely * base_config.nelx) @@ -252,6 +249,40 @@ def optimize( xnew.reshape(base_config.nelx * base_config.nely, 1) - x.reshape(base_config.nelx * base_config.nely, 1), np.inf, ) + + # Record the current state in optisteps_history, now that the move + # this design led to is known. + # + # x_sensitivities is the filtered objective sensitivity alone, one + # value per design variable, so that it has the same shape as the + # design it belongs to. The volume sensitivity dv is deliberately not + # stacked alongside it: consumers flatten this field, so an extra + # channel doubles its length with nothing recording that it did, and + # photonics2d -- the other 2D problem that reports sensitivities -- + # reports the objective gradient by itself. + # + # The move is measured in printed density, the space the recorded + # design is in, so that design + update is the next step's design; + # inner_opt also returns the raw density field, and differencing that + # instead would mix the two spaces. The objective delta needs the + # *next* step's objective, so it is filled in on the following pass, + # and the last step keeps None rather than costing an extra solve. + obj_values = np.array(c) + if optisteps_history: + previous = optisteps_history[-1] + previous.obj_values_update = obj_values - previous.obj_values + current_step = ExtendedOptiStep( + obj_values=obj_values, + step=loop, + x=design, + x_sensitivities=dc.copy(), + x_update=xPrint - design, + ) + current_step.design = design + optisteps_history.append(current_step) + + loop += 1 + x = deepcopy(xnew) return design_to_image(xPrint, base_config.nelx, base_config.nely), optisteps_history diff --git a/engibench/problems/beams2d/v1.py b/engibench/problems/beams2d/v1.py index 0961b524..988d94ec 100644 --- a/engibench/problems/beams2d/v1.py +++ b/engibench/problems/beams2d/v1.py @@ -76,12 +76,9 @@ def optimize( self.reset_called = True # override for multiple reset calls in optimize c = self.simulate(xPrint, ce=ce, config=dataclasses.asdict(simulate_config)) - # Record the current state in optisteps_history - current_step = ExtendedOptiStep(obj_values=np.array(c), step=loop) - current_step.design = np.array(xPrint) - optisteps_history.append(current_step) - - loop += 1 + # The design this step was evaluated at, taken before the overhang + # filter below rebinds xPrint to the next one. + design = np.array(xPrint) dc = (-base_config.penal * xPrint ** (base_config.penal - 1) * (self.__st.Emax - self.__st.Emin)) * ce dv = np.ones(base_config.nely * base_config.nelx) @@ -102,6 +99,40 @@ def optimize( xnew.reshape(base_config.nelx * base_config.nely, 1) - x.reshape(base_config.nelx * base_config.nely, 1), np.inf, ) + + # Record the current state in optisteps_history, now that the move + # this design led to is known. + # + # x_sensitivities is the filtered objective sensitivity alone, one + # value per design variable, so that it has the same shape as the + # design it belongs to. The volume sensitivity dv is deliberately not + # stacked alongside it: consumers flatten this field, so an extra + # channel doubles its length with nothing recording that it did, and + # photonics2d -- the other 2D problem that reports sensitivities -- + # reports the objective gradient by itself. + # + # The move is measured in printed density, the space the recorded + # design is in, so that design + update is the next step's design; + # inner_opt also returns the raw density field, and differencing that + # instead would mix the two spaces. The objective delta needs the + # *next* step's objective, so it is filled in on the following pass, + # and the last step keeps None rather than costing an extra solve. + obj_values = np.array(c) + if optisteps_history: + previous = optisteps_history[-1] + previous.obj_values_update = obj_values - previous.obj_values + current_step = ExtendedOptiStep( + obj_values=obj_values, + step=loop, + x=design, + x_sensitivities=dc.copy(), + x_update=xPrint - design, + ) + current_step.design = design + optisteps_history.append(current_step) + + loop += 1 + x = deepcopy(xnew) return design_to_image(xPrint, base_config.nelx, base_config.nely), optisteps_history diff --git a/tests/test_beams2d.py b/tests/test_beams2d.py new file mode 100644 index 00000000..5fb56486 --- /dev/null +++ b/tests/test_beams2d.py @@ -0,0 +1,79 @@ +"""Tests for the Beams2D problem. + +`optimize` records, for every step, the design it was evaluated at, the filtered +sensitivities there, and the move the optimizer made from it. These tests pin each +field to the step it belongs to. + +A small grid and few iterations keep this fast. +""" + +from itertools import pairwise + +import numpy as np +import pytest + +from engibench.problems.beams2d import Beams2D + +NELX = 20 +NELY = 10 +MAX_ITER = 4 +VOLFRAC = 0.5 +N_ELEMS = NELX * NELY + + +@pytest.fixture(scope="module") +def problem() -> Beams2D: + return Beams2D(seed=0, config={"nelx": NELX, "nely": NELY, "volfrac": VOLFRAC}) + + +@pytest.fixture(scope="module") +def optimized(problem: Beams2D) -> tuple[np.ndarray, list]: + """Run optimize once and reuse across tests (the expensive step).""" + problem.reset(seed=0) + start = np.full(problem.design_space.shape, VOLFRAC, dtype=np.float64) + # max_iter goes through optimize rather than the constructor: it lives on + # Config, not SimulateConfig, so optimize rebuilds it from the default. + return problem.optimize(start, config={"max_iter": MAX_ITER}) + + +def test_history_is_numbered_from_zero(optimized: tuple) -> None: + _opt, history = optimized + assert 1 <= len(history) <= MAX_ITER + assert [step.step for step in history] == list(range(len(history))) + + +def test_every_step_records_design_sensitivities_and_update(optimized: tuple) -> None: + _opt, history = optimized + for step in history: + assert step.design.shape == (N_ELEMS,) + # `x` is the core field every generic consumer reads; `design` is kept + # for the callers that already use it. + assert step.x is step.design + # One sensitivity per design variable, so it is shaped like the design. + assert step.x_sensitivities.shape == step.x.shape + assert step.x_update.shape == (N_ELEMS,) + + +def test_sensitivities_are_nonpositive(optimized: tuple) -> None: + """Compliance falls as material is added, and the optimizer clips dc at zero.""" + _opt, history = optimized + for step in history: + assert np.all(step.x_sensitivities <= 0.0) + assert np.all(np.isfinite(step.x_sensitivities)) + + +def test_update_is_the_move_to_the_next_design(optimized: tuple) -> None: + """x_update is the step the optimizer took, so it lands on the next design.""" + _opt, history = optimized + for step, following in pairwise(history): + np.testing.assert_allclose(step.x + step.x_update, following.x, atol=1e-12) + + +def test_objective_delta_is_filled_in_except_on_the_last_step(optimized: tuple) -> None: + """A step's delta needs the next step's objective, so the last one has none.""" + _opt, history = optimized + if len(history) < 2: # noqa: PLR2004 - a single step has no delta to check + pytest.skip("optimization converged in one step") + for step, following in pairwise(history): + np.testing.assert_allclose(step.obj_values_update, following.obj_values - step.obj_values) + assert history[-1].obj_values_update is None