From b33fe3e717c0ef9b0bbad880dec29bdd46364be0 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 27 Aug 2026 14:44:18 -0700 Subject: [PATCH 1/2] Document supported features in docs/features.md Adds a running Julia-vs-Python feature list, compiled from the current sources, covering the filters, SSP1/SSP2 projections, dilation/erosion, lengthscale constraints, AD support, and the low-level Julia API, plus known limitations and gaps. Also links the new page from the top-level and Julia READMEs and adds Julia installation instructions to the top-level README. --- README.md | 19 ++++++++++++++++++ docs/features.md | 43 +++++++++++++++++++++++++++++++++++++++++ src/julia/SSP/README.md | 3 +++ 3 files changed, 65 insertions(+) create mode 100644 docs/features.md diff --git a/README.md b/README.md index e259682..f61f7af 100644 --- a/README.md +++ b/README.md @@ -6,8 +6,15 @@ This is a repository for code implementing the **smoothed subpixel projection (S * G. Romano, R. Arrieta, and S. G. Johnson, [“Differentiating through binarized topology changes: Second-order subpixel-smoothed projection,”](http://arxiv.org/abs/2601.10737) arXiv.org e-Print archive, 2601.10737, January 2026. * R. Arrieta, G. Romano, and S. G. Johnson, [“Hyperparameter-free minimum-lengthscale constraints for topology optimization,”](http://arxiv.org/abs/2507.16108) arXiv.org e-Print archive, 2507.16108, July 2025. +## Documentation + +* [Supported features](docs/features.md) — a running list of what the Julia and Python + packages each implement, along with known limitations and gaps. + ## Installation +### Python + Install the PyPI distribution: ```bash @@ -25,3 +32,15 @@ For local development: ```bash python -m pip install -e ".[dev]" ``` + +### Julia + +The Julia package is not registered yet, so install it from this repository: + +```julia +using Pkg +Pkg.develop(path="src/julia/SSP") +``` + +See [`src/julia/SSP/README.md`](src/julia/SSP/README.md) for usage of both the high-level +and low-level Julia APIs. diff --git a/docs/features.md b/docs/features.md new file mode 100644 index 0000000..a8c2266 --- /dev/null +++ b/docs/features.md @@ -0,0 +1,43 @@ +# Supported Features + +The Julia package lives in [`src/julia/SSP`](../src/julia/SSP) and the Python package in +[`src/python/ssp_topopt`](../src/python/ssp_topopt). The table below tracks what each one +currently implements; **please keep it up to date when adding or removing functionality.** + +| Feature | Julia (`SSP`) | Python (`ssp_topopt`) | +| --- | --- | --- | +| Conic ("hat") filter | ✅ `conic_filter` | ✅ `conic_filter` | +| Filter radius from an eroded threshold point | ❌ | ✅ `get_conic_radius_from_eta_e` | +| Plain tanh projection | ❌ (internal only) | ✅ `tanh_projection` | +| First-order subpixel smoothing (SSP1), linear interpolation | ✅ `ssp1_linear` | ✅ `ssp1_bilinear` | +| First-order subpixel smoothing (SSP1), cubic interpolation | ✅ `ssp1` | ❌ | +| Second-order subpixel smoothing (SSP2), differentiable through topology changes | ✅ `ssp2` | ✅ `ssp2` | +| Finite and infinite projection strength (0 ≤ β ≤ ∞) | ✅ | ✅ | +| Dilation/erosion of the projected contour | ✅ `dilation_distance` argument | ❌ | +| Minimum-lengthscale constraints for solid and void | ✅ `constraint_solid`, `constraint_void` | ❌ | +| Lengthscale constraints compatible with any SSP order | ✅ (constraints act on `rho_filtered`/`rho_projected`) | ❌ | +| Reverse-mode automatic differentiation | ✅ hand-written adjoints, exposed to Zygote.jl and friends through a ChainRulesCore.jl extension | ✅ through JAX (`grad`, `jit`, `vmap`) | +| Dimensionality | N-dimensional code paths (only 2D is currently tested) | 2D only | +| Periodic filter axes | ❌ | ✅ `periodic_axes` argument of `conic_filter` | +| Low-level `init`/`solve!`/`adjoint_solve!` API with reduced allocations | ✅ | ❌ | +| Explicit control over padding/boundary conditions, kernels, interpolation, and projection target points | ✅ (low-level API) | ❌ | + +## Known limitations + +* The subpixel fill factor is the analytic expression for a *circular* smoothing kernel, so + both implementations assume an isotropic grid (`dx == dy`). Julia asserts that all grid + steps are equal; Python takes a single scalar `resolution` in the projection routines + (`conic_filter` does accept an anisotropic `resolution`). +* Only 2D usage is covered by the tests and examples in this repository, even though the + Julia routines are written generically over the number of dimensions. +* The high-level Julia `conic_filter` always pads by replicating the boundary values. + Other padding styles (`FillPadding`, `Inner`) are only reachable through the low-level API. + +## Not yet supported + +Contributions welcome — these are known gaps rather than fundamental limitations: + +* Python: cubic-interpolation SSP1, dilation/erosion, and minimum-lengthscale constraints. +* Python: a low-level API with reusable workspaces. +* Julia: `get_conic_radius_from_eta_e`-style helpers and periodic filter axes. +* Both: validated 3D usage and anisotropic grid spacings in the projection. diff --git a/src/julia/SSP/README.md b/src/julia/SSP/README.md index 2c40f38..e747cf9 100644 --- a/src/julia/SSP/README.md +++ b/src/julia/SSP/README.md @@ -3,6 +3,9 @@ A Smoothed Subpixel Projection (SSP) package for topology optimization in Julia. Supports N-dimensional data and reverse-mode automatic differentiation with minimal allocations. +For a feature-by-feature comparison of this package with its Python cousin, see +[`docs/features.md`](../../../docs/features.md). + ## Usage This package provides a high-level API nearly identical to its python cousin. From 44dd5fbd9266afa03118a7fb075b4f75d56ebbf4 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 27 Aug 2026 15:04:09 -0700 Subject: [PATCH 2/2] Add minimum-lengthscale constraints to the Python package Ports the geometric lengthscale constraints from the Julia package: constraint_solid and constraint_void, with the hyperparameter-free thresholds of Arrieta et al. (arXiv:2507.16108). Both are written in terms of rho_filtered and rho_projected only, so they compose with any projection order (ssp1_bilinear, ssp2, or a plain tanh projection). Includes a test suite covering the threshold functions, lengthscale detection on stripe/gap geometries, gradients against finite differences, composition with each projection, jit, and argument validation, plus a worked example that sweeps the constraint over feature widths and runs a two-stage constrained optimization. --- README.md | 1 + docs/features.md | 6 +- examples/python/ssp_constrained_example.py | 269 ++++++++++++++++ src/python/ssp_topopt/__init__.py | 3 + src/python/ssp_topopt/constraints.py | 271 ++++++++++++++++ tests/python/test_constraints.py | 352 +++++++++++++++++++++ 6 files changed, 899 insertions(+), 3 deletions(-) create mode 100644 examples/python/ssp_constrained_example.py create mode 100644 src/python/ssp_topopt/constraints.py create mode 100644 tests/python/test_constraints.py diff --git a/README.md b/README.md index f61f7af..f5dddb1 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,7 @@ The Python import package is `ssp_topopt`: ```python from ssp_topopt import conic_filter, get_conic_radius_from_eta_e, ssp1_bilinear,ssp2 +from ssp_topopt import constraint_solid, constraint_void ``` For local development: diff --git a/docs/features.md b/docs/features.md index a8c2266..c93a153 100644 --- a/docs/features.md +++ b/docs/features.md @@ -14,8 +14,8 @@ currently implements; **please keep it up to date when adding or removing functi | Second-order subpixel smoothing (SSP2), differentiable through topology changes | ✅ `ssp2` | ✅ `ssp2` | | Finite and infinite projection strength (0 ≤ β ≤ ∞) | ✅ | ✅ | | Dilation/erosion of the projected contour | ✅ `dilation_distance` argument | ❌ | -| Minimum-lengthscale constraints for solid and void | ✅ `constraint_solid`, `constraint_void` | ❌ | -| Lengthscale constraints compatible with any SSP order | ✅ (constraints act on `rho_filtered`/`rho_projected`) | ❌ | +| Minimum-lengthscale constraints for solid and void | ✅ `constraint_solid`, `constraint_void` | ✅ `constraint_solid`, `constraint_void` | +| Lengthscale constraints compatible with any SSP order | ✅ (constraints act on `rho_filtered`/`rho_projected`) | ✅ (constraints act on `rho_filtered`/`rho_projected`) | | Reverse-mode automatic differentiation | ✅ hand-written adjoints, exposed to Zygote.jl and friends through a ChainRulesCore.jl extension | ✅ through JAX (`grad`, `jit`, `vmap`) | | Dimensionality | N-dimensional code paths (only 2D is currently tested) | 2D only | | Periodic filter axes | ❌ | ✅ `periodic_axes` argument of `conic_filter` | @@ -37,7 +37,7 @@ currently implements; **please keep it up to date when adding or removing functi Contributions welcome — these are known gaps rather than fundamental limitations: -* Python: cubic-interpolation SSP1, dilation/erosion, and minimum-lengthscale constraints. +* Python: cubic-interpolation SSP1 and dilation/erosion. * Python: a low-level API with reusable workspaces. * Julia: `get_conic_radius_from_eta_e`-style helpers and periodic filter axes. * Both: validated 3D usage and anisotropic grid spacings in the projection. diff --git a/examples/python/ssp_constrained_example.py b/examples/python/ssp_constrained_example.py new file mode 100644 index 0000000..075a82a --- /dev/null +++ b/examples/python/ssp_constrained_example.py @@ -0,0 +1,269 @@ +"""Minimum-lengthscale (geometric) constraints combined with SSP. + +The constraints of `ssp_topopt.constraints` measure whether the solid or void +features of a design are thinner than a target lengthscale. They are formulated +in terms of the filtered and projected densities only, so they work with any +order of subpixel smoothing; this example uses `ssp2`. + +The example has two parts: + +1. A sweep over a simple stripe geometry, which shows that the solid constraint + changes sign at the target lengthscale (compare fig. 1c of Arrieta et al.). +2. A two-stage optimization that matches a target pattern containing bars that + are thinner than the target lengthscale. The first stage is unconstrained and + happily reproduces the sub-lengthscale bars; the second stage turns the + constraints on and repairs the design, at a small cost in the figure of merit. + +Ref: R. Arrieta, G. Romano, and S. G. Johnson, "Hyperparameter-free +minimum-lengthscale constraints for topology optimization," arXiv.org e-Print +archive, 2507.16108, July 2025. +""" + +import time + +import nlopt +import numpy as np +from jax import jit, value_and_grad +from jax import numpy as jnp +from matplotlib import pyplot as plt + +from ssp_topopt import conic_filter, constraint_solid, constraint_void, ssp2 + +BETA = np.inf +ETA_I = 0.5 + + +def filter_and_project(rho, filter_radius, lx, ly, resolution): + """The usual filter-then-project pipeline.""" + rho_filtered = conic_filter(rho, filter_radius, lx, ly, resolution) + rho_projected = ssp2(rho_filtered, BETA, ETA_I, resolution) + return rho_filtered, rho_projected + + +def feature_width(rho_projected, resolution, solid=True): + """Width of the central feature of a stripe design, in physical units.""" + profile = np.asarray(rho_projected)[:, rho_projected.shape[1] // 2] + mask = profile > 0.5 if solid else profile < 0.5 + return float(np.count_nonzero(mask)) / resolution + + +def constraint_sweep(): + """Evaluate the constraints for stripes of varying width.""" + lx = ly = 0.5 + resolution = 100 + target_length = 0.1 + filter_radius = target_length + + nx = int(np.round(lx * resolution)) + 1 + ny = int(np.round(ly * resolution)) + 1 + coords = np.linspace(-lx / 2, lx / 2, nx) + x = coords[:, None] * np.ones((1, ny)) + + latent_widths = np.linspace(0.04, 0.24, 21) + widths = [] + solid_values = [] + void_values = [] + + # Note that the void constraint also fires for very thin stripes: once the + # stripe is thin enough that it barely projects to solid, the crest of the + # filtered density is a nearly-solid feature sitting inside the void region. + print("Constraint sweep over stripe widths " f"(target lengthscale {target_length})") + for latent_width in latent_widths: + # A binary stripe of width `latent_width`, which the filter and the + # projection turn into a stripe of some (smaller) physical width. + rho = jnp.asarray((np.abs(x) <= latent_width / 2).astype(float)) + rho_filtered, rho_projected = filter_and_project( + rho, filter_radius, lx, ly, resolution + ) + + width = feature_width(rho_projected, resolution) + solid_value = float( + constraint_solid(rho_filtered, rho_projected, resolution, target_length) + ) + void_value = float( + constraint_void(rho_filtered, rho_projected, resolution, target_length) + ) + + widths.append(width) + solid_values.append(solid_value) + void_values.append(void_value) + print( + f" latent width={latent_width:.3f} physical width={width:.3f} " + f"solid={solid_value:+.3e} void={void_value:+.3e}" + ) + + # Feasible designs have an exactly vanishing violation, so clip to a floor to + # keep them visible on a log scale. + floor = 1e-2 + plt.figure(figsize=(5, 3.5)) + plt.axhspan(floor / 2, 1.0, color="tab:green", alpha=0.1, label="feasible") + plt.plot( + widths, np.maximum(np.asarray(solid_values) + 1, floor), marker="o", label="solid" + ) + plt.plot( + widths, np.maximum(np.asarray(void_values) + 1, floor), marker="s", label="void" + ) + plt.axhline(1.0, color="k", linestyle=":", label="threshold") + plt.axvline(target_length, color="r", linestyle="--", label="target lengthscale") + plt.yscale("log") + plt.ylim(bottom=floor / 2) + plt.xlabel("physical width of the stripe") + plt.ylabel("constraint / threshold") + plt.legend(fontsize=8) + plt.title("Constraint value vs feature width") + plt.tight_layout() + plt.savefig("constraint_sweep.png") + + +def optimization_demo(): + """Repair a design with sub-lengthscale features using the constraints.""" + lx, ly = 1.2, 0.6 + resolution = 50 + target_length = 0.12 + filter_radius = target_length + + nx = int(np.round(lx * resolution)) + 1 + ny = int(np.round(ly * resolution)) + 1 + num_vars = nx * ny + coords = np.linspace(-lx / 2, lx / 2, nx) + + # A target pattern of vertical bars, two of which are thinner than the + # target lengthscale and therefore cannot be manufactured. + target = np.zeros((nx, ny)) + position = -lx / 2 + 0.08 + for width in (0.24, 0.16, 0.08, 0.04): + target[(coords >= position) & (coords <= position + width), :] = 1.0 + position += width + 0.12 + target_jnp = jnp.asarray(target) + + def pipeline(rho_flat): + return filter_and_project( + rho_flat.reshape((nx, ny)), filter_radius, lx, ly, resolution + ) + + def figure_of_merit(rho_flat): + _, rho_projected = pipeline(rho_flat) + return jnp.mean((rho_projected - target_jnp) ** 2) + + def solid_constraint(rho_flat): + rho_filtered, rho_projected = pipeline(rho_flat) + return constraint_solid( + rho_filtered, rho_projected, resolution, target_length + ) + + def void_constraint(rho_flat): + rho_filtered, rho_projected = pipeline(rho_flat) + return constraint_void(rho_filtered, rho_projected, resolution, target_length) + + objective_and_grad = jit(value_and_grad(figure_of_merit)) + constraints_and_grad = { + "solid": jit(value_and_grad(solid_constraint)), + "void": jit(value_and_grad(void_constraint)), + } + + fom_history = [] + constraint_history = {"solid": [], "void": []} + + def run_stage(x_init, constrained, maxeval): + def nlopt_objective(x, grad_out): + value, gradient = objective_and_grad(jnp.asarray(x)) + if grad_out.size > 0: + grad_out[:] = np.asarray(gradient, dtype=float) + fom_history.append(float(value)) + for name, value_and_grad_fn in constraints_and_grad.items(): + constraint_history[name].append( + float(value_and_grad_fn(jnp.asarray(x))[0]) + ) + return float(value) + + opt = nlopt.opt(nlopt.LD_CCSAQ, num_vars) + opt.set_lower_bounds(np.zeros(num_vars)) + opt.set_upper_bounds(np.ones(num_vars)) + opt.set_min_objective(nlopt_objective) + + if constrained: + for value_and_grad_fn in constraints_and_grad.values(): + + def nlopt_constraint(x, grad_out, fn=value_and_grad_fn): + value, gradient = fn(jnp.asarray(x)) + if grad_out.size > 0: + grad_out[:] = np.asarray(gradient, dtype=float) + return float(value) + + # The constraints are normalized, so the feasible region is + # exactly where they are nonpositive. + opt.add_inequality_constraint(nlopt_constraint, 0.0) + + opt.set_maxeval(maxeval) + start = time.perf_counter() + x_opt = opt.optimize(np.asarray(x_init, dtype=float)) + elapsed = time.perf_counter() - start + + label = "constrained" if constrained else "unconstrained" + print( + f" {label}: {maxeval} evaluations in {elapsed:6.1f}s " + f"FOM={fom_history[-1]:.4e} " + f"solid={constraint_history['solid'][-1]:+.3e} " + f"void={constraint_history['void'][-1]:+.3e}" + ) + return x_opt + + print(f"\nTwo-stage optimization (target lengthscale {target_length})") + rng = np.random.default_rng(0) + x_init = 0.5 * np.ones(num_vars) + 0.01 * rng.standard_normal(num_vars) + + stage1_evals = 120 + x_stage1 = run_stage(x_init, constrained=False, maxeval=stage1_evals) + # The constrained stage needs a few hundred iterations: CCSA has to first + # drag the design back into the feasible region and then re-minimize. + x_stage2 = run_stage(x_stage1, constrained=True, maxeval=400) + + _, projected_stage1 = pipeline(jnp.asarray(x_stage1)) + _, projected_stage2 = pipeline(jnp.asarray(x_stage2)) + + plt.figure(figsize=(9, 3)) + for index, (image, title) in enumerate( + ( + (target, "target"), + (np.asarray(projected_stage1), "stage 1: unconstrained"), + (np.asarray(projected_stage2), "stage 2: lengthscale constrained"), + ) + ): + plt.subplot(1, 3, index + 1) + plt.imshow(image.T, vmin=0, vmax=1, cmap="binary", origin="lower") + plt.title(title, fontsize=9) + plt.axis("off") + plt.tight_layout() + plt.savefig("constraint_designs.png") + + plt.figure(figsize=(8, 3)) + plt.subplot(1, 2, 1) + plt.semilogy(fom_history, linewidth=1.5) + plt.axvline(stage1_evals, color="k", linestyle="--", linewidth=1) + plt.xlabel("iteration") + plt.ylabel("FOM") + plt.title("Objective history", fontsize=9) + + plt.subplot(1, 2, 2) + for name, values in constraint_history.items(): + # Shift by one so that the (normalized) constraint can be shown on a log + # scale: values below one are feasible. + plt.semilogy(np.asarray(values) + 1, linewidth=1.5, label=name) + plt.axhline(1.0, color="k", linestyle=":", label="threshold") + plt.axvline(stage1_evals, color="k", linestyle="--", linewidth=1) + plt.xlabel("iteration") + plt.ylabel("constraint / threshold") + plt.title("Constraint history", fontsize=9) + plt.legend(fontsize=8) + plt.tight_layout() + plt.savefig("constraint_history.png") + plt.show() + + +def main(): + constraint_sweep() + optimization_demo() + + +if __name__ == "__main__": + main() diff --git a/src/python/ssp_topopt/__init__.py b/src/python/ssp_topopt/__init__.py index e58537f..d0f57f7 100644 --- a/src/python/ssp_topopt/__init__.py +++ b/src/python/ssp_topopt/__init__.py @@ -1,5 +1,6 @@ """Public Python API for smoothed subpixel projection (SSP) for topology optimization.""" +from .constraints import constraint_solid, constraint_void from .core import ssp1_bilinear,ssp2 from .utils import conic_filter, get_conic_radius_from_eta_e, tanh_projection @@ -7,6 +8,8 @@ "ssp1_bilinear", "ssp2", "conic_filter", + "constraint_solid", + "constraint_void", "get_conic_radius_from_eta_e", "tanh_projection", ] diff --git a/src/python/ssp_topopt/constraints.py b/src/python/ssp_topopt/constraints.py new file mode 100644 index 0000000..8a98757 --- /dev/null +++ b/src/python/ssp_topopt/constraints.py @@ -0,0 +1,271 @@ +"""Minimum-lengthscale (geometric) constraints for topology optimization. + +These are the geometric constraints of Zhou et al. (2015) combined with the +hyperparameter-free thresholds derived by Arrieta et al. (2025). The constraints +penalize solid (or void) features whose lengthscale falls below a target value, +and are formulated purely in terms of the filtered density and the projected +density. Consequently they work with *any* order of subpixel smoothing -- +`ssp1_bilinear`, `ssp2`, or even a plain `tanh_projection` -- since the +projection only enters through `rho_projected`. + +The solid constraint reads + + g_s = (1/N) Σ_i I_s,i [min(rho_filtered_i - eta_e, 0)]^2 , + I_s = rho_projected * exp(-c |∇ rho_filtered|^2) , + +and the void constraint is the complementary expression + + g_v = (1/N) Σ_i I_v,i [min(eta_d - rho_filtered_i, 0)]^2 , + I_v = (1 - rho_projected) * exp(-c |∇ rho_filtered|^2) , + +where the structural functions I_s and I_v single out the "inflection regions" +of the design (the interior of a feature, where the filtered density is +stationary), and eta_e, eta_d are the eroded/dilated threshold points of the +conic filter. Following Arrieta et al., the decay rate is c = 64 R^2 and the +constraint threshold is eps = 1e-8, where R is the conic filter radius; the +constraint is well behaved for target_length / R roughly in [0.25, 1.5]. + +Refs: + +R. Arrieta, G. Romano, and S. G. Johnson, "Hyperparameter-free minimum-lengthscale +constraints for topology optimization," arXiv.org e-Print archive, 2507.16108, +July 2025. + +M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, "Minimum length scale in topology +optimization by geometric constraints," Computer Methods in Applied Mechanics and +Engineering, vol. 293, pp. 266-282, 2015. + +X. Qian and O. Sigmund, "Topological design of electromechanical actuators with +robustness toward over- and under-etching," Computer Methods in Applied Mechanics +and Engineering, vol. 253, pp. 237-251, 2013. +""" + +from typing import Optional + +from jax import numpy as jnp + +from .utils import ArrayLikeType, gradient + +# Hyperparameters derived in section 4.1 of Arrieta et al. (2025). The decay rate +# is expressed relative to the square of the conic filter radius, i.e. the actual +# decay rate is c = DEFAULT_CONSTRAINT_DECAYRATE * conic_radius**2. +DEFAULT_CONSTRAINT_THRESHOLD = 1e-8 +DEFAULT_CONSTRAINT_DECAYRATE = 64.0 + + +def solid_threshold(lengthscale_ratio: float): + """The eroded threshold point eta_e of a conic filter. + + Ref: Eq. (9) of Arrieta et al. (2025), originally from Qian and Sigmund (2013). + + Args: + lengthscale_ratio: the ratio of the target lengthscale to the conic + filter radius, which must be nonnegative. + + Returns: + The threshold point in the range [1/2, 1]. + """ + x = jnp.asarray(lengthscale_ratio) + return jnp.where( + x < 1, + x**2 / 4 + 1 / 2, + jnp.where(x < 2, -(x**2) / 4 + x, 1.0), + ) + + +def void_threshold(lengthscale_ratio: float): + """The dilated threshold point eta_d of a conic filter. + + Ref: Eq. (12) of Arrieta et al. (2025), originally from Qian and Sigmund (2013). + + Args: + lengthscale_ratio: the ratio of the target lengthscale to the conic + filter radius, which must be nonnegative. + + Returns: + The threshold point in the range [0, 1/2]. + """ + x = jnp.asarray(lengthscale_ratio) + return jnp.where( + x < 1, + 1 / 2 - x**2 / 4, + jnp.where(x < 2, 1 + x**2 / 4 - x, 0.0), + ) + + +def _geometric_constraint( + rho_filtered: ArrayLikeType, + rho_projected: ArrayLikeType, + resolution: float, + target_length: float, + conic_radius: Optional[float], + constraint_threshold: float, + constraint_decayrate: float, + solid: bool, +): + """Shared implementation of the solid and void lengthscale constraints.""" + if target_length < 0: + raise ValueError("The target lengthscale must be nonnegative.") + + if conic_radius is None: + conic_radius = target_length + if conic_radius <= 0: + raise ValueError("The conic filter radius must be positive.") + + rho_filtered = jnp.asarray(rho_filtered) + rho_projected = jnp.asarray(rho_projected) + if rho_filtered.ndim != 2: + raise ValueError( + f"Only 2D designs are supported, got {rho_filtered.ndim} dimensions." + ) + if rho_filtered.shape != rho_projected.shape: + raise ValueError( + "rho_filtered and rho_projected must have the same shape, got " + f"{rho_filtered.shape} and {rho_projected.shape}." + ) + + decayrate = constraint_decayrate * conic_radius**2 + + # The gradient of the filtered density vanishes in the interior of a feature, + # so this term restricts the constraint to those "inflection regions". + rho_filtered_grad = gradient(rho_filtered, resolution) + rho_filtered_grad_normsq = jnp.sum(rho_filtered_grad**2, axis=-1) + extremal_region = jnp.exp(-decayrate * rho_filtered_grad_normsq) + + if solid: + eta_m = solid_threshold(target_length / conic_radius) + inflection_region = rho_projected * extremal_region + # Only densities below the eroded threshold, i.e. features that are too + # thin to survive an erosion by the target lengthscale, are penalized. + beyond_threshold = jnp.minimum(rho_filtered - eta_m, 0.0) + else: + eta_m = void_threshold(target_length / conic_radius) + inflection_region = (1 - rho_projected) * extremal_region + beyond_threshold = jnp.minimum(eta_m - rho_filtered, 0.0) + + violation = jnp.mean(inflection_region * beyond_threshold**2) + + # Normalize so that the constraint is satisfied when it is nonpositive. + return violation / constraint_threshold - 1 + + +def constraint_solid( + rho_filtered: ArrayLikeType, + rho_projected: ArrayLikeType, + resolution: float, + target_length: float, + conic_radius: Optional[float] = None, + constraint_threshold: float = DEFAULT_CONSTRAINT_THRESHOLD, + constraint_decayrate: float = DEFAULT_CONSTRAINT_DECAYRATE, +): + """Calculate a solid minimum-lengthscale constraint function. + + This technique takes smoothed data, e.g. from filtering, `rho_filtered` and + binary data, e.g. from projection, `rho_projected` both defined on the same + grid, and measures whether features in the solid region, i.e. where + `rho_projected` takes values of 1, violate the minimum `target_length`. + + The returned value is normalized by the constraint threshold, so the + constraint is satisfied when the value is nonpositive and violated when it is + positive. It can therefore be handed directly to a nonlinear optimizer such + as `nlopt` as an inequality constraint. The unnormalized constraint value of + Eq. (7) of Arrieta et al. (2025) is `(value + 1) * constraint_threshold`. + + Any projection may be used to produce `rho_projected`, including + `ssp1_bilinear`, `ssp2`, and `tanh_projection`. + + Args: + rho_filtered: the (2D) filtered design parameters, e.g. from + `conic_filter`. + rho_projected: the (2D) projected design parameters, e.g. from `ssp2`. + resolution: resolution of the design grid. + target_length: the minimum lengthscale to impose on the solid region. + conic_radius: the radius of the conic filter used to obtain + `rho_filtered`. Defaults to `target_length`, which is the + recommended choice. + constraint_threshold: the threshold that separates feasible from + infeasible designs. May be tuned if feasible designs still don't + meet the target lengthscale. + constraint_decayrate: the decay rate of the structural function outside + of the inflection region, relative to `conic_radius**2`. + + Returns: + The normalized constraint value, which is nonpositive when the design + satisfies the minimum lengthscale. + + Example: + >>> rho_filtered = conic_filter(rho, filter_radius, lx, ly, resolution) + >>> rho_projected = ssp2(rho_filtered, beta, eta_i, resolution) + >>> constraint_solid(rho_filtered, rho_projected, resolution, filter_radius) + """ + return _geometric_constraint( + rho_filtered, + rho_projected, + resolution, + target_length, + conic_radius, + constraint_threshold, + constraint_decayrate, + solid=True, + ) + + +def constraint_void( + rho_filtered: ArrayLikeType, + rho_projected: ArrayLikeType, + resolution: float, + target_length: float, + conic_radius: Optional[float] = None, + constraint_threshold: float = DEFAULT_CONSTRAINT_THRESHOLD, + constraint_decayrate: float = DEFAULT_CONSTRAINT_DECAYRATE, +): + """Calculate a void minimum-lengthscale constraint function. + + This technique takes smoothed data, e.g. from filtering, `rho_filtered` and + binary data, e.g. from projection, `rho_projected` both defined on the same + grid, and measures whether features in the void region, i.e. where + `rho_projected` takes values of 0, violate the minimum `target_length`. + + The returned value is normalized by the constraint threshold, so the + constraint is satisfied when the value is nonpositive and violated when it is + positive. It can therefore be handed directly to a nonlinear optimizer such + as `nlopt` as an inequality constraint. The unnormalized constraint value of + Eq. (10) of Arrieta et al. (2025) is `(value + 1) * constraint_threshold`. + + Any projection may be used to produce `rho_projected`, including + `ssp1_bilinear`, `ssp2`, and `tanh_projection`. + + Args: + rho_filtered: the (2D) filtered design parameters, e.g. from + `conic_filter`. + rho_projected: the (2D) projected design parameters, e.g. from `ssp2`. + resolution: resolution of the design grid. + target_length: the minimum lengthscale to impose on the void region. + conic_radius: the radius of the conic filter used to obtain + `rho_filtered`. Defaults to `target_length`, which is the + recommended choice. + constraint_threshold: the threshold that separates feasible from + infeasible designs. May be tuned if feasible designs still don't + meet the target lengthscale. + constraint_decayrate: the decay rate of the structural function outside + of the inflection region, relative to `conic_radius**2`. + + Returns: + The normalized constraint value, which is nonpositive when the design + satisfies the minimum lengthscale. + + Example: + >>> rho_filtered = conic_filter(rho, filter_radius, lx, ly, resolution) + >>> rho_projected = ssp2(rho_filtered, beta, eta_i, resolution) + >>> constraint_void(rho_filtered, rho_projected, resolution, filter_radius) + """ + return _geometric_constraint( + rho_filtered, + rho_projected, + resolution, + target_length, + conic_radius, + constraint_threshold, + constraint_decayrate, + solid=False, + ) diff --git a/tests/python/test_constraints.py b/tests/python/test_constraints.py new file mode 100644 index 0000000..db5cc6d --- /dev/null +++ b/tests/python/test_constraints.py @@ -0,0 +1,352 @@ +"""Tests for the minimum-lengthscale (geometric) constraints. + +The reference behavior is the one described in R. Arrieta, G. Romano, and +S. G. Johnson, "Hyperparameter-free minimum-lengthscale constraints for topology +optimization," arXiv:2507.16108 (2025): with the derived hyperparameters, the +constraint is violated (positive) when the physical lengthscale of a feature +falls below the target lengthscale and satisfied (nonpositive) otherwise. +""" + +import unittest + +import numpy as np +from jax import grad, jit, value_and_grad +from jax import numpy as jnp +from jax.experimental import enable_x64 + +from ssp_topopt import ( + conic_filter, + constraint_solid, + constraint_void, + ssp1_bilinear, + ssp2, + tanh_projection, +) +from ssp_topopt.constraints import solid_threshold, void_threshold + + +class TestThresholdFunctions(unittest.TestCase): + """The conic-filter threshold points of Qian and Sigmund (2013).""" + + def test_known_values(self): + # Eqs. (9) and (12) of Arrieta et al. (2025). + self.assertAlmostEqual(float(solid_threshold(0.0)), 0.5) + self.assertAlmostEqual(float(solid_threshold(1.0)), 0.75) + self.assertAlmostEqual(float(solid_threshold(2.0)), 1.0) + self.assertAlmostEqual(float(solid_threshold(3.0)), 1.0) + + self.assertAlmostEqual(float(void_threshold(0.0)), 0.5) + self.assertAlmostEqual(float(void_threshold(1.0)), 0.25) + self.assertAlmostEqual(float(void_threshold(2.0)), 0.0) + self.assertAlmostEqual(float(void_threshold(3.0)), 0.0) + + def test_thresholds_are_complementary_and_monotonic(self): + ratios = np.linspace(0.0, 3.0, 61) + eta_e = np.asarray(solid_threshold(ratios)) + eta_d = np.asarray(void_threshold(ratios)) + + np.testing.assert_allclose(eta_e + eta_d, 1.0, atol=1e-6) + self.assertTrue(np.all(np.diff(eta_e) >= -1e-7)) + self.assertTrue(np.all(np.diff(eta_d) <= 1e-7)) + self.assertTrue(np.all((eta_e >= 0.5 - 1e-7) & (eta_e <= 1.0 + 1e-7))) + self.assertTrue(np.all((eta_d >= -1e-7) & (eta_d <= 0.5 + 1e-7))) + + +class LengthscaleFixture(unittest.TestCase): + """A stripe geometry whose physical lengthscale can be measured directly.""" + + def setUp(self): + self.lx = 0.5 + self.ly = 0.5 + self.resolution = 100 + self.target_length = 0.1 + self.filter_radius = self.target_length + self.beta = np.inf + self.eta_i = 0.5 + self.nx = int(np.round(self.lx * self.resolution)) + 1 + self.ny = int(np.round(self.ly * self.resolution)) + 1 + + coords = np.linspace(-self.lx / 2, self.lx / 2, self.nx) + self.x = coords[:, None] * np.ones((1, self.ny)) + + def _filter_and_project(self, rho, projection=ssp2): + rho_filtered = conic_filter( + jnp.asarray(rho), self.filter_radius, self.lx, self.ly, self.resolution + ) + rho_projected = projection( + rho_filtered, self.beta, self.eta_i, self.resolution + ) + return rho_filtered, rho_projected + + def _feature_width(self, rho_projected, solid=True): + """Width of the central feature of the projected design, in physical units.""" + profile = np.asarray(rho_projected)[:, self.ny // 2] + mask = profile > 0.5 if solid else profile < 0.5 + return float(np.count_nonzero(mask)) / self.resolution + + +class TestSolidConstraint(LengthscaleFixture): + def test_detects_thin_solid_features(self): + """A stripe thinner than the target violates the solid constraint.""" + for latent_width in (0.06, 0.08): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) <= latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected) + self.assertGreater(width, 0.0) + self.assertLess(width, self.target_length) + + value = float( + constraint_solid( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertGreater(value, 0.0) + + def test_accepts_thick_solid_features(self): + """A stripe at or above the target satisfies the solid constraint.""" + for latent_width in (0.10, 0.14, 0.20): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) <= latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected) + self.assertGreaterEqual(width, self.target_length) + + value = float( + constraint_solid( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertLessEqual(value, 0.0) + + def test_uniform_solid_design_is_feasible(self): + """A fully solid design has no interfaces and no lengthscale violation.""" + rho_filtered = jnp.ones((self.nx, self.ny)) + rho_projected = jnp.ones((self.nx, self.ny)) + + value = float( + constraint_solid( + rho_filtered, rho_projected, self.resolution, self.target_length + ) + ) + self.assertAlmostEqual(value, -1.0, places=6) + + +class TestVoidConstraint(LengthscaleFixture): + def test_detects_thin_void_features(self): + """A gap thinner than the target violates the void constraint.""" + for latent_width in (0.06, 0.08): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) > latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected, solid=False) + self.assertGreater(width, 0.0) + self.assertLess(width, self.target_length) + + value = float( + constraint_void( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertGreater(value, 0.0) + + def test_accepts_thick_void_features(self): + """A gap at or above the target satisfies the void constraint.""" + for latent_width in (0.10, 0.14, 0.20): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) > latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected, solid=False) + self.assertGreaterEqual(width, self.target_length) + + value = float( + constraint_void( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertLessEqual(value, 0.0) + + def test_uniform_void_design_is_feasible(self): + rho_filtered = jnp.zeros((self.nx, self.ny)) + rho_projected = jnp.zeros((self.nx, self.ny)) + + value = float( + constraint_void( + rho_filtered, rho_projected, self.resolution, self.target_length + ) + ) + self.assertAlmostEqual(value, -1.0, places=6) + + +class TestAnyProjectionOrder(LengthscaleFixture): + """The constraints only see `rho_projected`, so any SSP order works.""" + + def _design(self): + rng = np.random.default_rng(42) + return rng.random((self.nx, self.ny)) + + def test_gradient_flows_through_each_projection(self): + projections = { + "tanh_projection": tanh_projection, + "ssp1_bilinear": ssp1_bilinear, + "ssp2": ssp2, + } + rho = self._design() + + for name, projection in projections.items(): + for constraint in (constraint_solid, constraint_void): + with self.subTest(projection=name, constraint=constraint.__name__): + + def objective(rho_flat, projection=projection, constraint=constraint): + rho_design = rho_flat.reshape((self.nx, self.ny)) + rho_filtered = conic_filter( + rho_design, + self.filter_radius, + self.lx, + self.ly, + self.resolution, + ) + if projection is tanh_projection: + # A finite beta keeps the plain tanh projection smooth. + rho_projected = projection(rho_filtered, 8.0, self.eta_i) + else: + rho_projected = projection( + rho_filtered, self.beta, self.eta_i, self.resolution + ) + return constraint( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + + value, gradient = value_and_grad(objective)( + jnp.asarray(rho.ravel()) + ) + + self.assertTrue(np.isfinite(float(value))) + gradient = np.asarray(gradient) + self.assertTrue(np.isfinite(gradient).all()) + self.assertGreater(np.linalg.norm(gradient), 0.0) + + def test_jit_matches_eager(self): + rho = self._design() + rho_filtered, rho_projected = self._filter_and_project(rho) + + def constraints(rho_filtered, rho_projected): + return ( + constraint_solid( + rho_filtered, rho_projected, self.resolution, self.target_length + ), + constraint_void( + rho_filtered, rho_projected, self.resolution, self.target_length + ), + ) + + eager = [float(v) for v in constraints(rho_filtered, rho_projected)] + compiled = [ + float(v) for v in jit(constraints)(rho_filtered, rho_projected) + ] + + np.testing.assert_allclose(compiled, eager, rtol=1e-5) + + +class TestGradientsAgainstFiniteDifferences(unittest.TestCase): + """Compare reverse-mode gradients with central finite differences. + + Double precision is required here: the normalized constraint is O(1/epsilon), + so differencing it in single precision is dominated by roundoff. + """ + + def test_adjoints_match_finite_differences(self): + with enable_x64(): + lx = ly = 0.4 + resolution = 50 + target_length = 0.15 + nx = int(np.round(lx * resolution)) + 1 + ny = int(np.round(ly * resolution)) + 1 + + rng = np.random.default_rng(0) + rho = rng.random((nx, ny)) + rho_filtered = conic_filter( + jnp.asarray(rho), target_length, lx, ly, resolution + ) + rho_projected = ssp2(rho_filtered, np.inf, 0.5, resolution) + perturbation = rng.standard_normal((nx, ny)) + step = 1e-6 + + for constraint in (constraint_solid, constraint_void): + for argument in ("rho_filtered", "rho_projected"): + with self.subTest( + constraint=constraint.__name__, argument=argument + ): + + def scalar(value, argument=argument, constraint=constraint): + if argument == "rho_filtered": + return constraint( + value, rho_projected, resolution, target_length + ) + return constraint( + rho_filtered, value, resolution, target_length + ) + + base = ( + rho_filtered + if argument == "rho_filtered" + else rho_projected + ) + adjoint = float( + np.sum(np.asarray(grad(scalar)(base)) * perturbation) + ) + finite_difference = ( + float(scalar(base + step * perturbation)) + - float(scalar(base - step * perturbation)) + ) / (2 * step) + + self.assertAlmostEqual( + adjoint / finite_difference, 1.0, places=5 + ) + + +class TestArgumentValidation(unittest.TestCase): + def setUp(self): + self.rho_filtered = jnp.full((8, 8), 0.5) + self.rho_projected = jnp.zeros((8, 8)) + + def test_negative_target_length(self): + with self.assertRaises(ValueError): + constraint_solid(self.rho_filtered, self.rho_projected, 10, -0.1) + + def test_nonpositive_conic_radius(self): + with self.assertRaises(ValueError): + constraint_void( + self.rho_filtered, self.rho_projected, 10, 0.1, conic_radius=0.0 + ) + + def test_mismatched_shapes(self): + with self.assertRaises(ValueError): + constraint_solid(self.rho_filtered, jnp.zeros((8, 4)), 10, 0.1) + + def test_non_2d_input(self): + with self.assertRaises(ValueError): + constraint_solid(jnp.full((4, 4, 4), 0.5), jnp.zeros((4, 4, 4)), 10, 0.1) + + +if __name__ == "__main__": + unittest.main()