diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index dd44e76..a927af7 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -44,7 +44,8 @@ pipeline orchestration, and GUI stay in their own packages. - **Geometries:** a `BaseGeometry` subclass provides `name`, `stackup_xml`, `simconfig_filename`, `input_parameter_iterator`, `create_gds_file`, and `create_dataset`; `is_feasible(params)` is optional and rejects draws before - they are drawn, and `feasibility_constraints()` states the same rules as + they are drawn, `simulation_ports(params)` is optional and changes port fields + (never their number or order) per sample, and `feasibility_constraints()` states the same rules as expressions (`geometry/constraints.py` grammar) for the ONNX metadata — a test must keep the two in agreement. Never clamp or repair parameters inside `create_gds_file` — the parameter table records the requested values, so the model would learn a diff --git a/README.md b/README.md index ac2ab8d..ed1118e 100644 --- a/README.md +++ b/README.md @@ -110,7 +110,7 @@ Each stage can be left out, e.g. to retrain on existing simulation results. An i | Stage | Class | What it does | |-------|-------|--------------| | GDS generation | `GDSGenerator` | Samples the geometry's parameters and draws a GDS layout for each | -| Design-rule check | `DRCChecker` | Snaps layouts to the manufacturing grid and drops those that violate the IHP SG13G2 rules | +| Design-rule check | `DRCChecker` | Snaps layouts to the manufacturing grid and drops those that violate the IHP SG13G2 rules or have a port marker off its metal | | GDS conversion | `GDSConverter` | Meshes the layouts for Palace with [gds2palace](https://github.com/VolkerMuehlhaus/gds2palace_ihp_sg13g2) | | EM simulation | `PalaceSimulator` | Runs Palace and stores the S-parameters as Touchstone files | | Model training | `ModelTrainer` | Trains (and by default tunes) a PyTorch model from geometry and frequency to S-parameters | diff --git a/docs/custom_class.md b/docs/custom_class.md index 35cadfd..e7ed2a5 100644 --- a/docs/custom_class.md +++ b/docs/custom_class.md @@ -55,8 +55,10 @@ The Python class should be a `@dataclass` extending `orca.BaseGeometry` and must **Optional methods:** +- `electrical_parameters(ntwk) -> dict[str, np.ndarray]` — The figures of merit of a network of this geometry, which `ModelTester` reports the model's error for next to the S-parameter error: curves over `ntwk.f` (such as L, R and Q) and scalars as 0-d arrays (such as the self-resonance frequency `srf_f`). Which port is which is a property of the layout, so the geometry decides how to read its network. `orca.utils.postprocessing` has the usual ones, with the ports given explicitly: `inductor_parameters(ntwk, ends, shorted)` and `transformer_parameters(ntwk, primary, secondary, shorted)`, which drive each winding differentially with its center tap (`shorted`) AC-grounded. Import it inside the method. The default returns nothing, so only the S-parameter error is reported. Errors are relative, except for the names in the class attribute `absolute_error_parameters` (e.g. `frozenset({"k"})` for a coupling factor that is near zero for weak coupling). With `ModelTester(plot=True)` the same parameters are plotted, reference against prediction. - `feasibility_constraints() -> list[str]` — The same rules as `is_feasible()`, written as boolean expressions over the input parameter names, e.g. `"bottom_linewidth <= bottom_winding_diameter / 3"`. `OnnxExporter` stores them in the model's `input_constraints` metadata, so COBRA can refuse a query for a geometry that cannot be built instead of returning a prediction the model was never trained for. The grammar is a small subset of Python (arithmetic, comparisons, `and`/`or`/`not`, `a if c else b`, and `abs min max sqrt sin cos tan radians ceil floor round`, plus `pi` and `sqrt2`), documented in `orca.geometry.constraints`; anything else is rejected at export. Derive the strings from the same numbers as `is_feasible()` and add a test that they agree on random draws (see `tests/test_constraints.py` for the presets' version). - `is_feasible(params) -> bool` — Whether a parameter combination describes a layout that can be drawn (default: always `True`). `GDSGenerator` calls it for every draw of the iterator; rejected draws are counted and, with the `"sobol"`, `"lhs"` and `"random"` strategies, replaced by further draws, so the requested number of samples is met with buildable layouts only. Put cheap, closed-form constraints between parameters here — a winding that must fit its diameter, a feed gap that must fit the octagon's side. The presets derive it from the same check their cell code runs, so the two cannot disagree. +- `simulation_ports(params) -> list[dict]` — The gds2palace ports of one sample, as entries of the simconfig's `ports` list (default: the simconfig's ports unchanged). Which metal a port has to reach is a property of the layout, so it may depend on the parameters: `InductorOcta` feeds a single turn on TopMetal2 and more turns on TopMetal1, so it returns ports 1/2 with `"to_layername": "TopMetal2"` for `turns == 1`. Change port fields only; the simulation settings stay those of the simconfig, so every sample of one model is meshed and solved alike. Keep every port, with its number and in its order — the Touchstone files and the model's outputs depend on them — which `ports_for(params)`, the method the stages call, checks. `DRCChecker` checks each layout against its own ports and `GDSConverter` hands them to gds2palace. !!! warning "Reject, never clamp" @@ -114,7 +116,7 @@ class TransformerOcta(BaseGeometry): """ name: str = "tf_octa_c_ports" - stackup_xml: str = StackupXML.SG13G2_FEM_200um # from orca.geometry.presets + stackup_xml: str = StackupXML.SG13G2_FEM_200um_passi3D # from orca.geometry.presets simconfig_filename: str = os.path.join(os.path.dirname(__file__), "tf_octa_c_ports.simcfg") input_parameter_iterator: InputParameterIterator = field(default_factory=_input_parameters) @@ -132,6 +134,16 @@ class TransformerOcta(BaseGeometry): output_normalizer=StandardNormalizer(), ) + # Close to zero for weakly coupled windings, where a relative error says nothing + absolute_error_parameters: ClassVar[frozenset[str]] = frozenset({"k"}) + + def electrical_parameters(self, ntwk: "rf.Network") -> dict[str, "np.ndarray"]: + from orca.utils.postprocessing import transformer_parameters + + # Zero-based ports: 0 op, 1 on (top winding), 2 ip, 3 in (bottom winding), + # 4 and 5 the center taps, AC-grounded as in differential use + return transformer_parameters(ntwk, primary=(0, 1), secondary=(2, 3), shorted=(4, 5)) + @staticmethod def create_gds_file(name: str, output_path: str, params: dict[str, Any]) -> str: # @@ -333,3 +345,5 @@ Example: ] } ``` + +The `ports` list is the default for every sample. When the metal a port lands on changes with the parameters, override `simulation_ports(params)` in the geometry class rather than keeping a second simcfg: one file of simulation settings cannot drift apart between samples of the same model. Each port's marker on its `source_layernum` must touch the metals it connects; `DRCChecker` drops layouts where it does not, as such a port would simulate as an open circuit. diff --git a/docs/index.md b/docs/index.md index ffbf5bd..fa6181a 100644 --- a/docs/index.md +++ b/docs/index.md @@ -82,7 +82,7 @@ flowchart LR | Stage | Purpose | |---|---| | `GDSGenerator` | Samples geometry parameters and writes GDS layout files | -| `DRCChecker` | Snaps layouts to the manufacturing grid and drops those violating the SG13G2 design rules | +| `DRCChecker` | Snaps layouts to the manufacturing grid and drops those violating the SG13G2 design rules or with port markers off their metal | | `GDSConverter` | Converts GDS files to Palace-compatible mesh inputs | | `PalaceSimulator` | Runs full-wave EM simulations and stores Touchstone results | | `ModelTrainer` | Trains a PyTorch neural network on the simulation dataset | diff --git a/docs/pipeline.md b/docs/pipeline.md index 0816baa..3652af3 100644 --- a/docs/pipeline.md +++ b/docs/pipeline.md @@ -13,7 +13,7 @@ ORCA runs a linear pipeline. Each stage receives a context dictionary and adds i ```mermaid flowchart TB - A[GDSGenerator] -- GDS files --> X["DRCChecker
(grid snap + SG13G2 rules)"] + A[GDSGenerator] -- GDS files --> X["DRCChecker
(grid snap + SG13G2 rules + port contact)"] X -- clean layouts --> B["GDSConverter
(gds2palace mesh)"] B -- mesh files --> C["PalaceSimulator
(full-wave EM)"] C -- "Touchstone .sNp" --> D["ModelTrainer
(PyTorch MLP)"] @@ -33,11 +33,15 @@ Every generated layout is snapped to the manufacturing grid and checked against Off-grid vertices are repaired rather than reported: `snap_to_grid=True` (the default) moves every vertex onto the `grid_nm` grid (5 nm for SG13G2) and writes the GDS file back, so geometry code does not need its own snapping. Layouts with remaining violations are left out of the parameter table the later stages use (`drop_violations=True`); the stage writes `_drc_report.csv` with the per-layout counts and `_drc.csv` with the layouts that passed, both next to the GDS files, and logs a summary. This is where parameter combinations that draw unbuildable geometry — a winding folded over itself, a via clipped by a miter — are stopped before they cost simulation time or teach the model shapes that cannot be fabricated. +The stage also checks every layout against its ports (`check_ports=True`, the default) — the `.simcfg`'s, or the sample's own if the geometry overrides `simulation_ports(params)`: each port's marker — the shape on its `source_layernum` — has to overlap or touch the metals the port connects (`from_layername` and `to_layername`, or `target_layername` for an in-plane port, with GDS numbers from the stackup XML). gds2palace places a port wherever its marker is drawn and does not check that it meets a conductor, so a marker a few nanometres off its feed line still meshes and simulates, as an open circuit: the S-parameters then describe a different circuit than the parameter table says, and nothing fails. The report lists such layouts as `port.missing` or `port.open_`, and they are left out of the parameter table even with `drop_violations=False`. + ## Stage 3 — GDS conversion (`GDSConverter`) -Each GDS file — those that passed DRC when the stage ran, otherwise all of them — is converted to a Palace-ready simulation setup using [gds2palace](https://github.com/VolkerMuehlhaus/gds2palace_ihp_sg13g2). The geometry's `stackup_xml` defines the physical layer stackup and material properties; the `simconfig_filename` defines the simulation parameters (port positions, frequency sweep, mesh settings). +Each GDS file — those that passed DRC when the stage ran, otherwise all of them — is converted to a Palace-ready simulation setup using [gds2palace](https://github.com/VolkerMuehlhaus/gds2palace_ihp_sg13g2). The geometry's `stackup_xml` defines the physical layer stackup and material properties; the `simconfig_filename` defines the simulation parameters (ports, frequency sweep, mesh settings). Each sample gets the ports of the geometry's `simulation_ports(params)`, by default the simconfig's; `InductorOcta` uses this to land a single turn's ports on TopMetal2, where its feeds are. + +Conversions run in parallel worker processes, one sample per task. A sample whose geometry gds2palace/gmsh cannot mesh (e.g. `PLC Error: A segment and a facet intersect`) is logged and skipped. gmsh can also loop forever on degenerate geometry, so each conversion has a time limit — `GDSConverter(timeout=180)` seconds by default, which leaves room for fine meshes (a large inductor at `refined_cellsize` 2 µm takes up to about 50 s) — after which its worker is killed and the sample is skipped as well, instead of stalling the whole pipeline. Models converted by an earlier run for the same parameters are reused. -Conversions run in parallel worker processes, one sample per task. A sample whose geometry gds2palace/gmsh cannot mesh (e.g. `PLC Error: A segment and a facet intersect`) is logged and skipped. gmsh can also loop forever on degenerate geometry, so each conversion has a time limit — `GDSConverter(timeout=60)` seconds by default — after which its worker is killed and the sample is skipped as well, instead of stalling the whole pipeline. Models converted by an earlier run for the same parameters are reused. +Every new mesh is also checked for flat or inverted tetrahedra. gmsh occasionally produces such zero-volume elements where thin layers meet in one plane — a 0.85 µm port sheet flush with the ground ring's outer edge, as in both presets, at a coarse `refined_cellsize` — and Palace cannot solve such a mesh (`HasPositiveFiniteDiagonal(...) is false`). Rarely (about 1 in 300 small layouts with the `passi3D` stackup), gmsh also writes a corrupt mesh whose elements reference a node that does not exist; it cannot be read back and counts as quality −∞. A sample whose worst element quality (gmsh's `minSICN`: 1 for a regular tetrahedron, 0 for a flat one) is at or below `min_element_quality` (default `1e-6`) is left out of the Palace table, and the stage warns that the simcfg's `refined_cellsize` may be too coarse. Meshes below `warn_element_quality` (default `1e-4`) are kept but listed in a warning; meshes Palace has solved had a worst quality of 6·10⁻⁴ and above. The worst quality of every mesh is written to `palace_sims/_mesh_report.csv`. ## Stage 4 — EM simulation (`PalaceSimulator`) @@ -56,7 +60,7 @@ Each simulation that finishes adds its row to `results/.csv` immediately, A PyTorch MLP is trained on the simulation data. Inputs are geometry parameters and frequency; outputs are the real and imaginary parts of each S-parameter entry. Normalization is defined in the geometry's dataset and applied automatically. An optional basis expansion of the inputs — for example a Chebyshev expansion of frequency — is chosen on the stage itself with `ModelTrainer(basis="chebyshev")`; it lives inside the model, so it is tuned with it and exported into the ONNX graph. Hyperparameters such as learning rate, batch size, and network depth can be passed to `ModelTrainer`. -The result table is split by geometry, never by frequency point: `test_frac` of the geometries are held out for `ModelTester`, and `val_frac` of the rest select the best checkpoint during training. The split is recorded in `models/_split.csv`. Without `hyperparameters`, Optuna tunes them with `n_fold_cv`-fold cross-validation over the geometries, so every fold is scored on layouts the model has not seen (`n_trials` trials, or until `tuning_timeout` seconds have passed). `n_fold_cv=1` skips the cross-validation and scores each trial on the `val_frac` split of the final training instead: one training per trial rather than `n_fold_cv`, which is usually enough once there are a few thousand geometries; with a few hundred, the average over folds ranks trials more reliably. The number of epochs is not tuned: each fold runs up to `tuning_max_epochs` with early stopping, the final model up to `max_epochs`, and a trial that falls behind the others at the same fold and epoch is pruned after any epoch. The batch sizes searched are `batch_sizes` (32 to 512 if not given); for a per-point dataset with millions of samples, larger ones such as `[1024, 2048, 4096]` train much faster. A trial that diverges or runs out of memory is pruned; any other error stops the stage. Every training, in tuning and in the final run, starts with a linear learning-rate warmup over `warmup_epochs` (default 1) and then follows `lr_schedule`: `"cosine"` (default) decays the rate smoothly to 1% of its tuned value over the epoch limit, `"plateau"` halves it whenever the validation loss stalls. Gradients are clipped to a total norm of `grad_clip_norm` (default 1.0; `None` disables it), so one bad batch cannot undo the progress so far. With `allow_tf32=True`, tuning and the final training let matrix multiplies round their inputs to TF32 (10 instead of 23 mantissa bits) on GPUs with TF32 tensor cores (Ampere and newer). Testing and the exported ONNX model always run in full FP32. The MLP's width is searched on a log scale from 64 to 2048 and its depth from 2 to 8 layers, so small networks are tried as often as large ones. With `regularization=True` the weight decay and the model's own regularization (dropout for the MLP) are tuned too; otherwise AdamW's default weight decay and no dropout are used. The hyperparameters a run trained with are saved to `models/_hyperparameters.json`, and `hyperparameters` accepts the path of such a file as well as a dict, so a later run can retrain with them without tuning again. Each Touchstone file is parsed once per run, for tuning and the final training together. The run seed (`ORCA.run(seed=...)`) fixes the splits, the tuner and the weight initialisation, so two runs with the same seed and data train the same model. +The result table is split by geometry, never by frequency point: `test_frac` of the geometries are held out for `ModelTester`, and `val_frac` of the rest select the best checkpoint during training. The split is recorded in `models/_split.csv`. Without `hyperparameters`, Optuna tunes them with `n_fold_cv`-fold cross-validation over the geometries, so every fold is scored on layouts the model has not seen (`n_trials` trials, or until `tuning_timeout` seconds have passed). `n_fold_cv=1` skips the cross-validation and scores each trial on the `val_frac` split of the final training instead: one training per trial rather than `n_fold_cv`, which is usually enough once there are a few thousand geometries; with a few hundred, the average over folds ranks trials more reliably. The number of epochs is not tuned: each fold runs up to `tuning_max_epochs` with early stopping, the final model up to `max_epochs`, and a trial that falls behind the others at the same fold and epoch is pruned after any epoch. The batch sizes searched are `batch_sizes` (32 to 512 if not given); for a per-point dataset with millions of samples, larger ones such as `[1024, 2048, 4096]` train much faster. A trial that diverges or runs out of memory is pruned; any other error stops the stage. Every training, in tuning and in the final run, starts with a linear learning-rate warmup over `warmup_epochs` (default 1) and then follows `lr_schedule`: `"cosine"` (default) decays the rate smoothly to 1% of its tuned value over the epoch limit, `"plateau"` halves it whenever the validation loss stalls. Gradients are clipped to a total norm of `grad_clip_norm` (default 1.0; `None` disables it), so one bad batch cannot undo the progress so far. With `allow_tf32=True`, tuning and the final training let matrix multiplies round their inputs to TF32 (10 instead of 23 mantissa bits) on GPUs with TF32 tensor cores (Ampere and newer). Testing and the exported ONNX model always run in full FP32. The MLP's width is searched on a log scale from 64 to 2048 and its depth from 2 to 8 layers, so small networks are tried as often as large ones. With `regularization=True` the weight decay and the model's own regularization (dropout for the MLP) are tuned too; otherwise AdamW's default weight decay and no dropout are used. Three optional loss settings judge the prediction as a circuit, not only as a vector of normalized S-parameters. They need a per-point dataset, as the presets have. `admittance_weight` adds the relative error of the admittance matrix Y, as a whole and of its real part: a small error in S can be a large error in L and Q, because at low frequency an inductor is almost a short between its ports. 0.1 to 1 puts it on the scale of the S-parameter loss. `passivity_weight` penalises predictions with a largest singular value of S above 1, or above the simulated one where that is larger (the DC point that `dc_deembedded` extrapolates is slightly non-passive). `above_srf_weight` weights the frequency points above each geometry's first self-resonance relative to those below it. The resonance is found in the simulated S-parameters as the frequency where the first inductive mode turns capacitive, which for the center-tapped inductor is its differential self-resonance. In the InductorOcta results (0 to 500 GHz) about 95% of the frequency points lie above it, so an `above_srf_weight` of 0.1 still leaves the points below with about a third of the loss. All three default to off, and they apply to tuning and the final training alike; validation losses are only comparable between runs with the same settings. The hyperparameters a run trained with are saved to `models/_hyperparameters.json`, and `hyperparameters` accepts the path of such a file as well as a dict, so a later run can retrain with them without tuning again. Each Touchstone file is parsed once per run, for tuning and the final training together. The run seed (`ORCA.run(seed=...)`) fixes the splits, the tuner and the weight initialisation, so two runs with the same seed and data train the same model. ## Stage 6 — ONNX export (`OnnxExporter`) @@ -64,4 +68,4 @@ The trained PyTorch model is exported to ONNX format with a fixed frequency swee ## Stage 7 — Model testing (`ModelTester`) -The trained model (or, if training did not run in this pipeline, the exported ONNX model) is evaluated against the held-out geometries listed in `models/_split.csv`. Without that file, for example for a model tested against a fresh results folder, every row of the result table is used and a warning says so. Besides the mean absolute S-parameter error and the median relative error of each electrical parameter, the stage reports the spread: the median, 95th percentile and worst geometry, the error in each of `n_frequency_bands` frequency bands, and the 95th percentile of each electrical parameter's error. The errors of every test geometry, next to its parameters, are written to `models/_test_errors.csv`, for example to plot the error against each parameter and find under-sampled regions. The error is also resolved over frequency: `models/_errors_vs_frequency.png` shows, for the S-parameters and each electrical parameter, the median and the 25th–75th and 5th–95th percentiles over the test geometries at every frequency point, so you can see which frequency ranges the model gets right; the values are in `models/_errors_vs_frequency.csv`. The coupling factor k is reported as an absolute error, since a relative one explodes for weakly coupled layouts. Prediction errors are logged to help assess whether the surrogate is accurate enough for use in [COBRA](https://github.com/DI-PASSIONATE/COBRA). +The trained model (or, if training did not run in this pipeline, the exported ONNX model) is evaluated against the held-out geometries listed in `models/_split.csv`. Without that file, for example for a model tested against a fresh results folder, every row of the result table is used and a warning says so. Besides the mean absolute S-parameter error and the median relative error of each electrical parameter, the stage reports the spread: the median, 95th percentile and worst geometry, the error in each of `n_frequency_bands` frequency bands, and the 95th percentile of each electrical parameter's error. The errors of every test geometry, next to its parameters, are written to `models/_test_errors.csv`, for example to plot the error against each parameter and find under-sampled regions. The error is also resolved over frequency: `models/_errors_vs_frequency.png` shows, for the S-parameters and each electrical parameter, the median and the 25th–75th and 5th–95th percentiles over the test geometries at every frequency point, so you can see which frequency ranges the model gets right; the values are in `models/_errors_vs_frequency.csv`. The electrical parameters are the geometry's own (`electrical_parameters()`, see [Custom Classes](custom_class.md)), since only the geometry knows which port is which. `InductorOcta` reports its differential inductance `L`, resistance `R`, quality factor `Q` and self-resonance `srf_f`, with the center tap AC-grounded as in differential operation; `TransformerOcta` reports `Lp`, `Ls`, `Rp`, `Rs`, `Qp`, `Qs`, the coupling factor `k` and the primary's self-resonance, each winding driven differentially with both center taps AC-grounded. A geometry that declares none is tested on its S-parameters only. The coupling factor k is reported as an absolute error, since a relative one explodes for weakly coupled layouts; a geometry lists such parameters in `absolute_error_parameters`. Prediction errors are logged to help assess whether the surrogate is accurate enough for use in [COBRA](https://github.com/DI-PASSIONATE/COBRA). diff --git a/docs/running_orca.md b/docs/running_orca.md index 3358d93..23e91ac 100644 --- a/docs/running_orca.md +++ b/docs/running_orca.md @@ -73,6 +73,9 @@ orca_instance = ORCA( # warmup_epochs=1.0, # Linear learning-rate warmup at the start of training # grad_clip_norm=1.0, # Gradient-norm clipping; None disables it # allow_tf32=False, # TF32 matmuls while training (faster on A100 and newer) + # admittance_weight=0.0, # Loss on the relative Y-parameter error (tracks L and Q) + # passivity_weight=0.0, # Penalty on predicted S with a singular value above 1 + # above_srf_weight=1.0, # Loss weight of frequencies above each self-resonance ), orca.OnnxExporter(), orca.ModelTester(), diff --git a/examples/main.py b/examples/main.py index 3a31ba4..4c4c0a3 100644 --- a/examples/main.py +++ b/examples/main.py @@ -1,7 +1,5 @@ import orca -from orca.geometry.presets import TransformerOcta - -# from orca.geometry.presets import InductorOcta +from orca.geometry.presets import InductorOcta, TransformerOcta PLOT = False @@ -21,6 +19,7 @@ def main(): # Use predefined geometry from examples geometry = TransformerOcta(name="transformer_octa") + geometry = InductorOcta(name="inductor_octa") orca_instance = orca.ORCA( [ diff --git a/examples/slurm_runs/orca_slurm_train.sh b/examples/slurm_runs/orca_slurm_train.sh index e929caa..0aa9c7f 100644 --- a/examples/slurm_runs/orca_slurm_train.sh +++ b/examples/slurm_runs/orca_slurm_train.sh @@ -14,7 +14,7 @@ # --gres=gpu:rtx3080:1 --partition=rtx3080 or --gres=gpu:v100:1 --partition=v100 #SBATCH --gres=gpu:a100:1 #SBATCH --partition=a100 -#SBATCH --time=24:00:00 +#SBATCH --time=12:00:00 #SBATCH --job-name=ORCA-TRAIN #SBATCH --export=NONE diff --git a/examples/slurm_runs/train.py b/examples/slurm_runs/train.py index 9ae8193..0055922 100644 --- a/examples/slurm_runs/train.py +++ b/examples/slurm_runs/train.py @@ -19,19 +19,22 @@ def main(): orca_instance = orca.ORCA( [ orca.ModelTrainer( - n_fold_cv=3, - n_trials=100, + n_fold_cv=1, + n_trials=200, batch_sizes=[1024, 2048, 4096, 8192, 16384], tuning_max_epochs=30, basis="chebyshev", - # Stop tuning after 16 h, leaving the rest of the 24 h limit for the final + # Stop tuning after 10 h, leaving the rest of the 24 h limit for the final # training, the export and the test. A trial still running at that point is # pruned after its current epoch. - tuning_timeout=16 * 3600, - max_epochs=200, + tuning_timeout=10 * 3600, + max_epochs=300, # TF32 matrix multiplies while training: faster on the A100, ignored on the # V100. Testing and the exported model stay in FP32. allow_tf32=True, + admittance_weight=0.5, + above_srf_weight=0.6, + passivity_weight=1.0, ), orca.OnnxExporter(), orca.ModelTester(), diff --git a/src/orca/__init__.py b/src/orca/__init__.py index ca7a044..96df8d6 100644 --- a/src/orca/__init__.py +++ b/src/orca/__init__.py @@ -62,6 +62,7 @@ "HyperparameterTuner": ".training.tuner", "ComplexMSELoss": ".training.losses", "MSEPlusLogCoshLoss": ".training.losses", + "SParameterLoss": ".training.losses", "NetworkPredictor": ".training.predictors", "TorchNetworkPredictor": ".training.predictors", "OnnxNetworkPredictor": ".training.predictors", diff --git a/src/orca/geometry/base_geometry.py b/src/orca/geometry/base_geometry.py index 7760c6c..e10cc2a 100644 --- a/src/orca/geometry/base_geometry.py +++ b/src/orca/geometry/base_geometry.py @@ -1,11 +1,15 @@ from __future__ import annotations +import copy from abc import ABC, abstractmethod from dataclasses import dataclass from functools import cached_property -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING, Any, ClassVar if TYPE_CHECKING: + import numpy as np + import skrf as rf + from orca.geometry.input_parameters import InputParameterIterator from orca.training.datasets.base_dataset import BaseDataset @@ -22,6 +26,10 @@ class BaseGeometry(ABC): simconfig_filename: str input_parameter_iterator: InputParameterIterator + #: Electrical parameters whose test error is an absolute difference rather than a + #: relative one, e.g. a coupling factor that is close to zero for weak coupling. + absolute_error_parameters: ClassVar[frozenset[str]] = frozenset() + @property def input_iterator(self) -> InputParameterIterator: # Return the input parameter iterator, ensuring it is initialized with iter() @@ -44,9 +52,54 @@ def n_ports(self) -> int: the simulation config. This sets the ``.sNp`` extension of the Touchstone results and must match the codec of :attr:`dataset`. """ + return len(self._simconfig_ports) + + @cached_property + def _simconfig_ports(self) -> list[dict[str, Any]]: from orca.simulation.simulate import read_simconfig - return len(read_simconfig(self.simconfig_filename)["ports"]) + return read_simconfig(self.simconfig_filename)["ports"] + + def simulation_ports(self, params: dict[str, Any]) -> list[dict[str, Any]]: # noqa: ARG002 - hook with a default + """ + The gds2palace ports of one sample, as entries of the simconfig's ``ports`` list. + + Which metal a port has to reach is a property of the layout, so it may depend + on the parameters: a single-turn inductor feeds on another metal than a + multi-turn one. Override this to change port fields per sample, typically + ``to_layername``; the simulation settings stay those of the simconfig, so + every sample of one model is meshed and solved alike. Keep every port, with + its number and in its order: the Touchstone files and the model's outputs + depend on them, and :meth:`ports_for` checks it. + + The default returns the simconfig's ports unchanged. + + Args: + params (dict[str, Any]): The sample's input parameters, as in the + parameter table (integers may come back as floats). + + Returns: + list[dict[str, Any]]: A fresh copy of the port entries, safe to modify. + """ + return copy.deepcopy(self._simconfig_ports) + + def ports_for(self, params: dict[str, Any]) -> list[dict[str, Any]]: + """ + :meth:`simulation_ports` for one sample, checked against the simconfig. + + Raises: + ValueError: The ports are not the simconfig's port numbers in its order. + """ + ports = self.simulation_ports(params) + expected = [port["portnumber"] for port in self._simconfig_ports] + numbers = [port.get("portnumber") for port in ports] + if numbers != expected: + raise ValueError( + f"{type(self).__name__}.simulation_ports returned ports {numbers} for " + f"{params}; the simconfig defines {expected}, and every sample must keep " + "those numbers in that order." + ) + return ports def is_feasible(self, params: dict[str, Any]) -> bool: # noqa: ARG002 - hook with a default """ @@ -84,6 +137,30 @@ def feasibility_constraints(self) -> list[str]: """ return [] + def electrical_parameters(self, ntwk: rf.Network) -> dict[str, np.ndarray]: # noqa: ARG002 - hook with a default + """ + The figures of merit of a simulated or predicted network, which ``ModelTester`` + reports the model's error for, next to the S-parameter error. + + Which port is which is a property of the layout, so the geometry decides how to + read its network: for example the differential L, R and Q of an inductor whose + center tap is AC-grounded. :mod:`orca.utils.postprocessing` has the usual ones + (:func:`~orca.utils.postprocessing.inductor_parameters`, + :func:`~orca.utils.postprocessing.transformer_parameters`). Import it inside + this method, as the presets do. + + The default returns nothing, so only the S-parameter error is reported. + + Args: + ntwk (rf.Network): The network, simulated or predicted. + + Returns: + dict: Parameter name to a curve over ``ntwk.f`` or to a scalar (a 0-d + array), such as a self-resonance frequency. Errors are relative unless the + name is in :attr:`absolute_error_parameters`. + """ + return {} + @staticmethod @abstractmethod def create_gds_file(name: str, output_path: str, params: dict[str, Any]) -> str: diff --git a/src/orca/geometry/cells/inductor.py b/src/orca/geometry/cells/inductor.py index bed849c..af6c63e 100644 --- a/src/orca/geometry/cells/inductor.py +++ b/src/orca/geometry/cells/inductor.py @@ -22,8 +22,8 @@ TopMetal2 (134) spiral windings TopMetal1 (126) crossovers + feedlines TopVia2 (133) vias - Metal1 (8) ground frame (only when forEM=True; ``ground_layer`` picks - another layer, ``ground_style`` a strip instead of a ring) + Metal1 (8) ground frame, and ground under the feeds or port tabs (only + when forEM=True; ``ground_layer`` picks another layer) 201/202/203 EM ports (only when forEM=True) Standalone GDS build (just gdspy + matplotlib): @@ -85,8 +85,6 @@ VIA_MARGIN = 0.5 # IHP TopVia2 rule TV2.c, TV2.d DELTA = 0.1 # size of EM port perpendicular to width -GROUND_STRIP_OVERLAP = 2.0 # ground strip reaches this far past the feed ends -GROUND_SLOT_WIDTH = 10.0 # gap cut into the right bar of a "slotted_ring" ground MU0 = 4 * math.pi * 1e-7 @@ -212,31 +210,20 @@ def calculate_octa_diameter(N, w, s, Ltarget, K1=2.15522, K2=3.61868, L0=0): # ==================== def symmetric_octa_IHP(N, D, w, s, includeCenterTap=False, LBE=False, forEM=False, - include_nofill=True, ground_layer=None, ring_spacing=None, - ring_width=None, ground_style="ring", filename="inductor.gds", + include_nofill=True, ground_layer=FRAME_LAYER_NUM, ring_spacing=None, + ring_width=None, feeds_to_ring_edge=False, filename="inductor.gds", textlabel=""): # Drawing unit and parameter unit is micron - # ground_layer: GDS layer number for the EM ground (forEM=True). - # None -> FRAME_LAYER_NUM (8 = Metal1). Use 67 for Metal5, 250 for SUBGND. - # ground_style: "ring" draws a closed frame around the inductor with the - # ports on its outer edge. "strip" draws one plate below the feed ends - # only, so no closed loop surrounds the spiral; use it with a lossless - # ground layer (SUBGND) as the common port reference. "slotted_ring" draws - # the ring with a GROUND_SLOT_WIDTH gap in its right bar, so it is not a - # closed loop either, and the ports sit on its inner part like on the - # strip; use it when ports are on opposite sides (a center tap at the top - # for odd N) and need one connected ground. - # ring_spacing: gap [um] from the inductor outer radius (D/2) to the inner - # edge of the ground (forEM=True). None -> D/2 (original behaviour, - # scales with diameter). Set e.g. 10 or 20 for a fixed clearance. - # ring_width: thickness [um] of the ground-ring frame, or depth of the - # ground strip below the feed ends (forEM=True). - # None -> min(20, 5*w) (original behaviour). Set e.g. 10 for a fixed width. - if ground_style not in ("ring", "strip", "slotted_ring"): - raise ValueError( - f"ground_style must be 'ring', 'strip' or 'slotted_ring', not {ground_style!r}") - if ground_layer is None: - ground_layer = FRAME_LAYER_NUM + # ORCA additions (all keep upstream's layout at their defaults): + # include_nofill: also draw the OPDK nofill octagons when not forEM. + # ground_layer: GDS layer of the EM ground frame (upstream: FRAME_LAYER_NUM, Metal1). + # ring_spacing: gap [um] from the outer diameter to the frame (forEM). None: D/2 as + # upstream, but at least past the 30 um feeds (see the frame below). + # ring_width: width [um] of the frame bars (forEM). None: min(20, 5 w) as upstream. + # feeds_to_ring_edge: run the feeds out to the frame's outer edge and put the ports + # there (forEM), so the reference planes are the cell's boundary and layouts can + # be placed side by side; instead of upstream's 30 um feeds over the ground under + # the feedline. # GDSII setup lib = gdspy.GdsLibrary() @@ -246,10 +233,12 @@ def symmetric_octa_IHP(N, D, w, s, includeCenterTap=False, LBE=False, forEM=Fals else: cellname = f"inductor2_N{N}_Do{D}_w{w}_s{s}" - try: - cell = lib.new_cell(cellname, overwrite_duplicate=True) - except ValueError: - cell = lib.new_cell("final_" + cellname, overwrite_duplicate=True) + # ORCA fix: upstream's lib.new_cell also registers the cell in gdspy's process-wide + # library, so drawing the same parameters a third time in one process (a reused + # worker meeting a repeated corner draw) failed even after its "final_" fallback. + # The cell only needs to live in this function's own library. + cell = gdspy.Cell(cellname, exclude_from_current=True) + lib.add(cell) # list with all geometries that we created all_geometries_list = [] @@ -279,19 +268,22 @@ def symmetric_octa_IHP(N, D, w, s, includeCenterTap=False, LBE=False, forEM=Fals # Inner diameter Di = gridsnap(D - 2 * N * w - 2 * (N - 1) * s) - # Ground-ring geometry (only present when forEM). Computed up front so the - # feed length can reach the ring. + # Ground frame (forEM), computed up front because the feeds may run out to it frame_width = min(20, gridsnap(5 * w)) if ring_width is None else gridsnap(ring_width) - frame_margin = gridsnap(D / 2) if ring_spacing is None else gridsnap(ring_spacing) - - # Feed length: when forEM, extend the feedlines so the pins/ports always - # land on the OUTER edge of the ground ring, or just inside the ground - # strip / slotted ring; otherwise keep the default. - if not forEM: - feed_length = 30 - elif ground_style in ("strip", "slotted_ring"): - feed_length = gridsnap(frame_margin + GROUND_STRIP_OVERLAP) + if ring_spacing is not None: + frame_margin = gridsnap(ring_spacing) + elif feeds_to_ring_edge: + frame_margin = gridsnap(D / 2) else: + # ORCA fix: upstream places the frame D/2 outside the spiral. Below D = 64 um + # the 30 um feeds then cross the frame, so the ports end up over it or beyond it + # (open), and how much feed runs over ground depends on D. Keeping the frame + # 2 um beyond the feed ends gives every size upstream's large-D layout. + frame_margin = max(gridsnap(D / 2), 30 + 2) + + # Feed length + feed_length = 30 + if forEM and feeds_to_ring_edge: feed_length = gridsnap(frame_margin + frame_width) # --- Feedline drawing --- @@ -621,17 +613,8 @@ def symmetric_octa_IHP(N, D, w, s, includeCenterTap=False, LBE=False, forEM=Fals if LBE: add_poly(all_geometries_list, layer=LBE_LAYER_NUM, purpose=PURPOSE_DRAWING, points=points) - # --- ground for EM simulation using gds2palace ------- - # (frame_width / frame_margin were computed up front, near the feed length) - if forEM and ground_style == "strip": - # One plate below the feed ends, as in the gds2palace L6n2 study. It - # reaches GROUND_STRIP_OVERLAP past the feed ends so the port boxes sit - # on it, and its inner edge is frame_margin from the outer diameter. - y_feed_end = y0 - D / 2 - feed_length - add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, - p1=(gridsnap(x0 - D / 2), gridsnap(y_feed_end + GROUND_STRIP_OVERLAP)), - p2=(gridsnap(x0 + D / 2), gridsnap(y_feed_end - frame_width))) - elif forEM: + # --- ground frame for EM simulation using gds2palace ------- + if forEM: xmin_frame_inner = gridsnap(x0 - D / 2 - frame_margin) xmax_frame_inner = gridsnap(x0 + D / 2 + frame_margin) ymin_frame_inner = gridsnap(y0 - D / 2 - frame_margin) @@ -645,18 +628,9 @@ def symmetric_octa_IHP(N, D, w, s, includeCenterTap=False, LBE=False, forEM=Fals add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, p1=(xmin_frame_outer, ymin_frame_outer), p2=(xmin_frame_inner, ymax_frame_outer)) - if ground_style == "slotted_ring": - # right bar in two pieces, so the frame is not a closed loop - add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, - p1=(xmax_frame_inner, ymin_frame_outer), - p2=(xmax_frame_outer, gridsnap(y0 - GROUND_SLOT_WIDTH / 2))) - add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, - p1=(xmax_frame_inner, gridsnap(y0 + GROUND_SLOT_WIDTH / 2)), - p2=(xmax_frame_outer, ymax_frame_outer)) - else: - add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, - p1=(xmax_frame_inner, ymin_frame_outer), - p2=(xmax_frame_outer, ymax_frame_outer)) + add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, + p1=(xmax_frame_inner, ymin_frame_outer), + p2=(xmax_frame_outer, ymax_frame_outer)) add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, p1=(xmin_frame_inner, ymin_frame_inner), p2=(xmax_frame_inner, ymin_frame_outer)) @@ -664,12 +638,22 @@ def symmetric_octa_IHP(N, D, w, s, includeCenterTap=False, LBE=False, forEM=Fals p1=(xmin_frame_inner, ymax_frame_inner), p2=(xmax_frame_inner, ymax_frame_outer)) - # NOTE: the original gds2palace code added a "ground under feedline" - # rectangle here (filling the feed opening down to the ring). With the - # feed length now reaching the ring it overlapped the ring frame, so it - # is intentionally omitted — the ring is just the clean square frame. - - # add all created shapes to cell now + if not feeds_to_ring_edge: + # ground under feedline at pin LA,LB + add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, + p1=(x0 - feedline_spacing / 2 - w, y0 - D / 2 - feed_length + 2), + p2=(x0 + feedline_spacing / 2 + w, ymin_frame_inner)) + + if includeCenterTap and not is_even(N): + # ground under feedline at top side pin LC + add_box(all_geometries_list, layer=ground_layer, purpose=PURPOSE_DRAWING, + p1=(x0 - feedline_spacing / 2 - w, y0 + D / 2 + feed_length - 2), + p2=(x0 + feedline_spacing / 2 + w, ymax_frame_inner)) + + # add all created shapes to cell now. Upstream v4 snaps every polygon vertex to + # the 10 nm grid here; that narrows thin TopMetal2 traces below the 2 um minimum + # (TM2.a) at the miter joins, so ORCA leaves snapping to the DRCChecker stage + # (5 nm grid), as for every geometry. for geometry in all_geometries_list: cell.add(geometry) diff --git a/src/orca/geometry/cells/transformer.py b/src/orca/geometry/cells/transformer.py index 221d73f..2a4b9cf 100644 --- a/src/orca/geometry/cells/transformer.py +++ b/src/orca/geometry/cells/transformer.py @@ -20,10 +20,51 @@ def _ensure_active_pdk() -> None: except ValueError: gf.gpdk.PDK.activate() +def _outer_flat(diameter: float, width: float) -> float: + """Distance from a winding's centre to the outer edge of its flat sides, on the grid. + + The windings are octagons with flat sides facing +-x and +-y, so this is their + outermost metal in both directions (the feeds and center taps aside). + """ + return round((diameter / 2.0 * math.cos(math.radians(22.5)) + width / 2.0) / _GRID) * _GRID + + +def _ground_ring( + top_winding_diameter: float, + bottom_winding_diameter: float, + top_linewidth: float, + bottom_linewidth: float, + center_displacement: float, + gnd_upper_spacing: float, + gnd_lower_spacing: float, + gnd_side_spacing: float, + gnd_ring_width: float, +) -> tuple[float, float, float, float]: + """Where the windings and the ground ring sit. + + Returns: + tuple: The windings' offset from the origin (half the center displacement, on the + grid); the x of the ring's left and right outer edges, where the ports sit; and the + y of its top outer edge (the bottom one is at -y). The ring's inner edges keep + gnd_*_spacing from the windings' outermost metal. + """ + # On the grid like the windings (at most 2.5 nm from the requested offset): an off-grid + # shift would take every vertex off the grid, and snapping them afterwards could tilt + # the 45 degree sides. + half_offset = round(center_displacement / 2.0 / _GRID) * _GRID + top = _outer_flat(top_winding_diameter, top_linewidth) + bottom = _outer_flat(bottom_winding_diameter, bottom_linewidth) + port_xr = max(half_offset + top, -half_offset + bottom) + gnd_upper_spacing + gnd_ring_width + port_xl = min(half_offset - top, -half_offset - bottom) - gnd_lower_spacing - gnd_ring_width + tf_y = max(top, bottom) + gnd_side_spacing + gnd_ring_width + on_grid = lambda x: round(x / _GRID) * _GRID # noqa: E731 - one-line helper + return half_offset, on_grid(port_xl), on_grid(port_xr), on_grid(tf_y) + + def check_tf_octa_c_parameters( bottom_winding_diameter: float = 50.0, top_winding_diameter: float = 50.0, - center_displacement: float = 15.0, + center_displacement: float = 15.0, # noqa: ARG001 - same arguments as tf_octa_c bottom_linewidth: float = 5.0, bottom_center_tap_width: float = 0.0, lower_feed_type: int = 1, @@ -88,36 +129,20 @@ def check_tf_octa_c_parameters( f"a {diameter:g} octagon." ) - # The ports sit on the ring's inner edge, gnd_*_spacing - gnd_ring_width beyond the - # windings' vertices; closer than half a trace width, a feed would run back into its - # own winding instead of out of it. - for label, spacing in (("gnd_upper_spacing", gnd_upper_spacing), ("gnd_lower_spacing", gnd_lower_spacing)): - if spacing - gnd_ring_width < max(top_linewidth, bottom_linewidth) / 2.0: + # The ground ring keeps gnd_*_spacing from the windings' outermost metal; without a + # positive clearance it would run under the windings. + for label, spacing in ( + ("gnd_upper_spacing", gnd_upper_spacing), + ("gnd_lower_spacing", gnd_lower_spacing), + ("gnd_side_spacing", gnd_side_spacing), + ): + if spacing <= 0: raise ValueError( - f"{label} - gnd_ring_width = {spacing - gnd_ring_width:g} puts the ports inside " - "the windings; it must be at least half the widest trace." + f"{label} = {spacing:g} puts the ground ring under the windings; it is the " + "clearance between them and must be positive." ) - - # Ground ring: the ports on both sides and the ring bars must leave a positive opening. - tf_y = max(top_winding_diameter, bottom_winding_diameter) / 2.0 + gnd_side_spacing - port_xr = ( - max( - center_displacement / 2.0 + top_winding_diameter / 2.0, - -center_displacement / 2.0 + bottom_winding_diameter / 2.0, - ) - + gnd_upper_spacing - ) - port_xl = ( - min( - center_displacement / 2.0 - top_winding_diameter / 2.0, - -center_displacement / 2.0 - bottom_winding_diameter / 2.0, - ) - - gnd_lower_spacing - ) - if not (port_xr - gnd_ring_width > port_xl + gnd_ring_width and tf_y - gnd_ring_width > 0): - raise ValueError( - "Ground ring dimensions are invalid due to port spacing. Adjust parameters." - ) + if gnd_ring_width <= 0: + raise ValueError(f"gnd_ring_width = {gnd_ring_width:g} must be positive.") def tf_octa_c( @@ -156,10 +181,11 @@ def tf_octa_c( upper_feed_type: Center tap of the upper winding, same encoding as lower_feed_type (port ``oci`` on layer 205 when set to 1). feedline_spacing: Feedline spacing (gap between inner sides of feed lines). - gnd_upper_spacing: Ring spacing on the upper winding side. - gnd_lower_spacing: Ring spacing on the lower winding side. - gnd_side_spacing: Ring spacing at the side. - gnd_ring_width: Ring width. + gnd_upper_spacing: Clearance from the windings' outermost metal to the ground + ring's inner edge on the right, the side of the upper winding's feeds. + gnd_lower_spacing: The same on the left, the side of the lower winding's feeds. + gnd_side_spacing: The same at the top and bottom. + gnd_ring_width: Ring width. The ports sit on the ring's outer edge. textlabel: Text placed at the transformer's centre on the TEXT layer, to read the dimensions in a layout viewer. Empty lists every dimension drawn. """ @@ -207,26 +233,18 @@ def tf_octa_c( # Top Center Tap goes LEFT. It crosses Bot Gap. fs_bot = max(feedline_spacing, top_centertap_width) if draw_top_tap else feedline_spacing - # Geometry Limits - tf_y = max(top_winding_diameter, bottom_winding_diameter) / 2.0 + gnd_side_spacing - - # The windings are built on the manufacturing grid, so their centres are placed on it - # too (at most 2.5 nm from the requested offset): an off-grid shift would take every - # vertex off the grid, and snapping them afterwards could tilt the 45 degree sides. - half_offset = round(center_displacement / 2.0 / _GRID) * _GRID - - def on_grid(x: float) -> float: - return round(x / _GRID) * _GRID - - # X Limits for Ports, on the grid like the windings whose feeds end there - # Note: Winding edges are approx at center +/- diameter/2 - top_right_x = half_offset + (top_winding_diameter / 2.0) - bot_right_x = -half_offset + (bottom_winding_diameter / 2.0) - port_xr = on_grid(max(top_right_x, bot_right_x) + gnd_upper_spacing) - - top_left_x = half_offset - (top_winding_diameter / 2.0) - bot_left_x = -half_offset - (bottom_winding_diameter / 2.0) - port_xl = on_grid(min(top_left_x, bot_left_x) - gnd_lower_spacing) + # Windings and ground ring, the ring gnd_*_spacing clear of the windings' metal + half_offset, port_xl, port_xr, tf_y = _ground_ring( + top_winding_diameter, + bottom_winding_diameter, + top_linewidth, + bottom_linewidth, + center_displacement, + gnd_upper_spacing, + gnd_lower_spacing, + gnd_side_spacing, + gnd_ring_width, + ) # ------------------------------------------------- # 2. Helper: Winding Generator @@ -288,7 +306,7 @@ def region(shape): # which the SG13G2 angle rule rejects. Rounding the outer octagon's vertex offset # up and the inner one's down keeps the diagonal trace at least `width` wide. tan_22 = math.tan(math.radians(22.5)) - outer = round((diameter / 2.0 * math.cos(math.radians(22.5)) + width / 2.0) / _GRID) * _GRID + outer = _outer_flat(diameter, width) inner = outer - width apothem = outer - width / 2.0 # flat side of the trace's centre line winding = region(octagon(outer, math.ceil(outer * tan_22 / _GRID) * _GRID)) - region( @@ -333,8 +351,8 @@ def region(shape): center_x=half_offset, center_y=0, rotation_deg=0, - feed_target_x=port_xr - gnd_ring_width, - centertap_target_x=port_xl + gnd_ring_width if draw_top_tap else None, + feed_target_x=port_xr, + centertap_target_x=port_xl if draw_top_tap else None, centertap_width=top_centertap_width, ) @@ -347,8 +365,8 @@ def region(shape): center_x=-half_offset, center_y=0, rotation_deg=180, - feed_target_x=port_xl + gnd_ring_width, - centertap_target_x=port_xr - gnd_ring_width if draw_bottom_tap else None, + feed_target_x=port_xl, + centertap_target_x=port_xr if draw_bottom_tap else None, centertap_width=bottom_centertap_width, ) @@ -364,6 +382,9 @@ def region(shape): y_bot_n = -fs_bot / 2.0 - bottom_linewidth / 2.0 # Zero-width paths on port layers create Palace's 2D vertical port sheets. + # The feeds and center taps run across the ground ring to its outer edge, where the + # ports sit: the reference planes are the cell's boundary, so cells can be abutted + # with touching ports, as the inductor's. def add_port_marker(center, width, layer, orientation): dbu = c.layout().dbu @@ -383,13 +404,13 @@ def add_port_marker(center, width, layer, orientation): # OP: top winding, right side, upper port c.add_port( name="op", - center=(round(port_xr - gnd_ring_width, 2), round(y_top_p, 2)), + center=(port_xr, y_top_p), width=top_linewidth, orientation=0, layer=(201, 0), ) add_port_marker( - (round(port_xr - gnd_ring_width, 2), round(y_top_p, 2)), + (port_xr, y_top_p), top_linewidth, (201, 0), 0, @@ -397,13 +418,13 @@ def add_port_marker(center, width, layer, orientation): # ON: top winding, right side, lower port c.add_port( name="on", - center=(round(port_xr - gnd_ring_width, 2), round(y_top_n, 2)), + center=(port_xr, y_top_n), width=top_linewidth, orientation=0, layer=(202, 0), ) add_port_marker( - (round(port_xr - gnd_ring_width, 2), round(y_top_n, 2)), + (port_xr, y_top_n), top_linewidth, (202, 0), 0, @@ -412,26 +433,26 @@ def add_port_marker(center, width, layer, orientation): if draw_top_tap: c.add_port( name="oci", - center=(round(port_xl + gnd_ring_width, 2), 0.0), + center=(port_xl, 0.0), width=top_centertap_width, orientation=180, layer=(205, 0), ) add_port_marker( - (round(port_xl + gnd_ring_width, 2), 0.0), top_centertap_width, (205, 0), 180 + (port_xl, 0.0), top_centertap_width, (205, 0), 180 ) ### BOT LAYER (ports on the LEFT) -> Port 3 and 4 -> Layer 203, 204 # IP: bottom winding, left side, upper port c.add_port( name="ip", - center=(round(port_xl + gnd_ring_width, 2), round(y_bot_p, 2)), + center=(port_xl, y_bot_p), width=bottom_linewidth, orientation=180, layer=(203, 0), ) add_port_marker( - (round(port_xl + gnd_ring_width, 2), round(y_bot_p, 2)), + (port_xl, y_bot_p), bottom_linewidth, (203, 0), 180, @@ -439,13 +460,13 @@ def add_port_marker(center, width, layer, orientation): # IN: bottom winding, left side, lower port c.add_port( name="in", - center=(round(port_xl + gnd_ring_width, 2), round(y_bot_n, 2)), + center=(port_xl, y_bot_n), width=bottom_linewidth, orientation=180, layer=(204, 0), ) add_port_marker( - (round(port_xl + gnd_ring_width, 2), round(y_bot_n, 2)), + (port_xl, y_bot_n), bottom_linewidth, (204, 0), 180, @@ -454,13 +475,13 @@ def add_port_marker(center, width, layer, orientation): if draw_bottom_tap: c.add_port( name="ico", - center=(round(port_xr - gnd_ring_width, 2), 0.0), + center=(port_xr, 0.0), width=bottom_centertap_width, orientation=0, layer=(206, 0), ) add_port_marker( - (round(port_xr - gnd_ring_width, 2), 0.0), bottom_centertap_width, (206, 0), 0 + (port_xr, 0.0), bottom_centertap_width, (206, 0), 0 ) # ------------------------------------------------- @@ -479,12 +500,12 @@ def add_port_marker(center, width, layer, orientation): # Top bar top = gf.components.rectangle(size=(outer_w, gnd_ring_width), layer=LAYER_RING) top_ref = c << top - top_ref.move((round(port_xl, 2), round(tf_y - gnd_ring_width, 2))) + top_ref.move((port_xl, tf_y - gnd_ring_width)) # Bottom bar bot = gf.components.rectangle(size=(outer_w, gnd_ring_width), layer=LAYER_RING) bot_ref = c << bot - bot_ref.move((round(port_xl, 2), round(-tf_y, 2))) + bot_ref.move((port_xl, -tf_y)) # Left bar left_h = outer_h - 2 * gnd_ring_width @@ -493,7 +514,7 @@ def add_port_marker(center, width, layer, orientation): size=(gnd_ring_width, left_h), layer=LAYER_RING ) left_ref = c << left - left_ref.move((round(port_xl, 2), round(-tf_y + gnd_ring_width, 2))) + left_ref.move((port_xl, -tf_y + gnd_ring_width)) # Right bar if left_h > 0: @@ -502,7 +523,7 @@ def add_port_marker(center, width, layer, orientation): ) right_ref = c << right right_ref.move( - (round(port_xr - gnd_ring_width, 2), round(-tf_y + gnd_ring_width, 2)) + (port_xr - gnd_ring_width, -tf_y + gnd_ring_width) ) else: raise ValueError( diff --git a/src/orca/geometry/drc.py b/src/orca/geometry/drc.py index 9a1f92f..b66b957 100644 --- a/src/orca/geometry/drc.py +++ b/src/orca/geometry/drc.py @@ -8,16 +8,20 @@ and the rule names follow its report (``TM2.a``, ``TV2.c``, ...), so a finding here can be looked up in the IHP design rule manual directly. -Two entry points: :func:`snap_to_grid` moves every vertex of a layout onto -the grid, and :func:`check_layout` counts the violations per rule. Only the -*drawing* datatype of each layer is checked, as in the PDK deck; pin, text and -the gds2palace port layers are left alone. +Three entry points: :func:`snap_to_grid` moves every vertex of a layout onto +the grid, :func:`check_layout` counts the violations per rule, and +:func:`check_ports` checks that every gds2palace port marker touches the metals +its port connects (:func:`port_contacts` looks those up in the stackup). Only the +*drawing* datatype of each layer is checked, as in the PDK deck; pin and text +layers are left alone. """ from __future__ import annotations +import itertools +import xml.etree.ElementTree as ET from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Final +from typing import TYPE_CHECKING, Any, Final import klayout.db as kdb @@ -110,14 +114,34 @@ class DRCResult: """Vertices moved onto the grid before checking (0 when snapping was off).""" violations: dict[str, int] = field(default_factory=dict) """Violation count per rule name; rules without findings are absent.""" + port_findings: dict[str, int] = field(default_factory=dict) + """Port markers that do not touch their metal, per finding (see :func:`check_ports`).""" @property def clean(self) -> bool: - return not self.violations + return not self.violations and not self.port_findings + + @property + def ports_connected(self) -> bool: + return not self.port_findings @property def total(self) -> int: - return sum(self.violations.values()) + return sum(self.violations.values()) + sum(self.port_findings.values()) + + +@dataclass(frozen=True) +class PortContact: + """The metals one gds2palace port marker has to touch to excite the layout.""" + + number: int + """Port number, as in the simconfig and the Touchstone file.""" + marker_layer: int + """GDS layer number of the port marker (``source_layernum``).""" + metals: tuple[tuple[str, int], ...] + """``(layer name, GDS layer number)`` of every metal the port connects.""" + datatypes: tuple[int, ...] = (0,) + """Datatypes gds2palace reads the marker and the metals from (the simconfig's ``purpose``).""" def _snap(value: int, grid: int) -> int: @@ -270,7 +294,138 @@ def geometry_checks(layer: Layer, diagonal: bool) -> None: return violations -def check_gds_file(path: str, grid_nm: int = GRID_NM, snap: bool = True) -> DRCResult: +def port_contacts( + ports: list[dict[str, Any]], stackup_xml: str, datatypes: tuple[int, ...] = (0,) +) -> tuple[PortContact, ...]: + """The metals each port has to touch, with GDS numbers from the stackup. + + A port between two layers (``from_layername``/``to_layername``, a vertical sheet + in gds2palace) has to touch both; an in-plane port (``target_layername``) the one. + + Args: + ports: Port entries in the format of a simconfig's ``ports`` list, e.g. one + sample's from ``BaseGeometry.ports_for``. + stackup_xml: gds2palace stackup XML that maps the layer names to GDS numbers. + datatypes: Datatypes gds2palace reads the layout from (the simconfig's + ``purpose``). + + Returns: + tuple[PortContact, ...]: One entry per port, in the given order. + + Raises: + ValueError: A port names a layer the stackup does not define. + """ + # The stackup is a trusted file shipped with the geometry, not external input + root = ET.parse(stackup_xml).getroot() # noqa: S314 + layer_numbers = { + element.get("Name"): int(element.get("Layer", "")) + for element in root.iter("Layer") + if element.get("Name") and element.get("Layer", "").isdigit() + } + contacts = [] + for port in ports: + target = port.get("target_layername") + names = [target] if target else [port.get("from_layername"), port.get("to_layername")] + metals = [] + for name in filter(None, names): + if name not in layer_numbers: + raise ValueError( + f"Port {port['portnumber']} connects layer {name!r}, which the stackup " + f"{stackup_xml} does not define." + ) + metals.append((name, layer_numbers[name])) + contacts.append( + PortContact( + number=port["portnumber"], + marker_layer=port["source_layernum"], + metals=tuple(metals), + datatypes=tuple(datatypes), + ) + ) + return tuple(contacts) + + +def check_ports(layout: kdb.Layout, ports: tuple[PortContact, ...]) -> dict[str, int]: + """Count the port markers of *layout* that do not touch the metals their port connects. + + gds2palace places each port where its marker is drawn and does not check that + the port meets the conductors. A marker a few nanometres off its metal still + meshes and simulates, but as an open circuit: the S-parameters then describe + another circuit than the parameters say, with no error anywhere. Touching an + edge is enough, as for a vertical port sheet on the end face of a feed line. + + - ``port.missing``: no marker on the port's layer; + - ``port.open_``: markers that neither overlap nor touch ````. + + Args: + layout: Layout to check. It is not modified. + ports: The ports to check, from :func:`port_contacts`. + + Returns: + dict[str, int]: Finding counts keyed by name; empty when every port is connected. + """ + findings: dict[str, int] = {} + metals: dict[tuple[int, tuple[int, ...]], kdb.Region] = {} + for port in ports: + markers = _markers(layout, port.marker_layer, port.datatypes) + if not markers: + findings[f"port{port.number}.missing"] = 1 + continue + for name, layer_number in port.metals: + key = (layer_number, port.datatypes) + if key not in metals: + metals[key] = _region(layout, layer_number, port.datatypes) + open_markers = sum(marker.interacting(metals[key]).is_empty() for marker in markers) + if open_markers: + findings[f"port{port.number}.open_{name}"] = open_markers + return findings + + +def _region(layout: kdb.Layout, layer_number: int, datatypes: tuple[int, ...]) -> kdb.Region: + """Everything drawn on a layer in any of *datatypes*, flattened and merged.""" + region = kdb.Region() + for datatype in datatypes: + index = layout.find_layer(layer_number, datatype) + if index is not None: + for top in layout.top_cells(): + region += kdb.Region(top.begin_shapes_rec(index)) + return region.merged() + + +def _markers( + layout: kdb.Layout, layer_number: int, datatypes: tuple[int, ...] +) -> list[kdb.Region | kdb.Edges]: + """Each port marker on a layer, in top-cell coordinates. + + A zero-width path, gds2palace's vertical port sheet, has no area, so it is + returned as its centre line; any other shape as its polygon. + """ + markers: list[kdb.Region | kdb.Edges] = [] + for datatype in datatypes: + index = layout.find_layer(layer_number, datatype) + if index is None: + continue + for top in layout.top_cells(): + iterator = top.begin_shapes_rec(index) + while not iterator.at_end(): + shape, trans = iterator.shape(), iterator.trans() + if shape.is_path() and shape.path.width == 0: + points = [trans * p for p in shape.path.each_point()] + markers.append( + kdb.Edges([kdb.Edge(a, b) for a, b in itertools.pairwise(points)]) + ) + elif not shape.is_text(): + markers.append(kdb.Region(shape.polygon.transformed(trans))) + iterator.next() + return markers + + +def check_gds_file( + path: str, + grid_nm: int = GRID_NM, + snap: bool = True, + ports: tuple[PortContact, ...] = (), +) -> DRCResult: """Check one GDS file, optionally snapping it to the grid first. Args: @@ -278,6 +433,8 @@ def check_gds_file(path: str, grid_nm: int = GRID_NM, snap: bool = True) -> DRCR grid_nm: Manufacturing grid in nanometres. snap: Snap all vertices to the grid and write the file back before checking. Off-grid vertices are then repaired rather than reported. + ports: Ports whose markers must touch their metals (:func:`check_ports`), + checked after snapping; empty skips the port check. Returns: DRCResult: Vertices moved and the violations found. @@ -289,7 +446,11 @@ def check_gds_file(path: str, grid_nm: int = GRID_NM, snap: bool = True) -> DRCR moved = snap_to_grid(layout, grid_nm) if moved: layout.write(path) - return DRCResult(snapped_vertices=moved, violations=check_layout(layout, grid_nm)) + return DRCResult( + snapped_vertices=moved, + violations=check_layout(layout, grid_nm), + port_findings=check_ports(layout, ports), + ) def _grid_dbu(layout: kdb.Layout, grid_nm: int) -> int: diff --git a/src/orca/geometry/presets/inductor/inductor_octa.py b/src/orca/geometry/presets/inductor/inductor_octa.py index b1dbfdb..107e241 100644 --- a/src/orca/geometry/presets/inductor/inductor_octa.py +++ b/src/orca/geometry/presets/inductor/inductor_octa.py @@ -14,18 +14,26 @@ from orca.geometry.presets.paths import StackupXML if TYPE_CHECKING: + import numpy as np + import skrf as rf + from orca.training.datasets.base_dataset import BaseDataset -# All ports run up from a ground on Metal5 (matches "from_layername": "Metal5" in -# the simcfg): ports 1/2 to the feeds on TopMetal1, port 3 to the center tap on -# TopMetal2. The ground is a closed square ring around the spiral for every N, so -# one connected reference serves all three ports whether the center tap leaves -# between the feeds (even N) or at the top (odd N). The feeds and the center tap run out to the ring's outer edge, -# where the ports sit. The feed on TopMetal1 requires N >= 2 turns (see -# symmetric_octa_IHP: N == 1 feeds on TopMetal2 instead). +# The spiral is that of the upstream gds2palace example (synthesize_ihp_inductor_v4). +# Its EM ground is a closed square Metal5 ring at a fixed distance around it, rather +# than upstream's Metal1 frame D/2 away, which makes large inductors very large. The +# feeds run out to the ring's outer edge, where the ports sit, so each port's +# reference plane is the boundary of the cell and simulated cells can be placed side +# by side, ports touching. With the 20 µm gap and the 10 µm ring the feeds reach 30 µm +# beyond the outer diameter, the lead length of the upstream example. All ports +# run up from the Metal5 ground (matches "from_layername": "Metal5" in the simcfg): +# ports 1/2 to the feeds, port 3 to the center tap on TopMetal2. The feeds are on +# TopMetal1, as the simcfg says, except for a single turn, which feeds on TopMetal2; +# simulation_ports moves ports 1/2 there for N == 1, as the upstream script does. GROUND_LAYER = 67 # Metal5 -GROUND_SPACING = 20.0 # µm, gap between inductor outer edge and the ground -GROUND_DEPTH = 20.0 # µm, width of the ground ring bars +GROUND_SPACING = 20.0 # µm, gap between the inductor's outer diameter and the ring +GROUND_DEPTH = 10.0 # µm, width of the ground ring bars, as the transformer's +SINGLE_TURN_FEED_LAYER = "TopMetal2" # stackup name of the metal a single turn feeds on # Built per instance rather than shared as a class attribute - see the note in @@ -59,15 +67,18 @@ class InductorOcta(BaseGeometry): Represents a symmetric octagonal spiral inductor geometry (IHP SG13G2). 3-port spiral inductor (LA, LB and the center tap LC), ported from the - gds2palace IHP example by Volker Muehlhaus. Requires N >= 2 turns, since the feedline sits on - TopMetal1 (single-turn inductors feed on TopMetal2 instead). + gds2palace IHP example by Volker Muehlhaus. Multi-turn spirals feed on TopMetal1, + a single turn on TopMetal2; the ports follow the feeds (:meth:`simulation_ports`). """ name: str = "inductor_octa" # Conformal SiO2/passivation over TopMetal2 (gds2palace L6n2 study): the planar # stackup fills the gaps between turns with oxide and overstates the turn-to-turn - # capacitance. Paired with refined_cellsize = 5 in the simcfg, the study's fast - # "daily driver" setting; it meshes smaller than planar at 2 µm. + # capacitance. The simcfg meshes it with refined_cellsize = 2 rather than the study's + # fast 5: with the ports flush with the ring's outer edge, the ring face, the 0.85 µm + # port sheet and the feed's end face lie in one plane, and at 5 µm gmsh filled it with + # flat tetrahedra in about 40% of the 4-5 turn layouts, which Palace cannot solve. At + # 2 µm none of 70 such layouts had one, at about twice the elements. stackup_xml: str = StackupXML.SG13G2_FEM_200um_passi3D simconfig_filename: str = os.path.join(os.path.dirname(__file__), "inductor_octa.simcfg") input_parameter_iterator: InputParameterIterator = field( @@ -91,6 +102,22 @@ def create_dataset(self) -> "BaseDataset": output_normalizer=StandardNormalizer(), ) + def simulation_ports(self, params: dict[str, Any]) -> list[dict[str, Any]]: + ports = super().simulation_ports(params) + if round(params["turns"]) == 1: + # The feeds (LA, LB) of a single turn are on TopMetal2, not TopMetal1 + for port in ports: + if port["portnumber"] in (1, 2): + port["to_layername"] = SINGLE_TURN_FEED_LAYER + return ports + + def electrical_parameters(self, ntwk: "rf.Network") -> dict[str, "np.ndarray"]: + from orca.utils.postprocessing import inductor_parameters + + # Ports 1 and 2 feed the two ends (LA, LB); port 3, the center tap, is AC-grounded + # as in differential use + return inductor_parameters(ntwk, ends=(0, 1), shorted=(2,)) + def feasibility_constraints(self) -> list[str]: # get_min_outer_diameter as one expression; kept in step with it by a test. two_vias = 2 * VIA_SIZE + VIA_GAP + 2 * VIA_MARGIN @@ -137,9 +164,9 @@ def create_gds_file(name: str, output_path: str, params: dict[str, Any]) -> str: LBE=False, forEM=True, ground_layer=GROUND_LAYER, - ground_style="ring", ring_spacing=GROUND_SPACING, ring_width=GROUND_DEPTH, + feeds_to_ring_edge=True, filename=output_path, ) return output_path diff --git a/src/orca/geometry/presets/inductor/inductor_octa.simcfg b/src/orca/geometry/presets/inductor/inductor_octa.simcfg index ba8cbe8..e0be3fa 100644 --- a/src/orca/geometry/presets/inductor/inductor_octa.simcfg +++ b/src/orca/geometry/presets/inductor/inductor_octa.simcfg @@ -10,7 +10,7 @@ "fstart": 1.0, "fstop": 500.0, "fstep": 1.0, - "refined_cellsize": 5.0, + "refined_cellsize": 2.0, "order": 2, "cells_per_wavelength": 10.0, "meshsize_max": 100.0, diff --git a/src/orca/geometry/presets/transformer/tf_octa_c_ports.py b/src/orca/geometry/presets/transformer/tf_octa_c_ports.py index 140f1d8..904f183 100644 --- a/src/orca/geometry/presets/transformer/tf_octa_c_ports.py +++ b/src/orca/geometry/presets/transformer/tf_octa_c_ports.py @@ -1,6 +1,6 @@ import os from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING, Any, ClassVar from orca import BaseGeometry from orca.geometry.cells.transformer import check_tf_octa_c_parameters, tf_octa_c @@ -9,6 +9,9 @@ from orca.geometry.presets.paths import StackupXML if TYPE_CHECKING: + import numpy as np + import skrf as rf + from orca.training.datasets.base_dataset import BaseDataset # Layout choices fixed for every sample (µm). They are part of the device the model @@ -19,8 +22,11 @@ #: Gap between the two feed lines of a winding: the crossing 3 µm tap plus 1 µm on #: either side. Close feeds are a tight differential pair with a small loop inductance. FEED_GAP = 5.0 -#: The ports sit on the ground ring's inner edge, GROUND_SPACING - GROUND_RING_WIDTH = -#: 10 µm beyond the windings' vertices, so short feeds are part of every model. +#: The Metal5 ground ring keeps GROUND_SPACING = 20 µm clear of the windings' outermost +#: metal and is GROUND_RING_WIDTH = 10 µm wide, as the inductor's. The feeds run across it +#: to its outer edge, 30 µm beyond the windings, where the ports sit: the reference planes +#: are the boundary of the cell, so simulated cells (inductors too) can be abutted with +#: touching ports. The feeds and their crossing of the ring are part of every model. GROUND_SPACING = 20.0 GROUND_RING_WIDTH = 10.0 @@ -64,16 +70,32 @@ class TransformerOcta(BaseGeometry): One single-turn winding on TopMetal2 (ports ``op``/``on`` on the right, center tap ``oci`` to the left) over one on TopMetal1 (ports ``ip``/``in`` on the left, center - tap ``ico`` to the right), inside a Metal5 ground ring the ports refer to. + tap ``ico`` to the right), inside a Metal5 ground ring the ports refer to. The ports + sit on the ring's outer edge. """ name: str = "tf_octa_c_ports" - stackup_xml: str = StackupXML.SG13G2_FEM_200um + # Conformal SiO2/passivation over TopMetal2, as for the inductor (gds2palace L6n2 + # study). The simcfg meshes it with refined_cellsize = 2: the ports are flush with + # the ring's outer edge, and at coarser sizes gmsh can fill that plane with flat + # tetrahedra that Palace cannot solve (GDSConverter drops such meshes). + stackup_xml: str = StackupXML.SG13G2_FEM_200um_passi3D simconfig_filename: str = os.path.join(os.path.dirname(__file__), "tf_octa_c_ports.simcfg") input_parameter_iterator: InputParameterIterator = field( default_factory=_input_parameters ) + # Close to zero for weakly coupled windings, where a relative error says nothing + absolute_error_parameters: ClassVar[frozenset[str]] = frozenset({"k"}) + + def electrical_parameters(self, ntwk: "rf.Network") -> dict[str, "np.ndarray"]: + from orca.utils.postprocessing import transformer_parameters + + # Ports (simcfg order): 1 op, 2 on (top winding), 3 ip, 4 in (bottom winding), + # 5 oci and 6 ico (the center taps), which are AC-grounded as in differential use. + # The top winding (op/on) is the primary, as in COBRA's Lp/Qp goals. + return transformer_parameters(ntwk, primary=(0, 1), secondary=(2, 3), shorted=(4, 5)) + def create_dataset(self) -> "BaseDataset": # Imported here so the geometry can be drawn and simulated without the # "train" extra (PyTorch) installed. diff --git a/src/orca/pipeline/context.py b/src/orca/pipeline/context.py index b7f0720..5af6c0d 100644 --- a/src/orca/pipeline/context.py +++ b/src/orca/pipeline/context.py @@ -159,6 +159,11 @@ def palace_csv_path(self) -> str: """Where the GDS conversion stage writes its parameter table.""" return os.path.join(self.palace_sim_dir, f"{self.geometry.name}.csv") + @property + def mesh_report_path(self) -> str: + """Where the GDS conversion stage writes the worst element quality of each mesh.""" + return os.path.join(self.palace_sim_dir, f"{self.geometry.name}_mesh_report.csv") + @property def result_dir(self) -> str: """Directory holding the Touchstone results the model is trained on.""" @@ -224,6 +229,7 @@ def to_json_dict(self) -> dict[str, Any]: "gds_coverage_plot": self.gds_coverage_plot_path, "drc_report": self.drc_report_path, "palace_sim_dir": self.palace_sim_dir, + "mesh_report": self.mesh_report_path, "result_dir": self.result_dir, "result_csv": self.result_csv, "model_dir": self.model_dir, diff --git a/src/orca/pipeline/drc_stage.py b/src/orca/pipeline/drc_stage.py index 5483102..4412f20 100644 --- a/src/orca/pipeline/drc_stage.py +++ b/src/orca/pipeline/drc_stage.py @@ -1,3 +1,4 @@ +import json import os from collections import Counter from collections.abc import Callable @@ -7,9 +8,10 @@ import pandas as pd import tqdm -from orca.geometry.drc import GRID_NM, DRCResult, check_gds_file +from orca.geometry.drc import GRID_NM, DRCResult, PortContact, check_gds_file, port_contacts from orca.logger import logger from orca.pipeline.pipeline_stage import PipelineStage +from orca.simulation.simulate import read_simconfig if TYPE_CHECKING: from orca.pipeline.context import PipelineContext @@ -27,6 +29,11 @@ class DRCChecker(PipelineStage): cannot be built. This stage repairs what can be repaired in place (the grid) and records the rest, so the conversion stage only picks up clean layouts. + It also checks that every port marker touches the metals its port connects, with + each layout's own ports (``BaseGeometry.ports_for``). A marker that misses its metal by a few nanometres + still simulates, as an open circuit, so the result describes another circuit + than the parameters say. Such layouts are always left out. + Outputs, in the geometry folder: ``_drc_report.csv`` with the violation counts of every layout, and ``_drc.csv``, the parameter table of the layouts that passed, in the same layout as the GDS table. @@ -34,7 +41,11 @@ class DRCChecker(PipelineStage): """ def __init__( - self, grid_nm: int = GRID_NM, snap_to_grid: bool = True, drop_violations: bool = True + self, + grid_nm: int = GRID_NM, + snap_to_grid: bool = True, + drop_violations: bool = True, + check_ports: bool = True, ): """ Args: @@ -44,6 +55,10 @@ def __init__( drop_violations (bool): Leave layouts with violations out of the parameter table the later stages use. Off, every layout goes on and the violations are only reported. + check_ports (bool): Check that every port marker touches the metals its port + connects, with the ports of each sample. Layouts with an open port are left out + whatever ``drop_violations`` says: their simulation would describe a + circuit with that port disconnected. """ super().__init__(name="DRC Checker", index=1) if grid_nm <= 0: @@ -51,6 +66,7 @@ def __init__( self.grid_nm = grid_nm self.snap_to_grid = snap_to_grid self.drop_violations = drop_violations + self.check_ports = check_ports def run( self, @@ -67,6 +83,7 @@ def run( gds_data = pd.read_csv(gds_csv) gds_dir = os.path.dirname(gds_csv) + ports = self._port_contacts(context, gds_data) if self.check_ports else {} logger.info( f"Starting DRC of {len(gds_data)} GDS files on a {self.grid_nm} nm grid " f"using {context.num_processes} CPU cores." @@ -76,7 +93,11 @@ def run( with ProcessPoolExecutor(max_workers=context.num_processes) as executor: futures = { executor.submit( - check_gds_file, os.path.join(gds_dir, name), self.grid_nm, self.snap_to_grid + check_gds_file, + os.path.join(gds_dir, name), + self.grid_nm, + self.snap_to_grid, + ports.get(name, ()), ): name for name in gds_data["name"] } @@ -105,7 +126,10 @@ def run( "snapped_vertices": results[name].snapped_vertices, "violations": results[name].total, "rules": ";".join( - f"{rule}:{count}" for rule, count in sorted(results[name].violations.items()) + f"{rule}:{count}" + for rule, count in sorted( + (results[name].violations | results[name].port_findings).items() + ) ), } for name in gds_data["name"] @@ -114,12 +138,15 @@ def run( report.to_csv(context.drc_report_path, index=False) clean = gds_data["name"].map(lambda name: results[name].clean) - passed = gds_data if not self.drop_violations else gds_data[clean] + connected = gds_data["name"].map(lambda name: results[name].ports_connected) + passed = gds_data[clean] if self.drop_violations else gds_data[connected] passed.to_csv(context.drc_csv_path, index=False) summary: Counter[str] = Counter() for result in results.values(): summary.update(result.violations) + summary.update(result.port_findings) + self._log_open_ports(results, context.drc_report_path) self._log_summary( len(gds_data), int(clean.sum()), @@ -132,6 +159,40 @@ def run( context.drc_summary = dict(sorted(summary.items())) return context + @staticmethod + def _port_contacts( + context: "PipelineContext", gds_data: pd.DataFrame + ) -> dict[str, tuple[PortContact, ...]]: + """ + Per layout, the metals each of its ports has to touch: the sample's ports + (``BaseGeometry.ports_for``), with layer numbers from the stackup. + """ + geometry = context.geometry + datatypes = tuple( + read_simconfig(geometry.simconfig_filename)["saved_values"].get("purpose", [0]) + ) + # Most samples share one port set; look each distinct set up once + by_ports: dict[str, tuple[PortContact, ...]] = {} + contacts = {} + for row in gds_data.to_dict("records"): + name = row.pop("name") + ports = geometry.ports_for(row) + key = json.dumps(ports, sort_keys=True) + if key not in by_ports: + by_ports[key] = port_contacts(ports, geometry.stackup_xml, datatypes) + contacts[name] = by_ports[key] + return contacts + + @staticmethod + def _log_open_ports(results: dict[str, DRCResult], report_path: str) -> None: + open_ports = [name for name, result in results.items() if not result.ports_connected] + if open_ports: + logger.error( + f"{len(open_ports)} of {len(results)} layouts have port markers that do not " + "touch their metal; they would simulate as open circuits and were dropped. " + f"Examples: {', '.join(open_ports[:5])}. Details are in {report_path}." + ) + def _log_summary( self, total: int, n_clean: int, summary: Counter[str], snapped: int, report_path: str ) -> None: diff --git a/src/orca/pipeline/gds_conversion_stage.py b/src/orca/pipeline/gds_conversion_stage.py index 64a2638..cc49631 100644 --- a/src/orca/pipeline/gds_conversion_stage.py +++ b/src/orca/pipeline/gds_conversion_stage.py @@ -12,6 +12,7 @@ from orca.pipeline.pipeline_stage import PipelineStage from orca.pipeline.resume import append_row, resume_table from orca.simulation.gds_converter import create_palace_model_from_gds +from orca.simulation.simulate import read_simconfig if TYPE_CHECKING: from orca.geometry.base_geometry import BaseGeometry @@ -34,19 +35,39 @@ class GDSConverter(PipelineStage): Palace models of an earlier, possibly aborted run are kept: a layout is only converted if the Palace table has no row for it with the same parameters, or its Palace config is gone. + + Every new mesh is checked for flat, inverted or corrupt elements, which Palace cannot + solve (``HasPositiveFiniteDiagonal(...) is false``): such a sample is left out of + the Palace table instead of failing on the cluster. The worst element quality of + each mesh is written to ``_mesh_report.csv`` next to the table. """ - def __init__(self, timeout: float = 60.0, overwrite: bool = False): + def __init__( + self, + timeout: float = 180.0, + overwrite: bool = False, + min_element_quality: float = 1e-6, + warn_element_quality: float = 1e-4, + ): """ Args: timeout (float): Maximum seconds a single GDS conversion may take - before its worker is killed and the sample is skipped. + before its worker is killed and the sample is skipped. It only guards + against gmsh hanging: a large inductor at refined_cellsize 2 already + takes up to ~50 s, and a sample cut off here is lost silently. overwrite (bool): Delete the Palace model folder first instead of keeping the models an earlier run left there. + min_element_quality (float): Meshes whose worst element quality (gmsh's + minSICN: 1 regular, 0 flat) is at or below this are left out. Flat + elements come out at 0 up to rounding, so the default only catches those. + warn_element_quality (float): Meshes with a worst element quality below this + are kept but reported. Meshes that Palace has solved had 6e-4 and more. """ super().__init__(name="GDS Converter", index=2) self.timeout = timeout self.overwrite = overwrite + self.min_element_quality = min_element_quality + self.warn_element_quality = warn_element_quality def run( self, @@ -131,12 +152,14 @@ def run( "stackup_xml": geometry.stackup_xml, "simconfig_filename": geometry.simconfig_filename, "show_mesh_results": False, + "ports": geometry.ports_for(params), }, timeout=self.timeout, ) futures[future] = name # Collect finished results + qualities: dict[str, float] = {} # Print progress bar using tqdm for i, future in enumerate( tqdm.tqdm( @@ -145,7 +168,11 @@ def run( ): name = futures[future] try: - geo_name, params, config_name, sim_path, data_dir = future.result() + geo_name, params, config_name, sim_path, data_dir, quality = future.result() + qualities[name] = quality + if quality <= self.min_element_quality: + # Palace would stop on this mesh; keep it out of the table + continue # Save input parameters to CSV append_row( palace_csv, @@ -175,5 +202,49 @@ def run( logger.info("GDS conversion completed.") + self._report_mesh_quality(qualities, context) context.palace_csv = palace_csv return context + + def _report_mesh_quality(self, qualities: dict[str, float], context: "PipelineContext") -> None: + """Write the worst element quality of each new mesh and warn about poor ones.""" + if not qualities: + return + report_path = context.mesh_report_path + report = pd.DataFrame( + {"name": list(qualities), "worst_element_quality": list(qualities.values())} + ) + if os.path.exists(report_path): + # Meshes of earlier runs stay listed; a sample meshed again replaces its row + earlier = pd.read_csv(report_path) + report = pd.concat([earlier[~earlier["name"].isin(qualities)], report]) + report.to_csv(report_path, index=False) + + flat = sorted(n for n, q in qualities.items() if q <= self.min_element_quality) + poor = sorted( + n for n, q in qualities.items() + if self.min_element_quality < q < self.warn_element_quality + ) + cell_size = read_simconfig(context.geometry.simconfig_filename)["saved_values"].get( + "refined_cellsize" + ) + hint = ( + f"If this happens often, refined_cellsize = {cell_size:g} µm in the simcfg may be too " + "coarse for these layouts." + if cell_size is not None + else "If this happens often, the mesh may be too coarse for these layouts." + ) + if flat: + logger.warning( + f"{len(flat)} of {len(qualities)} meshes contain flat, inverted or corrupt " + f"elements (worst quality at or below {self.min_element_quality:g}; -inf: gmsh " + f"cannot read the mesh back) and were left out, as Palace cannot solve them. " + f"{hint} Samples: {', '.join(flat[:10])}" + f"{' ...' if len(flat) > 10 else ''}. Qualities are in {report_path}." + ) + if poor: + logger.warning( + f"{len(poor)} of {len(qualities)} meshes have poorly shaped elements (worst " + f"quality below {self.warn_element_quality:g}); they were kept. {hint} " + f"Samples: {', '.join(poor[:10])}{' ...' if len(poor) > 10 else ''}." + ) diff --git a/src/orca/pipeline/test_model_stage.py b/src/orca/pipeline/test_model_stage.py index 5e85f51..8621ae9 100644 --- a/src/orca/pipeline/test_model_stage.py +++ b/src/orca/pipeline/test_model_stage.py @@ -19,9 +19,8 @@ TorchNetworkPredictor, ) from orca.utils.postprocessing import ( - calculate_electrical_parameters, + plot_electrical_parameters, plot_errors_vs_frequency, - plot_rfic_transformer_metrics, pointwise_relative_error, ) @@ -32,11 +31,6 @@ from orca.pipeline.context import PipelineContext -#: Electrical parameters whose error is an absolute difference instead of a relative one. -#: The coupling factor k is close to zero for weakly coupled layouts, where a relative error -#: explodes without saying anything about the model. -ABSOLUTE_ERROR_PARAMETERS = frozenset({"k"}) - #: Name of the S-parameter error among the per-frequency errors S_PARAMETER_ERROR = "|S|" @@ -55,12 +49,15 @@ class EvaluationErrors: per_frequency: Per quantity (``|S|`` and each electrical parameter that is a curve), the error of every geometry at every frequency point, shape ``(n_geometries, n_freq)``. Relative errors are in percent, except for - :data:`ABSOLUTE_ERROR_PARAMETERS` and ``|S|``. + ``|S|`` and the ``absolute_parameters``. + absolute_parameters: Electrical parameters whose error is absolute, as the + geometry declared them (``BaseGeometry.absolute_error_parameters``). """ per_geometry: pd.DataFrame frequencies: np.ndarray | None = None per_frequency: dict[str, np.ndarray] = field(default_factory=dict) + absolute_parameters: frozenset[str] = frozenset() class ModelTester(PipelineStage): @@ -118,7 +115,7 @@ def run( logger.error(str(e)) return context - errors = self.evaluate(test_dataset, predictor, progress_callback) + errors = self.evaluate(test_dataset, predictor, progress_callback, context.geometry) results = self.summarize(errors.per_geometry) if results: @@ -190,6 +187,7 @@ def test_model( test_dataset: GeoToNtwkDataset, predictor: NetworkPredictor, progress_callback: Callable[[str, int, int, str], None] | None = None, + geometry: "BaseGeometry | None" = None, ) -> dict[str, Any]: """ Evaluates the predictor on the test dataset. @@ -197,23 +195,33 @@ def test_model( Returns: dict: The summary of :meth:`summarize`; empty if nothing was evaluated. """ - return self.summarize(self.evaluate(test_dataset, predictor, progress_callback).per_geometry) + errors = self.evaluate(test_dataset, predictor, progress_callback, geometry) + return self.summarize(errors.per_geometry) def evaluate( self, test_dataset: GeoToNtwkDataset, predictor: NetworkPredictor, progress_callback: Callable[[str, int, int, str], None] | None = None, + geometry: "BaseGeometry | None" = None, ) -> EvaluationErrors: """ The errors of every test geometry, overall and at every frequency point. + Args: + test_dataset: The held-out networks and their parameters. + predictor: The model to test. + progress_callback: Called after every geometry. + geometry: Supplies the electrical parameters to compare + (``BaseGeometry.electrical_parameters``); None compares S-parameters only. + Returns: EvaluationErrors: ``per_geometry`` holds one row per geometry: ``name``, the geometry parameters, the mean and maximum absolute S-parameter error, the mean absolute error per frequency band (``s_error ``), and the median error of each electrical parameter, relative in percent (`` error %``) - or, for :data:`ABSOLUTE_ERROR_PARAMETERS`, absolute (`` abs error``). + or, for the geometry's ``absolute_error_parameters``, absolute + (`` abs error``). Empty in plot mode. """ num_samples = len(test_dataset) @@ -231,8 +239,13 @@ def evaluate( ntwk_gt.name = "Ground Truth" if self.plot: - plot_rfic_transformer_metrics(ntwk_gt) - plot_rfic_transformer_metrics(ntwk_pred) + if geometry is not None: + plot_electrical_parameters( + ntwk_gt.f, + geometry.electrical_parameters(ntwk_gt), + geometry.electrical_parameters(ntwk_pred), + title=test_dataset.names[i], + ) continue abs_error = np.abs(ntwk_pred.s - ntwk_gt.s) # (n_freq, n_ports, n_ports) @@ -245,8 +258,10 @@ def evaluate( } for label, in_band in self._frequency_bands(ntwk_gt.f).items(): row[f"s_error {label}"] = float(per_frequency[in_band].mean()) - electrical, electrical_curves = self._electrical_errors( - ntwk_pred, ntwk_gt, test_dataset.names[i] + electrical, electrical_curves = ( + self._electrical_errors(ntwk_pred, ntwk_gt, test_dataset.names[i], geometry) + if geometry is not None + else ({}, {}) ) row |= electrical rows.append(row) @@ -281,6 +296,9 @@ def evaluate( for quantity, stack in curves.items() if len(stack) == n_on_grid }, + absolute_parameters=( + geometry.absolute_error_parameters if geometry is not None else frozenset() + ), ) @staticmethod @@ -360,7 +378,7 @@ def frequency_profile(errors: EvaluationErrors) -> pd.DataFrame: with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) levels = np.nanpercentile(stack, list(PERCENTILES.values()), axis=0) - absolute = quantity == S_PARAMETER_ERROR or quantity in ABSOLUTE_ERROR_PARAMETERS + absolute = quantity == S_PARAMETER_ERROR or quantity in errors.absolute_parameters tables.append( pd.DataFrame( { @@ -404,10 +422,13 @@ def _frequency_bands(self, frequencies: np.ndarray) -> dict[str, np.ndarray]: @staticmethod def _electrical_errors( - predicted_ntwk: "rf.Network", reference_ntwk: "rf.Network", name: str + predicted_ntwk: "rf.Network", + reference_ntwk: "rf.Network", + name: str, + geometry: "BaseGeometry", ) -> tuple[dict[str, float], dict[str, np.ndarray]]: """ - The error of each electrical parameter that can be derived. + The error of each electrical parameter the geometry declares, where it can be derived. Returns: tuple: The median error per parameter, keyed by its per-geometry column name; @@ -415,8 +436,8 @@ def _electrical_errors( (not, say, the self-resonance frequency). """ try: - predicted = calculate_electrical_parameters(predicted_ntwk) - reference = calculate_electrical_parameters(reference_ntwk) + predicted = geometry.electrical_parameters(predicted_ntwk) + reference = geometry.electrical_parameters(reference_ntwk) except Exception as e: # noqa: BLE001 - skip samples whose metrics cannot be derived logger.debug(f"Could not compute electrical parameters for {name}: {e}") return {}, {} @@ -424,7 +445,7 @@ def _electrical_errors( medians: dict[str, float] = {} curves: dict[str, np.ndarray] = {} for param, gt in reference.items(): - if param in ABSOLUTE_ERROR_PARAMETERS: + if param in geometry.absolute_error_parameters: curve = np.abs(np.atleast_1d(predicted[param]) - np.atleast_1d(gt)).astype(float) curve[~np.isfinite(curve)] = np.nan column = f"{param} abs error" diff --git a/src/orca/pipeline/training_stage.py b/src/orca/pipeline/training_stage.py index 8bd664d..d84ebe7 100644 --- a/src/orca/pipeline/training_stage.py +++ b/src/orca/pipeline/training_stage.py @@ -12,7 +12,9 @@ from orca.pipeline.pipeline_stage import PipelineStage from orca.training.basis_expansion import BasisExpansion, get_basis_class from orca.training.datasets.base_dataset import BaseDataset +from orca.training.losses import SParameterLoss, srf_sample_weights from orca.training.models.base_model import OrcaModel, get_model_class +from orca.training.spec import FrequencyMode from orca.training.trainer import Trainer, TrainingConfig from orca.training.tuner import HyperparameterTuner @@ -45,6 +47,9 @@ def __init__( warmup_epochs: float = 1.0, grad_clip_norm: float | None = 1.0, allow_tf32: bool = False, + admittance_weight: float = 0.0, + passivity_weight: float = 0.0, + above_srf_weight: float = 1.0, ): """ Initializes the ModelTrainer stage with the architecture to train and optional @@ -93,6 +98,23 @@ def __init__( (Ampere and newer, e.g. A100): faster, with 10 instead of 23 mantissa bits in the multiplies. Tuning and the final training use it; testing and the exported model stay in full FP32. + admittance_weight: Weight of a loss term on the relative error of the + predicted admittance matrix Y, of Y as a whole and of its real part, which + track L and Q far more closely than S does (see + orca.training.losses.admittance_error). 0 (default) leaves it out; + 0.1 to 1 puts it on the scale of the S-parameter loss. + passivity_weight: Weight of a penalty on predicted S-matrices whose largest + singular value exceeds 1 (or the simulated one, where that is larger). + 0 (default) leaves it out. + above_srf_weight: Loss weight of the frequency points above each geometry's + first self-resonance, relative to the points below it; e.g. 0.1 spends + the model's capacity on the band an inductor is used in. 1 (default) + weights all points alike. The resonance is found in the simulated + S-parameters (orca.training.losses.first_self_resonance). + + The three loss options need a per-point dataset. They change what the + validation loss measures, so losses are only comparable between runs with the + same options. """ super().__init__(name="Model Trainer", index=4) self.model_cls = get_model_class(model) @@ -114,8 +136,15 @@ def __init__( "grad_clip_norm": grad_clip_norm, "allow_tf32": allow_tf32, } + self.admittance_weight = admittance_weight + self.passivity_weight = passivity_weight + self.above_srf_weight = above_srf_weight # Checked now, so a typo fails before hours of tuning rather than after TrainingConfig.from_hyperparameters(self.training_defaults) + if admittance_weight < 0 or passivity_weight < 0: + raise ValueError("admittance_weight and passivity_weight must not be negative.") + if above_srf_weight <= 0: + raise ValueError(f"above_srf_weight must be positive, got {above_srf_weight}.") def run( self, @@ -172,18 +201,14 @@ def run( val_dataset = geometry.dataset.new_split(directory=result_dir, data_df=val_df) geometry.dataset.clear_cache() self._check_compatibility(train_dataset) + criterion_factory = self._criterion_factory(train_dataset) + self._weight_samples(train_dataset) + self._weight_samples(val_dataset) logger.info( f"Loaded {len(train_dataset)} training samples and {len(val_dataset)} validation samples for model training. Beginning training..." ) - trainer = Trainer( - config=TrainingConfig.from_hyperparameters( - {"epochs": self.max_epochs, **self.training_defaults, **hyperparameters} - ), - progress_callback=progress_callback, - stage_name=self.name, - ) spec = train_dataset.io_spec basis = self.basis_cls.from_spec(spec, hyperparameters) if self.basis_cls else None @@ -191,11 +216,16 @@ def run( if seed is not None: torch.manual_seed(seed) - result = trainer.fit( - model=self.model_cls.from_spec(spec, hyperparameters, basis), - train_dataset=train_dataset, - val_dataset=val_dataset, + model = self.model_cls.from_spec(spec, hyperparameters, basis) + trainer = Trainer( + config=TrainingConfig.from_hyperparameters( + {"epochs": self.max_epochs, **self.training_defaults, **hyperparameters} + ), + criterion=criterion_factory(model) if criterion_factory else None, + progress_callback=progress_callback, + stage_name=self.name, ) + result = trainer.fit(model=model, train_dataset=train_dataset, val_dataset=val_dataset) if not math.isfinite(result.best_loss): raise RuntimeError( "Training diverged: the validation loss was never finite. Lower the learning " @@ -251,6 +281,7 @@ def _tune( directory=result_dir, data_df=train_val_df, fit_normalizers=True ) self._check_compatibility(train_val_dataset) + self._weight_samples(train_val_dataset) tuner = HyperparameterTuner( model_cls=self.model_cls, dataset=train_val_dataset, @@ -265,11 +296,50 @@ def _tune( training_defaults=self.training_defaults, # The dataset labels each sample with its result file, the name column of the table holdout_groups=set(val_df["name"]) if self.n_fold_cv == 1 else None, + criterion_factory=self._criterion_factory(train_val_dataset), ) hyperparameters = tuner.tune() logger.info(f"Hyperparameter tuning completed. Best hyperparameters: {hyperparameters}") return hyperparameters + @property + def _custom_loss(self) -> bool: + """Whether any loss option asks for more than the model's default loss.""" + return self.admittance_weight > 0 or self.passivity_weight > 0 or self.above_srf_weight != 1 + + def _criterion_factory( + self, dataset: BaseDataset + ) -> Callable[[OrcaModel], SParameterLoss] | None: + """ + The loss of each model trained on ``dataset``'s outputs, or None for the model's + default loss. A factory, since the data term is the default loss of the model. + """ + if not self._custom_loss: + return None + if type(dataset).frequency_mode is not FrequencyMode.PER_POINT: + raise ValueError( + "ModelTrainer's admittance_weight, passivity_weight and above_srf_weight need " + f"a per-point dataset, but {type(dataset).__name__} lays out frequency as " + f"{type(dataset).frequency_mode.name}." + ) + codec, output_normalizer = dataset.codec, dataset.output_normalizer + + def build(model: OrcaModel) -> SParameterLoss: + return SParameterLoss( + model.default_loss(), + codec=codec, + output_normalizer=output_normalizer, + admittance_weight=self.admittance_weight, + passivity_weight=self.passivity_weight, + ) + + return build + + def _weight_samples(self, dataset: BaseDataset) -> None: + """Weight the frequency points above each geometry's self-resonance, if asked to.""" + if self.above_srf_weight != 1: + dataset.set_sample_weights(srf_sample_weights(dataset, self.above_srf_weight)) + def _save_split( self, context: "PipelineContext", diff --git a/src/orca/simulation/gds_converter.py b/src/orca/simulation/gds_converter.py index 4900b04..be0681d 100644 --- a/src/orca/simulation/gds_converter.py +++ b/src/orca/simulation/gds_converter.py @@ -1,3 +1,4 @@ +import glob import os from contextlib import ExitStack, redirect_stdout from typing import Any @@ -13,7 +14,8 @@ def create_palace_model_from_gds( stackup_xml: str, simconfig_filename: str, show_mesh_results: bool = False, -) -> tuple[str, dict[str, Any], str, str, str]: + ports: list[dict[str, Any]] | None = None, +) -> tuple[str, dict[str, Any], str, str, str, float]: """ Uses gds2palace to create a Palace model from a GDS file and simulation configuration. The simconfig is a json and can either be created manually or by using setupEM GUI and saving the configuration. @@ -28,10 +30,13 @@ def create_palace_model_from_gds( stackup_xml (str): Path to the XML file describing the layer stackup. simconfig_filename (str): Path to the simulation configuration file (json). show_mesh_results (bool): Show the gmsh GUI with the mesh and keep gds2palace's console output. + ports (list[dict[str, Any]] | None): The sample's ports, in the format of the simconfig's + ``ports`` list (``BaseGeometry.ports_for``). None uses the simconfig's own. Returns: - tuple[str, dict, str, str, str]: geometry_name, params, Palace config name, simulation - directory and data directory of the created Palace model. + tuple[str, dict, str, str, str, float]: geometry_name, params, Palace config name, + simulation directory and data directory of the created Palace model, and the worst + element quality of its mesh (:func:`worst_element_quality`). """ # Imported here: the PyPI gmsh wheel loads libGLU, which headless machines (CI runners, # HPC compute nodes) often lack, and `import orca` must not depend on it. @@ -57,9 +62,9 @@ def create_palace_model_from_gds( # The settings dictionary contains all simulation parameters (e.g. frequency range, mesh settings...) settings = simconfig["saved_values"] - # Add all ports from simconfig + # Add the sample's ports; everything else comes from the simconfig simulation_ports = simulation_setup.all_simulation_ports() - for port in simconfig["ports"]: + for port in simconfig["ports"] if ports is None else ports: simulation_ports.add_port( simulation_setup.simulation_port( portnumber=port["portnumber"], @@ -121,4 +126,43 @@ def create_palace_model_from_gds( # for convenience, write run script to model directory utilities.create_run_script(settings["sim_path"]) - return geometry_name, params, config_name, sim_path, data_dir + meshes = glob.glob(os.path.join(sim_path, "*.msh")) + quality = worst_element_quality(meshes[0]) if meshes else float("nan") + return geometry_name, params, config_name, sim_path, data_dir, quality + + +def worst_element_quality(mesh_filename: str) -> float: + """ + The worst shape quality of the volume elements of a gmsh mesh. + + The quality is gmsh's minimum scaled inverse condition number (``minSICN``): 1 for a + regular tetrahedron, 0 for a flat one, negative for an inverted one. A flat element + has no volume, which puts a zero on the diagonal of Palace's system matrix; Palace + then stops with ``HasPositiveFiniteDiagonal(...) is false``. + + Args: + mesh_filename (str): Path to the ``.msh`` file. + + Returns: + float: The lowest quality of any volume element; NaN for a mesh without one, and + -inf for a mesh gmsh cannot read back. gmsh occasionally writes such a corrupt mesh + (elements referencing node 0, which does not exist), which Palace cannot use either. + """ + import gmsh # imported here like in create_palace_model_from_gds (libGLU on HPC nodes) + + gmsh.initialize() + try: + gmsh.option.setNumber("General.Terminal", 0) + try: + gmsh.open(mesh_filename) + except Exception: # noqa: BLE001 - gmsh raises a bare Exception for an unreadable mesh + return float("-inf") + _, tags, _ = gmsh.model.mesh.getElements(3) + lowest = [ + min(gmsh.model.mesh.getElementQualities(element_tags, "minSICN")) + for element_tags in tags + if len(element_tags) + ] + return min(lowest) if lowest else float("nan") + finally: + gmsh.finalize() diff --git a/src/orca/training/README.md b/src/orca/training/README.md index a161d16..c091a0d 100644 --- a/src/orca/training/README.md +++ b/src/orca/training/README.md @@ -18,7 +18,7 @@ architecture lives in a geometry class. | `models/` | `OrcaModel` — architecture, hyperparameter space, loss, guarantees | | `trainer.py` | `Trainer` / `TrainingConfig` — optimizer, schedule, early stopping, history | | `tuner.py` | `HyperparameterTuner` — optuna study with k-fold cross-validation over geometries, or one fixed hold-out (`n_fold_cv=1`) | -| `losses.py` | Loss modules (`ComplexMSELoss`, `MSEPlusLogCoshLoss`) | +| `losses.py` | Loss modules (`ComplexMSELoss`, `MSEPlusLogCoshLoss`, `SParameterLoss`), self-resonance detection and sample weights | The model owns *what* is fitted (architecture and loss); the trainer owns *how* it is fitted. `TrainingConfig` holds the hyperparameters the trainer owns, and @@ -71,6 +71,29 @@ to be re-pointed. Codec and model both declare `PhysicsGuarantees`; the exporter writes their union into the ONNX metadata, so a model exported with the upper-triangle codec is marked `reciprocal` whatever the architecture is. +### Loss terms and sample weights + +A model trains with its `default_loss()` (L1 on the normalized outputs for the MLP). +`SParameterLoss` wraps that data loss with terms that judge the prediction as a circuit, +and `ModelTrainer` builds it from three settings: + +| Setting | Adds | +| --- | --- | +| `admittance_weight` | `admittance_error`: relative L1 error of Y = (I + S)⁻¹(I − S), of all of Y and of its real part, which track L and Q | +| `passivity_weight` | `passivity_violation`: how far σ_max of the predicted S exceeds max(1, σ_max of the simulated S) | +| `above_srf_weight` | per-sample weights from `srf_sample_weights`: points above each geometry's first self-resonance count this much relative to those below | + +The terms are computed per sample on the denormalized S-matrices, which the codec's +`decode_tensor` rebuilds, so the data loss has to be elementwise (a torch loss with a +`reduction` attribute). Sample weights live on the dataset (`set_sample_weights`), are +appended to its `tensors` and batched like the inputs, also through the tuner's `Subset` +folds; the trainer passes each batch's weights to the loss as a third argument. They are +normalized to average 1, so the epoch loss stays on the scale of the data loss. + +`first_self_resonance` needs no knowledge of the ports: it counts the negative eigenvalues +of the susceptance matrix Im(Y), one per inductive mode, and returns where their number +first drops. For InductorOcta this is the differential self-resonance. + ### Adding an output codec ```python @@ -84,6 +107,8 @@ class MyCodec(OutputCodec): def encode(self, ntwk): ... # skrf.Network -> (n_freq, output_dim) def decode(self, raw): ... # (n_freq, output_dim) -> (n_freq, N, N) complex + def decode_tensor(self, raw): ... # the same in torch, differentiable; optional, + # needed by the admittance and passivity losses ``` ### Basis expansions diff --git a/src/orca/training/codecs.py b/src/orca/training/codecs.py index be3d03d..521067b 100644 --- a/src/orca/training/codecs.py +++ b/src/orca/training/codecs.py @@ -15,13 +15,18 @@ from __future__ import annotations from abc import ABC, abstractmethod -from typing import ClassVar +from typing import TYPE_CHECKING, ClassVar import numpy as np import skrf as rf from orca.training.guarantees import PhysicsGuarantees +if TYPE_CHECKING: + # Only for annotations: `import orca` loads this module, and has to work in a + # simulation-only install without PyTorch. + import torch + class OutputCodec(ABC): """Bidirectional mapping between a network response and a model target. @@ -64,6 +69,24 @@ def decode(self, raw: np.ndarray) -> np.ndarray: np.ndarray: Complex array of shape (n_freq, n_ports, n_ports). """ + def decode_tensor(self, raw: torch.Tensor) -> torch.Tensor: + """Differentiable :meth:`decode` of a batch of denormalized outputs, in torch. + + The loss terms that work on the S-matrix (``ModelTrainer(admittance_weight=..., + passivity_weight=...)``) need it; a codec that does not implement it cannot be + trained with them. + + Args: + raw (torch.Tensor): Outputs of shape ``(batch, output_dim)``. + + Returns: + torch.Tensor: Complex S-matrices of shape ``(batch, n_ports, n_ports)``. + """ + raise NotImplementedError( + f"{type(self).__name__} does not implement decode_tensor, which the admittance " + "and passivity loss terms need." + ) + def to_network(self, raw: np.ndarray, frequencies: np.ndarray) -> rf.Network: """Decode model outputs and wrap them in a ``skrf.Network``.""" s = self.decode(np.asarray(raw)) @@ -100,6 +123,12 @@ def decode(self, raw: np.ndarray) -> np.ndarray: raw = np.asarray(raw).reshape(-1, n, n, 2) return (raw[..., 0] + 1j * raw[..., 1]).astype(np.complex64) + def decode_tensor(self, raw: torch.Tensor) -> torch.Tensor: + import torch # Local: see the TYPE_CHECKING import + + pairs = raw.reshape(-1, self.n_ports, self.n_ports, 2) + return torch.complex(pairs[..., 0], pairs[..., 1]) + class UpperTriangleReImCodec(OutputCodec): """The upper triangle of a reciprocal S-matrix, as interleaved real/imag pairs. @@ -167,3 +196,17 @@ def decode(self, raw: np.ndarray) -> np.ndarray: s[:, rows, cols] = entries s[:, cols, rows] = entries # mirror; the diagonal is written twice, harmlessly return s + + def decode_tensor(self, raw: torch.Tensor) -> torch.Tensor: + import torch # Local: see the TYPE_CHECKING import + + n = self.n_ports + pairs = raw.reshape(-1, self.n_entries, 2) + entries = torch.complex(pairs[..., 0], pairs[..., 1]) + # Gathered rather than written into a zero matrix, so autograd sees one plain + # indexing operation: entry (i, j) and entry (j, i) both read triangle entry k. + rows, cols = self._triangle_indices() + entry_of = np.empty((n, n), dtype=np.int64) + entry_of[rows, cols] = entry_of[cols, rows] = np.arange(self.n_entries) + index = torch.as_tensor(entry_of.reshape(-1), device=raw.device) + return entries[:, index].reshape(-1, n, n) diff --git a/src/orca/training/datasets/base_dataset.py b/src/orca/training/datasets/base_dataset.py index 728ec9f..6053d71 100644 --- a/src/orca/training/datasets/base_dataset.py +++ b/src/orca/training/datasets/base_dataset.py @@ -43,6 +43,8 @@ def __init__( # The normalized samples, stacked once loading is done self.inputs = torch.empty(0) self.targets = torch.empty(0) + # Optional loss weight of each sample, see set_sample_weights + self.sample_weights: torch.Tensor | None = None # Parsed Touchstone files, shared by every split made with new_split() self._touchstone_cache: dict[tuple[str, int], tuple[np.ndarray, np.ndarray]] = {} self.input_normalizer = input_normalizer @@ -95,10 +97,30 @@ def clear_cache(self) -> None: self._touchstone_cache.clear() @property - def tensors(self) -> tuple[torch.Tensor, torch.Tensor]: - """All normalized inputs and targets, stacked along the first dimension. + def tensors(self) -> tuple[torch.Tensor, ...]: + """All normalized inputs and targets, stacked along the first dimension, followed + by the sample weights if any are set. The trainer batches all of them alike. """ - return self.inputs, self.targets + if self.sample_weights is None: + return self.inputs, self.targets + return self.inputs, self.targets, self.sample_weights + + def set_sample_weights(self, weights: torch.Tensor | None) -> None: + """ + Weight each sample's contribution to the loss, e.g. with + :func:`~orca.training.losses.srf_sample_weights`. The trainer then hands the + weights of each batch to the loss as a third argument, so the loss has to + accept them (:class:`~orca.training.losses.SParameterLoss` does). + + Args: + weights (torch.Tensor | None): One weight per sample; None removes them. + """ + if weights is not None and weights.shape != (len(self.inputs),): + raise ValueError( + f"Expected one weight per sample, shape ({len(self.inputs)},), " + f"got {tuple(weights.shape)}." + ) + self.sample_weights = weights @property def io_spec(self) -> IOSpec: diff --git a/src/orca/training/losses.py b/src/orca/training/losses.py index b767895..789901a 100644 --- a/src/orca/training/losses.py +++ b/src/orca/training/losses.py @@ -3,15 +3,37 @@ Losses are ``nn.Module`` subclasses so they can be configured, moved between devices and returned from :meth:`~orca.training.models.base_model.OrcaModel.default_loss` like any other torch loss. + +:class:`SParameterLoss` adds optional terms that judge the prediction as a circuit +rather than as a vector of numbers (the relative error of its admittance matrix, a +passivity penalty) and per-sample weights, such as those of :func:`srf_sample_weights`, +which weight the frequency points above each geometry's self-resonance differently. """ from __future__ import annotations +import copy import math +from typing import TYPE_CHECKING, Any +import numpy as np import torch from torch import nn +from orca.logger import logger +from orca.training.spec import FrequencyMode + +if TYPE_CHECKING: + from collections.abc import Callable + + from orca.training.codecs import OutputCodec + from orca.training.datasets.base_dataset import BaseDataset + from orca.training.normalize import Normalizer + +#: Conductance floor of :func:`admittance_error`, relative to the magnitude of the +#: admittance matrix: keeps the conductance term finite where the conductance vanishes. +CONDUCTANCE_FLOOR = 1e-2 + class ComplexMSELoss(nn.Module): """Mean squared error over interleaved real/imaginary output pairs. @@ -51,3 +73,268 @@ def forward(self, pred: torch.Tensor, target: torch.Tensor) -> torch.Tensor: abs_error + torch.nn.functional.softplus(-2.0 * abs_error) - math.log(2.0) ) return self.log_cosh_weight * log_cosh + self.mse_weight * self.mse(pred, target) + + +def s_to_y(s: torch.Tensor) -> torch.Tensor: + """Admittance matrices of a batch of S-matrices, in units of the reference admittance. + + ``y = (I + S)^-1 (I - S)``, which is ``Z0 * Y`` for a common, real reference + impedance ``Z0``. The two factors commute, so the inverse is applied by a solve. + + Args: + s (torch.Tensor): Complex S-matrices of shape ``(batch, n_ports, n_ports)``. + + Returns: + torch.Tensor: Normalized admittance matrices of the same shape. + """ + eye = torch.eye(s.shape[-1], dtype=s.dtype, device=s.device) + return torch.linalg.solve(eye + s, eye - s) + + +def admittance_error(s_pred: torch.Tensor, s_true: torch.Tensor) -> torch.Tensor: + """Relative error of the predicted admittance matrix, per sample. + + A small error in S can be a large error in the inductance and quality factor: at + low frequency an inductor is close to a short between its ports, where Y = 1/(R + + jwL) reacts strongly to S. The sum of two relative L1 errors over all entries: + + * of Y as a whole, which follows the reactance (L, C, the self-resonance), and + * of its real part, the conductance, which sets the losses and with them Q. It is + a few percent of |Y| for a high-Q inductor, so the first term barely sees it. + + Both are scaled by the reference matrix of each sample, so every frequency point + counts alike whatever the magnitude of its Y. The conductance is floored at + :data:`CONDUCTANCE_FLOOR` of that magnitude. + + Args: + s_pred (torch.Tensor): Predicted complex S-matrices, ``(batch, n, n)``. + s_true (torch.Tensor): Reference complex S-matrices, same shape. + + Returns: + torch.Tensor: The error of each sample, shape ``(batch,)``. + """ + y_true = s_to_y(s_true) + error = s_to_y(s_pred) - y_true + + def entry_sum(x: torch.Tensor) -> torch.Tensor: + return x.abs().sum(dim=(-2, -1)) + + scale = entry_sum(y_true) + 1e-12 + conductance = entry_sum(y_true.real) + CONDUCTANCE_FLOOR * scale + return entry_sum(error) / scale + entry_sum(error.real) / conductance + + +def passivity_violation(s_pred: torch.Tensor, s_true: torch.Tensor) -> torch.Tensor: + """How far the largest singular value of each predicted S-matrix exceeds 1. + + A passive network has ``sigma_max(S) <= 1``. The reference may exceed that + slightly, as the DC point extrapolated by ``touchstone_type="dc_deembedded"`` + does; a prediction is only penalised beyond ``max(1, sigma_max(S_true))``, so the + penalty never pulls the model away from its own training data. + + Args: + s_pred (torch.Tensor): Predicted complex S-matrices, ``(batch, n, n)``. + s_true (torch.Tensor): Reference complex S-matrices, same shape. + + Returns: + torch.Tensor: The violation of each sample, zero where passive; ``(batch,)``. + """ + allowed = largest_singular_value(s_true).clamp(min=1.0) + return torch.relu(largest_singular_value(s_pred) - allowed) + + +def largest_singular_value(s: torch.Tensor) -> torch.Tensor: + """``sigma_max`` of each matrix of a batch, as the root of the top eigenvalue of S^H S. + + Equal to ``torch.linalg.matrix_norm(s, ord=2)``, but about five times faster for a + batch of small matrices on a GPU, where batched SVDs are slow. + """ + top = torch.linalg.eigvalsh(s.mH @ s)[..., -1] + # Floored, so the square root has a finite gradient for an all-zero matrix + return top.clamp(min=1e-12).sqrt() + + +class SParameterLoss(nn.Module): + """A model's data loss, plus optional circuit-level terms and per-sample weights. + + Every term is computed per sample and summed with its weight; the per-sample sum + is then averaged, weighted by ``sample_weights`` if given. Those weights should + average 1 over a dataset (as :func:`srf_sample_weights` returns them), so the + loss of a whole epoch is their weighted mean and stays on the scale of the data + loss alone. + + Args: + data_loss (Callable): Elementwise loss on the normalized outputs, usually + the model's :meth:`~orca.training.models.base_model.OrcaModel.default_loss`. + It needs a ``reduction`` attribute, as torch's own losses have; a copy + set to ``"none"`` is used. + codec (OutputCodec): Layout of the outputs; its + :meth:`~orca.training.codecs.OutputCodec.decode_tensor` rebuilds the S-matrices. + output_normalizer (Normalizer | None): Undoes the output normalization before + the S-matrices are rebuilt. + admittance_weight (float): Weight of :func:`admittance_error`; 0 leaves it out. + passivity_weight (float): Weight of :func:`passivity_violation`; 0 leaves it out. + """ + + def __init__( + self, + data_loss: Callable[..., torch.Tensor], + codec: OutputCodec, + output_normalizer: Normalizer | None = None, + admittance_weight: float = 0.0, + passivity_weight: float = 0.0, + ): + super().__init__() + if not hasattr(data_loss, "reduction"): + raise TypeError( + f"SParameterLoss needs an elementwise data loss with a 'reduction' attribute " + f"(such as nn.L1Loss), got {type(data_loss).__name__}." + ) + if admittance_weight < 0 or passivity_weight < 0: + raise ValueError("Loss term weights must not be negative.") + # Typed Any: the reduction attribute is torch's convention, not part of Callable + elementwise: Any = copy.copy(data_loss) + elementwise.reduction = "none" + self.data_loss = elementwise + self.codec = codec + self.output_normalizer = output_normalizer + self.admittance_weight = admittance_weight + self.passivity_weight = passivity_weight + + def forward( + self, + pred: torch.Tensor, + target: torch.Tensor, + sample_weights: torch.Tensor | None = None, + ) -> torch.Tensor: + loss = self.data_loss(pred, target).reshape(len(pred), -1).mean(dim=1) + if self.admittance_weight > 0 or self.passivity_weight > 0: + s_pred, s_true = self._s_matrices(pred), self._s_matrices(target) + if self.admittance_weight > 0: + loss = loss + self.admittance_weight * admittance_error(s_pred, s_true) + if self.passivity_weight > 0: + loss = loss + self.passivity_weight * passivity_violation(s_pred, s_true) + if sample_weights is not None: + loss = loss * sample_weights + return loss.mean() + + def _s_matrices(self, outputs: torch.Tensor) -> torch.Tensor: + if self.output_normalizer is not None: + outputs = self.output_normalizer.denormalize(outputs) + return self.codec.decode_tensor(outputs) + + +def first_self_resonance(frequencies: np.ndarray, s: np.ndarray) -> float: + """The first self-resonance frequency of a network, from its S-parameters. + + The frequency at which the first inductive mode of the network turns capacitive: + the susceptance matrix ``Im(Y)`` (its Hermitian part) has one negative + eigenvalue per inductive mode, and the resonance is where their number first + drops, interpolated linearly. This needs no knowledge of the ports: for a + center-tapped inductor it finds the differential self-resonance (within the + frequency step, on the InductorOcta results), since the differential mode + resonates first. The modes are counted at the lowest frequency above DC, where Y + is still almost real. + + Args: + frequencies (np.ndarray): Frequencies in Hz, shape ``(n_freq,)``, any order. + s (np.ndarray): Complex S-matrices at those frequencies, ``(n_freq, n, n)``. + + Returns: + float: The resonance in Hz; ``inf`` if there is none in the band, or if the + network has no inductive mode at its lowest frequency. + """ + order = np.argsort(frequencies) + frequencies, s = frequencies[order], s[order] + eye = np.eye(s.shape[-1]) + try: + y = np.linalg.solve(eye + s, eye - s) + except np.linalg.LinAlgError: + logger.warning("Cannot convert an S-matrix to Y (a short circuit); no resonance found.") + return math.inf + susceptance = (y - np.conj(np.swapaxes(y, -1, -2))) / 2j + eigenvalues = np.linalg.eigvalsh(susceptance) # ascending, per frequency + # Relative to the largest mode, so rounding noise around zero is not a mode + threshold = -1e-6 * np.abs(eigenvalues).max(axis=1) + n_inductive = (eigenvalues < threshold[:, None]).sum(axis=1) + + ac = np.flatnonzero(frequencies > 0) + if ac.size < 2 or n_inductive[ac[0]] == 0: + return math.inf + reference = n_inductive[ac[0]] + dropped = ac[n_inductive[ac] < reference] + if dropped.size == 0: + return math.inf + + k = dropped[0] + # The eigenvalue that turned is the reference-th smallest one + before, after = eigenvalues[k - 1, reference - 1], eigenvalues[k, reference - 1] + f0, f1 = frequencies[k - 1], frequencies[k] + if after <= before: + return float(f1) + return float(np.clip(f0 - before * (f1 - f0) / (after - before), f0, f1)) + + +def srf_sample_weights(dataset: BaseDataset, above_srf_weight: float) -> torch.Tensor: + """Per-sample loss weights that weight frequencies above the self-resonance apart. + + Below its first self-resonance (:func:`first_self_resonance`) a passive is used + as what it is, an inductor say; above it the response swings through resonances + that the network has to spend capacity on, but that a designer rarely needs. + Each frequency point of a geometry above that geometry's resonance is weighted + ``above_srf_weight`` relative to the points below it. Geometries that do not + resonate in the band keep weight 1 everywhere. + + Args: + dataset (BaseDataset): A loaded per-point dataset with a ``frequency`` input. + above_srf_weight (float): Weight of the points above the resonance, relative + to those below; positive. + + Returns: + torch.Tensor: One weight per sample, on the dataset's device, averaging 1. + """ + if type(dataset).frequency_mode is not FrequencyMode.PER_POINT: + raise ValueError( + f"Weighting by self-resonance needs a per-point dataset, but " + f"{type(dataset).__name__} lays out frequency as {type(dataset).frequency_mode.name}." + ) + if "frequency" not in dataset.input_param_names: + raise ValueError( + f"Weighting by self-resonance needs a 'frequency' input, but the inputs are " + f"{dataset.input_param_names}." + ) + if above_srf_weight <= 0: + raise ValueError(f"above_srf_weight must be positive, got {above_srf_weight}.") + + inputs, targets = dataset.inputs, dataset.targets + if dataset.input_normalizer is not None: + inputs = dataset.input_normalizer.denormalize(inputs) + if dataset.output_normalizer is not None: + targets = dataset.output_normalizer.denormalize(targets) + column = dataset.input_param_names.index("frequency") + frequencies = inputs[:, column].double().cpu().numpy() + s = dataset.codec.decode(targets.cpu().numpy()) + + # The samples of each geometry, found by one sort rather than a scan per geometry + _, group_of, counts = np.unique( + np.asarray(dataset.sample_groups), return_inverse=True, return_counts=True + ) + order = np.argsort(group_of, kind="stable") + resonance = np.array( + [ + first_self_resonance(frequencies[members], s[members]) + for members in np.split(order, np.cumsum(counts)[:-1]) + ] + ) + + above = frequencies > resonance[group_of] + weights = np.where(above, above_srf_weight, 1.0) + weights /= weights.mean() + in_band = resonance[np.isfinite(resonance)] + logger.info( + f"Self-resonance found in the band for {in_band.size} of {resonance.size} geometries" + + (f" (median {np.median(in_band) / 1e9:.3g} GHz)" if in_band.size else "") + + f"; {above.mean():.1%} of the frequency points lie above it and are weighted " + f"{above_srf_weight:g} relative to the rest." + ) + return torch.as_tensor(weights, dtype=torch.float32, device=dataset.inputs.device) diff --git a/src/orca/training/trainer.py b/src/orca/training/trainer.py index 6d9b960..bcea9d7 100644 --- a/src/orca/training/trainer.py +++ b/src/orca/training/trainer.py @@ -178,47 +178,47 @@ def epochs_run(self) -> int: return len(self.history) -def _stacked_tensors(dataset: Dataset) -> tuple[torch.Tensor, torch.Tensor] | None: - """The inputs and targets of ``dataset`` as two tensors, if it keeps them stacked. +def _stacked_tensors(dataset: Dataset) -> tuple[torch.Tensor, ...] | None: + """The inputs, targets and any per-sample extras of ``dataset`` as stacked tensors. Datasets exposing ``tensors`` (ORCA's datasets and ``TensorDataset``) and ``Subset`` - views of them qualify; anything else returns ``None``. + views of them qualify; anything else returns ``None``. ORCA's datasets append their + sample weights, if set, as a third tensor. """ if isinstance(dataset, Subset): stacked = _stacked_tensors(dataset.dataset) if stacked is None: return None indices = torch.as_tensor(dataset.indices, dtype=torch.long, device=stacked[0].device) - return stacked[0][indices], stacked[1][indices] + return tuple(tensor[indices] for tensor in stacked) tensors = getattr(dataset, "tensors", None) - if isinstance(tensors, tuple) and len(tensors) == 2: # (inputs, targets) + if isinstance(tensors, tuple) and len(tensors) >= 2: # (inputs, targets, *extras) return tensors return None class _TensorBatches: - """Mini-batches cut from two stacked tensors.""" + """Mini-batches cut from stacked tensors of equal length, all sliced alike.""" - def __init__(self, inputs: torch.Tensor, targets: torch.Tensor, batch_size: int, shuffle: bool): - self.inputs = inputs - self.targets = targets + def __init__(self, tensors: tuple[torch.Tensor, ...], batch_size: int, shuffle: bool): + self.tensors = tensors self.batch_size = batch_size self.shuffle = shuffle def __len__(self) -> int: - return math.ceil(len(self.inputs) / self.batch_size) + return math.ceil(len(self.tensors[0]) / self.batch_size) - def __iter__(self) -> Iterator[tuple[torch.Tensor, torch.Tensor]]: - n = len(self.inputs) + def __iter__(self) -> Iterator[tuple[torch.Tensor, ...]]: + n = len(self.tensors[0]) if self.shuffle: - order = torch.randperm(n, device=self.inputs.device) + order = torch.randperm(n, device=self.tensors[0].device) for start in range(0, n, self.batch_size): batch = order[start : start + self.batch_size] - yield self.inputs[batch], self.targets[batch] + yield tuple(tensor[batch] for tensor in self.tensors) else: for start in range(0, n, self.batch_size): end = start + self.batch_size - yield self.inputs[start:end], self.targets[start:end] + yield tuple(tensor[start:end] for tensor in self.tensors) @contextlib.contextmanager @@ -238,7 +238,7 @@ def make_batches( """Mini-batches of ``dataset``, sliced from stacked tensors where it has them.""" stacked = _stacked_tensors(dataset) if stacked is not None: - return _TensorBatches(*stacked, batch_size=batch_size, shuffle=shuffle) + return _TensorBatches(stacked, batch_size=batch_size, shuffle=shuffle) return DataLoader(dataset, batch_size=batch_size, shuffle=shuffle) @@ -296,7 +296,9 @@ class Trainer: Args: config (TrainingConfig | None): Training hyperparameters. Defaults are used if omitted. - criterion (Callable | None): Loss to optimize. If omitted, the model's + criterion (Callable | None): Loss to optimize, called as + ``criterion(prediction, target)``, or as ``criterion(prediction, target, + sample_weights)`` for a dataset with sample weights. If omitted, the model's :meth:`~orca.training.models.base_model.OrcaModel.default_loss` is used. progress_callback (Callable | None): Called as ``(stage_name, current_epoch, total_epochs, message)`` after every epoch. @@ -307,7 +309,7 @@ class Trainer: def __init__( self, config: TrainingConfig | None = None, - criterion: Callable[[torch.Tensor, torch.Tensor], torch.Tensor] | None = None, + criterion: Callable[..., torch.Tensor] | None = None, progress_callback: Callable[[str, int, int, str], None] | None = None, stage_name: str = "Training", verbose: bool = True, @@ -428,14 +430,12 @@ def _train_epoch(self, model, criterion, optimizer, schedule, loader) -> float: total = torch.zeros((), device=self.config.device) count = 0 - for batch_x, batch_y in tqdm.tqdm( - loader, desc="Training", leave=False, disable=not self.verbose - ): - x = batch_x.to(self.config.device) - y = batch_y.to(self.config.device) + for batch in tqdm.tqdm(loader, desc="Training", leave=False, disable=not self.verbose): + # Inputs, targets, and the sample weights if the dataset has them + x, y, *weights = (tensor.to(self.config.device) for tensor in batch) optimizer.zero_grad() - loss = criterion(model(x), y) + loss = criterion(model(x), y, *weights) loss.backward() if self.config.grad_clip_norm is not None: torch.nn.utils.clip_grad_norm_(model.parameters(), self.config.grad_clip_norm) @@ -461,10 +461,9 @@ def _run_eval(self, model, criterion, loader, desc: str, show_progress: bool = T iterator = tqdm.tqdm(loader, desc=desc, leave=False, disable=not self.verbose) with torch.no_grad(): - for batch_x, batch_y in iterator: - x = batch_x.to(self.config.device) - y = batch_y.to(self.config.device) - total += criterion(model(x), y) * len(x) + for batch in iterator: + x, y, *weights = (tensor.to(self.config.device) for tensor in batch) + total += criterion(model(x), y, *weights) * len(x) count += len(x) return total.item() / count diff --git a/src/orca/training/tuner.py b/src/orca/training/tuner.py index e787387..a7a4d9a 100644 --- a/src/orca/training/tuner.py +++ b/src/orca/training/tuner.py @@ -22,7 +22,7 @@ from orca.training.trainer import EpochResult, Trainer, TrainingConfig if TYPE_CHECKING: - from collections.abc import Collection + from collections.abc import Callable, Collection from orca.training.basis_expansion import BasisExpansion from orca.training.datasets.base_dataset import BaseDataset @@ -137,6 +137,9 @@ class HyperparameterTuner: (see :class:`~orca.training.trainer.TrainingConfig`). holdout_groups (Collection[str] | None): Groups (result files) to validate on when ``n_fold_cv`` is 1, e.g. the validation split of the final training. + criterion_factory (Callable | None): Builds the loss of each trial from its + freshly built model, e.g. an :class:`~orca.training.losses.SParameterLoss` + around the model's default loss. ``None`` trains with the model's default loss. """ def __init__( @@ -155,6 +158,7 @@ def __init__( regularization: bool = False, training_defaults: dict[str, Any] | None = None, holdout_groups: Collection[str] | None = None, + criterion_factory: Callable[[OrcaModel], Callable[..., torch.Tensor]] | None = None, ): self.model_cls = model_cls self.dataset = dataset @@ -170,6 +174,7 @@ def __init__( self.batch_sizes = batch_sizes self.regularization = regularization self.training_defaults = training_defaults or {} + self.criterion_factory = criterion_factory self._deadline = math.inf self.study: optuna.Study | None = None # The folds depend only on the data, so every trial is scored on the same split @@ -263,7 +268,13 @@ def report(epoch: EpochResult) -> None: basis = self.basis_cls.from_spec(self.spec, hyperparameters) if self.basis_cls else None model = self.model_cls.from_spec(self.spec, hyperparameters, basis) - trainer = Trainer(config=config, stage_name=f"Tuning ({fold_label})", verbose=False) + criterion = self.criterion_factory(model) if self.criterion_factory else None + trainer = Trainer( + config=config, + criterion=criterion, + stage_name=f"Tuning ({fold_label})", + verbose=False, + ) try: result = trainer.fit( model=model, diff --git a/src/orca/utils/postprocessing.py b/src/orca/utils/postprocessing.py index 24a6c5d..690c8ff 100644 --- a/src/orca/utils/postprocessing.py +++ b/src/orca/utils/postprocessing.py @@ -1,3 +1,5 @@ +from collections.abc import Sequence + import matplotlib.pyplot as plt import numpy as np import skrf as rf @@ -5,70 +7,116 @@ from orca.logger import logger -def to_mixed_mode(ntwk): - """Mixed-mode view of a single-ended network, as a copy. - - Returned separately rather than alongside the electrical parameters: callers - iterate over that dict computing per-curve errors, and a Network is not a - curve. +def differential_impedance( + ntwk: rf.Network, pairs: Sequence[tuple[int, int]], shorted: Sequence[int] = () +) -> np.ndarray: """ - mm_ntwk = ntwk.copy() - if ntwk.nports >= 4: - mm_ntwk.se2gmm(p=ntwk.nports // 2) - return mm_ntwk + The Z-matrix of port pairs driven differentially, each by a floating source. + The ports in ``shorted`` are AC-grounded (a center tap in differential operation); + every other port that is in no pair is left open. Entry ``(a, b)`` is the voltage + across pair ``a`` per current driven through pair ``b`` (into its first port, out of + its second): ``Z[pa, pb] - Z[pa, nb] - Z[na, pb] + Z[na, nb]`` of the network with + the shorted ports removed. -def calculate_electrical_parameters(ntwk): - # 1. Single-ended to Mixed-Mode Conversion - mm_ntwk = to_mixed_mode(ntwk) + Args: + ntwk (rf.Network): The single-ended network. + pairs: ``(positive, negative)`` port indices, zero-based, per differential port. + shorted: Zero-based indices of the ports to short. - freq_ghz = ntwk.f / 1e9 - omega = 2 * np.pi * ntwk.f + Returns: + np.ndarray: Complex, shape ``(n_freq, len(pairs), len(pairs))``. + """ + kept = [port for port in range(ntwk.nports) if port not in shorted] + position = {port: i for i, port in enumerate(kept)} + # Shorting a port is dropping its row and column of Y; the inverse is then the + # Z-matrix with that port grounded and the rest open. A pseudo-inverse, since a + # winding without a path to ground has a singular Y (its common mode is undefined), + # while its differential impedance, orthogonal to that mode, is not. + z = np.linalg.pinv(ntwk.y[:, kept][:, :, kept]) + z_diff = np.empty((len(ntwk.f), len(pairs), len(pairs)), dtype=complex) + for a, (pa, na) in enumerate(pairs): + for b, (pb, nb) in enumerate(pairs): + i, j, k, m = position[pa], position[na], position[pb], position[nb] + z_diff[:, a, b] = z[:, i, k] - z[:, i, m] - z[:, j, k] + z[:, j, m] + return z_diff + + +def inductor_parameters( + ntwk: rf.Network, ends: tuple[int, int] = (0, 1), shorted: Sequence[int] = () +) -> dict[str, np.ndarray]: + """ + Differential inductance, resistance, quality factor and self-resonance of an inductor. - # Extract Differential Z-parameters for lumped metrics - # Index 0 = Primary Diff (d1), Index 1 = Secondary Diff (d2) - z_d11 = mm_ntwk.z[:, 0, 0] - z_d22 = mm_ntwk.z[:, 1, 1] - z_d12 = mm_ntwk.z[:, 0, 1] + Args: + ntwk (rf.Network): The single-ended network. + ends: Zero-based ports at the two ends of the winding. + shorted: Ports to AC-ground, such as a center tap. - # Calculate Parameters + Returns: + dict: ``L`` [nH], ``R`` [Ohm] and ``Q`` over frequency, and ``srf_f`` [GHz], + where the reactance first turns capacitive (NaN if it stays inductive). + """ + z = differential_impedance(ntwk, [ends], shorted)[:, 0, 0] with np.errstate(divide="ignore", invalid="ignore"): - Lp = np.imag(z_d11) / omega * 1e9 - Ls = np.imag(z_d22) / omega * 1e9 - Rp, Rs = np.real(z_d11), np.real(z_d22) - Qp = np.imag(z_d11) / np.real(z_d11) - Qs = np.imag(z_d22) / np.real(z_d22) - k = np.abs(np.imag(z_d12)) / np.sqrt(np.abs(np.imag(z_d11) * np.imag(z_d22))) + L = np.imag(z) / (2 * np.pi * ntwk.f) * 1e9 + Q = np.imag(z) / np.real(z) + return { + "L": L, + "R": np.real(z), + "Q": Q, + "srf_f": np.array(_first_inductive_to_capacitive(ntwk.f / 1e9, np.imag(z))), + } - im = np.imag(z_d11) - cross = np.where(im[:-1] * im[1:] < 0)[0] # actual sign changes only - f_min = 20.0 - cross = cross[freq_ghz[cross] >= f_min] +def transformer_parameters( + ntwk: rf.Network, + primary: tuple[int, int], + secondary: tuple[int, int], + shorted: Sequence[int] = (), +) -> dict[str, np.ndarray]: + """ + Inductance, resistance and quality factor of both windings of a transformer, their + coupling factor, and the primary's self-resonance, all differential. - if cross.size == 0: - # NaN rather than None, so callers can do arithmetic on the result and - # filter it out with the usual isfinite() masks. - srf_f = float("nan") - else: - # Note: this index must not be called k - that name holds the coupling - # coefficient computed above. - cross_idx = cross[0] - f0, f1 = freq_ghz[cross_idx], freq_ghz[cross_idx + 1] - y0, y1 = im[cross_idx], im[cross_idx + 1] - srf_f = float(f0 - y0 * (f1 - f0) / (y1 - y0)) + Args: + ntwk (rf.Network): The single-ended network. + primary: Zero-based ports at the two ends of the primary winding. + secondary: The same for the secondary winding. + shorted: Ports to AC-ground, such as the center taps. - return { - "Lp": np.array(Lp), - "Ls": np.array(Ls), - "Rp": np.array(Rp), - "Rs": np.array(Rs), - "Qp": np.array(Qp), - "Qs": np.array(Qs), - "k": np.array(k), - "z_d11": np.array(z_d11), - "srf_f": np.array(srf_f), - } + Returns: + dict: ``Lp``, ``Ls`` [nH], ``Rp``, ``Rs`` [Ohm], ``Qp``, ``Qs`` and ``k`` over + frequency, and ``srf_f`` [GHz] of the primary (NaN if it stays inductive). + """ + z = differential_impedance(ntwk, [primary, secondary], shorted) + omega = 2 * np.pi * ntwk.f + zp, zs, zm = z[:, 0, 0], z[:, 1, 1], z[:, 0, 1] + with np.errstate(divide="ignore", invalid="ignore"): + return { + "Lp": np.imag(zp) / omega * 1e9, + "Ls": np.imag(zs) / omega * 1e9, + "Rp": np.real(zp), + "Rs": np.real(zs), + "Qp": np.imag(zp) / np.real(zp), + "Qs": np.imag(zs) / np.real(zs), + "k": np.abs(np.imag(zm)) / np.sqrt(np.abs(np.imag(zp) * np.imag(zs))), + "srf_f": np.array(_first_inductive_to_capacitive(ntwk.f / 1e9, np.imag(zp))), + } + + +def _first_inductive_to_capacitive(freq_ghz: np.ndarray, reactance: np.ndarray) -> float: + """ + Where a reactance first falls from positive to negative, interpolated linearly; NaN if + never. The DC point is skipped: its reactance is zero up to rounding, of either sign. + """ + cross = np.flatnonzero((freq_ghz[:-1] > 0) & (reactance[:-1] > 0) & (reactance[1:] <= 0)) + if cross.size == 0: + return float("nan") + i = cross[0] + f0, f1 = freq_ghz[i], freq_ghz[i + 1] + x0, x1 = reactance[i], reactance[i + 1] + return float(f0 - x0 * (f1 - f0) / (x1 - x0)) def pointwise_relative_error(pred, gt) -> np.ndarray: @@ -120,83 +168,55 @@ def median_relative_error(pred, gt) -> float: return float(np.median(errors)) if errors.size else float("nan") -def plot_rfic_transformer_metrics(ntwk): - metrics = calculate_electrical_parameters(ntwk) - mm_ntwk = to_mixed_mode(ntwk) - freq = ntwk.f / 1e9 - Lp, Ls = metrics["Lp"], metrics["Ls"] - Rp, Rs = metrics["Rp"], metrics["Rs"] - Qp, Qs = metrics["Qp"], metrics["Qs"] - k = metrics["k"] - z_d11 = metrics["z_d11"] - - # Setup Plot - fig, axes = plt.subplots(3, 2, figsize=(14, 10)) - fig.suptitle( - f"RFIC Transformer Report: {ntwk.name}", fontsize=16, fontweight="bold" - ) +def plot_electrical_parameters( + frequencies: np.ndarray, + reference: dict[str, np.ndarray], + predicted: dict[str, np.ndarray] | None = None, + title: str = "", +) -> None: + """ + Show each electrical parameter over frequency, the reference against a prediction. - # Subplot 1: S-Parameters (S11 & S21 Mixed-Mode) - axes[0, 0].plot( - freq, - mm_ntwk.s_db[:, 1, 0], - label="Sdd21 (Insertion Loss)", - color="teal", - lw=2.5, - ) - axes[0, 0].plot( - freq, - mm_ntwk.s_db[:, 0, 0], - label="Sdd11 (Return Loss)", - color="darkorange", - ls="--", - ) - axes[0, 0].set_title("Mixed-Mode S-Parameters", fontsize=14) - axes[0, 0].set_ylabel("Magnitude [dB]") - axes[0, 0].legend() - - # Subplot 2: Inductance (Lp & Ls) - axes[0, 1].plot(freq, Lp, label="Lp (Primary)", color="blue") - axes[0, 1].plot(freq, Ls, label="Ls (Secondary)", color="cyan") - axes[0, 1].set_title("Inductance [nH]", fontsize=14) - axes[0, 1].set_ylabel("L [nH]") - axes[0, 1].legend() - - # Subplot 3: Quality Factor (Qp & Qs) - axes[1, 0].plot(freq, Qp, label="Qp (Primary)", color="red") - axes[1, 0].plot(freq, Qs, label="Qs (Secondary)", color="magenta") - axes[1, 0].set_title("Quality Factor (Q)", fontsize=14) - axes[1, 0].set_ylabel("Q") - axes[1, 0].legend() - - # Subplot 4: Resistance (Rp & Rs) - axes[1, 1].plot(freq, Rp, label="Rp (Primary)", color="darkgreen") - axes[1, 1].plot(freq, Rs, label="Rs (Secondary)", color="lime") - axes[1, 1].set_title("Loss / Resistance [Ω]", fontsize=14) - axes[1, 1].set_ylabel("R [Ω]") - axes[1, 1].legend() - - # Subplot 5: Coupling Coefficient (k) - axes[2, 0].plot(freq, k, color="purple", lw=2) - axes[2, 0].set_title("Coupling Coefficient (k)", fontsize=14) - axes[2, 0].set_ylabel("k") - axes[2, 0].set_ylim(0, 1.1) - - # Subplot 6: Reactance & SRF Identification - axes[2, 1].plot(freq, np.imag(z_d11), label="Im(Zdd11)", color="brown") - axes[2, 1].axhline(0, color="black", lw=1) # y=0 line to find zero-crossing - srf_f = metrics["srf_f"] - if np.isfinite(srf_f): - axes[2, 1].axvline( - srf_f, color="red", linestyle=":", label=f"SRF: {srf_f:.2f} GHz" - ) - axes[2, 1].set_title("Primary Reactance & SRF", fontsize=14) - axes[2, 1].set_ylabel("Im(Z) [Ω]") - axes[2, 1].legend() + One panel per curve; scalar parameters (such as ``srf_f``) are listed in the title + and marked as a vertical line. Shown with pyplot, for interactive use. - for ax in axes.flat: + Args: + frequencies (np.ndarray): Frequencies in Hz. + reference (dict): Parameters of the reference, as a geometry's + ``electrical_parameters`` returns them. + predicted (dict | None): The same for the prediction. + title (str): Figure title. + """ + curves = [name for name, value in reference.items() if np.ndim(value) == 1] + scalars = [name for name, value in reference.items() if np.ndim(value) == 0] + if not curves: + return + freq = frequencies / 1e9 + n_columns = min(3, len(curves)) + n_rows = -(-len(curves) // n_columns) + fig, axes = plt.subplots( + n_rows, n_columns, figsize=(4.8 * n_columns, 3.4 * n_rows), squeeze=False + ) + labels = [ + f"{name}: {float(reference[name]):.3g}" + + (f" (predicted {float(predicted[name]):.3g})" if predicted and name in predicted else "") + for name in scalars + ] + fig.suptitle(" | ".join([title, *labels]) if labels else title, fontsize=13) + + for ax, name in zip(axes.flat, curves, strict=False): + ax.plot(freq, reference[name], color="black", lw=2, label="Reference") + if predicted is not None and name in predicted: + ax.plot(freq, predicted[name], color="tab:blue", ls="--", label="Predicted") + for scalar in scalars: + if np.isfinite(reference[scalar]): + ax.axvline(float(reference[scalar]), color="red", ls=":", lw=1) + ax.set_title(name) ax.set_xlabel("Frequency [GHz]") ax.grid(True, alpha=0.3) + axes.flat[0].legend() + for ax in list(axes.flat)[len(curves) :]: + ax.set_visible(False) plt.tight_layout() plt.show() diff --git a/tests/test_drc_ports.py b/tests/test_drc_ports.py new file mode 100644 index 0000000..7d2f698 --- /dev/null +++ b/tests/test_drc_ports.py @@ -0,0 +1,212 @@ +"""Tests for the port contact check: every gds2palace port marker must touch its metals.""" + +from __future__ import annotations + +from dataclasses import dataclass + +import klayout.db as kdb +import pandas as pd +import pytest + +from orca.geometry.drc import PortContact, check_gds_file, check_ports, port_contacts +from orca.geometry.layers import SG13G2 +from orca.geometry.presets import InductorOcta, StackupXML, TransformerOcta +from orca.pipeline.context import PipelineContext +from orca.pipeline.drc_stage import DRCChecker +from orca.simulation.simulate import read_simconfig + +#: A vertical port from Metal5 up to TopMetal2, marked on layer 201 +PORT = PortContact(number=1, marker_layer=201, metals=(("Metal5", 67), ("TopMetal2", 134))) + + +def _layout(marker_x_um: float, marker: str = "path") -> kdb.Layout: + """A TopMetal2 feed ending at x = 0 over a Metal5 bar starting there, plus a port marker.""" + layout = kdb.Layout() + top = layout.create_cell("top") + um = lambda value: round(value / layout.dbu) # noqa: E731 - one-line unit helper + top.shapes(layout.layer(*SG13G2.TopMetal2)).insert(kdb.Box(um(-20), um(-2), 0, um(2))) + top.shapes(layout.layer(*SG13G2.Metal5)).insert(kdb.Box(0, um(-10), um(10), um(10))) + x = um(marker_x_um) + if marker == "path": + line = kdb.Path([kdb.Point(x, um(-2)), kdb.Point(x, um(2))], 0) + top.shapes(layout.layer(201, 0)).insert(line) + else: + top.shapes(layout.layer(201, 0)).insert(kdb.Box(x - um(0.1), um(-2), x, um(2))) + return layout + + +@pytest.mark.parametrize("marker", ["path", "box"]) +def test_marker_on_the_feed_end_is_connected(marker): + assert check_ports(_layout(0.0, marker), (PORT,)) == {} + + +@pytest.mark.parametrize("marker", ["path", "box"]) +def test_marker_5nm_off_the_feed_end_is_open(marker): + # Beyond the feed end, over the ground bar: it touches Metal5 but not the feed + assert check_ports(_layout(0.105, marker), (PORT,)) == {"port1.open_TopMetal2": 1} + + +def test_missing_marker_is_reported(): + layout = _layout(0.0) + layout.clear_layer(layout.find_layer(201, 0)) + + assert check_ports(layout, (PORT,)) == {"port1.missing": 1} + + +def test_port_contacts_follow_the_simconfig_and_stackup(): + geometry = TransformerOcta() + contacts = port_contacts( + read_simconfig(geometry.simconfig_filename)["ports"], geometry.stackup_xml + ) + + assert [c.number for c in contacts] == [1, 2, 3, 4, 5, 6] + assert contacts[0] == PORT + assert contacts[2].metals == (("Metal5", 67), ("TopMetal1", 126)) + + +def test_port_on_a_layer_the_stackup_lacks_is_rejected(): + ports = [{"portnumber": 1, "source_layernum": 201, + "from_layername": "Metal5", "to_layername": "Metal9"}] + + with pytest.raises(ValueError, match="Metal9"): + port_contacts(ports, StackupXML.SG13G2_FEM_200um) + + +@pytest.mark.parametrize("geometry", [InductorOcta(), TransformerOcta()], ids=lambda g: g.name) +def test_sampled_preset_layouts_have_every_port_connected(tmp_path, geometry): + iterator = geometry.input_parameter_iterator + iterator.set_sample_count(30, seed=5, feasible=geometry.is_feasible) + + failures = {} + for i, params in enumerate(iterator): + path = str(tmp_path / f"{geometry.name}_{i}.gds") + geometry.create_gds_file(f"{geometry.name}_{i}", path, params) + ports = port_contacts(geometry.ports_for(params), geometry.stackup_xml) + result = check_gds_file(path, ports=ports) + if not result.ports_connected: + failures[str(params)] = result.port_findings + + assert not failures + + +def test_drc_stage_drops_open_ports_even_when_violations_are_kept(tmp_path): + context = PipelineContext(geometry=TransformerOcta(), base_dir=str(tmp_path), num_processes=1) + gds_dir = context.geometry_dir + tmp_path.joinpath(gds_dir).mkdir(parents=True, exist_ok=True) + good = {"bottom_winding_diameter": 21.4, "top_winding_diameter": 20.1, + "relative_displacement": 0.135, "bottom_linewidth": 7.0, "top_linewidth": 5.2} + rows = [] + for name, x in (("connected.gds", 0.0), ("open.gds", 0.105)): + if name == "connected.gds": + TransformerOcta.create_gds_file("connected", f"{gds_dir}/{name}", good) + else: + _layout(x).write(f"{gds_dir}/{name}") + rows.append({"name": name} | good) + pd.DataFrame(rows).to_csv(context.gds_csv_path, index=False) + + context = DRCChecker(drop_violations=False).run(context) + + passed = pd.read_csv(context.drc_csv_path) + assert list(passed["name"]) == ["connected.gds"] + report = pd.read_csv(context.drc_report_path).set_index("name") + assert "port1.open_TopMetal2:1" in report.loc["open.gds", "rules"] + + +def test_default_ports_are_the_simconfigs(): + geometry = TransformerOcta() + + ports = geometry.ports_for({}) + + assert ports == read_simconfig(geometry.simconfig_filename)["ports"] + ports[0]["to_layername"] = "Metal1" # a copy: changing it leaves the geometry alone + assert geometry.ports_for({})[0]["to_layername"] == "TopMetal2" + + +@pytest.mark.parametrize(("turns", "feed_layer"), [(1, "TopMetal2"), (1.0, "TopMetal2"), + (2, "TopMetal1"), (5, "TopMetal1")]) +def test_inductor_ports_follow_the_feed_layer(turns, feed_layer): + ports = InductorOcta().ports_for({"turns": turns, "width": 4.0, "space": 3.0, + "diameter": 100.0}) + + assert [p["to_layername"] for p in ports] == [feed_layer, feed_layer, "TopMetal2"] + assert {p["from_layername"] for p in ports} == {"Metal5"} + + +@dataclass +class _RenumberedPorts(TransformerOcta): + def simulation_ports(self, params): + return super().simulation_ports(params)[::-1] + + +@dataclass +class _DroppedPort(TransformerOcta): + def simulation_ports(self, params): + return super().simulation_ports(params)[:-1] + + +@pytest.mark.parametrize("geometry", [_RenumberedPorts(), _DroppedPort()], + ids=["renumbered", "dropped"]) +def test_ports_must_keep_the_simconfigs_numbers_and_order(geometry): + with pytest.raises(ValueError, match="simulation_ports"): + geometry.ports_for({}) + + +def test_drc_stage_checks_each_layout_with_its_own_ports(tmp_path): + geometry = InductorOcta() + context = PipelineContext(geometry=geometry, base_dir=str(tmp_path), num_processes=1) + gds_dir = tmp_path.joinpath(context.geometry_dir) + gds_dir.mkdir(parents=True, exist_ok=True) + rows = [] + for name, turns in (("one_turn.gds", 1), ("three_turns.gds", 3)): + params = {"turns": turns, "width": 4.0, "space": 3.0, "diameter": 160.0} + geometry.create_gds_file(name, str(gds_dir / name), params) + rows.append({"name": name} | params) + pd.DataFrame(rows).to_csv(context.gds_csv_path, index=False) + + context = DRCChecker().run(context) + + report = pd.read_csv(context.drc_report_path).set_index("name") + assert not report["rules"].fillna("").str.contains("port").any(), report["rules"] + + +@pytest.mark.parametrize( + "params", + [ + # the smallest single turns + {"turns": 1, "width": 2.02, "space": 2.32, "diameter": 32.0}, + {"turns": 1, "width": 2.1, "space": 2.98, "diameter": 34.0}, + # a larger multi-turn spiral + {"turns": 3, "width": 6.0, "space": 3.0, "diameter": 200.0}, + ], + ids=["small-D32", "small-D34", "large"], +) +def test_inductor_ports_sit_on_ground_at_every_size(tmp_path, params): + geometry = InductorOcta() + path = str(tmp_path / "inductor.gds") + geometry.create_gds_file("inductor", path, params) + + result = check_gds_file(path, ports=port_contacts(geometry.ports_for(params), geometry.stackup_xml)) + + assert result.port_findings == {} + + +@pytest.mark.parametrize("params", [{"turns": 1, "width": 2.46, "space": 2.36, "diameter": 32.0}, + {"turns": 2, "width": 5.0, "space": 3.0, "diameter": 120.0}, + {"turns": 3, "width": 8.0, "space": 4.1, "diameter": 260.0}], + ids=["one-turn", "two-turns", "three-turns"]) +def test_inductor_ports_sit_on_the_ring_outer_edge(tmp_path, params): + # The reference planes are the cell's boundary, so simulated cells can be abutted + # with their ports touching: each port is flush with the ring's outer edge + path = str(tmp_path / "inductor.gds") + InductorOcta().create_gds_file("inductor", path, params) + layout = kdb.Layout() + layout.read(path) + top = layout.top_cell() + ring = kdb.Region(top.begin_shapes_rec(layout.find_layer(*SG13G2.Metal5))).bbox() + + for marker_layer in (201, 202, 203): + marker = kdb.Region(top.begin_shapes_rec(layout.find_layer(marker_layer, 0))).bbox() + if marker.center().y < 0: + assert marker.bottom == ring.bottom + else: + assert marker.top == ring.top diff --git a/tests/test_losses.py b/tests/test_losses.py new file mode 100644 index 0000000..a9e8fc9 --- /dev/null +++ b/tests/test_losses.py @@ -0,0 +1,249 @@ +"""Loss terms, self-resonance detection and sample weighting (no Palace needed).""" + +import math +import os +from dataclasses import dataclass, field + +import numpy as np +import pandas as pd +import pytest + +torch = pytest.importorskip("torch") +skrf = pytest.importorskip("skrf") + +from orca.geometry.base_geometry import BaseGeometry +from orca.geometry.input_parameters import InputParameterIterator, RangeParameter +from orca.pipeline.context import PipelineContext +from orca.pipeline.training_stage import ModelTrainer +from orca.training.codecs import FlatReImCodec, UpperTriangleReImCodec +from orca.training.datasets.base_dataset import BaseDataset +from orca.training.datasets.geo_to_s_param_single_f import ( + GeoToSParamDatasetSingleFrequency, +) +from orca.training.losses import ( + SParameterLoss, + admittance_error, + first_self_resonance, + passivity_violation, + s_to_y, + srf_sample_weights, +) +from orca.training.normalize import MinMaxNormalizer, StandardNormalizer +from orca.training.trainer import make_batches + +# From DC, so the DC point that dc_deembedded results start with is covered too +FREQUENCIES = np.linspace(0.0, 10e9, 41) +L_NH = [1.0, 2.0, 3.0, 4.0] +C_PF = [0.5, 1.0, 1.5] +R_SERIES = 1.0 +Z0 = 50.0 + + +def _pi_network(inductance: float, capacitance: float, frequencies=FREQUENCIES) -> np.ndarray: + """S-matrices of a series R-L between two ports, with a shunt C at each port.""" + omega = 2 * np.pi * frequencies + y_series = 1 / (R_SERIES + 1j * omega * inductance) + y_shunt = 1j * omega * capacitance + y = np.empty((len(frequencies), 2, 2), dtype=complex) + y[:, 0, 0] = y[:, 1, 1] = y_series + y_shunt + y[:, 0, 1] = y[:, 1, 0] = -y_series + eye = np.eye(2) + return (eye - Z0 * y) @ np.linalg.inv(eye + Z0 * y) + + +def _resonance(inductance: float, capacitance: float) -> float: + """Where the differential susceptance Im(2 Y_series + Y_shunt) crosses zero.""" + return math.sqrt(2 / (inductance * capacitance) - (R_SERIES / inductance) ** 2) / (2 * math.pi) + + +def _iterator() -> InputParameterIterator: + return InputParameterIterator( + RangeParameter("a", L_NH[0], L_NH[-1], step=1.0), + RangeParameter("b", C_PF[0], C_PF[-1], step=0.5), + frequency=[FREQUENCIES[0], FREQUENCIES[-1]], + ) + + +def _iterator_normalizer() -> MinMaxNormalizer: + return MinMaxNormalizer(_iterator()) + + +def _dataset() -> GeoToSParamDatasetSingleFrequency: + return GeoToSParamDatasetSingleFrequency( + codec=FlatReImCodec(n_ports=2), + input_normalizer=_iterator_normalizer(), + output_normalizer=StandardNormalizer(), + ) + + +@dataclass +class ResonatorGeometry(BaseGeometry): + name: str = "resonator" + stackup_xml: str = "" + simconfig_filename: str = "" + input_parameter_iterator: InputParameterIterator = field(default_factory=_iterator) + + @staticmethod + def create_gds_file(name: str, output_path: str, params: dict) -> str: + raise NotImplementedError + + def create_dataset(self) -> GeoToSParamDatasetSingleFrequency: + return _dataset() + + +@pytest.fixture +def result_dir(tmp_path) -> str: + """One pi-network per (L, C) pair; the smallest L*C resonates just above the band.""" + rows = [] + for i, (a, b) in enumerate((a, b) for a in L_NH for b in C_PF): + s = _pi_network(a * 1e-9, b * 1e-12) + ntwk = skrf.Network(frequency=skrf.Frequency.from_f(FREQUENCIES, unit="hz"), s=s) + ntwk.write_touchstone(filename=f"res_{i}", dir=str(tmp_path)) + rows.append({"name": f"res_{i}.s2p", "a": a, "b": b}) + pd.DataFrame(rows).to_csv(tmp_path / "resonator.csv", index=False) + return str(tmp_path) + + +def _loaded(result_dir: str) -> BaseDataset: + table = pd.read_csv(os.path.join(result_dir, "resonator.csv")) + return _dataset().new_split(result_dir, table, fit_normalizers=True) + + +@pytest.mark.parametrize("codec", [FlatReImCodec(n_ports=3), UpperTriangleReImCodec(n_ports=3)]) +def test_tensor_decode_matches_numpy_decode(codec): + raw = np.random.default_rng(0).normal(size=(5, codec.output_dim)).astype(np.float32) + + decoded = codec.decode_tensor(torch.from_numpy(raw)).numpy() + + np.testing.assert_array_equal(decoded, codec.decode(raw)) + + +def test_admittance_matches_scikit_rf(): + s = _pi_network(2e-9, 1e-12) + ntwk = skrf.Network(frequency=skrf.Frequency.from_f(FREQUENCIES, unit="hz"), s=s, z0=Z0) + + y = s_to_y(torch.from_numpy(s)).numpy() + + np.testing.assert_allclose(y, ntwk.y * Z0, rtol=1e-9, atol=1e-12) + + +@pytest.mark.parametrize(("inductance", "capacitance"), [(1e-9, 1e-12), (4e-9, 1.5e-12)]) +def test_self_resonance_of_a_pi_network(inductance, capacitance): + s = _pi_network(inductance, capacitance) + shuffled = np.random.default_rng(1).permutation(len(FREQUENCIES)) + + found = first_self_resonance(FREQUENCIES[shuffled], s[shuffled]) + + assert found == pytest.approx(_resonance(inductance, capacitance), rel=1e-2) + + +def test_no_self_resonance_in_the_band_is_infinite(): + # 1 nH and 0.5 pF resonate at about 10.07 GHz, just above the band + assert first_self_resonance(FREQUENCIES, _pi_network(1e-9, 0.5e-12)) == math.inf + + +def test_points_above_the_resonance_are_weighted_apart(result_dir): + dataset = _loaded(result_dir) + + weights = srf_sample_weights(dataset, above_srf_weight=0.1).cpu().numpy() + + assert weights.mean() == pytest.approx(1.0, rel=1e-6) + frequencies = _iterator_normalizer().denormalize(dataset.inputs.cpu())[:, -1].numpy() + table = pd.read_csv(os.path.join(result_dir, "resonator.csv")).set_index("name") + groups = np.asarray(dataset.sample_groups) + below_weight = weights[frequencies <= 1e9].max() + for name, (a, b) in table[["a", "b"]].iterrows(): + mine = groups == name + above = frequencies[mine] > _resonance(a * 1e-9, b * 1e-12) + expected = np.where(above, 0.1 * below_weight, below_weight) + np.testing.assert_allclose(weights[mine], expected, rtol=1e-5) + # 1 nH, 0.5 pF has no resonance in the band and keeps full weight everywhere + assert np.all(weights[groups == "res_0.s2p"] == below_weight) + + +def test_without_terms_the_loss_is_the_data_loss(): + pred, target = torch.randn(6, 8), torch.randn(6, 8) + loss = SParameterLoss(torch.nn.L1Loss(), FlatReImCodec(n_ports=2)) + + assert loss(pred, target) == pytest.approx(torch.nn.functional.l1_loss(pred, target).item()) + + weights = torch.tensor([3.0, 0.0, 0.0, 1.0, 1.0, 1.0]) + per_sample = (pred - target).abs().mean(dim=1) + assert loss(pred, target, weights) == pytest.approx((per_sample * weights).mean().item()) + + +def test_loss_needs_an_elementwise_data_loss(): + with pytest.raises(TypeError, match="reduction"): + SParameterLoss(lambda pred, target: (pred - target).abs().mean(), FlatReImCodec(2)) + + +def test_admittance_error_sees_a_loss_error_that_s_hides(): + s_true = torch.from_numpy(_pi_network(2e-9, 1e-12)[1:5]) # 0.25 to 1 GHz + # Twice the series resistance: Q halves, but S barely moves + eye = np.eye(2) + omega = 2 * np.pi * FREQUENCIES[1:5] + y_series = 1 / (2 * R_SERIES + 1j * omega * 2e-9) + y = np.empty((4, 2, 2), dtype=complex) + y[:, 0, 0] = y[:, 1, 1] = y_series + 1j * omega * 1e-12 + y[:, 0, 1] = y[:, 1, 0] = -y_series + s_lossier = torch.from_numpy((eye - Z0 * y) @ np.linalg.inv(eye + Z0 * y)) + + assert torch.all(admittance_error(s_true, s_true) == 0) + assert torch.all((s_lossier - s_true).abs().amax(dim=(1, 2)) < 0.05) + assert torch.all(admittance_error(s_lossier, s_true) > 0.5) + + +def test_passivity_is_only_penalised_beyond_the_reference(): + eye = torch.eye(2, dtype=torch.complex64).expand(3, 2, 2) + s_true = eye * torch.tensor([0.5, 1.0, 1.1]).view(3, 1, 1) + s_pred = eye * torch.tensor([0.9, 1.2, 1.1]).view(3, 1, 1) + + violation = passivity_violation(s_pred, s_true) + + np.testing.assert_allclose(violation.numpy(), [0.0, 0.2, 0.0], atol=1e-6) + + +def test_sample_weights_are_batched_with_their_samples(result_dir): + dataset = _loaded(result_dir) + dataset.set_sample_weights(torch.arange(len(dataset), dtype=torch.float32, device=dataset.device)) + subset = torch.utils.data.Subset(dataset, [5, 2, 7]) + + (x, _, w), *_ = list(make_batches(subset, batch_size=10, shuffle=False)) + + assert torch.equal(x, dataset.inputs[[5, 2, 7]]) + assert w.tolist() == [5.0, 2.0, 7.0] + with pytest.raises(ValueError, match="one weight per sample"): + dataset.set_sample_weights(torch.ones(3)) + + +def test_stage_trains_with_every_loss_option(result_dir, tmp_path): + context = PipelineContext( + geometry=ResonatorGeometry(), + base_dir=str(tmp_path / "run"), + result_dir_override=result_dir, + seed=3, + ) + trainer = ModelTrainer( + n_fold_cv=1, + n_trials=1, + tuning_max_epochs=2, + max_epochs=2, + batch_sizes=[16], + admittance_weight=0.5, + passivity_weight=1.0, + above_srf_weight=0.1, + ) + + context = trainer.run(context) + + assert context.final_val_loss is not None + assert math.isfinite(context.final_val_loss) + + +@pytest.mark.parametrize( + "options", + [{"admittance_weight": -1.0}, {"passivity_weight": -0.1}, {"above_srf_weight": 0.0}], +) +def test_invalid_loss_options_fail_when_the_stage_is_built(options): + with pytest.raises(ValueError, match="must"): + ModelTrainer(**options) diff --git a/tests/test_mesh_quality.py b/tests/test_mesh_quality.py new file mode 100644 index 0000000..fff2db6 --- /dev/null +++ b/tests/test_mesh_quality.py @@ -0,0 +1,59 @@ +"""Tests for the mesh quality check of the GDS conversion stage.""" + +from __future__ import annotations + +import math +import os + +import pandas as pd + +import orca.pipeline.gds_conversion_stage as conversion_stage +from orca.geometry.presets import InductorOcta +from orca.pipeline.context import PipelineContext +from orca.pipeline.gds_conversion_stage import GDSConverter + + +def _tetrahedron_mesh(path, apex_height: float) -> str: + """A one-tetrahedron mesh in gmsh's MSH 2.2 format; height 0 makes it flat.""" + with open(path, "w") as f: + f.write( + "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n" + "$Nodes\n4\n1 0 0 0\n2 1 0 0\n3 0 1 0\n" + f"4 0.25 0.25 {apex_height}\n$EndNodes\n" + "$Elements\n1\n1 4 2 1 1 1 2 3 4\n$EndElements\n" + ) + return str(path) + + +QUALITY = {"good.gds": 0.01, "poor.gds": 5e-5, "flat.gds": -1.7e-14, "corrupt.gds": float("-inf")} + + +def _fake_conversion(geometry_name, params, output_dir, gds_filename, stackup_xml, + simconfig_filename, show_mesh_results=False, ports=None): + sim_path = os.path.join(output_dir, "palace_sims", f"{geometry_name[:-4]}_data") + os.makedirs(sim_path, exist_ok=True) + open(os.path.join(sim_path, "config.json"), "w").close() + return geometry_name, params, "config.json", sim_path, sim_path, QUALITY[geometry_name] + + +def test_flat_meshes_are_left_out_and_reported(tmp_path, monkeypatch, caplog): + monkeypatch.setattr(conversion_stage, "create_palace_model_from_gds", _fake_conversion) + context = PipelineContext(geometry=InductorOcta(), base_dir=str(tmp_path), num_processes=1) + os.makedirs(context.geometry_dir, exist_ok=True) + params = {"turns": 2, "width": 4.0, "space": 3.0, "diameter": 120.0} + pd.DataFrame([{"name": name} | params for name in QUALITY]).to_csv( + context.gds_csv_path, index=False + ) + + context = GDSConverter().run(context) + + assert sorted(pd.read_csv(context.palace_csv_path)["name"]) == ["good.gds", "poor.gds"] + report = pd.read_csv(context.mesh_report_path).set_index("name")["worst_element_quality"] + assert all( + report[name] == q or math.isclose(report[name], q, abs_tol=1e-15) + for name, q in QUALITY.items() + ) + messages = " ".join(r.getMessage() for r in caplog.records) + assert "2 of 4 meshes contain flat, inverted or corrupt elements" in messages + assert "refined_cellsize = 2 µm" in messages + assert "1 of 4 meshes have poorly shaped elements" in messages diff --git a/tests/test_postprocessing.py b/tests/test_postprocessing.py new file mode 100644 index 0000000..fdda52a --- /dev/null +++ b/tests/test_postprocessing.py @@ -0,0 +1,125 @@ +"""Electrical parameters of the presets' devices, on circuits with known values.""" + +import math + +import numpy as np +import pytest + +skrf = pytest.importorskip("skrf") + +from orca.geometry.presets import InductorOcta, TransformerOcta +from orca.utils.postprocessing import _first_inductive_to_capacitive, differential_impedance + +FREQUENCIES = np.linspace(0.0, 20e9, 81) +L_DIFF = 2e-9 +C_SHUNT = 0.5e-12 +R_DIFF = 2.0 +Z0 = 50.0 + + +def _network(y: np.ndarray, frequencies: np.ndarray) -> "skrf.Network": + n = y.shape[-1] + eye = np.eye(n) + s = (eye - Z0 * y) @ np.linalg.inv(eye + Z0 * y) + return skrf.Network(frequency=skrf.Frequency.from_f(frequencies, unit="hz"), s=s, z0=Z0) + + +def _center_tapped_inductor(frequencies=FREQUENCIES) -> "skrf.Network": + """Two half windings from ports 1 and 2 to the center tap (port 3), with a shunt C at the ends.""" + omega = 2 * np.pi * frequencies + y_half = 1 / (R_DIFF / 2 + 1j * omega * L_DIFF / 2) + y_shunt = 1j * omega * C_SHUNT + y = np.zeros((len(frequencies), 3, 3), dtype=complex) + y[:, 0, 0] = y[:, 1, 1] = y_half + y_shunt + y[:, 0, 2] = y[:, 2, 0] = y[:, 1, 2] = y[:, 2, 1] = -y_half + y[:, 2, 2] = 2 * y_half + return _network(y, frequencies) + + +#: Branch inductances [nH] of four coupled half windings: top op->oci, oci->on and bottom +#: ip->ico, ico->in, so a current through either winding runs forward through both halves +HALF_WINDINGS = np.array( + [ + [0.50, 0.10, 0.15, 0.15], + [0.10, 0.50, 0.15, 0.15], + [0.15, 0.15, 0.50, 0.10], + [0.15, 0.15, 0.10, 0.50], + ] +) +L_WINDING = 2 * 0.5 + 2 * 0.1 # nH, both halves and their mutual inductance +M_WINDINGS = 4 * 0.15 # nH, every half of one winding with every half of the other + + +def _transformer(frequency: float = 1e9) -> "skrf.Network": + """TransformerOcta's six ports: 1 op, 2 on, 3 ip, 4 in, 5 oci, 6 ico.""" + omega = 2 * np.pi * frequency + z_branch = 0.5 * np.eye(4) + 1j * omega * HALF_WINDINGS * 1e-9 + # Node-branch incidence: each branch leaves one port and enters another + incidence = np.zeros((6, 4)) + for branch, (start, end) in enumerate([(0, 4), (4, 1), (2, 5), (5, 3)]): + incidence[start, branch], incidence[end, branch] = 1, -1 + y = incidence @ np.linalg.inv(z_branch) @ incidence.T + 1j * omega * 5e-15 * np.eye(6) + return _network(y[None], np.array([frequency])) + + +def test_inductor_preset_reads_the_center_tapped_inductor(): + metrics = InductorOcta().electrical_parameters(_center_tapped_inductor()) + + assert set(metrics) == {"L", "R", "Q", "srf_f"} + low = 1 # 0.25 GHz, well below resonance + omega = 2 * np.pi * FREQUENCIES[low] + assert metrics["L"][low] == pytest.approx(L_DIFF * 1e9, rel=1e-2) + assert metrics["R"][low] == pytest.approx(R_DIFF, rel=1e-2) + assert metrics["Q"][low] == pytest.approx(omega * L_DIFF / R_DIFF, rel=2e-2) + # Each half winding resonates with its shunt C: w^2 = 2 / (L C) - (R / L)^2 + resonance = math.sqrt(2 / (L_DIFF * C_SHUNT) - (R_DIFF / L_DIFF) ** 2) / (2 * math.pi) + assert float(metrics["srf_f"]) == pytest.approx(resonance / 1e9, rel=1e-2) + + +def test_inductor_below_resonance_has_no_srf(): + metrics = InductorOcta().electrical_parameters(_center_tapped_inductor(FREQUENCIES[:20])) + + assert np.isnan(metrics["srf_f"]) + + +def test_transformer_preset_pairs_the_ports_of_each_winding(): + metrics = TransformerOcta().electrical_parameters(_transformer()) + + assert metrics["Lp"][0] == pytest.approx(L_WINDING, rel=1e-2) + assert metrics["Ls"][0] == pytest.approx(L_WINDING, rel=1e-2) + assert metrics["Rp"][0] == pytest.approx(1.0, rel=1e-2) + assert metrics["k"][0] == pytest.approx(M_WINDINGS / L_WINDING, rel=1e-2) + assert TransformerOcta.absolute_error_parameters == {"k"} + + +def test_differential_impedance_of_a_floating_source_across_a_series_element(): + # A lone series impedance between two ports: the floating source sees exactly it + z_series = 3 + 4j + y = np.array([[1, -1], [-1, 1]]) / z_series + ntwk = _network(y[None], np.array([1e9])) + + assert differential_impedance(ntwk, [(0, 1)])[0, 0, 0] == pytest.approx(z_series) + + +def test_tester_reports_the_geometry_errors(): + from orca.pipeline.test_model_stage import ModelTester + + reference = _center_tapped_inductor() + predicted = reference.copy() + predicted.s = predicted.s * 0.99 + + medians, curves = ModelTester._electrical_errors( + predicted, reference, "inductor", InductorOcta() + ) + + assert set(medians) == {"L error %", "R error %", "Q error %", "srf_f error %"} + assert set(curves) == {"L", "R", "Q"} + assert all(np.isfinite(error) and error > 0 for error in medians.values()) + + +def test_srf_ignores_the_sign_of_the_reactance_at_dc(): + # An open winding is capacitive from DC on; rounding can leave +0 reactance at DC + frequencies = np.array([0.0, 1e9, 2e9]) + reactance = np.array([1e-12, -500.0, -250.0]) + + assert np.isnan(_first_inductive_to_capacitive(frequencies / 1e9, reactance)) diff --git a/tests/test_training.py b/tests/test_training.py index bd73c09..62d64ab 100644 --- a/tests/test_training.py +++ b/tests/test_training.py @@ -6,6 +6,7 @@ import os import random from dataclasses import dataclass, field +from typing import ClassVar import numpy as np import pandas as pd @@ -67,6 +68,11 @@ class ToyGeometry(BaseGeometry): simconfig_filename: str = "" input_parameter_iterator: InputParameterIterator = field(default_factory=_iterator) + absolute_error_parameters: ClassVar[frozenset[str]] = frozenset({"|S21|"}) + + def electrical_parameters(self, ntwk) -> dict[str, np.ndarray]: + return {"|S11|": np.abs(ntwk.s[:, 0, 0]), "|S21|": np.abs(ntwk.s[:, 1, 0])} + @staticmethod def create_gds_file(name: str, output_path: str, params: dict) -> str: raise NotImplementedError @@ -297,12 +303,14 @@ def test_tester_evaluates_the_split_the_trainer_recorded(result_dir, tmp_path): per_geometry = pd.read_csv(later.test_errors_csv_path) assert set(per_geometry["name"]) == test_names assert {"a", "b", "mean_abs_s_error", "max_abs_s_error"} <= set(per_geometry.columns) - # k near zero would make a relative error meaningless, so it is reported as absolute - assert "k abs error" in per_geometry.columns - assert "k error %" not in per_geometry.columns + # The geometry's electrical parameters, absolute where it says so + assert "|S11| error %" in per_geometry.columns + assert "|S21| abs error" in per_geometry.columns + assert "|S21| error %" not in per_geometry.columns profile = pd.read_csv(later.errors_vs_frequency_csv_path) - assert {"|S|", "Lp", "k"} <= set(profile["parameter"]) + assert {"|S|", "|S11|", "|S21|"} <= set(profile["parameter"]) + assert set(profile.loc[profile["parameter"] == "|S21|", "unit"]) == {"abs"} s_rows = profile[profile["parameter"] == "|S|"] assert len(s_rows) == len(FREQUENCIES) assert (s_rows["p5"] <= s_rows["median"]).all() diff --git a/tests/test_transformer.py b/tests/test_transformer.py index a2a20eb..11555a0 100644 --- a/tests/test_transformer.py +++ b/tests/test_transformer.py @@ -8,7 +8,7 @@ import pytest from orca.geometry.cells.transformer import check_tf_octa_c_parameters, tf_octa_c -from orca.geometry.drc import check_gds_file +from orca.geometry.drc import check_gds_file, port_contacts from orca.geometry.layers import SG13G2 from orca.geometry.presets import TransformerOcta @@ -80,6 +80,29 @@ def test_each_winding_with_its_feeds_and_tap_is_one_polygon(tmp_path): assert _metal(path, layer).count() == 1 +@pytest.mark.parametrize( + "params", + [ + # Draws whose port markers were rounded 5 nm past the feed ends: all six ports, + # the right-hand ones (op, on, ico) and the left-hand ones (ip, in, oci) open + {"bottom_winding_diameter": 28.5, "top_winding_diameter": 23.6, + "relative_displacement": 0.17, "bottom_linewidth": 7.8, "top_linewidth": 6.2}, + {"bottom_winding_diameter": 28.6, "top_winding_diameter": 24.7, + "relative_displacement": 0.02, "bottom_linewidth": 7.6, "top_linewidth": 5.2}, + {"bottom_winding_diameter": 29.2, "top_winding_diameter": 28.8, + "relative_displacement": 0.07, "bottom_linewidth": 3.5, "top_linewidth": 3.1}, + ], + ids=["all-open", "right-open", "left-open"], +) +def test_port_markers_touch_the_feed_ends(tmp_path, params): + geometry = TransformerOcta() + ports = port_contacts(geometry.ports_for(params), geometry.stackup_xml) + + result = check_gds_file(_draw(tmp_path, params), ports=ports) + + assert result.port_findings == {} + + def test_feed_gap_wider_than_the_flat_side_is_rejected(): # A 20 µm octagon's flat side is 20 * sin(22.5 deg) = 7.65 µm along the centre line check_tf_octa_c_parameters(bottom_winding_diameter=20.0, top_winding_diameter=20.0, @@ -107,9 +130,30 @@ def test_gap_beyond_the_inner_flat_side_draws_without_folding(tmp_path): assert check_gds_file(path).clean -def test_ports_inside_the_windings_are_rejected(): - with pytest.raises(ValueError, match="inside the windings"): - check_tf_octa_c_parameters(top_linewidth=6.0, gnd_upper_spacing=12.0, gnd_ring_width=10.0) +def test_ring_without_clearance_is_rejected(): + for no_clearance in ( + lambda: check_tf_octa_c_parameters(gnd_upper_spacing=0.0), + lambda: check_tf_octa_c_parameters(gnd_lower_spacing=0.0), + lambda: check_tf_octa_c_parameters(gnd_side_spacing=0.0), + ): + with pytest.raises(ValueError, match="under the windings"): + no_clearance() + + +def test_ring_keeps_its_clearance_from_the_winding_metal(tmp_path): + # The ring's inner edges keep GROUND_SPACING from the windings' outermost metal, as + # the inductor's ring does from its spiral + path = _draw(tmp_path, SMALLEST | {"relative_displacement": 0.2, "bottom_linewidth": 6.0}) + windings = (_metal(path, SG13G2.TopMetal1) + _metal(path, SG13G2.TopMetal2)).merged() + ring = _metal(path, SG13G2.Metal5) + opening = (kdb.Region(ring.bbox()) - ring).merged() + + # the feeds and center taps cross the ring; the octagons alone set the clearance + windings_body = windings - (kdb.Region(ring.bbox()) - opening) + spacing = 20.0 / 0.001 # GROUND_SPACING in database units (1 nm) + body, inner = windings_body.bbox(), opening.bbox() + assert inner.top - body.top == pytest.approx(spacing, abs=5) + assert body.bottom - inner.bottom == pytest.approx(spacing, abs=5) @pytest.mark.parametrize( @@ -151,3 +195,25 @@ def test_offset_scales_with_the_windings_and_lands_on_the_grid(tmp_path): (text,) = [s.text.string for s in layout.top_cell().shapes(layout.find_layer(*SG13G2.TEXT)).each()] assert "center displacement: 1.240" in text # each winding centre on the 5 nm grid assert check_gds_file(path).snapped_vertices == 0 + + +@pytest.mark.parametrize( + "params", + [SMALLEST, SMALLEST | {"bottom_winding_diameter": 90.0, "top_winding_diameter": 80.0, + "relative_displacement": 0.2, "bottom_linewidth": 8.0}], + ids=["smallest", "large-offset"], +) +def test_ports_sit_on_the_ring_outer_edge(tmp_path, params): + # The reference planes are the cell's boundary, so cells can be abutted with touching + # ports, as the inductor's + path = _draw(tmp_path, params) + layout = kdb.Layout() + layout.read(path) + top = layout.top_cell() + ring = kdb.Region(top.begin_shapes_rec(layout.find_layer(*SG13G2.Metal5))).bbox() + + for marker_layer in range(201, 207): + # each marker is a zero-width path across its feed, at constant x + (marker,) = top.shapes(layout.find_layer(marker_layer, 0)).each() + xs = {point.x for point in marker.path.each_point()} + assert xs in ({ring.left}, {ring.right}), marker_layer