diff --git a/docs/tutorial.md b/docs/tutorial.md index 54c2a1d..e637e6e 100644 --- a/docs/tutorial.md +++ b/docs/tutorial.md @@ -131,7 +131,59 @@ See [`python/examples/`](../python/examples/) for a runnable script for each. ## Optimization -Design optimization (resolution / noise PSD goal functions with parameter and nonlinear constraints) is implemented in the MATLAB version but has not yet been ported to Python. Track [GitHub Issues](https://github.com/MicrosystemsLab/PiezoD/issues) for progress, or use the MATLAB implementation when optimization is required. +Design optimization is built around three pieces: + +- A **goal callable** that maps a cantilever to the scalar to minimize. PiezoD ships factories that match the MATLAB tutorial's units (pN, nm, pN/sqrt(Hz), etc.): + +```python +from piezod import ( + force_resolution_goal, + displacement_resolution_goal, + force_noise_density_goal, + surface_stress_resolution_goal, +) +``` + +- Optional **parameter_constraints** that override the default state-variable bounds (e.g. clamp thickness, cap bias voltage). Keys are `min_` / `max_` for any state variable returned by `c.optimization_state_vars()`. + +- Optional **metric_constraints** -- inequality constraints on derived quantities like power dissipation, resonant frequency, or stiffness. Each constraint is a `CantileverMetricConstraint` referencing a `CantileverMetric` enum value. + +```python +from piezod import ( + CantileverMetric, + CantileverMetricConstraint, + optimize_performance, +) + +constraints = [ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + CantileverMetricConstraint(CantileverMetric.OMEGA_VACUUM_HZ, minimum=5 * c.freq_max), + CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=1e-3, maximum=1e1), +] + +result = optimize_performance( + c, + force_resolution_goal(), + parameter_constraints={"max_v_bridge": 10.0}, + metric_constraints=constraints, + n_starts=5, + max_iterations=10, + random_seed=0, +) +``` + +`optimize_performance` runs SciPy's SLSQP (or L-BFGS-B when there are no nonlinear constraints) from random initial conditions and returns the best result once two starts agree within `convergence_tolerance` (1% by default), capped at `max_iterations` total runs. Use `optimize_performance_from_current` for a single-shot refinement of the existing design. + +```python +print(f"force_resolution: {c.force_resolution() * 1e12:.1f} pN -> " + f"{result.optimized.force_resolution() * 1e12:.2f} pN") +``` + +Default geometric sanity constraints (`l/w >= 2`, `w/t >= 2`, `l_pr/w_pr >= 2`, `l_pr >= 2 um`) are added automatically; pass `default_aspect_constraints=False` to opt out. + +For a runnable end-to-end example with full output, see [`python/examples/optimization.py`](../python/examples/optimization.py). The example mirrors the MATLAB tutorial flow and typically improves Harley-1999-style force resolution by several hundred times. + +For implant-process-only optimization (geometry fixed, only `annealing_time/temp` and `implantation_energy/dose` varied), `CantileverImplantation.optimize_doping_for_hooge_noise` is also available with a Hooge-noise-limited default objective. ## Next Steps diff --git a/python/examples/optimization.py b/python/examples/optimization.py new file mode 100644 index 0000000..f539ffe --- /dev/null +++ b/python/examples/optimization.py @@ -0,0 +1,93 @@ +"""Cantilever optimization example. + +Mirrors the optimization section of the MATLAB tutorial: an epitaxial +piezoresistive cantilever is optimized for force resolution subject to +power-dissipation and resonant-frequency constraints, with the bridge +voltage capped at 10 V and the thickness floor at 1 um. + +Run: + uv run python examples/optimization.py +""" + +from __future__ import annotations + +import warnings + +from piezod import ( + CantileverEpitaxy, + CantileverMetric, + CantileverMetricConstraint, + force_resolution_goal, + optimize_performance, +) + + +def main() -> None: + # Starting design: 200 um x 20 um x 2 um epitaxial boron piezoresistor. + c = CantileverEpitaxy( + freq_min=1, + freq_max=1000, + l=200e-6, + w=20e-6, + t=2e-6, + l_pr_ratio=0.3, + v_bridge=2.0, + doping_type="boron", + dopant_concentration=1e19, + t_pr_ratio=0.3, + ) + c.fluid = "vacuum" + c.number_of_piezoresistors = 4 + + print("=== Starting design ===") + print(f" geometry: l={c.l * 1e6:.1f} um, w={c.w * 1e6:.1f} um, t={c.t * 1e6:.2f} um") + print(f" v_bridge: {c.v_bridge:.2f} V") + print(f" N_dopant: {c.dopant_concentration:.2e} cm^-3") + print(f" t_pr_ratio: {c.t_pr_ratio:.3f}") + print(f" force_resolution: {c.force_resolution() * 1e12:.2f} pN") + print(f" power: {c.power_dissipation() * 1e3:.3f} mW") + print(f" vacuum freq: {c.omega_vacuum_hz() / 1e3:.1f} kHz") + print(f" stiffness: {c.stiffness():.3g} N/m") + print() + + constraints = [ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + CantileverMetricConstraint(CantileverMetric.OMEGA_VACUUM_HZ, minimum=5 * 1000), + CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=1e-3, maximum=1e1), + ] + + print("Optimizing for force resolution with multi-start (this can take ~10 s)...") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + result = optimize_performance( + c, + force_resolution_goal(), + parameter_constraints={"max_v_bridge": 10.0}, + metric_constraints=constraints, + n_starts=5, + max_iterations=10, + random_seed=0, + ) + + opt = result.optimized + print(f" ran {len(result.all_results)} starts; best objective {result.objective_value:.3f} pN") + print() + + print("=== Optimized design ===") + print(f" geometry: l={opt.l * 1e6:.1f} um, w={opt.w * 1e6:.1f} um, t={opt.t * 1e6:.3f} um") + print(f" v_bridge: {opt.v_bridge:.2f} V") + print(f" N_dopant: {opt.dopant_concentration:.2e} cm^-3") + print(f" t_pr_ratio: {opt.t_pr_ratio:.3f}") + print(f" force_resolution: {opt.force_resolution() * 1e12:.3f} pN") + print(f" power: {opt.power_dissipation() * 1e3:.3f} mW (max {2.0:.1f})") + print(f" vacuum freq: {opt.omega_vacuum_hz() / 1e3:.2f} kHz (min {5.0:.1f})") + print(f" stiffness: {opt.stiffness():.3g} N/m") + print(f" l/w ratio: {opt.l / opt.w:.2f} (default >= 2)") + print(f" w/t ratio: {opt.w / opt.t:.2f} (default >= 2)") + print() + improvement = c.force_resolution() / opt.force_resolution() + print(f"Force resolution improved by {improvement:.1f}x") + + +if __name__ == "__main__": + main() diff --git a/python/src/piezod/__init__.py b/python/src/piezod/__init__.py index 6a0456e..c253fdc 100644 --- a/python/src/piezod/__init__.py +++ b/python/src/piezod/__init__.py @@ -20,6 +20,22 @@ PiezoMaterial, ) from piezod.cantilever_poly import CantileverPoly, Material +from piezod.optimization import ( + CantileverMetric, + CantileverMetricConstraint, + OptimizationResult, + StateVar, + charge_displacement_resolution_goal, + charge_force_resolution_goal, + displacement_resolution_goal, + force_noise_density_goal, + force_resolution_goal, + optimize_performance, + optimize_performance_from_current, + surface_stress_resolution_goal, + voltage_displacement_resolution_goal, + voltage_force_resolution_goal, +) from piezod.piezoresistance import ( CrystalOrientation, default_orientation, @@ -33,6 +49,8 @@ "CantileverDiffusion", "CantileverEpitaxy", "CantileverImplantation", + "CantileverMetric", + "CantileverMetricConstraint", "CantileverPiezoelectric", "CantileverPoly", "CrystalOrientation", @@ -42,13 +60,25 @@ "FluidType", "Material", "MetricConstraint", + "OptimizationResult", "PiezoMaterial", "PiezoresistorFromProfile", + "StateVar", + "charge_displacement_resolution_goal", + "charge_force_resolution_goal", "default_orientation", "diffusion_length_cm", + "displacement_resolution_goal", + "force_noise_density_goal", + "force_resolution_goal", "hooge_alpha_from_anneal", + "optimize_performance", + "optimize_performance_from_current", "pi_low_doping", "rotate_in_plane_stress", + "surface_stress_resolution_goal", + "voltage_displacement_resolution_goal", + "voltage_force_resolution_goal", ] __version__ = version("piezod") diff --git a/python/src/piezod/cantilever.py b/python/src/piezod/cantilever.py index bcbcedb..cae5609 100644 --- a/python/src/piezod/cantilever.py +++ b/python/src/piezod/cantilever.py @@ -758,7 +758,9 @@ def print_performance(self) -> None: omega_damped_hz, Q = self.omega_damped_hz_and_Q() x, active_doping, total_doping = self.doping_profile() TMax_approx, TTip_approx = self.approxTempRise() - TMax, TTip = self.calculateMaxAndTipTemp() + # NOTE: the FD temperature solver (calculateMaxAndTipTemp / + # calculateTempProfile) has unported MATLAB-isms that raise on the + # default cantilever_type. Skip the F-D line until that is ported. thermoLimit = self.thermo_integrated() / self.force_sensitivity() print("=======================") @@ -769,7 +771,9 @@ def print_performance(self) -> None: print(f"PR Length Ratio: {self.l_pr_ratio:g}") print(f"Force resolution (N): {self.force_resolution():g}") - print(f"Force noise at 1 kHz (fN): {self.force_noise_density(1e3):g}") + # force_noise_density returns a (1,1) array even for a scalar + # frequency input; squeeze to a Python float for the format spec. + print(f"Force noise at 1 kHz (fN): {float(np.squeeze(self.force_noise_density(1e3))) * 1e15:g}") print(f"Displacement resolution (m): {self.displacement_resolution():g}") print(f"Force sensitivity (V/N): {self.force_sensitivity():g}") print(f"Displacement sensitivity (V/m): {self.displacement_sensitivity():g}") @@ -786,7 +790,6 @@ def print_performance(self) -> None: print(f"Sheet Resistance: {self.sheet_resistance()}") print(f"Power dissipation (mW): {self.power_dissipation() * 1e3:g}") print(f"Approx. Temp Rises (C) - Tip: {TTip_approx} Max: {TMax_approx}") - print(f"F-D Temp Rises (C) - Tip: {TTip} Max: {TMax}") print(f"Integrated noise (V): {self.integrated_noise():g}") print(f"Integrated johnson noise (V): {self.johnson_integrated():g}") @@ -802,7 +805,7 @@ def print_performance(self) -> None: print(f"Number of silicon resistors: {self.number_of_piezoresistors}") print(f"Si Thermal Conductivity (W/m-K): {self.k_base()}") print(f"E (GPa): {self.modulus() * 1e-9}") - print(f"Alpha: {self.alpha:g}") + print(f"Alpha: {self.alpha():g}") if self.cantilever_type == "step": print("=======================") @@ -970,6 +973,7 @@ def hooge_integrated(self): # Johnson noise PSD from the entire Wheatstone bridge (V^2/Hz) def johnson_PSD(self, freq): + freq = np.atleast_1d(np.asarray(freq, dtype=float)) resistance = self.resistance() # resistance() includes contacts TPR = self.piezoresistor_temp() resistance *= 1 + TPR * self.TCR @@ -1001,6 +1005,7 @@ def piezoresistor_temp(self): # Thermomechanical noise PSD # Units: V^2/Hz def thermo_PSD(self, freq): + freq = np.atleast_1d(np.asarray(freq, dtype=float)) omega_damped_hz, Q_M = self.omega_damped_hz_and_Q() TPR = self.piezoresistor_temp() return ( @@ -1134,11 +1139,11 @@ def integrated_noise(self): + amplifier_integrated**2 ) - # Calculate the noise at a given frequency (V/rtHz) + # Calculate the noise at a given frequency (V/rtHz). Accepts either a + # scalar or an array; the result has the same shape semantics. np.sqrt + # is used (not math.sqrt) so array inputs work. def voltage_noise(self, freq): - return math.sqrt( - self.johnson_PSD(freq) + self.hooge_PSD(freq) + self.thermo_PSD(freq) + self.amplifier_PSD(freq) - ) + return np.sqrt(self.johnson_PSD(freq) + self.hooge_PSD(freq) + self.thermo_PSD(freq) + self.amplifier_PSD(freq)) def f_min_cumulative(self): frequency = np.logspace(math.log10(self.freq_min), math.log10(self.freq_max), Cantilever.numFrequencyPoints) @@ -1747,10 +1752,12 @@ def calculateTempProfileTempDependent(self): T_final = T_current return x, Q, T_final - # Calculate the approximate thermal conductivity of the cantilever (W/m-K) + # Calculate the approximate thermal conductivity of the cantilever (W/m-K). + # k_x() already returns the scalar effective conductivity; the legacy + # ``k_x[0]`` indexing was a MATLAB-port leftover and raised on the float + # returned by the "approx" branch. def k_base(self): - k_x = self.k_x() - return k_x[0] + return self.k_x() # Calculate k(x) as a function of cantilever thickness # None/approx = Asheghi (1997) model @@ -1884,9 +1891,9 @@ def force_noise_density(self, freq): return self.voltage_noise(freq) / self.force_sensitivity() def resonant_force_noise_density(self): - [omega_damped_hz, Q] = self.omega_damped_hz_and_Q() ##ok<*NASGU> - freq = omega_damped_hz - return self.force_noise_density(freq) * 1e15 + omega_damped_hz, _ = self.omega_damped_hz_and_Q() + psd = self.force_noise_density(np.atleast_1d(float(omega_damped_hz))) + return float(np.squeeze(psd)) * 1e15 def surface_stress_resolution(self): return self.integrated_noise() / self.surface_stress_sensitivity() @@ -2292,11 +2299,14 @@ def plot_noise_spectrum(self): freq = np.logspace(math.log10(self.freq_min), math.log10(self.freq_max), Cantilever.numFrequencyPoints) plt.figure() - plt.plot(freq, math.sqrt(self.johnson_PSD(freq))) - plt.plot(freq, math.sqrt(self.hooge_PSD(freq))) - plt.plot(freq, math.sqrt(self.thermo_PSD(freq))) - plt.plot(freq, math.sqrt(self.amplifier_PSD(freq))) - plt.plot(freq, self.voltage_noise(freq)) + # np.sqrt (not math.sqrt) so each PSD array is rooted element-wise. + # Squeeze converts the (1, N) shape returned by johnson_PSD/thermo_PSD + # into (N,) so matplotlib gets matched 1-D x and y. + plt.plot(freq, np.squeeze(np.sqrt(self.johnson_PSD(freq)))) + plt.plot(freq, np.squeeze(np.sqrt(self.hooge_PSD(freq)))) + plt.plot(freq, np.squeeze(np.sqrt(self.thermo_PSD(freq)))) + plt.plot(freq, np.squeeze(np.sqrt(self.amplifier_PSD(freq)))) + plt.plot(freq, np.squeeze(self.voltage_noise(freq))) plt.gca().set_xscale("log") plt.gca().set_yscale("log") plt.ylabel("Noise Voltage Spectral Density (V/rtHz)") diff --git a/python/src/piezod/cantilever_diffusion.py b/python/src/piezod/cantilever_diffusion.py index 6f5f2ed..ff1e0ff 100644 --- a/python/src/piezod/cantilever_diffusion.py +++ b/python/src/piezod/cantilever_diffusion.py @@ -17,6 +17,7 @@ from scipy.special import erfc from piezod.cantilever import Cantilever +from piezod.optimization.state import StateVar class CantileverDiffusion(Cantilever): @@ -340,3 +341,21 @@ def doping_initial_conditions_random(self, parameter_constraints: dict[str, floa diffusion_temp_random = diffusion_temp_min + np.random.rand() * (diffusion_temp_max - diffusion_temp_min) return [diffusion_time_random, diffusion_temp_random] + + def optimization_state_vars(self) -> tuple[StateVar, ...]: + """Declarative state spec for joint geometry + diffusion optimization. + + Seven state variables: cantilever length, width, thickness, + piezoresistor length ratio, bridge bias voltage, diffusion time, + and diffusion temperature. Default bounds match MATLAB + ``cantileverDiffusion``. + """ + return ( + StateVar("l", 1e5, 10e-6, 3e-3), + StateVar("w", 1e7, 2e-6, 100e-6), + StateVar("t", 1e8, 1e-6, 100e-6), + StateVar("l_pr_ratio", 1e2, 0.01, 0.99), + StateVar("v_bridge", 1e1, 0.1, 10.0), + StateVar("diffusion_time", 1e-3, 5 * 60, 90 * 60), + StateVar("diffusion_temp", 1e-3, 273 + 800, 273 + 1000), + ) diff --git a/python/src/piezod/cantilever_epitaxy.py b/python/src/piezod/cantilever_epitaxy.py index 728e3ce..86265a9 100644 --- a/python/src/piezod/cantilever_epitaxy.py +++ b/python/src/piezod/cantilever_epitaxy.py @@ -11,6 +11,16 @@ from numpy.typing import NDArray from piezod.cantilever import Cantilever +from piezod.optimization.state import StateVar + +# Solid-solubility approximations at 1000C used as upper bounds for the +# dopant concentration during optimization. Mirrors MATLAB +# cantileverEpitaxy.doping_optimization_bounds. +_MAX_DOPANT_CONCENTRATION_BY_TYPE: dict[str, float] = { + "boron": 2e20, + "phosphorus": 4e20, + "arsenic": 8e20, +} class CantileverEpitaxy(Cantilever): @@ -245,3 +255,24 @@ def doping_initial_conditions_random(self, parameter_constraints: Optional[dict] x0 = np.array([dopant_concentration_random, t_pr_ratio_random]) return x0 + + def optimization_state_vars(self) -> tuple[StateVar, ...]: + """Declarative state spec for joint geometry + epitaxial optimization. + + Returns the seven optimizable state variables: cantilever length, + width, thickness, piezoresistor length ratio, bridge bias voltage, + epitaxial dopant concentration (log-scale), and piezoresistor + thickness ratio. Default bounds match MATLAB ``cantileverEpitaxy``; + the dopant-concentration upper bound is solid-solubility-based and + depends on ``self.doping_type``. + """ + max_concentration = _MAX_DOPANT_CONCENTRATION_BY_TYPE.get(self.doping_type, 1e20) + return ( + StateVar("l", 1e5, 10e-6, 3e-3), + StateVar("w", 1e7, 2e-6, 100e-6), + StateVar("t", 1e8, 1e-6, 100e-6), + StateVar("l_pr_ratio", 1e2, 0.01, 0.99), + StateVar("v_bridge", 1e1, 0.1, 10.0), + StateVar("dopant_concentration", 1.0, 1e17, max_concentration, log_scale=True), + StateVar("t_pr_ratio", 10.0, 0.01, 0.99), + ) diff --git a/python/src/piezod/cantilever_implantation.py b/python/src/piezod/cantilever_implantation.py index b049776..18f2663 100644 --- a/python/src/piezod/cantilever_implantation.py +++ b/python/src/piezod/cantilever_implantation.py @@ -42,6 +42,7 @@ from scipy.optimize import OptimizeResult, minimize from piezod.cantilever import Cantilever +from piezod.optimization.state import StateVar # Module-level cache for lookup table data, keyed by source name _LOOKUP_TABLE_CACHE: dict[str, dict] = {} @@ -939,6 +940,50 @@ def doping_initial_conditions_random(self, parameter_constraints: dict | None = lb, ub = self.doping_optimization_bounds(parameter_constraints) return lb + np.random.rand(4) * (ub - lb) + def optimization_state_vars(self) -> tuple[StateVar, ...]: + """Declarative state spec for joint geometry + implant optimization. + + Nine state variables: cantilever length, width, thickness, + piezoresistor length ratio, bridge bias voltage, annealing time, + annealing temperature, implant energy, and implant dose + (log-scale). Default doping bounds come from + ``_DEFAULT_BOUNDS[self.lookup_source]`` so the bounds match + whichever lookup table the cantilever was constructed with. + """ + bounds = _DEFAULT_BOUNDS[self.lookup_source] + return ( + StateVar("l", 1e5, 10e-6, 3e-3), + StateVar("w", 1e7, 2e-6, 100e-6), + StateVar("t", 1e8, 1e-6, 100e-6), + StateVar("l_pr_ratio", 1e2, 0.01, 0.99), + StateVar("v_bridge", 1e1, 0.1, 10.0), + StateVar( + "annealing_time", + 1e-2, + bounds["min_annealing_time"], + bounds["max_annealing_time"], + ), + StateVar( + "annealing_temp", + 1e-2, + bounds["min_annealing_temp"], + bounds["max_annealing_temp"], + ), + StateVar( + "implantation_energy", + 1.0, + bounds["min_implantation_energy"], + bounds["max_implantation_energy"], + ), + StateVar( + "implantation_dose", + 1.0, + bounds["min_implantation_dose"], + bounds["max_implantation_dose"], + log_scale=True, + ), + ) + def optimize_doping_for_hooge_noise( self, objective: Callable[[DopingProcessMetrics], float] | None = None, diff --git a/python/src/piezod/cantilever_piezoelectric.py b/python/src/piezod/cantilever_piezoelectric.py index 0444348..6ecb447 100644 --- a/python/src/piezod/cantilever_piezoelectric.py +++ b/python/src/piezod/cantilever_piezoelectric.py @@ -13,6 +13,8 @@ from scipy.interpolate import RegularGridInterpolator from scipy.optimize import minimize_scalar +from piezod.optimization.state import StateVar + class PiezoMaterial(Enum): """Piezoelectric material options.""" @@ -520,3 +522,35 @@ def Qn_total(self) -> float: freq = self.freq_for_plot() Qn = self.Qn(freq) return float(np.sqrt(np.trapezoid(Qn**2, freq))) + + def optimization_post_apply(self) -> None: + """Mirror silicon length/width onto the piezoelectric layer. + + MATLAB ``cantileverPiezoelectric.cantilever_from_state`` performs + the same coupling. Methods like ``stiffness`` and ``tipDeflection`` + read ``l_pe``/``w_pe``, so they must track ``l_si``/``w_si``. + """ + self.l_pe = self.l_si + self.w_pe = self.w_si + + def optimization_state_vars(self) -> tuple[StateVar, ...]: + """Declarative state spec for piezoelectric cantilever optimization. + + Five state variables: silicon length, width, thickness, + piezoelectric layer thickness, and shunt resistance. Default + bounds match MATLAB ``cantileverPiezoelectric``. The piezoelectric + layer geometry tracks the silicon length/width and is not + independently optimized. Material choice and fluid environment are + also not optimized. + + Note: MATLAB's default pins ``r_shunt`` at 1e12 (min == max). To + actually vary it, pass ``parameter_constraints={"min_r_shunt": ..., + "max_r_shunt": ...}``. + """ + return ( + StateVar("l_si", 1e6, 10e-6, 10e-3), + StateVar("w_si", 1e6, 1e-6, 500e-6), + StateVar("t_si", 1e9, 1e-6, 10e-6), + StateVar("t_pe", 1e9, 200e-9, 1e-6), + StateVar("r_shunt", 1e-9, 1e12, 1e12), + ) diff --git a/python/src/piezod/cantilever_poly.py b/python/src/piezod/cantilever_poly.py index 1e2515a..735edf7 100644 --- a/python/src/piezod/cantilever_poly.py +++ b/python/src/piezod/cantilever_poly.py @@ -13,6 +13,8 @@ import numpy as np from numpy.typing import NDArray +from piezod.optimization.state import StateVar + class Material(Enum): """Material types for cantilever layers.""" @@ -598,3 +600,35 @@ def print_performance(self) -> None: print(f"Number of carriers: {self.number_of_carriers():.3g}") print(f"Stiffness: {self.stiffness():.3g} N/m") print(f"Resonant frequency: {self.resonant_frequency():.1f} Hz") + + def optimization_post_apply(self) -> None: + """Enforce top/bottom thickness symmetry for two-piezoresistor layouts. + + ``CantileverPoly.__init__`` forces ``t_bot = t_top`` whenever + ``number_of_piezoresistors_on_cantilever == 2``. The same constraint + must be re-applied during optimization, otherwise the optimizer + could find an asymmetric design that violates the configured + layout. + """ + if self.number_of_piezoresistors_on_cantilever == 2: + self.t_bot = self.t_top + + def optimization_state_vars(self) -> tuple[StateVar, ...]: + """Declarative state spec for joint geometry + multi-layer optimization. + + Eight state variables: cantilever length, width, piezoresistor + length ratio, bridge bias voltage, top/mid/bottom layer + thicknesses, and dopant concentration (log-scale). Default bounds + match MATLAB ``cantileverPoly``. Layer materials and the + ``number_of_piezoresistors`` settings are not optimized. + """ + return ( + StateVar("l", 1e6, 20e-6, 1e-3), + StateVar("w", 1e6, 5e-6, 100e-6), + StateVar("l_pr_ratio", 1.0, 0.01, 1.0), + StateVar("v_bridge", 1.0, 0.1, 10.0), + StateVar("t_top", 1e9, 50e-9, 2e-6), + StateVar("t_mid", 1e9, 10e-9, 2e-6), + StateVar("t_bot", 1e9, 50e-9, 2e-6), + StateVar("dopant_concentration", 1.0, 1e17, 1e20, log_scale=True), + ) diff --git a/python/src/piezod/optimization/__init__.py b/python/src/piezod/optimization/__init__.py new file mode 100644 index 0000000..65ef18d --- /dev/null +++ b/python/src/piezod/optimization/__init__.py @@ -0,0 +1,71 @@ +"""Cantilever optimization framework. + +Public API mirrors the MATLAB optimization toolbox: + +- :func:`optimize_performance` -- multi-start with convergence check +- :func:`optimize_performance_from_current` -- single-shot from current state +- :class:`OptimizationResult` -- structured return value +- :class:`StateVar` -- declarative state-variable spec +- :class:`CantileverMetric` -- enum of constrainable metrics +- :class:`CantileverMetricConstraint` -- inequality constraint dataclass +- Goal factories: :func:`force_resolution_goal`, etc. + +Example: + >>> from piezod import CantileverEpitaxy + >>> from piezod.optimization import ( + ... optimize_performance, + ... force_resolution_goal, + ... CantileverMetric, + ... CantileverMetricConstraint, + ... ) + >>> c = CantileverEpitaxy() # default values + >>> result = optimize_performance( + ... c, + ... force_resolution_goal(), + ... metric_constraints=[ + ... CantileverMetricConstraint( + ... CantileverMetric.POWER_DISSIPATION, maximum=2e-3 + ... ), + ... ], + ... random_seed=0, + ... ) + >>> result.optimized.force_resolution() # doctest: +SKIP +""" + +from piezod.optimization.constraints import ( + CantileverMetric, + CantileverMetricConstraint, +) +from piezod.optimization.goals import ( + charge_displacement_resolution_goal, + charge_force_resolution_goal, + displacement_resolution_goal, + force_noise_density_goal, + force_resolution_goal, + surface_stress_resolution_goal, + voltage_displacement_resolution_goal, + voltage_force_resolution_goal, +) +from piezod.optimization.optimizer import ( + OptimizationResult, + optimize_performance, + optimize_performance_from_current, +) +from piezod.optimization.state import StateVar + +__all__ = [ + "CantileverMetric", + "CantileverMetricConstraint", + "OptimizationResult", + "StateVar", + "charge_displacement_resolution_goal", + "charge_force_resolution_goal", + "displacement_resolution_goal", + "force_noise_density_goal", + "force_resolution_goal", + "optimize_performance", + "optimize_performance_from_current", + "surface_stress_resolution_goal", + "voltage_displacement_resolution_goal", + "voltage_force_resolution_goal", +] diff --git a/python/src/piezod/optimization/constraints.py b/python/src/piezod/optimization/constraints.py new file mode 100644 index 0000000..e7932fc --- /dev/null +++ b/python/src/piezod/optimization/constraints.py @@ -0,0 +1,284 @@ +"""Nonlinear constraint scaffolding for cantilever optimization. + +User-supplied :class:`CantileverMetricConstraint` entries plus a small set of +always-on geometric sanity constraints are translated into the +``scipy.optimize.minimize`` SLSQP constraint format. The MATLAB +``optimization_constraints`` method dispatched on string keys +(``'omega_min_hz'``, ``'min_k'``, etc.); the Python equivalent uses an enum +so typos fail at construction time. +""" + +from __future__ import annotations + +import copy +from collections.abc import Callable, Sequence +from dataclasses import dataclass +from enum import Enum +from typing import Any + +import numpy as np + +from piezod.optimization.state import ( + StateVar, + apply_physical_state, + to_physical_state, +) + +_RESIDUAL_SCALE_FLOOR = 1.0 +_DEFAULT_MIN_PR_LENGTH_M = 2e-6 +_DEFAULT_MIN_L_W_RATIO = 2.0 +_DEFAULT_MIN_W_T_RATIO = 2.0 +_DEFAULT_MIN_PR_L_W_RATIO = 2.0 + + +class CantileverMetric(str, Enum): + """Derived metrics that can be constrained during optimization. + + Each value resolves to a method (or a method plus tuple index) on a + :class:`Cantilever`. Constraints reference metrics by enum so + misspellings fail at construction time rather than silently being + skipped, the way MATLAB's ``exist('omega_min_hz', 'var')`` pattern fails. + """ + + FORCE_RESOLUTION = "force_resolution" + DISPLACEMENT_RESOLUTION = "displacement_resolution" + POWER_DISSIPATION = "power_dissipation" + OMEGA_VACUUM_HZ = "omega_vacuum_hz" + OMEGA_DAMPED_HZ = "omega_damped_hz" + STIFFNESS = "stiffness" + # Approximate (lumped circuit) temperature rises. The FD-based exact + # variants (MATLAB calculateMaxAndTipTemp) are not yet ported -- their + # underlying calculateTempProfile has MATLAB-style 1-indexed/sequence + # bugs. Use the APPROX metrics below as a temperature constraint until + # the FD solver is ported. + TEMP_TIP_APPROX = "temp_tip_approx" + TEMP_MAX_APPROX = "temp_max_approx" + TIP_DEFLECTION = "tip_deflection" + SURFACE_STRESS_RESOLUTION = "surface_stress_resolution" + + +@dataclass(frozen=True) +class CantileverMetricConstraint: + """Inequality constraint on a derived cantilever metric. + + Either ``minimum`` or ``maximum`` (or both) must be set. With both, the + constraint becomes a closed interval. To pin a metric to a specific + value, set ``minimum == maximum``. + """ + + metric: CantileverMetric + minimum: float | None = None + maximum: float | None = None + + def __post_init__(self) -> None: + if self.minimum is None and self.maximum is None: + raise ValueError( + f"CantileverMetricConstraint for {self.metric.value!r} must specify at least one of minimum or maximum." + ) + if self.minimum is not None and not np.isfinite(self.minimum): + raise ValueError(f"minimum must be finite, got {self.minimum!r}.") + if self.maximum is not None and not np.isfinite(self.maximum): + raise ValueError(f"maximum must be finite, got {self.maximum!r}.") + if self.minimum is not None and self.maximum is not None and self.minimum > self.maximum: + raise ValueError(f"minimum ({self.minimum}) must be <= maximum ({self.maximum}).") + + +def _call_or_attr(obj: Any, name: str) -> float: + """Read a quantity that might be a method, property, or plain attribute.""" + value = getattr(obj, name) + if callable(value): + value = value() + return float(value) + + +def evaluate_metric(cantilever: Any, metric: CantileverMetric) -> float: + """Evaluate a metric on a cantilever instance, handling tuple-returning methods.""" + if metric == CantileverMetric.OMEGA_DAMPED_HZ: + omega_damped_hz, _ = cantilever.omega_damped_hz_and_Q() + return float(omega_damped_hz) + if metric == CantileverMetric.TEMP_MAX_APPROX: + t_max, _ = cantilever.approxTempRise() + return float(t_max) + if metric == CantileverMetric.TEMP_TIP_APPROX: + _, t_tip = cantilever.approxTempRise() + return float(t_tip) + if metric == CantileverMetric.TIP_DEFLECTION: + return float(cantilever.tipDeflection()) + return _call_or_attr(cantilever, metric.value) + + +def _residual_scale(limit: float) -> float: + return max(abs(limit), _RESIDUAL_SCALE_FLOOR) + + +def _make_state_applier( + cantilever_template: Any, + state_vars: Sequence[StateVar], +) -> Callable[[np.ndarray], Any]: + """Build a closure that turns an optimizer-space state into a fresh cantilever copy. + + Uses ``copy.copy`` (shallow) so attached lookup tables / arrays are + shared across iterations rather than re-allocated. Mirrors the pattern + used by ``CantileverImplantation._copy_with_doping_state``. + """ + + def apply(optimizer_state: np.ndarray) -> Any: + physical = to_physical_state(optimizer_state, state_vars) + candidate = copy.copy(cantilever_template) + apply_physical_state(candidate, state_vars, physical) + return candidate + + return apply + + +def build_metric_constraint_funcs( + cantilever_template: Any, + state_vars: Sequence[StateVar], + metric_constraints: Sequence[CantileverMetricConstraint], +) -> list[dict]: + """Translate :class:`CantileverMetricConstraint` entries into SLSQP ``ineq`` dicts. + + Each constraint becomes one or two SLSQP entries (one for min, one for max). + SLSQP treats inequality constraints as ``g(x) >= 0`` for feasible, so the + residual returned is positive when the bound is satisfied. + """ + apply_state = _make_state_applier(cantilever_template, state_vars) + constraints: list[dict] = [] + + for constraint in metric_constraints: + metric = constraint.metric + + if constraint.minimum is not None: + lo = constraint.minimum + scale = _residual_scale(lo) + + def make_min_fun(metric=metric, lo=lo, scale=scale): + def fun(x: np.ndarray) -> float: + candidate = apply_state(x) + try: + value = evaluate_metric(candidate, metric) + except (ValueError, FloatingPointError): + return -1.0 + if not np.isfinite(value): + return -1.0 + return (value - lo) / scale + + return fun + + constraints.append({"type": "ineq", "fun": make_min_fun()}) + + if constraint.maximum is not None: + hi = constraint.maximum + scale = _residual_scale(hi) + + def make_max_fun(metric=metric, hi=hi, scale=scale): + def fun(x: np.ndarray) -> float: + candidate = apply_state(x) + try: + value = evaluate_metric(candidate, metric) + except (ValueError, FloatingPointError): + return -1.0 + if not np.isfinite(value): + return -1.0 + return (hi - value) / scale + + return fun + + constraints.append({"type": "ineq", "fun": make_max_fun()}) + + return constraints + + +def _has_attr(obj: Any, name: str) -> bool: + try: + getattr(obj, name) + except (AttributeError, NotImplementedError): + return False + return True + + +def build_default_aspect_ratio_constraints( + cantilever_template: Any, + state_vars: Sequence[StateVar], +) -> list[dict]: + """Add the always-on geometric sanity constraints from MATLAB. + + MATLAB enforces ``l/w >= 2``, ``w/t >= 2``, ``l_pr/w_pr >= 2``, and + ``l_pr >= 2 um`` on every standard optimization. These prevent the + optimizer from picking geometrically silly cantilevers. We add them + only when the cantilever exposes the relevant attributes -- e.g. + Piezoelectric uses ``l_si``/``w_si``/``t_si`` and has no piezoresistor, + so its applicable defaults are different. + """ + apply_state = _make_state_applier(cantilever_template, state_vars) + constraints: list[dict] = [] + + has_l = _has_attr(cantilever_template, "l") + has_w = _has_attr(cantilever_template, "w") + has_t = _has_attr(cantilever_template, "t") + has_l_si = _has_attr(cantilever_template, "l_si") + has_w_si = _has_attr(cantilever_template, "w_si") + has_t_si = _has_attr(cantilever_template, "t_si") + has_l_pr = _has_attr(cantilever_template, "l_pr") + has_w_pr = _has_attr(cantilever_template, "w_pr") + + def safe_eval(x: np.ndarray, fn: Callable[[Any], float]) -> float: + try: + value = fn(apply_state(x)) + except (ValueError, FloatingPointError, ZeroDivisionError): + return -1.0 + if not np.isfinite(value): + return -1.0 + return float(value) + + if has_l and has_w: + constraints.append( + { + "type": "ineq", + "fun": lambda x: safe_eval(x, lambda c: c.l / c.w - _DEFAULT_MIN_L_W_RATIO), + } + ) + elif has_l_si and has_w_si: + constraints.append( + { + "type": "ineq", + "fun": lambda x: safe_eval(x, lambda c: c.l_si / c.w_si - _DEFAULT_MIN_L_W_RATIO), + } + ) + + if has_w and has_t: + constraints.append( + { + "type": "ineq", + "fun": lambda x: safe_eval(x, lambda c: c.w / c.t - _DEFAULT_MIN_W_T_RATIO), + } + ) + elif has_w_si and has_t_si: + constraints.append( + { + "type": "ineq", + "fun": lambda x: safe_eval(x, lambda c: c.w_si / c.t_si - _DEFAULT_MIN_W_T_RATIO), + } + ) + + if has_l_pr and has_w_pr: + constraints.append( + { + "type": "ineq", + "fun": lambda x: safe_eval( + x, + lambda c: _call_or_attr(c, "l_pr") / _call_or_attr(c, "w_pr") - _DEFAULT_MIN_PR_L_W_RATIO, + ), + } + ) + constraints.append( + { + "type": "ineq", + "fun": lambda x: safe_eval( + x, + lambda c: (_call_or_attr(c, "l_pr") - _DEFAULT_MIN_PR_LENGTH_M) * 1e6, + ), + } + ) + + return constraints diff --git a/python/src/piezod/optimization/goals.py b/python/src/piezod/optimization/goals.py new file mode 100644 index 0000000..e3cc127 --- /dev/null +++ b/python/src/piezod/optimization/goals.py @@ -0,0 +1,82 @@ +"""Named objective factories matching MATLAB's ``goal*`` constants. + +MATLAB exposed four goals (``goalForceResolution``, +``goalDisplacementResolution``, ``goalForceNoiseDensity``, +``goalSurfaceStress``) and the piezoelectric subclass added voltage-mode +and charge-mode variants. Each MATLAB ``optimize_*`` wrapper applied a +unit-scaling factor purely for fmincon's numerical conditioning: pN, nm, +pN/sqrt(Hz), MPa. The Python factories preserve those scaling choices so +optimizer-space objective values stay in a sensible range and so test +expectations match the MATLAB tutorial's reported numbers. +""" + +from __future__ import annotations + +from collections.abc import Callable +from typing import Any + +import numpy as np + +# All goals share the signature ``Callable[[Cantilever], float]`` and are +# typed as ``Callable[[Any], float]`` here because the cantilever +# subclasses (Epitaxy/Diffusion/Implantation, Poly, Piezoelectric) do not +# share a common base class and this avoids an import cycle with +# cantilever.py. +Goal = Callable[[Any], float] + + +def force_resolution_goal() -> Goal: + """Goal: minimum detectable force in pN (force_resolution * 1e12).""" + return lambda c: float(c.force_resolution()) * 1e12 + + +def displacement_resolution_goal() -> Goal: + """Goal: minimum detectable displacement in nm (displacement_resolution * 1e9).""" + return lambda c: float(c.displacement_resolution()) * 1e9 + + +def force_noise_density_goal() -> Goal: + """Goal: force noise density at the damped resonance, in pN/sqrt(Hz). + + Mirrors MATLAB's ``optimize_resonant_force_noise_density``. Evaluates + ``Cantilever.force_noise_density`` at the damped resonance and applies + the MATLAB 1e12 (pN) scaling. ``force_noise_density`` returns a 2-D + ``(1, 1)`` array for a single frequency, so the result is squeezed to + a scalar. + """ + + def goal(c: Any) -> float: + omega_damped_hz, _ = c.omega_damped_hz_and_Q() + psd = c.force_noise_density(np.atleast_1d(float(omega_damped_hz))) + return float(np.squeeze(psd)) * 1e12 + + return goal + + +def surface_stress_resolution_goal() -> Goal: + """Goal: minimum detectable surface stress in N/mm (surface_stress_resolution * 1e6). + + Matches MATLAB's ``optimize_surface_stress_resolution`` scaling. The + underlying SI value is in N/m, so 1e6 converts to N/mm. + """ + return lambda c: float(c.surface_stress_resolution()) * 1e6 + + +def voltage_force_resolution_goal() -> Goal: + """Piezoelectric voltage-mode force resolution in pN (Fminv * 1e12).""" + return lambda c: float(c.Fminv()) * 1e12 + + +def voltage_displacement_resolution_goal() -> Goal: + """Piezoelectric voltage-mode displacement resolution in nm (Xminv * 1e9).""" + return lambda c: float(c.Xminv()) * 1e9 + + +def charge_force_resolution_goal() -> Goal: + """Piezoelectric charge-mode force resolution in pN (Fminq * 1e12).""" + return lambda c: float(c.Fminq()) * 1e12 + + +def charge_displacement_resolution_goal() -> Goal: + """Piezoelectric charge-mode displacement resolution in nm (Xminq * 1e9).""" + return lambda c: float(c.Xminq()) * 1e9 diff --git a/python/src/piezod/optimization/optimizer.py b/python/src/piezod/optimization/optimizer.py new file mode 100644 index 0000000..b4d8720 --- /dev/null +++ b/python/src/piezod/optimization/optimizer.py @@ -0,0 +1,349 @@ +"""Multi-start cantilever optimization built on ``scipy.optimize.minimize``. + +Two entry points mirror MATLAB: + +- :func:`optimize_performance_from_current` -- single shot from the + cantilever's current state, equivalent to MATLAB + ``optimize_performance_from_current``. +- :func:`optimize_performance` -- multi-start with random initial + conditions and a convergence check that stops early when the two best + results agree within ``convergence_tolerance``. Equivalent to MATLAB + ``optimize_performance`` driven by ``cantilever.numOptimizationIterations``. + +The state machinery in :mod:`piezod.optimization.state` and constraint +builders in :mod:`piezod.optimization.constraints` keep this module +focused on the optimizer dispatch and convergence logic. +""" + +from __future__ import annotations + +import copy +import math +from collections.abc import Callable, Sequence +from dataclasses import dataclass +from typing import Any + +import numpy as np +from scipy.optimize import OptimizeResult, minimize + +from piezod.optimization.constraints import ( + CantileverMetricConstraint, + build_default_aspect_ratio_constraints, + build_metric_constraint_funcs, +) +from piezod.optimization.state import ( + StateVar, + apply_physical_state, + optimizer_bounds, + physical_state_from_cantilever, + random_physical_state, + to_optimizer_state, + to_physical_state, +) + +_DEFAULT_TOL = 1e-9 +_DEFAULT_MAXITER = 2000 +_INFEASIBLE_OBJECTIVE = 1e30 + + +@dataclass(frozen=True) +class OptimizationResult: + """Return value of :func:`optimize_performance` and friends. + + Attributes: + optimized: A copy of the input cantilever with the optimal state + applied. Original cantilever is left untouched. + objective_value: Objective function value at ``optimized``, in the + same units as the goal callable returned (e.g. pN for + :func:`force_resolution_goal`). + physical_state: Final state vector in physical units, in the order + given by ``cantilever.optimization_state_vars()``. + scipy_result: The ``OptimizeResult`` from the best run. + all_results: Tuple of ``OptimizeResult`` for every multi-start. + """ + + optimized: object + objective_value: float + physical_state: np.ndarray + scipy_result: OptimizeResult + all_results: tuple[OptimizeResult, ...] + + +def _resolve_method( + method: str | None, + has_constraints: bool, +) -> str: + if method is not None: + if method == "L-BFGS-B" and has_constraints: + raise ValueError("L-BFGS-B cannot enforce nonlinear constraints; use SLSQP.") + return method + return "SLSQP" if has_constraints else "L-BFGS-B" + + +def _run_one( + cantilever_template: Any, + state_vars: Sequence[StateVar], + objective: Callable[[object], float], + optimizer_start: np.ndarray, + bounds: list[tuple[float, float]], + constraints: list[dict], + method: str, + options: dict, +) -> OptimizeResult: + def objective_from_optimizer_state(x: np.ndarray) -> float: + physical = to_physical_state(x, state_vars) + candidate = copy.copy(cantilever_template) + try: + apply_physical_state(candidate, state_vars, physical) + value = float(objective(candidate)) + except (ValueError, FloatingPointError, ZeroDivisionError): + return _INFEASIBLE_OBJECTIVE + if not math.isfinite(value): + return _INFEASIBLE_OBJECTIVE + return value + + minimize_kwargs: dict = { + "method": method, + "bounds": bounds, + "options": options, + } + if constraints: + minimize_kwargs["constraints"] = constraints + + return minimize(objective_from_optimizer_state, optimizer_start, **minimize_kwargs) + + +def _build_optimized_cantilever( + cantilever_template: Any, + state_vars: Sequence[StateVar], + optimizer_state: np.ndarray, +) -> tuple[object, np.ndarray]: + physical = to_physical_state(optimizer_state, state_vars) + optimized = copy.copy(cantilever_template) + apply_physical_state(optimized, state_vars, physical) + return optimized, physical + + +def optimize_performance_from_current( + cantilever: Any, + objective: Callable[[object], float], + *, + parameter_constraints: dict[str, float] | None = None, + metric_constraints: Sequence[CantileverMetricConstraint] | None = None, + default_aspect_constraints: bool = True, + method: str | None = None, + options: dict | None = None, +) -> OptimizationResult: + """Single-shot optimization starting from the cantilever's current state. + + Mirrors MATLAB ``optimize_performance_from_current``. Useful for + refining an existing design rather than searching the full bound box. + + Args: + cantilever: Instance of a :class:`Cantilever` subclass that + implements ``optimization_state_vars()``. The instance is not + mutated; the result is a copy. + objective: Callable mapping a cantilever to a scalar to minimize. + Use one of the named goal factories from + :mod:`piezod.optimization.goals` or supply your own. + parameter_constraints: Optional dict overriding default state + bounds. Keys are ``"min_"`` / ``"max_"`` for any + ``StateVar.name`` returned by ``optimization_state_vars()``. + metric_constraints: Sequence of :class:`CantileverMetricConstraint` + inequality constraints on derived metrics + (``power_dissipation``, ``stiffness``, etc.). + default_aspect_constraints: When True, add the standard + ``l/w >= 2``, ``w/t >= 2``, ``l_pr/w_pr >= 2``, ``l_pr >= 2 um`` + sanity constraints. Mirrors MATLAB's always-on aspect ratio + checks. Disable for advanced cases (sub-2 aspect ratios). + method: SciPy ``minimize`` method. Default chooses SLSQP if any + nonlinear constraints are present, L-BFGS-B otherwise. + options: Override SciPy ``minimize`` options. Defaults: ``ftol``, + ``maxiter``, set to match MATLAB's ``TolFun``/``MaxIter``. + + Returns: + :class:`OptimizationResult` with the optimized cantilever and full + SciPy diagnostic information. + """ + state_vars = cantilever.optimization_state_vars() + bounds = optimizer_bounds(state_vars, parameter_constraints) + physical_start = physical_state_from_cantilever(cantilever, state_vars) + physical_start = _clip_to_bounds(physical_start, state_vars, parameter_constraints) + optimizer_start = to_optimizer_state(physical_start, state_vars) + + constraints = list(metric_constraints) if metric_constraints else [] + slsqp_constraints = build_metric_constraint_funcs(cantilever, state_vars, constraints) + if default_aspect_constraints: + slsqp_constraints.extend(build_default_aspect_ratio_constraints(cantilever, state_vars)) + + resolved_method = _resolve_method(method, has_constraints=bool(slsqp_constraints)) + resolved_options = _build_options(options) + + scipy_result = _run_one( + cantilever_template=cantilever, + state_vars=state_vars, + objective=objective, + optimizer_start=optimizer_start, + bounds=bounds, + constraints=slsqp_constraints, + method=resolved_method, + options=resolved_options, + ) + + optimized, physical = _build_optimized_cantilever(cantilever, state_vars, scipy_result.x) + return OptimizationResult( + optimized=optimized, + objective_value=float(scipy_result.fun), + physical_state=physical, + scipy_result=scipy_result, + all_results=(scipy_result,), + ) + + +def optimize_performance( + cantilever: Any, + objective: Callable[[object], float], + *, + parameter_constraints: dict[str, float] | None = None, + metric_constraints: Sequence[CantileverMetricConstraint] | None = None, + default_aspect_constraints: bool = True, + n_starts: int = 5, + max_iterations: int = 20, + convergence_tolerance: float = 0.01, + include_current_state: bool = False, + random_seed: int | None = None, + method: str | None = None, + options: dict | None = None, +) -> OptimizationResult: + """Multi-start optimization with early-stop convergence check. + + Mirrors MATLAB ``optimize_performance``: run optimization from random + initial conditions until the two best objective values agree within + ``convergence_tolerance``, capped at ``max_iterations`` runs total. + + Args: + cantilever: Instance to optimize (not mutated). + objective: Callable returning the scalar to minimize. + parameter_constraints: Optional state-bound overrides. + metric_constraints: Optional derived-metric inequality constraints. + default_aspect_constraints: Add standard geometric sanity checks. + n_starts: Number of random starts to run before checking convergence. + After that point, additional starts run one at a time until + convergence or ``max_iterations`` is reached. Mirrors MATLAB's + iterative ``ii > 1`` check. + max_iterations: Hard upper bound on total optimization runs. + Equivalent to MATLAB's ``cantilever.numOptimizationIterations``. + convergence_tolerance: Two best objective values must agree within + this relative tolerance for early stop. ``0.01`` = 1%, the + MATLAB default. + include_current_state: When True, the cantilever's current state + is used as the first start before any random ones. + random_seed: Optional seed for reproducibility. + method: SciPy method override. + options: SciPy options override. + + Returns: + :class:`OptimizationResult` whose ``optimized`` cantilever is the + best of all starts. ``all_results`` contains every SciPy run. + """ + if n_starts < 1: + raise ValueError(f"n_starts must be >= 1, got {n_starts}.") + if max_iterations < n_starts: + raise ValueError(f"max_iterations ({max_iterations}) must be >= n_starts ({n_starts}).") + + state_vars = cantilever.optimization_state_vars() + bounds = optimizer_bounds(state_vars, parameter_constraints) + + user_constraints = list(metric_constraints) if metric_constraints else [] + slsqp_constraints = build_metric_constraint_funcs(cantilever, state_vars, user_constraints) + if default_aspect_constraints: + slsqp_constraints.extend(build_default_aspect_ratio_constraints(cantilever, state_vars)) + + resolved_method = _resolve_method(method, has_constraints=bool(slsqp_constraints)) + resolved_options = _build_options(options) + + rng = np.random.default_rng(random_seed) + all_results: list[OptimizeResult] = [] + objective_values: list[float] = [] + + def run_start(physical_start: np.ndarray) -> None: + clipped = _clip_to_bounds(physical_start, state_vars, parameter_constraints) + optimizer_start = to_optimizer_state(clipped, state_vars) + result = _run_one( + cantilever_template=cantilever, + state_vars=state_vars, + objective=objective, + optimizer_start=optimizer_start, + bounds=bounds, + constraints=slsqp_constraints, + method=resolved_method, + options=resolved_options, + ) + all_results.append(result) + # SciPy SLSQP often reports success=False on constrained problems + # ("inequality constraints incompatible") even when it found a + # useful local minimum. Filter only on the objective being finite + # so we don't discard genuinely good solutions. MATLAB filtered on + # exitflag in [1, 2]; SciPy's analog is too coarse. + if math.isfinite(result.fun) and result.fun < _INFEASIBLE_OBJECTIVE / 2: + objective_values.append(float(result.fun)) + else: + objective_values.append(math.inf) + + if include_current_state: + run_start(physical_state_from_cantilever(cantilever, state_vars)) + + while len(all_results) < n_starts: + run_start(random_physical_state(state_vars, parameter_constraints, rng)) + + while len(all_results) < max_iterations: + finite = sorted(v for v in objective_values if math.isfinite(v)) + if len(finite) >= 2: + best, second = finite[0], finite[1] + if best > 0 and abs(1.0 - best / second) < convergence_tolerance: + break + run_start(random_physical_state(state_vars, parameter_constraints, rng)) + + finite_indices = [i for i, v in enumerate(objective_values) if math.isfinite(v)] + if not finite_indices: + raise RuntimeError( + f"All {len(all_results)} optimization starts failed to produce a finite objective. " + "Check parameter_constraints, metric_constraints, and the objective callable." + ) + + best_index = min(finite_indices, key=lambda i: objective_values[i]) + best_result = all_results[best_index] + + optimized, physical = _build_optimized_cantilever(cantilever, state_vars, best_result.x) + return OptimizationResult( + optimized=optimized, + objective_value=float(best_result.fun), + physical_state=physical, + scipy_result=best_result, + all_results=tuple(all_results), + ) + + +def _clip_to_bounds( + physical_state: np.ndarray, + state_vars: Sequence[StateVar], + parameter_constraints: dict[str, float] | None, +) -> np.ndarray: + """Project a physical state vector inside the bound box. + + SciPy's SLSQP requires the initial guess to satisfy bounds (otherwise + it can produce NaNs in the first step). Pure-clipping is sufficient + here -- we don't need a feasible-with-respect-to-nonlinear-constraints + start, just a bound-feasible one. + """ + from piezod.optimization.state import physical_bounds + + lower, upper = physical_bounds(state_vars, parameter_constraints) + return np.clip(physical_state, lower, upper) + + +def _build_options(user_options: dict | None) -> dict: + options = {"ftol": _DEFAULT_TOL, "maxiter": _DEFAULT_MAXITER} + if user_options: + options.update(user_options) + return options diff --git a/python/src/piezod/optimization/state.py b/python/src/piezod/optimization/state.py new file mode 100644 index 0000000..ceea1ca --- /dev/null +++ b/python/src/piezod/optimization/state.py @@ -0,0 +1,198 @@ +"""State-vector machinery for cantilever optimization. + +Each cantilever subclass declares the variables it wants to optimize over by +returning a tuple of :class:`StateVar` from ``optimization_state_vars()``. +This module provides the helpers that convert between physical attribute +values on a cantilever, optimizer-space state vectors (log-transformed and +scaled to O(1)), and bounds derived from defaults plus user overrides. + +The design mirrors the MATLAB implementation's pattern of a scaling vector +plus per-class bound/initial-condition hooks, but expressed as a single +declarative spec instead of five parallel methods. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from dataclasses import dataclass +from typing import Any + +import numpy as np + + +@dataclass(frozen=True) +class StateVar: + """Declarative spec for one optimizable state variable. + + Attributes: + name: Cantilever attribute name (e.g. ``"l"``, ``"dopant_concentration"``). + The optimizer reads/writes this attribute via getattr/setattr. + scale: Multiplicative factor applied in optimizer space so the scaled + state lands near O(1). Mirrors MATLAB's ``optimization_scaling`` + vector. Example: ``"l"`` (meters, ~1e-4) uses ``scale = 1e5``. + default_min: Lower bound in physical units, used unless overridden by + ``parameter_constraints`` keyed ``"min_"``. + default_max: Upper bound in physical units, used unless overridden by + ``parameter_constraints`` keyed ``"max_"``. + log_scale: If True, optimize log10(value) instead of value. Used for + quantities that span many orders of magnitude (e.g. dopant + concentration, implant dose). Bounds are interpreted in physical + space; the log transform is applied internally. + """ + + name: str + scale: float + default_min: float + default_max: float + log_scale: bool = False + + +def physical_state_from_cantilever( + cantilever: Any, + state_vars: Sequence[StateVar], +) -> np.ndarray: + """Read the current physical state vector from a cantilever instance.""" + return np.array([float(getattr(cantilever, var.name)) for var in state_vars], dtype=float) + + +def apply_physical_state( + cantilever: Any, + state_vars: Sequence[StateVar], + physical_state: np.ndarray, +) -> None: + """Apply a physical state vector by setting the corresponding attributes. + + After all state attributes are set, calls + ``cantilever.optimization_post_apply()`` if the method exists. Subclasses + use this to keep derived attributes in sync (e.g. + ``CantileverPiezoelectric`` mirrors ``l_si``/``w_si`` to ``l_pe``/``w_pe``, + ``CantileverPoly`` enforces top/bottom thickness symmetry). + """ + if len(physical_state) != len(state_vars): + raise ValueError(f"physical_state has length {len(physical_state)}, expected {len(state_vars)}.") + for var, value in zip(state_vars, physical_state, strict=True): + setattr(cantilever, var.name, float(value)) + post_apply = getattr(cantilever, "optimization_post_apply", None) + if post_apply is not None: + post_apply() + + +def to_optimizer_state( + physical_state: np.ndarray, + state_vars: Sequence[StateVar], +) -> np.ndarray: + """Convert a physical state vector to optimizer (scaled, log-where-applicable) space.""" + physical_state = np.asarray(physical_state, dtype=float) + if len(physical_state) != len(state_vars): + raise ValueError(f"physical_state has length {len(physical_state)}, expected {len(state_vars)}.") + out = np.empty_like(physical_state) + for i, var in enumerate(state_vars): + value = physical_state[i] + if var.log_scale: + if value <= 0: + raise ValueError(f"State variable {var.name!r} has log_scale=True but value {value} is not positive.") + value = np.log10(value) + out[i] = value * var.scale + return out + + +def to_physical_state( + optimizer_state: np.ndarray, + state_vars: Sequence[StateVar], +) -> np.ndarray: + """Invert the optimizer-space transform back to physical units.""" + optimizer_state = np.asarray(optimizer_state, dtype=float) + if len(optimizer_state) != len(state_vars): + raise ValueError(f"optimizer_state has length {len(optimizer_state)}, expected {len(state_vars)}.") + out = np.empty_like(optimizer_state) + for i, var in enumerate(state_vars): + value = optimizer_state[i] / var.scale + if var.log_scale: + value = 10.0**value + out[i] = value + return out + + +def physical_bounds( + state_vars: Sequence[StateVar], + parameter_constraints: dict[str, float] | None, +) -> tuple[np.ndarray, np.ndarray]: + """Build (lower, upper) bound vectors in physical units. + + ``parameter_constraints`` overrides defaults using keys ``"min_"`` + and ``"max_"``. Unknown keys raise :class:`ValueError` so typos + fail loudly instead of silently being ignored. + """ + lower = np.array([var.default_min for var in state_vars], dtype=float) + upper = np.array([var.default_max for var in state_vars], dtype=float) + + if parameter_constraints is None: + return lower, upper + + name_to_index = {var.name: i for i, var in enumerate(state_vars)} + for key, value in parameter_constraints.items(): + if key.startswith("min_"): + target = key[4:] + if target not in name_to_index: + raise ValueError( + f"parameter_constraints key {key!r} does not match any state variable. " + f"Known names: {sorted(name_to_index)}" + ) + lower[name_to_index[target]] = float(value) + elif key.startswith("max_"): + target = key[4:] + if target not in name_to_index: + raise ValueError( + f"parameter_constraints key {key!r} does not match any state variable. " + f"Known names: {sorted(name_to_index)}" + ) + upper[name_to_index[target]] = float(value) + else: + raise ValueError(f"parameter_constraints key {key!r} must start with 'min_' or 'max_'.") + + if np.any(lower > upper): + bad = [(var.name, float(lower[i]), float(upper[i])) for i, var in enumerate(state_vars) if lower[i] > upper[i]] + raise ValueError(f"parameter_constraints produced lower > upper for: {bad}") + + return lower, upper + + +def optimizer_bounds( + state_vars: Sequence[StateVar], + parameter_constraints: dict[str, float] | None, +) -> list[tuple[float, float]]: + """Build SciPy-style bound list in optimizer (scaled, log) space.""" + lower, upper = physical_bounds(state_vars, parameter_constraints) + lo_opt = to_optimizer_state(lower, state_vars) + hi_opt = to_optimizer_state(upper, state_vars) + return list(zip(lo_opt.tolist(), hi_opt.tolist(), strict=True)) + + +def random_physical_state( + state_vars: Sequence[StateVar], + parameter_constraints: dict[str, float] | None, + rng: np.random.Generator, +) -> np.ndarray: + """Sample a uniformly-random physical state within bounds. + + For ``log_scale`` variables, samples uniformly in log10 space so that + e.g. dopant concentration spans the full bound range rather than + clustering near the upper end. Mirrors MATLAB's + ``doping_initial_conditions_random`` for epitaxy. + """ + lower, upper = physical_bounds(state_vars, parameter_constraints) + out = np.empty(len(state_vars), dtype=float) + for i, var in enumerate(state_vars): + lo = lower[i] + hi = upper[i] + if var.log_scale: + if lo <= 0 or hi <= 0: + raise ValueError( + f"State variable {var.name!r} has log_scale=True but bounds ({lo}, {hi}) are not strictly positive." + ) + log_lo = np.log10(lo) + log_hi = np.log10(hi) + out[i] = 10.0 ** (log_lo + rng.random() * (log_hi - log_lo)) + else: + out[i] = lo + rng.random() * (hi - lo) + return out diff --git a/python/tests/test_optimization.py b/python/tests/test_optimization.py new file mode 100644 index 0000000..9584e53 --- /dev/null +++ b/python/tests/test_optimization.py @@ -0,0 +1,693 @@ +"""Tests for the generic geometry+process cantilever optimizer.""" + +from __future__ import annotations + +import warnings + +import numpy as np +import pytest + +from piezod import ( + CantileverDiffusion, + CantileverEpitaxy, + CantileverImplantation, + CantileverMetric, + CantileverMetricConstraint, + CantileverPiezoelectric, + CantileverPoly, + Material, + OptimizationResult, + PiezoMaterial, + StateVar, + charge_force_resolution_goal, + displacement_resolution_goal, + force_noise_density_goal, + force_resolution_goal, + optimize_performance, + optimize_performance_from_current, + surface_stress_resolution_goal, + voltage_force_resolution_goal, +) +from piezod.optimization.state import ( + physical_bounds, + physical_state_from_cantilever, + random_physical_state, + to_optimizer_state, + to_physical_state, +) + + +@pytest.fixture(autouse=True) +def _silence_warnings(): + """SciPy and the cantilever model emit RuntimeWarnings during normal optimizer + iteration (singular interpolation, divide-by-zero). They are harmless but + noisy in test output.""" + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + yield + + +def _epitaxy_default() -> CantileverEpitaxy: + """Return an epitaxial cantilever inside the default optimization bound box.""" + c = CantileverEpitaxy( + freq_min=10, + freq_max=1000, + l=200e-6, + w=20e-6, + t=2e-6, + l_pr_ratio=0.3, + v_bridge=2.0, + doping_type="boron", + dopant_concentration=1e19, + t_pr_ratio=0.3, + ) + c.fluid = "vacuum" + c.number_of_piezoresistors = 4 + return c + + +def _diffusion_default() -> CantileverDiffusion: + c = CantileverDiffusion( + freq_min=10, + freq_max=1000, + l=200e-6, + w=20e-6, + t=2e-6, + l_pr_ratio=0.3, + v_bridge=2.0, + doping_type="phosphorus", + diffusion_time=30 * 60, + diffusion_temp=273 + 850, + ) + c.fluid = "vacuum" + c.number_of_piezoresistors = 4 + return c + + +def _implantation_default() -> CantileverImplantation: + c = CantileverImplantation( + freq_min=10, + freq_max=1000, + l=200e-6, + w=20e-6, + t=2e-6, + l_pr_ratio=0.3, + v_bridge=2.0, + doping_type="boron", + annealing_time=30 * 60, + annealing_temp=273.15 + 1000, + annealing_type="inert", + implantation_energy=40, + implantation_dose=1e15, + ) + c.fluid = "vacuum" + c.number_of_piezoresistors = 4 + return c + + +def _poly_default() -> CantileverPoly: + return CantileverPoly( + freq_min=10, + freq_max=1000, + l=200e-6, + w=20e-6, + t_top=200e-9, + t_mid=1e-6, + t_bot=200e-9, + matl_top=Material.POLY, + matl_mid=Material.SI, + matl_bot=Material.POLY, + l_pr_ratio=0.3, + v_bridge=2.0, + dopant_concentration=1e19, + number_of_piezoresistors=2, + number_of_piezoresistors_on_cantilever=2, + ) + + +def _piezoelectric_default() -> CantileverPiezoelectric: + return CantileverPiezoelectric( + freq_min=10, + freq_max=1000, + l_si=200e-6, + w_si=20e-6, + t_si=2e-6, + t_pe=500e-9, + material=PiezoMaterial.ALN, + r_shunt=1e12, + ) + + +# ----------------------------------------------------------------------------- +# State machinery +# ----------------------------------------------------------------------------- + + +class TestStateMachinery: + def test_state_var_round_trip_linear(self): + var = StateVar("l", 1e5, 10e-6, 3e-3) + physical = np.array([1e-4]) + opt = to_optimizer_state(physical, [var]) + back = to_physical_state(opt, [var]) + assert np.allclose(back, physical) + + def test_state_var_round_trip_log(self): + var = StateVar("dopant_concentration", 1.0, 1e17, 1e20, log_scale=True) + physical = np.array([1e19]) + opt = to_optimizer_state(physical, [var]) + # log10(1e19) * 1.0 = 19.0 + assert np.isclose(opt[0], 19.0) + back = to_physical_state(opt, [var]) + assert np.isclose(back[0], physical[0]) + + def test_log_scale_rejects_nonpositive(self): + var = StateVar("x", 1.0, 1e-3, 1.0, log_scale=True) + with pytest.raises(ValueError, match="log_scale"): + to_optimizer_state(np.array([0.0]), [var]) + + def test_physical_bounds_default(self): + c = _epitaxy_default() + lower, upper = physical_bounds(c.optimization_state_vars(), None) + assert lower[0] == 10e-6 # min_l + assert upper[0] == 3e-3 # max_l + + def test_physical_bounds_override(self): + c = _epitaxy_default() + lower, upper = physical_bounds( + c.optimization_state_vars(), + {"min_l": 50e-6, "max_v_bridge": 5.0}, + ) + assert lower[0] == 50e-6 + assert upper[4] == 5.0 # v_bridge upper + + def test_physical_bounds_unknown_key(self): + c = _epitaxy_default() + with pytest.raises(ValueError, match="not match any state variable"): + physical_bounds(c.optimization_state_vars(), {"min_nonexistent": 1.0}) + + def test_physical_bounds_bad_prefix(self): + c = _epitaxy_default() + with pytest.raises(ValueError, match="must start with"): + physical_bounds(c.optimization_state_vars(), {"length": 1.0}) + + def test_physical_bounds_lower_above_upper(self): + c = _epitaxy_default() + with pytest.raises(ValueError, match="lower > upper"): + physical_bounds( + c.optimization_state_vars(), + {"min_l": 1e-3, "max_l": 1e-4}, + ) + + def test_random_physical_state_inside_bounds(self): + c = _epitaxy_default() + rng = np.random.default_rng(0) + state_vars = c.optimization_state_vars() + lower, upper = physical_bounds(state_vars, None) + for _ in range(10): + sample = random_physical_state(state_vars, None, rng) + assert np.all(sample >= lower) + assert np.all(sample <= upper) + + def test_random_physical_state_log_scale_distribution(self): + c = _epitaxy_default() + rng = np.random.default_rng(0) + state_vars = c.optimization_state_vars() + # Sample many times; the geometric mean of dopant_concentration + # should be near 10**((17+log10(2e20))/2) ~= 10**18.65 = 4.5e18 + # for boron (max 2e20) + samples = np.array([random_physical_state(state_vars, None, rng) for _ in range(500)]) + dopant = samples[:, 5] + log_mean = np.mean(np.log10(dopant)) + expected_log_mean = (np.log10(1e17) + np.log10(2e20)) / 2 + assert abs(log_mean - expected_log_mean) < 0.2 + + +# ----------------------------------------------------------------------------- +# Per-subclass state-var declarations +# ----------------------------------------------------------------------------- + + +class TestStateVarDeclarations: + @pytest.mark.parametrize( + "factory,expected_count,expected_first", + [ + (_epitaxy_default, 7, "l"), + (_diffusion_default, 7, "l"), + (_implantation_default, 9, "l"), + (_poly_default, 8, "l"), + (_piezoelectric_default, 5, "l_si"), + ], + ) + def test_subclass_state_vars(self, factory, expected_count, expected_first): + c = factory() + state_vars = c.optimization_state_vars() + assert len(state_vars) == expected_count + assert state_vars[0].name == expected_first + # All declared attributes exist on the cantilever + for var in state_vars: + assert hasattr(c, var.name), f"{type(c).__name__} missing attribute {var.name!r}" + # All bounds finite and lower < upper (or equal for pinned) + for var in state_vars: + assert np.isfinite(var.default_min) + assert np.isfinite(var.default_max) + assert var.default_min <= var.default_max + + def test_epitaxy_concentration_bound_depends_on_dopant(self): + c = _epitaxy_default() + c.doping_type = "boron" + boron_max = next(v for v in c.optimization_state_vars() if v.name == "dopant_concentration").default_max + c.doping_type = "arsenic" + arsenic_max = next(v for v in c.optimization_state_vars() if v.name == "dopant_concentration").default_max + assert arsenic_max > boron_max # 8e20 > 2e20 per MATLAB defaults + + def test_implantation_bounds_depend_on_lookup_source(self): + c_tsuprem = CantileverImplantation( + freq_min=10, + freq_max=1000, + l=200e-6, + w=20e-6, + t=2e-6, + l_pr_ratio=0.3, + v_bridge=2.0, + doping_type="boron", + annealing_time=30 * 60, + annealing_temp=273.15 + 1000, + annealing_type="inert", + implantation_energy=40, + implantation_dose=1e15, + lookup_source="tsuprem4", + ) + c_dope = CantileverImplantation( + freq_min=10, + freq_max=1000, + l=200e-6, + w=20e-6, + t=2e-6, + l_pr_ratio=0.3, + v_bridge=2.0, + doping_type="boron", + annealing_time=30 * 60, + annealing_temp=273.15 + 1000, + annealing_type="inert", + implantation_energy=40, + implantation_dose=1e15, + lookup_source="dopedealer", + ) + ts = next(v for v in c_tsuprem.optimization_state_vars() if v.name == "implantation_energy") + dd = next(v for v in c_dope.optimization_state_vars() if v.name == "implantation_energy") + # dopedealer table covers a wider energy range + assert dd.default_min < ts.default_min + assert dd.default_max > ts.default_max + + +# ----------------------------------------------------------------------------- +# Goal factories +# ----------------------------------------------------------------------------- + + +class TestGoalFactories: + def test_force_resolution_goal_units(self): + c = _epitaxy_default() + goal = force_resolution_goal() + # MATLAB-style units: pN + assert np.isclose(goal(c), c.force_resolution() * 1e12) + + def test_displacement_resolution_goal_units(self): + c = _epitaxy_default() + goal = displacement_resolution_goal() + assert np.isclose(goal(c), c.displacement_resolution() * 1e9) + + def test_surface_stress_goal_units(self): + c = _epitaxy_default() + goal = surface_stress_resolution_goal() + assert np.isclose(goal(c), c.surface_stress_resolution() * 1e6) + + def test_force_noise_density_goal_returns_finite_scalar(self): + c = _epitaxy_default() + goal = force_noise_density_goal() + value = goal(c) + assert isinstance(value, float) + assert np.isfinite(value) + assert value > 0 + # Equivalent path: resonant_force_noise_density returns fN/sqrt(Hz) + # and the goal returns pN/sqrt(Hz); they differ by a factor of 1e-3. + assert np.isclose(value, c.resonant_force_noise_density() * 1e-3) + + def test_force_noise_density_goal_drives_optimization(self): + c = _epitaxy_default() + before = force_noise_density_goal()(c) + result = optimize_performance_from_current( + c, + force_noise_density_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + assert np.isfinite(result.objective_value) + assert result.objective_value <= before + + def test_piezoelectric_voltage_goal(self): + c = _piezoelectric_default() + goal = voltage_force_resolution_goal() + assert np.isclose(goal(c), c.Fminv() * 1e12) + + def test_piezoelectric_charge_goal(self): + c = _piezoelectric_default() + goal = charge_force_resolution_goal() + assert np.isclose(goal(c), c.Fminq() * 1e12) + + +# ----------------------------------------------------------------------------- +# Constraint validation +# ----------------------------------------------------------------------------- + + +class TestConstraintValidation: + def test_metric_constraint_requires_min_or_max(self): + with pytest.raises(ValueError, match="at least one"): + CantileverMetricConstraint(CantileverMetric.STIFFNESS) + + def test_metric_constraint_min_above_max(self): + with pytest.raises(ValueError, match="<= maximum"): + CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=10.0, maximum=1.0) + + def test_metric_constraint_rejects_inf(self): + with pytest.raises(ValueError, match="finite"): + CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=float("inf")) + + +# ----------------------------------------------------------------------------- +# End-to-end optimization smoke tests +# ----------------------------------------------------------------------------- + + +class TestOptimizeFromCurrent: + def test_epitaxy_improves(self): + c = _epitaxy_default() + before = c.force_resolution() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + assert isinstance(result, OptimizationResult) + assert result.optimized is not c + # Original cantilever untouched + assert np.isclose(c.force_resolution(), before) + assert result.optimized.force_resolution() <= before + + def test_diffusion_runs(self): + c = _diffusion_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + assert np.isfinite(result.objective_value) + assert result.optimized.diffusion_time > 0 + + def test_implantation_runs(self): + c = _implantation_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + assert np.isfinite(result.objective_value) + assert result.optimized.implantation_dose > 0 + + def test_poly_runs(self): + c = _poly_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + assert np.isfinite(result.objective_value) + # Symmetric layout enforces t_bot == t_top via optimization_post_apply + assert result.optimized.t_bot == result.optimized.t_top + + def test_piezoelectric_runs(self): + c = _piezoelectric_default() + result = optimize_performance_from_current( + c, + voltage_force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=1e-3, maximum=1e2), + ], + ) + assert np.isfinite(result.objective_value) + # l_pe / w_pe must track l_si / w_si via optimization_post_apply + assert result.optimized.l_pe == result.optimized.l_si + assert result.optimized.w_pe == result.optimized.w_si + + +# ----------------------------------------------------------------------------- +# Multi-start + reproducibility +# ----------------------------------------------------------------------------- + + +class TestOptimizePerformance: + def test_multi_start_reproducible_with_seed(self): + c = _epitaxy_default() + kwargs = { + "cantilever": c, + "objective": force_resolution_goal(), + "metric_constraints": [ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + "n_starts": 2, + "max_iterations": 2, + "random_seed": 42, + } + result_a = optimize_performance(**kwargs) + result_b = optimize_performance(**kwargs) + assert np.isclose(result_a.objective_value, result_b.objective_value) + assert np.allclose(result_a.physical_state, result_b.physical_state) + + def test_multi_start_finds_better_than_single_start(self): + c = _epitaxy_default() + constraints = [ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + CantileverMetricConstraint(CantileverMetric.OMEGA_VACUUM_HZ, minimum=5000), + ] + single = optimize_performance_from_current(c, force_resolution_goal(), metric_constraints=constraints) + multi = optimize_performance( + c, + force_resolution_goal(), + metric_constraints=constraints, + n_starts=4, + max_iterations=8, + random_seed=0, + ) + # Multi-start should be at least as good as single-start (or close) + # within a generous tolerance since SLSQP is local and may converge + # to comparable basins. + assert multi.objective_value <= single.objective_value * 1.5 + + def test_n_starts_must_be_at_least_one(self): + c = _epitaxy_default() + with pytest.raises(ValueError, match="n_starts"): + optimize_performance(c, force_resolution_goal(), n_starts=0) + + def test_max_iterations_below_n_starts_rejected(self): + c = _epitaxy_default() + with pytest.raises(ValueError, match="max_iterations"): + optimize_performance( + c, + force_resolution_goal(), + n_starts=5, + max_iterations=2, + ) + + +# ----------------------------------------------------------------------------- +# Constraint enforcement +# ----------------------------------------------------------------------------- + + +class TestConstraintEnforcement: + def test_power_dissipation_respected(self): + c = _epitaxy_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=1e-3), + ], + ) + # SLSQP can violate constraints by ~1e-6 of scale; allow 1% tolerance + assert result.optimized.power_dissipation() <= 1e-3 * 1.01 + + def test_stiffness_minimum_respected(self): + c = _epitaxy_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=1e-3), + ], + ) + assert result.optimized.stiffness() >= 1e-3 * 0.99 + + def test_aspect_ratio_default_constraints_respected(self): + c = _epitaxy_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + opt = result.optimized + assert opt.l / opt.w >= 2.0 * 0.99 + assert opt.w / opt.t >= 2.0 * 0.99 + assert opt.l_pr() / opt.w_pr() >= 2.0 * 0.99 + assert opt.l_pr() >= 2e-6 * 0.99 + + def test_disabling_aspect_constraints(self): + c = _epitaxy_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + default_aspect_constraints=False, + ) + assert np.isfinite(result.objective_value) + + def test_parameter_constraint_override(self): + c = _epitaxy_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + parameter_constraints={"max_v_bridge": 3.0}, + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + assert result.optimized.v_bridge <= 3.0 + 1e-9 + + +# ----------------------------------------------------------------------------- +# Custom objectives +# ----------------------------------------------------------------------------- + + +class TestCustomObjective: + def test_lambda_objective_runs(self): + c = _epitaxy_default() + result = optimize_performance_from_current( + c, + lambda c: c.force_resolution() * 1e12 + c.power_dissipation() * 1e3, + ) + assert np.isfinite(result.objective_value) + + def test_l_bfgs_b_with_no_constraints(self): + c = _epitaxy_default() + # Disable defaults so we can use unconstrained L-BFGS-B + result = optimize_performance_from_current( + c, + force_resolution_goal(), + default_aspect_constraints=False, + method="L-BFGS-B", + ) + assert np.isfinite(result.objective_value) + + def test_l_bfgs_b_rejected_with_metric_constraints(self): + c = _epitaxy_default() + with pytest.raises(ValueError, match="L-BFGS-B"): + optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.STIFFNESS, minimum=1e-3), + ], + method="L-BFGS-B", + ) + + +# ----------------------------------------------------------------------------- +# Original cantilever immutability +# ----------------------------------------------------------------------------- + + +class TestNoiseChainAndPrintPerformance: + """Coverage for the cantilever methods that were broken pre-PR.""" + + def test_voltage_noise_accepts_scalar_and_array(self): + c = _epitaxy_default() + # Scalar input was broken by johnson_PSD's freq.size pattern + scalar_value = float(np.squeeze(c.voltage_noise(1e3))) + assert np.isfinite(scalar_value) and scalar_value > 0 + # Array input was broken by math.sqrt + array_values = np.squeeze(c.voltage_noise(np.logspace(1, 3, 5))) + assert array_values.shape == (5,) + assert np.all(np.isfinite(array_values)) + + def test_resonant_force_noise_density_returns_scalar(self): + c = _epitaxy_default() + value = c.resonant_force_noise_density() + assert isinstance(value, float) + assert np.isfinite(value) and value > 0 + + def test_k_base_returns_scalar(self): + c = _epitaxy_default() + value = c.k_base() + assert isinstance(value, float) + assert value > 0 + + def test_approx_temp_rise_returns_finite(self): + c = _epitaxy_default() + t_max, t_tip = c.approxTempRise() + assert np.isfinite(t_max) and np.isfinite(t_tip) + assert t_max >= 0 and t_tip >= 0 + + def test_print_performance_runs_end_to_end(self, capsys): + c = _epitaxy_default() + c.print_performance() + captured = capsys.readouterr() + # Spot-check that key sections rendered + assert "Force resolution" in captured.out + assert "Force noise at 1 kHz" in captured.out + assert "Approx. Temp Rises" in captured.out + + def test_temp_max_approx_constraint_works(self): + c = _epitaxy_default() + result = optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + CantileverMetricConstraint(CantileverMetric.TEMP_MAX_APPROX, maximum=20.0), + ], + ) + # The constraint should bound the optimized cantilever's max temp rise. + # SLSQP can violate by ~1% of scale. + t_max, _ = result.optimized.approxTempRise() + assert t_max <= 20.0 * 1.05 + + +def test_original_cantilever_not_mutated(): + c = _epitaxy_default() + state_vars = c.optimization_state_vars() + before = physical_state_from_cantilever(c, state_vars) + optimize_performance_from_current( + c, + force_resolution_goal(), + metric_constraints=[ + CantileverMetricConstraint(CantileverMetric.POWER_DISSIPATION, maximum=2e-3), + ], + ) + after = physical_state_from_cantilever(c, state_vars) + assert np.allclose(before, after)