From 4ab3265a3806146c911ff5e06e5704bcf7dad777 Mon Sep 17 00:00:00 2001 From: PaulSharp <44529197+DrPaulSharp@users.noreply.github.com> Date: Mon, 16 Mar 2026 18:39:04 +0000 Subject: [PATCH 1/4] Add trace and determinant to quantity operations --- requirements.txt | 1 + sasdata/quantities/quantity.py | 99 ++++++++++++++- test/quantities/utest_math_operations.py | 151 +++++++++++------------ test/quantities/utest_operations.py | 21 +++- 4 files changed, 191 insertions(+), 81 deletions(-) diff --git a/requirements.txt b/requirements.txt index 64bfe29a4..623acaac1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,6 +5,7 @@ lxml # Calculation numpy scipy +sympy # Unit testing pytest diff --git a/sasdata/quantities/quantity.py b/sasdata/quantities/quantity.py index 08e0be6fe..8088d488a 100644 --- a/sasdata/quantities/quantity.py +++ b/sasdata/quantities/quantity.py @@ -5,6 +5,7 @@ import h5py import numpy as np +import sympy as sp from numpy._typing import ArrayLike from sasdata.quantities import units @@ -39,6 +40,30 @@ def transpose(a: Union["Quantity[ArrayLike]", ArrayLike], axes: tuple | None = N return np.transpose(a, axes=axes) +def trace(a: Union["Quantity[ArrayLike]", ArrayLike], offset: int = 0, axis1: int = 0, axis2: int = 1): + """Find the trace of an array or an array based quantity.""" + if isinstance(a, Quantity): + return DerivedQuantity( + value=np.trace(a.value, offset, axis1, axis2), + units=a.units, + history=QuantityHistory.apply_operation(Trace, a.history, offset=offset, axis1=axis1, axis2=axis2), + ) + else: + return np.trace(a, offset, axis1, axis2) + + +def determinant(a: Union["Quantity[ArrayLike]", ArrayLike]): + """Find the determinant of an array or an array based quantity.""" + if isinstance(a, Quantity): + return DerivedQuantity( + value=np.linalg.det(a.value), + units=a.units, # TODO: hmmm, need to think about this + history=QuantityHistory.apply_operation(Determinant, a.history), + ) + else: + return np.linalg.det(a) + + def dot(a: Union["Quantity[ArrayLike]", ArrayLike], b: Union["Quantity[ArrayLike]", ArrayLike]): """Dot product of two arrays or two array based quantities""" a_is_quantity = isinstance(a, Quantity) @@ -999,14 +1024,14 @@ def __init__(self, a: Operation, axes: tuple[int] | None = None): self.axes = axes def evaluate(self, variables: dict[int, T]) -> T: - return np.transpose(self.a.evaluate(variables)) + return np.transpose(self.a.evaluate(variables), self.axes) def _derivative(self, hash_value: int) -> Operation: - return Transpose(self.a.derivative(hash_value)) # TODO: Check! + return Transpose(self.a.derivative(hash_value), self.axes) # TODO: Check! def _clean(self): clean_a = self.a._clean() - return Transpose(clean_a) + return Transpose(clean_a, self.axes) def _serialise_parameters(self) -> dict[str, Any]: if self.axes is None: @@ -1044,6 +1069,72 @@ def __eq__(self, other): return False +class Trace(Operation): + """Trace operation - as per numpy""" + + serialisation_name = "trace" + + def __init__(self, a: Operation, offset: int = 0, axis1: int = 0, axis2: int = 1): + self.a = a + self.offset = offset + self.axis1 = axis1 + self.axis2 = axis2 + + def evaluate(self, variables: dict[int, T]) -> T: + return np.trace(self.a.evaluate(variables), self.offset, self.axis1, self.axis2) + + def _derivative(self, hash_value: int) -> Operation: + return Trace(self.a.derivative(hash_value), self.offset, self.axis1, self.axis2) + + def _clean(self): + clean_a = self.a._clean() + return Trace(clean_a, self.offset, self.axis1, self.axis2) + + def _serialise_parameters(self) -> dict[str, Any]: + return {"a": self.a._serialise_json(), "offset": self.offset, "axis1": self.axis1, "axis2": self.axis2} + + @staticmethod + def _deserialise(parameters: dict) -> "Operation": + return Trace( + a=Operation.deserialise_json(parameters["a"]), + offset=parameters["offset"], + axis1=parameters["axis1"], + axis2=parameters["axis2"], + ) + + def summary(self, indent_amount: int = 0, indent=" "): + return ( + f"{indent_amount * indent}Trace(\n" + + self.a.summary(indent_amount + 1, indent) + + "\n" + + f"{(indent_amount + 1) * indent}{self.offset}\n" + + f"{(indent_amount + 1) * indent}{self.axis1}\n" + + f"{(indent_amount + 1) * indent}{self.axis2}\n" + + f"{indent_amount * indent})" + ) + + def __eq__(self, other): + if isinstance(other, Trace): + return other.a == self.a + return False + + +class Determinant(UnaryOperation): + """Array Determinant - backed by numpy's det method""" + + serialisation_name = "determinant" + + def evaluate(self, variables: dict[int, T]) -> T: + return np.linalg.det(self.a.evaluate(variables)) + + def _derivative(self, hash_value: int) -> Operation: + return Trace(sp.adjugate(self.a) * self.a._derivative(hash_value)) + + def _clean(self): + clean_a = self.a._clean() + return Determinant(clean_a) # Do nothing for now - can consider identity matrix and common columns + + class Dot(BinaryOperation): """Dot product - backed by numpy's dot method""" @@ -1141,6 +1232,8 @@ def _deserialise(parameters: dict) -> "Operation": Pow, Log, Transpose, + Trace, + Determinant, Dot, MatMul, TensorDot, diff --git a/test/quantities/utest_math_operations.py b/test/quantities/utest_math_operations.py index 85e9ba1d6..3b38bad55 100644 --- a/test/quantities/utest_math_operations.py +++ b/test/quantities/utest_math_operations.py @@ -1,62 +1,83 @@ -""" Tests for math operations """ +"""Tests for math operations""" import numpy as np import pytest from sasdata.quantities import units -from sasdata.quantities.quantity import NamedQuantity, tensordot, transpose +from sasdata.quantities.quantity import NamedQuantity, tensordot, trace, transpose -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 +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_with_quantity(order: list[int]): - """ Check that the transpose operation changes the order of indices correctly - uses sizes as way of tracking""" +def test_transpose(order: list[int]): + """Check that the transpose operation changes the order of indices correctly for raw data and quantities - 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) + input_mat = np.zeros(input_shape) + input_quantity = NamedQuantity("testmat", np.zeros(input_shape), units=units.none) measured_mat = transpose(input_mat, axes=tuple(order)) + measured_quantity = transpose(input_quantity, axes=tuple(order)) - assert measured_mat.value.shape == expected_shape + assert measured_mat.shape == expected_shape + assert measured_quantity.value.shape == expected_shape + + +@pytest.mark.parametrize( + "matrix, offset, expected_trace", + [ + (np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), 0, 15), + (np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), 1, 8), + (np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), 2, 3), + (np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), -1, 12), + ], +) +def test_trace_offset(matrix, offset, expected_trace): + """Check that the trace operation correctly identifies the offset value for raw data and quantities.""" + assert (trace(matrix, offset=offset) == expected_trace).all() + assert (trace(NamedQuantity("testmat", matrix, units=units.none), offset).value == expected_trace).all() + + +@pytest.mark.parametrize( + "matrix, axis1, axis2, expected_trace", + [ + (np.array([[[1, 2], [3, 4]], [[1, 2], [3, 4]], [[1, 2], [3, 4]]]), 0, 1, np.array([4, 6])), + (np.array([[[1, 2], [3, 4]], [[1, 2], [3, 4]], [[1, 2], [3, 4]]]), 1, 2, np.array([5, 5, 5])), + (np.array([[[1, 2], [3, 4]], [[1, 2], [3, 4]], [[1, 2], [3, 4]]]), 0, 2, np.array([3, 7])), + ], +) +def test_trace_axes(matrix, axis1, axis2, expected_trace): + """Check that the trace operation correctly identifies the offset value for raw data and quantities.""" + assert (trace(matrix, axis1=axis1, axis2=axis2) == expected_trace).all() + assert ( + trace(NamedQuantity("testmat", matrix, units=units.none), axis1=axis1, axis2=axis2).value == expected_trace + ).all() rng_seed = 1979 -tensor_product_with_identity_sizes = (4,6,5) +tensor_product_with_identity_sizes = (4, 6, 5) + +@pytest.mark.parametrize( + "x, x_unit", + [ + (NamedQuantity("x", np.random.rand(*tensor_product_with_identity_sizes), units=units.meters), units.meters), + ((np.random.rand(*tensor_product_with_identity_sizes), units.none)), + ], +) @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)""" +def test_tensor_product_with_identity_quantities(x, x_unit, 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 + assert z.units == x_unit * units.seconds # Expected sizes - last index gets moved to end output_order = [i for i in (0, 1, 2) if i != index] + [index] @@ -69,29 +90,37 @@ def test_tensor_product_with_identity_quantities(index, size): for to_index, from_index in enumerate(output_order): reverse_order[from_index] = to_index - z_reordered = transpose(z, axes = tuple(reverse_order)) + 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() + try: + mat_in = x.in_si() + except AttributeError: + mat_in = x + mat_out = transpose(z, axes=tuple(reverse_order)).in_si() assert np.all(np.abs(mat_in - mat_out) < 1e-10) +@pytest.mark.parametrize( + "x, x_unit", + [ + (NamedQuantity("x", np.random.rand(*tensor_product_with_identity_sizes), units=units.meters), units.meters), + ], +) @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)""" +def test_tensor_product_with_identity_quantity_matrix(x, x_unit, 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 + assert z.units == x_unit # Expected sizes - last index gets moved to end output_order = [i for i in (0, 1, 2) if i != index] + [index] @@ -104,49 +133,17 @@ def test_tensor_product_with_identity_quantity_matrix(index, size): for to_index, from_index in enumerate(output_order): reverse_order[from_index] = to_index - z_reordered = transpose(z, axes = tuple(reverse_order)) + 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 + try: + mat_in = x.in_si() + except AttributeError: + mat_in = x - mat_in = x mat_out = transpose(z, axes=tuple(reverse_order)).in_si() assert np.all(np.abs(mat_in - mat_out) < 1e-10) diff --git a/test/quantities/utest_operations.py b/test/quantities/utest_operations.py index 5b64eab0e..586c1ba26 100644 --- a/test/quantities/utest_operations.py +++ b/test/quantities/utest_operations.py @@ -11,6 +11,7 @@ ArcTan, Constant, Cos, + Determinant, Div, Dot, Exp, @@ -26,6 +27,7 @@ Sin, Sub, Tan, + Trace, Transpose, Variable, ) @@ -46,7 +48,7 @@ ) -@pytest.fixture(params=[Inv, Exp, Ln, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose]) +@pytest.fixture(params=[Determinant, Inv, Exp, Ln, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose]) def unary_operation(request): return request.param(x) @@ -93,6 +95,14 @@ def test_log_pow_serialise_deserialise(log_pow_operation): assert serialised == reserialised +def test_trace_serialise_deserialise(): + serialised = Trace(x).serialise() + deserialised = Operation.deserialise(serialised) + reserialised = deserialised.serialise() + + assert serialised == reserialised + + @pytest.mark.parametrize( "op, summary", [(AdditiveIdentity, "0 [Add.Id.]"), (MultiplicativeIdentity, "1 [Mul.Id.]"), (Operation, "Operation(\n)")], @@ -118,6 +128,11 @@ def test_log_pow_summary(log_pow_operation): assert log_pow_operation.summary() == f"{log_pow_operation.__class__.__name__}(\n x\n 2\n)" +def test_trace_summary(): + op = Trace(x) + assert op.summary() == f"{op.__class__.__name__}(\n x\n 0\n 0\n 1\n)" + + @pytest.mark.parametrize("op, result", [(AdditiveIdentity, 0), (MultiplicativeIdentity, 1), (Operation, None)]) def test_evaluation(op, result): f = op() @@ -147,6 +162,10 @@ def test_evaluation(op, result): (ArcCos, -1.0, math.pi), (ArcTan, 0.0, 0.0), (ArcTan, -1.0, -0.25 * math.pi), + (Determinant, np.array([[1, 2], [3, 4]]), pytest.approx(-2.0)), + (Determinant, np.array([[1, 4, 1], [2, 4, 2], [3, 4, 3]]), pytest.approx(0.0)), + (Trace, np.array([[1, 2], [3, 4]]), pytest.approx(5.0)), + (Trace, np.array([[1, 4, 1], [2, 4, 2], [3, 4, 3]]), pytest.approx(8.0)), ], ) def test_unary_evaluation(op, a, result): From 933f5829b7ff80747c290a80e4b8f0399d021b32 Mon Sep 17 00:00:00 2001 From: PaulSharp <44529197+DrPaulSharp@users.noreply.github.com> Date: Tue, 17 Mar 2026 11:34:09 +0000 Subject: [PATCH 2/4] Adds matrix inverse operation --- sasdata/quantities/quantity.py | 39 ++++++++++++++++++++++++ test/quantities/utest_math_operations.py | 17 ++++++++++- test/quantities/utest_operations.py | 4 ++- 3 files changed, 58 insertions(+), 2 deletions(-) diff --git a/sasdata/quantities/quantity.py b/sasdata/quantities/quantity.py index 8088d488a..6f0ab86e5 100644 --- a/sasdata/quantities/quantity.py +++ b/sasdata/quantities/quantity.py @@ -52,6 +52,18 @@ def trace(a: Union["Quantity[ArrayLike]", ArrayLike], offset: int = 0, axis1: in return np.trace(a, offset, axis1, axis2) +def matinv(a: Union["Quantity[ArrayLike]", ArrayLike]): + """Find the inverse of a matrix.""" + if isinstance(a, Quantity): + return DerivedQuantity( + value=np.linalg.inv(a.value), + units=a.units, + history=QuantityHistory.apply_operation(MatInv, a.history), + ) + else: + return np.linalg.inv(a) + + def determinant(a: Union["Quantity[ArrayLike]", ArrayLike]): """Find the determinant of an array or an array based quantity.""" if isinstance(a, Quantity): @@ -1180,6 +1192,32 @@ def _clean_ab(self, a, b): return MatMul(a, b) +class MatInv(UnaryOperation): + """Matrix inversion, using numpy""" + + serialisation_name = "matinv" + + def evaluate(self, variables: dict[int, T]) -> T: + return np.linalg.inv(self.a.evaluate(variables)) + + def _derivative(self, hash_value: int) -> Operation: + return Neg(Matmul(MatMul(MatInv(self.a), self.a._derivative(hash_value)), MatInv(self.a))) + + def _clean(self): + clean_a = self.a._clean() + + if isinstance(clean_a, MatInv): + # Removes double inversions + return clean_a.a + + elif isinstance(clean_a, Determinant): + # Inverse of determinant is determinant of inverse + return Determinant(MatInv(clean_a.a)) + + else: + return MatInv(clean_a) + + class TensorDot(Operation): serialisation_name = "tensor_product" @@ -1236,6 +1274,7 @@ def _deserialise(parameters: dict) -> "Operation": Determinant, Dot, MatMul, + MatInv, TensorDot, ] diff --git a/test/quantities/utest_math_operations.py b/test/quantities/utest_math_operations.py index 3b38bad55..11d0b11b3 100644 --- a/test/quantities/utest_math_operations.py +++ b/test/quantities/utest_math_operations.py @@ -4,7 +4,7 @@ import pytest from sasdata.quantities import units -from sasdata.quantities.quantity import NamedQuantity, tensordot, trace, transpose +from sasdata.quantities.quantity import NamedQuantity, matinv, tensordot, trace, transpose order_list = [[0, 1, 2, 3], [0, 2, 1], [1, 0], [0, 1], [2, 0, 1], [3, 1, 2, 0]] @@ -57,6 +57,21 @@ def test_trace_axes(matrix, axis1, axis2, expected_trace): ).all() +@pytest.mark.parametrize( + "matrix, expected_inverse", + [ + (np.array([[1]]), np.array([[1]])), + (np.array([[-2.0, 1.0], [1.5, -0.5]]), np.array([[1, 2], [3, 4]])), + ], +) +def test_inverse(matrix, expected_inverse): + """Check that the matinv operation correctly inverse for raw data and quantities.""" + print(matinv(matrix)) + print(expected_inverse) + assert (matinv(matrix) == expected_inverse).all() + assert (matinv(NamedQuantity("testmat", matrix, units=units.none)).value == expected_inverse).all() + + rng_seed = 1979 tensor_product_with_identity_sizes = (4, 6, 5) diff --git a/test/quantities/utest_operations.py b/test/quantities/utest_operations.py index 586c1ba26..d564b044b 100644 --- a/test/quantities/utest_operations.py +++ b/test/quantities/utest_operations.py @@ -18,6 +18,7 @@ Inv, Ln, Log, + MatInv, MatMul, Mul, MultiplicativeIdentity, @@ -48,7 +49,7 @@ ) -@pytest.fixture(params=[Determinant, Inv, Exp, Ln, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose]) +@pytest.fixture(params=[Determinant, Inv, Exp, Ln, MatInv, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose]) def unary_operation(request): return request.param(x) @@ -232,6 +233,7 @@ def test_derivative(op, result): [ (Neg(Neg(x))), (Inv(Inv(x))), + (MatInv(MatInv(x))), ], ) def test_clean_double_applications(op): From 51bffaf2327bcaf17133afd3b6311bede8230f8b Mon Sep 17 00:00:00 2001 From: PaulSharp <44529197+DrPaulSharp@users.noreply.github.com> Date: Thu, 19 Mar 2026 14:13:18 +0000 Subject: [PATCH 3/4] Adds norm operations --- sasdata/quantities/quantity.py | 184 +++++++++++++++++++++++++++- test/quantities/utest_operations.py | 67 +++++++++- 2 files changed, 248 insertions(+), 3 deletions(-) diff --git a/sasdata/quantities/quantity.py b/sasdata/quantities/quantity.py index 6f0ab86e5..7f56900c2 100644 --- a/sasdata/quantities/quantity.py +++ b/sasdata/quantities/quantity.py @@ -76,6 +76,48 @@ def determinant(a: Union["Quantity[ArrayLike]", ArrayLike]): return np.linalg.det(a) +def norm_1(a: Union["Quantity[ArrayLike]", ArrayLike], axes: int | tuple[int] | None = None): + """Caculate the 1-norm of an array or an array based quantity.""" + if isinstance(a, Quantity): + if axes is None: + return DerivedQuantity( + value=np.linalg.norm(a.value, ord=1, axes=axes), + units=a.units, + history=QuantityHistory.apply_operation(Norm_1, a.history), + ) + + else: + return DerivedQuantity( + value=np.linalg.norm(a.value, ord=1, axes=axes), + units=a.units, + history=QuantityHistory.apply_operation(Norm_1, a.history, axes=axes), + ) + + else: + return np.linalg.norm(a.value, ord=1, axes=axes) + + +def norm_2(a: Union["Quantity[ArrayLike]", ArrayLike], axes: int | tuple[int] | None = None): + """Caculate the 2-norm of an array or an array based quantity.""" + if isinstance(a, Quantity): + if axes is None: + return DerivedQuantity( + value=np.linalg.norm(a.value, axes=axes), + units=a.units, + history=QuantityHistory.apply_operation(Norm_2, a.history), + ) + + else: + return DerivedQuantity( + value=np.linalg.norm(a.value, axes=axes), + units=a.units, + history=QuantityHistory.apply_operation(Norm_2, a.history, axes=axes), + ) + + else: + return np.linalg.norm(a.value, axes=axes) + + def dot(a: Union["Quantity[ArrayLike]", ArrayLike], b: Union["Quantity[ArrayLike]", ArrayLike]): """Dot product of two arrays or two array based quantities""" a_is_quantity = isinstance(a, Quantity) @@ -331,6 +373,34 @@ def __eq__(self, other): return False +class MatrixIdentity(ConstantBase): + serialisation_name = "identity" + + def __init__(self, size): + self.size = size + + def evaluate(self, variables: dict[int, T]) -> T: + return np.eye(self.size) + + def _derivative(self, hash_value: int): + return Constant(np.zeros((self.size, self.size))) + + @staticmethod + def _deserialise(parameters: dict) -> "Operation": + return MatrixIdentity(self.size) + + def _serialise_parameters(self) -> dict[str, Any]: + return {"size": numerical_encode(self.size)} + + def summary(self, indent_amount: int = 0, indent=" "): + return f"{indent_amount * indent}{self.size} [Matrix Id.]" + + def __eq__(self, other): + if isinstance(other, MatrixIdentity): + return self.size == other.size + return False + + class Constant(ConstantBase): serialisation_name = "constant" @@ -1071,7 +1141,7 @@ def summary(self, indent_amount: int = 0, indent=" "): f"{indent_amount * indent}Transpose(\n" + self.a.summary(indent_amount + 1, indent) + "\n" - + f"{(indent_amount + 1) * indent}{self.axes}\n" + + f"{(indent_amount + 1) * indent}{list(self.axes)}\n" + f"{indent_amount * indent})" ) @@ -1248,6 +1318,116 @@ def _deserialise(parameters: dict) -> "Operation": ) +class Norm_1(Operation): + """1-norm of a matrix from numpy""" + + serialisation_name = "norm_1" + + def __init__(self, a: Operation, axes: int | tuple[int] | None = None): + self.a = a + self.axes = axes + + def evaluate(self, variables: dict[int, T]) -> T: + return np.linalg.norm(self.a.evaluate(variables), ord=1, axis=self.axes) + + def _derivative(self, hash_value: int) -> Operation: + return np.sign(self.a) + + def _clean(self): + clean_a = self.a._clean() + return Norm_1(clean_a, self.axes) + + def _serialise_parameters(self) -> dict[str, Any]: + if self.axes is None: + return {"a": self.a._serialise_json()} + else: + return {"a": self.a._serialise_json(), "axes": list(self.axes)} + + @staticmethod + def _deserialise(parameters: dict) -> "Operation": + if "axes" in parameters: + return Norm_1(a=Operation.deserialise_json(parameters["a"]), axes=tuple(parameters["axes"])) + else: + return Norm_1(a=Operation.deserialise_json(parameters["a"])) + + def summary(self, indent_amount: int = 0, indent=" "): + if self.axes is None: + return ( + f"{indent_amount * indent}Norm_1(\n" + + self.a.summary(indent_amount + 1, indent) + + "\n" + + f"{indent_amount * indent})" + ) + else: + return ( + f"{indent_amount * indent}Norm_1(\n" + + self.a.summary(indent_amount + 1, indent) + + "\n" + + f"{(indent_amount + 1) * indent}{list(self.axes)}\n" + + f"{indent_amount * indent})" + ) + + def __eq__(self, other): + if isinstance(other, Norm_1): + return other.a == self.a + return False + + +class Norm_2(Operation): + """2-norm of a matrix from numpy""" + + serialisation_name = "norm_2" + + def __init__(self, a: Operation, axes: int | tuple[int] | None = None): + self.a = a + self.axes = axes + + def evaluate(self, variables: dict[int, T]) -> T: + return np.linalg.norm(self.a.evaluate(variables), axis=self.axes) + + def _derivative(self, hash_value: int) -> Operation: + return Transpose(Div(self.a, Norm_2(self.a, self.axes))) + + def _clean(self): + clean_a = self.a._clean() + return Norm_2(clean_a, self.axes) + + def _serialise_parameters(self) -> dict[str, Any]: + if self.axes is None: + return {"a": self.a._serialise_json()} + else: + return {"a": self.a._serialise_json(), "axes": list(self.axes)} + + @staticmethod + def _deserialise(parameters: dict) -> "Operation": + if "axes" in parameters: + return Norm_2(a=Operation.deserialise_json(parameters["a"]), axes=tuple(parameters["axes"])) + else: + return Norm_2(a=Operation.deserialise_json(parameters["a"])) + + def summary(self, indent_amount: int = 0, indent=" "): + if self.axes is None: + return ( + f"{indent_amount * indent}Norm_2(\n" + + self.a.summary(indent_amount + 1, indent) + + "\n" + + f"{indent_amount * indent})" + ) + else: + return ( + f"{indent_amount * indent}Norm_2(\n" + + self.a.summary(indent_amount + 1, indent) + + "\n" + + f"{(indent_amount + 1) * indent}{list(self.axes)}\n" + + f"{indent_amount * indent})" + ) + + def __eq__(self, other): + if isinstance(other, Norm_2): + return other.a == self.a + return False + + _serialisable_classes = [ AdditiveIdentity, MultiplicativeIdentity, @@ -1276,6 +1456,8 @@ def _deserialise(parameters: dict) -> "Operation": MatMul, MatInv, TensorDot, + Norm_1, + Norm_2, ] _serialisation_lookup = {class_.serialisation_name: class_ for class_ in _serialisable_classes} diff --git a/test/quantities/utest_operations.py b/test/quantities/utest_operations.py index d564b044b..46a746c07 100644 --- a/test/quantities/utest_operations.py +++ b/test/quantities/utest_operations.py @@ -20,9 +20,12 @@ Log, MatInv, MatMul, + MatrixIdentity, Mul, MultiplicativeIdentity, Neg, + Norm_1, + Norm_2, Operation, Pow, Sin, @@ -64,6 +67,16 @@ def log_pow_operation(request): return request.param(x, 2) +@pytest.fixture(params=[Transpose, Norm_1, Norm_2]) +def axis_operation(request): + return request.param(x, (0,)) + + +@pytest.fixture(params=[Transpose, Norm_1, Norm_2]) +def axis_none_operation(request): + return request.param(x, None) + + def test_serialise_deserialise(): serialised = operation_with_everything.serialise() deserialised = Operation.deserialise(serialised) @@ -96,6 +109,22 @@ def test_log_pow_serialise_deserialise(log_pow_operation): assert serialised == reserialised +def test_axis_serialise_deserialise(axis_operation): + serialised = axis_operation.serialise() + deserialised = Operation.deserialise(serialised) + reserialised = deserialised.serialise() + + assert serialised == reserialised + + +def test_axis_none_serialise_deserialise(axis_none_operation): + serialised = axis_none_operation.serialise() + deserialised = Operation.deserialise(serialised) + reserialised = deserialised.serialise() + + assert serialised == reserialised + + def test_trace_serialise_deserialise(): serialised = Trace(x).serialise() deserialised = Operation.deserialise(serialised) @@ -113,6 +142,11 @@ def test_summary(op, summary): assert f.summary() == summary +def test_matrix_id_summary(): + f = MatrixIdentity(1) + assert f.summary() == "1 [Matrix Id.]" + + def test_variable_summary(): assert x.summary() == "x" @@ -129,6 +163,14 @@ def test_log_pow_summary(log_pow_operation): assert log_pow_operation.summary() == f"{log_pow_operation.__class__.__name__}(\n x\n 2\n)" +def test_axis_summary(axis_operation): + assert axis_operation.summary() == f"{axis_operation.__class__.__name__}(\n x\n [0]\n)" + + +def test_axis_none_summary(axis_none_operation): + assert axis_none_operation.summary() == f"{axis_none_operation.__class__.__name__}(\n x\n)" + + def test_trace_summary(): op = Trace(x) assert op.summary() == f"{op.__class__.__name__}(\n x\n 0\n 0\n 1\n)" @@ -212,10 +254,31 @@ def test_matmul_evaluation(op, a, b, result): @pytest.mark.parametrize( "op, a, result", - [(Transpose, np.array([[1, 2]]), np.array([[1], [2]])), (Transpose, [[1, 2], [3, 4]], [[1, 3], [2, 4]])], + [ + (Transpose, np.array([[1, 2]]), np.array([[1], [2]])), + (Transpose, [[1, 2], [3, 4]], [[1, 3], [2, 4]]), + (Norm_1, [[1, 2], [3, 4]], np.float64(6.0)), + (Norm_2, [[1, 2], [3, 4]], np.float64(np.sqrt(30.0))), + (Norm_2, [[1.0, 2.5], [3.0, 4.0]], np.float64(np.sqrt(32.25))), + ], ) -def test_transpose_evaluation(op, a, result): +def test_axis_none_evaluation(op, a, result): f = op(Constant(a)) + print(f.evaluate({})) + print(result) + assert (f.evaluate({}) == result).all() + + +@pytest.mark.parametrize( + "op, a, axes, result", + [ + (Transpose, [[1, 2], [3, 4]], (1, 0), [[1, 3], [2, 4]]), + (Norm_1, np.array([[1, 2], [3, 4]]), 1, np.array([3.0, 7.0])), + (Norm_2, np.array([[1, 2], [3, 4]]), 1, np.array([np.sqrt(5.0), 5.0])), + ], +) +def test_axis_evaluation(op, a, axes, result): + f = op(Constant(a), axes=axes) assert (f.evaluate({}) == result).all() From b72f6fa82c50fba6e3fc6eddaa68a12262f05f66 Mon Sep 17 00:00:00 2001 From: PaulSharp <44529197+DrPaulSharp@users.noreply.github.com> Date: Fri, 20 Mar 2026 19:22:22 +0000 Subject: [PATCH 4/4] Removes determinant and sympy dependency --- requirements.txt | 1 - sasdata/quantities/quantity.py | 34 ----------------------------- test/quantities/utest_operations.py | 7 +----- 3 files changed, 1 insertion(+), 41 deletions(-) diff --git a/requirements.txt b/requirements.txt index 623acaac1..64bfe29a4 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,7 +5,6 @@ lxml # Calculation numpy scipy -sympy # Unit testing pytest diff --git a/sasdata/quantities/quantity.py b/sasdata/quantities/quantity.py index 7f56900c2..6f3f26410 100644 --- a/sasdata/quantities/quantity.py +++ b/sasdata/quantities/quantity.py @@ -5,7 +5,6 @@ import h5py import numpy as np -import sympy as sp from numpy._typing import ArrayLike from sasdata.quantities import units @@ -64,18 +63,6 @@ def matinv(a: Union["Quantity[ArrayLike]", ArrayLike]): return np.linalg.inv(a) -def determinant(a: Union["Quantity[ArrayLike]", ArrayLike]): - """Find the determinant of an array or an array based quantity.""" - if isinstance(a, Quantity): - return DerivedQuantity( - value=np.linalg.det(a.value), - units=a.units, # TODO: hmmm, need to think about this - history=QuantityHistory.apply_operation(Determinant, a.history), - ) - else: - return np.linalg.det(a) - - def norm_1(a: Union["Quantity[ArrayLike]", ArrayLike], axes: int | tuple[int] | None = None): """Caculate the 1-norm of an array or an array based quantity.""" if isinstance(a, Quantity): @@ -1201,22 +1188,6 @@ def __eq__(self, other): return False -class Determinant(UnaryOperation): - """Array Determinant - backed by numpy's det method""" - - serialisation_name = "determinant" - - def evaluate(self, variables: dict[int, T]) -> T: - return np.linalg.det(self.a.evaluate(variables)) - - def _derivative(self, hash_value: int) -> Operation: - return Trace(sp.adjugate(self.a) * self.a._derivative(hash_value)) - - def _clean(self): - clean_a = self.a._clean() - return Determinant(clean_a) # Do nothing for now - can consider identity matrix and common columns - - class Dot(BinaryOperation): """Dot product - backed by numpy's dot method""" @@ -1280,10 +1251,6 @@ def _clean(self): # Removes double inversions return clean_a.a - elif isinstance(clean_a, Determinant): - # Inverse of determinant is determinant of inverse - return Determinant(MatInv(clean_a.a)) - else: return MatInv(clean_a) @@ -1451,7 +1418,6 @@ def __eq__(self, other): Log, Transpose, Trace, - Determinant, Dot, MatMul, MatInv, diff --git a/test/quantities/utest_operations.py b/test/quantities/utest_operations.py index 46a746c07..77a9faa54 100644 --- a/test/quantities/utest_operations.py +++ b/test/quantities/utest_operations.py @@ -11,7 +11,6 @@ ArcTan, Constant, Cos, - Determinant, Div, Dot, Exp, @@ -52,7 +51,7 @@ ) -@pytest.fixture(params=[Determinant, Inv, Exp, Ln, MatInv, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose]) +@pytest.fixture(params=[Inv, Exp, Ln, MatInv, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose]) def unary_operation(request): return request.param(x) @@ -205,8 +204,6 @@ def test_evaluation(op, result): (ArcCos, -1.0, math.pi), (ArcTan, 0.0, 0.0), (ArcTan, -1.0, -0.25 * math.pi), - (Determinant, np.array([[1, 2], [3, 4]]), pytest.approx(-2.0)), - (Determinant, np.array([[1, 4, 1], [2, 4, 2], [3, 4, 3]]), pytest.approx(0.0)), (Trace, np.array([[1, 2], [3, 4]]), pytest.approx(5.0)), (Trace, np.array([[1, 4, 1], [2, 4, 2], [3, 4, 3]]), pytest.approx(8.0)), ], @@ -264,8 +261,6 @@ def test_matmul_evaluation(op, a, b, result): ) def test_axis_none_evaluation(op, a, result): f = op(Constant(a)) - print(f.evaluate({})) - print(result) assert (f.evaluate({}) == result).all()