Repository: CUQI-DTU/CUQIpy

Issue:
Can we make linear model with 2D Image geometry more intuitive
Description:
------------
As suggested by @chaozg, it seems to be a recurring issue with users to find the right setting for 2D geometry when working with linear models (in particular setting the flag "visual_only")

There is a need to clear some confusion the user might encounter when preparing their forward model: If their forward model works on 2D images, they probably need `visual_only=False` (the default) when defining the domain geometry, and they probably need `visual_only=True` when the model works on vectors.

DoD
-----
- [x] A solution is found to make this setup more intuitive to the user.
    - It is decided to add a how-to to demonstrate/differentiate the use of Image2D
- [x] In a how-to file, demonstrate when `visual_only` is used and when it is not. 
    - The how-to should be concise and the examples should be intuitive.
    - The effect of using `visual_only` should be properly demonstrated and documented.
- [x] Guide the user to the how-to file by putting links at the following places
    - It is decided at stand-up that we only put a link at the docstring of Image2D
- [x] The solution is reviewed and approved

Discussion:
@amal-ghamdi: fyi @jakobsj, @chaozg and @nabriis , issue is created@chaozg: Hi @jakobsj, could you please help review the DoD here?@jakobsj: @chaozg I think you suggested in a recent meeting to add some more details on this issue and I think it is a good idea to do this, before revieving the DoD.  

So it would be great if you could add something on: what is the confusion, ie short examples of how the code is to be used, how a user might expect it instead and initial proposed suggestion(s) for how to make the behaviour more intuitive.@chaozg: Hi @jakobsj, could you please review the DoD here since it was updated a bit after our stand-up?@jakobsj: @chaozg A howto seems a useful solution to help guide the user. My only comment is how we make sure the user is made aware of this howto in their workflow when they face the challenges of understanding Image2D. Would it be suitable to refer the user to the howto in the documentation/docstring of Image2D geometry? Are there other places a user might look for help during that situation and where we could add a link or similar?  Please can you think of good way(s) and update the DoD and PR accordingly.@chaozg: Thank you @jakobsj: that's a very nice suggestion. I agree it would be useful to add a link to the how-to file at the docstrings of relevant CUQI classes. For now, CUQI model, Image2D, array and sample are in my mind.

Solution patch:
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 5654f5c11..0cf002c9c 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -22,6 +22,15 @@ jobs:
       - name: Sphinx build
         run: |
           sphinx-build docs ./public
+
+      # If PR-based build. Upload artifact of the documentation for download and review
+      - name: Upload artifact
+        if: ${{ github.event_name == 'pull_request' }}
+        uses: actions/upload-artifact@v3
+        with:
+          name: docs_build
+          path: ./public
+
       - name: Deploy
         uses: peaceiris/actions-gh-pages@v3
         if: ${{ github.event_name == 'push' && github.ref == 'refs/heads/main' }}
diff --git a/cuqi/experimental/mcmc/_langevin_algorithm.py b/cuqi/experimental/mcmc/_langevin_algorithm.py
index 9c1e9742c..5759a8548 100644
--- a/cuqi/experimental/mcmc/_langevin_algorithm.py
+++ b/cuqi/experimental/mcmc/_langevin_algorithm.py
@@ -3,7 +3,7 @@
 from cuqi.experimental.mcmc import Sampler
 from cuqi.implicitprior import RestorationPrior, MoreauYoshidaPrior
 from cuqi.array import CUQIarray
