Add Abscissa class (#220)

* Adds the "Abscissa" class * Implements new SasData structure and fixes tests * Implements abscissa object * Refactors "deduce_qz" routine * Reformatting files * Adds test cases for SasMeasurement class --------- Co-authored-by: lucas-wilkins <lucas.wilkins@stfc.ac.uk>
lucas-wilkins · Mar 10, 2025 · Mar 10, 2025 · Mar 10, 2025 · Mar 10, 2025 · Mar 10, 2025
commit dc510234638dd2b2fb8fe62178b7d1397fcc9338
diff --git a/conftest.py b/conftest.py
@@ -1,11 +1,33 @@
 import pytest
 
+from sasdata.metadata import Instrument, Metadata, Source
+from sasdata.quantities.quantity import Quantity
+from sasdata.quantities.units import angstroms
+
 
 def pytest_addoption(parser):
-    parser.addoption(
-        "--show_plots", action="store_true", default=False, help="Display diagnostic plots during tests"
-    )
+    parser.addoption("--show_plots", action="store_true", default=False, help="Display diagnostic plots during tests")
+
 
 @pytest.fixture
 def show_plots(request):
     return request.config.getoption("--show_plots")
+
+
+@pytest.fixture
+def basic_metadata():
+    """A minimal metadata object."""
+    wavelength = Quantity(1.0, angstroms)
+    source = Source(
+        radiation=None,
+        beam_shape=None,
+        beam_size=None,
+        wavelength=wavelength,
+        wavelength_max=None,
+        wavelength_min=None,
+        wavelength_spread=None,
+    )
+    instrument = Instrument(collimations=[], source=source, detector=[])
+    metadata = Metadata(title=None, run=[], definition=None, process=[], sample=None, instrument=instrument, raw=None)
+
+    return metadata
diff --git a/sasdata/abscissa.py b/sasdata/abscissa.py
@@ -0,0 +1,121 @@
+from abc import ABC, abstractmethod
+
+import numpy as np
+from numpy._typing import ArrayLike
+
+from sasdata.exceptions import InterpretationError
+from sasdata.quantities.quantity import Quantity
+from sasdata.util import is_increasing
+
+
+class Abscissa(ABC):
+    def __init__(self, axes: list[Quantity]):
+        self._axes = axes
+        self._dimensionality = len(axes)
+
+    @property
+    def dimensionality(self) -> int:
+        """Dimensionality of this data"""
+        return self._dimensionality
+
+    @property
+    @abstractmethod
+    def is_grid(self) -> bool:
+        """Are these coordinates using a grid representation
+        (as opposed to a general list representation)
+
+        is_grid = True: implies that the corresponding ordinate is n-dimensional tensor
+        is_grid = False: implies that the corresponding ordinate is a 1D list
+
+        If the data is one dimensional, is_grid=True
+        """
+
+    @property
+    def axes(self) -> list[Quantity]:
+        """Axes of the data:
+
+        If it's an (n1-by-n2-by-n3...) grid (is_grid=True): give the values for each axis, returning a list like
+          [Quantity(length n1), Quantity(length n2), Quantity(length n3) ... ]
+
+        If it is not grid data (is_grid=False), but n points on a general mesh, give one array for each dimension
+          [Quantity(length n), Quantity(length n), Quantity(length n) ... ]
+        """
+        return self._axes
+
+    @staticmethod
+    def _determine_error_message(axis_arrays: list[np.ndarray], ordinate_shape: tuple):
+        """Error message for the `.determine` function"""
+
+        shape_string = ", ".join([str(axis.shape) for axis in axis_arrays])
+
+        return f"Cannot interpret array shapes axis: [{shape_string}], ordinate: {ordinate_shape}"
+
+    @staticmethod
+    def determine(axis_data: list[Quantity[ArrayLike]], ordinate_data: Quantity[ArrayLike]) -> "Abscissa":
+        """Get an Abscissa object that fits the combination of axes and data"""
+
+        # Different posibilites:
+        #   1: axes_data[i].shape == axes_data[j].shape == ordinate_data.shape
+        #    1a: axis_data[i] is 1D =>
+        #      1a-i:  len(axes_data) == 1 => Grid type or Scatter type depending on sortedness
+        #      1a-ii: len(axes_data) != 1 => Scatter type
+        #    1b: axis_data[i] dimensionality matches len(axis_data) => Meshgrid type
+        #    1c: other => Error
+        #   2: (len(axes_data[0]), len(axes_data[1])... ) == ordinate_data.shape => Grid type
+        #   3: None of the above => Error
+
+        ordinate_shape = np.array(ordinate_data.value).shape
+        axis_arrays = [np.array(axis.value) for axis in axis_data]
+
+        # 1:
+        if all([axis.shape == ordinate_shape for axis in axis_arrays]):
+            # 1a:
+            if all([len(axis.shape) == 1 for axis in axis_arrays]):
+                # 1a-i:
+                if len(axis_arrays) == 1:
+                    # Is it sorted
+                    if is_increasing(axis_arrays[0]):
+                        return GridAbscissa(axis_data)
+                    else:
+                        return ScatterAbscissa(axis_data)
+                # 1a-ii
+                else:
+                    return ScatterAbscissa(axis_data)
+            # 1b
+            elif all([len(axis.shape) == len(axis_arrays) for axis in axis_arrays]):
+                return MeshgridAbscissa(axis_data)
+
+            else:
+                raise InterpretationError(Abscissa._determine_error_message(axis_arrays, ordinate_shape))
+
+        elif (
+            all([len(axis.shape) == 1 for axis in axis_arrays])
+            and tuple([axis.shape[0] for axis in axis_arrays]) == ordinate_shape
+        ):
+            # Require that they are sorted
+            if all([is_increasing(axis) for axis in axis_arrays]):
+                return GridAbscissa(axis_data)
+
+            else:
+                raise InterpretationError("Grid axes are not sorted")
+
+        else:
+            raise InterpretationError(Abscissa._determine_error_message(axis_arrays, ordinate_shape))
+
+
+class GridAbscissa(Abscissa):
+    @property
+    def is_grid(self):
+        return True
+
+
+class MeshgridAbscissa(Abscissa):
+    @property
+    def is_grid(self):
+        return True
+
+
+class ScatterAbscissa(Abscissa):
+    @property
+    def is_grid(self):
+        return False
diff --git a/sasdata/data.py b/sasdata/data.py
@@ -7,12 +7,66 @@
 from h5py._hl.group import Group as HDF5Group
 
 from sasdata import dataset_types
+from sasdata.abscissa import Abscissa
 from sasdata.dataset_types import DatasetType
-from sasdata.metadata import Metadata, MetadataEncoder
+from sasdata.metadata import DerivedMetadata, Metadata, MetadataEncoder
 from sasdata.quantities.quantity import Quantity
 
 
 class SasData:
+    """General object containing data in the SasView ecosystem"""
+
+    def __init__(
+        self,
+        name: str,
+        ordinate: Quantity,
+        abscissae: Abscissa,
+        mask: Quantity,
+        dependents: list["SasData"],
+        metadata: Metadata,
+    ):
+        self.name = name
+        self._ordinate = ordinate
+        self._abscissae = abscissae
+        self._mask = mask
+        self.dependents = dependents
+        self.metadata = metadata
+
+    @property
+    def ordinate(self) -> Quantity:
+        return self._ordinate
+
+    @property
+    def abscissae(self) -> Abscissa:
+        return self._abscissae
+
+    @property
+    def mask(self) -> Quantity:
+        return self._mask
+
+    def scatter_data(self):
+        """Return data in the coordinate/value form [(x1, x2, x3, y)...]"""
+
+
+class SasDerivedMeasurement(SasData):
+    """General object sas measurement that has not come directly from a file,
+    for example, the difference between two datasets"""
+
+    def __init__(
+        self,
+        name: str,
+        ordinate: Quantity,
+        abscissae: Abscissa,
+        mask: Quantity,
+        dependents: list["SasData"],
+        metadata: DerivedMetadata,
+    ):
+        super().__init__(
+            name=name, ordinate=ordinate, abscissae=abscissae, mask=mask, dependents=dependents, metadata=metadata
+        )
+
+
+class SasMeasurement(SasData):
     def __init__(
         self,
         name: str,
@@ -34,52 +88,37 @@ def __init__(
         self.dataset_type: DatasetType = dataset_type
 
         # Components that need to be organised after creation
-        self.mask = None  # TODO: fill out
+        self._mask = None  # TODO: fill out
         self.model_requirements = None  # TODO: fill out
 
     # TODO: Handle the other data types.
     @property
     def ordinate(self) -> Quantity:
         match self.dataset_type:
-            case (dataset_types.one_dim |
-                  dataset_types.two_dim |
-                  dataset_types.angle_dim):
+            case dataset_types.one_dim | dataset_types.two_dim | dataset_types.three_dim | dataset_types.angle_dim:
                 return self._data_contents["I"]
             case dataset_types.sesans:
                 return self._data_contents["Depolarisation"]
             case _:
                 return None
 
     @property
-    def abscissae(self) -> Quantity:
+    def abscissae(self) -> Abscissa:
         match self.dataset_type:
             case dataset_types.one_dim:
-                return self._data_contents["Q"]
+                return Abscissa.determine([self._data_contents["Q"]], self.ordinate)
             case dataset_types.two_dim:
-                # Type hinting is a bit lacking. Assume each part of the zip is a scalar value.
-                data_contents = np.array(
-                    list(
-                        zip(
-                            self._data_contents["Qx"].value,
-                            self._data_contents["Qy"].value,
-                        )
-                    )
-                )
-                # Use this value to extract units etc. Assume they will be the same for Qy.
-                reference_data_content = self._data_contents["Qx"]
-                # TODO: If this is a derived quantity then we are going to lose that
-                # information.
-                #
-                # TODO: Won't work when there's errors involved. On reflection, we
-                # probably want to avoid creating a new Quantity but at the moment I
-                # can't see a way around it.
-                return Quantity(data_contents, reference_data_content.units, name=self._data_contents["Qx"].name, id_header=self._data_contents["Qx"]._id_header)
+                return Abscissa.determine([self._data_contents["Qx"], self._data_contents["Qy"]], self.ordinate)
             case dataset_types.angle_dim:
-                return self._data_contents["Phi"]
+                return Abscissa.determine([self._data_contents["Phi"]], self.ordinate)
+            case dataset_types.three_dim:
+                return Abscissa.determine(
+                    [self._data_contents["Qx"], self._data_contents["Qy"], self._data_contents["Qz"]], self.ordinate
+                )
             case dataset_types.sesans:
-                return self._data_contents["SpinEchoLength"]
+                return Abscissa.determine([self._data_contents["SpinEchoLength"]], self.ordinate)
             case _:
-                None
+                return None
 
     def __getitem__(self, item: str):
         return self._data_contents[item]
@@ -98,7 +137,7 @@ def summary(self, indent="  "):
 
     @staticmethod
     def from_json(obj):
-        return SasData(
+        return SasMeasurement(
             name=obj["name"],
             dataset_type=DatasetType(
                 name=obj["type"]["name"],
@@ -121,12 +160,11 @@ def _save_h5(self, sasentry: HDF5Group):
         for idx, (key, sasdata) in enumerate(self._data_contents.items()):
             sasdata.as_h5(group, key)
 
-
     @staticmethod
     def save_h5(data: dict[str, typing.Self], path: str | typing.BinaryIO):
         with h5py.File(path, "w") as f:
             for idx, (key, data) in enumerate(data.items()):
-                sasentry = f.create_group(f"sasentry{idx+1:02d}")
+                sasentry = f.create_group(f"sasentry{idx + 1:02d}")
                 if not key.startswith("sasentry"):
                     sasentry.attrs["sasview_key"] = key
                 data._save_h5(sasentry)
@@ -140,7 +178,7 @@ def deserialise(data: str) -> "SasData":
     def deserialise_json(json_data: dict) -> "SasData":
         name = json_data["name"]
         data_contents = {}
-        dataset_type = json_data["dataset_type"] # TODO: update when DatasetType is more finalized
+        dataset_type = json_data["dataset_type"]  # TODO: update when DatasetType is more finalized
         metadata = json_data["metadata"].deserialise_json()
         for quantity in json_data["data_contents"]:
             data_contents[quantity["label"]] = Quantity.deserialise_json(quantity)
@@ -159,11 +197,11 @@ def _serialise_json(self) -> dict[str, Any]:
         return {
             "name": self.name,
             "data_contents": data,
-            "dataset_type": None, # TODO: update when DatasetType is more finalized
+            "dataset_type": None,  # TODO: update when DatasetType is more finalized
             "verbose": self._verbose,
             "metadata": self.metadata.serialise_json(),
             "mask": {},
-            "model_requirements": {}
+            "model_requirements": {},
         }
 
 
@@ -210,16 +248,16 @@ def sasdata_reader2D_converter(data2d: SasData | None = None) -> SasData:
     qx_data = new_x.flatten()
     qy_data = new_y.flatten()
     err_data = np.sqrt(data2d._data_contents["I"].variance.value)
-    if not data2d._data_contents["I"].has_variance or np.any(err_data <= 0):
+    if not data2d._data_contents["I"].has_error or np.any(err_data <= 0):
         new_err_data = np.sqrt(np.abs(new_data))
     else:
         new_err_data = err_data.flatten()
     mask = np.ones(len(new_data), dtype=bool)
 
     data2d._data_contents["I"].value = new_data
-    data2d._data_contents["I"].variance.value = new_err_data ** 2
+    data2d._data_contents["I"].variance.value = new_err_data**2
     data2d._data_contents["Qx"].value = qx_data
     data2d._data_contents["Qy"].value = qy_data
-    data2d.mask = mask
+    data2d._mask = mask
 
     return data2d