From 1128c83867145557af9c17125f687976d35e8fbe Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Thu, 29 Jan 2026 16:15:00 +0000 Subject: [PATCH 1/3] initial draft --- .../operators/DiagonalOperator.py | 54 ++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py b/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py index 93c9cd9eda..8e406f5fe4 100644 --- a/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py +++ b/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py @@ -18,7 +18,8 @@ import numpy as np from cil.framework import ImageData -from cil.optimisation.operators import LinearOperator +from cil.optimisation.operators import LinearOperator, BlockOperator +from cil.framework import BlockDataContainer class DiagonalOperator(LinearOperator): @@ -34,6 +35,56 @@ class DiagonalOperator(LinearOperator): If the :class:`DataContainer` `x` is also flattened, we have a :math:`M*N` vector. Now, matrix-vector multiplcation is allowed and results to a :math:`(M*N,1)` vector. After reshaping we recover a :math:`M\times N` :class:`DataContainer`. + Parameters + ---------- + diagonal : DataContainer + DataContainer with the same dimensions as the data to be operated on. + domain_geometry : ImageGeometry + Specifies the geometry of the operator domain. If 'None' will use the diagonal geometry directly. default=None . + + """ + def __init__(self, diagonal, domain_geometry=None): + if domain_geometry is None: + domain_geometry = diagonal.geometry.copy() + super(DiagonalOperator, self).__init__(domain_geometry=domain_geometry, + range_geometry=domain_geometry) + if isinstance(diagonal, BlockDataContainer): + operator_list=[] + for i in range(len(diagonal)): + operator_list.append(DiagonalOperator(diagonal[i], diagonal[i].geometry)) + self.operator = BlockOperator(operator_list) + else: + self.operator = _DiagonalOperator(diagonal, domain_geometry) + def direct(self,x,out=None): + "Returns :math:`D\circ x` " + return self.operator.direct(x,out=out) + + def adjoint(self,x, out=None): + "Returns :math:`D^*\circ x` " + return self.operator.adjoint(x,out=out) + + def calculate_norm(self, **kwargs): + r""" Returns the operator norm of DiagonalOperator which is the :math:`\infty` norm of `diagonal` + + .. math:: \|D\|_{\infty} = \max_{i}\{|D_{i}|\} + """ + return self.operator.calculate_norm(**kwargs) + + +class _DiagonalOperator(LinearOperator): + + r"""DiagonalOperator + + Performs an element-wise multiplication, i.e., `Hadamard Product `_ + of a :class:`DataContainer` `x` and :class:`DataContainer` `diagonal`, `d` . + + .. math:: (D\circ x) = \sum_{i,j}^{M,N} D_{i,j} x_{i, j} + + In matrix-vector interpretation, if `D` is a :math:`M\times N` dense matrix and is flattened, we have a :math:`M*N \times M*N` vector. + A sparse diagonal matrix, i.e., :class:`DigaonalOperator` can be created if we add the vector above to the main diagonal. + If the :class:`DataContainer` `x` is also flattened, we have a :math:`M*N` vector. + Now, matrix-vector multiplcation is allowed and results to a :math:`(M*N,1)` vector. After reshaping we recover a :math:`M\times N` :class:`DataContainer`. + Parameters ---------- diagonal : DataContainer @@ -67,3 +118,4 @@ def calculate_norm(self, **kwargs): .. math:: \|D\|_{\infty} = \max_{i}\{|D_{i}|\} """ return self.diagonal.abs().max() + From 488dd169ff5da0143aed1327ae29a24d4c57f318 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Thu, 29 Jan 2026 16:43:03 +0000 Subject: [PATCH 2/3] attempt 2 with a seperate _BlockDiagonalOperator --- .../operators/DiagonalOperator.py | 60 +++++++++++++++++-- 1 file changed, 55 insertions(+), 5 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py b/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py index 8e406f5fe4..442dc5ae2a 100644 --- a/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py +++ b/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py @@ -49,10 +49,7 @@ def __init__(self, diagonal, domain_geometry=None): super(DiagonalOperator, self).__init__(domain_geometry=domain_geometry, range_geometry=domain_geometry) if isinstance(diagonal, BlockDataContainer): - operator_list=[] - for i in range(len(diagonal)): - operator_list.append(DiagonalOperator(diagonal[i], diagonal[i].geometry)) - self.operator = BlockOperator(operator_list) + self.operator = _BlockDiagonalOperator(diagonal) else: self.operator = _DiagonalOperator(diagonal, domain_geometry) def direct(self,x,out=None): @@ -96,7 +93,7 @@ class _DiagonalOperator(LinearOperator): def __init__(self, diagonal, domain_geometry=None): if domain_geometry is None: domain_geometry = diagonal.geometry.copy() - super(DiagonalOperator, self).__init__(domain_geometry=domain_geometry, + super(_DiagonalOperator, self).__init__(domain_geometry=domain_geometry, range_geometry=domain_geometry) self.diagonal = diagonal @@ -119,3 +116,56 @@ def calculate_norm(self, **kwargs): """ return self.diagonal.abs().max() +class _BlockDiagonalOperator(LinearOperator): + + r"""BlockDiagonalOperator + + Performs an element-wise multiplication, i.e., `Hadamard Product `_ + of a :class:`DataContainer` `x` and :class:`DataContainer` `diagonal`, `d` . + + .. math:: (D\circ x) = \sum_{i,j}^{M,N} D_{i,j} x_{i, j} + + In matrix-vector interpretation, if `D` is a :math:`M\times N` dense matrix and is flattened, we have a :math:`M*N \times M*N` vector. + A sparse diagonal matrix, i.e., :class:`DigaonalOperator` can be created if we add the vector above to the main diagonal. + If the :class:`DataContainer` `x` is also flattened, we have a :math:`M*N` vector. + Now, matrix-vector multiplcation is allowed and results to a :math:`(M*N,1)` vector. After reshaping we recover a :math:`M\times N` :class:`DataContainer`. + + Parameters + ---------- + diagonal : DataContainer + DataContainer with the same dimensions as the data to be operated on. + domain_geometry : ImageGeometry + Specifies the geometry of the operator domain. If 'None' will use the diagonal geometry directly. default=None . + + """ + def __init__(self, diagonal, domain_geometry=None): + if domain_geometry is None: + domain_geometry = diagonal.geometry.copy() + super(_DiagonalOperator, self).__init__(domain_geometry=domain_geometry, + range_geometry=domain_geometry) + self.diagonal = diagonal + self.diagonal_operator_list = [ DiagonalOperator(diagonal[i]) for i in range(len(diagonal)) ] + + def direct(self,x,out=None): + "Returns :math:`D\circ x` " + if out is None: + out = x.copy() + for i in range(len(self.diagonal)): + self.diagonal_operator_list[i].direct(x[i], out=out[i]) + return out + + def adjoint(self,x, out=None): + "Returns :math:`D^*\circ x` " + if out is None: + out = x.copy() + for i in range(len(self.diagonal)): + self.diagonal_operator_list[i].adjoint(x[i], out=out[i]) + return out + + def calculate_norm(self, **kwargs): + r""" Returns the operator norm of DiagonalOperator which is the :math:`\infty` norm of `diagonal` + + .. math:: \|D\|_{\infty} = \max_{i}\{|D_{i}|\} + """ + norms = [ op.calculate_norm(**kwargs) for op in self.diagonal_operator_list ] + return max(norms) \ No newline at end of file From 89c9ab2ce30d00bbdc19ca35175ab016ea356f35 Mon Sep 17 00:00:00 2001 From: Margaret Duff Date: Thu, 29 Jan 2026 16:50:14 +0000 Subject: [PATCH 3/3] Added in diagonal property and fixed a couple of bugs --- .../cil/optimisation/operators/DiagonalOperator.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py b/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py index 442dc5ae2a..cab5298fb2 100644 --- a/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py +++ b/Wrappers/Python/cil/optimisation/operators/DiagonalOperator.py @@ -45,13 +45,14 @@ class DiagonalOperator(LinearOperator): """ def __init__(self, diagonal, domain_geometry=None): if domain_geometry is None: - domain_geometry = diagonal.geometry.copy() + domain_geometry = diagonal.geometry super(DiagonalOperator, self).__init__(domain_geometry=domain_geometry, range_geometry=domain_geometry) if isinstance(diagonal, BlockDataContainer): self.operator = _BlockDiagonalOperator(diagonal) else: self.operator = _DiagonalOperator(diagonal, domain_geometry) + self.diagonal = diagonal def direct(self,x,out=None): "Returns :math:`D\circ x` " return self.operator.direct(x,out=out) @@ -132,8 +133,8 @@ class _BlockDiagonalOperator(LinearOperator): Parameters ---------- - diagonal : DataContainer - DataContainer with the same dimensions as the data to be operated on. + diagonal : BlockDataContainer + BlockDataContainer with the same dimensions as the data to be operated on. domain_geometry : ImageGeometry Specifies the geometry of the operator domain. If 'None' will use the diagonal geometry directly. default=None . @@ -141,7 +142,7 @@ class _BlockDiagonalOperator(LinearOperator): def __init__(self, diagonal, domain_geometry=None): if domain_geometry is None: domain_geometry = diagonal.geometry.copy() - super(_DiagonalOperator, self).__init__(domain_geometry=domain_geometry, + super(_BlockDiagonalOperator, self).__init__(domain_geometry=domain_geometry, range_geometry=domain_geometry) self.diagonal = diagonal self.diagonal_operator_list = [ DiagonalOperator(diagonal[i]) for i in range(len(diagonal)) ]