Comparison

Compare the same field between two runs, across codes or across parameter studies.

Cross-model comparison converts to SI at the comparison boundary — different normalizations make code units incomparable — so units defaults to "si". Same-model comparisons with identical normalization can pass units="code", and dimensionless quantities (beta, Mach numbers, entropy) need no conversion at all.

Datasets on different grids must be aligned first; see align_grids in Regridding.

nan_policy is forwarded unchanged to the underlying diagnostics, so masked regions from SphereSelection, FieldDataset.where(), or out-of-domain regrid fills behave the same whether you pass arrays or datasets — see Conventions § NaN handling.

comparison

Grid-aware cross-model comparison utilities.

These three functions are the only place the pure pypic.diagnostics norms are wrapped as a FieldDataset-level public API (pypic.reductions is the other module that runs NumPy reductions over dataset-held arrays). The norms themselves stay pure (NumPy in, NumPy out); the functions here add the glue layer — alignment via pypic.regridding, alias resolution through both datasets, and SI conversion at the comparison boundary — then delegate the actual norm evaluation back to the pure helpers.

Cross-model comparisons (e.g. iPIC3D vs BATSRUS) default to SI because different normalizations are incomparable in code units; same-model runs can opt in to units="code" to skip the conversion.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4,), spacing=(1.0,), origin=(0.0,))
>>> a = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
...     grid, Normalization.identity(),
... )
>>> b = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
...     grid, Normalization.identity(),
... )
>>> float(compare_fields(a, b, "B_1"))
0.0

compare_fields(a, b, field, *, metric='l2', units='si', method='linear', nan_policy='omit', frame=None)

Compute an error norm between one field of two datasets.

Aligns a and b onto their common grid via pypic.regridding.align_grids, resolves field through both datasets' aliases to a shared canonical name, converts to SI (by default) or leaves in code units, and delegates to the pure diagnostic in pypic.diagnostics.

When the two datasets are in different frames, b is transformed to a's frame via FieldDataset.transform_to before alignment. Pass an explicit frame to compare in a third reference frame — both inputs are then transformed to that frame instead. A clear ValueError is raised if any required transform is missing.

Parameters:

Name	Type	Description	Default
a	FieldDataset	Datasets to compare. Grids may differ in resolution, extent, or both — the overlap is computed automatically.	required
b	FieldDataset	Datasets to compare. Grids may differ in resolution, extent, or both — the overlap is computed automatically.	required
field	str	Field name; canonical or alias. Must resolve to the same canonical name in both datasets.	required
metric	('l2', 'linf')	Error norm. "l2" uses l2_relative_error (relative to b); "linf" uses linf_error (absolute max).	"l2"
units	('si', 'code')	"si" converts both fields through FieldDataset.in_si before comparing — the default, safe for cross-model runs. "code" compares raw code-unit values; valid only when both datasets share a normalization.	"si"
method	str	Interpolation method passed through to align_grids (e.g. "linear", "nearest", "cubic"). Default "linear".	'linear'
nan_policy	('omit', 'propagate', 'raise')	Forwarded to the pure diagnostic. Default "omit" masks NaN cells from the metric (with a warning naming the dropped count) — useful for sphere selections, masked regions, and other upstream sources of NaN. Use "propagate" for strict verification where any NaN should poison the result.	"omit"
frame	str | None	Reference frame to compare in. None (default) uses a's frame, transforming b if needed. A non-None value transforms both a and b to that frame first — useful when neither dataset lives natively in the frame you want to plot or report in. Each dataset must have a transform registered to the target (or already be in it).	None

Returns:

Type	Description
float	Scalar error metric.

Raises:

Type	Description
KeyError	If field is missing in either dataset (message from FieldDataset.resolve_key includes close-match suggestions).
ValueError	If metric or units is unknown, if field resolves to different canonical names in the two datasets, or if the grids do not overlap.
NotImplementedError	If either grid is non-Cartesian (propagated from align_grids).

Warns:

Type	Description
UserWarning	If any axis' spacing ratio exceeds 10×.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4,), spacing=(1.0,), origin=(0.0,))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
...     grid, Normalization.identity(),
... )
>>> float(compare_fields(ds, ds, "B_1"))
0.0

Source code in src/pypic/comparison.py

