Interpolation and Rebinning

requirements.txt

@@ -4,10 +4,14 @@ lxml
 
 # Calculation
 numpy==1.*
+scipy
 
 # Unit testing
 pytest
 unittest-xml-reporting
 
 # Documentation
 sphinx
+
+# Other stuff
+matplotlib

sasdata/data.py

@@ -2,29 +2,32 @@
 from typing import TypeVar, Any, Self
 from dataclasses import dataclass
 
+import numpy as np
+
 from quantities.quantity import NamedQuantity
 from sasdata.metadata import Metadata
 from sasdata.quantities.accessors import AccessorTarget
 from sasdata.data_backing import Group, key_tree
 
 
 class SasData:
-    def __init__(self, name: str, data_contents: list[NamedQuantity], raw_metadata: Group, verbose: bool=False):
+    def __init__(self, name: str,
+                 data_contents: list[NamedQuantity],
+                 raw_metadata: Group,
+                 verbose: bool=False):
+
         self.name = name
         self._data_contents = data_contents
         self._raw_metadata = raw_metadata
         self._verbose = verbose
 
         self.metadata = Metadata(AccessorTarget(raw_metadata, verbose=verbose))
 
-        # TO IMPLEMENT
-
-        # abscissae: list[NamedQuantity[np.ndarray]]
-        # ordinate: NamedQuantity[np.ndarray]
-        # other: list[NamedQuantity[np.ndarray]]
-        #
-        # metadata: Metadata
-        # model_requirements: ModellingRequirements
+        # Components that need to be organised after creation
+        self.ordinate: NamedQuantity[np.ndarray] = None # TODO: fill out
+        self.abscissae: list[NamedQuantity[np.ndarray]] = None # TODO: fill out
+        self.mask = None # TODO: fill out
+        self.model_requirements = None # TODO: fill out
 
     def summary(self, indent = "  ", include_raw=False):
         s = f"{self.name}\n"

sasdata/manual_tests/interpolation.py

@@ -0,0 +1,44 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from sasdata.quantities.quantity import NamedQuantity
+from sasdata.quantities.plotting import quantity_plot
+from sasdata.quantities import units
+
+from sasdata.transforms.rebinning import calculate_interpolation_matrix_1d
+from sasdata.transforms.rebinning import InterpolationOptions
+
+def linear_interpolation_check():
+
+    for from_bins in [(-10, 10, 10),
+                      (-10, 10, 1000),
+                      (-15, 5, 10),
+                      (15,5, 10)]:
+        for to_bins in [
+            (-15, 0, 10),
+            (-15, 15, 10),
+            (0, 20, 100)]:
+
+            plt.figure()
+
+            x = NamedQuantity("x", np.linspace(*from_bins), units=units.meters)
+            y = x**2
+
+            quantity_plot(x, y)
+
+            new_x = NamedQuantity("x_new", np.linspace(*to_bins), units=units.meters)
+
+            rebin_mat = calculate_interpolation_matrix_1d(x, new_x, order=InterpolationOptions.LINEAR)
+
+            new_y = y @ rebin_mat
+
+            quantity_plot(new_x, new_y)
+
+            print(new_y.history.summary())
+
+    plt.show()
+
+
+
+
+linear_interpolation_check()

sasdata/model_requirements.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from sasdata.metadata import Metadata
-from transforms.operation import Operation
+from sasdata.quantities.quantity import Operation
 
 
 @dataclass

sasdata/quantities/math_operations_test.py

