Source code for superscreen.device.device

import logging
import numbers
import os
from collections import defaultdict
from contextlib import contextmanager, nullcontext
from typing import Dict, List, Literal, Optional, Sequence, Tuple, Union

import dill
import h5py
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import HTML
from matplotlib.patches import PathPatch
from matplotlib.path import Path
from shapely import geometry as geo
from tqdm import tqdm

from .. import fem
from ..geometry import ensure_unique
from ..units import ureg
from . import utils
from .layer import Layer
from .mesh import Mesh
from .polygon import Polygon

logger = logging.getLogger("device")


[docs]class Device: """An object representing a device composed of one or more layers of thin film superconductor. Args: name: Name of the device. layers: ``Layers`` making up the device. films: ``Polygons`` representing regions of superconductor. holes: ``Holes`` representing holes in superconducting films. terminals: A dict of ``{film_name: terminals}`` representing transport terminals in the device. abstract_regions: ``Polygons`` representing abstract regions in a device. length_units: Distance units for the coordinate system. solve_dtype: The float data type to use when solving the device. """ ureg = ureg def __init__( self, name: str, *, layers: Union[Sequence[Layer], Dict[str, Layer]], films: Union[Sequence[Polygon], Dict[str, Polygon]], holes: Optional[Union[Sequence[Polygon], Dict[str, Polygon]]] = None, terminals: Optional[Dict[str, List[Polygon]]] = None, abstract_regions: Optional[Union[Sequence[Polygon], Dict[str, Polygon]]] = None, length_units: str = "um", solve_dtype: Union[str, np.dtype] = "float32", ): self.name = name if isinstance(layers, dict): layers = list(layers.values()) self.layers = {layer.name: layer for layer in layers} if isinstance(films, dict): films = list(films.values()) self.films = {film.name: film for film in films} if holes is None: holes = [] if isinstance(holes, dict): holes = list(holes.values()) self.holes = {hole.name: hole for hole in holes} if terminals is None: terminals = {} self.terminals = terminals if not set(self.terminals).issubset(self.films): raise ValueError( f"terminals.keys() must be a subset of films.keys() ({list(self.films)!r})." ) for film, terminals in self.terminals.items(): for terminal in terminals: terminal.layer = self.films[film].layer if abstract_regions is None: abstract_regions = [] if isinstance(abstract_regions, dict): abstract_regions = list(abstract_regions.values()) self.abstract_regions = {region.name: region for region in abstract_regions} for polygons, label in [ (self.films.values(), "film"), (self.holes.values(), "hole"), ]: for polygon in polygons: if not polygon.is_valid: raise ValueError(f"The following {label} is not valid: {polygon}.") if polygon.layer not in self.layers: raise ValueError( f"The following {label} is assigned to a layer that doesn not " f"exist in the device: {polygon}." ) # Make units a "read-only" attribute. # It should never be changed after instantiation. self._length_units = length_units self.solve_dtype = solve_dtype self.meshes: Union[Dict[str, Mesh], None] = None @property def length_units(self) -> str: """Length units used for the device geometry.""" return self._length_units @property def solve_dtype(self) -> np.dtype: """Numpy dtype to use for floating point numbers.""" return self._solve_dtype @solve_dtype.setter def solve_dtype(self, dtype) -> None: try: _ = np.finfo(dtype) except ValueError as e: raise ValueError(f"Invalid float dtype: {dtype}") from e self._solve_dtype = np.dtype(dtype)
[docs] def get_polygons(self, include_terminals: bool = True) -> List[Polygon]: """Returns list of all Polygons in the device. Args: include_terminals: Include transport terminals in the list. Returns: A list of all Polygons in the device. """ polygons = [] for attr_name in ("films", "holes", "abstract_regions"): polygons.extend(list(getattr(self, attr_name).values())) if include_terminals: for terminals in self.terminals.values(): polygons.extend(terminals) return polygons
@property def poly_points(self) -> np.ndarray: """Shape (n, 2) array of (x, y) coordinates of all Polygons in the Device.""" points = np.concatenate( [poly.points for poly in self.get_polygons(include_terminals=False)] ) # Remove duplicate points to avoid meshing issues. # If you don't do this and there are duplicate points, # meshpy.triangle will segfault. return ensure_unique(points)
[docs] def polygons_by_layer( self, polygon_type: Optional[ Literal["film", "hole", "abstract", "terminal", "all"] ] = None, ) -> Dict[str, List[Polygon]]: """Returns a dict of ``{layer_name: list of polygons in layer}``. Args: polygon_type: One of 'film', 'hole', 'abstract', 'terminal', or 'all', specifying which types of polygons to include. Returns: A dict of ``{layer_name: list of polygons in layer of the given type}``. """ valid_types = ("film", "hole", "abstract", "terminal", "all") if polygon_type is None: polygon_type = "all" polygon_type = polygon_type.lower() if polygon_type not in valid_types: raise ValueError( f"Invalid polygon type ({polygon_type}). Expected one of {valid_types!r}." ) if polygon_type == "film": all_polygons = self.films.values() elif polygon_type == "hole": all_polygons = self.holes.values() elif polygon_type == "abstract": all_polygons = self.abstract_regions.values() elif polygon_type == "terminal": all_polygons = [] for terminals in self.terminals.values(): all_polygons.extend(terminals) else: # Films + holes + terminals + abstract regions all_polygons = self.get_polygons() polygons = {} for layer in self.layers: polygons[layer] = [p for p in all_polygons if p.layer == layer] return polygons
[docs] def holes_by_film(self) -> Dict[str, List[Polygon]]: """Generates a mapping of films to holes contained in the film. Returns: A dict of ``{film_name: list of holes in the film}``. """ holes_by_layer = self.polygons_by_layer("hole") holes_by_film = {} for film in self.films.values(): holes_by_film[film.name] = [] for hole in holes_by_layer[film.layer]: if film.contains_points(hole.points).all(): holes_by_film[film.name].append(hole) return holes_by_film
[docs] def copy(self, with_mesh: bool = True, copy_mesh: bool = False) -> "Device": """Copy this Device to create a new one. Args: with_mesh: Whether to shallow copy the ``meshes`` dictionary. copy_mesh: Whether to deepcopy the arrays defining the mesh. Returns: A new Device instance, copied from self """ layers = [layer.copy() for layer in self.layers.values()] films = [film.copy() for film in self.films.values()] holes = [hole.copy() for hole in self.holes.values()] terminals = { film: [term.copy() for term in film_terms] for film, film_terms in self.terminals.items() } abstract_regions = [region.copy() for region in self.abstract_regions.values()] device = Device( self.name, layers=layers, films=films, holes=holes, terminals=terminals, abstract_regions=abstract_regions, length_units=self.length_units, ) if with_mesh and self.meshes is not None: meshes = self.meshes if copy_mesh: meshes = {name: mesh.copy() for name, mesh in meshes.items()} device.meshes = meshes return device
def __copy__(self) -> "Device": # Shallow copy (create new references to existing arrays). return self.copy(with_mesh=True, copy_mesh=False) def __deepcopy__(self, memo) -> "Device": # Deep copy (copy all arrays and return references to copies) return self.copy(with_mesh=True, copy_mesh=True) def _warn_if_mesh_exist(self, method: str) -> None: if not self.meshes: return message = ( f"Calling device.{method} on a device whose mesh already exists returns " f"a new device with no mesh. Call new_device.make_mesh() to generate the mesh " f"for the new device." ) logger.warning(message)
[docs] def scale( self, xfact: float = 1, yfact: float = 1, origin: Tuple[float, float] = (0, 0) ) -> "Device": """Returns a new device with polygons scaled horizontally and/or vertically. Negative ``xfact`` (``yfact``) can be used to reflect the device horizontally (vertically) about the ``origin``. Args: xfact: Factor by which to scale the device horizontally. yfact: Factor by which to scale the device vertically. origin: (x, y) coorindates of the origin. Returns: The scaled :class:`superscreen.Device`. """ if not ( isinstance(origin, tuple) and len(origin) == 2 and all(isinstance(val, numbers.Real) for val in origin) ): raise TypeError("Origin must be a tuple of floats (x, y).") self._warn_if_mesh_exist("scale()") device = self.copy(with_mesh=False) for polygon in device.get_polygons(): polygon.scale(xfact=xfact, yfact=yfact, origin=origin, inplace=True) return device
[docs] def rotate(self, degrees: float, origin: Tuple[float, float] = (0, 0)) -> "Device": """Returns a new device with polygons rotated a given amount counterclockwise about specified origin. Args: degrees: The amount by which to rotate the polygons. origin: (x, y) coorindates of the origin. Returns: The rotated :class:`superscreen.Device`. """ if not ( isinstance(origin, tuple) and len(origin) == 2 and all(isinstance(val, numbers.Real) for val in origin) ): raise TypeError("Origin must be a tuple of floats (x, y).") self._warn_if_mesh_exist("rotate()") device = self.copy(with_mesh=False) for polygon in device.get_polygons(): polygon.rotate(degrees, origin=origin, inplace=True) return device
[docs] def mirror_layers(self, about_z: float = 0.0) -> "Device": """Returns a new device with its layers mirrored about the plane ``z = about_z``. Args: about_z: The z-position of the plane (parallel to the x-y plane) about which to mirror the layers. Returns: The mirrored :class:`superscreen.Device`. """ self._warn_if_mesh_exist("mirror_layers()") device = self.copy(with_mesh=False) for layer in device.layers.values(): layer.z0 = about_z - layer.z0 return device
[docs] def translate( self, dx: float = 0, dy: float = 0, dz: float = 0, inplace: bool = False, ) -> "Device": """Translates the device polygons, layers, and meshes in space by a given amount. Args: dx: Distance by which to translate along the x-axis. dy: Distance by which to translate along the y-axis. dz: Distance by which to translate layers along the z-axis. inplace: If True, modifies the device (``self``) in-place and returns None, otherwise, creates a new device, translates it, and returns it. Returns: The translated device. """ if inplace: device = self else: device = self.copy(with_mesh=True, copy_mesh=True) for polygon in device.get_polygons(): polygon.translate(dx, dy, inplace=True) if device.meshes: for mesh in device.meshes.values(): mesh.sites += np.array([[dx, dy]]) if dz: for layer in device.layers.values(): layer.z0 += dz return device
[docs] @contextmanager def translation(self, dx: float, dy: float, dz: float = 0): """A context manager that temporarily translates a device in-place, then returns it to its original position. Args: dx: Distance by which to translate polygons along the x-axis. dy: Distance by which to translate polygons along the y-axis. dz: Distance by which to translate layers along the z-axis. """ try: self.translate(dx, dy, dz=dz, inplace=True) yield finally: self.translate(-dx, -dy, dz=-dz, inplace=True)
[docs] def make_mesh( self, buffer_factor: Union[float, Dict[str, float], None] = 0.05, buffer: Union[float, Dict[str, float], None] = None, join_style: str = "round", min_points: Union[int, Dict[str, int], None] = None, max_edge_length: Union[float, Dict[str, float], None] = None, preserve_boundary: bool = False, smooth: Union[int, Dict[str, int]] = 0, **meshpy_kwargs, ) -> None: """Generates the triangular mesh for each film and stores them in the ``self.meshes`` dictionary. The arguments ``buffer_factor``, ``buffer``, ``min_points``, ``max_edge_length``, and ``smooth`` can be specified either as a single value for all films or as a dict of ``{film_name: argument_value}``. Args: buffer_factor: Buffer for the film bounding box(es), in units of the maximum film dimension. This argument is ignored if ``buffer`` is not None. buffer: Buffer for the film bounding box(es), in ``length_units``. join_style: The join style for the buffered region (see :meth:`superscreen.Polygon.buffer`). min_points: Minimum number of vertices in the mesh. If None, then the number of vertices will be determined by meshpy_kwargs and the number of vertices in the underlying polygons. max_edge_length: The maximum distance between vertices in the resulting mesh. preserve_boundary: Do not add any mesh sites to the boundary. smooth: Number of Laplacian smoothing iterations to perform. **meshpy_kwargs: Passed to meshpy.triangle.build(). """ films = self.films meshes = {} if not isinstance(buffer_factor, dict): buffer_factor = {name: buffer_factor for name in films} if not isinstance(buffer, dict): buffer = {name: buffer for name in films} if not isinstance(min_points, dict): min_points = {name: min_points for name in films} if not isinstance(max_edge_length, dict): max_edge_length = {name: max_edge_length for name in films} if not isinstance(smooth, dict): smooth = {name: smooth for name in films} holes_by_layer = self.polygons_by_layer("hole") abs_regions_by_layer = self.polygons_by_layer("abstract") for name, film in films.items(): film_terminals = self.terminals.get(name) holes = holes_by_layer[film.layer] abs_regions = abs_regions_by_layer[film.layer] coords = [film.points] for poly in holes + abs_regions: if film.contains_points(poly.points).all(): coords.append(poly.points) if ( film_terminals is not None or buffer[name] == 0 or (buffer_factor[name] is None and buffer[name] is None) ): boundary = film.points else: boundary = Polygon(points=film.points) if buffer[name] is None: buffer_size = buffer_factor[name] * max(boundary.extents) else: buffer_size = buffer[name] buffered = boundary.polygon.buffer( buffer_size, join_style=getattr(geo.JOIN_STYLE, join_style), single_sided=True, mitre_limit=5.0, ).exterior.coords boundary = Polygon(points=buffered).resample(len(film.points)).points coords.append(boundary) points, triangles = utils.generate_mesh( ensure_unique(np.concatenate(coords, axis=0)), min_points=min_points[name], max_edge_length=max_edge_length[name], boundary=boundary, convex_hull=False, preserve_boundary=preserve_boundary or (film_terminals is not None), **meshpy_kwargs, ) if smooth[name]: meshes[name] = Mesh.from_triangulation( points, triangles, build_operators=False ).smooth(smooth[name]) else: meshes[name] = Mesh.from_triangulation(points, triangles) self.meshes = meshes
[docs] def boundary_vertices(self, film: str) -> np.ndarray: """An array of boundary vertex indices, ordered counterclockwise. Args: film: The name of the film for which to find boundary indices. Returns: An array of indices for vertices that are on the film boundary, ordered counterclockwise. """ if self.meshes is None: return None mesh = self.meshes[film] points = mesh.sites triangles = mesh.elements indices = utils.boundary_vertices(points, triangles) if film not in self.terminals: return indices # Ensure that the indices wrap around outside of any terminals. for terminal in self.terminals[film]: boundary_points = points[indices] terminal_indices = terminal.contains_points(boundary_points, index=True) discont = np.diff(terminal_indices) != 1 if np.any(discont): i_discont = np.where(discont)[0][0] indices = np.roll(indices, -(i_discont + 1)) break return indices
[docs] def mesh_stats_dict(self) -> Optional[Dict[str, Dict[str, Union[int, float]]]]: """Returns a dictionary of information about all meshes.""" if self.meshes is None: return None return {name: mesh.stats() for name, mesh in self.meshes.items()}
[docs] def mesh_stats(self, precision: int = 3) -> Optional[HTML]: """When called with in Jupyter notebook, displays a table of information about the mesh. Args: precision: Number of digits after the decimal for float values. Returns: An HTML table of mesh statistics. """ all_stats = self.mesh_stats_dict() if all_stats is None: return None def make_row(*cols): return "<tr>" + "".join([f"<td>{c}</td>" for c in cols]) + "</tr>" html = ["<table>", "<tr><h2>Mesh Statistics</h2></tr>"] html.append(make_row("", "<b>length_units</b>", repr(self.length_units))) for name, stats in all_stats.items(): for i, (key, value) in enumerate(stats.items()): if isinstance(value, float): value = f"{value:.{precision}e}" if i == 0: html.append(make_row(f"<b>{name!r}</b>", f"<b>{key}</b>", value)) else: html.append(make_row("", f"<b>{key}</b>", value)) html.append("</table>") return HTML("".join(html))
[docs] def mutual_inductance_matrix( self, hole_polygon_mapping: Optional[Dict[str, np.ndarray]] = None, units: str = "pH", all_iterations: bool = False, progress_bar: bool = False, **solve_kwargs, ) -> Union[np.ndarray, List[np.ndarray]]: """Calculates the mutual inductance matrix :math:`\\mathbf{M}` for the Device. :math:`\\mathbf{M}` is defined such that element :math:`M_{ij}` is the mutual inductance :math:`\\Phi^f_{S_i} / I_j` between hole :math:`j` and polygon :math:`S_{i}`, which encloses hole :math:`i`. :math:`\\Phi^f_{S_i}` is the fluxoid for polygon :math:`S_i` and :math:`I_j` is the current circulating around hole :math:`j`. Args: hole_polygon_mapping: A dict of ``{hole_name: polygon_coordinates}`` specifying a mapping between holes in the device and polygons enclosing those holes, for which the fluxoid will be calculated. The length of this dict, ``n_holes``, will determine the dimension of the square mutual inductance matrix :math:`M`. units: The units in which to report the mutual inductance. all_iterations: Whether to return mutual inductance matrices for all ``iterations + 1`` solutions, or just the final solution. progress_bar: Display a progress bar. solve_kwargs: Keyword arguments passed to :func:`superscreen.solve.solve`, e.g. ``iterations``. Returns: If all_iterations is False, returns a shape ``(n_holes, n_holes)`` mutual inductance matrix from the final iteration. Otherwise, returns a list of mutual inductance matrices, each with shape ``(n_holes, n_holes)``. The length of the list is ``1`` if the device has a single layer, or ``iterations + 1`` if the device has multiple layers. """ from ..solver import factorize_model, solve holes = self.holes if hole_polygon_mapping is None: from ..fluxoid import make_fluxoid_polygons hole_polygon_mapping = make_fluxoid_polygons(self) n_holes = len(hole_polygon_mapping) for hole_name, polygon in hole_polygon_mapping.items(): if hole_name not in holes: raise ValueError(f"Hole '{hole_name}' does not exist in the device.") if not fem.in_polygon(polygon, holes[hole_name].points).all(): raise ValueError( f"Hole '{hole_name}' is not completely contained " f"within the given polygon." ) solve_kwargs = solve_kwargs.copy() iterations = solve_kwargs.get("iterations", 1) solve_kwargs["current_units"] = None solve_kwargs["progress_bar"] = False # The magnitude of this current is not important I_circ = self.ureg("1 mA") if all_iterations: n_iter = 1 if len(self.layers) == 1 else iterations + 1 solution_slice = slice(None) else: n_iter = 1 solution_slice = slice(-1, None) mutual_inductance = np.zeros((n_iter, n_holes, n_holes)) films_by_hole = {} for film, holes in self.holes_by_film().items(): for hole in holes: films_by_hole[hole.name] = film model = None for j, hole_name in enumerate( tqdm(hole_polygon_mapping, desc="Holes", disable=(not progress_bar)) ): logger.info( f"Evaluating {self.name!r} mutual inductance matrix " f"column ({j+1}/{len(hole_polygon_mapping)}), " f"source = {hole_name!r}." ) if model is None: model = factorize_model( device=self, current_units="mA", circulating_currents={hole_name: str(I_circ)}, ) I_circ_val = model.circulating_currents[hole_name] else: model.set_circulating_currents({hole_name: I_circ_val}) solutions = solve(device=None, model=model, **solve_kwargs)[solution_slice] for n, solution in enumerate(solutions): logger.info( f"Evaluating fluxoids for solution {n + 1}/{len(solutions)}." ) for i, (name, polygon) in enumerate(hole_polygon_mapping.items()): fluxoid = solution.polygon_fluxoid( polygon, film=films_by_hole[name] ) mutual_inductance[n, i, j] = ( (sum(fluxoid) / I_circ).to(units).magnitude ) mutual_inductance = mutual_inductance * self.ureg(units) # Return a list to make it clear that we are returning n_iter distinct # matrices. You can convert back to an ndarray using # np.stack(result, axis=0). result = [m for m in mutual_inductance] if not all_iterations: assert len(result) == 1 result = result[0] return result
[docs] def plot_polygons( self, ax: Optional[plt.Axes] = None, subplots: bool = False, legend: bool = False, figsize: Optional[Tuple[float, float]] = None, **kwargs, ) -> Tuple[plt.Figure, plt.Axes]: """Plot all of the Device's polygons. Args: ax: matplotlib axis on which to plot. If None, a new figure is created. subplots: If True, plots each film on a different subplot. legend: Whether to add a legend. figsize: matplotlib figsize, only used if ax is None. kwargs: Passed to ``ax.plot()`` for the polygon boundaries. Returns: Matplotlib Figure and Axes """ if len(self.films) > 1 and subplots and ax is not None: raise ValueError( "Axes may not be provided if subplots is True and the device has " "multiple films." ) if ax is None: if subplots: from ..visualization import auto_grid fig, axes = auto_grid( len(self.films), max_cols=2, figsize=figsize, constrained_layout=True, ) else: fig, axes = plt.subplots(figsize=figsize, constrained_layout=True) axes = np.array([axes for _ in self.films]) else: subplots = False fig = ax.get_figure() axes = np.array([ax for _ in self.films]) holes_by_film = self.holes_by_film() terminals = self.terminals for ax, (name, film) in zip(axes.flat, self.films.items()): _ = film.plot(ax=ax, **kwargs) for hole in holes_by_film[name]: _ = hole.plot(ax=ax, **kwargs) if name in terminals: for terminal in terminals[name]: _ = terminal.plot(ax=ax, **kwargs) if subplots: ax.set_title(name) if legend: ax.legend(bbox_to_anchor=(1, 1), loc="upper left") units = self.ureg(self.length_units).units ax.set_xlabel(f"$x$ $[{units:~L}]$") ax.set_ylabel(f"$y$ $[{units:~L}]$") ax.set_aspect("equal") if not subplots: axes = axes[0] return fig, axes
[docs] def plot_mesh( self, ax: Optional[plt.Axes] = None, subplots: bool = False, figsize: Optional[Tuple[float, float]] = None, show_sites: bool = False, show_edges: bool = True, site_color: Union[str, Sequence[float], None] = None, edge_color: Union[str, Sequence[float], None] = None, linewidth: float = 0.75, linestyle: str = "-", marker: str = ".", ) -> Tuple[plt.Figure, plt.Axes]: """Plot all of the Device's meshes. Args: ax: matplotlib axis on which to plot. If None, a new figure is created. subplots: If True, plots each film on a different subplot. figsize: matplotlib figsize, only used if ax is None. show_sites: Whether to show the mesh sites. show_edges: Whether to show the mesh edges. site_color: The color for the sites. edge_color: The color for the edges. linewidth: The line width for all edges. linestyle: The line style for all edges. marker: The marker to use for the mesh sites. Returns: Matplotlib Figure and Axes """ if len(self.films) > 1 and subplots and ax is not None: raise ValueError( "Axes may not be provided if subplots is True and the device has " "multiple films." ) if self.meshes is None: raise ValueError( "Mesh doesn't exist. Run Device.make_mesh() to generate one." ) if ax is None: if subplots: from ..visualization import auto_grid fig, axes = auto_grid( len(self.films), max_cols=2, figsize=figsize, constrained_layout=True, ) else: fig, axes = plt.subplots(figsize=figsize, constrained_layout=True) axes = np.array([axes for _ in self.films]) else: subplots = False fig = ax.get_figure() axes = np.array([ax for _ in self.films]) for i, (ax, (name, mesh)) in enumerate(zip(axes.flat, self.meshes.items())): sc = f"C{i}" if site_color is None else site_color ec = f"C{i}" if edge_color is None else edge_color ax = mesh.plot( ax=ax, show_sites=show_sites, show_edges=show_edges, site_color=sc, edge_color=ec, linestyle=linestyle, linewidth=linewidth, marker=marker, ) if subplots: ax.set_title(name) units = self.ureg(self.length_units).units ax.set_xlabel(f"$x$ $[{units:~L}]$") ax.set_ylabel(f"$y$ $[{units:~L}]$") ax.set_aspect("equal") if not subplots: axes = axes[0] return fig, axes
[docs] def patches(self) -> Dict[str, Dict[str, PathPatch]]: """Returns a dict of ``{layer_name: {film_name: PathPatch}}`` for visualizing the device. `""" abstract_regions = self.abstract_regions polygons_by_layer = self.polygons_by_layer() holes_by_layer = self.polygons_by_layer(polygon_type="hole") patches = defaultdict(dict) for layer, regions in polygons_by_layer.items(): for region in regions: if region.name in holes_by_layer[layer]: continue coords = region.points.tolist() codes = [Path.LINETO for _ in coords] codes[0] = Path.MOVETO codes[-1] = Path.CLOSEPOLY poly = region.polygon for hole in holes_by_layer[layer]: if region.name not in abstract_regions and poly.contains( hole.polygon ): hole_coords = hole.points.tolist()[::-1] hole_codes = [Path.LINETO for _ in hole_coords] hole_codes[0] = Path.MOVETO hole_codes[-1] = Path.CLOSEPOLY coords.extend(hole_coords) codes.extend(hole_codes) patches[layer][region.name] = PathPatch(Path(coords, codes)) return dict(patches)
[docs] def draw( self, ax: Optional[plt.Axes] = None, subplots: bool = False, max_cols: int = 3, legend: bool = False, figsize: Optional[Tuple[float, float]] = None, alpha: float = 0.5, exclude: Optional[Union[str, List[str]]] = None, layer_order: str = "increasing", ) -> Tuple[plt.Figure, Union[plt.Axes, np.ndarray]]: """Draws all polygons in the device as matplotlib patches. Args: ax: matplotlib axis on which to plot. If None, a new figure is created. subplots: If True, plots each layer on a different subplot. max_cols: The maximum number of columns to create if ``subplots`` is True. legend: Whether to add a legend. figsize: matplotlib figsize, only used if ax is None. alpha: The alpha (opacity) value for the patches (0 <= alpha <= 1). exclude: A polygon name or list of polygon names to exclude from the figure. layer_order: If ``"increasing"`` (``"decreasing"``) draw polygons in increasing (decreasing) order by layer height ``layer.z0``. Returns: Matplotlib Figre and Axes. """ if len(self.layers) > 1 and subplots and ax is not None: raise ValueError( "Axes may not be provided if subplots is True and the device has " "multiple layers." ) layer_order = layer_order.lower() layer_orders = ("increasing", "decreasing") if layer_order not in layer_orders: raise ValueError( f"Invalid layer_order: {layer_order}. " f"Valid layer orders are {layer_orders}." ) if ax is None: if subplots: from ..visualization import auto_grid fig, axes = auto_grid( len(self.layers), max_cols=max_cols, figsize=figsize, constrained_layout=True, ) else: fig, ax = plt.subplots(figsize=figsize, constrained_layout=True) axes = np.array([ax for _ in self.layers]) else: subplots = False fig = ax.get_figure() axes = np.array([ax for _ in self.layers]) exclude = exclude or [] if isinstance(exclude, str): exclude = [exclude] layers = [ layer.name for layer in sorted(self.layers.values(), key=lambda x: x.z0) ] if layer_order == "decreasing": layers = layers[::-1] patches = self.patches() used_axes = set() units = self.ureg(self.length_units).units x, y = self.poly_points.T margin = 0.1 dx = np.ptp(x) dy = np.ptp(y) x0 = x.min() + dx / 2 y0 = y.min() + dy / 2 dx *= 1 + margin dy *= 1 + margin labels = [] handles = [] for i, (layer, ax) in enumerate(zip(layers, axes.flat)): ax.set_aspect("equal") ax.grid(False) ax.set_xlim(x0 - dx / 2, x0 + dx / 2) ax.set_ylim(y0 - dy / 2, y0 + dy / 2) ax.set_xlabel(f"$x$ $[{units:~L}]$") ax.set_ylabel(f"$y$ $[{units:~L}]$") if subplots: labels = [] handles = [] j = 0 for name, patch in patches[layer].items(): if name in exclude or name in self.holes: continue patch.set_facecolor(f"C{i}") patch.set_alpha(alpha) ax.add_artist(patch) used_axes.add(ax) if j == 0: labels.append(layer) handles.append(patch) j += 1 if subplots: ax.set_title(layer) if legend: ax.legend(handles, labels, bbox_to_anchor=(1, 1), loc="upper left") if subplots: for ax in fig.axes: if ax not in used_axes: fig.delaxes(ax) else: axes = axes[0] if legend: axes.legend(handles, labels, bbox_to_anchor=(1, 1), loc="upper left") return fig, axes
[docs] def to_hdf5( self, path_or_group: Union[os.PathLike, h5py.Group], save_mesh: bool = True, compress: bool = True, ) -> None: """Serializes the Device to and HDF5 file. Args: path_or_group: Path to an HDF5 file to create, or an open HD5F file. save_mesh: Whether to save the full mesh to file. compress: Save the minimum amount of data needed to recreate the mesh. """ if isinstance(path_or_group, h5py.Group): save_context = nullcontext(path_or_group) else: save_context = h5py.File(path_or_group, "x") with save_context as h5group: h5group.attrs["name"] = self.name h5group.attrs["length_units"] = self.length_units h5group.attrs["solve_dtype"] = str(self.solve_dtype) layer_grp = h5group.create_group("layers") film_grp = h5group.create_group("films") hole_grp = h5group.create_group("holes") terminals_grp = h5group.create_group("terminals") abs_grp = h5group.create_group("abstract_regions") for name, layer in self.layers.items(): layer.to_hdf5(layer_grp.create_group(name)) for name, polygon in self.films.items(): polygon.to_hdf5(film_grp.create_group(name)) for name, polygon in self.holes.items(): polygon.to_hdf5(hole_grp.create_group(name)) for name, polygon in self.abstract_regions.items(): polygon.to_hdf5(abs_grp.create_group(name)) for film_name, terminals in self.terminals.items(): grp = terminals_grp.create_group(film_name) for i, terminal in enumerate(terminals): terminal.to_hdf5(grp.create_group(str(i))) if save_mesh and self.meshes: mesh_grp = h5group.create_group("mesh") for name, mesh in self.meshes.items(): mesh.to_hdf5(mesh_grp.create_group(name), compress=compress)
[docs] @staticmethod def from_hdf5(path_or_group: Union[os.PathLike, h5py.Group]) -> "Device": """Loads a Device from an HDF5 file. Args: path_or_group: Path to an HDF5 file to read, or an open HD5F file. Returns: The deserialized Device. """ if isinstance(path_or_group, h5py.Group): read_context = nullcontext(path_or_group) else: read_context = h5py.File(path_or_group, "r") with read_context as h5group: terminals = {} for film, grp in h5group["terminals"].items(): terminals[film] = [] for i in range(len(grp)): terminals[film].append(Polygon.from_hdf5(grp[str(i)])) device = Device( name=h5group.attrs["name"], layers=[Layer.from_hdf5(grp) for grp in h5group["layers"].values()], films=[Polygon.from_hdf5(grp) for grp in h5group["films"].values()], holes=[Polygon.from_hdf5(grp) for grp in h5group["holes"].values()], terminals=terminals, abstract_regions=[ Polygon.from_hdf5(grp) for grp in h5group["abstract_regions"].values() ], length_units=h5group.attrs["length_units"], solve_dtype=h5group.attrs["solve_dtype"], ) if "mesh" in h5group: device.meshes = { name: Mesh.from_hdf5(grp) for name, grp in h5group["mesh"].items() } return device
def __repr__(self) -> str: # Normal tab "\t" renders a bit too big in jupyter if you ask me. indent = 4 t = " " * indent nt = "\n" + t def format_list(L): if not L: return None items = [f"{t}{value}" for value in L] return "[" + nt + (", " + nt).join(items) + "," + nt + "]" def format_dict(D): if not D: return None items = [f"{t}{key!r}: {value}" for key, value in D.items()] return "{" + nt + (", " + nt).join(items) + "," + nt + "}" args = [ f'"{self.name}"', f"layers={format_list(self.layers.values())}", f"films={format_list(self.films.values())}", f"holes={format_list(self.holes.values())}", f"terminals={format_dict(self.terminals)}", f"abstract_regions={format_list(self.abstract_regions.values())}", f'length_units="{self.length_units}"', ] return f"{self.__class__.__name__}(" + nt + (", " + nt).join(args) + ",\n)" def __eq__(self, other) -> bool: if other is self: return True if not isinstance(other, Device): return False def equals_sorted(first, second): def key(x): return x.name return sorted(first, key=key) == sorted(second, key=key) return ( self.name == other.name and equals_sorted(self.layers.values(), other.layers.values()) and equals_sorted(self.films.values(), other.films.values()) and equals_sorted(self.holes.values(), other.holes.values()) and self.terminals == other.terminals and equals_sorted( self.abstract_regions.values(), other.abstract_regions.values() ) and self.length_units == other.length_units ) def __getstate__(self): state = self.__dict__.copy() # Use dill for Layer objects because Layer.Lambda could be a Parameter state["layers"] = dill.dumps(self.layers) return state def __setstate__(self, state): # Use dill for Layer objects because Layer.Lambda could be a Parameter state["layers"] = dill.loads(state["layers"]) self.__dict__.update(state)