resolved merge conflicts

Author: Salman Naqvi

.github/workflows/build.yaml

@@ -4,25 +4,30 @@ on: [push, pull_request]
 
 jobs:
   build:
-
-    runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8]
+        platform: [ ubuntu-latest, macos-latest ]
+        python-version: ["3.8", "3.9", "3.10"]
 
+    runs-on: ${{ matrix.platform }}
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
+    - name: Get history and tags for SCM versioning to work
+      run: |
+        git fetch --prune --unshallow
+        git fetch --depth=1 origin +refs/tags/*:refs/tags/*
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
-        python -m pip install --upgrade pip
+        python -m pip install --upgrade pip setuptools
         pip install flake8 pytest
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
     - name: Install pyproximal
       run: |
+        python -m setuptools_scm
         pip install .
     - name: Test with pytest
       run: |

.readthedocs.yaml

@@ -21,5 +21,5 @@ sphinx:
 python:
   install:
     - requirements: requirements-dev.txt
-    - method: setuptools
+    - method: pip
       path: .

Makefile

@@ -31,7 +31,7 @@ dev-install_conda:
 
 tests:
 	make pythoncheck
-	$(PYTHON) setup.py test
+	pytest
 
 doc:
 	cd docs  && rm -rf source/api/generated && rm -rf source/gallery &&\

azure-pipelines.yml

@@ -64,7 +64,7 @@ jobs:
     displayName: 'Install prerequisites and library'
 
   - script: |
-      python setup.py test
+      pytest
     condition: succeededOrFailed()
     displayName: 'Run tests'
 
@@ -93,7 +93,7 @@ jobs:
     displayName: 'Install prerequisites and library'
 
   - script: |
-      python setup.py test
+      pytest
     condition: succeededOrFailed()
     displayName: 'Run tests'
 

docs/source/api/index.rst

@@ -92,6 +92,7 @@ Non-Convex
     Log1
     QuadraticEnvelopeCard
     QuadraticEnvelopeCardIndicator
+    RelaxedMumfordShah
     SCAD
 
 

docs/source/conf.py

@@ -2,7 +2,6 @@
 import sys
 import os
 import datetime
-import sphinx_rtd_theme
 import sphinx_gallery
 from sphinx_gallery.sorting import ExampleTitleSortKey
 from pyproximal import __version__

environment-dev.yml

@@ -20,7 +20,7 @@ dependencies:
     - bm3d
     - pytest-runner
     - setuptools_scm
-    - sphinx-rtd-theme
+    - pydata-sphinx-theme
     - sphinx-gallery
     - nbsphinx
     - image

pyproject.toml

@@ -0,0 +1,49 @@
+[build-system]
+requires = [
+    "setuptools >= 65",
+    "setuptools_scm[toml]",
+    "wheel",
+]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "pyproximal"
+description = "Python library implementing proximal operators to solve non-smooth, constrained convex problems with proximal algorithms"
+readme = "README.md"
+authors = [
+    {name = "Matteo Ravasi", email = "[email protected]"},
+]
+license = {file = "LICENSE.md"}
+keywords = ["algebra", "inverse problems", "proximal", "convex optimization", "large-scale optimization"]
+classifiers = [
+    "Development Status :: 5 - Production/Stable",
+    "Intended Audience :: Developers",
+    "Intended Audience :: Science/Research",
+    "Intended Audience :: Education",
+    "License :: OSI Approved :: GNU Lesser General Public License v3 (LGPLv3)",
+    "Natural Language :: English",
+    "Operating System :: OS Independent",
+    "Programming Language :: Python :: 3 :: Only",
+    "Programming Language :: Python :: 3.8",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10",
+    "Topic :: Scientific/Engineering : Mathematics",
+]
+dependencies = [
+    "numpy >= 1.15.0",
+    "scipy >= 1.8.0",
+    "pylops >= 2.0.0",
+]
+dynamic = ["version"]
+
+[project.optional-dependencies]
+advanced = [
+    "llvmlite",
+    "numba",
+]
+
+[tool.setuptools.packages.find]
+exclude = ["pytests"]
+
+[tool.setuptools_scm]
+version_file = "pyproximal/version.py"

pyproximal/optimization/primal.py

@@ -309,7 +309,7 @@ def GeneralizedProximalGradient(proxfs, proxgs, x0, tau=None,
     .. math::
 
         \mathbf{x} = \argmin_\mathbf{x} \sum_{i=1}^n f_i(\mathbf{x}) 
-        + \sum_{j=1}^m \tau_j g_j(\mathbf{x}),~~n,m \in \mathbb{N}^+
+        + \sum_{j=1}^m \epsilon_j g_j(\mathbf{x}),~~n,m \in \mathbb{N}^+
 
     where the :math:`f_i(\mathbf{x})` are smooth convex functions with a uniquely
     defined gradient and the :math:`g_j(\mathbf{x})` are any convex function that
@@ -330,7 +330,7 @@ def GeneralizedProximalGradient(proxfs, proxgs, x0, tau=None,
         backtracking is used to adaptively estimate the best tau at each
         iteration.
     epsg : :obj:`float` or :obj:`np.ndarray`, optional
-        Scaling factor of g function
+        Scaling factor(s) of ``g`` function(s)
     niter : :obj:`int`, optional
         Number of iterations of iterative scheme
     acceleration:  :obj:`str`, optional
@@ -353,11 +353,12 @@ def GeneralizedProximalGradient(proxfs, proxgs, x0, tau=None,
 
     .. math::
         \text{for } j=1,\cdots,n, \\
-        ~~~~\mathbf z_j^{k+1} = \mathbf z_j^{k} + \eta_k (prox_{\frac{\tau^k}{\omega_j} g_j}(2 \mathbf{x}^{k} - z_j^{k}) 
-        - \tau^k \sum_{i=1}^n \nabla f_i(\mathbf{x}^{k})) - \mathbf{x}^{k} \\
+        ~~~~\mathbf z_j^{k+1} = \mathbf z_j^{k} + \epsilon_j
+        \left[prox_{\frac{\tau^k}{\omega_j} g_j}\left(2 \mathbf{x}^{k} - \mathbf{z}_j^{k}
+        - \tau^k \sum_{i=1}^n \nabla f_i(\mathbf{x}^{k})\right) - \mathbf{x}^{k} \right] \\
         \mathbf{x}^{k+1} = \sum_{j=1}^n \omega_j f_j \\
         
-    where :math:`\sum_{j=1}^n \omega_j=1`. 
+    where :math:`\sum_{j=1}^n \omega_j=1`. In the current implementation :math:`\omega_j=1/n`.
     """
     # check if epgs is a vector
     if np.asarray(epsg).size == 1.:
@@ -394,23 +395,23 @@ def GeneralizedProximalGradient(proxfs, proxgs, x0, tau=None,
     for iiter in range(niter):
         xold = x.copy()
 
-        # proximal step
+        # gradient
         grad = np.zeros_like(x)
         for i, proxf in enumerate(proxfs):
             grad += proxf.grad(x)
 
-        sol = np.zeros_like(x)
+        # proximal step
+        x = np.zeros_like(x)
         for i, proxg in enumerate(proxgs):
-            tmp = 2 * y - zs[i] - tau * grad
-            tmp[:] = proxg.prox(tmp, epsg *tau *len(proxgs) )
-            zs[i] += (tmp - y)
-            sol += zs[i] / len(proxgs)
-        x[:] = sol.copy()
+            ztmp = 2 * y - zs[i] - tau * grad
+            ztmp = proxg.prox(ztmp, epsg * tau * len(proxgs))
+            zs[i] += (ztmp - y)
+            x += zs[i] / len(proxgs)
 
         # update y
         if acceleration == 'vandenberghe':
             omega = iiter / (iiter + 3)
-        elif acceleration== 'fista':
+        elif acceleration == 'fista':
             told = t
             t = (1. + np.sqrt(1. + 4. * t ** 2)) / 2.
             omega = ((told - 1.) / t)
@@ -782,7 +783,7 @@ def ADMML2(proxg, Op, b, A, x0, tau, niter=10, callback=None, show=False, **kwar
 
         if show:
             if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
-                pf, pg = np.linalg.norm(Op @ x - b), proxg(Ax)
+                pf, pg = 0.5 * np.linalg.norm(Op @ x - b) ** 2, proxg(Ax)
                 msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e' % \
                       (iiter + 1, x[0], pf, pg, pf + pg)
                 print(msg)

pyproximal/proximal/RelaxedMS.py

@@ -0,0 +1,99 @@
+import numpy as np
+
+from pyproximal.ProxOperator import _check_tau
+from pyproximal import ProxOperator
+from pyproximal.proximal.L1 import _current_sigma
+
+
+def _l2(x, alpha):
+    r"""Scaling operation.
+
+    Applies the proximal of ``alpha||y - x||_2^2`` which is essentially a scaling operation.
+
+    Parameters
+    ----------
+    x : :obj:`numpy.ndarray`
+        Vector
+    alpha : :obj:`float`
+        Scaling parameter
+
+    Returns
+    -------
+    y : :obj:`numpy.ndarray`
+        Scaled vector
+
+    """
+    y = 1 / (1 + 2 * alpha) * x
+    return y
+
+
+def _current_kappa(kappa, count):
+    if not callable(kappa):
+        return kappa
+    else:
+        return kappa(count)
+
+
+class RelaxedMumfordShah(ProxOperator):
+    r"""Relaxed Mumford-Shah norm proximal operator.
+
+    Proximal operator of the relaxed Mumford-Shah norm:
+    :math:`\text{rMS}(x) = \min (\alpha\Vert x\Vert_2^2, \kappa)`.
+
+    Parameters
+    ----------
+    sigma : :obj:`float` or :obj:`list` or :obj:`np.ndarray` or :obj:`func`, optional
+        Multiplicative coefficient of L2 norm that controls the smoothness of the solutuon.
+        This can be a constant number, a list of values (for multidimensional inputs, acting
+        on the second dimension) or a function that is called passing a counter which keeps
+        track of how many times the ``prox`` method has been invoked before and returns a
+        scalar (or a list of) ``sigma`` to be used.
+    kappa : :obj:`float` or :obj:`list` or :obj:`np.ndarray` or :obj:`func`, optional
+        Constant value in the rMS norm which essentially controls when the norm allows a jump. This can be a
+        constant number, a list of values (for multidimensional inputs, acting on the second dimension) or
+        a function that is called passing a counter which keeps track of how many
+        times the ``prox`` method has been invoked before and returns a scalar (or a list of)
+        ``kappa`` to be used.
+
+    Notes
+    -----
+    The :math:`rMS` proximal operator is defined as [1]_:
+
+    .. math::
+        \text{prox}_{\tau \text{rMS}}(x) =
+        \begin{cases}
+        \frac{1}{1+2\tau\alpha}x & \text{ if } & \vert x\vert \leq \sqrt{\frac{\kappa}{\alpha}(1 + 2\tau\alpha)} \\
+        \kappa & \text{ else }
+        \end{cases}.
+
+    .. [1] Strekalovskiy, E., and D. Cremers, 2014, Real-time minimization of the piecewise smooth
+            Mumford-Shah functional: European Conference on Computer Vision, 127–141.
+
+    """
+    def __init__(self, sigma=1., kappa=1.):
+        super().__init__(None, False)
+        self.sigma = sigma
+        self.kappa = kappa
+        self.count = 0
+
+    def __call__(self, x):
+        sigma = _current_sigma(self.sigma, self.count)
+        kappa = _current_sigma(self.kappa, self.count)
+        return np.minimum(sigma * np.linalg.norm(x) ** 2, kappa)
+
+    def _increment_count(func):
+        """Increment counter
+        """
+        def wrapped(self, *args, **kwargs):
+            self.count += 1
+            return func(self, *args, **kwargs)
+        return wrapped
+
+    @_increment_count
+    @_check_tau
+    def prox(self, x, tau):
+        sigma = _current_sigma(self.sigma, self.count)
+        kappa = _current_sigma(self.kappa, self.count)
+
+        x = np.where(np.abs(x) <= np.sqrt(kappa / sigma * (1 + 2 * tau * sigma)), _l2(x, tau * sigma), x)
+        return x