Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sirt documentation improvements #1959

Merged
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 16 additions & 16 deletions Wrappers/Python/cil/optimisation/algorithms/SIRT.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,23 +27,22 @@


class SIRT(Algorithm):

r"""Simultaneous Iterative Reconstruction Technique, see :cite:`Kak2001`.

Simultaneous Iterative Reconstruction Technique (SIRT) solves
the following problem

.. math:: A x = b

The SIRT algorithm is
The SIRT update step for iteration :math:`k` is given by

.. math:: x^{k+1} = \mathrm{proj}_{C}( x^{k} + \omega * D ( A^{T} ( M * (b - Ax^{k}) ) ) ),
.. math:: x^{k+1} = \mathrm{proj}_{C}( x^{k} + \omega D ( A^{T} ( M (b - Ax^{k}) ) ) ),
hrobarts marked this conversation as resolved.
Show resolved Hide resolved

where,
:math:`M = \frac{1}{A*\mathbb{1}}`,
:math:`M = \frac{1}{A\mathbb{1}}`,
:math:`D = \frac{1}{A^{T}\mathbb{1}}`,
:math:`\mathbb{1}` is a :code:`DataContainer` of ones,
:math:`\mathrm{prox}_{C}` is the projection over a set :math:`C`,
:math:`\mathrm{proj}_{C}` is the projection over a set :math:`C`,
and :math:`\omega` is the relaxation parameter.

Parameters
Expand All @@ -63,9 +62,6 @@ class SIRT(Algorithm):
A function with :code:`proximal` method, e.g., :class:`.IndicatorBox` function and :meth:`.IndicatorBox.proximal`,
or :class:`.TotalVariation` function and :meth:`.TotalVariation.proximal`.

kwargs:
Keyword arguments used from the base class :class:`.Algorithm`.

Note
----
If :code:`constraint` is not passed, :code:`lower` and :code:`upper` are used to create an :class:`.IndicatorBox` and apply its :code:`proximal`.
Expand All @@ -77,21 +73,22 @@ class SIRT(Algorithm):

The preconditioning arrays (weights) :code:`M` and :code:`D` used in SIRT are defined as

.. math:: M = \frac{1}{A*\mathbb{1}} = \frac{1}{\sum_{j}a_{i,j}}

.. math:: D = \frac{1}{A*\mathbb{1}} = \frac{1}{\sum_{i}a_{i,j}}
.. math:: M = \frac{1}{A\mathbb{1}}
.. math:: D = \frac{1}{A^T\mathbb{1}}


Examples
--------
.. math:: \underset{x}{\mathrm{argmin}} \frac{1}{2}\| x - d\|^{2}
.. math:: \underset{x}{\mathrm{argmin}} \frac{1}{2}\| Ax - d\|^{2}

>>> sirt = SIRT(initial = ig.allocate(0), operator = A, data = d, max_iteration = 5)
>>> sirt = SIRT(initial = ig.allocate(0), operator = A, data = d)

"""


def __init__(self, initial=None, operator=None, data=None, lower=None, upper=None, constraint=None, **kwargs):
"""Constructor of SIRT algorithm"""

super(SIRT, self).__init__(**kwargs)

Expand Down Expand Up @@ -140,10 +137,12 @@ def set_up(self, initial, operator, data, lower=None, upper=None, constraint=Non

@property
def relaxation_parameter(self):
"""Get the relaxation parameter :math:`\omega`"""
return self._relaxation_parameter

@property
def D(self):
"""Get the preconditioning array :math:`D`"""
return self._Dscaled / self._relaxation_parameter

def set_relaxation_parameter(self, value=1.0):
Expand All @@ -164,6 +163,7 @@ def set_relaxation_parameter(self, value=1.0):


def _set_up_weights(self):
"""Set up the preconditioning arrays M and D"""
self.M = 1./self.operator.direct(self.operator.domain_geometry().allocate(value=1.0))
self._Dscaled = 1./self.operator.adjoint(self.operator.range_geometry().allocate(value=1.0))

Expand Down Expand Up @@ -196,9 +196,9 @@ def _remove_nan_or_inf(self, datacontainer, replace_with=1.0):

def update(self):

r""" Performs a single iteration of the SIRT algorithm
r""" Performs a single iteration of the SIRT algorithm. The update step for iteration :math:`k` is given by

.. math:: x^{k+1} = \mathrm{proj}_{C}( x^{k} + \omega * D ( A^{T} ( M * (b - Ax) ) ) )
.. math:: x^{k+1} = \mathrm{proj}_{C}( x^{k} + \omega D ( A^{T} ( M (b - Ax^{k}) ) ) )

"""

Expand All @@ -218,7 +218,7 @@ def update(self):
self.x=self.constraint.proximal(self.x, tau=1)

def update_objective(self):
r"""Returns the objective
r""" Appends the current objective value to the list of previous objective values

.. math:: \frac{1}{2}\|A x - b\|^{2}

Expand Down
Loading