From 64d5882d2aa36cdd2572c904db969d4806d7234d Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 17:27:29 -0800
Subject: [PATCH 01/34] feat: implement PyTensorOperator

---
 pylops/__init__.py               |  1 +
 pylops/pytensoroperator.py       | 85 ++++++++++++++++++++++++++++++++
 pylops/utils/deps.py             | 19 +++++++
 pytests/test_pytensoroperator.py | 51 +++++++++++++++++++
 4 files changed, 156 insertions(+)
 create mode 100644 pylops/pytensoroperator.py
 create mode 100755 pytests/test_pytensoroperator.py

diff --git a/pylops/__init__.py b/pylops/__init__.py
index 7672fda4..b1ee04c8 100755
--- a/pylops/__init__.py
+++ b/pylops/__init__.py
@@ -48,6 +48,7 @@
 from .config import *
 from .linearoperator import *
 from .torchoperator import *
+from .pytensoroperator import *
 from .jaxoperator import *
 from .basicoperators import *
 from . import (
diff --git a/pylops/pytensoroperator.py b/pylops/pytensoroperator.py
new file mode 100644
index 00000000..cdd04c04
--- /dev/null
+++ b/pylops/pytensoroperator.py
@@ -0,0 +1,85 @@
+import pylops
+from pylops.utils import deps
+
+pytensor_message = deps.pytensor_import("the pytensor module")
+
+if pytensor_message is not None:
+
+    class PyTensorOperator:
+        """PyTensor Op which applies a PyLops Linear Operator, including gradient support.
+
+        This class "converts" a PyLops `LinearOperator` class into a PyTensor `Op`.
+        This applies the `LinearOperator` in "forward-mode" in `self.perform`, and applies
+        its adjoint when computing the vector-Jacobian product (`self.grad`), as that is
+        the analytically correct gradient for linear operators. This class should pass
+        `pytensor.gradient.verify_grad`.
+
+        Parameters
+        ----------
+        LOp : pylops.LinearOperator
+        """
+
+        def __init__(self, LOp: pylops.LinearOperator) -> None:
+            if not deps.pytensor_enabled:
+                raise NotImplementedError(pytensor_message)
+
+else:
+    import pytensor.tensor as pt
+    from pytensor.graph.basic import Apply
+    from pytensor.graph.op import Op
+
+    class _PyTensorOperatorNoGrad(Op):
+        """PyTensor Op which applies a PyLops Linear Operator, excluding gradient support.
+
+        This class "converts" a PyLops `LinearOperator` class into a PyTensor `Op`.
+        This applies the `LinearOperator` in "forward-mode" in `self.perform`.
+
+        Parameters
+        ----------
+        LOp : pylops.LinearOperator
+        """
+
+        __props__ = ("dims", "dimsd", "shape")
+
+        def __init__(self, LOp: pylops.LinearOperator) -> None:
+            self._LOp = LOp
+            self.dims = self._LOp.dims
+            self.dimsd = self._LOp.dimsd
+            self.shape = self._LOp.shape
+            super().__init__()
+
+        def make_node(self, x) -> Apply:
+            x = pt.as_tensor_variable(x)
+            inputs = [x]
+            outputs = [pt.tensor(dtype=x.type.dtype, shape=self._LOp.dimsd)]
+            return Apply(self, inputs, outputs)
+
+        def perform(
+            self, node: Apply, inputs: list, output_storage: list[list[None]]
+        ) -> None:
+            (x,) = inputs
+            (yt,) = output_storage
+            yt[0] = self._LOp @ x
+
+    class PyTensorOperator(_PyTensorOperatorNoGrad):
+        """PyTensor Op which applies a PyLops Linear Operator, including gradient support.
+
+        This class "converts" a PyLops `LinearOperator` class into a PyTensor `Op`.
+        This applies the `LinearOperator` in "forward-mode" in `self.perform`, and applies
+        its adjoint when computing the vector-Jacobian product (`self.grad`), as that is
+        the analytically correct gradient for linear operators. This class should pass
+        `pytensor.gradient.verify_grad`.
+
+        Parameters
+        ----------
+        LOp : pylops.LinearOperator
+        """
+
+        def __init__(self, LOp: pylops.LinearOperator) -> None:
+            super().__init__(LOp)
+            self._gradient_op = _PyTensorOperatorNoGrad(self._LOp.H)
+
+        def grad(
+            self, inputs: list[pt.TensorVariable], output_grads: list[pt.TensorVariable]
+        ):
+            return [self._gradient_op(output_grads[0])]
diff --git a/pylops/utils/deps.py b/pylops/utils/deps.py
index ecf69a95..df4a0d9e 100644
--- a/pylops/utils/deps.py
+++ b/pylops/utils/deps.py
@@ -10,6 +10,7 @@
     "spgl1_enabled",
     "sympy_enabled",
     "torch_enabled",
+    "pytensor_enabled",
 ]
 
 import os
@@ -223,6 +224,23 @@ def sympy_import(message: Optional[str] = None) -> str:
     return sympy_message
 
 
+def pytensor_import(message: Optional[str] = None) -> str:
+    if pytensor_enabled:
+        try:
+            import_module("pytensor")  # noqa: F401
+
+            pytensor_message = None
+        except Exception as e:
+            pytensor_message = f"Failed to import pytensor (error:{e})."
+    else:
+        pytensor_message = (
+            f"pytensor package not installed. In order to be able to use "
+            f"{message} run "
+            f'"pip install pytensor" or "conda install -c conda-forge pytensor".'
+        )
+    return pytensor_message
+
+
 # Set package availability booleans
 # cupy and jax: the package is imported to check everything is working correctly,
 # if not the package is disabled. We do this here as these libraries are used as drop-in
@@ -245,3 +263,4 @@ def sympy_import(message: Optional[str] = None) -> str:
 spgl1_enabled = util.find_spec("spgl1") is not None
 sympy_enabled = util.find_spec("sympy") is not None
 torch_enabled = util.find_spec("torch") is not None
+pytensor_enabled = util.find_spec("pytensor") is not None
diff --git a/pytests/test_pytensoroperator.py b/pytests/test_pytensoroperator.py
new file mode 100755
index 00000000..e86ddac3
--- /dev/null
+++ b/pytests/test_pytensoroperator.py
@@ -0,0 +1,51 @@
+import numpy as np
+import pytensor
+import pytest
+from numpy.testing import assert_array_equal
+
+from pylops import MatrixMult, PyTensorOperator
+
+par1 = {"ny": 11, "nx": 11, "dtype": np.float32}  # square
+par2 = {"ny": 21, "nx": 11, "dtype": np.float32}  # overdetermined
+
+np.random.seed(0)
+rng = np.random.default_rng()
+
+
+@pytest.mark.parametrize("par", [(par1)])
+def test_PyTensorOperator(par):
+    """Verify output and gradient of PyTensor function obtained from a LinearOperator."""
+    Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"])))
+    pytensor_op = PyTensorOperator(Dop)
+
+    # Check gradient
+    inp = np.random.randn(*pytensor_op.dims)
+    pytensor.gradient.verify_grad(pytensor_op, (inp,), rng=rng)
+
+    # Check value
+    x = pytensor.tensor.dvector()
+    f = pytensor.function([x], pytensor_op(x))
+    out = f(inp)
+    assert_array_equal(out, Dop @ inp)
+
+
+@pytest.mark.parametrize("par", [(par1)])
+def test_PyTensorOperator_nd(par):
+    """Verify output and gradient of PyTensor function obtained from a LinearOperator."""
+
+    otherdims = rng.choice(range(1, 3), size=rng.choice(range(2, 8)))
+    Dop = MatrixMult(
+        np.random.normal(0.0, 1.0, (par["ny"], par["nx"])), otherdims=otherdims
+    )
+    pytensor_op = PyTensorOperator(Dop)
+
+    # Check gradient
+    inp = np.random.randn(*pytensor_op.dims)
+    pytensor.gradient.verify_grad(pytensor_op, (inp,), rng=rng)
+
+    # Check value
+    tensor = pytensor.tensor.TensorType(dtype="float64", shape=(None,) * inp.ndim)
+    x = tensor()
+    f = pytensor.function([x], pytensor_op(x))
+    out = f(inp)
+    assert_array_equal(out, Dop @ inp)

From 434796f928b106ee5d1c7a0b7002666d6f7754f0 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 17:27:44 -0800
Subject: [PATCH 02/34] fix: add __pycache__ to gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 42f9cc53..0365f10d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -17,6 +17,7 @@ build
 dist
 pylops.egg-info/
 .eggs/
+__pycache__
 
 # setuptools_scm generated  #
 pylops/version.py

From 1ed3efeb4483228f9b8a60d6447d22f9589228d0 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 17:29:40 -0800
Subject: [PATCH 03/34] fix: docs

---
 pytests/test_pytensoroperator.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pytests/test_pytensoroperator.py b/pytests/test_pytensoroperator.py