def compare_fields(
    a: FieldDataset,
    b: FieldDataset,
    field: str,
    *,
    metric: str = "l2",
    units: str = "si",
    method: str = "linear",
    nan_policy: NanPolicy = "omit",
    frame: str | None = None,
) -> float:
    r"""Compute an error norm between one field of two datasets.

    Aligns *a* and *b* onto their common grid via
    [`pypic.regridding.align_grids`][pypic.regridding.align_grids],
    resolves *field* through both datasets' aliases to a shared canonical
    name, converts to SI (by default) or leaves in code units, and
    delegates to the pure diagnostic in [`pypic.diagnostics`][pypic.diagnostics].

    When the two datasets are in different frames, *b* is transformed to
    *a*'s frame via `FieldDataset.transform_to` before alignment.
    Pass an explicit *frame* to compare in a third reference frame —
    both inputs are then transformed to that frame instead. A clear
    `ValueError` is raised if any required transform is missing.

    Parameters
    ----------
    a, b : FieldDataset
        Datasets to compare. Grids may differ in resolution, extent,
        or both — the overlap is computed automatically.
    field : str
        Field name; canonical or alias. Must resolve to the same
        canonical name in both datasets.
    metric : {"l2", "linf"}
        Error norm. ``"l2"`` uses
        [`l2_relative_error`][pypic.diagnostics.l2_relative_error]
        (relative to *b*); ``"linf"`` uses [`linf_error`][pypic.diagnostics.linf_error]
        (absolute max).
    units : {"si", "code"}
        ``"si"`` converts both fields through `FieldDataset.in_si`
        before comparing — the default, safe for cross-model runs.
        ``"code"`` compares raw code-unit values; valid only when both
        datasets share a normalization.
    method : str
        Interpolation method passed through to
        [`align_grids`][pypic.regridding.align_grids] (e.g. ``"linear"``,
        ``"nearest"``, ``"cubic"``). Default ``"linear"``.
    nan_policy : {"omit", "propagate", "raise"}
        Forwarded to the pure diagnostic. Default ``"omit"`` masks NaN
        cells from the metric (with a warning naming the dropped count)
        — useful for sphere selections, masked regions, and other
        upstream sources of NaN. Use ``"propagate"`` for strict
        verification where any NaN should poison the result.
    frame : str | None
        Reference frame to compare in. ``None`` (default) uses *a*'s
        frame, transforming *b* if needed. A non-``None`` value
        transforms both *a* and *b* to that frame first — useful when
        neither dataset lives natively in the frame you want to plot or
        report in. Each dataset must have a transform registered to the
        target (or already be in it).

    Returns
    -------
    float
        Scalar error metric.

    Raises
    ------
    KeyError
        If *field* is missing in either dataset (message from
        `FieldDataset.resolve_key` includes close-match suggestions).
    ValueError
        If *metric* or *units* is unknown, if *field* resolves to
        different canonical names in the two datasets, or if the grids
        do not overlap.
    NotImplementedError
        If either grid is non-Cartesian (propagated from `align_grids`).

    Warns
    -----
    UserWarning
        If any axis' spacing ratio exceeds 10×.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4,), spacing=(1.0,), origin=(0.0,))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
    ...     grid, Normalization.identity(),
    ... )
    >>> float(compare_fields(ds, ds, "B_1"))
    0.0
    """
    _validate_choice(metric, _ALLOWED_METRICS, "metric")
    _validate_choice(units, _ALLOWED_UNITS, "units")
    # Duplicates the check inside ``_apply_nan_policy`` intentionally:
    # a bad policy caught *here* fails before the expensive alignment
    # step, turning a wasted multi-field regrid into an instant error.
    _validate_choice(nan_policy, _VALID_NAN_POLICIES, "nan_policy")
    _validate_code_units_compatible(a, b, units)
    # Resolve against the *original* datasets: alignment rebuilds them
    # without custom aliases (``transform_to`` drops them, ``regrid`` keeps
    # only geometry defaults), so a post-alignment lookup would lose anything
    # from ``from_arrays(aliases=...)``.  Canonical names survive both.
    canonical = _resolve_common_field(a, b, field)
    a, b = _align_frames(a, b, frame=frame)
    _warn_if_coarse_mismatch(a.grid, b.grid)
    # Regrid only the requested field — 50× cheaper than full-dataset
    # alignment on a multi-moment PIC dump.
    a_aligned, b_aligned = align_grids(a, b, fields=[canonical], method=method)
    va = _extract_values(a_aligned, canonical, units)
    vb = _extract_values(b_aligned, canonical, units)
    if metric == "l2":
        return float(l2_relative_error(va, vb, nan_policy=nan_policy))
    return float(linf_error(va, vb, nan_policy=nan_policy))