-from copy import deepcopy
+from copy import copy
 
 class ULA(Sampler): # Refactor to Proposal-based sampler?
     """Unadjusted Langevin algorithm (ULA) (Roberts and Tweedie, 1996)
@@ -316,11 +316,14 @@ def target(self, value):
 
     def _create_smoothed_target(self, value):
         """ Create a smoothed target using a Moreau-Yoshida envelope. """
-        copied_value = deepcopy(value)
+        copied_value = copy(value)
         if isinstance(copied_value.prior, RestorationPrior):
+            # Acceess the prior name
+            name = value.prior.name
             copied_value.prior = MoreauYoshidaPrior(
                 copied_value.prior,
-                self.smoothing_strength)
+                self.smoothing_strength,
+                name=name)
         return copied_value
 
     def validate_target(self):
diff --git a/cuqi/geometry/_geometry.py b/cuqi/geometry/_geometry.py
index 6361d357f..742ccf6ed 100644
--- a/cuqi/geometry/_geometry.py
+++ b/cuqi/geometry/_geometry.py
@@ -522,6 +522,9 @@ class Image2D(Geometry):
     Plotting is handled via matplotlib.pyplot.imshow.
     Colormap is defaulted to grayscale.
 
+    A How-To guide on the use of Image2D as the domain/range geometry of a CUQI 
+    :class:`Model` is available `here <https://cuqi-dtu.github.io/CUQIpy/user/_auto_howtos/myula.html>`_.
+
     Parameters
     -----------
     im_shape : tuple
diff --git a/demos/howtos/Image2D.py b/demos/howtos/Image2D.py
new file mode 100644
index 000000000..0a1ac893a
--- /dev/null
+++ b/demos/howtos/Image2D.py
@@ -0,0 +1,70 @@
+"""
+The use of Image2D
+==================
+
+This script demonstrates the use of :class:`~cuqi.geometry.Image2D` for domain and range geometries in CUQIpy forward models. The settings of this geometry, in particular the so-called `visual_only` flag, depend on how the forward model input is operated on in the forward model. `visual_only` flag tells the geometry to interpret the variable as a 2D image only when visualizing the variable but keep it as a vector when being operated on in the forward model. Generally, if the forward operator expects an image input, set `visual_only` to `False`, which is the default in CUQIpy. Otherwise, if it expects a vector input, set `visual_only` to `True` in the domain geometry. In the same way, the range geometry `visual_only` setting should be consistent with your forward operator output.
+"""
+
+# %%
+# Import necessary libraries
+import cuqi
+import numpy as np
+import matplotlib.pyplot as plt
+
+# %%
+# Define dimensions
+dim_x, dim_y = 2, 4
+
+# %%
+# Example illustrating the default behavior: `visual_only=False`
+# With `visual_only=False`, the underlying structure of a CUQI array or sample will be reshaped to the
+# shape of the specified geometry, i.e., an image of size (dim_x, dim_y), before being passed to the model.
+# Here we intend to define a forward operator that downsamples the input image by a factor of 2. It's easier
+# to define such a forward operator in terms of an image. We do not need to explicitly set `visual_only` to 
+#`False` because it is the default behavior.
+
+def forward_func_a(image):
+    assert image.shape == (dim_x, dim_y)
+    return image[::2, ::2]
+
+model_A = cuqi.model.LinearModel(forward_func_a,
+                                 domain_geometry=cuqi.geometry.Image2D((dim_x, dim_y)),
+                                 range_geometry=cuqi.geometry.Image2D((dim_x // 2, dim_y // 2)))
+
+print(model_A)
+
+x_sample = cuqi.array.CUQIarray(np.linspace(0, 1, dim_x * dim_y),
+                                geometry=cuqi.geometry.Image2D((dim_x, dim_y)))
+
+# Plot the original sample
+plt.figure()
+x_sample.plot(cmap='Greens')
+plt.title("input image x")
+
+# Apply the model and plot the result
+plt.figure()
+y_sample = model_A @ x_sample
+y_sample.plot(cmap='Greens')
+plt.title("output y=model_A(x)")
+
+# %%
+# Example with `visual_only=True`
+# With `visual_only=True`, the underlying structure of a cuqi array or sample will not be changed before
+# being passed to the model. This is useful when the model expects a vector input. Here we define a 
+# forward operator that reverses the input vector, which can be easily defined in terms of a vector.
+def forward_func_b(x):
+    assert x.shape == (dim_x * dim_y, )
+    return x[::-1]
+
+model_B = cuqi.model.LinearModel(forward_func_b,
+                                 domain_geometry=cuqi.geometry.Image2D((dim_x, dim_y), visual_only=True),
+                                 range_geometry=cuqi.geometry.Image2D((dim_x, dim_y), visual_only=True))
+
+# Note we explicitely set the `visual_only` flag to be `True` here.
+print(model_B)
+
+# Apply the model and plot the result
+plt.figure()
+y_sample = model_B @ x_sample
+y_sample.plot(cmap='Greens')
+plt.title("output y=model_B(x)")
\ No newline at end of file


Test patch:
diff --git a/tests/zexperimental/test_mcmc.py b/tests/zexperimental/test_mcmc.py
index ad31b2387..177962dcd 100644
--- a/tests/zexperimental/test_mcmc.py
+++ b/tests/zexperimental/test_mcmc.py
@@ -330,12 +330,15 @@ def create_myula_target(dim=16):
     """Create a target for MYULA."""
     def func(x, restoration_strength=1):
         return x, True
-    likelihood = cuqi.testproblem.Deconvolution1D(
+    y = cuqi.testproblem.Deconvolution1D(
         dim=dim).posterior.likelihood 
-    restoration_prior = cuqi.implicitprior.RestorationPrior(
-        func, geometry=likelihood.model.domain_geometry)
+    x = cuqi.implicitprior.RestorationPrior(
+        func, geometry=y.model.domain_geometry)
+    # access name property to ensure RestorationPrior
+    # name is set 
+    x.name
     posterior = cuqi.distribution.Posterior(
-        likelihood, restoration_prior)
+        y, x)
     return posterior
 
 def create_myula_smoothed_target(dim=16):
@@ -368,6 +371,10 @@ def test_myula():
     assert samples_ula.Ns == 10
     assert samples_myula.Ns == 10
     assert np.allclose(samples_myula.samples, samples_ula.samples)
+    # Assert the name of the MoreauYoshidaPrior created by MYULA is
+    # the same name as the corresponding RestorationPrior (both are named x)
+    assert myula._smoothed_target.prior.name == 'x'
+    assert myula.target.prior.name == 'x'
 
 def test_myula_object_creation_fails_with_target_without_restore_method():
     """ Test that MYULA object creation fails with target that does not