Clean up example scripts

README.md

@@ -8,8 +8,9 @@
 [![Test coverage](https://codecov.io/gh/lanl/scico/branch/main/graph/badge.svg?token=wQimmjnzFf)](https://codecov.io/gh/lanl/scico)
 [![CodeFactor](https://www.codefactor.io/repository/github/lanl/scico/badge/main)](https://www.codefactor.io/repository/github/lanl/scico/overview/main)\
 [![PyPI package version](https://badge.fury.io/py/scico.svg)](https://badge.fury.io/py/scico)
-[![PyPI download statistics](https://static.pepy.tech/personalized-badge/scico?period=month&left_color=grey&right_color=brightgreen)](https://pepy.tech/project/scico)
-[![Conda Forge Release](https://img.shields.io/conda/vn/conda-forge/scico.svg)](https://anaconda.org/conda-forge/scico)\
+[![PyPI download statistics](https://static.pepy.tech/personalized-badge/scico?period=total&left_color=grey&right_color=brightgreen&left_text=downloads)](https://pepy.tech/project/scico)
+[![Conda Forge Release](https://img.shields.io/conda/vn/conda-forge/scico.svg)](https://anaconda.org/conda-forge/scico)
+[![Conda Forge Downloads](https://img.shields.io/conda/dn/conda-forge/scico.svg)](https://anaconda.org/conda-forge/scico)\
 [![View notebooks at nbviewer](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg)](https://nbviewer.jupyter.org/github/lanl/scico-data/tree/main/notebooks/index.ipynb)
 [![Run notebooks on binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/lanl/scico-data/binder?labpath=notebooks%2Findex.ipynb)
 [![Run notebooks on google colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/lanl/scico-data/blob/colab/notebooks/index.ipynb)

examples/scripts/ct_abel_tv_admm.py

@@ -8,16 +8,16 @@
 TV-Regularized Abel Inversion
 =============================
 
-This example demonstrates a TV-regularized Abel inversion by solving the
-problem
+This example demonstrates a total variation (TV) regularized Abel
+inversion by solving the problem
 
   $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x}
   \|_2^2 + \lambda \| C \mathbf{x} \|_1 \;,$$
 
 where $A$ is the Abel projector (with an implementation based on a
 projector from PyAbel :cite:`pyabel-2022`), $\mathbf{y}$ is the measured
 data, $C$ is a 2D finite difference operator, and $\mathbf{x}$ is the
-desired image.
+solution.
 """
 
 

examples/scripts/ct_abel_tv_admm_tune.py

@@ -23,7 +23,7 @@
 because `JAX_PLATFORMS` was intended to replace `JAX_PLATFORM_NAME` but
 this change has yet to be correctly implemented
 (see [google/jax#6805](https://github.com/google/jax/issues/6805) and
-[google/jax#10272](https://github.com/google/jax/pull/10272).
+[google/jax#10272](https://github.com/google/jax/pull/10272)).
 """
 
 # isort: off

examples/scripts/ct_astra_3d_tv_admm.py

@@ -17,7 +17,7 @@
 
 where $C$ is the X-ray transform (the CT forward projection operator),
 $\mathbf{y}$ is the sinogram, $D$ is a 3D finite difference operator,
-and $\mathbf{x}$ is the desired image.
+and $\mathbf{x}$ is the reconstructed image.
 
 In this example the problem is solved via ADMM, while proximal
 ADMM is used in a [companion example](ct_astra_3d_tv_padmm.rst).

examples/scripts/ct_astra_3d_tv_padmm.py

@@ -17,7 +17,7 @@
 
 where $C$ is the X-ray transform (the CT forward projection operator),
 $\mathbf{y}$ is the sinogram, $D$ is a 3D finite difference operator,
-and $\mathbf{x}$ is the desired image.
+and $\mathbf{x}$ is the reconstructed image.
 
 In this example the problem is solved via proximal ADMM, while standard
 ADMM is used in a [companion example](ct_astra_3d_tv_admm.rst).

examples/scripts/ct_astra_noreg_pcg.py

@@ -58,8 +58,8 @@
 
 r"""
 Invert in the Fourier domain to form a preconditioner $\mathbf{M}
-\approx (\mathbf{A}^T \mathbf{A})^{-1}$. See
-:cite:`clinthorne-1993-preconditioning` Section V.A. for more details.
+\approx (\mathbf{A}^T \mathbf{A})^{-1}$ (see
+:cite:`clinthorne-1993-preconditioning` Section V.A. for more details).
 """
 # γ limits the gain of the preconditioner; higher gives a weaker filter.
 γ = 1e-2
@@ -76,7 +76,9 @@
 r"""
 Check that $\mathbf{M}$ does approximately invert $\mathbf{A}^T \mathbf{A}$.
 """
-plot_args = dict(norm=plot.matplotlib.colors.Normalize(vmin=0, vmax=1.5))
+plot_args = dict(
+    norm=plot.matplotlib.colors.Normalize(vmin=0, vmax=1.5), cmap=plot.matplotlib.cm.Blues_r
+)
 
 fig, axes = plot.subplots(nrows=1, ncols=3, figsize=(12, 4.5))
 plot.imview(x_gt, title="Ground truth, $x_{gt}$", fig=fig, ax=axes[0], **plot_args)
@@ -96,8 +98,8 @@
     axes[2].get_images()[0],
     ax=axes,
     location="right",
-    shrink=1.0,
-    pad=0.05,
+    shrink=0.82,
+    pad=0.02,
     label="Arbitrary Units",
 )
 fig.show()
@@ -134,13 +136,13 @@
 f_cg = loss.SquaredL2Loss(y=A.T @ y, A=A.T @ A)
 f_data = loss.SquaredL2Loss(y=y, A=A)
 print(
-    f"{'Method':10s}{'Iterations':>15s}{'Time (s)':>15s}{'||ATAx - ATy||':>15s}{'||Ax - y||':>15s}"
+    f"{'Method':8s}{'Iterations':>11s}{'Time (s)':>12s}{'||ATAx - ATy||':>17s}{'||Ax - y||':>15s}"
 )
 print(
-    f"{'CG':10s}{info_cg['num_iter']:>15d}{time_cg:>15.2f}{f_cg(x_cg):>15.2e}{f_data(x_cg):>15.2e}"
+    f"{'CG':8s}{info_cg['num_iter']:>11d}{time_cg:>12.2f}{f_cg(x_cg):>17.2e}{f_data(x_cg):>15.2e}"
 )
 print(
-    f"{'PCG':10s}{info_pcg['num_iter']:>15d}{time_pcg:>15.2f}{f_cg(x_pcg):>15.2e}"
+    f"{'PCG':8s}{info_pcg['num_iter']:>11d}{time_pcg:>12.2f}{f_cg(x_pcg):>17.2e}"
     f"{f_data(x_pcg):>15.2e}"
 )
 

examples/scripts/ct_astra_tv_admm.py

@@ -16,8 +16,8 @@
 
 where $A$ is the X-ray transform (the CT forward projection operator),
 $\mathbf{y}$ is the sinogram, $C$ is a 2D finite difference operator, and
-$\mathbf{x}$ is the desired image. This example uses the CT projector
-provided by the astra package, while the companion
+$\mathbf{x}$ is the reconstructed image. This example uses the CT
+projector provided by the astra package, while the companion
 [example script](ct_tv_admm.rst) uses the projector integrated into
 scico.
 """

examples/scripts/ct_astra_weighted_tv_admm.py

@@ -8,8 +8,8 @@
 TV-Regularized Low-Dose CT Reconstruction
 =========================================
 
-This example demonstrates solution of a low-dose CT reconstruction problem
-with isotropic total variation (TV) regularization
+This example demonstrates solution of a low-dose CT reconstruction
+problem with isotropic total variation (TV) regularization
 
   $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x}
   \|_W^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$
@@ -18,7 +18,7 @@
 $\mathbf{y}$ is the sinogram, the norm weighting $W$ is chosen so that
 the weighted norm is an approximation to the Poisson negative log
 likelihood :cite:`sauer-1993-local`, $C$ is a 2D finite difference
-operator, and $\mathbf{x}$ is the desired image.
+operator, and $\mathbf{x}$ is the reconstructed image.
 """
 
 import numpy as np
@@ -56,11 +56,11 @@
 r"""
 Add Poisson noise to projections according to
 
-$$\mathrm{counts} \sim \mathrm{Poi}\left(I_0 exp\left\{- \alpha A
-\mathbf{x} \right\}\right)$$
+$$\mathrm{counts} \sim \mathrm{Poi}\left(I_0 \exp (- \alpha A
+\mathbf{x} ) \right)$$
 
 $$\mathbf{y} = - \frac{1}{\alpha} \log\left(\mathrm{counts} /
-I_0\right).$$
+I_0\right) \;.$$
 
 We use the NumPy random functionality so we can generate using 64-bit
 numbers.
@@ -124,7 +124,7 @@ def postprocess(x):
 
 where
 
-  $$W = \mathrm{diag}\left\{ \mathrm{counts} / I_0 \right\} \;.$$
+  $$W = \mathrm{diag}( \mathrm{counts} / I_0 ) \;.$$
 
 The data fidelity term in this formulation follows
 :cite:`sauer-1993-local` (9) except for the scaling by $I_0$, which we

examples/scripts/ct_fan_svmbir_ppp_bm3d_admm_prox.py

@@ -46,7 +46,9 @@
 density = 0.025  # attenuation density of the image
 np.random.seed(1234)
 pad_len = 5
-x_gt = discrete_phantom(Foam(size_range=[0.05, 0.02], gap=0.02, porosity=0.3), size=N - 2 * pad_len)
+x_gt = discrete_phantom(
+    Foam(size_range=[0.075, 0.005], gap=2e-3, porosity=1.0), size=N - 2 * pad_len
+)
 x_gt = x_gt / np.max(x_gt) * density
 x_gt = np.pad(x_gt, pad_len)
 x_gt[x_gt < 0] = 0
@@ -267,7 +269,7 @@ def add_poisson_noise(sino, max_intensity):
     fig=fig,
     ax=ax[0],
 )
-ax[0].set_ylim([5e-3, 5e0])
+ax[0].set_ylim([1e-1, 1e1])
 ax[0].xaxis.set_major_locator(MaxNLocator(integer=True))
 plot.plot(
     snp.vstack((hist_extloss_fan.Prml_Rsdl, hist_extloss_fan.Dual_Rsdl)).T,
@@ -278,7 +280,7 @@ def add_poisson_noise(sino, max_intensity):
     fig=fig,
     ax=ax[1],
 )
-ax[1].set_ylim([5e-3, 5e0])
+ax[1].set_ylim([1e-1, 1e1])
 ax[1].xaxis.set_major_locator(MaxNLocator(integer=True))
 fig.show()
 

examples/scripts/ct_multi_tv_admm.py

@@ -16,9 +16,9 @@
 
 where $A$ is the X-ray transform (the CT forward projection operator),
 $\mathbf{y}$ is the sinogram, $C$ is a 2D finite difference operator, and
-$\mathbf{x}$ is the desired image. The solution is computed and compared
-for all three 2D CT projectors available in scico, using a sinogram
-computed with the astra projector.
+$\mathbf{x}$ is the reconstructed image. The solution is computed and
+compared for all three 2D CT projectors available in scico, using a
+sinogram computed with the astra projector.
 """
 
 import numpy as np

examples/scripts/ct_svmbir_ppp_bm3d_admm_cg.py

@@ -41,7 +41,7 @@
 N = 256  # image size
 density = 0.025  # attenuation density of the image
 np.random.seed(1234)
-x_gt = discrete_phantom(Foam(size_range=[0.05, 0.02], gap=0.02, porosity=0.3), size=N - 10)
+x_gt = discrete_phantom(Foam(size_range=[0.075, 0.005], gap=2e-3, porosity=1.0), size=N - 10)
 x_gt = x_gt / np.max(x_gt) * density
 x_gt = np.pad(x_gt, 5)
 x_gt[x_gt < 0] = 0
@@ -105,7 +105,7 @@
     x0=x0,
     maxiter=20,
     subproblem_solver=LinearSubproblemSolver(cg_kwargs={"tol": 1e-4, "maxiter": 100}),
-    itstat_options={"display": True, "period": 1},
+    itstat_options={"display": True, "period": 5},
 )
 
 

examples/scripts/ct_svmbir_ppp_bm3d_admm_prox.py

@@ -55,7 +55,7 @@
 N = 256  # image size
 density = 0.025  # attenuation density of the image
 np.random.seed(1234)
-x_gt = discrete_phantom(Foam(size_range=[0.05, 0.02], gap=0.02, porosity=0.3), size=N - 10)
+x_gt = discrete_phantom(Foam(size_range=[0.075, 0.005], gap=2e-3, porosity=1.0), size=N - 10)
 x_gt = x_gt / np.max(x_gt) * density
 x_gt = np.pad(x_gt, 5)
 x_gt[x_gt < 0] = 0
@@ -219,7 +219,7 @@
     fig=fig,
     ax=ax[0],
 )
-ax[0].set_ylim([5e-3, 5e0])
+ax[0].set_ylim([1e-1, 5e0])
 ax[0].xaxis.set_major_locator(MaxNLocator(integer=True))
 plot.plot(
     snp.vstack((hist_extloss.Prml_Rsdl, hist_extloss.Dual_Rsdl)).T,
@@ -230,7 +230,7 @@
     fig=fig,
     ax=ax[1],
 )
-ax[1].set_ylim([5e-3, 5e0])
+ax[1].set_ylim([1e-1, 5e0])
 ax[1].xaxis.set_major_locator(MaxNLocator(integer=True))
 fig.show()
 

examples/scripts/ct_svmbir_tv_multi.py

@@ -16,8 +16,8 @@
 
 where $A$ is the X-ray transform (implemented using the SVMBIR
 :cite:`svmbir-2020` tomographic projection), $\mathbf{y}$ is the sinogram,
-$C$ is a 2D finite difference operator, and $\mathbf{x}$ is the desired
-image.
+$C$ is a 2D finite difference operator, and $\mathbf{x}$ is the
+reconstructed image.
 """
 
 import numpy as np
@@ -40,7 +40,7 @@
 N = 256  # image size
 density = 0.025  # attenuation density of the image
 np.random.seed(1234)
-x_gt = discrete_phantom(Foam(size_range=[0.05, 0.02], gap=0.02, porosity=0.3), size=N - 10)
+x_gt = discrete_phantom(Foam(size_range=[0.075, 0.005], gap=2e-3, porosity=1.0), size=N - 10)
 x_gt = x_gt / np.max(x_gt) * density
 x_gt = np.pad(x_gt, 5)
 x_gt[x_gt < 0] = 0

examples/scripts/ct_tv_admm.py

@@ -16,8 +16,8 @@
 
 where $A$ is the X-ray transform (the CT forward projection operator),
 $\mathbf{y}$ is the sinogram, $C$ is a 2D finite difference operator, and
-$\mathbf{x}$ is the desired image. This example uses the CT projector
-integrated into scico, while the companion
+$\mathbf{x}$ is the reconstructed image. This example uses the CT
+projector integrated into scico, while the companion
 [example script](ct_astra_tv_admm.rst) uses the projector provided by
 the astra package.
 """

examples/scripts/deconv_datagen_foam1.py

@@ -9,8 +9,8 @@
 ===============================================
 
 This example demonstrates how to generate blurred image data for
-training neural network models for deconvolution (deblurring). Foam
-phantoms from xdesign are used to generate the clean images.
+training neural network models for deconvolution (deblurring), using foam
+phantoms generated by `xdesign`.
 """
 
 # isort: off

examples/scripts/deconv_microscopy_allchn_tv_admm.py

@@ -25,7 +25,7 @@
 where $M$ is a mask operator, $A$ is circular convolution,
 $\mathbf{y}$ is the blurred image, $C$ is a convolutional gradient
 operator, $\iota_{\mathrm{NN}}$ is the indicator function of the
-non-negativity constraint, and $\mathbf{x}$ is the desired image.
+non-negativity constraint, and $\mathbf{x}$ is the deconvolved image.
 """
 
 # isort: off
@@ -48,8 +48,8 @@
 example on a GPU it may be necessary to set environment variables
 `XLA_PYTHON_CLIENT_ALLOCATOR=platform` and
 `XLA_PYTHON_CLIENT_PREALLOCATE=false`. If your GPU does not have enough
-memory, you can try setting the environment variable
-`JAX_PLATFORM_NAME=cpu` to run on CPU.
+memory, try setting the environment variable `JAX_PLATFORM_NAME=cpu` to
+run on CPU.
 """
 downsampling_rate = 2
 

examples/scripts/deconv_microscopy_tv_admm.py

@@ -25,7 +25,7 @@
 where $M$ is a mask operator, $A$ is circular convolution,
 $\mathbf{y}$ is the blurred image, $C$ is a convolutional gradient
 operator, $\iota_{\mathrm{NN}}$ is the indicator function of the
-non-negativity constraint, and $\mathbf{x}$ is the desired image.
+non-negativity constraint, and $\mathbf{x}$ is the deconvolved image.
 """
 
 import scico.numpy as snp
@@ -40,8 +40,8 @@
 example on a GPU it may be necessary to set environment variables
 `XLA_PYTHON_CLIENT_ALLOCATOR=platform` and
 `XLA_PYTHON_CLIENT_PREALLOCATE=false`. If your GPU does not have enough
-memory, you can try setting the environment variable
-`JAX_PLATFORM_NAME=cpu` to run on CPU.
+memory, try setting the environment variable `JAX_PLATFORM_NAME=cpu` to
+run on CPU.
 """
 channel = 0
 downsampling_rate = 2

examples/scripts/deconv_tv_padmm.py

@@ -15,7 +15,7 @@
   \|_2^2 + \lambda \| D \mathbf{x} \|_{2,1} \;,$$
 
 where $C$ is a convolution operator, $\mathbf{y}$ is the blurred image,
-$D$ is a 2D finite fifference operator, and $\mathbf{x}$ is the
+$D$ is a 2D finite difference operator, and $\mathbf{x}$ is the
 deconvolved image.
 
 In this example the problem is solved via proximal ADMM, while standard

examples/scripts/denoise_tv_admm.py

@@ -16,7 +16,7 @@
   \|_2^2 + \lambda R(\mathbf{x}) \;,$$
 
 where $R$ is either the isotropic or anisotropic TV regularizer.
-In SCICO, switching between these two regularizers is a one-line
+In SCICO, switching between these two regularizers involves a one-line
 change: replacing an
 [L1Norm](../_autosummary/scico.functional.rst#scico.functional.L1Norm)
 with a

examples/scripts/denoise_tv_apgm.py

@@ -121,8 +121,8 @@ def prox(self, v: Array, lam: float, **kwargs) -> Array:
 
 
 """
-Use RobustLineSearchStepSize object and set up AcceleratedPGM solver
-object. Run the solver.
+Set up `AcceleratedPGM` solver object using `RobustLineSearchStepSize`
+step size policy. Run the solver.
 """
 reg_weight_iso = 1.4e0
 f_iso = DualTVLoss(y=y, A=A, lmbda=reg_weight_iso)
@@ -165,9 +165,9 @@ def prox(self, v: Array, lam: float, **kwargs) -> Array:
 
 
 """
-Use RobustLineSearchStepSize object and set up AcceleratedPGM solver
-object. Weight was tuned to give the same data fidelity as the
-isotropic case. Run the solver.
+Set up `AcceleratedPGM` solver object using `RobustLineSearchStepSize`
+step size policy. (Weight was tuned to give the same data fidelity as the
+isotropic case.) Run the solver.
 """
 
 reg_weight_aniso = 1.2e0

examples/scripts/denoise_tv_multi.py

@@ -65,7 +65,7 @@
 the timing results. The precise cause of the remaining differences in
 time required to compute the first step of each algorithm is unknown,
 but it is worth noting that this difference becomes negligible when
-just-in-time compilation is disabled (e.g. via the JAX_DISABLE_JIT
+just-in-time compilation is disabled (e.g. via the `JAX_DISABLE_JIT`
 environment variable).
 """
 solver_admm = ADMM(