"""Field-probe orientation and coordinate transformation helpers.
The functions in this module transform measurements from probe coordinates
into TEM-cell coordinates. They deliberately distinguish signed vectors and
synchronous waveforms from component magnitudes: a general three-dimensional
rotation cannot be reconstructed from three unsigned component magnitudes.
"""
import ast
import math
import numpy as np
from scuq.quantities import Quantity
from scuq.si import METER, VOLT
EFIELD = VOLT / METER
CELL_AXES = ("cell_x", "cell_y", "cell_z")
PROBE_AXES = ("probe_x", "probe_y", "probe_z")
PROBE_AXIS_INDEX = {axis: index for index, axis in enumerate(PROBE_AXES)}
PROBE_ORIENTATION_TOLERANCE = 1e-9
SIGNED_VECTOR = "signed_vector"
COMPONENT_MAGNITUDES = "component_magnitudes"
WAVEFORM = "waveform"
PROBE_DATA_KINDS = (SIGNED_VECTOR, COMPONENT_MAGNITUDES, WAVEFORM)
DEFAULT_PROBE_AXIS_MAP = {
"cell_x": "+probe_x",
"cell_y": "+probe_y",
"cell_z": "+probe_z",
}
ORIENTATION_KEYS = (
"probe_axis_map",
"probe_rotation_matrix",
"probe_rotation_angles_deg",
)
class ProbeOrientationConflictError(ValueError):
"""Report contradictory orientation definitions for one field probe.
Parameters
----------
graph_name : str
Physical field-probe node name.
candidates : iterable of mapping
Normalized source and orientation entries participating in the
conflict.
Attributes
----------
graph_name : str
Physical field-probe node name.
candidates : tuple of dict
Normalized source and orientation entries participating in the
conflict.
"""
def __init__(self, graph_name, candidates):
self.graph_name = str(graph_name)
self.candidates = tuple(
{
"source": str(candidate["source"]),
"orientation": candidate["orientation"],
}
for candidate in candidates
)
detail = "; ".join(
"%s=%r" % (candidate["source"], candidate["orientation"])
for candidate in self.candidates
)
super().__init__(
"conflicting probe orientation for field probe %r: %s"
% (self.graph_name, detail)
)
[docs]
def parse_probe_axis_map(axis_map=None):
"""Parse a complete signed mapping from probe axes to cell axes.
Parameters
----------
axis_map : mapping of str to str, optional
Mapping for ``cell_x``, ``cell_y``, and ``cell_z``. Values use signed
probe-axis expressions such as ``+probe_y`` or ``-probe_x``.
Returns
-------
dict
Parsed probe-axis index, sign, and normalized expression for every
cell axis.
"""
raw = dict(DEFAULT_PROBE_AXIS_MAP)
if axis_map is not None:
raw.update(axis_map)
parsed = {}
for cell_axis in CELL_AXES:
if cell_axis not in raw:
raise ValueError("probe_axis_map needs mapping for %s" % cell_axis)
token = str(raw[cell_axis]).strip()
sign = 1
if token.startswith("+"):
token = token[1:]
elif token.startswith("-"):
token = token[1:]
sign = -1
if token in ("x", "y", "z"):
token = "probe_" + token
if token not in PROBE_AXIS_INDEX:
raise ValueError(
"invalid probe axis mapping for %s: %r"
% (cell_axis, raw[cell_axis])
)
parsed[cell_axis] = {
"probe_axis": token,
"probe_index": PROBE_AXIS_INDEX[token],
"sign": sign,
"expression": ("+" if sign > 0 else "-") + token,
}
unknown = sorted(set(raw) - set(CELL_AXES))
if unknown:
raise ValueError(
"probe_axis_map contains unsupported cell axes: %s"
% ", ".join(unknown)
)
return parsed
def _probe_axis_map_to_matrix(axis_map=None):
parsed = parse_probe_axis_map(axis_map)
matrix = []
for cell_axis in CELL_AXES:
row = [0.0, 0.0, 0.0]
entry = parsed[cell_axis]
row[entry["probe_index"]] = float(entry["sign"])
matrix.append(row)
return matrix
def _coerce_probe_rotation_matrix(value):
if value is None:
return None
if isinstance(value, str):
text = value.strip().strip('"\'')
if not text:
return None
try:
value = ast.literal_eval(text)
except (ValueError, SyntaxError):
value = [
[item.strip() for item in row_text.split(",")]
for row_text in text.split(";")
if row_text.strip()
]
rows = [list(row) for row in value]
if len(rows) != 3 or any(len(row) != 3 for row in rows):
raise ValueError("probe_rotation_matrix must be a 3x3 matrix")
matrix = []
for row in rows:
matrix_row = []
for item in row:
try:
number = float(item)
except (TypeError, ValueError) as exc:
raise ValueError(
"probe_rotation_matrix entries must be real numbers"
) from exc
if not math.isfinite(number):
raise ValueError("probe_rotation_matrix entries must be finite")
matrix_row.append(number)
matrix.append(matrix_row)
_validate_probe_rotation_matrix(matrix)
return matrix
def _coerce_probe_rotation_angles_deg(value):
if value is None:
return None
if isinstance(value, str):
text = value.strip().strip('"\'')
if not text:
return None
try:
value = ast.literal_eval(text)
except (ValueError, SyntaxError) as exc:
raise ValueError(
"probe_rotation_angles_deg must be a mapping"
) from exc
if not isinstance(value, dict):
raise ValueError("probe_rotation_angles_deg must be a mapping")
allowed = ("about_cell_z", "about_cell_x", "about_cell_y")
unknown = sorted(set(value) - set(allowed))
if unknown:
raise ValueError(
"probe_rotation_angles_deg contains unsupported keys: %s"
% ", ".join(unknown)
)
angles = {}
for key in allowed:
try:
angle = float(value.get(key, 0.0))
except (TypeError, ValueError) as exc:
raise ValueError(
"probe_rotation_angles_deg entries must be real numbers"
) from exc
if not math.isfinite(angle):
raise ValueError("probe_rotation_angles_deg entries must be finite")
angles[key] = angle
return angles
def _matrix_multiply(left, right):
return [
[
sum(left[row][index] * right[index][col] for index in range(3))
for col in range(3)
]
for row in range(3)
]
def _rotation_matrix_about_cell_axis(axis, angle_deg):
angle = math.radians(angle_deg)
c = math.cos(angle)
s = math.sin(angle)
if axis == "cell_x":
return [[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]]
if axis == "cell_y":
return [[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]]
if axis == "cell_z":
return [[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]]
raise ValueError("unsupported cell axis %r" % axis)
def _probe_rotation_matrix_from_angles_deg(value):
angles = _coerce_probe_rotation_angles_deg(value)
if angles is None:
return None
matrix = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
# Fixed-cell-axis convention: rotate about cell z, then x, then y.
for axis, key in (
("cell_z", "about_cell_z"),
("cell_x", "about_cell_x"),
("cell_y", "about_cell_y"),
):
matrix = _matrix_multiply(
_rotation_matrix_about_cell_axis(axis, angles[key]), matrix
)
_validate_probe_rotation_matrix(matrix)
return matrix
def _matrix_product_left_transpose(matrix):
return [
[
sum(matrix[row][i] * matrix[row][j] for row in range(3))
for j in range(3)
]
for i in range(3)
]
def _matrix_determinant(matrix):
return (
matrix[0][0]
* (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])
- matrix[0][1]
* (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])
+ matrix[0][2]
* (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0])
)
def _validate_probe_rotation_matrix(matrix):
product = _matrix_product_left_transpose(matrix)
for row_index, row in enumerate(product):
for col_index, value in enumerate(row):
expected = 1.0 if row_index == col_index else 0.0
if abs(value - expected) > PROBE_ORIENTATION_TOLERANCE:
raise ValueError("probe_rotation_matrix must be orthonormal")
determinant = _matrix_determinant(matrix)
if abs(determinant - 1.0) > PROBE_ORIENTATION_TOLERANCE:
raise ValueError(
"probe_rotation_matrix must describe a proper rotation with "
"determinant +1"
)
def _axis_map_from_rotation_matrix(matrix):
used_columns = set()
axis_map = {}
for row_index, cell_axis in enumerate(CELL_AXES):
row = matrix[row_index]
non_zero = [
col_index
for col_index, value in enumerate(row)
if abs(value) > PROBE_ORIENTATION_TOLERANCE
]
if len(non_zero) != 1:
return None
col_index = non_zero[0]
if col_index in used_columns:
return None
value = row[col_index]
if abs(abs(value) - 1.0) > PROBE_ORIENTATION_TOLERANCE:
return None
used_columns.add(col_index)
axis_map[cell_axis] = ("+" if value > 0 else "-") + PROBE_AXES[col_index]
return axis_map
def _normalized_probe_rotation_matrix(matrix):
return [
[
0.0 if abs(value) < PROBE_ORIENTATION_TOLERANCE else float(value)
for value in row
]
for row in matrix
]
def _probe_axis_map_from_config_value(value):
if value is None:
return None
if isinstance(value, dict):
return value
if isinstance(value, str):
text = value.strip().strip('"\'')
if not text:
return None
if text.startswith("{"):
try:
parsed = ast.literal_eval(text)
except (ValueError, SyntaxError) as exc:
raise ValueError("probe_axis_map must be a mapping") from exc
if not isinstance(parsed, dict):
raise ValueError("probe_axis_map must be a mapping")
return parsed
parsed = {}
for item in text.split(","):
key, separator, expression = item.partition(":")
if not separator:
raise ValueError(
"probe_axis_map item %r must use 'cell_axis:probe_axis' "
"syntax" % item
)
parsed[key.strip()] = expression.strip()
return parsed
raise TypeError("probe_axis_map must be a mapping or string")
def _probe_orientation_from_config_value(value):
if value is None:
return None
if isinstance(value, dict):
orientation_keys = [
key
for key in (
"probe_rotation_matrix",
"rotation_matrix",
"probe_rotation_angles_deg",
"probe_axis_map",
)
if key in value and value.get(key) is not None
]
if len(orientation_keys) > 1:
normalized = [
parse_probe_orientation({key: value[key]})["normalized"]
for key in orientation_keys
]
if any(item != normalized[0] for item in normalized[1:]):
raise ValueError("conflicting probe orientation entries")
return {
"probe_rotation_matrix": parse_probe_orientation(
{orientation_keys[0]: value[orientation_keys[0]]}
)["matrix"]
}
if "probe_rotation_matrix" in value:
return {
"probe_rotation_matrix": _coerce_probe_rotation_matrix(
value["probe_rotation_matrix"]
)
}
if "rotation_matrix" in value:
return {
"probe_rotation_matrix": _coerce_probe_rotation_matrix(
value["rotation_matrix"]
)
}
if "probe_rotation_angles_deg" in value:
return {
"probe_rotation_matrix": _probe_rotation_matrix_from_angles_deg(
value["probe_rotation_angles_deg"]
)
}
if "probe_axis_map" in value:
return {
"probe_axis_map": _probe_axis_map_from_config_value(
value["probe_axis_map"]
)
}
return {"probe_axis_map": _probe_axis_map_from_config_value(value)}
[docs]
def parse_probe_orientation(orientation=None):
"""Parse a field-probe orientation into a cell-from-probe matrix.
Parameters
----------
orientation : mapping, str, or None, optional
Signed axis map, rotation matrix, or fixed-cell-axis rotation angles.
``None`` selects :data:`DEFAULT_PROBE_AXIS_MAP`.
Returns
-------
dict
Rotation ``matrix``, an ``axis_map`` when the rotation is a signed
permutation, and its ``normalized`` representation.
"""
orientation = _probe_orientation_from_config_value(orientation)
if orientation is None:
orientation = {"probe_axis_map": None}
if "probe_rotation_matrix" in orientation:
matrix = _coerce_probe_rotation_matrix(
orientation["probe_rotation_matrix"]
)
else:
matrix = _probe_axis_map_to_matrix(orientation.get("probe_axis_map"))
axis_map = _axis_map_from_rotation_matrix(matrix)
return {
"matrix": matrix,
"axis_map": axis_map,
"normalized": (
axis_map
if axis_map is not None
else {"probe_rotation_matrix": _normalized_probe_rotation_matrix(matrix)}
),
}
def normalize_probe_orientation(orientation=None):
"""Return the normalized representation of a probe orientation.
Parameters
----------
orientation : mapping, str, or None, optional
Signed axis map, rotation matrix, or fixed-cell-axis rotation angles.
Returns
-------
dict
Signed axis map for a permutation, otherwise a normalized rotation
matrix mapping probe coordinates to cell coordinates.
"""
return parse_probe_orientation(orientation)["normalized"]
def _coerce_field_components(data):
if not isinstance(data, (list, tuple)):
raise ValueError("probe field data must contain three components")
values = []
for index, probe_axis in enumerate(PROBE_AXES):
try:
value = data[index]
except IndexError as exc:
raise ValueError("probe field data is missing %s" % probe_axis) from exc
if not isinstance(value, Quantity):
value = Quantity(EFIELD, value)
values.append(value.reduce_to(EFIELD))
return values
def _map_probe_field_vector_parsed(data, parsed, data_kind=SIGNED_VECTOR):
if data_kind not in (SIGNED_VECTOR, COMPONENT_MAGNITUDES):
raise ValueError(
"vector data_kind must be %r or %r"
% (SIGNED_VECTOR, COMPONENT_MAGNITUDES)
)
probe_values = _coerce_field_components(data)
if data_kind == COMPONENT_MAGNITUDES:
axis_map = parsed.get("axis_map")
if axis_map is None:
raise ValueError(
"general probe rotations require signed vector or waveform "
"data; component magnitudes only support axis permutations"
)
parsed_axis_map = parse_probe_axis_map(axis_map)
return {
cell_axis: abs(probe_values[entry["probe_index"]]).reduce_to(EFIELD)
for cell_axis, entry in parsed_axis_map.items()
}
matrix = np.asarray(parsed["matrix"], dtype=float)
cell_values = matrix @ np.asarray(probe_values, dtype=object)
return {
cell_axis: cell_values[row_index].reduce_to(EFIELD)
for row_index, cell_axis in enumerate(CELL_AXES)
}
[docs]
def map_probe_field_vector(data, orientation=None, *, data_kind=SIGNED_VECTOR):
"""Map a three-axis field reading into cell coordinates.
Parameters
----------
data : sequence of Quantity or float
Probe x-, y-, and z-axis field components.
orientation : mapping, str, or None, optional
Probe orientation accepted by :func:`parse_probe_orientation`.
data_kind : {"signed_vector", "component_magnitudes"}, optional
``signed_vector`` applies the complete rotation matrix. For
``component_magnitudes`` only an axis permutation is physically
identifiable; signs are ignored and general rotations are rejected.
Returns
-------
dict
``cell_x``, ``cell_y``, and ``cell_z`` field quantities.
"""
parsed = parse_probe_orientation(orientation)
return _map_probe_field_vector_parsed(data, parsed, data_kind=data_kind)
def map_probe_field_waveform(data, orientation=None):
"""Rotate synchronous signed probe waveforms into cell coordinates.
Parameters
----------
data : sequence of three array-like objects
Synchronous signed samples for probe x, y, and z. All arrays must have
equal one-dimensional shapes.
orientation : mapping, str, or None, optional
Probe orientation accepted by :func:`parse_probe_orientation`.
Returns
-------
dict
Arrays for ``cell_x``, ``cell_y``, and ``cell_z``.
"""
if not isinstance(data, (list, tuple)) or len(data) != 3:
raise ValueError("probe waveform data must contain exactly three axes")
probe_values = tuple(np.asarray(axis, dtype=float) for axis in data)
if any(axis.ndim != 1 for axis in probe_values):
raise ValueError("probe waveform axes must be one-dimensional")
if any(axis.shape != probe_values[0].shape for axis in probe_values[1:]):
raise ValueError("probe waveform axes must have equal lengths")
matrix = np.asarray(parse_probe_orientation(orientation)["matrix"], dtype=float)
cell_values = matrix @ np.asarray(probe_values)
return {
cell_axis: cell_values[row_index]
for row_index, cell_axis in enumerate(CELL_AXES)
}
def mapped_probe_reading(
value,
orientation,
*,
data_kind=SIGNED_VECTOR,
):
"""Return cell-axis and unchanged probe-axis values for one reading.
Parameters
----------
value : sequence of Quantity or float
Probe x-, y-, and z-axis field components.
orientation : mapping, str, or None
Probe orientation accepted by :func:`parse_probe_orientation`.
data_kind : {"signed_vector", "component_magnitudes"}, optional
Physical interpretation of the three input components.
Returns
-------
tuple of (list, list)
Transformed cell-coordinate vector followed by the unchanged probe
vector.
"""
mapped = map_probe_field_vector(value, orientation, data_kind=data_kind)
return [mapped[axis] for axis in CELL_AXES], list(value)
def _orientation_from_section(section, source):
if not isinstance(section, dict):
return None
entries = {
key: section.get(key)
for key in ORIENTATION_KEYS
if key in section and section.get(key) is not None
}
if not entries:
return None
try:
normalized = normalize_probe_orientation(entries)
except ValueError as exc:
if len(entries) > 1:
raise ValueError(
"conflicting probe orientation in %s" % source
) from exc
raise
return {"source": source, "orientation": normalized}
def _node_attributes(measurement_graph, graph_name):
try:
node = measurement_graph.nodes[graph_name]
except (AttributeError, KeyError, TypeError):
return {}
attributes = {}
graph_node = node.get("gnode") if isinstance(node, dict) else None
get_attributes = getattr(graph_node, "get_attributes", None)
if callable(get_attributes):
attributes.update(get_attributes())
if isinstance(node, dict):
attributes.update(
{key: value for key, value in node.items() if key != "gnode"}
)
return attributes
def _orientation_from_dot(measurement_graph, graph_name):
return _orientation_from_section(
_node_attributes(measurement_graph, graph_name),
"dot:%s" % graph_name,
)
def _orientation_from_ini(measurement_graph, graph_name):
try:
node = measurement_graph.nodes[graph_name]
except (AttributeError, KeyError, TypeError):
return None
if not isinstance(node, dict):
return None
ini_data = node.get("inidic") or {}
description = ini_data.get("description") or {}
if str(description.get("type", "")).strip().lower() != "fieldprobe":
return None
for section_name in ("description", "init_value"):
candidate = _orientation_from_section(
ini_data.get(section_name) or {},
"ini:%s:%s" % (graph_name, section_name),
)
if candidate is not None:
return candidate
return None
def _configuration_for_probe(config_probe_orientations, graph_name):
if not config_probe_orientations:
return None
if isinstance(config_probe_orientations, dict):
if graph_name in config_probe_orientations:
return config_probe_orientations[graph_name]
if set(config_probe_orientations) & set(ORIENTATION_KEYS):
return {
key: config_probe_orientations[key]
for key in ORIENTATION_KEYS
if key in config_probe_orientations
}
if set(config_probe_orientations) & set(CELL_AXES):
return {"probe_axis_map": config_probe_orientations}
return config_probe_orientations
def _probe_orientation_candidates(
measurement_graph,
graph_name,
config_probe_orientations=None,
):
candidates = []
configured = _configuration_for_probe(config_probe_orientations, graph_name)
if configured is not None:
if isinstance(configured, dict) and set(configured) & set(ORIENTATION_KEYS):
config_section = configured
else:
config_section = {"probe_axis_map": configured}
candidate = _orientation_from_section(
config_section,
"config:%s" % graph_name,
)
if candidate is not None:
candidates.append(candidate)
for candidate in (
_orientation_from_dot(measurement_graph, graph_name),
_orientation_from_ini(measurement_graph, graph_name),
):
if candidate is not None:
candidates.append(candidate)
return candidates
[docs]
def resolve_probe_orientation(
measurement_graph,
graph_name,
config_probe_orientations=None,
):
"""Resolve one physical field probe's orientation and source.
Parameters
----------
measurement_graph : mpylab.tools.mgraph.MGraph or compatible object
Graph containing the field-probe node, DOT attributes, and parsed INI
sections.
graph_name : str
Physical field-probe node name.
config_probe_orientations : mapping, optional
Per-node or direct orientation supplied by the application.
Returns
-------
dict
``orientation`` contains the normalized mapping or matrix, ``source``
identifies the selected source, and ``candidates`` retains all
agreeing explicit definitions.
Raises
------
ProbeOrientationConflictError
If explicit sources define different normalized orientations.
"""
candidates = _probe_orientation_candidates(
measurement_graph,
graph_name,
config_probe_orientations=config_probe_orientations,
)
unique = []
for candidate in candidates:
if candidate["orientation"] not in [
entry["orientation"] for entry in unique
]:
unique.append(candidate)
if len(unique) > 1:
raise ProbeOrientationConflictError(graph_name, unique)
if unique:
orientation = unique[0]["orientation"]
source = unique[0]["source"]
else:
orientation = normalize_probe_orientation(None)
source = "default"
return {
"orientation": orientation,
"source": source,
"candidates": tuple(candidates),
}
def resolve_field_probe_orientations(
measurement_graph,
probe_names,
config_probe_orientations=None,
):
"""Resolve orientations independently for multiple physical probes.
Parameters
----------
measurement_graph : mpylab.tools.mgraph.MGraph or compatible object
Graph containing all requested field-probe nodes.
probe_names : iterable of str or str
Physical field-probe node names.
config_probe_orientations : mapping, optional
Per-node or direct application orientation.
Returns
-------
dict
Resolution result from :func:`resolve_probe_orientation` for every
physical probe node.
"""
if isinstance(probe_names, str):
probe_names = (probe_names,)
return {
graph_name: resolve_probe_orientation(
measurement_graph,
graph_name,
config_probe_orientations=config_probe_orientations,
)
for graph_name in probe_names
}
def _format_orientation_matrix(matrix):
return [
" [% .3f, % .3f, % .3f]" % (row[0], row[1], row[2])
for row in matrix
]
def _format_orientation_axis_expression(cell_axis, matrix):
row = matrix[CELL_AXES.index(cell_axis)]
terms = []
for factor, probe_axis in zip(row, PROBE_AXES):
if abs(factor) <= PROBE_ORIENTATION_TOLERANCE:
continue
terms.append("%+.3f*%s" % (factor, probe_axis))
return "%s = %s" % (cell_axis, " ".join(terms) if terms else "0")
__all__ = [
"CELL_AXES",
"COMPONENT_MAGNITUDES",
"DEFAULT_PROBE_AXIS_MAP",
"PROBE_AXES",
"PROBE_DATA_KINDS",
"ProbeOrientationConflictError",
"SIGNED_VECTOR",
"WAVEFORM",
"format_probe_orientation_preflight",
"map_probe_field_vector",
"map_probe_field_waveform",
"mapped_probe_reading",
"normalize_probe_orientation",
"parse_probe_axis_map",
"parse_probe_orientation",
"resolve_field_probe_orientations",
"resolve_probe_orientation",
]