field_comparison_report(a, b, *, fields=None, units='si', method='linear', nan_policy='omit', frame=None)

Compute L2 and L∞ errors for every common field, plus grid context.

Aligns the datasets once, then loops over the requested (or shared) canonical field names, computing both norms per field. The grid context captures domain extent and resolution ratio so a reader can judge whether the comparison is physically meaningful — e.g. a 100× spacing mismatch between a kinetic-scale PIC run and an MHD-scale run is numerically computable but suspect.

Parameters:

Name	Type	Description	Default
a	FieldDataset	Datasets to compare.	required
b	FieldDataset	Datasets to compare.	required
fields	Iterable[str] | None	Field names to report on. None reports on the intersection of canonical field names present in both datasets. Explicit names may be aliases; they resolve through both datasets.	None
units	('si', 'code')	Unit convention; see compare_fields.	"si"
method	str	Interpolation method passed through to align_grids. Default "linear".	'linear'
nan_policy	('omit', 'propagate', 'raise')	Forwarded to the pure diagnostics; see compare_fields.	"omit"
frame	str | None	Reference frame to compare in; see compare_fields.	None

Returns:

Type	Description
dict	Report with keys "fields" — {canonical_name: {"l2": float, "linf": float}} "grid" — common-grid and per-axis resolution-ratio info "units" — echo of the units argument

Raises:

Type	Description
ValueError	If there are no common fields, or on invalid units.
KeyError	If an explicit fields entry is missing in either dataset.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4,), spacing=(1.0,), origin=(0.0,))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
...     grid, Normalization.identity(),
... )
>>> report = field_comparison_report(ds, ds)
>>> report["fields"]["B_1"]["l2"]
0.0
>>> report["units"]
'si'

Source code in src/pypic/comparison.py

def field_comparison_report(
    a: FieldDataset,
    b: FieldDataset,
    *,
    fields: Iterable[str] | None = None,
    units: str = "si",
    method: str = "linear",
    nan_policy: NanPolicy = "omit",
    frame: str | None = None,
) -> dict[str, Any]:
    r"""Compute L2 and L∞ errors for every common field, plus grid context.

    Aligns the datasets once, then loops over the requested (or shared)
    canonical field names, computing both norms per field. The grid
    context captures domain extent and resolution ratio so a reader can
    judge whether the comparison is physically meaningful — e.g. a
    100× spacing mismatch between a kinetic-scale PIC run and an
    MHD-scale run is numerically computable but suspect.

    Parameters
    ----------
    a, b : FieldDataset
        Datasets to compare.
    fields : Iterable[str] | None
        Field names to report on. ``None`` reports on the intersection
        of canonical field names present in both datasets. Explicit names
        may be aliases; they resolve through both datasets.
    units : {"si", "code"}
        Unit convention; see `compare_fields`.
    method : str
        Interpolation method passed through to
        [`align_grids`][pypic.regridding.align_grids]. Default ``"linear"``.
    nan_policy : {"omit", "propagate", "raise"}
        Forwarded to the pure diagnostics; see `compare_fields`.
    frame : str | None
        Reference frame to compare in; see `compare_fields`.

    Returns
    -------
    dict
        Report with keys

        - ``"fields"`` — ``{canonical_name: {"l2": float, "linf": float}}``
        - ``"grid"`` — common-grid and per-axis resolution-ratio info
        - ``"units"`` — echo of the *units* argument

    Raises
    ------
    ValueError
        If there are no common fields, or on invalid *units*.
    KeyError
        If an explicit *fields* entry is missing in either dataset.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4,), spacing=(1.0,), origin=(0.0,))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
    ...     grid, Normalization.identity(),
    ... )
    >>> report = field_comparison_report(ds, ds)
    >>> report["fields"]["B_1"]["l2"]
    0.0
    >>> report["units"]
    'si'
    """
    _validate_choice(units, _ALLOWED_UNITS, "units")
    _validate_choice(nan_policy, _VALID_NAN_POLICIES, "nan_policy")
    _validate_code_units_compatible(a, b, units)
    # Resolve names against the *originals* so custom aliases from
    # ``from_arrays(aliases=...)`` survive — both transform_to and
    # regrid drop user aliases, so resolving post-alignment would lose
    # them. Bad names also raise *before* the expensive alignment step.
    names = _resolve_field_list(a, b, fields)
    a, b = _align_frames(a, b, frame=frame)
    _warn_if_coarse_mismatch(a.grid, b.grid)
    a_aligned, b_aligned = align_grids(a, b, fields=names, method=method)

    per_field: dict[str, dict[str, float]] = {}
    for name in names:
        va = _extract_values(a_aligned, name, units)
        vb = _extract_values(b_aligned, name, units)
        per_field[name] = {
            "l2": float(l2_relative_error(va, vb, nan_policy=nan_policy)),
            "linf": float(linf_error(va, vb, nan_policy=nan_policy)),
        }

    grid_context: dict[str, Any] = {
        "common_dimensions": a_aligned.grid.dimensions,
        "common_spacing": a_aligned.grid.spacing,
        "common_origin": a_aligned.grid.origin,
        "resolution_ratio": _resolution_ratio(a.grid, b.grid),
        "source_a_dimensions": a.grid.dimensions,
        "source_b_dimensions": b.grid.dimensions,
    }
    return {"fields": per_field, "grid": grid_context, "units": units}