index e86ddac3..968bc768 100755
--- a/pytests/test_pytensoroperator.py
+++ b/pytests/test_pytensoroperator.py
@@ -31,7 +31,8 @@ def test_PyTensorOperator(par):
 
 @pytest.mark.parametrize("par", [(par1)])
 def test_PyTensorOperator_nd(par):
-    """Verify output and gradient of PyTensor function obtained from a LinearOperator."""
+    """Verify output and gradient of PyTensor function obtained from a LinearOperator
+    using an ND-array."""
 
     otherdims = rng.choice(range(1, 3), size=rng.choice(range(2, 8)))
     Dop = MatrixMult(

From 68c415d96602c1d3872f6aab2a5bc15fa3e829b7 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 21:39:57 -0800
Subject: [PATCH 04/34] fix: add to docs

---
 environment-dev-arm.yml | 2 ++
 environment-dev.yml     | 2 ++
 requirements-dev.txt    | 2 ++
 requirements-doc.txt    | 2 ++
 4 files changed, 8 insertions(+)

diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml
index 0081a0ad..4fbc35b4 100755
--- a/environment-dev-arm.yml
+++ b/environment-dev-arm.yml
@@ -25,6 +25,8 @@ dependencies:
   - autopep8
   - isort
   - black
+  - pymc
+  - python-graphviz
   - pip:
     - devito
     - dtcwt
diff --git a/environment-dev.yml b/environment-dev.yml
index eb51c4dc..132d745c 100755
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -26,6 +26,8 @@ dependencies:
   - autopep8
   - isort
   - black
+  - pymc
+  - python-graphviz
   - pip:
     - devito
     - dtcwt
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 2fca40dd..4e2b0f01 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -28,3 +28,5 @@ isort
 black
 flake8
 mypy
+pymc
+graphviz
diff --git a/requirements-doc.txt b/requirements-doc.txt
index 2e62dea7..067ee37c 100644
--- a/requirements-doc.txt
+++ b/requirements-doc.txt
@@ -34,3 +34,5 @@ isort
 black
 flake8
 mypy
+pymc
+graphviz

From 66668048f1a663073f900fb21d86df541583990a Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 21:41:02 -0800
Subject: [PATCH 05/34] fix: formatting

---
 docs/source/gpu.rst                  |  4 ++--
 examples/plot_slopeest.py            |  6 +++---
 pylops/basicoperators/matrixmult.py  | 14 ++++++++------
 pylops/basicoperators/regression.py  |  2 +-
 pylops/basicoperators/restriction.py |  2 +-
 pylops/waveeqprocessing/marchenko.py | 17 +++++++++--------
 6 files changed, 24 insertions(+), 21 deletions(-)

diff --git a/docs/source/gpu.rst b/docs/source/gpu.rst
index ecc87120..a37597e9 100755
--- a/docs/source/gpu.rst
+++ b/docs/source/gpu.rst
@@ -462,5 +462,5 @@ Again, the code is almost unchanged apart from the fact that we now use ``jax``
 
 .. note::
 
-   More examples for the CuPy and JAX backends be found `here <https://github.com/PyLops/pylops_notebooks/tree/master/developement-cupy>`_
-   and `here <https://github.com/PyLops/pylops_notebooks/tree/master/developement/Basic_JAX.ipynb>`_.
\ No newline at end of file
+   More examples for the CuPy and JAX backends be found `here <https://github.com/PyLops/pylops_notebooks/tree/master/developement-cupy>`__
+   and `here <https://github.com/PyLops/pylops_notebooks/tree/master/developement/Basic_JAX.ipynb>`__.
\ No newline at end of file
diff --git a/examples/plot_slopeest.py b/examples/plot_slopeest.py
index 89938a4a..22cf9096 100755
--- a/examples/plot_slopeest.py
+++ b/examples/plot_slopeest.py
@@ -14,7 +14,7 @@
 precondition sparsity-promoting inverse problems.
 
 We will show examples of a variety of different settings, including a comparison
-with the original implementation in [1].
+with the original implementation in [1]_.
 
 .. [1] van Vliet, L. J.,  Verbeek, P. W., "Estimators for orientation and
     anisotropy in digitized images", Journal ASCI Imaging Workshop. 1995.
@@ -145,7 +145,7 @@
 ###############################################################################
 # Concentric circles
 # ------------------
-# The original paper by van Vliet and Verbeek [1] has an example with concentric
+# The original paper by van Vliet and Verbeek [1]_ has an example with concentric
 # circles. We recover their original images and compare our implementation with
 # theirs.
 
@@ -215,7 +215,7 @@ def rgb2gray(rgb):
 ###############################################################################
 # Core samples
 # ------------------
-# The original paper by van Vliet and Verbeek [1] also has an example with images
+# The original paper by van Vliet and Verbeek [1]_ also has an example with images
 # of core samples. Since the original paper does not have a scale with which to
 # plot the angles, we have chosen ours it to match their image as closely as
 # possible.
diff --git a/pylops/basicoperators/matrixmult.py b/pylops/basicoperators/matrixmult.py
index a5f713ab..66b96c70 100644
--- a/pylops/basicoperators/matrixmult.py
+++ b/pylops/basicoperators/matrixmult.py
@@ -79,12 +79,14 @@ def __init__(
         else:
             otherdims = _value_or_sized_to_array(otherdims)
             self.otherdims = np.array(otherdims, dtype=int)
-            dims, dimsd = np.insert(self.otherdims, 0, self.A.shape[1]), np.insert(
-                self.otherdims, 0, self.A.shape[0]
+            dims, dimsd = (
+                np.insert(self.otherdims, 0, self.A.shape[1]),
+                np.insert(self.otherdims, 0, self.A.shape[0]),
+            )
+            self.dimsflatten, self.dimsdflatten = (
+                np.insert([np.prod(self.otherdims)], 0, self.A.shape[1]),
+                np.insert([np.prod(self.otherdims)], 0, self.A.shape[0]),
             )
-            self.dimsflatten, self.dimsdflatten = np.insert(
-                [np.prod(self.otherdims)], 0, self.A.shape[1]
-            ), np.insert([np.prod(self.otherdims)], 0, self.A.shape[0])
             self.reshape = True
             explicit = False
 
@@ -138,7 +140,7 @@ def inv(self) -> NDArray:
         r"""Return the inverse of :math:`\mathbf{A}`.
 
         Returns
-        ----------
+        -------
         Ainv : :obj:`numpy.ndarray`
             Inverse matrix.
 
diff --git a/pylops/basicoperators/regression.py b/pylops/basicoperators/regression.py
index dc51ded5..1160fe9b 100644
--- a/pylops/basicoperators/regression.py
+++ b/pylops/basicoperators/regression.py
@@ -124,7 +124,7 @@ def apply(self, t: npt.ArrayLike, x: NDArray) -> NDArray:
             Regression coefficients
 
         Returns
-        ----------
+        -------
         y : :obj:`numpy.ndarray`
             Values along y-axis
 
diff --git a/pylops/basicoperators/restriction.py b/pylops/basicoperators/restriction.py
index c2fceab7..e27610bf 100644
--- a/pylops/basicoperators/restriction.py
+++ b/pylops/basicoperators/restriction.py
@@ -189,7 +189,7 @@ def mask(self, x: NDArray) -> NDArray:
             Input array (can be either flattened or not)
 
         Returns
-        ----------
+        -------
         y : :obj:`numpy.ma.core.MaskedArray`
             Masked array.
 
diff --git a/pylops/waveeqprocessing/marchenko.py b/pylops/waveeqprocessing/marchenko.py
index bf8e3305..d810c70b 100644
--- a/pylops/waveeqprocessing/marchenko.py
+++ b/pylops/waveeqprocessing/marchenko.py
@@ -301,7 +301,7 @@ def apply_onepoint(
         greens: bool = False,
         dottest: bool = False,
         usematmul: bool = False,
-        **kwargs_solver
+        **kwargs_solver,
     ) -> Union[
         Tuple[NDArray, NDArray, NDArray, NDArray, NDArray],
         Tuple[NDArray, NDArray, NDArray, NDArray],
@@ -341,7 +341,7 @@ def apply_onepoint(
             for numpy and cupy `data`, respectively)
 
         Returns
-        ----------
+        -------
         f1_inv_minus : :obj:`numpy.ndarray`
             Inverted upgoing focusing function of size :math:`[n_r \times n_t]`
         f1_inv_plus : :obj:`numpy.ndarray`
@@ -473,7 +473,7 @@ def apply_onepoint(
                 Mop,
                 d.ravel(),
                 x0=self.ncp.zeros(2 * (2 * self.nt - 1) * self.nr, dtype=self.dtype),
-                **kwargs_solver
+                **kwargs_solver,
             )[0]
 
         f1_inv = f1_inv.reshape(2 * self.nt2, self.nr)
@@ -486,8 +486,9 @@ def apply_onepoint(
             # Create Green's functions
             g_inv = Gop * f1_inv_tot.ravel()
             g_inv = g_inv.reshape(2 * self.nt2, self.ns)
-            g_inv_minus, g_inv_plus = -g_inv[: self.nt2].T, np.fliplr(
-                g_inv[self.nt2 :].T
+            g_inv_minus, g_inv_plus = (
+                -g_inv[: self.nt2].T,
+                np.fliplr(g_inv[self.nt2 :].T),
             )
         if rtm and greens:
             return f1_inv_minus, f1_inv_plus, p0_minus, g_inv_minus, g_inv_plus
@@ -507,7 +508,7 @@ def apply_multiplepoints(
         greens: bool = False,
         dottest: bool = False,
         usematmul: bool = False,
-        **kwargs_solver
+        **kwargs_solver,
     ) -> Union[
         Tuple[NDArray, NDArray, NDArray, NDArray, NDArray],
         Tuple[NDArray, NDArray, NDArray, NDArray],
@@ -548,7 +549,7 @@ def apply_multiplepoints(
             for numpy and cupy `data`, respectively)
 
         Returns
-        ----------
+        -------
         f1_inv_minus : :obj:`numpy.ndarray`
             Inverted upgoing focusing function of size
             :math:`[n_r \times n_{vs} \times n_t]`
@@ -695,7 +696,7 @@ def apply_multiplepoints(
                 x0=self.ncp.zeros(
                     2 * (2 * self.nt - 1) * self.nr * nvs, dtype=self.dtype
                 ),
-                **kwargs_solver
+                **kwargs_solver,
             )[0]
 
         f1_inv = f1_inv.reshape(2 * self.nt2, self.nr, nvs)

From edc0377de644e8a9ad47a735de99f2b5c06bafff Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 21:41:59 -0800
Subject: [PATCH 06/34] fix: remove options not available to
 pydata_sphinx_theme

---
 docs/source/conf.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/docs/source/conf.py b/docs/source/conf.py
index c5e6536d..7f2620a6 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -1,8 +1,10 @@
 # -*- coding: utf-8 -*-
-import sys
-import os
 import datetime
+import os
+import sys
+
 from sphinx_gallery.sorting import ExampleTitleSortKey
+
 from pylops import __version__
 
 # Sphinx needs to be able to import the package to use autodoc and get the version number
@@ -103,9 +105,7 @@
 # These enable substitutions using |variable| in the rst files
 rst_epilog = """
 .. |year| replace:: {year}
-""".format(
-    year=year
-)
+""".format(year=year)
 html_static_path = ["_static"]
 html_last_updated_fmt = "%b %d, %Y"
 html_title = "PyLops"
@@ -122,15 +122,15 @@
 # Theme config
 html_theme = "pydata_sphinx_theme"
 html_theme_options = {
-    "logo_only": True,
-    "display_version": True,
+    # "logo_only": True,
+    # "display_version": True,
     "logo": {
         "image_light": "pylops_b.png",
         "image_dark": "pylops.png",
     }
 }
 html_css_files = [
-    'css/custom.css',
+    "css/custom.css",
 ]
 
 html_context = {

From 06654008d7cffb540f5f5f1eeee74d23b840b4fd Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 21:48:30 -0800
Subject: [PATCH 07/34] fix: math formatting

---
 tutorials/torchop.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/tutorials/torchop.py b/tutorials/torchop.py
index 9e73d7b3..d13a764d 100755
--- a/tutorials/torchop.py
+++ b/tutorials/torchop.py
@@ -14,6 +14,7 @@
 modelling operators.
 
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 import torch
@@ -30,7 +31,7 @@
 # In this example we consider a simple multidimensional functional:
 #
 # .. math::
-#   \mathbf{y} = \mathbf{A} sin(\mathbf{x})
+#   \mathbf{y} = \mathbf{A} \sin(\mathbf{x})
 #
 # and we use AD to compute the gradient with respect to the input vector
 # evaluated at :math:`\mathbf{x}=\mathbf{x}_0` :
@@ -44,10 +45,10 @@
 #        ... & ... & ... \\
 #        dy_N / dx_1 & ... & dy_N / dx_M
 #        \end{bmatrix} = \begin{bmatrix}
-#        a_{11} cos(x_1) & ... & a_{1M} cos(x_M) \\
+#        a_{11} \cos(x_1) & ... & a_{1M} \cos(x_M) \\
 #        ... & ... & ... \\
-#        a_{N1} cos(x_1) & ... & a_{NM} cos(x_M)
-#        \end{bmatrix} = \textbf{A} cos(\mathbf{x})
+#        a_{N1} \cos(x_1) & ... & a_{NM} \cos(x_M)
+#        \end{bmatrix} = \textbf{A} \cos(\mathbf{x})
 #
 # Since both input and output are multidimensional,
 # PyTorch ``backward`` actually computes the product between the transposed

From 77a8cb28b5a205b5c23cfb293afda5a6b41516e7 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 21:52:56 -0800
Subject: [PATCH 08/34] fix: math formatting

---
 tutorials/torchop.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/tutorials/torchop.py b/tutorials/torchop.py
index d13a764d..a18ff2ef 100755
--- a/tutorials/torchop.py
+++ b/tutorials/torchop.py
@@ -35,19 +35,19 @@
 #
 # and we use AD to compute the gradient with respect to the input vector
 # evaluated at :math:`\mathbf{x}=\mathbf{x}_0` :
-# :math:`\mathbf{g} = d\mathbf{y} / d\mathbf{x} |_{\mathbf{x}=\mathbf{x}_0}`.
+# :math:`\mathbf{g} = \partial\mathbf{y} / \partial\mathbf{x} |_{\mathbf{x}=\mathbf{x}_0}`.
 #
 # Let's start by defining the Jacobian:
 #
 #   .. math::
 #        \textbf{J} = \begin{bmatrix}
-#        dy_1 / dx_1 & ... & dy_1 / dx_M \\
-#        ... & ... & ... \\
-#        dy_N / dx_1 & ... & dy_N / dx_M
+#        \frac{\partial y_1}{\partial x_1} & \cdots & \frac{\partial y_1}{\partial x_M} \\
+#        \vdots & \ddots & \vdots \\
+#        \frac{\partial y_N}{\partial x_1} & \cdots & \frac{\partial y_N}{\partial x_M}
 #        \end{bmatrix} = \begin{bmatrix}
-#        a_{11} \cos(x_1) & ... & a_{1M} \cos(x_M) \\
-#        ... & ... & ... \\
-#        a_{N1} \cos(x_1) & ... & a_{NM} \cos(x_M)
+#        a_{11} \cos(x_1) & \cdots & a_{1M} \cos(x_M) \\
+#        \vdots & \ddots & \vdots \\
+#        a_{N1} \cos(x_1) & \cdots & a_{NM} \cos(x_M)
 #        \end{bmatrix} = \textbf{A} \cos(\mathbf{x})
 #
 # Since both input and output are multidimensional,

From 50d8caa8419fc66a0c0ba73877666a565a2c61a4 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 23:00:44 -0800
Subject: [PATCH 09/34] fix: add intersphinx

---
 docs/source/conf.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/source/conf.py b/docs/source/conf.py
index 7f2620a6..75680cac 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -39,6 +39,8 @@
     "matplotlib": ("https://matplotlib.org/", None),
     "pyfftw": ("https://pyfftw.readthedocs.io/en/latest/", None),
     "spgl1": ("https://spgl1.readthedocs.io/en/latest/", None),
+    "pymc": ("https://www.pymc.io/", None),
+    "arviz": ("https://python.arviz.org/en/latest/", None),
 }
 
 # Generate autodoc stubs with summaries from code

From 59e40e67244e2c35e914cbf856e6aeb0b4909595 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 23:00:56 -0800
Subject: [PATCH 10/34] feature: Bayesian Linear Regression

---
 examples/plot_bayeslinearregr.py | 230 +++++++++++++++++++++++++++++++
 1 file changed, 230 insertions(+)
 create mode 100644 examples/plot_bayeslinearregr.py

diff --git a/examples/plot_bayeslinearregr.py b/examples/plot_bayeslinearregr.py
new file mode 100644
index 00000000..87d11081
--- /dev/null
+++ b/examples/plot_bayeslinearregr.py
@@ -0,0 +1,230 @@
+r"""
+Bayesian Linear Regression
+==========================
+
+In the :ref:`sphx_glr_gallery_plot_linearregr.py` example, we
+performed linear regression by applying a variety of solvers to the
+:py:class:`pylops.LinearRegression` operator.
+
+In this example, we will apply linear regression the Bayesian way.
+In Bayesian inference, we are not looking for a "best" estimate
+of the linear regression parameters; rather, we are looking for
+all possible parameters and their associated (posterior) probability,
+that is, how likely that those are the parameters that generated our data.
+
+To do this, we will leverage the probabilistic programming library
+`PyMC <https://www.pm.io>`_.
+
+In the Bayesian formulation, we write the problem in the following manner:
+
+    .. math::
+        y_i \sim N(x_0 + x_1 t_i, \sigma)  \qquad \forall i=0,1,\ldots,N-1
+
+where :math:`x_0` is the intercept and :math:`x_1` is the gradient.
+This notation means that the obtained measurements :math:`y_i` are normally distributed around
+mean :math:`x_0 + x_1 t_i` with a standard deviation of :math:`\sigma`.
+We can also express this problem in a matrix form, which makes it clear that we
+can use a PyLops operator to describe this relationship.
+
+    .. math::
+        \mathbf{y} \sim N(\mathbf{A} \mathbf{x}, \sigma)
+
+In this example, we will combine the Bayesian power of PyMC with the linear language of
+PyLops.
+"""
+
+from io import BytesIO
+
+import arviz as az
+import matplotlib.pyplot as plt
+import numpy as np
+import pymc as pm
+from PIL import Image
+
+import pylops
+
+plt.close("all")
+np.random.seed(10)
+
+###############################################################################
+# Define the input parameters: number of samples along the t-axis (``N``),
+# linear regression coefficients (``x``), and standard deviation of noise
+# to be added to data (``sigma``).
+N = 30
+x = np.array([1.0, 0.5])
+sigma = 0.25
+
+###############################################################################
+# Let's create the time axis and initialize the
+# :py:class:`pylops.LinearRegression` operator
+t = np.linspace(0, 1, N)
+LRop = pylops.LinearRegression(t, dtype=t.dtype)
+
+###############################################################################
+# We can then apply the operator in forward mode to compute our data points
+# along the x-axis (``y``). We will also generate some random gaussian noise
+# and create a noisy version of the data (``yn``).
+y = LRop @ x
+yn = y + np.random.normal(0, sigma, N)
+
+###############################################################################
+# The deterministic solution is to solve the
+# :math:`\mathbf{y} =  \mathbf{A} \mathbf{x}` in a least-squares sense.
+# Using PyLops, the ``/`` operator solves the iteratively (i.e.,
+# :py:func:`scipy.sparse.linalg.lsqr`).
+xnest = LRop / yn
+noise_est = np.sqrt(np.sum((yn - LRop @ xnest) ** 2) / (N - 1))
+
+###############################################################################
+# Let's plot the best fitting line for the case of noise free and noisy data
+fig, ax = plt.subplots(figsize=(8, 4))
+ax.plot(
+    np.array([t.min(), t.max()]),
+    np.array([t.min(), t.max()]) * x[1] + x[0],
+    "k",
+    lw=4,
+    label=rf"True: $x_0$ = {x[0]:.2f}, $x_1$ = {x[1]:.2f}",
+)
+ax.plot(
+    np.array([t.min(), t.max()]),
+    np.array([t.min(), t.max()]) * xnest[1] + xnest[0],
+    "--g",
+    lw=4,
+    label=rf"MAP Estimated: $x_0$ = {xnest[0]:.2f}, $x_1$ = {xnest[1]:.2f}",
+)
+ax.scatter(t, y, c="r", s=70)
+ax.scatter(t, yn, c="g", s=70)
+ax.legend()
+fig.tight_layout()
+
+###############################################################################
+# Let's solve this problem the Bayesian way, which consists in obtaining the
+# posterior probability :math:`p(\mathbf{x}\,|\,\mathbf{y})` via Bayes theorem:
+#
+# .. math::
+#    \underbrace{p(\mathbf{x} \,|\, \mathbf{y})}_{\text{posterior}}
+#    \propto \overbrace{p(\mathbf{y} \,|\, \mathbf{x})}^{\text{likelihood}}\;
+#    \overbrace{p(\mathbf{x})}^{\text{prior}}
+#
+# To do so, we need to define the priors and the likelihood.
+#
+# Priors in Bayesian analysis can be interpreted as the probabilistic
+# equivalent to regularization. Finding the maximum a posteriori (MAP) estimate
+# to a least-squares problem with a Gaussian prior on the parameters is
+# is equivalent to applying a Tikhonov (L2) regularization to these parameters.
+# A Laplace prior is equivalent to a sparse (L1) regularization. In addition,
+# the weight of the regularization is controlled by the "scale" of the
+# distribution of the prior; the standard deviation (in the case of a Gaussian)
+# is inversely proportional strength of the regularization.
+#
+# In this problem we will use weak, not very informative priors, by settings
+# their "weights" to be small (high scale parameters):
+#
+# .. math::
+#     x_0 \sim N(0, 20)
+#
+#     x_1 \sim N(0, 20)
+#
+#     \sigma \sim \text{HalfCauchy}(10)
+#
+# The (log) likelihood in Bayesian analysis is the equivalent of the cost
+# function in deterministic inverse problems. In this case we have already
+# seen this likelihood:
+#
+# .. math::
+#       p(\mathbf{y}\,|\,\mathbf{x}) \sim N(\mathbf{A}\mathbf{x}, \sigma)
+#
+
+# Construct a PyTensor `Op` which can be used in a PyMC model.
+pytensor_lrop = pylops.PyTensorOperator(LRop)
+dims = pytensor_lrop.dims  # Inherits dims, dimsd and shape from LRop
+
+# Construct the PyMC model
+with pm.Model() as model:
+    y_data = pm.Data("y_data", yn)
+
+    # Define priors
+    sp = pm.HalfCauchy("σ", beta=10)
+    xp = pm.Normal("x", 0, sigma=20, shape=dims)
+    mu = pm.Deterministic("mu", pytensor_lrop(xp))
+
+    # Define likelihood
+    likelihood = pm.Normal("y", mu=mu, sigma=sp, observed=y_data)
+
+    # Inference!
+    idata = pm.sample(2000, tune=1000, chains=2)
+
+###############################################################################
+# The PyMC library offers a way to visualize the Bayesian model which is
+# helpful in mapping dependencies.
+dot = pm.model_to_graphviz(model)
+
+# Some magic to display in docs
+plt.imshow(Image.open(BytesIO(dot.pipe("png"))))
+plt.axis("off")
+
+###############################################################################
+axes = az.plot_trace(idata, figsize=(10, 7), var_names=["~mu"])
+axes[0, 0].axvline(x[0], label="True Intercept", lw=2, color="k")
+axes[0, 0].axvline(xnest[0], label="Intercept MAP", lw=2, color="C0", ls="--")
+axes[0, 0].axvline(x[1], label="True Slope", lw=2, color="darkgray")
+axes[0, 0].axvline(xnest[1], label="Slope MAP", lw=2, color="C1", ls="--")
+axes[0, 1].axhline(x[0], label="True Intercept", lw=2, color="k")
+axes[0, 1].axhline(xnest[0], label="Intercept MAP", lw=2, color="C0", ls="--")
+axes[0, 1].axhline(x[1], label="True Slope", lw=2, color="darkgray")
+axes[0, 1].axhline(xnest[1], label="Slope MAP", lw=2, color="C1", ls="--")
+axes[1, 0].axvline(sigma, label="True Sigma", lw=2, color="k")
+axes[1, 0].axvline(noise_est, label="Sigma MAP", lw=2, color="C0", ls="--")
+axes[1, 1].axhline(sigma, label="True Sigma", lw=2, color="k")
+axes[1, 1].axhline(noise_est, label="Sigma MAP", lw=2, color="C0", ls="--")
+for ax in axes.ravel():
+    ax.legend()
+ax.get_figure().tight_layout()
+
+################################################################################
+# The plot above is known as the "trace" plot. The plots on the left display the
+# posterior distribution of all latent variables in the model. The top-left
+# has multiple colored, one for each parameter in the latent vector
+# :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
+# estimated noise.
+#
+# In these plots there are multiple distributions of the same color. Each
+# represents a "chain". A chain is a single run of a Monte Carlo algorithm.
+# Generally, Monte Carlo methods run various chains to ensure that all regions
+# of the posterior distribution are sampled. These chains are shown on the right
+# hand plots.
+#
+# With this model, we can obtain an uncertainty measurement via the High Density
+# Interval. To do that, we need to sample the "preditive posterior", that is,
+# the posterior distribution of the data, given the model.
+
+with model:
+    pm.sample_posterior_predictive(idata, extend_inferencedata=True)
+
+###############################################################################
+fig, ax = plt.subplots(figsize=(8, 4))
+az.plot_hdi(
+    t,
+    idata.posterior_predictive["y"],
+    fill_kwargs={"label": "95% HDI"},
+    hdi_prob=0.95,
+    ax=ax,
+)
+ax.plot(
+    np.array([t.min(), t.max()]),
+    np.array([t.min(), t.max()]) * x[1] + x[0],
+    "k",
+    lw=4,
+    label=rf"True: $x_0$ = {x[0]:.2f}, $x_1$ = {x[1]:.2f}",
+)
+ax.plot(
+    np.array([t.min(), t.max()]),
+    np.array([t.min(), t.max()]) * xnest[1] + xnest[0],
+    "--g",
+    lw=4,
+    label=rf"MAP Estimated: $x_0$ = {xnest[0]:.2f}, $x_1$ = {xnest[1]:.2f}",
+)
+ax.scatter(t, y, c="r", s=70)
+ax.scatter(t, yn, c="g", s=70)
+ax.legend()
+fig.tight_layout()

From d68ee794964978852974458f00f9d55888b2880b Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sat, 16 Nov 2024 23:13:09 -0800
Subject: [PATCH 11/34] fix: change thumbnail

---
 examples/plot_bayeslinearregr.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/plot_bayeslinearregr.py b/examples/plot_bayeslinearregr.py
index 87d11081..5a0123f6 100644
--- a/examples/plot_bayeslinearregr.py
+++ b/examples/plot_bayeslinearregr.py
@@ -202,6 +202,7 @@
     pm.sample_posterior_predictive(idata, extend_inferencedata=True)
 
 ###############################################################################
+# sphinx_gallery_thumbnail_number = 4
 fig, ax = plt.subplots(figsize=(8, 4))
 az.plot_hdi(
     t,

From 571f6d5828613e006e45ecf3c6628e7f4dc66a4d Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:14:46 -0800
Subject: [PATCH 12/34] fix: use == to compare string literal

---
 pytests/test_oneway.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pytests/test_oneway.py b/pytests/test_oneway.py
index 48f73a9e..1f1b936c 100755
--- a/pytests/test_oneway.py
+++ b/pytests/test_oneway.py
@@ -96,7 +96,7 @@ def test_PhaseShift_3dsignal(par):
 @pytest.mark.parametrize("par", [(par1), (par2), (par1v), (par2v)])
 def test_Deghosting_2dsignal(par, create_data2D):
     """Deghosting of 2d data"""
-    p2d, p2d_minus = create_data2D(1 if par["kind"] is "p" else -1)
+    p2d, p2d_minus = create_data2D(1 if par["kind"] == "p" else -1)
 
     p2d_minus_inv, p2d_plus_inv = Deghosting(
         p2d,
@@ -111,7 +111,7 @@ def test_Deghosting_2dsignal(par, create_data2D):
         npad=0,
         ntaper=0,
         dtype=par["dtype"],
-        **dict(damp=1e-10, iter_lim=60)
+        **dict(damp=1e-10, iter_lim=60),
     )
 
     assert np.linalg.norm(p2d_minus_inv - p2d_minus) / np.linalg.norm(p2d_minus) < 3e-1

From 004bee8153cee6bf01c4d6555e328c2dbad8837d Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:24:02 -0800
Subject: [PATCH 13/34] fix: do not test on mac

---
 pytests/test_pytensoroperator.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/pytests/test_pytensoroperator.py b/pytests/test_pytensoroperator.py
index 968bc768..05cf30b4 100755
--- a/pytests/test_pytensoroperator.py
+++ b/pytests/test_pytensoroperator.py
@@ -1,3 +1,5 @@
+import platform
+
 import numpy as np
 import pytensor
 import pytest
@@ -15,6 +17,10 @@
 @pytest.mark.parametrize("par", [(par1)])
 def test_PyTensorOperator(par):
     """Verify output and gradient of PyTensor function obtained from a LinearOperator."""
+    # See: https://discourse.pymc.io/t/installation-issues-v5-9-macos/13094
+    if platform.system() == "Darwin":
+        return
+
     Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"])))
     pytensor_op = PyTensorOperator(Dop)
 
@@ -33,6 +39,9 @@ def test_PyTensorOperator(par):
 def test_PyTensorOperator_nd(par):
     """Verify output and gradient of PyTensor function obtained from a LinearOperator
     using an ND-array."""
+    # See: https://discourse.pymc.io/t/installation-issues-v5-9-macos/13094
+    if platform.system() == "Darwin":
+        return
 
     otherdims = rng.choice(range(1, 3), size=rng.choice(range(2, 8)))
     Dop = MatrixMult(

From df9f8e5660b56dfa2680c442c647c3fe3666135f Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:36:56 -0800
Subject: [PATCH 14/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 4e2b0f01..f9a78517 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -14,7 +14,7 @@ ipython
 pytest
 pytest-runner
 setuptools_scm
-docutils<0.18
+docutils>=0.18
 Sphinx
 pydata-sphinx-theme
 sphinx-gallery
@@ -28,5 +28,6 @@ isort
 black
 flake8
 mypy
-pymc
+pymc>=5
+pandas>=2
 graphviz

From e8a713d76860b55bfdc8187a2b071d6ef639bf23 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:40:11 -0800
Subject: [PATCH 15/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index f9a78517..df84cef4 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -30,4 +30,6 @@ flake8
 mypy
 pymc>=5
 pandas>=2
+arviz>=0.20
+contourpy>=1.3
 graphviz

From 8ca8ece3bed84a6d037e32faa30680ba1d783578 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:44:29 -0800
Subject: [PATCH 16/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index df84cef4..57fb442b 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,6 +1,6 @@
-numpy>=1.21.0
+numpy>=1.21.0,<2
 scipy>=1.11.0
-jax
+jax>=0.4
 numba
 pyfftw
 PyWavelets

From 49b25459c753c8c16e47f641af10b8b99f4d30ea Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:45:57 -0800
Subject: [PATCH 17/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 57fb442b..dd43be22 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -29,6 +29,7 @@ black
 flake8
 mypy
 pymc>=5
+pytensor>=2.23
 pandas>=2
 arviz>=0.20
 contourpy>=1.3

From 24ad574dd19009c229ccf2aad0cd85a0334a19ca Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:50:03 -0800
Subject: [PATCH 18/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index dd43be22..3e7f7b30 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -6,7 +6,7 @@ pyfftw
 PyWavelets
 spgl1
 scikit-fmm
-sympy
+sympy=0.12
 devito
 dtcwt
 matplotlib

From 2dca9d9d29ccf558e5d5e8cf421086c00b8bf1ef Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:50:19 -0800
Subject: [PATCH 19/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 3e7f7b30..b8136d26 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -6,7 +6,7 @@ pyfftw
 PyWavelets
 spgl1
 scikit-fmm
-sympy=0.12
+sympy==0.12
 devito
 dtcwt
 matplotlib

From aa230dd974d703c72b46e1dbc0feb5931e9e11cf Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:51:01 -0800
Subject: [PATCH 20/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index b8136d26..4a2924a1 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -6,7 +6,7 @@ pyfftw
 PyWavelets
 spgl1
 scikit-fmm
-sympy==0.12
+sympy<1.13
 devito
 dtcwt
 matplotlib

From 8117fca2b2e6c986dd1dc804aa317a1143dc183e Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 11:58:29 -0800
Subject: [PATCH 21/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 4a2924a1..80456f11 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -6,8 +6,8 @@ pyfftw
 PyWavelets
 spgl1
 scikit-fmm
-sympy<1.13
-devito
+sympy
+devito>=4.8.11
 dtcwt
 matplotlib
 ipython

From 4720b5edb4a86c9cb9a1cc55317b4d637d633bec Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 12:02:05 -0800
Subject: [PATCH 22/34] fix: test new versions for dev deps

---
 requirements-dev.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 80456f11..16eb77ec 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -3,7 +3,7 @@ scipy>=1.11.0
 jax>=0.4
 numba
 pyfftw
-PyWavelets
+PyWavelets>=1.7
 spgl1
 scikit-fmm
 sympy

From 038351f00ab37603d6cd0c8a0f8375bca0296e3d Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 12:06:25 -0800
Subject: [PATCH 23/34] lock all

---
 requirements-dev.txt | 72 ++++++++++++++++++++++----------------------
 1 file changed, 36 insertions(+), 36 deletions(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 16eb77ec..b4d547e4 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,36 +1,36 @@
-numpy>=1.21.0,<2
-scipy>=1.11.0
-jax>=0.4
-numba
-pyfftw
-PyWavelets>=1.7
-spgl1
-scikit-fmm
-sympy
-devito>=4.8.11
-dtcwt
-matplotlib
-ipython
-pytest
-pytest-runner
-setuptools_scm
-docutils>=0.18
-Sphinx
-pydata-sphinx-theme
-sphinx-gallery
-sphinxemoji
-numpydoc
-nbsphinx
-image
-pre-commit
-autopep8
-isort
-black
-flake8
-mypy
-pymc>=5
-pytensor>=2.23
-pandas>=2
-arviz>=0.20
-contourpy>=1.3
-graphviz
+numpy==1.26.3
+scipy==1.14.1
+jax==0.4.35
+numba==0.60.0
+pyfftw==0.15.0
+PyWavelets==1.7.0
+spgl1==0.0.3
+scikit-fmm==2024.9.16
+sympy==1.13.3
+devito==4.8.11
+dtcwt==0.14.0
+matplotlib==3.10.0rc1
+ipython==8.29.0
+pytest==8.3.3
+pytest-runner==6.0.1
+setuptools_scm==8.1.0
+docutils==0.21.2
+Sphinx==8.1.3
+pydata-sphinx-theme==0.16.0
+sphinx-gallery==0.18.0
+sphinxemoji==0.3.1
+numpydoc==1.8.0
+nbsphinx==0.9.5
+image==1.5.33
+pre-commit==4.0.1
+autopep8==2.3.1
+isort==6.0.0b2
+black==24.10.0
+flake8==7.1.1
+mypy==1.13.0
+pymc==5.18.1
+pytensor==2.26.3
+pandas==2.2.3
+arviz==0.20.0
+contourpy==1.3.1
+graphviz==0.20.3

From d5c772c9cdaede77abc1568f0305528717762b11 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 14:40:18 -0800
Subject: [PATCH 24/34] fix: test new versions for dev deps

---
 requirements-dev.txt   | 68 ++++++++++++++++++++----------------------
 requirements-torch.txt |  2 +-
 2 files changed, 33 insertions(+), 37 deletions(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index b4d547e4..bf5bd639 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,36 +1,32 @@
-numpy==1.26.3
-scipy==1.14.1
-jax==0.4.35
-numba==0.60.0
-pyfftw==0.15.0
-PyWavelets==1.7.0
-spgl1==0.0.3
-scikit-fmm==2024.9.16
-sympy==1.13.3
-devito==4.8.11
-dtcwt==0.14.0
-matplotlib==3.10.0rc1
-ipython==8.29.0
-pytest==8.3.3
-pytest-runner==6.0.1
-setuptools_scm==8.1.0
-docutils==0.21.2
-Sphinx==8.1.3
-pydata-sphinx-theme==0.16.0
-sphinx-gallery==0.18.0
-sphinxemoji==0.3.1
-numpydoc==1.8.0
-nbsphinx==0.9.5
-image==1.5.33
-pre-commit==4.0.1
-autopep8==2.3.1
-isort==6.0.0b2
-black==24.10.0
-flake8==7.1.1
-mypy==1.13.0
-pymc==5.18.1
-pytensor==2.26.3
-pandas==2.2.3
-arviz==0.20.0
-contourpy==1.3.1
-graphviz==0.20.3
+numpy>=1.21.0,<2
+scipy>=1.11.0
+jax
+numba
+pyfftw
+PyWavelets
+spgl1
+scikit-fmm
+sympy
+devito
+dtcwt
+matplotlib
+ipython
+pytest
+pytest-runner
+setuptools_scm
+docutils
+Sphinx
+pydata-sphinx-theme
+sphinx-gallery
+sphinxemoji
+numpydoc
+nbsphinx
+image
+pre-commit
+autopep8
+isort
+black
+flake8
+mypy
+pymc
+graphviz
diff --git a/requirements-torch.txt b/requirements-torch.txt
index f2c3b105..3f94f19f 100644
--- a/requirements-torch.txt
+++ b/requirements-torch.txt
@@ -1,2 +1,2 @@
 --index-url https://download.pytorch.org/whl/cpu
-torch>=1.2.0
+torch>=1.2.0,<2.5

From 48477a5e58997f74a2737f6d89c49c702dd8ae25 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 14:52:21 -0800
Subject: [PATCH 25/34] fix: test on darwin

---
 pytests/test_pytensoroperator.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pytests/test_pytensoroperator.py b/pytests/test_pytensoroperator.py
index 05cf30b4..f5e9bf97 100755
--- a/pytests/test_pytensoroperator.py
+++ b/pytests/test_pytensoroperator.py
@@ -18,8 +18,8 @@
 def test_PyTensorOperator(par):
     """Verify output and gradient of PyTensor function obtained from a LinearOperator."""
     # See: https://discourse.pymc.io/t/installation-issues-v5-9-macos/13094
-    if platform.system() == "Darwin":
-        return
+    # if platform.system() == "Darwin":
+    #     return
 
     Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"])))
     pytensor_op = PyTensorOperator(Dop)
@@ -40,8 +40,8 @@ def test_PyTensorOperator_nd(par):
     """Verify output and gradient of PyTensor function obtained from a LinearOperator
     using an ND-array."""
     # See: https://discourse.pymc.io/t/installation-issues-v5-9-macos/13094
-    if platform.system() == "Darwin":
-        return
+    # if platform.system() == "Darwin":
+    #     return
 
     otherdims = rng.choice(range(1, 3), size=rng.choice(range(2, 8)))
     Dop = MatrixMult(

From 85e82532ee5408c9f0adefd32bfa88ea8dfde3d0 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 14:53:36 -0800
Subject: [PATCH 26/34] fix: revert to old docutils

---
 requirements-dev.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index bf5bd639..0b924c1c 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -14,7 +14,7 @@ ipython
 pytest
 pytest-runner
 setuptools_scm
-docutils
+docutils<0.18
 Sphinx
 pydata-sphinx-theme
 sphinx-gallery

From 23fea350300a629c497329108065fdb7b0f4f3a1 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 15:08:41 -0800
Subject: [PATCH 27/34] fix: passing on mac

---
 pytests/test_pytensoroperator.py | 10 ----------
 1 file changed, 10 deletions(-)

diff --git a/pytests/test_pytensoroperator.py b/pytests/test_pytensoroperator.py
index f5e9bf97..9a59bc88 100755
--- a/pytests/test_pytensoroperator.py
+++ b/pytests/test_pytensoroperator.py
@@ -1,5 +1,3 @@
-import platform
-
 import numpy as np
 import pytensor
 import pytest
@@ -17,10 +15,6 @@
 @pytest.mark.parametrize("par", [(par1)])
 def test_PyTensorOperator(par):
     """Verify output and gradient of PyTensor function obtained from a LinearOperator."""
-    # See: https://discourse.pymc.io/t/installation-issues-v5-9-macos/13094
-    # if platform.system() == "Darwin":
-    #     return
-
     Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"])))
     pytensor_op = PyTensorOperator(Dop)
 
@@ -39,10 +33,6 @@ def test_PyTensorOperator(par):
 def test_PyTensorOperator_nd(par):
     """Verify output and gradient of PyTensor function obtained from a LinearOperator
     using an ND-array."""
-    # See: https://discourse.pymc.io/t/installation-issues-v5-9-macos/13094
-    # if platform.system() == "Darwin":
-    #     return
-
     otherdims = rng.choice(range(1, 3), size=rng.choice(range(2, 8)))
     Dop = MatrixMult(
         np.random.normal(0.0, 1.0, (par["ny"], par["nx"])), otherdims=otherdims

From 33677d6c170e263a30355163cda87ff8f8ac3e3b Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 16:30:08 -0800
Subject: [PATCH 28/34] fix: bump arviz version

---
 requirements-doc.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/requirements-doc.txt b/requirements-doc.txt
index 067ee37c..5905ff1d 100644
--- a/requirements-doc.txt
+++ b/requirements-doc.txt
@@ -35,4 +35,5 @@ black
 flake8
 mypy
 pymc
+arviz>=0.18
 graphviz

From fd6d703ac3ea64a00e26acf23c239fe55aa249dd Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 16:33:14 -0800
Subject: [PATCH 29/34] fix: use old scipy

---
 requirements-doc.txt | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/requirements-doc.txt b/requirements-doc.txt
index 5905ff1d..ed570d2b 100644
--- a/requirements-doc.txt
+++ b/requirements-doc.txt
@@ -3,7 +3,7 @@
 # same reason, we force devito==4.8.7 as later versions of devito
 # require numpy>=2.0.0
 numpy>=1.21.0,<2.0.0
-scipy>=1.11.0
+scipy>=1.11.0,<1.13
 jax
 --extra-index-url https://download.pytorch.org/whl/cpu
 torch>=1.2.0
@@ -35,5 +35,5 @@ black
 flake8
 mypy
 pymc
-arviz>=0.18
+arviz
 graphviz

From 88648c4e1c7b59a364df7d46891e7cd23d666d23 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 16:52:06 -0800
Subject: [PATCH 30/34] fix: remove gtraphviz (requies binary), use pytensor
 instead of pymc for dev dependencies, include pytensor as doc dependency

---
 environment-dev-arm.yml          |  3 +--
 environment-dev.yml              |  3 +--
 examples/plot_bayeslinearregr.py | 40 ++++++++++++++------------------
 requirements-dev.txt             |  3 +--
 requirements-doc.txt             |  3 +--
 5 files changed, 21 insertions(+), 31 deletions(-)

diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml
index 4fbc35b4..acda12aa 100755
--- a/environment-dev-arm.yml
+++ b/environment-dev-arm.yml
@@ -25,8 +25,7 @@ dependencies:
   - autopep8
   - isort
   - black
-  - pymc
-  - python-graphviz
+  - pytensor
   - pip:
     - devito
     - dtcwt
diff --git a/environment-dev.yml b/environment-dev.yml
index 132d745c..d29d25b5 100755
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -26,8 +26,7 @@ dependencies:
   - autopep8
   - isort
   - black
-  - pymc
-  - python-graphviz
+  - pytensor
   - pip:
     - devito
     - dtcwt
diff --git a/examples/plot_bayeslinearregr.py b/examples/plot_bayeslinearregr.py
index 5a0123f6..27e76e05 100644
--- a/examples/plot_bayeslinearregr.py
+++ b/examples/plot_bayeslinearregr.py
@@ -33,13 +33,10 @@
 PyLops.
 """
 
-from io import BytesIO
-
 import arviz as az
 import matplotlib.pyplot as plt
 import numpy as np
 import pymc as pm
-from PIL import Image
 
 import pylops
 
@@ -155,15 +152,24 @@
     idata = pm.sample(2000, tune=1000, chains=2)
 
 ###############################################################################
-# The PyMC library offers a way to visualize the Bayesian model which is
-# helpful in mapping dependencies.
-dot = pm.model_to_graphviz(model)
-
-# Some magic to display in docs
-plt.imshow(Image.open(BytesIO(dot.pipe("png"))))
-plt.axis("off")
+# The plot above is known as the "trace" plot. The plots on the left display the
+# posterior distribution of all latent variables in the model. The top-left
+# has multiple colored, one for each parameter in the latent vector
+# :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
+# estimated noise.
+#
+# The plot above is known as the "trace" plot. The plots on the left display the
+# posterior distribution of all latent variables in the model. The top-left
+# has multiple colored, one for each parameter in the latent vector
+# :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
+# estimated noise.
+#
+# In these plots there are multiple distributions of the same color. Each
+# represents a "chain". A chain is a single run of a Monte Carlo algorithm.
+# Generally, Monte Carlo methods run various chains to ensure that all regions
+# of the posterior distribution are sampled. These chains are shown on the right
+# hand plots.
 
-###############################################################################
 axes = az.plot_trace(idata, figsize=(10, 7), var_names=["~mu"])
 axes[0, 0].axvline(x[0], label="True Intercept", lw=2, color="k")
 axes[0, 0].axvline(xnest[0], label="Intercept MAP", lw=2, color="C0", ls="--")
@@ -182,18 +188,6 @@
 ax.get_figure().tight_layout()
 
 ################################################################################
-# The plot above is known as the "trace" plot. The plots on the left display the
-# posterior distribution of all latent variables in the model. The top-left
-# has multiple colored, one for each parameter in the latent vector
-# :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
-# estimated noise.
-#
-# In these plots there are multiple distributions of the same color. Each
-# represents a "chain". A chain is a single run of a Monte Carlo algorithm.
-# Generally, Monte Carlo methods run various chains to ensure that all regions
-# of the posterior distribution are sampled. These chains are shown on the right
-# hand plots.
-#
 # With this model, we can obtain an uncertainty measurement via the High Density
 # Interval. To do that, we need to sample the "preditive posterior", that is,
 # the posterior distribution of the data, given the model.
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 0b924c1c..8eb9f87d 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -28,5 +28,4 @@ isort
 black
 flake8
 mypy
-pymc
-graphviz
+pytensor
diff --git a/requirements-doc.txt b/requirements-doc.txt
index ed570d2b..e13ad06e 100644
--- a/requirements-doc.txt
+++ b/requirements-doc.txt
@@ -34,6 +34,5 @@ isort
 black
 flake8
 mypy
+pytensor
 pymc
-arviz
-graphviz

From 9a83bf52a95e2275b0a9b4b69953374b00489068 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 17:10:10 -0800
Subject: [PATCH 31/34] fix: improve descriptions, fix thumbnail

---
 examples/plot_bayeslinearregr.py | 34 +++++++++++++++-----------------
 1 file changed, 16 insertions(+), 18 deletions(-)

diff --git a/examples/plot_bayeslinearregr.py b/examples/plot_bayeslinearregr.py
index 27e76e05..5faacd5e 100644
--- a/examples/plot_bayeslinearregr.py
+++ b/examples/plot_bayeslinearregr.py
@@ -152,23 +152,17 @@
     idata = pm.sample(2000, tune=1000, chains=2)
 
 ###############################################################################
-# The plot above is known as the "trace" plot. The plots on the left display the
-# posterior distribution of all latent variables in the model. The top-left
-# has multiple colored, one for each parameter in the latent vector
-# :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
-# estimated noise.
+# The plot below is known as the "trace" plot. The left column displays the
+# posterior distributions of all latent variables in the model. The top-left
+# plot has multiple colored posteriors, one for each parameter of the latent
+# vector :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
+# estimated noise :math:`\sigma`.
 #
-# The plot above is known as the "trace" plot. The plots on the left display the
-# posterior distribution of all latent variables in the model. The top-left
-# has multiple colored, one for each parameter in the latent vector
-# :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
-# estimated noise.
-#
-# In these plots there are multiple distributions of the same color. Each
-# represents a "chain". A chain is a single run of a Monte Carlo algorithm.
-# Generally, Monte Carlo methods run various chains to ensure that all regions
-# of the posterior distribution are sampled. These chains are shown on the right
-# hand plots.
+# In these plots there are multiple distributions of the same color and multipl
+# line styles. Each of these represents a "chain". A chain is a single run of
+# a Monte Carlo algorithm. Generally, Monte Carlo methods run various chains
+# to ensure that all regions of the posterior distribution are sampled. These
+# chains are shown on the right hand plots.
 
 axes = az.plot_trace(idata, figsize=(10, 7), var_names=["~mu"])
 axes[0, 0].axvline(x[0], label="True Intercept", lw=2, color="k")
@@ -190,13 +184,17 @@
 ################################################################################
 # With this model, we can obtain an uncertainty measurement via the High Density
 # Interval. To do that, we need to sample the "preditive posterior", that is,
-# the posterior distribution of the data, given the model.
+# the posterior distribution of the data, given the model. What this does is
+# sample the latent vetors from their posteriors (above), and use the model
+# to construct realizations of the data given these realizations. They represent
+# what the model thinks the data should look like, given everything it has
+# already seen.
 
 with model:
     pm.sample_posterior_predictive(idata, extend_inferencedata=True)
 
 ###############################################################################
-# sphinx_gallery_thumbnail_number = 4
+# sphinx_gallery_thumbnail_number = 3
 fig, ax = plt.subplots(figsize=(8, 4))
 az.plot_hdi(
     t,

From 80bf466a5e9cc23c9f40bf0a4ef7d40d3500c7c4 Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Sun, 17 Nov 2024 17:12:53 -0800
Subject: [PATCH 32/34] fix: typo

---
 examples/plot_bayeslinearregr.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/examples/plot_bayeslinearregr.py b/examples/plot_bayeslinearregr.py
index 5faacd5e..b6f215aa 100644
--- a/examples/plot_bayeslinearregr.py
+++ b/examples/plot_bayeslinearregr.py
@@ -158,11 +158,11 @@
 # vector :math:`\mathbf{x}`. The bottom left plot displays the posterior of the
 # estimated noise :math:`\sigma`.
 #
-# In these plots there are multiple distributions of the same color and multipl
-# line styles. Each of these represents a "chain". A chain is a single run of
-# a Monte Carlo algorithm. Generally, Monte Carlo methods run various chains
-# to ensure that all regions of the posterior distribution are sampled. These
-# chains are shown on the right hand plots.
+# In these plots there are multiple distributions of the same color and
+# multiple line styles. Each of these represents a "chain". A chain is a single
+# run of a Monte Carlo algorithm. Generally, Monte Carlo methods run various
+# chains to ensure that all regions of the posterior distribution are sampled.
+# These chains are shown on the right hand plots.
 
 axes = az.plot_trace(idata, figsize=(10, 7), var_names=["~mu"])
 axes[0, 0].axvline(x[0], label="True Intercept", lw=2, color="k")

From 956e74eea3312d284038a6687e28f9cbad6d96cc Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Mon, 18 Nov 2024 18:41:23 -0800
Subject: [PATCH 33/34] fix: improve workinng, add MAP, make chains shorter

---
 examples/plot_bayeslinearregr.py | 108 +++++++++++++++++--------------
 1 file changed, 60 insertions(+), 48 deletions(-)

diff --git a/examples/plot_bayeslinearregr.py b/examples/plot_bayeslinearregr.py
index b6f215aa..42417eb1 100644
--- a/examples/plot_bayeslinearregr.py
+++ b/examples/plot_bayeslinearregr.py
@@ -69,26 +69,35 @@
 # :math:`\mathbf{y} =  \mathbf{A} \mathbf{x}` in a least-squares sense.
 # Using PyLops, the ``/`` operator solves the iteratively (i.e.,
 # :py:func:`scipy.sparse.linalg.lsqr`).
-xnest = LRop / yn
-noise_est = np.sqrt(np.sum((yn - LRop @ xnest) ** 2) / (N - 1))
+# In Bayesian terminology, this estimator is known as the maximulum likelihood
+# estimation (MLE).
+x_mle = LRop / yn
+noise_mle = np.sqrt(np.sum((yn - LRop @ x_mle) ** 2) / (N - 1))
+
+###############################################################################
+# Alternatively, we may regularize the problem. In this case we will condition
+# the solution towards smaller magnitude parameters, we can use a regularized
+# least squares approach. Since the weight is pretty small, we expect the
+# result to be very similar to the one above.
+sigma_prior = 20
+eps = 1 / np.sqrt(2) / sigma_prior
+x_map, *_ = pylops.optimization.basic.lsqr(LRop, yn, damp=eps)
+noise_map = np.sqrt(np.sum((yn - LRop @ x_map) ** 2) / (N - 1))
 
 ###############################################################################
 # Let's plot the best fitting line for the case of noise free and noisy data
 fig, ax = plt.subplots(figsize=(8, 4))
-ax.plot(
-    np.array([t.min(), t.max()]),
-    np.array([t.min(), t.max()]) * x[1] + x[0],
-    "k",
-    lw=4,
-    label=rf"True: $x_0$ = {x[0]:.2f}, $x_1$ = {x[1]:.2f}",
-)
-ax.plot(
-    np.array([t.min(), t.max()]),
-    np.array([t.min(), t.max()]) * xnest[1] + xnest[0],
-    "--g",
-    lw=4,
-    label=rf"MAP Estimated: $x_0$ = {xnest[0]:.2f}, $x_1$ = {xnest[1]:.2f}",
-)
+for est, est_label, c in zip(
+    [x, x_mle, x_map], ["True", "MLE", "MAP"], ["k", "C0", "C1"]
+):
+    ax.plot(
+        np.array([t.min(), t.max()]),
+        np.array([t.min(), t.max()]) * est[1] + est[0],
+        color=c,
+        ls="--" if est_label == "MAP" else "-",
+        lw=4,
+        label=rf"{est_label}: $x_0$ = {est[0]:.2f}, $x_1$ = {est[1]:.2f}",
+    )
 ax.scatter(t, y, c="r", s=70)
 ax.scatter(t, yn, c="g", s=70)
 ax.legend()
@@ -105,17 +114,23 @@
 #
 # To do so, we need to define the priors and the likelihood.
 #
-# Priors in Bayesian analysis can be interpreted as the probabilistic
-# equivalent to regularization. Finding the maximum a posteriori (MAP) estimate
-# to a least-squares problem with a Gaussian prior on the parameters is
-# is equivalent to applying a Tikhonov (L2) regularization to these parameters.
-# A Laplace prior is equivalent to a sparse (L1) regularization. In addition,
-# the weight of the regularization is controlled by the "scale" of the
-# distribution of the prior; the standard deviation (in the case of a Gaussian)
-# is inversely proportional strength of the regularization.
+# As hinted above, priors in Bayesian analysis can be interpreted as the
+# probabilistic equivalent to regularization. Finding the maximum a posteriori
+# (MAP) estimate to a least-squares problem with a Gaussian prior on the
+# parameters is equivalent to applying a Tikhonov (L2) regularization to these
+# parameters. A Laplace prior is equivalent to a sparse (L1) regularization.
+# In addition, the weight of the regularization is controlled by the "scale" of
+# the distribution of the prior; the standard deviation (in the case of a Gaussian)
+# is inversely proportional strength of the regularization. So if we use the same
+# sigma_prior above as the standard deviation of our prior distribition, we
+# should get the same MAP out of them. In practice, in Bayesian analysis we are
+# not only interested in point estimates like MAP, but rather, the whole
+# posterior distribution. If you want the MAP only, there are better,
+# methods to obtain them, such as the one shown above.
 #
-# In this problem we will use weak, not very informative priors, by settings
-# their "weights" to be small (high scale parameters):
+# In this problem we will use weak, not very informative priors, by setting
+# their prior to accept a wide range of probable values. This is equivalent to
+# setting the "weights" to be small, as shown above:
 #
 # .. math::
 #     x_0 \sim N(0, 20)
@@ -142,14 +157,14 @@
 
     # Define priors
     sp = pm.HalfCauchy("σ", beta=10)
-    xp = pm.Normal("x", 0, sigma=20, shape=dims)
+    xp = pm.Normal("x", 0, sigma=sigma_prior, shape=dims)
     mu = pm.Deterministic("mu", pytensor_lrop(xp))
 
     # Define likelihood
     likelihood = pm.Normal("y", mu=mu, sigma=sp, observed=y_data)
 
     # Inference!
-    idata = pm.sample(2000, tune=1000, chains=2)
+    idata = pm.sample(500, tune=200, chains=2)
 
 ###############################################################################
 # The plot below is known as the "trace" plot. The left column displays the
@@ -166,17 +181,17 @@
 
 axes = az.plot_trace(idata, figsize=(10, 7), var_names=["~mu"])
 axes[0, 0].axvline(x[0], label="True Intercept", lw=2, color="k")
-axes[0, 0].axvline(xnest[0], label="Intercept MAP", lw=2, color="C0", ls="--")
+axes[0, 0].axvline(x_map[0], label="Intercept MAP", lw=2, color="C0", ls="--")
 axes[0, 0].axvline(x[1], label="True Slope", lw=2, color="darkgray")
-axes[0, 0].axvline(xnest[1], label="Slope MAP", lw=2, color="C1", ls="--")
+axes[0, 0].axvline(x_map[1], label="Slope MAP", lw=2, color="C1", ls="--")
 axes[0, 1].axhline(x[0], label="True Intercept", lw=2, color="k")
-axes[0, 1].axhline(xnest[0], label="Intercept MAP", lw=2, color="C0", ls="--")
+axes[0, 1].axhline(x_map[0], label="Intercept MAP", lw=2, color="C0", ls="--")
 axes[0, 1].axhline(x[1], label="True Slope", lw=2, color="darkgray")
-axes[0, 1].axhline(xnest[1], label="Slope MAP", lw=2, color="C1", ls="--")
+axes[0, 1].axhline(x_map[1], label="Slope MAP", lw=2, color="C1", ls="--")
 axes[1, 0].axvline(sigma, label="True Sigma", lw=2, color="k")
-axes[1, 0].axvline(noise_est, label="Sigma MAP", lw=2, color="C0", ls="--")
+axes[1, 0].axvline(noise_map, label="Sigma MAP", lw=2, color="C0", ls="--")
 axes[1, 1].axhline(sigma, label="True Sigma", lw=2, color="k")
-axes[1, 1].axhline(noise_est, label="Sigma MAP", lw=2, color="C0", ls="--")
+axes[1, 1].axhline(noise_map, label="Sigma MAP", lw=2, color="C0", ls="--")
 for ax in axes.ravel():
     ax.legend()
 ax.get_figure().tight_layout()
@@ -203,20 +218,17 @@
     hdi_prob=0.95,
     ax=ax,
 )
-ax.plot(
-    np.array([t.min(), t.max()]),
-    np.array([t.min(), t.max()]) * x[1] + x[0],
-    "k",
-    lw=4,
-    label=rf"True: $x_0$ = {x[0]:.2f}, $x_1$ = {x[1]:.2f}",
-)
-ax.plot(
-    np.array([t.min(), t.max()]),
-    np.array([t.min(), t.max()]) * xnest[1] + xnest[0],
-    "--g",
-    lw=4,
-    label=rf"MAP Estimated: $x_0$ = {xnest[0]:.2f}, $x_1$ = {xnest[1]:.2f}",
-)
+for est, est_label, c in zip(
+    [x, x_mle, x_map], ["True", "MLE", "MAP"], ["k", "C0", "C1"]
+):
+    ax.plot(
+        np.array([t.min(), t.max()]),
+        np.array([t.min(), t.max()]) * est[1] + est[0],
+        color=c,
+        ls="--" if est_label == "MAP" else "-",
+        lw=4,
+        label=rf"{est_label}: $x_0$ = {est[0]:.2f}, $x_1$ = {est[1]:.2f}",
+    )
 ax.scatter(t, y, c="r", s=70)
 ax.scatter(t, yn, c="g", s=70)
 ax.legend()

From 417fc24bcca12c148787bb444d7b31150d67a69a Mon Sep 17 00:00:00 2001
From: Carlos da Costa <cako@users.noreply.github.com>
Date: Tue, 19 Nov 2024 17:13:24 -0800
Subject: [PATCH 34/34] fix: remove pytensor from yamls

---
 environment-dev-arm.yml | 1 -
 environment-dev.yml     | 1 -
 2 files changed, 2 deletions(-)

diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml
index acda12aa..0081a0ad 100755
--- a/environment-dev-arm.yml
+++ b/environment-dev-arm.yml
@@ -25,7 +25,6 @@ dependencies:
   - autopep8
   - isort
   - black
-  - pytensor
   - pip:
     - devito
     - dtcwt
diff --git a/environment-dev.yml b/environment-dev.yml
index d29d25b5..eb51c4dc 100755
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -26,7 +26,6 @@ dependencies:
   - autopep8
   - isort
   - black
-  - pytensor
   - pip:
     - devito
     - dtcwt