From 0e93020b15ca4371478890af1f7912d31969af21 Mon Sep 17 00:00:00 2001 From: Michael-T-McCann Date: Mon, 18 Oct 2021 10:57:57 -0600 Subject: [PATCH 1/5] Copy edit --- docs/source/overview.rst | 154 +++++++++++++++++++++++++++++++-------- 1 file changed, 125 insertions(+), 29 deletions(-) diff --git a/docs/source/overview.rst b/docs/source/overview.rst index 2ea26a1cb..0b30f7875 100644 --- a/docs/source/overview.rst +++ b/docs/source/overview.rst @@ -25,17 +25,40 @@ Overview -`SCICO `__ is a Python package for solving the inverse problems that arise in scientific imaging applications. Its primary focus is providing methods for solving ill-posed inverse problems by using an appropriate prior model of the reconstruction space. SCICO includes a growing suite of operators, cost functionals, regularizers, and optimization routines that may be combined to solve a wide range of problems, and is designed so that it is easy to add new building blocks. +`SCICO `__ is a Python package +for solving the inverse problems that arise +in scientific imaging applications. +Its primary focus is providing methods +for solving ill-posed inverse problems +by using an appropriate prior model of the reconstruction space. +SCICO includes a growing suite of operators, cost +functionals, regularizers, and optimization routines that may be +combined to solve a wide range of problems and is designed so that it is +easy to add new building blocks. SCICO's main advantages are: - - Operators in SCICO are defined using a simple, NumPy-like syntax, yet automatically benefit from GPU/TPU acceleration and `just-in-time compilation `__. - - SCICO can compute the adjoint of linear operators automatically, which saves time when defining new operators. - - SCICO's operator calculus makes code for optimization routines look like the pseudocode in scientific papers. - - SCICO provides state-of-the-art optimization routines, including projected gradients and the Alternating Direction Method of Multipliers, with the flexibility of plug-and-play priors including BM3D and DnCNN denoisers. - - -If you use this library for published work, please cite :cite:`pfister-2021-scico` (see bibtex entry ``pfister-2021-scico`` in `docs/source/references.bib `_ in the source distribution). + - Operators in SCICO are defined using a simple, NumPy-like syntax, + yet automatically benefit from GPU/TPU acceleration and + `just-in-time compilation + `__ + - SCICO can compute the adjoint of linear operators automatically, + which saves time when defining new operators. + - SCICO's operator calculus makes code for optimization routines + look like the pseudocode in scientific papers. + - SCICO provides state-of-the-art optimization routines, + including projected gradients + and the alternating direction method of multipliers (ADMM) + with the flexibility of plug-and-play priors + including BM3D and DnCNN denoisers. + + +If you use this library for published work, +please cite :cite:`pfister-2021-scico` +(see bibtex entry ``pfister-2021-scico`` in +`docs/source/references.bib +`_ +in the source distribution). Anatomy of an Inverse Problem @@ -47,35 +70,84 @@ SCICO is designed to solve problems of the form \argmin_{\mb{x}} \; f(\mb{x}) + \sum_{r=1}^R g_r(C_r (\mb{x})), -with :math:`\mb{x} \in \mathbb{R}^{n}` representing the reconstruction, e.g., a 1D, 2D, 3D, or 3D+time image, :math:`C_r: \mathbb{R}^{n} \to \mathbb{R}^{m_r}` being regularization operators, and :math:`f: \mathbb{R}^{n} \to \mathbb{R}` and :math:`g_r: \mathbb{R}^{m_r} \to \mathbb{R}` being functionals associated with the data fidelity and regularization, respectively. - -In a typical computational imaging problem, we have a `forward model` that models the data acquisition process. We denote this forward model by the (possibly nonlinear) operator :math:`A`. - -We want to find :math:`\mb{x}` such that :math:`\mb{y} \approx A(\mb{x})` -in some appropriate sense. This discrepency is measured in our data fidelity, or `loss`, function - -.. math:: +with :math:`\mb{x} \in \mathbb{R}^{n}` representing the reconstruction, +e.g., a 1D, 2D, 3D, or 3D+time image, +:math:`C_r: \mathbb{R}^{n} \to \mathbb{R}^{m_r}` +being regularization operators, +and :math:`f: \mathbb{R}^{n} \to \mathbb{R}` +and :math:`g_r: \mathbb{R}^{m_r} \to \mathbb{R}` +being functionals associated with the data fidelity +and regularization, respectively. + +In a typical computational imaging problem, +we have a `forward model` that models the data acquisition process. +We denote this forward model +by the (possibly nonlinear) operator :math:`A`. +We want to find :math:`\mb{x}` +such that :math:`\mb{y} \approx A(\mb{x})` +in some appropriate sense. +This discrepency is measured in our data fidelity, or `loss`, function + +.. Math:: f(\mb{x}) = L(A(\mb{x}), \mb{y}) \,, -where :math:`L` is a `functional` that maps a vector (or array) into a scalar. A common choice is :math:`f(\mb{x}) = \norm{\mb{y} - A(\mb{x})}_2^2`. The SCICO :class:`.Loss` object encapsulates the data :math:`\mb{y}`, the forward model :math:`A`, and the functional :math:`L` in a single object. A library of functionals and losses are implemented in :mod:`.functional` and :mod:`.loss`, respectively. - -Note that :math:`f(\mb{x})` can often be interpreted as the negative log likelihood of the candidate :math:`\mb{x}`, given the measurements :math:`\mb{y}` and an underlying noise model. - -The functionals :math:`g_r(C_r (\mb{x}))` are regularization functionals. The :math:`C_r` are operators, usually linear operators. Together, these terms encourage solutions that arex "simple" in some sense. A popular choice is :math:`g = \norm{ \cdot }_1` and :math:`C` to be a :class:`.FiniteDifferece` linear operator; this combination promotes piecewise smooth solutions to the inverse problem. +where :math:`L` is a `functional` that maps a vector (or array) +into a scalar. +A common choice is :math:`f(\mb{x}) = \norm{\mb{y} - A(\mb{x})}_2^2`. +Note that :math:`f(\mb{x})` can often be interpreted +as the negative log likelihood of the candidate :math:`\mb{x}`, +given the measurements :math:`\mb{y}` and an underlying noise model. +The SCICO :class:`.Loss` object encapsulates +the data :math:`\mb{y}`, +the forward model :math:`A`, +and the functional :math:`L` in a single object. +A library of functionals and losses are implemented +in :mod:`.functional` and :mod:`.loss`, respectively. + +The functionals :math:`g_r(C_r (\cdot))` are regularization functionals. +The :math:`C_r` are operators, usually linear operators. +Together, +these terms encourage solutions that are "simple" in some sense. +A popular choice is to let :math:`g = \norm{ \cdot }_1` +and :math:`C` be a :class:`.FiniteDifferece` linear operator; +this combination promotes piecewise smooth solutions +to the inverse problem. Usage Examples -------------- -Usage examples are available as Python scripts and Jupyter Notebooks. Example scripts are located in ``examples/scripts``. The corresponding Jupyter Notebooks are provided in the ``scico-data`` submodule and symlinked to ``examples/notebooks``. They are also viewable on `GitHub `_ and in the documentation :ref:`example_notebooks`. +Usage examples are available as Python scripts and Jupyter Notebooks. +Example scripts are located in ``examples/scripts``. +The corresponding Jupyter Notebooks +are provided in the ``scico-data`` submodule +and symlinked to ``examples/notebooks``. +They are also viewable on +`GitHub `_ +and in the documentation under :ref:`example_notebooks`. Related Projects ---------------- -The SCICO library is inspired by the `GlobalBiolm `_ MATLAB package, which provides a similar object-oriented design for solving computational imaging problems. `Pycsou `_ is a similar Python library for inverse problems that is also inspired by GlobalBioIm. - -A key advantage of SCICO over these libraries is the usage of JAX, which provides automatic hardware acceleration, automatic differentiation, and automatic adjoint calculations. Moreover, as JAX is a machine learning library, state of the art Plug-and-Play regularizers such as DnCNN can specified, trained, and implemented in the same software package. +The SCICO library is inspired by the +`GlobalBiolm `_ +MATLAB package, +which provides a similar object-oriented design +for solving computational imaging problems. +`Pycsou `_ +is a similar Python library for inverse problems +that is also inspired by GlobalBioIm. + +A key advantage of SCICO over these libraries is the usage of +`JAX `_, +which provides +automatic hardware acceleration, +automatic differentiation, +and automatic adjoint calculations. +Moreover, beause JAX is a machine learning library, +it enables implementation and training +of state of the art Plug-and-Play regularizers such as DnCNN. Other related projects that may be of interest include: @@ -85,7 +157,8 @@ Other related projects that may be of interest include: - `ProxImaL `_ - `ProxMin `_ - `ToMoBAR `_ - - `CCPi-Regularisation Toolkit `_ + - `CCPi-Regularisation Toolkit + `_ - `SPORCO `_ - `SigPy `_ - `MIRT `_ @@ -95,14 +168,37 @@ Other related projects that may be of interest include: Contributing ------------ -Bug reports, feature requests, and general suggestions are welcome, and should be submitted via the `github issue system `__. More substantial contributions are also welcome; please see :ref:`scico_dev_contributing`. +Bug reports, feature requests, and general suggestions are welcome, +and should be submitted via the +`github issue system `__. +More substantial contributions are also welcome; +please see :ref:`scico_dev_contributing`. License ------- -SCICO is distributed as open-source software under a BSD 3-Clause License (see the `LICENSE `__ file for details). LANL open source approval reference C20091. +SCICO is distributed as open-source software +under a BSD 3-Clause License +(see the +`LICENSE `__ file +for details). +LANL open source approval reference C20091. © 2020-2021. Triad National Security, LLC. All rights reserved. -This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National Nuclear Security Administration. All rights in the program are reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear Security Administration. The Government has granted for itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare derivative works, distribute copies to the public, perform publicly and display publicly, and to permit others to do so. +This program was produced under +U.S. Government contract 89233218CNA000001 +for Los Alamos National Laboratory (LANL), +which is operated by Triad National Security, LLC for the +U.S. Department of Energy/National Nuclear Security Administration. +All rights in the program are reserved by Triad National Security, LLC, +and the +U.S. Department of Energy/National Nuclear Security Administration. +The Government has granted for itself and others acting on its behalf +a nonexclusive, paid-up, irrevocable worldwide license in this material +to reproduce, +prepare derivative works, +distribute copies to the public, +perform publicly and display publicly, +and to permit others to do so. From 7613f9ed48e98141183952c27fbcd39053c4891b Mon Sep 17 00:00:00 2001 From: Michael-T-McCann Date: Wed, 20 Oct 2021 11:51:03 -0600 Subject: [PATCH 2/5] Read through Important Classes --- INSTALL.rst | 6 ++-- docs/source/functional.rst | 64 +++++++++++++++++++++++++------------- docs/source/overview.rst | 4 +-- 3 files changed, 49 insertions(+), 25 deletions(-) diff --git a/INSTALL.rst b/INSTALL.rst index a89632685..8879bddda 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -32,8 +32,10 @@ Automatic Installation The recommended way to install SCICO and its dependencies is via `conda `_ using the scripts in ``misc/conda``: - - ``install_conda.sh``: install ``miniconda`` (needed if conda is not already installed on your system) - - ``conda_env.sh``: install a ``conda`` environment with all SCICO dependencies + - ``install_conda.sh``: install ``miniconda`` + (needed if conda is not already installed on your system). + - ``conda_env.sh``: install a ``conda`` environment + with all SCICO dependencies. For GPU installation, see ``conda_env.sh -h``. Manual Installation diff --git a/docs/source/functional.rst b/docs/source/functional.rst index 772c118c4..1a0d4ed0a 100644 --- a/docs/source/functional.rst +++ b/docs/source/functional.rst @@ -22,25 +22,31 @@ Functionals and Losses } - -A functional maps an :code:`array-like` to a scalar; abstractly, a functional is +A functional is a mapping from :math:`\mathbb{R}^n` or :math:`\mathbb{C}^n` to :math:`\mathbb{R}`. - -A functional ``f`` can have three core operations. +In SCICO, functionals are +primarily used to represent a cost to be minimized +and are represented by instances of the :class:`.Functional` class. +An instance of :class:`.Functional`, ``f``, may provide three core operations. * Evaluation - - ``f(x)`` returns the value of the functional evaluated at :math:`\mb{x}`. - - A functional that can be evaluated has the attribute ``f.has_eval == True``. + - ``f(x)`` returns the value of the functional + evaluated at the point ``x``. + - A functional that can be evaluated + has the attribute ``f.has_eval == True``. - Not all functionals can be evaluated: see `Plug-and-Play`_. - * Gradient - - ``f.grad(x)`` returns the gradient of the functional evaluated at :math:`\mb{x}`. - - Calculated using JAX reverse-mode automatic differentiation, exposed through :func:`scico.grad`. + - ``f.grad(x)`` returns the gradient of the functional evaluated at ``x``. + - Gradients are calculated using JAX reverse-mode automatic differentiation, + exposed through :func:`scico.grad`. - A functional that is smooth has the attribute ``f.is_smooth == True``. - - NOTE: The gradient of a functional ``f`` can be evaluated even if ``f.is_smooth == False``. All that is required is that the functional can be evaluated, ``f.has_eval == True``. However, the result may not be a valid gradient (or subgradient) for all inputs :math:`\mb{x}`. - + - NOTE: The gradient of a functional ``f`` can be evaluated even if ``f.is_smooth == False``. + All that is required is that the functional can be evaluated, ``f.has_eval == True``. + However, the result may not be a valid gradient (or subgradient) for all inputs. * Proximal operator - - The proximal operator of a functional :math:`f : \mathbb{R}^n \to \mathbb{R}` is the mapping + - ``f.prox(v, lam)`` returns the result of the scaled proximal operator + at ``v`` with scale ``lam``. + - The proximal operator of a functional :math:`f : \mathbb{R}^n \to \mathbb{R}` is the mapping :math:`\mathrm{prox}_f : \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}^n` defined as .. math:: @@ -51,14 +57,18 @@ A functional ``f`` can have three core operations. Plug-and-Play ------------- -* For the Plug-and-Play framework :cite:`sreehari-2016-plug`, we encapsulate denoisers/CNNs in a Functional object that **cannot be evaluated**. -* Only the proximal operator is exposed. +For the plug-and-play framework :cite:`sreehari-2016-plug`, +we encapsulate generic denoisers including CNNs +in :class:`.Functional` objects that **cannot be evaluated**. +The denoiser is applied via the the proximal operator. +For examples, see :ref:`example_notebooks`. + Proximal Calculus ----------------- -We support a limited subset of proximal calculus rules. +We support a limited subset of proximal calculus rules: Scaled Functionals @@ -75,9 +85,15 @@ determine the proximal method of ``c * f`` as &= \mathrm{prox}_{f} (v, c \lambda) \end{align} -Note that we have made no assumptions regarding homogeneity of ``f``; rather, only that the proximal method of ``f`` is given in the parameterized form :math:`\mathrm{prox}_{c f}`. +Note that we have made no assumptions regarding homogeneity of ``f``; +rather, only that the proximal method of ``f`` is given +in the parameterized form :math:`\mathrm{prox}_{c f}`. -In SCICO, multiplying a :class:`.Functional` by a scalar will return a :class:`.ScaledFunctional`. This :class:`.ScaledFunctional` retains the ``has_eval, is_smooth``, and ``has_prox`` attributes from the original :class:`.Functional`, but the proximal method is modified to accomodate the additional scalar. +In SCICO, multiplying a :class:`.Functional` by a scalar +will return a :class:`.ScaledFunctional`. +This :class:`.ScaledFunctional` retains the ``has_eval``, ``is_smooth``, and ``has_prox`` attributes +from the original :class:`.Functional`, +but the proximal method is modified to accomodate the additional scalar. Separable Functionals @@ -89,7 +105,9 @@ of functionals :math:`f_i : \mathbb{C}^{N_i} \to \mathbb{R}` with :math:`\sum_i .. math:: f(\mb{x}) = f(\mb{x}_1, \dots, \mb{x}_N) = f_1(\mb{x}_1) + \dots + f_N(\mb{x}_N) -The proximal operator of a separable :math:`f` can be written in terms of the proximal operators of the :math:`f_i` (see Theorem 6.6 of :cite:`beck-2017-first`): +The proximal operator of a separable :math:`f` can be written +in terms of the proximal operators of the :math:`f_i` +(see Theorem 6.6 of :cite:`beck-2017-first`): .. math:: \mathrm{prox}_f(\mb{x}, \lambda) @@ -106,10 +124,14 @@ Separable Functionals are implemented in the :class:`.SeparableFunctional` class Adding New Functionals ---------------------- +To add a new functional, +create a class which + +1. inherits from base :class:`.Functional`; +2. has ``has_eval``, ``is_smooth``, and ``has_prox`` flags; +3. has ``_eval`` and ``prox`` methods, as necessary. -1. Inherit from base functional -2. Set ``has_eval``, ``is_smooth``, and ``has_prox`` flags. -3. Add ``_eval`` and ``prox`` methods, as necessary. +For example, :: diff --git a/docs/source/overview.rst b/docs/source/overview.rst index 0b30f7875..856ae8982 100644 --- a/docs/source/overview.rst +++ b/docs/source/overview.rst @@ -50,7 +50,7 @@ SCICO's main advantages are: including projected gradients and the alternating direction method of multipliers (ADMM) with the flexibility of plug-and-play priors - including BM3D and DnCNN denoisers. + including BM3D :cite:`dabov-2008-image` and DnCNN :cite:`zhang-2017-dncnn` denoisers. If you use this library for published work, @@ -147,7 +147,7 @@ automatic differentiation, and automatic adjoint calculations. Moreover, beause JAX is a machine learning library, it enables implementation and training -of state of the art Plug-and-Play regularizers such as DnCNN. +of state of the art plug-and-play regularizers such as DnCNN. Other related projects that may be of interest include: From 766750b847fdb1a9ade408c0a2318554e6a835ac Mon Sep 17 00:00:00 2001 From: Michael-T-McCann Date: Wed, 20 Oct 2021 13:12:33 -0600 Subject: [PATCH 3/5] Handle TODOs in Operator page --- docs/source/functional.rst | 1 - docs/source/op_linop.rst | 211 ++++++++++++++++--------------------- scico/operator/__init__.py | 4 + 3 files changed, 92 insertions(+), 124 deletions(-) diff --git a/docs/source/functional.rst b/docs/source/functional.rst index 1a0d4ed0a..e5751a4fe 100644 --- a/docs/source/functional.rst +++ b/docs/source/functional.rst @@ -64,7 +64,6 @@ The denoiser is applied via the the proximal operator. For examples, see :ref:`example_notebooks`. - Proximal Calculus ----------------- diff --git a/docs/source/op_linop.rst b/docs/source/op_linop.rst index b40a1d550..54ad98419 100644 --- a/docs/source/op_linop.rst +++ b/docs/source/op_linop.rst @@ -2,31 +2,29 @@ Operators ========= .. todo:: - * link to PyLops, scipy abstract linear operators * Document :func:`._wrap_mul_div_scalar`, :func:`._wrap_sum` -Matrix-free representation of operators. +An operator is a map from :math:`\mathbb{R}^n` or :math:`\mathbb{C}^n` +to :math:`\mathbb{R}^m` or :math:`\mathbb{C}^m`. +In SCICO, operators are primarily used to represent imaging systems +and provide regularization. +SCICO operators are represented by instances of the :class:`.Operator` class. -* Operator: generic operator -* LinearOperator: a linear map -SCICO :class:`.LinearOperator` is a linear operator that is designed to work with :class:`.BlockArray`. +SCICO :class:`.Operator` objects extend the notion of "shape" and "size" from the usual NumPy ``ndarray`` class. +Each :class:`.Operator` object has an ``input_shape`` and ``output_shape``; these shapes can be either tuples or a tuple of tuples +(in the case of a :class:`.BlockArray`). +The ``matrix_shape`` attribute describes the shape of the :class:`.LinearOperator` if it were to act on vectorized, or flattened, inputs. - -Consider a two dimensional array :math:`\mb{x} \in \mathbb{R}^{n \times m}`. +For example, consider a two dimensional array :math:`\mb{x} \in \mathbb{R}^{n \times m}`. We compute the discrete differences of :math:`\mb{x}` in the horizontal and vertical directions, generating two new arrays: :math:`\mb{x}_h \in \mathbb{R}^{n \times (m-1)}` and :math:`\mb{x}_v \in \mathbb{R}^{(n-1) \times m}`. We represent this linear operator by :math:`\mb{A} : \mathbb{R}^{n \times m} \to \mathbb{R}^{n \times (m-1)} \otimes \mathbb{R}^{(n-1) \times m}`. - In SCICO, this linear operator will return a :class:`.BlockArray` with the horizontal and vertical differences -stored as blocks. Letting :math:`y = \mb{A} x`, we have ``y.shape = ((n, m-1), (n-1, m))``. - - -SCICO :class:`.LinearOperator` objects extend the notion of "shape" and "size" from the usual NumPy ``ndarray`` class. -Each :class:`.LinearOperator` object has an ``input_shape`` and ``output_shape``; these shapes can be either tuples or a tuple of tuples -(in the case of a :class:`.BlockArray`). For the finite difference operator above, +stored as blocks. Letting :math:`y = \mb{A} x`, we have ``y.shape = ((n, m-1), (n-1, m))`` +and :: @@ -38,16 +36,11 @@ Each :class:`.LinearOperator` object has an ``input_shape`` and ``output_shape`` A.matrix_shape = (n*(n-1)*m*(m-1), n*m) # (output_size, input_size) -The ``matrix_shape`` attribute describes the shape of the :class:`.LinearOperator` if it were to act on vectorized, or flattened, -inputs. - - - - -Using An Operator +Operator Calculus ----------------- -Call: ``A(x)`` - +SCICO supports a variety of operator calculus rules, +allowing new operators to be defined in terms of old ones. +The following table summarizes the available operations. +----------------+-----------------+ | Operation | Result | @@ -70,8 +63,8 @@ Call: ``A(x)`` Defining A New Operator ----------------------- - -Pass a callable to the Operator constructor: +To define a new operator, +pass a callable to the :class:`.Operator` constructor: :: @@ -79,9 +72,7 @@ Pass a callable to the Operator constructor: eval_fn = lambda x: 2 * x) -Or via subclassing: - -At a minimum, the ``_eval`` function must be overridden: +Or use subclassing: :: @@ -93,74 +84,45 @@ At a minimum, the ``_eval`` function must be overridden: >>> A = MyOp(input_shape=(32,)) - - -* If either ``output_shape`` or ``output_dtype`` are unspecified, they are determined by evaluating - the Operator on an input of appropriate shape and dtype. - +At a minimum, the ``_eval`` function must be overridden. +If either ``output_shape`` or ``output_dtype`` are unspecified, they are determined by evaluating +the operator on an input of appropriate shape and dtype. Linear Operators ================ -Specialization of :class:`.Operator` to linear operators. - +Linear operators are those for which + .. math:: + H(a \mb{x} + b \mb{y}) = a H(\mb{x}) + b H(\mb{y}). - -Operations using LinearOperator -------------------------------- - -``A`` and ``B`` are :class:`.LinearOperator` of the same shape, -``x`` is an array of appropriate shape, -``c`` is a scalar, and -``O`` is :class:`.Operator` - -+----------------+----------------------------+ -| Operation | Result | -+----------------+----------------------------+ -| ``(A+B)(x)`` | ``A(x) + B(x)`` | -+----------------+----------------------------+ -| ``(A-B)(x)`` | ``A(x) - B(x)`` | -+----------------+----------------------------+ -| ``(c * A)(x)`` | ``c * A(x)`` | -+----------------+----------------------------+ -| ``(A/c)(x)`` | ``A(x)/c`` | -+----------------+----------------------------+ -| ``(-A)(x)`` | ``-A(x)`` | -+----------------+----------------------------+ -| ``(A@B)(x)`` | ``A@B@x`` | -+----------------+----------------------------+ -| ``A @ B`` | ``ComposedLinearOperator`` | -+----------------+----------------------------+ -| ``A @ O`` | ``Operator`` | -+----------------+----------------------------+ -| ``O(A)`` | ``Operator`` | -+----------------+----------------------------+ +SCICO represents linear operators as instances of the class :class:`.LinearOperator`. +While finite-dimensional linear operators +can always be associated with a matrix, +it is often useful to represent them in a matrix-free manner. +Most of SCICO's linear operators are implemented matrix-free. Using A LinearOperator ---------------------- -Evaluating a LinearOperator - - -We implement two ways to evaluate the LinearOperator. The first is using standard +We implement two ways to evaluate a :class:`.LinearOperator`. The first is using standard callable syntax: ``A(x)``. The second mimics the NumPy matrix multiplication syntax: ``A @ x``. Both methods perform shape and type checks to validate the -input before ultimately either calling `A._eval` or generating a new LinearOperator. +input before ultimately either calling `A._eval` or generating a new :class:`.LinearOperator`. -For LinearOperators that map real-valued inputs to real-valued outputs, there are two ways to apply the adjoint: +For linear operators that map real-valued inputs to real-valued outputs, there are two ways to apply the adjoint: ``A.adj(y)`` and ``A.T @ y``. -For complex-valued LinearOperators, there are three ways to apply the adjoint ``A.adj(y)``, ``A.H @ y``, and ``A.conj().T @ y``. +For complex-valued linear operators, there are three ways to apply the adjoint ``A.adj(y)``, ``A.H @ y``, and ``A.conj().T @ y``. Note that in this case, ``A.T`` returns the non-conjugated transpose of the LinearOperator. -While the cost of evaluating the LinearOperator is virtually identical for ``A(x)`` and ``A @ x``, +While the cost of evaluating the linear operator is virtually identical for ``A(x)`` and ``A @ x``, the ``A.H`` and ``A.conj().T`` methods are somewhat slower; especially the latter. This is because two -intermediate LinearOperators must be created before the function is evaluated. Evaluating ``A.conj().T @ y`` +intermediate linear operators must be created before the function is evaluated. Evaluating ``A.conj().T @ y`` is equivalent to: :: @@ -191,18 +153,46 @@ For instance: The public methods perform shape and type checking to validate the input before either calling the corresponding private method or returning a composite LinearOperator. -Jit Options ------------ -.. todo:: +Linear Operator Calculus +------------------------ +SCICO supports several linear operator calculus rules. +Given +``A`` and ``B`` of class :class:`.LinearOperator` and of appropriate shape, +``x`` an array of appropriate shape, +``c`` a scalar, and +``O`` an :class:`.Operator`, +we have + ++----------------+----------------------------+ +| Operation | Result | ++----------------+----------------------------+ +| ``(A+B)(x)`` | ``A(x) + B(x)`` | ++----------------+----------------------------+ +| ``(A-B)(x)`` | ``A(x) - B(x)`` | ++----------------+----------------------------+ +| ``(c * A)(x)`` | ``c * A(x)`` | ++----------------+----------------------------+ +| ``(A/c)(x)`` | ``A(x)/c`` | ++----------------+----------------------------+ +| ``(-A)(x)`` | ``-A(x)`` | ++----------------+----------------------------+ +| ``(A@B)(x)`` | ``A@B@x`` | ++----------------+----------------------------+ +| ``A @ B`` | ``ComposedLinearOperator`` | ++----------------+----------------------------+ +| ``A @ O`` | ``Operator`` | ++----------------+----------------------------+ +| ``O(A)`` | ``Operator`` | ++----------------+----------------------------+ - details -Defining A New LinearOperator +Defining A New Linear Operator ----------------------------- -Pass a callable to the LinearOperator constructor +To define a new linear operator, +pass a callable to the :class:`.LinearOperator` constructor :: @@ -211,9 +201,7 @@ Pass a callable to the LinearOperator constructor ... eval_fn = lambda x: 2 * x) -Subclassing: - -At a minimum, the ``_eval`` method must be overridden: +Or, use subclassing: :: @@ -223,54 +211,31 @@ At a minimum, the ``_eval`` method must be overridden: >>> A = MyLinearOperator(input_shape=(32,)) - -* If the ``_adj`` method is not overriden, the adjoint is determined using :func:`scico.linear_adjoint`. -* If either ``output_shape`` or ``output_dtype`` are unspecified, they are determined by evaluating - the Operator on an input of appropriate shape and dtype. - +At a minimum, the ``_eval`` method must be overridden. +If the ``_adj`` method is not overriden, the adjoint is determined using :func:`scico.linear_adjoint`. +If either ``output_shape`` or ``output_dtype`` are unspecified, they are determined by evaluating +the Operator on an input of appropriate shape and dtype. 🔪 Sharp Edges 🔪 ------------------ -Strict types in adjoint +Strict Types in Adjoint *********************** -.. todo:: - - We silently promote real->complex types in forward application, but have strict type checking in the adjoint. - This is due to the strict type-safe nature of jax adjoints +SCICO silently promotes real types to complex types in forward application, +but enforces strict type checking in the adjoint. +This is due to the strict type-safe nature of jax adjoints. -LinearOperators from External Code +LinearOperators From External Code ********************************** -.. todo:: - - Fill this out! - -* Pain point: adjoint and defining VJP for gradient computations - For example, might want to compute :math:`\nabla_x \norm{y - A x}_2^2` where :math:`A` is not a pure jax function -* Discuss VJP framework -* Can't use ``jax.linear_transpose``; must use VJP framework to determine adjoint -* Complexities for complex functions - -`Vector-Jacobian Product `_ - .. math:: - \begin{aligned} - &f : \mathbb{R}^n \to \mathbb{R}^m \\ - &\partial f(x) : \mathbb{R}^n \times \mathbb{R}^m \\ - &v \in \mathbb{R}^m \\ - &\mathrm{vjp}_f (x, v) \to v \partial f(x) - \end{aligned} - - .. math:: - \begin{aligned} - &A \in \mathbb{R}^{m \times n} \\ - &f(x) = A x \\ - & \partial f(x) = A \\ - &\mathrm{vjp}_f (x, v) \to v \partial f = v A = A^T v - \end{aligned} - - -`Custom VJP rules `_ +External code may be wrapped as a subclass of :class:`.Operator` or :class:`.LinearOperator` +and used in SCICO optimization routines; +however this process can be complicated and error-prone. +As a starting point, +look at the source for :class:`.radon_svmbir.ParallelBeamProjector` or :class:`.radon_astra.ParallelBeamProjector` +and the JAX documentation for the +`vector-jacobian product `_ +and `ustom VJP rules `_. diff --git a/scico/operator/__init__.py b/scico/operator/__init__.py index c8142cbd6..85acff459 100644 --- a/scico/operator/__init__.py +++ b/scico/operator/__init__.py @@ -17,3 +17,7 @@ "Operator", "BiConvolve", ] + +# Imported items in __all__ appear to originate in top-level linop module +for name in __all__: + getattr(sys.modules[__name__], name).__module__ = __name__ From c43e82902b99e0cee32781398fead9710f433bdc Mon Sep 17 00:00:00 2001 From: Michael-T-McCann Date: Wed, 20 Oct 2021 13:53:34 -0600 Subject: [PATCH 4/5] Fill in Losses section --- docs/source/functional.rst | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/docs/source/functional.rst b/docs/source/functional.rst index e5751a4fe..6bf31f129 100644 --- a/docs/source/functional.rst +++ b/docs/source/functional.rst @@ -1,5 +1,5 @@ -Functionals and Losses -====================== +Functionals +=========== .. raw:: html @@ -150,6 +150,18 @@ For example, Losses ------ -.. todo:: +In SCICO, a loss is a special type of functional - Content missing here + .. math:: + f(\mb{x}) = a l( \mb{y}, A(\mb{x}) ) + +where :math:`a` is a scale parameter, +:math:`l` is a functional, +:math:`\mb{y}` is a set of measurements, +and :math:`A` is an operator. +SCICO uses the class :class:`.Loss` to represent losses. +Loss functionals commonly arrise in the context of solving +inverse problems in scientific imaging, +where they are used to represent the mismatch +between predicted measurements :math:`A(\mb{x})` +and actual ones :math:`\mb{y}`. From aca86ec567e9c2a4c2c7096ff44d9476a9586e35 Mon Sep 17 00:00:00 2001 From: Michael-T-McCann Date: Wed, 20 Oct 2021 15:45:16 -0600 Subject: [PATCH 5/5] Minor corrections, reformatting --- docs/source/contributing.rst | 16 ++++---- docs/source/notes.rst | 22 +++++----- docs/source/style.rst | 79 +++++++++++++++++------------------- 3 files changed, 56 insertions(+), 61 deletions(-) diff --git a/docs/source/contributing.rst b/docs/source/contributing.rst index d84bc0120..974d13851 100644 --- a/docs/source/contributing.rst +++ b/docs/source/contributing.rst @@ -41,7 +41,7 @@ Installing a Development Version `the jax example `_) -1. Create a conda environment using Python >= 3.8. +1. Create a conda environment using Python >= 3.8: :: @@ -76,7 +76,7 @@ Installing a Development Version pip install -r examples/examples_requirements.txt # Installs example requirements pip install -e . # Installs SCICO from the current directory in editable mode. -6. Set up ``black`` and ``isort`` pre-commit hooks +6. Set up ``black`` and ``isort`` pre-commit hooks: :: @@ -150,7 +150,7 @@ NOTE: If you have added or modified an example script, see `Adding Usage Exampl git push --set-upstream origin name-of-change -9. Create a new pull request to the ``main`` branch; see `the GitHub instructions `_ +9. Create a new pull request to the ``main`` branch; see `the GitHub instructions `_. 10. Delete the branch after it has been merged. @@ -163,7 +163,7 @@ existing examples to ensure that the mechanism for automatically generating corresponding Jupyter notebooks functions correctly. In particular: -1. The initial lines of the script should consist of a comment block, followed by a blank line, followed by a multiline string with an RST heading on the first line, e.g. +1. The initial lines of the script should consist of a comment block, followed by a blank line, followed by a multiline string with an RST heading on the first line, e.g., :: @@ -218,7 +218,7 @@ and ``scico-data`` repositories must be updated and kept in sync. 3. Convert your new example to a Jupyter notebook by changing directory to the ``scico/examples`` directory and following the instructions in ``scico/examples/README.rst``. -4. Change directory to the ``data`` directory and add/commit the new Jupyter Notebook +4. Change directory to the ``data`` directory and add/commit the new Jupyter Notebook: :: @@ -254,7 +254,7 @@ scico-data repositories must be updated and kept in sync. 1. Add the ``new_data.npz`` file to the ``scico/data`` directory. -2. Navigate to the ``data`` directory and add/commit the new data file +2. Navigate to the ``data`` directory and add/commit the new data file: :: @@ -288,7 +288,7 @@ Running Tests ------------- -To be able to run the tests, install `pytest` and, optionally, `pytest-runner` +To be able to run the tests, install `pytest` and, optionally, `pytest-runner`: :: @@ -342,7 +342,7 @@ Test Coverage Test coverage is a measure of the fraction of the package code that is exercised by the tests. While this should not be the primary criterion in designing tests, it is a useful tool for finding obvious areas of omission. -To be able to check test coverage, install `coverage` +To be able to check test coverage, install `coverage`: :: diff --git a/docs/source/notes.rst b/docs/source/notes.rst index 2fd66b2f8..2b2b9fdba 100644 --- a/docs/source/notes.rst +++ b/docs/source/notes.rst @@ -81,7 +81,7 @@ Double Precision By default, JAX enforces single-precision numbers. Double precision can be enabled in one of two ways: 1. Setting the environment variable ``JAX_ENABLE_X64=TRUE`` before launching python. -2. Manually set the ``jax_enable_x64`` flag **at program startup**; that is, **before** importing SCICO. +2. Manually setting the ``jax_enable_x64`` flag **at program startup**; that is, **before** importing SCICO. :: @@ -90,7 +90,7 @@ By default, JAX enforces single-precision numbers. Double precision can be enabl import scico # continue as usual -For more information, `see the JAX notes on double precision `_ +For more information, see the `JAX notes on double precision `_. Random Number Generation @@ -146,11 +146,11 @@ The function :func:`scico.grad` returns the expected gradient, that is, the conj JAX gradient. For further discussion, see this `JAX issue `_. -As a concrete example, consider the function :math:`f(x) = \frac{1}{2}\norm{A -x}_2^2` where :math:`A` is a complex matrix. The gradient of :math:`f` is -usually given :math:`(\nabla f)(x) = A^H A x`, where :math:`A^H` is the -conjugate transpose of :math:`A`. Applying ``jax.grad`` to :math:`f` will yield -:math:`(A^H A x)^*`, where :math:`*` denotes complex conjugation. +As a concrete example, consider the function :math:`f(x) = \frac{1}{2}\norm{\mb{A} +\mb{x}}_2^2` where :math:`\mb{A}` is a complex matrix. The gradient of :math:`f` is +usually given :math:`(\nabla f)(\mb{x}) = \mb{A}^H \mb{A} \mb{x}`, where :math:`\mb{A}^H` is the +conjugate transpose of :math:`\mb{A}`. Applying ``jax.grad`` to :math:`f` will yield +:math:`(\mb{A}^H \mb{A} \mb{x})^*`, where :math:`*` denotes complex conjugation. The following code demonstrates the use of ``jax.grad`` and :func:`scico.grad`: @@ -173,11 +173,11 @@ The following code demonstrates the use of ``jax.grad`` and :func:`scico.grad`: Non-differentiable Functionals and scico.grad --------------------------------------------- -* :func:`scico.grad` can be applied to any function, but has undefined behavior for - non-differentiable functions. -* For non-differerentiable functions, :func:`scico.grad` may or may not return a valid subgradient. As an example, ``scico.grad(snp.abs)(0.) = 0``, which is a valid subgradient. However, ``scico.grad(snp.linalg.norm)([0., 0.]) = [nan, nan]``, which is not a valid subgradient of this function. -* Differentiable functions that are written as the composition of a differentiable and non-differentiable function should be avoided. As an example, :math:`f(x) = \norm{x}_2^2` can be implemented in as ``f = lambda x: snp.linalg.norm(x)**2``. This involves first calculating the non-squared :math:`\ell_2` norm, then squaring it. The un-squared :math:`\ell_2` norm is not differentiable at zero. +:func:`scico.grad` can be applied to any function, but has undefined behavior for +non-differentiable functions. +For non-differerentiable functions, :func:`scico.grad` may or may not return a valid subgradient. As an example, ``scico.grad(snp.abs)(0.) = 0``, which is a valid subgradient. However, ``scico.grad(snp.linalg.norm)([0., 0.]) = [nan, nan]``. +Differentiable functions that are written as the composition of a differentiable and non-differentiable function should be avoided. As an example, :math:`f(x) = \norm{x}_2^2` can be implemented in as ``f = lambda x: snp.linalg.norm(x)**2``. This involves first calculating the non-squared :math:`\ell_2` norm, then squaring it. The un-squared :math:`\ell_2` norm is not differentiable at zero. When evaluating the gradient of ``f`` at 0, :func:`scico.grad` returns ``nan``: :: diff --git a/docs/source/style.rst b/docs/source/style.rst index 20cfe20a2..7eee6af14 100644 --- a/docs/source/style.rst +++ b/docs/source/style.rst @@ -90,15 +90,15 @@ Things to avoid: - Single character names except for the following special cases: - - Counters or iterators (``i``, ``j``) - - `e` as an exception identifier (``Exception e``) - - `f` as a file in ``with`` statements - - Mathematical notation in which a reference to the paper or algorithm with said notation is preferred if not clear from the intended purpose. + - counters or iterators (``i``, ``j``); + - `e` as an exception identifier (``Exception e``); + - `f` as a file in ``with`` statements; + - mathematical notation in which a reference to the paper or algorithm with said notation is preferred if not clear from the intended purpose. - Trailing underscores unless the component is meant to be protected or private: - - Protected: Use a single underscore, ``_``, for protected access - - Pseudo-private: Use double underscores, ``_``, for pseudo-private access via name mangling. + - protected: Use a single underscore, ``_``, for protected access; and + - pseudo-private: Use double underscores, ``_``, for pseudo-private access via name mangling. | @@ -133,18 +133,14 @@ Variables Apart from naming conventions there are a few extra documentation and coding practices that can be applied to variables such as: -- Typing Variables: - - - Use a ``: type`` before the function value is assigned - - | Example: +- One may type a variables by using a ``: type`` before the function value is assigned, e.g., .. code-block:: python a : Foo = SomeDecoratedFunction() - - Avoid global variables. - - A function can refer to variables defined in enclosing functions but cannot assign to them. +- Avoid global variables. +- A function can refer to variables defined in enclosing functions but cannot assign to them. | @@ -153,10 +149,11 @@ Parameters There are three important stlyle components for parameters: -- We use type annotations meaning we specify the types of the inputs and outputs of any method. - - From the ``typing`` module we can use more types such as ``Optional``, ``Union``, and ``Any``. +1. Typing - | Example: + We use type annotations meaning we specify the types of the inputs and outputs of any method. + From the ``typing`` module we can use more types such as ``Optional``, ``Union``, and ``Any``. + For example, .. code-block:: python @@ -164,13 +161,13 @@ There are three important stlyle components for parameters: """Takes an input of type string and returns a value of type string""" ... -- Default Values: - - - Parameter should include ``parameter_name = value`` where value is the default for that particular parameter. - - If the parameter has a type then the format is ``parameter_name: Type = value`` - - Additionally when documenting parameters, a parameter can only assume one of a fixed set of values, those values can be listed in braces, with the default appearing first. +2. Default Values - | Example: + Parameters should include ``parameter_name = value`` where value is the default for that particular parameter. + If the parameter has a type then the format is ``parameter_name: Type = value``. + When documenting parameters, if a parameter can only assume one of a fixed set of values, + those values can be listed in braces, with the default appearing first. + For example, .. code-block:: python @@ -179,19 +176,21 @@ There are three important stlyle components for parameters: Description of `letters`. """ -- NoneType: +3. NoneType - - In Python a ``NoneType`` is a "first-class" type. - - ``None`` is the most commonly used alias. \* If any of the parameters/arguments can be ``None`` then it has to be declared. ``Optional[T]`` is preferred over ``Union[T, None]``. - - | Example: + In Python, ``NoneType`` is a first-class type, meaning the type itself + can be passed into and returned from functions. + ``None`` is the most commonly used alias for ``NoneType``. + If any of a function's parameters can be ``None`` then it has to be declared. + ``Optional[T]`` is preferred over ``Union[T, None]``. + For example, .. code-block:: python def foo(a: Optional[str], b: Optional[Union[str, int]]) -> str: ... - - For documentation purposes, ``NoneType`` or ``None`` should be written with double backticks + For documentation purposes, ``NoneType`` or ``None`` should be written with double backticks. | @@ -207,15 +206,14 @@ Typing The following are docstring-specific usages that must be explained before going into the creation of said docstrings: - Always enclose variables in single backticks. -- For the parameter types, be as precise as possible, no need for backticks. +- For the parameter types, be as precise as possible, do not use backticks. Modules ~~~~~~~ Files must start with a docstring that describes the functionality of the module. - -Example: +For example, .. code-block:: python @@ -233,14 +231,12 @@ Example: Functions ~~~~~~~~~ -| The word *function* encompasses functions, methods, or generators in this section. - +The word *function* encompasses functions, methods, or generators in this section. The docstring should give enough information to make calls to the function without needing to read the functions code. Functions should contain docstrings unless: - -- not externally visible (the function name is prefaced with an underscore) -- very short +- not externally visible (the function name is prefaced with an underscore) or +- very short. The docstring should be imperative-style ``"""Fetch rows from a Table"""`` instead of the descriptive-style ``"""Fetches rows from a Table"""``. If the method overrides a method from a base class then it may use a simple docstring referencing that base class such as ``"""See base class"""``, unless the behavior is different from the overridden method or there are extra details that need to be documented. @@ -353,9 +349,8 @@ The following are sections that can be added to functions, modules, classes, or - Notes - - Provides additional information about the code. May include mathematical equations in LaTeX format: - - | Example: + - Provides additional information about the code. May include mathematical equations in LaTeX format. + For example, .. code-block:: python @@ -381,8 +376,7 @@ The following are sections that can be added to functions, modules, classes, or - Uses the doctest format and is meant to showcase usage. - If there are multiple examples include blank lines before and after each example. - - | Example: + For example, .. code-block:: python @@ -416,7 +410,8 @@ The following are sections that can be added to functions, modules, classes, or Comments ~~~~~~~~ -There are two types of comments: *block* and *inline*. A good rule of thumb to follow for when to include a comment in your code is: if you have to explain it or is too hard to figure out at first glance, then comment it. An example of this is complicated operations which most likely require a block of comments beforehand. +There are two types of comments: *block* and *inline*. A good rule of thumb to follow for when to include a comment in your code is *if you have to explain it or is too hard to figure out at first glance, then comment it*. +An example of this is complicated operations which most likely require a block of comments beforehand. .. code-block:: Python