field_difference_dataset(a, b, *, fields=None, units='si', method='linear', frame=None)

Build a FieldDataset of pointwise differences on the common grid.

Each requested field is computed as a[name] - b[name] after the two datasets are aligned. The returned dataset inherits a's species, physics, and frame metadata, so the result plugs directly into plot_field_slice. For the three-panel A | B | diff layout, run align_grids yourself and pass the pair to plot_comparison — that path does not need this helper.

Unlike regrid, which preserves a's original metadata dict verbatim, this function replaces .metadata with a fresh {"comparison": {"source_frames": ..., "units": ...}} record — the diff is a new artifact, not a regrid of a, and any per-step provenance on the sources would be misleading if copied.

NaN cells in either input pass through the difference array unchanged (NaN minus anything = NaN). There is no nan_policy parameter because field_difference itself is pure subtraction; use compare_fields with nan_policy=... or l2_relative_error directly if you need masked reductions.

When units="si", the stored arrays carry SI values but the dataset's normalization is set to Normalization.identity so that FieldDataset.in_si returns the same values instead of re-applying the SI factor. The actual unit choice is recorded in metadata["comparison"]["units"] for provenance.

Parameters:

Name	Type	Description	Default
a	FieldDataset	Datasets to subtract. Grids may differ.	required
b	FieldDataset	Datasets to subtract. Grids may differ.	required
fields	Iterable[str] | None	Field names to include. None uses the full intersection of canonical names.	None
units	('si', 'code')	Unit convention; see compare_fields.	"si"
method	str	Interpolation method passed through to align_grids. Default "linear".	'linear'
frame	str | None	Reference frame for the result; see compare_fields. The returned dataset's frame attribute reflects this choice (a.frame when None, otherwise the requested frame).	None

Returns:

Type	Description
FieldDataset	Dataset on common_grid(a.grid, b.grid) containing a - b for each selected field.

Raises:

Type	Description
ValueError	If there are no common fields, or on invalid units.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4,), spacing=(1.0,), origin=(0.0,))
>>> a = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
...     grid, Normalization.identity(),
... )
>>> b = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 1.5, 2.5, 4.0])},
...     grid, Normalization.identity(),
... )
>>> diff = field_difference_dataset(a, b)
>>> diff["B_1"]
array([0. , 0.5, 0.5, 0. ])

Source code in src/pypic/comparison.py