@@ -0,0 +1,152 @@
+""" Tests for math operations """
+
+import pytest
+
+import numpy as np
+from sasdata.quantities.quantity import NamedQuantity, tensordot, transpose
+from sasdata.quantities import units
+
+order_list = [
+    [0, 1, 2, 3],
+    [0, 2, 1],
+    [1, 0],
+    [0, 1],
+    [2, 0, 1],
+    [3, 1, 2, 0]
+]
+
+@pytest.mark.parametrize("order", order_list)
+def test_transpose_raw(order: list[int]):
+    """ Check that the transpose operation changes the order of indices correctly - uses sizes as way of tracking"""
+
+    input_shape = tuple([i+1 for i in range(len(order))])
+    expected_shape = tuple([i+1 for i in order])
+
+    input_mat = np.zeros(input_shape)
+
+    measured_mat = transpose(input_mat, axes=tuple(order))
+
+    assert measured_mat.shape == expected_shape
+
+
+@pytest.mark.parametrize("order", order_list)
+def test_transpose_raw(order: list[int]):
+    """ Check that the transpose operation changes the order of indices correctly - uses sizes as way of tracking"""
+    input_shape = tuple([i + 1 for i in range(len(order))])
+    expected_shape = tuple([i + 1 for i in order])
+
+    input_mat = NamedQuantity("testmat", np.zeros(input_shape), units=units.none)
+
+    measured_mat = transpose(input_mat, axes=tuple(order))
+
+    assert measured_mat.value.shape == expected_shape
+
+
+rng_seed = 1979
+tensor_product_with_identity_sizes = (4,6,5)
+
+@pytest.mark.parametrize("index, size", [tup for tup in enumerate(tensor_product_with_identity_sizes)])
+def test_tensor_product_with_identity_quantities(index, size):
+    """ Check the correctness of the tensor product by multiplying by the identity (quantity, quantity)"""
+    np.random.seed(rng_seed)
+
+    x = NamedQuantity("x", np.random.rand(*tensor_product_with_identity_sizes), units=units.meters)
+    y = NamedQuantity("y", np.eye(size), units.seconds)
+
+    z = tensordot(x, y, index, 0)
+
+    # Check units
+    assert z.units == units.meters * units.seconds
+
+    # Expected sizes - last index gets moved to end
+    output_order = [i for i in (0, 1, 2) if i != index] + [index]
+    output_sizes = [tensor_product_with_identity_sizes[i] for i in output_order]
+
+    assert z.value.shape == tuple(output_sizes)
+
+    # Restore original order and check
+    reverse_order = [-1, -1, -1]
+    for to_index, from_index in enumerate(output_order):
+        reverse_order[from_index] = to_index
+
+    z_reordered = transpose(z, axes = tuple(reverse_order))
+
+    assert z_reordered.value.shape == tensor_product_with_identity_sizes
+
+    # Check values
+
+    mat_in = x.in_si()
+    mat_out = transpose(z, axes=tuple(reverse_order)).in_si()
+
+    assert np.all(np.abs(mat_in - mat_out) < 1e-10)
+
+
+@pytest.mark.parametrize("index, size", [tup for tup in enumerate(tensor_product_with_identity_sizes)])
+def test_tensor_product_with_identity_quantity_matrix(index, size):
+    """ Check the correctness of the tensor product by multiplying by the identity (quantity, matrix)"""
+    np.random.seed(rng_seed)
+
+    x = NamedQuantity("x", np.random.rand(*tensor_product_with_identity_sizes), units.meters)
+    y = np.eye(size)
+
+    z = tensordot(x, y, index, 0)
+
+    assert z.units == units.meters
+
+    # Expected sizes - last index gets moved to end
+    output_order = [i for i in (0, 1, 2) if i != index] + [index]
+    output_sizes = [tensor_product_with_identity_sizes[i] for i in output_order]
+
+    assert z.value.shape == tuple(output_sizes)
+
+    # Restore original order and check
+    reverse_order = [-1, -1, -1]
+    for to_index, from_index in enumerate(output_order):
+        reverse_order[from_index] = to_index
+
+    z_reordered = transpose(z, axes = tuple(reverse_order))
+
+    assert z_reordered.value.shape == tensor_product_with_identity_sizes
+
+    # Check values
+
+    mat_in = x.in_si()
+    mat_out = transpose(z, axes=tuple(reverse_order)).in_si()
+
+    assert np.all(np.abs(mat_in - mat_out) < 1e-10)
+
+
+@pytest.mark.parametrize("index, size", [tup for tup in enumerate(tensor_product_with_identity_sizes)])
+def test_tensor_product_with_identity_matrix_quantity(index, size):
+    """ Check the correctness of the tensor product by multiplying by the identity (matrix, quantity)"""
+    np.random.seed(rng_seed)
+
+    x = np.random.rand(*tensor_product_with_identity_sizes)
+    y = NamedQuantity("y", np.eye(size), units.seconds)
+
+    z = tensordot(x, y, index, 0)
+
+    assert z.units == units.seconds
+
+
+    # Expected sizes - last index gets moved to end
+    output_order = [i for i in (0, 1, 2) if i != index] + [index]
+    output_sizes = [tensor_product_with_identity_sizes[i] for i in output_order]
+
+    assert z.value.shape == tuple(output_sizes)
+
+    # Restore original order and check
+    reverse_order = [-1, -1, -1]
+    for to_index, from_index in enumerate(output_order):
+        reverse_order[from_index] = to_index
+
+    z_reordered = transpose(z, axes = tuple(reverse_order))
+
+    assert z_reordered.value.shape == tensor_product_with_identity_sizes
+
+    # Check values
+
+    mat_in = x
+    mat_out = transpose(z, axes=tuple(reverse_order)).in_si()
+
+    assert np.all(np.abs(mat_in - mat_out) < 1e-10)

