AbsFunctions added

NOTICE.txt

@@ -70,8 +70,10 @@ Sam Porter (2024) - 5
 Joshua Hellier (2024) - 3
 Nicholas Whyatt (2024) - 1
 Rasmia Kulan (2024) - 1
+Francis M Watson (2024) - 3 
 Emmanuel Ferdman (2025) - 14
 
+
 CIL Advisory Board:
 Llion Evans - 9
 William Lionheart - 3

Wrappers/Python/cil/optimisation/functions/AbsFunction.py

@@ -0,0 +1,227 @@
+#  Copyright 2024 United Kingdom Research and Innovation
+#  Copyright 2024 The University of Manchester
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+#
+# Authors:
+# CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
+#
+# This work has been supported by the Royal Academy of Engineering and the
+# Office of the Chief Science Adviser for National Security under the UK
+# Intelligence Community Postdoctoral Research Fellowship programme.
+#
+# Francis M Watson, University of Manchester 2024
+#
+
+
+import numpy as np
+from cil.optimisation.functions import Function
+from cil.framework import DataContainer
+from typing import Optional
+import warnings
+import logging
+
+log = logging.getLogger(__name__)
+
+class FunctionOfAbs(Function):
+    r'''A function which acts on the absolute value of the complex-valued input, 
+
+    .. math::  G(z) = H(abs(z))
+
+    This function is initialised with another CIL function, :math:`H` in the above formula. When this function is called, first the absolute value of the input is taken, and then the input is passed to the provided function.
+
+    Included in this class is the proximal map for FunctionOfAbs.  From this, the proximal conjugate is also available from the parent CIL Function class, which is valid for this function. 
+    In the case that :math:`H` is lower semi-continuous, convex, non-decreasing and finite at the origin, (and thus `assume_lower_semi` is set to `True` in the `init`) the convex conjugate is also defined.  
+    The gradient is not defined for this function.
+
+
+
+    Parameters
+    ----------
+    function : Function
+        Function acting on a real input, :math:`H` in the above formula.
+    assume_lower_semi : bool, default False
+        If True, assume that the function is lower semi-continuous, convex, non-decreasing and finite at the origin.
+        This allows the convex conjugate to be calculated as the monotone conjugate, which is less than or equal to the convex conjugate.
+        If False, the convex conjugate returned as 0. This is to ensure compatibility with Algorithms such as PDHG.
+    precision : str, default 'double'
+        Precision of the calculation, 'single' or 'double'
+        Some complex-valued imaging problems involve high dynamic range and/or require fine phase accuracy, necessitating the use of double precision (default).
+        For example, in synthetic aperture radar imagery 100dB+ dynamic range may be encountered, which is more than single precision allows.
+        In other cases use of single precision will reduce memory reservation and may improve performance, depending on the compute architecture.
+
+
+
+
+    Reference
+    ---------
+    For further details see https://doi.org/10.48550/arXiv.2410.22161
+
+    '''
+
+    def __init__(self, function: Function, assume_lower_semi: bool=False, precision: str='double'):
+        self._function = function
+        self._lower_semi = assume_lower_semi
+
+        if precision=='single' or precision=='double':
+            self.precision = precision
+        else: 
+            raise ValueError('Precision must be `single` or `double`')
+
+        super().__init__(L=function.L)
+
+    def __call__(self, x: DataContainer) -> float:
+        call_abs = (_take_abs_input(self.precision))(self._function.__call__)
+        return call_abs(self._function, x)
+
+    def proximal(self, x: DataContainer, tau: float, out: Optional[DataContainer]=None) -> DataContainer:
+        r'''Returns the proximal map of function :math:`\tau G`  evaluated at x
+
+        .. math:: \text{prox}_{\tau G}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\|z - x\|^{2} + \tau G(z)
+
+        This is accomplished by calculating a bounded proximal map and making a change of phase,
+        :math:`prox_G(z) = prox^+_H(r) \circ \Phi` where  :math:`z = r \circ \Phi`, :math:`r = abs(z)`, :math:`\Phi = \exp(i angl(z))`,
+        and :math:`\circ` is element-wise product.  Also define :math:`prox^+` to be the proximal map of :math:`H`  in which the minimisation carried out over the positive orthant.
+
+
+        Parameters
+        ----------
+        x : DataContainer
+            The input to the function
+        tau: scalar
+            The scalar multiplying the function in the proximal map
+        out: return DataContainer, if None a new DataContainer is returned, default None.
+            DataContainer to store the result of the proximal map
+
+
+        Returns
+        -------
+        DataContainer, the proximal map of the function at x with scalar :math:`\tau`.
+
+        '''
+        prox_abs = _abs_and_project(self.precision)(self._function.proximal)
+        return prox_abs(self._function, x, tau=tau, out=out)
+
+    def convex_conjugate(self, x: DataContainer) -> float:
+        r'''
+        Evaluation of the function G* at x, where G* is the convex conjugate of function G,
+
+        .. math:: G^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle - G(x)
+
+        If :math:`H` is lower semi-continuous, convex, non-decreasing 
+        finite at the origin, then :math:`G^*(z*) = H^+(|z*|)`, where the monotone conjugate :math:`g^+` is
+
+        .. math:: H^+(z^*) =sup {(z, z^*) - H(z) : z >= O}
+
+        The monotone conjugate will therefore be less than or equal to the convex conjugate, 
+        since it is taken over a smaller set.  It is not available directly, but may coincide with
+        the convex conjugate, which is therefore the best estimate we have.  This is only valid for
+        real x. In other cases, a general convex conjugate is not available or defined.      
+
+
+        For reference see:  Convex Analysis, R. Tyrrell Rocakfellar, pp110-111.
+
+
+        Parameters
+        ----------
+        x : DataContainer
+            The input to the function
+
+        Returns
+        -------
+        float:
+
+        '''
+
+        if self._lower_semi:
+            conv_abs = (_take_abs_input(self.precision))(self._function.convex_conjugate)
+            return conv_abs(self._function, x)
+        else:
+            warnings.warn('Convex conjugate is not properly for this function, returning 0 for compatibility with optimisation algorithms')
+            return 0.0
+
+    def gradient(self, x):
+        '''Gradient of the function at x is not defined for this function.
+        '''
+        raise NotImplementedError('Gradient not available for this function')
+
+def _take_abs_input(precision='double'):
+    def _take_abs_input_inner(func):
+            '''Decorator for function to act on abs of input of a method'''
+
+            def _take_abs_decorator(self, x: DataContainer, *args, **kwargs):
+                if precision == 'single':
+                    real_dtype = np.float32
+                elif precision == 'double':
+                    real_dtype = np.float64
+                else:
+                    raise ValueError('Precision must be `single` or `double`')
+                rgeo = x.geometry.copy()
+                rgeo.dtype = real_dtype
+                r = rgeo.allocate(0)
+                r.fill(np.abs(x.array).astype(real_dtype))
+                fval = func(r, *args, **kwargs)
+                return fval
+            return _take_abs_decorator
+    return _take_abs_input_inner
+
+def _abs_and_project(precision='double'):
+    def _abs_and_project_inner(func):
+            '''Decorator for function to act on abs of input, 
+            with return being projected to the angle of the input.
+            Requires function return to have the same shape as input,
+            such as prox.'''
+
+
+            def _abs_project_decorator(self, x: DataContainer, *args, **kwargs):
+                if precision == 'single':
+                    real_dtype = np.float32
+                    complex_dtype = np.complex64
+                elif precision == 'double':
+                    real_dtype = np.float64
+                    complex_dtype = np.complex128
+                else:
+                    raise ValueError('Precision must be `single` or `double`')
+                rgeo = x.geometry.copy()
+                rgeo.dtype = real_dtype
+                r = rgeo.allocate(None)
+                r.fill( np.abs(x.array).astype(real_dtype))
+                Phi = np.exp((1j*np.angle(x.array)))
+                out = kwargs.pop('out', None)
+
+                fvals = func(r, *args, **kwargs)
+
+                # Douglas-Rachford splitting to find solution in positive orthant
+                if np.any(fvals.array < 0):
+                    log.info('AbsFunctions: projection to +ve orthant triggered')
+                    cts = 0
+                    y = r.copy()
+                    while np.any(fvals.array < 0):
+                        tmp = fvals.array - 0.5*y.array + 0.5*r.array
+                        tmp[tmp < 0] = 0.
+                        y.array += tmp - fvals.array
+                        fvals = func(y, *args, **kwargs)
+                        cts += 1
+                        if cts > 10:
+                            fvals.array[fvals.array < 0] = 0.
+                            break
+
+                if out is not None:
+                    out.array = fvals.array.astype(complex_dtype)*Phi
+
+                else:
+                    out = x.geometry.allocate(None)
+                    out.array = fvals.array.astype(complex_dtype)*Phi
+                return out
+            return _abs_project_decorator
+    return _abs_and_project_inner