def field_difference_dataset(
    a: FieldDataset,
    b: FieldDataset,
    *,
    fields: Iterable[str] | None = None,
    units: str = "si",
    method: str = "linear",
    frame: str | None = None,
) -> FieldDataset:
    r"""Build a FieldDataset of pointwise differences on the common grid.

    Each requested field is computed as ``a[name] - b[name]`` after the
    two datasets are aligned. The returned dataset inherits *a*'s
    species, physics, and frame metadata, so the result plugs directly
    into [`plot_field_slice`][pypic.plotting.plot_field_slice]. For the
    three-panel A | B | diff layout, run
    [`align_grids`][pypic.regridding.align_grids] yourself and pass the pair to
    [`plot_comparison`][pypic.plotting.plot_comparison] — that path does not
    need this helper.

    Unlike [`regrid`][pypic.regridding.regrid], which preserves *a*'s original
    metadata dict verbatim, this function **replaces** ``.metadata``
    with a fresh ``{"comparison": {"source_frames": ..., "units": ...}}``
    record — the diff is a new artifact, not a regrid of *a*, and any
    per-step provenance on the sources would be misleading if copied.

    NaN cells in either input pass through the difference array
    unchanged (NaN minus anything = NaN). There is no ``nan_policy``
    parameter because `field_difference` itself is pure
    subtraction; use `compare_fields` with ``nan_policy=...`` or
    [`l2_relative_error`][pypic.diagnostics.l2_relative_error] directly if you need
    masked reductions.

    When ``units="si"``, the stored arrays carry SI values but the
    dataset's normalization is set to `Normalization.identity`
    so that `FieldDataset.in_si` returns the same values instead
    of re-applying the SI factor. The actual unit choice is recorded
    in ``metadata["comparison"]["units"]`` for provenance.

    Parameters
    ----------
    a, b : FieldDataset
        Datasets to subtract. Grids may differ.
    fields : Iterable[str] | None
        Field names to include. ``None`` uses the full intersection of
        canonical names.
    units : {"si", "code"}
        Unit convention; see `compare_fields`.
    method : str
        Interpolation method passed through to
        [`align_grids`][pypic.regridding.align_grids]. Default ``"linear"``.
    frame : str | None
        Reference frame for the result; see `compare_fields`. The
        returned dataset's ``frame`` attribute reflects this choice
        (``a.frame`` when ``None``, otherwise the requested frame).

    Returns
    -------
    FieldDataset
        Dataset on ``common_grid(a.grid, b.grid)`` containing
        ``a - b`` for each selected field.

    Raises
    ------
    ValueError
        If there are no common fields, or on invalid *units*.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4,), spacing=(1.0,), origin=(0.0,))
    >>> a = FieldDataset.from_arrays(
    ...     {"B_1": np.array([1.0, 2.0, 3.0, 4.0])},
    ...     grid, Normalization.identity(),
    ... )
    >>> b = FieldDataset.from_arrays(
    ...     {"B_1": np.array([1.0, 1.5, 2.5, 4.0])},
    ...     grid, Normalization.identity(),
    ... )
    >>> diff = field_difference_dataset(a, b)
    >>> diff["B_1"]
    array([0. , 0.5, 0.5, 0. ])
    """
    _validate_choice(units, _ALLOWED_UNITS, "units")
    _validate_code_units_compatible(a, b, units)
    # Capture original frames *before* _align_frames for the provenance
    # record below; the metadata should reflect what the user passed in,
    # not the post-transform frame on B.
    source_frames = (a.frame, b.frame)
    # Resolve against originals — see field_comparison_report for the
    # rationale (custom aliases survive transform_to/regrid only when
    # resolved up front; bad names raise pre-alignment).
    names = _resolve_field_list(a, b, fields)
    a, b = _align_frames(a, b, frame=frame)
    _warn_if_coarse_mismatch(a.grid, b.grid)
    a_aligned, b_aligned = align_grids(a, b, fields=names, method=method)

    diff_fields: dict[str, FloatArray] = {}
    for name in names:
        va = _extract_values(a_aligned, name, units)
        vb = _extract_values(b_aligned, name, units)
        diff_fields[name] = field_difference(va, vb)

    # With units="si" the stored arrays are already SI, so the result needs
    # an identity normalization — otherwise ``in_si()`` would re-apply the
    # factor and silently double-convert.
    result_norm = Normalization.identity() if units == "si" else a_aligned.normalization

    return FieldDataset.from_arrays(
        diff_fields,
        a_aligned.grid,
        result_norm,
        species=list(a_aligned.species),
        physics=a_aligned.physics,
        frame=a_aligned.frame,
        transforms=dict(a_aligned.transforms),
        metadata={
            "comparison": {
                "source_frames": source_frames,
                "units": units,
            }
        },
        strict_fields=False,
    )