sasdata/quantities/numerical_encoding.py

@@ -0,0 +1,72 @@
+import numpy as np
+from scipy.sparse import coo_matrix, csr_matrix, csc_matrix, coo_array, csr_array, csc_array
+
+import base64
+import struct
+
+
+def numerical_encode(obj: int | float | np.ndarray | coo_matrix | coo_array | csr_matrix | csr_array | csc_matrix | csc_array):
+
+    if isinstance(obj, int):
+        return {"type": "int",
+                "value": obj}
+
+    elif isinstance(obj, float):
+        return {"type": "float",
+                "value": base64.b64encode(bytearray(struct.pack('d', obj)))}
+
+    elif isinstance(obj, np.ndarray):
+        return {
+            "type": "numpy",
+            "value": base64.b64encode(obj.tobytes()),
+            "dtype": obj.dtype.str,
+            "shape": list(obj.shape)
+        }
+
+    elif isinstance(obj, (coo_matrix, coo_array, csr_matrix, csr_array, csc_matrix, csc_array)):
+
+        output = {
+            "type": obj.__class__.__name__, # not robust to name changes, but more concise
+            "dtype": obj.dtype.str,
+            "shape": list(obj.shape)
+        }
+
+        if isinstance(obj, (coo_array, coo_matrix)):
+
+            output["data"] = numerical_encode(obj.data)
+            output["coords"] = [numerical_encode(coord) for coord in obj.coords]
+
+
+        elif isinstance(obj, (csr_array, csr_matrix)):
+            pass
+
+
+        elif isinstance(obj, (csc_array, csc_matrix)):
+
+            pass
+
+
+        return output
+
+    else:
+        raise TypeError(f"Cannot serialise object of type: {type(obj)}")
+
+def numerical_decode(data: dict[str, str | int | list[int]]) -> int | float | np.ndarray | coo_matrix | coo_array | csr_matrix | csr_array | csc_matrix | csc_array:
+    obj_type = data["type"]
+
+    match obj_type:
+        case "int":
+            return int(data["value"])
+
+        case "float":
+            return struct.unpack('d', base64.b64decode(data["value"]))[0]
+
+        case "numpy":
+            value = base64.b64decode(data["value"])
+            dtype = np.dtype(data["dtype"])
+            shape = tuple(data["shape"])
+            return np.frombuffer(value, dtype=dtype).reshape(*shape)
+
+        case _:
+            raise ValueError(f"Cannot decode objects of type '{obj_type}'")
+