Port MATLAB geometry+process optimization to Python#8
Merged
Conversation
Port the MATLAB optimization framework (cantilever.m optimize_performance and the per-subclass optimize_performance overrides on cantileverPoly / cantileverPiezoelectric) to a single declarative Python module. The MATLAB design used a scaling vector plus five parallel hook methods per subclass (doping_optimization_scaling, doping_cantilever_from_state, doping_current_state, doping_optimization_bounds, doping_initial_conditions_random). The Python port replaces them with a single per-subclass declarative method, optimization_state_vars(), that returns a tuple of StateVar(name, scale, default_min, default_max, log_scale). The optimizer then assembles bounds, scaling, log-transform, and the apply/read functions automatically. This handles the standard geometry+doping pattern (Epitaxy, Diffusion, Implantation), the multi-layer Poly state vector, and the fully custom Piezoelectric state vector through one mechanism. Public API in piezod.optimization (re-exported from piezod): optimize_performance, optimize_performance_from_current, OptimizationResult, StateVar, CantileverMetric, CantileverMetricConstraint, force_resolution_goal, displacement_resolution_goal, surface_stress_resolution_goal, voltage_/charge_force_resolution_goal, voltage_/charge_displacement_resolution_goal. Implementation notes: - SciPy SLSQP is used when nonlinear constraints are present, L-BFGS-B when not. SLSQP often reports success=False on constrained problems even when it found a useful local minimum (the well-known "inequality constraints incompatible" behaviour); the optimizer filters on whether the objective is finite rather than result.success so good solutions are not discarded. - MATLAB's always-on aspect-ratio constraints (l/w >= 2, w/t >= 2, l_pr/w_pr >= 2, l_pr >= 2 um) are added automatically when the cantilever exposes the relevant attributes. Disable with default_aspect_constraints=False. - The shallow-copy + setattr apply path needs a hook for derived attributes. CantileverPiezoelectric mirrors l_si/w_si onto l_pe/w_pe; CantileverPoly enforces t_bot == t_top for symmetric two-piezoresistor layouts. Both are handled by overriding optimization_post_apply. - CantileverImplantation reuses the existing _DEFAULT_BOUNDS table-aware bounds so the new optimizer respects whichever lookup table the cantilever was constructed with. - The existing CantileverImplantation.optimize_doping_for_hooge_noise (implant-process-only optimizer with Hooge-noise default objective) is left untouched; it now coexists with the new generic optimizer rather than being refactored on top of it. Out of scope and deliberately not included: - force_noise_density_goal: the underlying force_noise_density / voltage_noise methods have a pre-existing bug (math.sqrt of an array) that affects callers from this PR. Skipping the goal until that is fixed in a separate PR. Tests: 43 new tests in tests/test_optimization.py covering state machinery, per-subclass state-var declarations, every named goal, constraint validation, end-to-end optimization on all five subclasses, multi-start reproducibility with seeded RNG, parameter and metric constraint enforcement, default aspect-ratio constraints, custom objective callables, and immutability of the input cantilever. Example: python/examples/optimization.py mirrors the MATLAB tutorial's Harley-1999-style optimization flow and improves force resolution from ~115 pN to ~0.24 pN under power, frequency, and stiffness constraints. Tutorial: the optimization section of docs/tutorial.md is rewritten with the actual API.
…_goal
Three layered defects made Cantilever.force_noise_density raise on every
input shape, blocking the MATLAB-equivalent goalForceNoiseDensity from
being ported in this PR:
1. johnson_PSD and thermo_PSD ended with `... * np.ones((1, freq.size))`,
broadcasting a constant across every frequency bin. A Python float has
no `.size`, so passing a scalar raised AttributeError before the
broadcast even ran.
2. voltage_noise summed the four PSD components and applied math.sqrt.
math.sqrt only accepts a scalar (or 0-d array), so any 1-D input
raised TypeError("only 0-dimensional arrays can be converted to Python
scalars").
3. resonant_force_noise_density passed a Python scalar straight into
force_noise_density and returned a 1-element array instead of a
scalar.
Together: force_noise_density(scalar) failed at #1, force_noise_density
(array) failed at #2, and resonant_force_noise_density failed at both #1
and #3. f_min_cumulative (which calls voltage_noise on a logspace array)
and the "Force noise at 1 kHz" line in print_performance were also
silently broken; neither has test coverage.
Fix:
- johnson_PSD and thermo_PSD now promote freq via
`np.atleast_1d(np.asarray(freq, dtype=float))` so scalar input is
accepted. The (1, N) output shape is unchanged for any caller that
was relying on it.
- voltage_noise switches math.sqrt -> np.sqrt, which handles scalar,
0-d, 1-d, and N-d inputs uniformly.
- resonant_force_noise_density wraps the damped frequency in
np.atleast_1d, then squeezes the (1,1) result back to a Python float.
The fix also unblocks the existing f_min_cumulative and
print_performance code paths, which were already broken; these are
strictly additive, not regressions.
With the noise chain fixed, the previously-omitted
force_noise_density_goal is restored. It evaluates
force_noise_density at the damped resonance and applies the MATLAB
1e12 (pN/sqrt(Hz)) scaling. Two new tests cover the goal: one asserts
finite scalar output with the right unit relationship to
resonant_force_noise_density (factor of 1e-3, fN -> pN), and one drives
end-to-end optimization with this goal under a power-dissipation
constraint to confirm it improves.
… temp metrics While auditing the new optimization paths I confirmed three more pre-existing defects that affected this PR's surface: 1. k_base() did `return k_x[0]` on a Python float. The "approx" branch of k_x() returns a scalar; the "exact" branch returns a numpy scalar. Indexing either with [0] raises. The legacy bracket was a MATLAB-port leftover. Fix: return self.k_x() directly. Unblocks approxTempRise, approxTempRiseAnalytical, and through them the TEMP_TIP_APPROX / TEMP_MAX_APPROX optimizer constraints. 2. plot_noise_spectrum used math.sqrt on the four PSD arrays (same pattern as the voltage_noise fix). Switched to np.sqrt and squeezed the (1, N) outputs of johnson_PSD/thermo_PSD to (N,) so matplotlib gets matched 1-D arrays. 3. print_performance() raised in two places: the "Force noise at 1 kHz" line dropped a (1,1) array into a `:g` format spec, and the "Alpha" line formatted a bound method instead of calling it. Both now squeeze/call to a scalar before formatting. The FD-based exact temperature solver (calculateMaxAndTipTemp / calculateTempProfile) has additional unported MATLAB-isms (1-indexed loops, sequence-into-scalar slot assignments) and would need a substantial rewrite. Out of scope here. Consequently this commit also drops TEMP_TIP_EXACT and TEMP_MAX_EXACT from the public CantileverMetric enum -- they would dispatch to calculateMaxAndTipTemp and never produce a finite value. The TEMP_TIP_APPROX and TEMP_MAX_APPROX variants remain and now actually work; the docstring on CantileverMetric notes this. The "F-D Temp Rises" line in print_performance is also dropped pending an FD-solver port; the lumped-circuit "Approx. Temp Rises" line still prints. Six new tests in TestNoiseChainAndPrintPerformance cover voltage_noise on scalar+array input, resonant_force_noise_density returning a Python float, k_base returning a finite scalar, approxTempRise returning finite numbers, print_performance running end-to-end with output spot-checks, and optimize_performance_from_current accepting a TEMP_MAX_APPROX constraint and respecting it on the optimized cantilever.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Ports the MATLAB optimization framework (
optimize_performance/optimize_performance_from_currentplus per-subclass overrides) to asingle declarative Python module. Brings Python optimization to feature
parity with MATLAB across all five cantilever subclasses (Epitaxy,
Diffusion, Implantation, Poly, Piezoelectric).
The MATLAB design used a scaling vector plus five parallel hook methods
per subclass. The Python port replaces them with a single per-subclass
`optimization_state_vars()` returning a tuple of
`StateVar(name, scale, default_min, default_max, log_scale)`. The
optimizer assembles bounds, scaling, log-transform, and the apply/read
functions automatically. Handles the standard geometry+doping pattern,
the multi-layer Poly state vector, and the fully custom Piezoelectric
state vector through one mechanism.
What's new
and `max_iterations` cap, mirroring MATLAB's
`numOptimizationIterations`.
state.
`CantileverMetricConstraint` dataclasses/enums.
`displacement_resolution_goal`, `force_noise_density_goal`,
`surface_stress_resolution_goal`, plus `voltage_` / `charge_`
variants for piezoelectric.
needed) on all five cantilever subclasses.
`python/examples/optimization.py`, rewritten optimization section in
`docs/tutorial.md`.
Implementation notes
`L-BFGS-B` when not. SLSQP often reports `success=False` even on
usable solutions; the optimizer filters on objective finiteness
instead.
`l_pr/w_pr >= 2`, `l_pr >= 2 um`) are added automatically when the
cantilever exposes the relevant attributes. Disable via
`default_aspect_constraints=False`.
and `CantileverPoly` enforces `t_bot == t_top` for symmetric
two-piezoresistor layouts, both via `optimization_post_apply`.
table-aware bounds so the new optimizer respects whichever lookup
table the cantilever was constructed with.
`CantileverImplantation.optimize_doping_for_hooge_noise`
(implant-process-only Hooge-noise optimizer) is left untouched; it now
co-exists with the new generic optimizer.
Side fixes to pre-existing defects
Auditing while wiring up the optimizer surfaced several defects in the
cantilever module that the optimizer paths run into. Fixing them was
in-scope here:
`johnson_PSD` / `thermo_PSD` requiring `freq.size`,
`voltage_noise` calling `math.sqrt` on arrays, and
`resonant_force_noise_density` returning a 1-element array. Fixed by
promoting `freq` to `np.atleast_1d`, swapping to `np.sqrt`, and
squeezing to a scalar.
returning `self.k_x()` directly, which unblocks `approxTempRise`,
`approxTempRiseAnalytical`, and the `TEMP_TIP_APPROX` /
`TEMP_MAX_APPROX` optimizer constraints.
`np.sqrt` plus a squeeze for matplotlib.
((1,1) array vs `:g` format) and the "Alpha" line
(`self.alpha:g` formatted a bound method). Both fixed.
Known follow-up (not addressed here)
`calculateMaxAndTipTemp`) has unported MATLAB-style 1-indexed loops
and sequence-into-scalar slot assignments. Restoring the FD-based
exact temperature solver requires a careful rewrite that would
significantly expand this PR. `TEMP_TIP_EXACT` /
`TEMP_MAX_EXACT` are removed from the `CantileverMetric` enum
pending that work; the `APPROX` lumped-circuit variants remain and
now work. The "F-D Temp Rises" line in `print_performance` is
similarly dropped for now.
Test plan
Harley-1999-style force resolution by several hundred times under
power/frequency/stiffness constraints
epitaxial cantilever (regression-tested)