Wrappers/Python/cil/optimisation/functions/__init__.py

@@ -39,4 +39,5 @@
 from .SGFunction import SGFunction
 from .SVRGFunction import SVRGFunction, LSVRGFunction
 from .SAGFunction import SAGFunction, SAGAFunction
+from .AbsFunction import FunctionOfAbs
 

Wrappers/Python/test/test_functions.py

@@ -32,7 +32,7 @@
                                          L1Norm, MixedL21Norm, LeastSquares, \
                                          SmoothMixedL21Norm, OperatorCompositionFunction,\
                                          Rosenbrock, IndicatorBox, TotalVariation, ScaledFunction, SumFunction, SumScalarFunction, \
-                                         WeightedL2NormSquared, MixedL11Norm, ZeroFunction, L1Sparsity
+                                         WeightedL2NormSquared, MixedL11Norm, ZeroFunction, L1Sparsity, FunctionOfAbs
 
 from cil.optimisation.functions import BlockFunction
 
@@ -2124,3 +2124,73 @@ def test_set_num_threads(self):
             N = 10
             ib.set_num_threads(N)
             assert ib.num_threads == N
+
+
+class MockFunction(Function):
+    """Mock function to test FunctionOfAbs"""
+    def __init__(self, L=1.0):
+        super().__init__(L=L)
+
+    def __call__(self, x):
+        return x.array**2  # Simple squared function
+
+    def proximal(self, x, tau, out=None):
+        if out is None:
+            out = x.geometry.allocate(None)
+        out.array = x.array / (1 + tau)  # Simple proximal operator
+        return out
+
+    def convex_conjugate(self, x):
+        return 0.5 * x.array**2  # Quadratic conjugate function
+
+
+class TestFunctionOfAbs(unittest.TestCase):
+    def setUp(self):
+        self.data = VectorData(np.array([1+2j, -3+4j, -5-6j]))
+        self.mock_function = MockFunction()
+        self.abs_function = FunctionOfAbs(self.mock_function)
+        self.abs_function_double = FunctionOfAbs(self.mock_function, precision='double')
+
+
+    def test_initialization(self):
+        self.assertIsInstance(self.abs_function._function, MockFunction)
+        self.assertFalse(self.abs_function._lower_semi)
+
+        self.assertEqual(self.abs_function.real_dtype, np.float32)
+        self.assertEqual(self.abs_function_double.real_dtype, np.float64)
+
+    def test_call(self):
+        result = self.abs_function(self.data)
+        expected = np.abs(self.data.array)**2
+        np.testing.assert_array_almost_equal(result, expected, decimal=4)
+
+        result = self.abs_function_double(self.data)
+        np.testing.assert_array_almost_equal(result, expected, decimal=6)
+
+
+    def test_proximal(self):
+        tau = 0.5
+        result = self.abs_function.proximal(self.data, tau)
+        expected = (np.abs(self.data.array) / (1 + tau)) * np.exp(1j * np.angle(self.data.array))
+        np.testing.assert_array_almost_equal(result.array, expected)
+
+        result = self.abs_function_double.proximal(self.data, tau)
+        np.testing.assert_array_almost_equal(result.array, expected, decimal=4)
+
+    def test_convex_conjugate_lower_semi(self):
+        self.abs_function._lower_semi = True
+        result = self.abs_function.convex_conjugate(self.data)
+        expected = 0.5 * np.abs(self.data.array) ** 2
+        np.testing.assert_array_almost_equal(result, expected, decimal=4)
+
+        self.abs_function_double._lower_semi = True
+        result = self.abs_function_double.convex_conjugate(self.data)
+        np.testing.assert_array_almost_equal(result, expected, decimal=6)
+
+
+    def test_convex_conjugate_not_implemented(self):
+        self.abs_function._lower_semi = False
+        with self.assertRaises(NotImplementedError):
+            self.abs_function.convex_conjugate(self.data)
+
+

