diff --git a/README.md b/README.md index 05a94bb..9ebf864 100644 --- a/README.md +++ b/README.md @@ -111,7 +111,7 @@ A template saved with File → Save as becomes a configuration that can be edite | jScope | idscope | |---|---| | Shot field + Apply | **Sources** field + Apply: comma-separated sources, optionally `label=source` | -| Expression (TDI) | Signal **path**, e.g. `equilibrium/time_slice/global_quantities/ip`, with scale and offset | +| Expression (TDI) | Signal **path**, e.g. `equilibrium/time_slice/global_quantities/ip`, with scale and offset, or arithmetic of paths (`a - b`, `a / b`, `deriv(a)`) | | Browse signals | **IDS browser** (left): double-click or drag a quantity onto a panel | | Wave setup (double-click) | **Setup** dialog: double-click a panel (or Enter) | | Zoom / Point / Pan modes | **Pan** (default) / **Zoom** / **Point** in the toolbar (P, Z, C) | @@ -297,6 +297,26 @@ distribution_sources/source[:]/ggd/particles[0]/values a map per source, on a `psi` go from the magnetic axis to the boundary, `norm(.../psi)` is psi_n. As x path, `core_profiles/profiles_1d/grid/rho_tor_norm | norm(core_profiles/profiles_1d/grid/psi)` draws against rho_tor_norm, or against psi_n when rho_tor_norm is not filled. +- Signals can be computed from others: `+`, `-`, `*` and `/` combine signals and + numbers, with parentheses, and `deriv(x)` (along time for traces, along the + coordinate for profiles), `abs(x)` and `norm(x)` apply to them. Since `/` also + separates the nodes of a path, a division has spaces around it: `a / b`. + + ``` + equilibrium/time_slice/global_quantities/ip - summary/global_quantities/ip/value + core_profiles/profiles_1d/t_i_average / core_profiles/profiles_1d/electrons/temperature + core_profiles/profiles_1d/electrons/temperature / core_profiles/profiles_1d/electrons/temperature[0] + deriv(summary/global_quantities/ip/value) * 1e-6 + abs(pf_active/coil[:]/current/data) + ``` + + Two signals are combined on the times and the grid of the first: the second is + interpolated in time for traces, and taken at its closest time slice (within its + time range) and interpolated on the grid of the first for profiles and maps. A + time trace applies to profiles and maps at their times, e.g. a profile divided by + its central value. Signals with wildcards combine member by member (by name), or + each member with the same signal. The result is labeled with the expression, and + its units follow (a ratio of like quantities has none). Data stored with another major Data Dictionary version than the one idscope uses (e.g. DD 3 data with DD 4) is shown in its own version, so paths use that version's @@ -323,6 +343,16 @@ When a panel compares several sources, curves take the color of their signal (or member, e.g. coil) and the line style of their source; the legend has an entry for each, instead of one for every combination. +A panel can also show every source but the first relative to the first, signal by +signal (*sources* in the setup dialog, or *Sources* in the menu of the panel, or +`relative:` in the configuration): `difference` shows B − A and `ratio` B / A, on +the times and the grid of B, as the arithmetic of signals does. The first source is +then not drawn, and the legend says what is compared (`[B − A]`). The scale of a +signal still applies to a difference; its offset, which cancels out, does not, and +a ratio has neither. Shapes cannot be compared: a map of the difference keeps the +boundary of each source, for one. With a single source, the panel shows it as it +is. + A graph against time says nothing of data that is all at one and the same time, such as what a code gives for the time it was run for. A trace panel then shows bars of the values at that time (a bar per signal, or per member of a signal with @@ -365,7 +395,8 @@ columns: ``` Panel options: `title`, `kind`, `signals`, `x_label`, `y_label`, `x_range`, `y_range`, -`log_x`, `log_y`, `legend`, `equal_aspect`, `contours`, `height`. Signal options: +`log_x`, `log_y`, `legend`, `equal_aspect`, `contours`, `height`, `relative` +(`difference` or `ratio`: the sources relative to the first). Signal options: `path`, `label`, `scale`, `offset`, `x`, `color`, `sources`. Other top-level options: `share_x` (share x axes and stack panels tightly, default true) and `dd_version`. diff --git a/src/idscope/config.py b/src/idscope/config.py index 4bf198c..661a2e1 100644 --- a/src/idscope/config.py +++ b/src/idscope/config.py @@ -32,6 +32,8 @@ from idscope.errors import ConfigError, IdscopeError KINDS = ("auto", "trace", "profile", "map", "spacetime") +RELATIVE = ("difference", "ratio") +"""How a panel can show the sources relative to the first one.""" FORMAT_VERSION = 1 @@ -77,6 +79,9 @@ class PanelConfig: (0: image).""" height: float = 1.0 """Relative height within the column.""" + relative: Optional[str] = None + """Show every source but the first relative to the first: ``difference`` + (B - A) or ``ratio`` (B / A), signal by signal (default: as they are).""" @dataclass @@ -211,6 +216,9 @@ def _panel(data, where: str) -> PanelConfig: kind = data.get("kind", "auto") if kind not in KINDS: raise ConfigError(f"{where}: 'kind' must be one of {', '.join(KINDS)}.") + relative = data.get("relative") + if relative is not None and relative not in RELATIVE: + raise ConfigError(f"{where}: 'relative' must be one of {', '.join(RELATIVE)}.") signals = data.get("signals") or [] if not isinstance(signals, list): raise ConfigError(f"{where}: 'signals' must be a list.") @@ -228,6 +236,7 @@ def _panel(data, where: str) -> PanelConfig: equal_aspect=_bool(data, "equal_aspect", False, where), contours=_count(data, "contours", where), height=_number(data, "height", 1.0, where), + relative=relative, ) diff --git a/src/idscope/data/__init__.py b/src/idscope/data/__init__.py index e44a048..b585cb2 100644 --- a/src/idscope/data/__init__.py +++ b/src/idscope/data/__init__.py @@ -1,5 +1,6 @@ """Data access layer of idscope: sources, signal paths and extraction.""" +from idscope.data.expressions import Expression, parse_signal from idscope.data.extract import ( MAP, PROFILE, @@ -39,6 +40,7 @@ "Axis", "Channels", "DataSource", + "Expression", "MapSeries", "ProfileSeries", "Range", @@ -51,6 +53,7 @@ "extract_signal", "format_source_list", "parse_alternatives", + "parse_signal", "parse_signal_path", "parse_source", "parse_source_list", diff --git a/src/idscope/data/expressions.py b/src/idscope/data/expressions.py new file mode 100644 index 0000000..1d1a4dc --- /dev/null +++ b/src/idscope/data/expressions.py @@ -0,0 +1,727 @@ +"""Signal expressions: arithmetic on signals. + +A signal can be computed from others and from numbers, e.g.:: + + equilibrium/time_slice/global_quantities/ip - summary/global_quantities/ip/value + core_profiles/profiles_1d/t_i_average / core_profiles/profiles_1d/electrons/temperature + deriv(core_profiles/global_quantities/ip) + abs(pf_active/coil[:]/current/data) * 1e-3 + +Rules: + +- The operators ``+``, ``-``, ``*`` and ``/`` combine signals and numbers, with the + usual precedence and parentheses. Since ``/`` also separates the nodes of a path, + a division is written with spaces around it: ``a / b``, while ``a/b`` is a path + (a ``/`` right after a path, with no space, belongs to the path). +- Functions: ``deriv(x)`` is the derivative along time for time traces and along + the coordinate for profiles, ``abs(x)`` the absolute value, and ``norm(x)`` + normalizes profiles (see :mod:`idscope.data.paths`, also for ``transpose``). +- Two signals are combined on the times and the grid of the first: the second is + interpolated in time for traces, and taken at its closest time slice (within its + time range) and interpolated on the grid of the first for profiles and maps. A + time trace combined with profiles or maps applies to each slice at its time, + e.g. a profile divided by its central value. +- A signal with wildcards (e.g. a trace per coil) is combined member by member with + another one, by name, or with the same signal for all its members. + +The same arithmetic compares data sources in a panel: :func:`relative_to` gives the +difference or the ratio of the data of two sources. +""" + +import operator +import re +from typing import Callable, List, Optional, Union + +import numpy as np + +from idscope.data.extract import ( + Axis, + Channels, + MapSeries, + ProfileSeries, + Shapes, + SignalData, + Trace, + _interp, + normalize, +) +from idscope.data.paths import SignalPath, _split_arguments, parse_signal_path +from idscope.errors import PathError + +FUNCTIONS = ("deriv", "abs", "norm", "transpose") +"""Functions that expressions can use.""" + +_OPERATORS = {"+": operator.add, "-": operator.sub, "*": operator.mul, + "/": operator.truediv} +_PRECEDENCE = {"+": 1, "-": 1, "*": 2, "/": 2} +_NUMBER = re.compile(r"(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?") +_NAME = re.compile(r"[a-z][a-z0-9_]*", re.ASCII) +_CLOSING = {"[": "]", "(": ")"} + +Value = Union[SignalData, float] + + +# --------------------------------------------------------------------------------- +# Parsing +# --------------------------------------------------------------------------------- + + +class _Node: + def canonical(self) -> str: + raise NotImplementedError + + def label(self) -> str: + raise NotImplementedError + + def paths(self) -> List[SignalPath]: + return [] + + def evaluate(self, read: Callable[[SignalPath], SignalData]) -> Value: + raise NotImplementedError + + +class _Path(_Node): + def __init__(self, path: SignalPath): + self.path = path + + def canonical(self) -> str: + return self.path.canonical() + + def label(self) -> str: + return self.path.label + + def paths(self) -> List[SignalPath]: + return [self.path] + + def evaluate(self, read): + return read(self.path) + + +class _Number(_Node): + def __init__(self, value: float, text: str): + self.value, self.text = value, text + + def canonical(self) -> str: + return repr(self.value) + + def label(self) -> str: + return self.text + + def evaluate(self, read): + return self.value + + +class _Negate(_Node): + def __init__(self, operand: _Node): + self.operand = operand + + def canonical(self) -> str: + return f"(-{self.operand.canonical()})" + + def label(self) -> str: + inner = self.operand.label() + return f"-({inner})" if isinstance(self.operand, _Binary) else f"-{inner}" + + def paths(self): + return self.operand.paths() + + def evaluate(self, read): + value = self.operand.evaluate(read) + return -value if isinstance(value, float) else _elementwise( + value, np.negative, Axis(f"-{value.value_axis.label}", value.value_axis.units)) + + +class _Binary(_Node): + def __init__(self, op: str, left: _Node, right: _Node): + self.op, self.left, self.right = op, left, right + + def canonical(self) -> str: + return f"({self.left.canonical()} {self.op} {self.right.canonical()})" + + def label(self) -> str: + def side(node, right): + text = node.label() + if isinstance(node, _Binary) and ( + _PRECEDENCE[node.op] < _PRECEDENCE[self.op] + or (right and _PRECEDENCE[node.op] == _PRECEDENCE[self.op] + and self.op in "-/")): + return f"({text})" + return text + + return f"{side(self.left, False)} {self.op} {side(self.right, True)}" + + def paths(self): + return self.left.paths() + self.right.paths() + + def evaluate(self, read): + return combine(self.op, self.left.evaluate(read), self.right.evaluate(read)) + + +class _Call(_Node): + def __init__(self, function: str, argument: _Node): + self.function, self.argument = function, argument + + def canonical(self) -> str: + return f"{self.function}({self.argument.canonical()})" + + def label(self) -> str: + return f"{self.function}({self.argument.label()})" + + def paths(self): + return self.argument.paths() + + def evaluate(self, read): + value = self.argument.evaluate(read) + if isinstance(value, float): + raise PathError(f"{self.function}() applies to a signal, not to a number.") + if self.function == "deriv": + return derivative(value) + if self.function == "abs": + return _elementwise(value, np.abs, Axis(f"abs({value.value_axis.label})", + value.value_axis.units)) + return normalize(value) # norm + + +class Expression: + """A signal computed from signal paths and numbers (see the module + documentation). It stands where a :class:`SignalPath` does: evaluated by a data + source, it gives :class:`SignalData` whose ``path`` is the expression.""" + + def __init__(self, text: str, node: _Node): + self.text = text + self._node = node + + def canonical(self) -> str: + return self._node.canonical() + + @property + def label(self) -> str: + """A label for legends, from the labels of its paths.""" + return self._node.label() + + @property + def name(self) -> str: + return self.label + + def paths(self) -> List[SignalPath]: + """The signal paths it reads.""" + return self._node.paths() + + def evaluate(self, read: Callable[[SignalPath], SignalData]) -> SignalData: + """The signal, reading the signal of each path with ``read``.""" + value = self._node.evaluate(read) + if isinstance(value, float): + raise PathError(f"'{self.text}' is a number, not a signal.") + return _renamed(value, self) + + def __str__(self) -> str: + return self.text + + def __repr__(self) -> str: + return f"Expression({self.text!r})" + + +def parse_signal(text: str) -> Union[SignalPath, Expression]: + """Parse a signal: a path (as :func:`parse_signal_path`) or an expression. + + Raises: + PathError: when it is neither. + """ + stripped = text.strip() + try: + node = _Parser(stripped).parse() + except PathError as expression_error: + try: # e.g. a path written with spaces around its "/" + return parse_signal_path(stripped) + except PathError as path_error: + chosen = expression_error if _looks_like_expression(stripped) else path_error + raise chosen from None + if isinstance(node, _Path): + return node.path + return Expression(stripped, node) + + +def _looks_like_expression(text: str) -> bool: + """Whether a text that is neither a path nor an expression was meant as an + expression: it starts with a function, or has an operator outside the indices + of its paths.""" + outside = re.sub(r"\[[^\[\]]*\]", "", text) + return bool(re.match(r"\s*[a-z]\w*\s*\(", text) + or re.search(r"[+*]|(?:^|[\s(])-|\s/\s", outside)) + + +class _Parser: + """Recursive descent: sum = product (('+' | '-') product)*, product = unary + (('*' | '/') unary)*, unary = '-' unary | atom, atom = number | path | + function '(' sum ')' | '(' sum ')'.""" + + def __init__(self, text: str): + self.text = text + self.pos = 0 + self.path_end = -1 # where the last path ends: a "/" right after it is its own + + def error(self, message: str) -> PathError: + return PathError(f"Invalid expression '{self.text}': {message}.") + + def parse(self) -> _Node: + if not self.text: + raise self.error("it is empty") + node = self.sum() + self.skip() + if self.pos != len(self.text): + raise self.error(f"unexpected '{self.text[self.pos:]}'") + return node + + def skip(self) -> None: + while self.pos < len(self.text) and self.text[self.pos].isspace(): + self.pos += 1 + + def peek(self) -> str: + self.skip() + return self.text[self.pos] if self.pos < len(self.text) else "" + + def sum(self) -> _Node: + node = self.product() + while self.peek() in ("+", "-"): + op = self.text[self.pos] + self.pos += 1 + node = _Binary(op, node, self.product()) + return node + + def product(self) -> _Node: + node = self.unary() + while self.peek() == "*" or (self.peek() == "/" and self.pos != self.path_end): + op = self.text[self.pos] + self.pos += 1 + node = _Binary(op, node, self.unary()) + return node + + def unary(self) -> _Node: + char = self.peek() + if char == "-": + self.pos += 1 + return _Negate(self.unary()) + if char == "+": + self.pos += 1 + return self.unary() + return self.atom() + + def atom(self) -> _Node: + char = self.peek() + if char == "(": + self.pos += 1 + node = self.sum() + self.expect(")") + return node + number = _NUMBER.match(self.text, self.pos) + if number: + self.pos = number.end() + return _Number(float(number.group()), number.group()) + name = _NAME.match(self.text, self.pos + (char == "#")) + if not name: + found = f"'{self.text[self.pos:]}'" if char else "the end" + raise self.error(f"expected a signal path, a number or '(' at {found}") + after = name.end() + while after < len(self.text) and self.text[after].isspace(): + after += 1 + if after < len(self.text) and self.text[after] == "(" and char != "#": + return self.call(name.group(), after) + return self.path() + + def expect(self, char: str) -> None: + if self.peek() != char: + found = f"'{self.text[self.pos:]}'" if self.pos < len(self.text) else "the end" + raise self.error(f"expected '{char}' at {found}") + self.pos += 1 + + def call(self, function: str, opening: int) -> _Node: + if function not in FUNCTIONS: + raise self.error(f"unknown function '{function}' (available: " + f"{', '.join(FUNCTIONS)})") + start = self.pos + closing = self.closing(opening) + if function in ("norm", "transpose"): + # Of a path, a function of the path itself (as transpose has to be) + arguments = _split_arguments(self.text[opening + 1: closing - 1]) + argument = _Parser(arguments[0].strip()).parse() + if isinstance(argument, _Path): + self.pos = closing + return _Path(parse_signal_path(self.text[start:closing])) + if function == "transpose": + raise self.error("transpose() takes a signal path, not an expression") + if len(arguments) > 1: + raise self.error("norm() takes one argument") + self.pos = closing + return _Call(function, argument) + self.pos = opening + 1 + argument = self.sum() + self.expect(")") + return _Call(function, argument) + + def path(self) -> _Node: + """A path: the IDS name, its occurrence, and nodes joined by '/' or '.' with + no space, each with an index right after its name.""" + start = self.pos + pos = self.pos + (self.text[self.pos] == "#") + pos = _NAME.match(self.text, pos).end() + occurrence = re.compile(r":\d+").match(self.text, pos) + if occurrence: + pos = occurrence.end() + while pos < len(self.text) and self.text[pos] in "/.": + name = _NAME.match(self.text, pos + 1) + if not name: + break + pos = name.end() + if pos < len(self.text) and self.text[pos] in _CLOSING: + pos = self.closing(pos) + self.pos = self.path_end = pos + return _Path(parse_signal_path(self.text[start:pos])) + + def closing(self, opening: int) -> int: + """The position after the bracket that closes the one at ``opening``.""" + depth = [] + for pos in range(opening, len(self.text)): + char = self.text[pos] + if char in _CLOSING: + depth.append(_CLOSING[char]) + elif char in ")]": + if not depth or depth.pop() != char: + break + if not depth: + return pos + 1 + raise self.error("unbalanced brackets") + + +# --------------------------------------------------------------------------------- +# Arithmetic on signals +# --------------------------------------------------------------------------------- + + +def combine(op: str, a: Value, b: Value, label: Optional[str] = None) -> Value: + """``a op b`` for signals and numbers (``op`` is one of ``+ - * /``). + + The result is on the times and the grid of ``a`` (of ``b`` when ``a`` is a + number or does not depend on time and ``b`` does); its values are labeled + ``label`` (by default, the labels of ``a`` and ``b`` with the operator). + + Raises: + PathError: when the two cannot be combined (shapes, maps with profiles, + profiles against different coordinates, ...). + """ + function = _OPERATORS[op] + if isinstance(a, float) and isinstance(b, float): + return function(a, b) + for value in (a, b): + if isinstance(value, Shapes): + raise PathError(f"'{value.path.text}' is a shape; shapes cannot be computed " + "with.") + axis = Axis(label or f"{_label(a)} {op} {_label(b)}", _units(op, a, b)) + if isinstance(a, Channels) or isinstance(b, Channels): + return _members(op, a, b, axis) + if isinstance(a, float) or isinstance(b, float): + number = a if isinstance(a, float) else b + data = b if isinstance(a, float) else a + + def with_number(values): + with np.errstate(divide="ignore", invalid="ignore"): + return function(number, values) if data is b else function(values, number) + + return _elementwise(data, with_number, axis) + kinds = (type(a), type(b)) + if kinds == (Trace, Trace): + return _traces(function, a, b, axis) + if isinstance(a, MapSeries) or isinstance(b, MapSeries): + return _maps(function, a, b, axis) + return _profiles(function, a, b, axis) + + +def derivative(data: SignalData) -> SignalData: + """d/dt of a time trace, d/dx of every profile (member by member).""" + label = f"deriv({data.value_axis.label})" + if isinstance(data, Channels): + members = [derivative(member) for member in data.members] + return Channels(data.path, members, data.names, members[0].value_axis, + data.documentation, data.axis_label) + x_units = data.x_axis.units if hasattr(data, "x_axis") else "" + axis = Axis(label, _ratio(data.value_axis.units, x_units)) + if isinstance(data, Trace): + if data.time is None or data.values.size < 2: + raise PathError(f"deriv() needs a time trace; '{data.path.text}' has " + "one value.") + return Trace(data.path, data.time, _gradient(data.values, data.time), axis, + data.documentation, data.fixed_time) + if isinstance(data, ProfileSeries): + ys = [_gradient(y, x) for x, y in zip(data.xs, data.ys)] + return ProfileSeries(data.path, data.time, data.xs, ys, data.x_axis, axis, + data.documentation, data.fixed_time) + raise PathError(f"deriv() applies to time traces and profiles, but " + f"'{data.path.text}' is a {data.kind}.") + + +def _gradient(y: np.ndarray, x: np.ndarray) -> np.ndarray: + """dy/dx along x (in the order of x), NaN where it cannot be told.""" + result = np.full(len(y), np.nan) + valid = np.isfinite(x) & np.isfinite(y) + if valid.sum() < 2: + return result + order = np.flatnonzero(valid)[np.argsort(x[valid], kind="stable")] + xs, ys = x[order], y[order] + keep = np.concatenate([[True], np.diff(xs) > 0]) # the same x twice: once + order, xs, ys = order[keep], xs[keep], ys[keep] + if len(xs) >= 2: + result[order] = np.gradient(ys, xs) + return result + + +def _label(value: Value) -> str: + if isinstance(value, float): + return f"{value:g}" + label = value.value_axis.label + return f"({label})" if " " in label and not label.startswith(("abs(", "deriv(")) \ + else label + + +def _unit(units: str) -> str: + """Units, empty for dimensionless ones ("1" and "-" in the Data Dictionary).""" + return "" if units in ("1", "-") else units + + +def _ratio(numerator: str, denominator: str) -> str: + numerator, denominator = _unit(numerator), _unit(denominator) + if not denominator: + return numerator + if not numerator: + return f"{denominator}^-1" + return f"{numerator}/{denominator}" + + +def _units(op: str, a: Value, b: Value) -> str: + ua = "" if isinstance(a, float) else _unit(a.value_axis.units) + ub = "" if isinstance(b, float) else _unit(b.value_axis.units) + if op in "+-": + if isinstance(b, float) or ua == ub: + return ua + return ub if isinstance(a, float) else "" + if op == "*": + return ".".join(unit for unit in (ua, ub) if unit) + if isinstance(b, float): + return ua + if ua == ub: + return "" + return _ratio(ua, ub) + + +def _elementwise(data: SignalData, function, axis: Axis) -> SignalData: + """``function`` applied to the values of every slice (and member).""" + if isinstance(data, Channels): + members = [_elementwise(member, function, axis) for member in data.members] + return Channels(data.path, members, data.names, axis, data.documentation, + data.axis_label) + if isinstance(data, Trace): + return Trace(data.path, data.time, function(data.values), axis, + data.documentation, data.fixed_time) + if isinstance(data, ProfileSeries): + return ProfileSeries(data.path, data.time, data.xs, [function(y) for y in data.ys], + data.x_axis, axis, data.documentation, data.fixed_time) + if isinstance(data, MapSeries): + def load(index): + x1, x2, z = data.at_index(index) + return x1, x2, function(z) + + return MapSeries(data.path, data.time, data.count, load, data.x_axis, + data.x2_axis, axis, data.documentation, data.fixed_time) + raise PathError(f"'{data.path.text}' is a {data.kind}; it cannot be computed with.") + + +def _members(op: str, a: Value, b: Value, axis: Axis) -> Channels: + """Member by member: by name, else in order when there are as many.""" + template = a if isinstance(a, Channels) else b + names = template.names + if isinstance(a, Channels) and isinstance(b, Channels): + if a.names == b.names or (len(a.names) == len(b.names) + and not set(a.names) & set(b.names)): + pairs = list(zip(a.members, b.members)) # the same, or all renamed + else: # those of the first that the second has too + others = dict(zip(b.names, b.members)) + names = [name for name in a.names if name in others] + if not names: + raise PathError( + f"'{a.path.text}' and '{b.path.text}' have no element of the same " + "name: they cannot be combined one by one.") + pairs = [(member, others[name]) for member, name in zip(a.members, a.names) + if name in others] + elif isinstance(a, Channels): + pairs = [(member, b) for member in a.members] + else: + pairs = [(a, member) for member in b.members] + members = [combine(op, x, y, axis.label) for x, y in pairs] + return Channels(template.path, members, names, axis, template.documentation, + template.axis_label) + + +def _traces(function, a: Trace, b: Trace, axis: Axis) -> Trace: + base = a if a.time is not None or b.time is None else b + with np.errstate(divide="ignore", invalid="ignore"): + values = function(_at_times(a, base.time), _at_times(b, base.time)) + return Trace(base.path, base.time, values, axis, base.documentation, + base.fixed_time) + + +def _at_times(trace: Trace, times: Optional[np.ndarray]) -> np.ndarray: + """The values of a trace at ``times`` (NaN outside of its time range; a lone + sample only at its time).""" + values = np.atleast_1d(trace.values) + if trace.time is None: + return values[0] if times is None else np.full(len(times), values[0]) + if times is None: + return values + times = np.asarray(times, dtype=float) + if len(trace.time) == len(times) and np.array_equal(trace.time, times): + return values + if len(trace.time) == 1: + return np.where(_same_time(times, trace.time[0]), values[0], np.nan) + return _interp(times, trace.time, values) + + +def _same_time(t, reference: float): + return np.isclose(t, reference, rtol=1e-9, atol=1e-12) + + +def _base(a: SignalData, b: SignalData, kind) -> tuple: + """The signal of type ``kind`` whose grid the result is on, and the times of + the result: those of the first signal, unless only the second depends on time + (e.g. a static profile times a time trace).""" + base = a if isinstance(a, kind) else b + other = b if base is a else a + if base.time is None and other.time is not None: + if isinstance(other, kind): + base = other + return base, other.time + return base, base.time + + +def _slice_near(data: SignalData, t: Optional[float]) -> Optional[int]: + """The slice of ``data`` closest to ``t``, None when ``t`` is out of its time + range (by more than half a step; a single slice only at its time).""" + if data.time is None or t is None: + return 0 + index = data.time_index(t) + if index is None: + return None + times = np.sort(data.time[np.isfinite(data.time)]) + if len(times) == 1: + return index if _same_time(t, times[0]) else None + step = float(np.median(np.diff(times))) + if t < times[0] - step / 2 or t > times[-1] + step / 2: + return None + return index + + +def _profiles(function, a: SignalData, b: SignalData, axis: Axis) -> ProfileSeries: + """Profiles with profiles, or with the value of a trace at the time of each.""" + if isinstance(a, ProfileSeries) and isinstance(b, ProfileSeries) and \ + a.x_axis.label != b.x_axis.label: + raise PathError( + f"'{a.path.text}' and '{b.path.text}' are profiles against " + f"{a.x_axis.label or 'their index'} and {b.x_axis.label or 'their index'}; " + "they cannot be combined.") + base, time = _base(a, b, ProfileSeries) + times = [None] if time is None else [float(t) for t in time] + xs, ys = [], [] + for t in times: + index = 0 if base.time is None else _slice_near(base, t) + x = base.xs[index] if index is not None and index < len(base.xs) else np.empty(0) + with np.errstate(divide="ignore", invalid="ignore"): + ys.append(function(_profile_at(a, t, x), _profile_at(b, t, x))) + xs.append(x) + return ProfileSeries(base.path, time, xs, ys, base.x_axis, axis, + base.documentation, base.fixed_time) + + +def _profile_at(data: SignalData, t: Optional[float], x: np.ndarray) -> np.ndarray: + """A profile (or the value of a trace) at time ``t``, on the grid ``x``.""" + if isinstance(data, Trace): + value = _at_times(data, None if t is None else np.array([t])) + return np.full(len(x), float(np.atleast_1d(value)[0])) + index = _slice_near(data, t) + if index is None or index >= len(data.ys): + return np.full(len(x), np.nan) + own_x, y = data.xs[index], data.ys[index] + if len(own_x) == len(x) and np.array_equal(own_x, x): + return y + return _interp(x, own_x, y) + + +def _maps(function, a: SignalData, b: SignalData, axis: Axis) -> MapSeries: + """Maps with maps (on the grid of the first), or with the value of a trace at + the time of each.""" + for data in (a, b): + if not isinstance(data, (MapSeries, Trace)): + raise PathError(f"'{a.path.text}' and '{b.path.text}' cannot be combined: " + "a map combines with maps, numbers and time traces.") + base, time = _base(a, b, MapSeries) + + def load(index): + t = None if time is None else float(time[index]) + x1, x2, _ = base.at_index(0 if base.time is None else _slice_near(base, t)) + with np.errstate(divide="ignore", invalid="ignore"): + return x1, x2, function(_map_at(a, t, x1, x2), _map_at(b, t, x1, x2)) + + count = base.count if time is base.time else len(time) + return MapSeries(base.path, time, count, load, base.x_axis, base.x2_axis, axis, + base.documentation, base.fixed_time) + + +def _map_at(data: SignalData, t, x1: np.ndarray, x2: np.ndarray) -> np.ndarray: + if isinstance(data, Trace): + value = _at_times(data, None if t is None else np.array([t])) + return np.full((len(x1), len(x2)), float(np.atleast_1d(value)[0])) + index = _slice_near(data, t) + if index is None: + return np.full((len(x1), len(x2)), np.nan) + own1, own2, z = data.at_index(index) + if own1.shape == x1.shape and own2.shape == x2.shape and \ + np.array_equal(own1, x1) and np.array_equal(own2, x2): + return z + # On the grid of the other: along x2 for every row, then along x1 + rows = np.array([_interp(x2, own2, row) for row in np.asarray(z, float)]) + if rows.size == 0: + return np.full((len(x1), len(x2)), np.nan) + return np.array([_interp(x1, own1, column) for column in rows.T]).T + + +def _renamed(data: SignalData, path) -> SignalData: + """``data`` as the signal of ``path`` (an expression), and its members too.""" + renamed = _shallow_copy(data) + renamed.path = path + if isinstance(data, Channels): + renamed.members = [_renamed(member, path) for member in data.members] + renamed.value_axis = Axis(path.label, data.value_axis.units) + return renamed + + +def _shallow_copy(data: SignalData) -> SignalData: + copy = object.__new__(type(data)) + copy.__dict__.update(data.__dict__) + return copy + + +# --------------------------------------------------------------------------------- +# Comparing sources +# --------------------------------------------------------------------------------- + +RELATIVE = {"difference": "-", "ratio": "/"} +"""How a panel can show the data of a source relative to the first one.""" + + +def relative_to(data: SignalData, reference: SignalData, relative: str) -> SignalData: + """The data of a source relative to that of the reference source: their + difference (``data - reference``) or their ratio (``data / reference``), on the + times and grid of ``data``. + + Raises: + PathError: when they cannot be combined (e.g. shapes). + """ + label = data.value_axis.label + return combine(RELATIVE[relative], data, reference, + f"Δ {label}" if relative == "difference" else f"{label} ratio") diff --git a/src/idscope/data/paths.py b/src/idscope/data/paths.py index 22a14e3..bf47531 100644 --- a/src/idscope/data/paths.py +++ b/src/idscope/data/paths.py @@ -43,6 +43,9 @@ its index), as ``numpy.transpose``: ``transpose(path, 1, 0, 2)`` puts the second dimension first. Without axes, their order is reversed (``2, 1, 0`` for three). Dimensions are otherwise never reordered, and time is one like the others. + +Signals computed from paths, such as ``a - b`` or ``deriv(a)``, are parsed and +evaluated by :mod:`idscope.data.expressions`; :func:`parse_alternatives` reads both. """ import dataclasses @@ -169,20 +172,23 @@ def _format_selector(selector: Selector) -> str: ALTERNATIVE_SEPARATOR = "|" -def parse_alternatives(text: str) -> List[SignalPath]: - """Parse alternative paths ``a | b | ...``, to be tried in order. +def parse_alternatives(text: str) -> list: + """Parse alternative signals ``a | b | ...``, to be tried in order: paths, or + expressions of paths (:mod:`idscope.data.expressions`). Alternatives let one configuration work with data that uses different names, e.g. ``core_profiles/profiles_1d/j_phi | core_profiles/profiles_1d/j_tor`` for Data Dictionary 4 and 3. Raises: - PathError: when one of the paths is malformed. + PathError: when one of the alternatives is malformed. """ + from idscope.data.expressions import parse_signal + parts = text.split(ALTERNATIVE_SEPARATOR) if len(parts) > 1 and any(not part.strip() for part in parts): raise PathError(f"Empty alternative in '{text.strip()}'.") - return [parse_signal_path(part) for part in parts] + return [parse_signal(part) for part in parts] def parse_signal_path(text: str) -> SignalPath: diff --git a/src/idscope/data/sources.py b/src/idscope/data/sources.py index 6600f06..fd615c8 100644 --- a/src/idscope/data/sources.py +++ b/src/idscope/data/sources.py @@ -16,6 +16,7 @@ import imas from imas.ids_data_type import IDSDataType +from idscope.data.expressions import Expression from idscope.data.extract import SignalData, apply_x, extract_signal, is_time_aos from idscope.data.paths import SignalPath from idscope.errors import IdscopeError, NoDataError, PathError, SourceError @@ -299,7 +300,9 @@ def signal( cached = self._signals.get(key) if cached is None: try: - if x_path is None: + if x_path is None and isinstance(path, Expression): + cached = path.evaluate(self.signal) # of paths, which are cached + elif x_path is None: ids = self.get_ids(path.ids_name, path.occurrence) cached = extract_signal(ids, path) else: # both are cached: x is often shared by many signals diff --git a/src/idscope/gui/panel.py b/src/idscope/gui/panel.py index cb5b2b1..b36ae83 100644 --- a/src/idscope/gui/panel.py +++ b/src/idscope/gui/panel.py @@ -1,5 +1,6 @@ """A plot panel: one cell of the scope grid (a "wave" in jScope).""" +import dataclasses import html import logging from typing import List, NamedTuple, Optional, Tuple @@ -12,10 +13,12 @@ from idscope.config import PanelConfig, SignalConfig from idscope.data import MAP, PROFILE, SHAPE, TRACE, Channels, MapSeries, ProfileSeries from idscope.data import Shapes, SignalData, parse_alternatives +from idscope.data.expressions import relative_to from idscope.data.extract import ANNULAR, GEOMETRY, MARKERS, POINT from idscope.data.sources import DataSource from idscope.errors import IdscopeError, NoDataError from idscope.gui import style +from idscope.gui.setup_dialog import RELATIVE_CHOICES logger = logging.getLogger(__name__) @@ -182,6 +185,17 @@ def _mark_single_point(item, x, y) -> None: item.setSymbolPen(color) +def _sources_of(entries, kept) -> List[Tuple[int, str]]: + """``(index, label)`` of the sources of the loaded entries whose data is + ``kept``, the label of their data rather than that of their shapes (e.g. the + "B − A" of a source compared with the first).""" + labels = {} + for _, _, source_index, label, data in entries: + if kept(data) and (source_index not in labels or not isinstance(data, Shapes)): + labels[source_index] = label + return sorted(labels.items()) + + def _format(value: float) -> str: return f"{value:.5g}" if np.isfinite(value) else "nan" @@ -271,6 +285,7 @@ def __init__(self, config: PanelConfig, parent=None, member_colors=None): """The kind of a panel of time traces or time-space images whose data is all at one time, and which then shows it as it is at that time (see load).""" self._value_labels: List[object] = [] # values written on the bars + self._compared = False # sources shown relative to the first (config.relative) self._cursor_items = None # the crosshair of the Point mode, once needed self._placing_cursor = False self._labeled_heights: List[float] = [] # heights of the bars labeled @@ -514,6 +529,7 @@ def clear(self) -> None: self._time_line = None self.kind = None self._one_time = None + self._compared = False self.errors = [] self.notes = [] self._image_curves = [] @@ -553,6 +569,8 @@ def load(self, sources: List[DataSource], time: Optional[float]) -> List[str]: errors.append(f"{source.label}: {exc}") continue loaded.append((signal_index, signal, source_index, source.label, data)) + if self.config.relative and len({entry[2] for entry in loaded}) > 1: + loaded = self._relative(loaded, errors) self.kind = self._resolve_kind(loaded) if self.kind in ("trace", "spacetime") and _at_one_time(loaded): @@ -570,18 +588,23 @@ def load(self, sources: List[DataSource], time: Optional[float]) -> List[str]: errors.append(f"{entry[3]}: {problem}") else: accepted.append(entry) - sources_shown = sorted({(entry[2], entry[3]) for entry in accepted}) + # A source compared with the first is named as that ("B − A") rather than + # as itself for its shapes; the line styles stand for the sources of the + # lines (the data of a map is an image) + sources_shown = _sources_of(accepted, lambda data: True) + line_sources = _sources_of(accepted, lambda data: not isinstance(data, MapSeries)) + # Compared sources are named even alone ("B − A") + n_sources = max(len(sources_shown), 2 if self._compared else 1) n_signals = len(self.config.signals) if self._one_time == "trace": self._add_bars_at_one_time(accepted, sources_shown) else: for entry in accepted: - self._add_curves(*entry, n_sources=len(sources_shown), - n_signals=n_signals) + self._add_curves(*entry, n_sources=n_sources, n_signals=n_signals) if len(sources_shown) > 1 and self._styled_sources: # Colors tell the signals (or members) apart: say which line style is # which source - for source_index, source in sources_shown: + for source_index, source in line_sources: proxy = pg.PlotDataItem(pen=style.source_pen(source_index)) self._legend_proxies.append(proxy) self._legend_entries.append((proxy, source)) @@ -612,6 +635,35 @@ def load(self, sources: List[DataSource], time: Optional[float]) -> List[str]: self.apply_ranges() return errors + def _relative(self, loaded, errors: List[str]) -> list: + """Every source but the first as its difference with (or ratio to) the + first, signal by signal (see :attr:`PanelConfig.relative`). The first one is + then not shown; shapes, which cannot be compared, stay as they are.""" + relative = self.config.relative + symbol = "−" if relative == "difference" else "/" + compared = [] + for signal_index in sorted({entry[0] for entry in loaded}): + entries = [entry for entry in loaded if entry[0] == signal_index] + compared += [entry for entry in entries if isinstance(entry[4], Shapes)] + data = [entry for entry in entries if not isinstance(entry[4], Shapes)] + if len(data) < 2: + continue # nothing to compare with: the others do not have it + reference = data[0] + for _, signal, source_index, label, values in data[1:]: + try: + values = relative_to(values, reference[4], relative) + except IdscopeError as exc: + errors.append(f"{label}: {exc}") + continue + # The offset of the signal cancels out, and a ratio has no scale + signal = dataclasses.replace( + signal, offset=0.0, scale=signal.scale if relative == "difference" + else 1.0) + compared.append((signal_index, signal, source_index, + f"{label} {symbol} {reference[3]}", values)) + self._compared = True + return compared + def _resolve_kind(self, loaded) -> Optional[str]: if self.config.kind != "auto": return self.config.kind @@ -815,7 +867,7 @@ def _add_bars_at_one_time(self, accepted, sources_shown) -> None: voltages of the coils), the bars have a color per signal and source instead, as those of signals with wildcards otherwise do. """ - n_sources = len(sources_shown) + n_sources = max(len(sources_shown), 2 if self._compared else 1) rank = {index: k for k, (index, _) in enumerate(sources_shown)} bars, signals_of = [], {} for signal_index, signal, source_index, source, data in accepted: @@ -1469,6 +1521,31 @@ def _install_menu(self) -> None: menu.addSeparator() else: menu.addAction(text).triggered.connect(slot) + # The sources as they are, or relative to the first + compare = menu.addMenu("Sources") + self._relative_actions = [] + for text, value in RELATIVE_CHOICES: + action = compare.addAction(text) + action.setCheckable(True) + action.triggered.connect(lambda checked=False, value=value: + self.set_relative(value)) + self._relative_actions.append((action, value)) + compare.aboutToShow.connect(self._check_relative) + self._check_relative() + + def _check_relative(self) -> None: + for action, value in self._relative_actions: + action.setChecked(value == self.config.relative) + + def set_relative(self, relative: Optional[str]) -> None: + """Show the sources as they are (``None``), or every one but the first as its + ``difference`` with the first or its ``ratio`` to it.""" + self._check_relative() + if relative == self.config.relative: + return + self.config.relative = relative + self._check_relative() + self.signalsChanged.emit(self, "compare sources") # ------------------------------------------------------------- drag & drop diff --git a/src/idscope/gui/setup_dialog.py b/src/idscope/gui/setup_dialog.py index dcc3788..403b6d8 100644 --- a/src/idscope/gui/setup_dialog.py +++ b/src/idscope/gui/setup_dialog.py @@ -11,6 +11,12 @@ COLUMNS = ["Path", "Label", "Scale", "Offset", "X path", "Color", "Sources"] _PATH, _LABEL, _SCALE, _OFFSET, _X, _COLOR, _SOURCES = range(len(COLUMNS)) +RELATIVE_CHOICES = [ + ("as they are", None), + ("minus the first (B − A)", "difference"), + ("over the first (B / A)", "ratio"), +] +"""How the panel shows the sources: the texts of the choices and their values.""" class SetupDialog(QtWidgets.QDialog): @@ -59,6 +65,18 @@ def __init__(self, config: PanelConfig, source_labels: List[str], parent=None): "Draw maps and time-space images as this many contour lines instead of an " "image (0: image)" ) + self.relative = QtWidgets.QComboBox() + for text, value in RELATIVE_CHOICES: + self.relative.addItem(text, value) + self.relative.setCurrentIndex( + next((i for i, (_, value) in enumerate(RELATIVE_CHOICES) + if value == config.relative), 0)) + self.relative.setToolTip( + "With several sources, show every source but the first relative to the\n" + "first, signal by signal: their difference or their ratio, on the times\n" + "and the grid of each source (the first is interpolated). Shapes stay as\n" + "they are. The scale still applies to a difference; the offset does not." + ) form.addWidget(QtWidgets.QLabel("Title"), 0, 0) form.addWidget(self.title, 0, 1, 1, 3) @@ -78,6 +96,9 @@ def __init__(self, config: PanelConfig, source_labels: List[str], parent=None): checks.addSpacing(12) checks.addWidget(QtWidgets.QLabel("contours")) checks.addWidget(self.contours) + checks.addSpacing(12) + checks.addWidget(QtWidgets.QLabel("sources")) + checks.addWidget(self.relative) checks.addStretch(1) form.addLayout(checks, 3, 1, 1, 5) layout.addLayout(form) @@ -102,7 +123,10 @@ def __init__(self, config: PanelConfig, source_labels: List[str], parent=None): "Alternatives, tried in order: path_a | path_b\n" "norm(path): profiles normalized from their first to last point\n" "transpose(path, 1, 0, 2): the dimensions in another order\n" - "(numpy style; reversed without axes)", + "(numpy style; reversed without axes)\n" + "Arithmetic: path_a - path_b, path_a / path_b (spaces around /),\n" + "2 * path, deriv(path), abs(path); the second signal is taken on\n" + "the times and the grid of the first", _SCALE: "plotted value = value × scale + offset", _OFFSET: "plotted value = value × scale + offset", _X: "x values of a profile (default: its coordinate), " @@ -276,6 +300,7 @@ def read_config(self) -> PanelConfig: equal_aspect=self.equal_aspect.isChecked(), contours=self.contours.value(), height=self.config.height, + relative=self.relative.currentData(), ) def apply(self) -> bool: diff --git a/tests/test_config.py b/tests/test_config.py index 32a9945..4ae2e86 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -106,6 +106,8 @@ def test_default_values_are_omitted(): ({"time": True}, "time"), ({"columns": [{"panels": [{"contours": -1}]}]}, "must not be negative"), ({"columns": [{"panels": [{"contours": 2.5}]}]}, "contours"), + ({"columns": [{"panels": [{"relative": "minus"}]}]}, + "'relative' must be one of difference, ratio"), ], ) def test_invalid(data, message): @@ -144,3 +146,18 @@ def test_load_errors(tmp_path): def test_default_layout_is_a_single_panel(): config = default_config() assert len(config.columns) == 1 and len(config.panels()) == 1 + + +def test_sources_relative_to_the_first(tmp_path): + config = config_from_dict({"columns": [{"panels": [ + {"relative": "difference", "signals": ["summary/global_quantities/ip/value - " + "equilibrium/time_slice/global_quantities/ip"]}, + {"relative": "ratio"}, {}]}]}) + assert [panel.relative for panel in config.panels()] == ["difference", "ratio", None] + path = tmp_path / "relative.yaml" + save_config(config, path) + assert "relative: difference" in path.read_text() + again = load_config(path) + assert [panel.relative for panel in again.panels()] == ["difference", "ratio", None] + assert again.panels()[0].signals[0].path.endswith(" - equilibrium/time_slice/" + "global_quantities/ip") diff --git a/tests/test_expressions.py b/tests/test_expressions.py new file mode 100644 index 0000000..eee7f76 --- /dev/null +++ b/tests/test_expressions.py @@ -0,0 +1,279 @@ +import numpy as np +import pytest + +from idscope.data import ( + Axis, + Channels, + MapSeries, + ProfileSeries, + SignalPath, + Trace, + parse_alternatives, + parse_signal_path, +) +from idscope.data.expressions import ( + Expression, + combine, + derivative, + parse_signal, + relative_to, +) +from idscope.errors import NoDataError, PathError + +IP = "equilibrium/time_slice/global_quantities/ip" +SUMMARY_IP = "summary/global_quantities/ip/value" +TE = "core_profiles/profiles_1d/electrons/temperature" +TI = "core_profiles/profiles_1d/t_i_average" +TIME = np.linspace(0.0, 10.0, 41) +RAMP = np.minimum(TIME / 2.0, 1.0) + + +def value(source, text): + (signal,) = parse_alternatives(text) + return source.signal(signal) + + +# --------------------------------------------------------------------------------- +# Parsing +# --------------------------------------------------------------------------------- + + +@pytest.mark.parametrize("text", [ + IP, "magnetics/flux_loop(2)/flux/data", f"{TE}[-1]", "equilibrium:1/time_slice/q", + "norm(equilibrium/time_slice/profiles_1d/psi)", + "transpose(equilibrium/time_slice/profiles_2d[0]/psi, 1, 0)", + "core_profiles/profiles_1d/ electrons/temperature", # spaces, as before +]) +def test_paths_stay_paths(text): + parsed = parse_signal(text) + assert isinstance(parsed, SignalPath) + assert parsed == parse_signal_path(text) + + +def test_expressions(): + difference = parse_signal(f"{IP} - {SUMMARY_IP}") + assert isinstance(difference, Expression) + assert difference.label == "global_quantities/ip - ip/value" + assert [path.text for path in difference.paths()] == [IP, SUMMARY_IP] + # "/" between nodes is a path; with spaces around it, a division + assert isinstance(parse_signal("a/b"), SignalPath) + ratio = parse_signal(f"{TI} / {TE}") + assert ratio.label == "profiles_1d/t_i_average / electrons/temperature" + # precedence, parentheses and unary minus, in the labels + assert parse_signal(f"{TE} - {TI} * 2").label == \ + "electrons/temperature - profiles_1d/t_i_average * 2" + assert parse_signal(f"({TE} - {TI}) / {TE}").label == \ + "(electrons/temperature - profiles_1d/t_i_average) / electrons/temperature" + assert parse_signal(f"{TE} - ({TI} - {TE})").label == \ + "electrons/temperature - (profiles_1d/t_i_average - electrons/temperature)" + assert parse_signal(f"-({IP} - {SUMMARY_IP})").label == \ + "-(global_quantities/ip - ip/value)" + assert parse_signal(f"deriv({IP}) * 1e-6").label == "deriv(global_quantities/ip) * 1e-6" + # norm() of a path is still a path; of an expression, an expression + assert isinstance(parse_signal(f"norm({TE}) * 2"), Expression) + assert isinstance(parse_signal(f"norm({TE} / {TI})"), Expression) + # an index with a minus sign is not a subtraction + assert parse_signal(f"{TE}[-1] + 1").label == "electrons/temperature[-1] + 1" + # the same expression, however written, is read once by a source + assert parse_signal(f"{IP}-{SUMMARY_IP}").canonical() == difference.canonical() + # alternatives of expressions + first, second = parse_alternatives(f"{IP} * 2 | {SUMMARY_IP}") + assert isinstance(first, Expression) and isinstance(second, SignalPath) + + +@pytest.mark.parametrize("text, message", [ + (f"{TE} +", "expected a signal path, a number or '\\(' at the end"), + (f"deriv({TE}", "unbalanced brackets"), + (f"sqrt({TE})", "unknown function 'sqrt' \\(available: deriv, abs, norm, transpose\\)"), + (f"transpose({TE} * 2)", "transpose\\(\\) takes a signal path, not an expression"), + (f"{TE} ** 2", "expected a signal path"), + (f"{TE} - core_profiles", "only names the IDS"), +]) +def test_invalid_expressions(text, message): + with pytest.raises(PathError, match=message): + parse_signal(text) + + +def test_a_number_is_not_a_signal(source_a): + with pytest.raises(PathError, match="is a number, not a signal"): + value(source_a, "2 * 3") + + +# --------------------------------------------------------------------------------- +# Evaluation by a data source +# --------------------------------------------------------------------------------- + + +def test_traces(source_a): + difference = value(source_a, f"{IP} - {SUMMARY_IP}") + assert isinstance(difference, Trace) + np.testing.assert_allclose(difference.values, -2.0e6 * RAMP) + assert difference.value_axis == Axis("global_quantities/ip - ip/value", "A") + assert difference.path.text == f"{IP} - {SUMMARY_IP}" + scaled = value(source_a, f"-{IP} * 1e-6") + np.testing.assert_allclose(scaled.values, RAMP) + assert scaled.value_axis.units == "A" # scaled by a number: the same quantity + rate = value(source_a, f"deriv({SUMMARY_IP})") + assert rate.value_axis == Axis("deriv(ip/value)", "A/s") + np.testing.assert_allclose(rate.values[:7], 0.5e6) + assert value(source_a, f"abs({IP})").values.min() >= 0 + ratio = value(source_a, f"{SUMMARY_IP} / {SUMMARY_IP}") + assert ratio.value_axis.units == "" # dimensionless + product = value(source_a, f"{SUMMARY_IP} * {IP}") + assert product.value_axis.units == "A.A" + # read once, then from the cache of the source + assert value(source_a, f"{IP}-{SUMMARY_IP}") is difference + + +def test_profiles(source_a): + ratio = value(source_a, f"{TI} / {TE}") + assert isinstance(ratio, ProfileSeries) and len(ratio.ys) == 41 + np.testing.assert_allclose(ratio.ys[20][:10], 0.8) + assert ratio.value_axis.units == "" + # divided by a time trace: the central value, at the time of each profile + peaked = value(source_a, f"{TE} / {TE}[0]") + np.testing.assert_allclose([y[0] for y in peaked.ys], 1.0) + slope = value(source_a, f"deriv({TE})") + x, y, _ = slope.at(5.0) + _, te, _ = value(source_a, TE).at(5.0) + np.testing.assert_allclose(y, np.gradient(te, x)) + assert slope.value_axis.units == "eV" # rho_tor_norm has no units + normalized = value(source_a, f"norm({TE} * 2)") + np.testing.assert_allclose([y[0] for y in normalized.ys], 0.0) + with pytest.raises(PathError, match="profiles against rho_tor_norm and psi"): + value(source_a, f"{TE} - equilibrium/time_slice/profiles_1d/q") + + +def test_maps_and_members(source_a): + zero = value(source_a, "equilibrium/time_slice/profiles_2d[0]/psi * 2 - " + "equilibrium/time_slice/profiles_2d[0]/psi * 2") + assert isinstance(zero, MapSeries) and zero.count == 41 + x1, x2, z, _ = zero.at(5.0) + assert z.shape == (len(x1), len(x2)) and np.all(z == 0) + currents = value(source_a, "abs(pf_active/coil[:]/current/data) * 1e-3") + assert isinstance(currents, Channels) and currents.names[6] == "PF1" + np.testing.assert_allclose(currents.members[6].values, 4.0 * RAMP) + # member by member, and with one signal for all + zero = value(source_a, "pf_active/coil[:]/current/data - pf_active/coil[:]/current/data") + assert all(np.all(member.values == 0) for member in zero.members) + minus_ip = value(source_a, f"pf_active/coil[:]/current/data - {SUMMARY_IP}") + np.testing.assert_allclose(minus_ip.members[6].values, (-4.0e3 - 1.0e6) * RAMP) + with pytest.raises(PathError, match="applies to time traces and profiles"): + value(source_a, "deriv(equilibrium/time_slice/profiles_2d[0]/psi)") + with pytest.raises(PathError, match="shapes cannot be computed with"): + value(source_a, "wall/description_2d[0]/limiter/unit[:]/outline * 2") + + +def test_what_a_source_does_not_have(source_a): + """Not filled in one of the signals: not filled (a note, not an error).""" + with pytest.raises(NoDataError): + value(source_a, f"{IP} - summary/global_quantities/beta_pol/value") + # the first alternative that can be read + data = source_a.first_signal(parse_alternatives( + f"{IP} - summary/global_quantities/beta_pol/value | {IP} * 2")) + np.testing.assert_allclose(data.values, -2.0e6 * RAMP) + + +# --------------------------------------------------------------------------------- +# The arithmetic +# --------------------------------------------------------------------------------- + +A = parse_signal_path("summary/global_quantities/ip/value") + + +def trace(time, values, units="A"): + return Trace(A, None if time is None else np.asarray(time, float), + np.asarray(values, float), Axis("ip", units)) + + +def profiles(time, xs, ys, x_label="rho_tor_norm"): + return ProfileSeries(A, None if time is None else np.asarray(time, float), xs, ys, + Axis(x_label, ""), Axis("te", "eV")) + + +def test_traces_on_the_times_of_the_first(): + a = trace(np.arange(0.0, 11.0), np.arange(0.0, 11.0)) + b = trace(np.arange(0.0, 5.5, 0.5), 2.0 * np.arange(0.0, 5.5, 0.5)) + difference = combine("-", a, b) + np.testing.assert_array_equal(difference.time, a.time) + np.testing.assert_allclose(difference.values[:6], -np.arange(0.0, 6.0)) + assert np.isnan(difference.values[6:]).all() # where the second has no data + # a constant, and a single sample (only at its time) + np.testing.assert_allclose(combine("+", a, trace(None, [1.0])).values, a.values + 1) + alone = combine("-", a, trace([3.0], [1.0])) + assert np.isnan(alone.values[np.arange(11) != 3]).all() and alone.values[3] == 2.0 + # the second without time, the first with: on the times of the first + assert combine("*", trace(None, [2.0]), a).time is a.time + + +def test_profiles_on_the_grid_and_times_of_the_first(): + grid, coarse = np.linspace(0.0, 1.0, 11), np.linspace(0.0, 1.0, 3) + a = profiles([1.0, 2.0, 3.0], [grid] * 3, [grid * 0 + k for k in (1.0, 2.0, 3.0)]) + b = profiles([2.0, 2.8], [coarse] * 2, [coarse, 2 * coarse]) + difference = combine("-", a, b) + assert len(difference.ys) == 3 and difference.xs[0] is grid + # before the second (by more than half its step): none of it + assert np.isnan(difference.ys[0]).all() + np.testing.assert_allclose(difference.ys[1], 2.0 - grid) # interpolated + np.testing.assert_allclose(difference.ys[2], 3.0 - 2 * grid) # its closest slice + # a static profile times a time trace: a profile at each time of the trace + scaled = combine("*", profiles(None, [grid], [grid]), trace([0.0, 1.0], [1.0, 2.0])) + np.testing.assert_array_equal(scaled.time, [0.0, 1.0]) + np.testing.assert_allclose(scaled.ys[1], 2 * grid) + with pytest.raises(PathError, match="against rho_tor_norm and psi_n"): + combine("-", a, profiles(None, [grid], [grid], x_label="psi_n")) + + +def test_maps_on_the_grid_of_the_first(): + x1, x2 = np.linspace(0.0, 1.0, 5), np.linspace(0.0, 2.0, 9) + fine = np.linspace(0.0, 1.0, 11) + plane = MapSeries(A, None, 1, lambda i: (x1, x2, np.add.outer(x1, x2)), + Axis("R", "m"), Axis("Z", "m"), Axis("psi", "Wb")) + other = MapSeries(A, None, 1, lambda i: (fine, x2, np.add.outer(fine, x2)), + Axis("R", "m"), Axis("Z", "m"), Axis("psi", "Wb")) + _, _, z = combine("-", plane, other).at_index(0) + np.testing.assert_allclose(z, 0.0, atol=1e-12) + ramp = trace([0.0, 1.0], [1.0, 3.0]) + scaled = combine("*", plane, ramp) # a static map at the times of a trace + assert scaled.count == 2 + np.testing.assert_allclose(scaled.at_index(1)[2], 3 * np.add.outer(x1, x2)) + with pytest.raises(PathError, match="a map combines with maps, numbers and time"): + combine("-", plane, profiles(None, [x1], [x1])) + + +def test_members_by_name(): + def members(names, values): + data = [trace(None, [v]) for v in values] + return Channels(A, data, names, Axis("current", "A"), axis_label="coil") + + a = members(["PF1", "PF2", "PF3"], [1.0, 2.0, 3.0]) + by_name = combine("-", a, members(["PF3", "PF1"], [30.0, 10.0])) + assert by_name.names == ["PF1", "PF3"] # those both have + np.testing.assert_allclose([float(m.values) for m in by_name.members], [-9.0, -27.0]) + renamed = combine("-", a, members(["P1", "P2", "P3"], [1.0, 1.0, 1.0])) + assert renamed.names == a.names # all renamed: in order + with pytest.raises(PathError, match="no element of the same name"): + combine("-", a, members(["X"], [1.0])) + + +def test_derivatives(): + t = np.array([0.0, 1.0, 1.0, 3.0, np.nan]) + rate = derivative(trace(t, [0.0, 2.0, 2.0, 6.0, 1.0])) + np.testing.assert_allclose(rate.values[[0, 1, 3]], 2.0) # the same time once + assert np.isnan(rate.values[4]) and rate.value_axis.units == "A/s" + with pytest.raises(PathError, match="has one value"): + derivative(trace(None, [1.0])) + x = np.linspace(0.0, 1.0, 21) + slope = derivative(profiles(None, [x[::-1]], [x[::-1] ** 2])) # x decreasing + np.testing.assert_allclose(slope.ys[0][1:-1], 2 * x[::-1][1:-1], rtol=1e-12) + + +def test_relative_to_a_reference(): + a = trace([0.0, 1.0, 2.0], [1.0, 2.0, 4.0]) + b = trace([0.0, 2.0], [2.0, 8.0]) + difference = relative_to(b, a, "difference") + np.testing.assert_allclose(difference.values, [1.0, 4.0]) # on the times of B + assert difference.value_axis == Axis("Δ ip", "A") + ratio = relative_to(b, a, "ratio") + np.testing.assert_allclose(ratio.values, [2.0, 2.0]) + assert ratio.value_axis == Axis("ip ratio", "") diff --git a/tests/test_gui.py b/tests/test_gui.py index 7a5da9b..b4bcf65 100644 --- a/tests/test_gui.py +++ b/tests/test_gui.py @@ -1151,3 +1151,71 @@ def test_coordinates_are_rounded_to_the_view(): # on a log axis, the value, to what the decades shown resolve assert _round_to_view(np.log10(0.39994), 1.8, log=True) == pytest.approx(0.4) assert _round_to_view(np.log10(123456.0), 3.0, log=True) == pytest.approx(123500.0) + + +def test_expressions_in_a_panel(make_window, data_dir): + window = show(make_window, [f"A={data_dir / 'a.nc'}", f"B={data_dir / 'b.nc'}"], { + "signals": ["equilibrium/time_slice/global_quantities/ip - " + "summary/global_quantities/ip/value", "sqrt(summary/global_quantities/ip/value)"] + }) + panel = window.panels[0] + assert [c.name for c in panel.curves] == [ + "global_quantities/ip - ip/value [A]", "global_quantities/ip - ip/value [B]"] + np.testing.assert_allclose(panel.curves[1].item.yData[-1], -1.6e6) + assert panel.errors == ["Invalid expression 'sqrt(summary/global_quantities/ip/value)': " + "unknown function 'sqrt' (available: deriv, abs, norm, " + "transpose)."] + + +def test_sources_relative_to_the_first(window, make_window, data_dir): + """Every source but the first as its difference with the first, or its ratio to it.""" + ip, te, psi = window.panels[0], window.panels[3], window.panels[6] + ip.set_relative("difference") + (curve,) = ip.curves + assert curve.name == "Ip (equilibrium) [B − A]" + # B has 80 % of A: B − A is 0.2 MA, and the scale of the signal (-1e-6) applies + np.testing.assert_allclose(curve.item.yData[-1], -0.2) + # which the crosshair of the Point mode reads + window.set_mode("point") + ip._on_move(ip.plot.getViewBox().mapViewToScene(QtCore.QPointF(5.0, -0.1))) + assert ip._cursor_items[2].textItem.toPlainText().splitlines()[1:] == [ + "■ Ip (equilibrium) [B − A] = -0.2"] + window.set_mode("pan") + # a ratio has neither the scale (1e-3) nor the offset of the signal + te.config.signals[0].offset = 5.0 + te.set_relative("ratio") + (curve,) = te.curves + assert curve.name == "electrons/temperature [B / A]" + np.testing.assert_allclose(curve.item.yData[:10], 0.8, rtol=0.01) + assert te.plot.getAxis("left").labelText == "Te [keV]" # as configured + # a map: the image of B − A, with the boundary of each source as it is + psi.set_relative("difference") + assert [c.name for c in psi.curves] == [ + "profiles_2d[0]/psi [B − A]", "boundary [A]", "boundary [B]"] + assert [label.text for _, label in psi.legend.items][-2:] == ["A", "B"] + assert psi._colorbar.getAxis("left").labelText.startswith("Δ profiles_2d[0]/psi") + # from the menu of the panel, and back, which can be undone + actions = dict((value, action) for action, value in ip._relative_actions) + actions[None].trigger() + assert window.panels[0].config.relative is None + assert [c.name for c in window.panels[0].curves] == [ + "Ip (equilibrium) [A]", "Ip (equilibrium) [B]"] + window.undo() + assert window.panels[0].config.relative == "difference" + # with one source, the data as it is + alone = show(make_window, [f"A={data_dir / 'a.nc'}"], + {"relative": "difference", "signals": ["summary/global_quantities/ip/value"]}) + assert [c.name for c in alone.panels[0].curves] == ["ip/value"] + + +def test_setup_dialog_relative(window): + from idscope.gui.setup_dialog import SetupDialog + + panel = window.panels[3] + panel.config.relative = "ratio" + dialog = SetupDialog(panel.config, ["A", "B"], window) + assert dialog.relative.currentData() == "ratio" + dialog.relative.setCurrentIndex(1) + assert dialog.apply() and panel.config.relative == "difference" + dialog.relative.setCurrentIndex(0) + assert dialog.apply() and panel.config.relative is None