Wrappers/Python/test/test_out_in_place.py

@@ -38,7 +38,7 @@
     L1Norm, L2NormSquared, MixedL21Norm, LeastSquares, \
     SmoothMixedL21Norm, OperatorCompositionFunction, \
      IndicatorBox, TotalVariation,  SumFunction, SumScalarFunction, \
-    WeightedL2NormSquared, MixedL11Norm, ZeroFunction
+    WeightedL2NormSquared, MixedL11Norm, ZeroFunction, FunctionOfAbs
 
 from cil.processors import AbsorptionTransmissionConverter, Binner, CentreOfRotationCorrector, MaskGenerator, Masker, Normaliser, Padder, \
 RingRemover, Slicer, TransmissionAbsorptionConverter, PaganinProcessor, FluxNormaliser
@@ -127,9 +127,10 @@ def setUp(self):
             (MixedL11Norm(), bg, True, True, False),
             (BlockFunction(L1Norm(),L2NormSquared()), bg, True, True, False),
             (BlockFunction(L2NormSquared(),L2NormSquared()), bg, True, True, True),
-            (L1Sparsity(WaveletOperator(ig)), ig, True, True, False)
-
-
+            (L1Sparsity(WaveletOperator(ig)), ig, True, True, False),
+            (FunctionOfAbs(TotalVariation(backend='numpy', warm_start= False), assume_lower_semi=True), ig , True, True, False),
+            (FunctionOfAbs(TotalVariation(backend='numpy', warm_start=False), assume_lower_semi=True, precision='double'), ig , True, True, False),
+
         ]
 
         np.random.seed(5)

docs/source/optimisation.rst

@@ -512,6 +512,14 @@ Total variation
    :members:
    :inherited-members:
 
+Function of Absolute Value 
+--------------------------
+
+.. autoclass:: cil.optimisation.functions.FunctionOfAbs
+   :members:
+   :inherited-members:
+
+
 Approximate Gradient base class 
 --------------------------------