diff --git a/CHANGES.rst b/CHANGES.rst index 2e1b4b5a..45e979d6 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -6,7 +6,8 @@ SPORCO Release Notes Version 0.2.2 (unreleased) ---------------------------- -• No changes yet +• Removed deprecated sporco.fista modules +• Removed deprecation warning redirects for functions renamed in 0.2.0 Version 0.2.1 (2022-02-17) diff --git a/sporco/__init__.py b/sporco/__init__.py index 003f2724..8d9c2529 100644 --- a/sporco/__init__.py +++ b/sporco/__init__.py @@ -1,6 +1,6 @@ from __future__ import absolute_import -__version__ = '0.2.2a1' +__version__ = '0.2.2a2' # This is a temporary solution to the circular imports resulting from the diff --git a/sporco/array.py b/sporco/array.py index 6717ed84..6d1c662e 100644 --- a/sporco/array.py +++ b/sporco/array.py @@ -15,8 +15,6 @@ import numpy as np -from sporco._util import renamed_function - __author__ = """Brendt Wohlberg """ @@ -94,7 +92,6 @@ def transpose_ntpl_list(lst): -@renamed_function(depname='zpad', depmod='sporco.linalg') def zpad(x, pd, ax): """Zero-pad array `x` with `pd = (leading, trailing)` zeros on axis `ax`. @@ -119,7 +116,6 @@ def zpad(x, pd, ax): -@renamed_function(depname='zdivide', depmod='sporco.linalg') def zdivide(x, y): """Return `x`/`y`, with 0 instead of NaN where `y` is 0. @@ -141,7 +137,6 @@ def zdivide(x, y): -@renamed_function(depname='promote16', depmod='sporco.linalg') def promote16(u, fn=None, *args, **kwargs): r"""Promote ``np.float16`` arguments to ``np.float32`` dtype. @@ -191,7 +186,6 @@ def promote16(u, fn=None, *args, **kwargs): -@renamed_function(depname='atleast_nd', depmod='sporco.linalg') def atleast_nd(n, u): """Append axes to an array so that it is ``n`` dimensional. @@ -219,7 +213,6 @@ def atleast_nd(n, u): -@renamed_function(depname='split', depmod='sporco.linalg') def split(u, axis=0): """Split an array into a list of arrays on the specified axis. diff --git a/sporco/fft.py b/sporco/fft.py index 4b5fb8cc..ab1ed280 100644 --- a/sporco/fft.py +++ b/sporco/fft.py @@ -25,7 +25,6 @@ else: have_pyfftw = True -from sporco._util import renamed_function __author__ = """Brendt Wohlberg """ @@ -60,7 +59,6 @@ def is_complex_dtype(dtype): -@renamed_function(depname='complex_dtype', depmod='sporco.linalg') def complex_dtype(dtype): """Construct the corresponding complex dtype for a given real dtype. @@ -105,7 +103,6 @@ def real_dtype(dtype): -@renamed_function(depname='pyfftw_byte_aligned', depmod='sporco.linalg') def byte_aligned(array, dtype=None, n=None): """Construct a byte-aligned array for FFTs. @@ -131,7 +128,6 @@ def byte_aligned(array, dtype=None, n=None): -@renamed_function(depname='pyfftw_empty_aligned', depmod='sporco.linalg') def empty_aligned(shape, dtype, order='C', n=None): """Construct an empty byte-aligned array for FFTs. @@ -160,7 +156,6 @@ def empty_aligned(shape, dtype, order='C', n=None): -@renamed_function(depname='pyfftw_rfftn_empty_aligned', depmod='sporco.linalg') def rfftn_empty_aligned(shape, axes, dtype, order='C', n=None): """Construct an empty byte-aligned array for real FFTs. @@ -201,7 +196,6 @@ def rfftn_empty_aligned(shape, axes, dtype, order='C', n=None): -@renamed_function(depname='fftn', depmod='sporco.linalg') def fftn(a, s=None, axes=None): """Multi-dimensional discrete Fourier transform. @@ -231,7 +225,6 @@ def fftn(a, s=None, axes=None): -@renamed_function(depname='ifftn', depmod='sporco.linalg') def ifftn(a, s=None, axes=None): """Multi-dimensional inverse discrete Fourier transform. @@ -261,7 +254,6 @@ def ifftn(a, s=None, axes=None): -@renamed_function(depname='rfftn', depmod='sporco.linalg') def rfftn(a, s=None, axes=None): """Multi-dimensional discrete Fourier transform for real input. @@ -291,7 +283,6 @@ def rfftn(a, s=None, axes=None): -@renamed_function(depname='irfftn', depmod='sporco.linalg') def irfftn(a, s, axes=None): """Multi-dimensional inverse discrete Fourier transform for real input. @@ -324,7 +315,6 @@ def irfftn(a, s, axes=None): -@renamed_function(depname='dctii', depmod='sporco.linalg') def dctii(x, axes=None): """Multi-dimensional DCT-II. @@ -354,7 +344,6 @@ def dctii(x, axes=None): -@renamed_function(depname='idctii', depmod='sporco.linalg') def idctii(x, axes=None): """Multi-dimensional inverse DCT-II. @@ -384,7 +373,6 @@ def idctii(x, axes=None): -@renamed_function(depname='fftconv', depmod='sporco.linalg') def fftconv(a, b, axes=None, origin=None): """Multi-dimensional convolution via the Discrete Fourier Transform. @@ -430,7 +418,6 @@ def fftconv(a, b, axes=None, origin=None): -@renamed_function(depname='fl2norm2', depmod='sporco.linalg') def fl2norm2(xf, axis=(0, 1)): r"""Compute the squared :math:`\ell_2` norm in the DFT domain. @@ -459,7 +446,6 @@ def fl2norm2(xf, axis=(0, 1)): -@renamed_function(depname='rfl2norm2', depmod='sporco.linalg') def rfl2norm2(xf, xs, axis=(0, 1)): r"""Compute the squared :math:`\ell_2` norm in the real DFT domain. diff --git a/sporco/fista/__init__.py b/sporco/fista/__init__.py deleted file mode 100644 index 4fbf6dc6..00000000 --- a/sporco/fista/__init__.py +++ /dev/null @@ -1,9 +0,0 @@ -from __future__ import absolute_import -import warnings - -import sporco.fista.fista - -warnings.simplefilter('always', DeprecationWarning) -warnings.warn('Module sporco.fista is deprecated; please use sporco.pgm', - DeprecationWarning) -warnings.simplefilter('default', DeprecationWarning) diff --git a/sporco/fista/bpdn.py b/sporco/fista/bpdn.py deleted file mode 100644 index c27b34f9..00000000 --- a/sporco/fista/bpdn.py +++ /dev/null @@ -1,229 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright (C) 2018-2019 by Brendt Wohlberg -# Cristina Garcia-Cardona -# All rights reserved. BSD 3-clause License. -# This file is part of the SPORCO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -"""Classes for FISTA algorithm for the BPDN problem""" - -from __future__ import division, absolute_import, print_function - -import copy -import numpy as np - -from sporco.util import u -from sporco.prox import prox_l1 -from sporco.fista import fista - - -__author__ = """Cristina Garcia-Cardona """ - - - -class BPDN(fista.FISTA): - r""" - Class for FISTA algorithm for the Basis Pursuit DeNoising (BPDN) - :cite:`chen-1998-atomic` problem. - - | - - .. inheritance-diagram:: BPDN - :parts: 2 - - | - - The problem form is - - .. math:: - \mathrm{argmin}_\mathbf{x} \; (1/2) \| D \mathbf{x} - \mathbf{s} - \|_2^2 + \lambda \| \mathbf{x} \|_1 - - where :math:`\mathbf{s}` is the input vector/matrix, :math:`D` is - the dictionary, and :math:`\mathbf{x}` is the sparse representation. - - After termination of the :meth:`solve` method, attribute - :attr:`itstat` is a list of tuples representing statistics of each - iteration. The fields of the named tuple ``IterationStats`` are: - - ``Iter`` : Iteration number - - ``ObjFun`` : Objective function value - - ``DFid`` : Value of data fidelity term :math:`(1/2) \| D - \mathbf{x} - \mathbf{s} \|_2^2` - - ``RegL1`` : Value of regularisation term :math:`\lambda \| - \mathbf{x} \|_1` - - ``Rsdl`` : Residual - - ``L`` : Inverse of gradient step parameter - - ``Time`` : Cumulative run time - """ - - - class Options(fista.FISTA.Options): - r"""BPDN algorithm options - - Options include all of those defined in - :class:`.fista.FISTA.Options`, together with - additional options: - - ``L1Weight`` : An array of weights for the :math:`\ell_1` - norm. The array shape must be such that the array is - compatible for multiplication with the X/Y variables. If this - option is defined, the regularization term is :math:`\lambda - \| \mathbf{w} \odot \mathbf{x} \|_1` where - :math:`\mathbf{w}` denotes the weighting array. - - """ - - defaults = copy.deepcopy(fista.FISTADFT.Options.defaults) - defaults.update({'L1Weight': 1.0}) - defaults.update({'L': 500.0}) - - - def __init__(self, opt=None): - """ - Parameters - ---------- - opt : dict or None, optional (default None) - BPDN algorithm options - """ - - if opt is None: - opt = {} - fista.FISTA.Options.__init__(self, opt) - - - def __setitem__(self, key, value): - """Set options.""" - - fista.FISTA.Options.__setitem__(self, key, value) - - - itstat_fields_objfn = ('ObjFun', 'DFid', 'RegL1') - hdrtxt_objfn = ('Fnc', 'DFid', u('Regℓ1')) - hdrval_objfun = {'Fnc': 'ObjFun', 'DFid': 'DFid', u('Regℓ1'): 'RegL1'} - - - - def __init__(self, D, S, lmbda=None, opt=None): - """ - Parameters - ---------- - D : array_like - Dictionary array (2d) - S : array_like - Signal array (1d or 2d) - lmbda : float - Regularisation parameter - opt : :class:`BPDN.Options` object - Algorithm options - """ - - # Set default options if none specified - if opt is None: - opt = BPDN.Options() - - # Set dtype attribute based on S.dtype and opt['DataType'] - self.set_dtype(opt, S.dtype) - - # Set default lambda value if not specified - if lmbda is None: - DTS = D.T.dot(S) - lmbda = 0.1 * abs(DTS).max() - - # Set l1 term scaling and weight array - self.lmbda = self.dtype.type(lmbda) - self.wl1 = np.asarray(opt['L1Weight'], dtype=self.dtype) - - # Call parent class __init__ - Nc = D.shape[1] - Nm = S.shape[1] - - xshape = (Nc, Nm) - super(BPDN, self).__init__(xshape, S.dtype, opt) - - self.S = np.asarray(S, dtype=self.dtype) - - self.store_prev() - self.Y = self.X.copy() - self.Yprv = self.Y.copy() + 1e5 - - self.setdict(D) - - - - def setdict(self, D): - """Set dictionary array.""" - - self.D = np.asarray(D, dtype=self.dtype) - - - - def getcoef(self): - """Get final coefficient array.""" - - return self.X - - - - def eval_grad(self): - """Compute gradient in spatial domain for variable Y.""" - - # Compute D^T(D Y - S) - return self.D.T.dot(self.D.dot(self.Y) - self.S) - - - - def eval_proxop(self, V): - """Compute proximal operator of :math:`g`.""" - - return np.asarray(prox_l1(V, (self.lmbda / self.L) * self.wl1), - dtype=self.dtype) - - - - def eval_objfn(self): - """Compute components of objective function as well as total - contribution to objective function. - """ - - dfd = self.obfn_f() - reg = self.obfn_reg() - obj = dfd + reg[0] - return (obj, dfd) + reg[1:] - - - - def obfn_reg(self): - """Compute regularisation term and contribution to objective - function. - """ - - rl1 = np.linalg.norm((self.wl1 * self.X).ravel(), 1) - return (self.lmbda * rl1, rl1) - - - - def obfn_f(self, X=None): - r"""Compute data fidelity term :math:`(1/2) \| D \mathbf{x} - - \mathbf{s} \|_2^2`. - """ - if X is None: - X = self.X - - return 0.5 * np.linalg.norm((self.D.dot(X) - self.S).ravel())**2 - - - - def reconstruct(self, X=None): - """Reconstruct representation.""" - - if X is None: - X = self.X - return self.D.dot(self.X) diff --git a/sporco/fista/cbpdn.py b/sporco/fista/cbpdn.py deleted file mode 100644 index 7acd73e0..00000000 --- a/sporco/fista/cbpdn.py +++ /dev/null @@ -1,474 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright (C) 2016-2019 by Brendt Wohlberg -# Cristina Garcia-Cardona -# All rights reserved. BSD 3-clause License. -# This file is part of the SPORCO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -"""Classes for FISTA algorithm for the Convolutional BPDN problem""" - -from __future__ import division, absolute_import, print_function - -import copy -import numpy as np - -from sporco.util import u -from sporco.cnvrep import CSC_ConvRepIndexing, mskWshape -from sporco.linalg import inner -from sporco.fft import (rfftn, irfftn, empty_aligned, rfftn_empty_aligned, - rfl2norm2) -from sporco.prox import prox_l1 -from sporco.fista import fista - - -__author__ = """Cristina Garcia-Cardona """ - - - -class ConvBPDN(fista.FISTADFT): - r""" - Base class for FISTA algorithm for the Convolutional BPDN (CBPDN) - :cite:`garcia-2018-convolutional1` problem. - - | - - .. inheritance-diagram:: ConvBPDN - :parts: 2 - - | - - The generic problem form is - - .. math:: - \mathrm{argmin}_\mathbf{x} \; - f( \{ \mathbf{x}_m \} ) + \lambda g( \{ \mathbf{x}_m \} ) - - where :math:`f = (1/2) \left\| \sum_m \mathbf{d}_m * \mathbf{x}_m - - \mathbf{s} \right\|_2^2`, and :math:`g(\cdot)` is a penalty - term or the indicator function of a constraint; with input - image :math:`\mathbf{s}`, dictionary filters :math:`\mathbf{d}_m`, - and coefficient maps :math:`\mathbf{x}_m`. It is solved via the - FISTA formulation - - Proximal step - - .. math:: - \mathbf{x}_k = \mathrm{prox}_{t_k}(g) (\mathbf{y}_k - 1/L \nabla - f(\mathbf{y}_k) ) \;\;. - - Combination step - - .. math:: - \mathbf{y}_{k+1} = \mathbf{x}_k + \left( \frac{t_k - 1}{t_{k+1}} - \right) (\mathbf{x}_k - \mathbf{x}_{k-1}) \;\;, - - with :math:`t_{k+1} = \frac{1 + \sqrt{1 + 4 t_k^2}}{2}`. - - - After termination of the :meth:`solve` method, attribute - :attr:`itstat` is a list of tuples representing statistics of each - iteration. The fields of the named tuple ``IterationStats`` are: - - ``Iter`` : Iteration number - - ``ObjFun`` : Objective function value - - ``DFid`` : Value of data fidelity term :math:`(1/2) \| \sum_m - \mathbf{d}_m * \mathbf{x}_m - \mathbf{s} \|_2^2` - - ``RegL1`` : Value of regularisation term :math:`\sum_m \| - \mathbf{x}_m \|_1` - - ``Rsdl`` : Residual - - ``L`` : Inverse of gradient step parameter - - ``Time`` : Cumulative run time - """ - - - class Options(fista.FISTADFT.Options): - r"""ConvBPDN algorithm options - - Options include all of those defined in - :class:`.fista.FISTADFT.Options`, together with - additional options: - - ``NonNegCoef`` : Flag indicating whether to force solution to - be non-negative. - - ``NoBndryCross`` : Flag indicating whether all solution - coefficients corresponding to filters crossing the image - boundary should be forced to zero. - - ``L1Weight`` : An array of weights for the :math:`\ell_1` - norm. The array shape must be such that the array is - compatible for multiplication with the X/Y variables. If this - option is defined, the regularization term is :math:`\lambda - \sum_m \| \mathbf{w}_m \odot \mathbf{x}_m \|_1` where - :math:`\mathbf{w}_m` denotes slices of the weighting array on - the filter index axis. - - """ - - defaults = copy.deepcopy(fista.FISTADFT.Options.defaults) - defaults.update({'NonNegCoef': False, 'NoBndryCross': False}) - defaults.update({'L1Weight': 1.0}) - defaults.update({'L': 500.0}) - - - def __init__(self, opt=None): - """ - Parameters - ---------- - opt : dict or None, optional (default None) - ConvBPDN algorithm options - """ - - if opt is None: - opt = {} - fista.FISTADFT.Options.__init__(self, opt) - - - - def __setitem__(self, key, value): - """Set options.""" - - fista.FISTADFT.Options.__setitem__(self, key, value) - - - itstat_fields_objfn = ('ObjFun', 'DFid', 'RegL1') - hdrtxt_objfn = ('Fnc', 'DFid', u('Regℓ1')) - hdrval_objfun = {'Fnc': 'ObjFun', 'DFid': 'DFid', u('Regℓ1'): 'RegL1'} - - - - def __init__(self, D, S, lmbda=None, opt=None, dimK=None, dimN=2): - """ - This class supports an arbitrary number of spatial dimensions, - `dimN`, with a default of 2. The input dictionary `D` is either - `dimN` + 1 dimensional, in which case each spatial component - (image in the default case) is assumed to consist of a single - channel, or `dimN` + 2 dimensional, in which case the final - dimension is assumed to contain the channels (e.g. colour - channels in the case of images). The input signal set `S` is - either `dimN` dimensional (no channels, only one signal), - `dimN` + 1 dimensional (either multiple channels or multiple - signals), or `dimN` + 2 dimensional (multiple channels and - multiple signals). Determination of problem dimensions is - handled by :class:`.cnvrep.CSC_ConvRepIndexing`. - - - Parameters - ---------- - D : array_like - Dictionary array - S : array_like - Signal array - lmbda : float - Regularisation parameter - opt : :class:`ConvBPDN.Options` object - Algorithm options - dimK : 0, 1, or None, optional (default None) - Number of dimensions in input signal corresponding to multiple - independent signals - dimN : int, optional (default 2) - Number of spatial/temporal dimensions - """ - - # Set default options if none specified - if opt is None: - opt = ConvBPDN.Options() - - # Infer problem dimensions and set relevant attributes of self - if not hasattr(self, 'cri'): - self.cri = CSC_ConvRepIndexing(D, S, dimK=dimK, dimN=dimN) - - # Set dtype attribute based on S.dtype and opt['DataType'] - self.set_dtype(opt, S.dtype) - - # Set default lambda value if not specified - if lmbda is None: - cri = CSC_ConvRepIndexing(D, S, dimK=dimK, dimN=dimN) - Df = rfftn(D.reshape(cri.shpD), cri.Nv, axes=cri.axisN) - Sf = rfftn(S.reshape(cri.shpS), axes=cri.axisN) - b = np.conj(Df) * Sf - lmbda = 0.1 * abs(b).max() - - # Set l1 term scaling and weight array - self.lmbda = self.dtype.type(lmbda) - self.wl1 = np.asarray(opt['L1Weight'], dtype=self.dtype) - - # Call parent class __init__ - self.Xf = None - xshape = self.cri.shpX - super(ConvBPDN, self).__init__(xshape, S.dtype, opt) - - # Reshape D and S to standard layout - self.D = np.asarray(D.reshape(self.cri.shpD), dtype=self.dtype) - self.S = np.asarray(S.reshape(self.cri.shpS), dtype=self.dtype) - - # Compute signal in DFT domain - self.Sf = rfftn(self.S, None, self.cri.axisN) - - # Create byte aligned arrays for FFT calls - self.Y = self.X.copy() - self.X = empty_aligned(self.Y.shape, dtype=self.dtype) - self.X[:] = self.Y - - # Initialise auxiliary variable Vf: Create byte aligned arrays - # for FFT calls - self.Vf = rfftn_empty_aligned(self.X.shape, self.cri.axisN, - self.dtype) - - - self.Xf = rfftn(self.X, None, self.cri.axisN) - self.Yf = self.Xf.copy() - self.store_prev() - self.Yfprv = self.Yf.copy() + 1e5 - - self.setdict() - - # Initialization needed for back tracking (if selected) - self.postinitialization_backtracking_DFT() - - - - def setdict(self, D=None): - """Set dictionary array.""" - - if D is not None: - self.D = np.asarray(D, dtype=self.dtype) - self.Df = rfftn(self.D, self.cri.Nv, self.cri.axisN) - - - - def getcoef(self): - """Get final coefficient array.""" - - return self.X - - - - def eval_grad(self): - """Compute gradient in Fourier domain.""" - - # Compute D X - S - Ryf = self.eval_Rf(self.Yf) - # Compute D^H Ryf - gradf = np.conj(self.Df) * Ryf - - # Multiple channel signal, multiple channel dictionary - if self.cri.Cd > 1: - gradf = np.sum(gradf, axis=self.cri.axisC, keepdims=True) - - return gradf - - - - def eval_Rf(self, Vf): - """Evaluate smooth term in Vf.""" - - return inner(self.Df, Vf, axis=self.cri.axisM) - self.Sf - - - - def eval_proxop(self, V): - """Compute proximal operator of :math:`g`.""" - - return prox_l1(V, (self.lmbda / self.L) * self.wl1) - - - - def rsdl(self): - """Compute fixed point residual in Fourier domain.""" - - diff = self.Xf - self.Yfprv - return rfl2norm2(diff, self.X.shape, axis=self.cri.axisN) - - - - def eval_objfn(self): - """Compute components of objective function as well as total - contribution to objective function. - """ - - dfd = self.obfn_dfd() - reg = self.obfn_reg() - obj = dfd + reg[0] - return (obj, dfd) + reg[1:] - - - - def obfn_dfd(self): - r"""Compute data fidelity term :math:`(1/2) \| \sum_m - \mathbf{d}_m * \mathbf{x}_m - \mathbf{s} \|_2^2`. - This function takes into account the unnormalised DFT scaling, - i.e. given that the variables are the DFT of multi-dimensional - arrays computed via :func:`.rfftn`, this returns the data - fidelity term in the original (spatial) domain. - """ - - Ef = self.eval_Rf(self.Xf) - return rfl2norm2(Ef, self.S.shape, axis=self.cri.axisN) / 2.0 - - - - def obfn_reg(self): - """Compute regularisation term and contribution to objective - function. - """ - - rl1 = np.linalg.norm((self.wl1 * self.X).ravel(), 1) - return (self.lmbda * rl1, rl1) - - - - def obfn_f(self, Xf=None): - r"""Compute data fidelity term :math:`(1/2) \| \sum_m - \mathbf{d}_m * \mathbf{x}_m - \mathbf{s} \|_2^2` - This is used for backtracking. Since the backtracking is - computed in the DFT, it is important to preserve the - DFT scaling. - """ - - if Xf is None: - Xf = self.Xf - - Rf = self.eval_Rf(Xf) - return 0.5 * np.linalg.norm(Rf.flatten(), 2)**2 - - - - def reconstruct(self, X=None): - """Reconstruct representation.""" - - if X is None: - X = self.X - Xf = rfftn(X, None, self.cri.axisN) - Sf = np.sum(self.Df * Xf, axis=self.cri.axisM) - return irfftn(Sf, self.cri.Nv, self.cri.axisN) - - - - - -class ConvBPDNMask(ConvBPDN): - r""" - FISTA algorithm for Convolutional BPDN with a spatial mask. - - | - - .. inheritance-diagram:: ConvBPDNMask - :parts: 2 - - | - - Solve the optimisation problem - - .. math:: - \mathrm{argmin}_\mathbf{x} \; - (1/2) \left\| W \left(\sum_m \mathbf{d}_m * \mathbf{x}_m - - \mathbf{s}\right) \right\|_2^2 + \lambda \sum_m - \| \mathbf{x}_m \|_1 \;\;, - - where :math:`W` is a mask array. - - See :class:`ConvBPDN` for interface details. - """ - - - def __init__(self, D, S, lmbda, W=None, opt=None, dimK=None, dimN=2): - """ - - | - - - Parameters - ---------- - D : array_like - Dictionary matrix - S : array_like - Signal vector or matrix - lmbda : float - Regularisation parameter - W : array_like - Mask array. The array shape must be such that the array is - compatible for multiplication with input array S (see - :func:`.cnvrep.mskWshape` for more details). - opt : :class:`ConvBPDNMask.Options` object - Algorithm options - dimK : 0, 1, or None, optional (default None) - Number of dimensions in input signal corresponding to multiple - independent signals - dimN : int, optional (default 2) - Number of spatial dimensions - """ - - super(ConvBPDNMask, self).__init__(D, S, lmbda, opt, dimK=dimK, - dimN=dimN) - - if W is None: - W = np.array([1.0], dtype=self.dtype) - self.W = np.asarray(W.reshape(mskWshape(W, self.cri)), - dtype=self.dtype) - - # Create byte aligned arrays for FFT calls - self.WRy = empty_aligned(self.S.shape, dtype=self.dtype) - self.Ryf = rfftn_empty_aligned(self.S.shape, self.cri.axisN, - self.dtype) - - - - def eval_grad(self): - """Compute gradient in Fourier domain.""" - - # Compute D X - S - self.Ryf[:] = self.eval_Rf(self.Yf) - - # Map to spatial domain to multiply by mask - Ry = irfftn(self.Ryf, self.cri.Nv, self.cri.axisN) - # Multiply by mask - self.WRy[:] = (self.W**2) * Ry - # Map back to frequency domain - WRyf = rfftn(self.WRy, self.cri.Nv, self.cri.axisN) - - gradf = np.conj(self.Df) * WRyf - - # Multiple channel signal, multiple channel dictionary - if self.cri.Cd > 1: - gradf = np.sum(gradf, axis=self.cri.axisC, keepdims=True) - - return gradf - - - - def obfn_dfd(self): - r"""Compute data fidelity term :math:`(1/2) \| W (\sum_m - \mathbf{d}_m * \mathbf{x}_{m} - \mathbf{s}) \|_2^2` - """ - - Ef = self.eval_Rf(self.Xf) - E = irfftn(Ef, self.cri.Nv, self.cri.axisN) - - return (np.linalg.norm(self.W * E)**2) / 2.0 - - - - def obfn_f(self, Xf=None): - r"""Compute data fidelity term :math:`(1/2) \| W (\sum_m - \mathbf{d}_m * \mathbf{x}_{m} - \mathbf{s}) \|_2^2`. - This is used for backtracking. Since the backtracking is - computed in the DFT, it is important to preserve the - DFT scaling. - """ - - if Xf is None: - Xf = self.Xf - - Rf = self.eval_Rf(Xf) - R = irfftn(Rf, self.cri.Nv, self.cri.axisN) - WRf = rfftn(self.W * R, self.cri.Nv, self.cri.axisN) - - return 0.5 * np.linalg.norm(WRf.flatten(), 2)**2 diff --git a/sporco/fista/ccmod.py b/sporco/fista/ccmod.py deleted file mode 100644 index b07da802..00000000 --- a/sporco/fista/ccmod.py +++ /dev/null @@ -1,564 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright (C) 2016-2019 by Brendt Wohlberg -# Cristina Garcia-Cardona -# All rights reserved. BSD 3-clause License. -# This file is part of the SPORCO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -"""FISTA algorithms for the CCMOD problem""" - -from __future__ import division, absolute_import - -import copy -import numpy as np - -from sporco.fista import fista -from sporco.array import atleast_nd -from sporco.fft import (rfftn, irfftn, empty_aligned, rfftn_empty_aligned, - rfl2norm2) -from sporco.linalg import inner -from sporco.cnvrep import CDU_ConvRepIndexing, getPcn, bcrop - - -__author__ = """Cristina Garcia-Cardona """ - - - -class ConvCnstrMOD(fista.FISTADFT): - r""" - Base class for FISTA algorithm for Convolutional Constrained MOD - problem :cite:`garcia-2018-convolutional1`. - - | - - .. inheritance-diagram:: ConvCnstrMOD - :parts: 2 - - | - - Solve the optimisation problem - - .. math:: - \mathrm{argmin}_\mathbf{d} \; - (1/2) \sum_k \left\| \sum_m \mathbf{d}_m * \mathbf{x}_{k,m} - - \mathbf{s}_k \right\|_2^2 \quad \text{such that} \quad - \mathbf{d}_m \in C - - via the FISTA problem - - .. math:: - \mathrm{argmin}_\mathbf{d} \; - (1/2) \sum_k \left\| \sum_m \mathbf{d}_m * \mathbf{x}_{k,m} - - \mathbf{s}_k \right\|_2^2 + \sum_m \iota_C(\mathbf{d}_m) \;\;, - - where :math:`\iota_C(\cdot)` is the indicator function of feasible - set :math:`C` consisting of filters with unit norm and constrained - support. Multi-channel problems with input image channels - :math:`\mathbf{s}_{c,k}` are also supported, either as - - .. math:: - \mathrm{argmin}_\mathbf{d} \; (1/2) \sum_c \sum_k \left\| \sum_m - \mathbf{d}_m * \mathbf{x}_{c,k,m} - \mathbf{s}_{c,k} \right\|_2^2 - \quad \text{such that} \quad \mathbf{d}_m \in C - - with single-channel dictionary filters :math:`\mathbf{d}_m` and - multi-channel coefficient maps :math:`\mathbf{x}_{c,k,m}`, or - - .. math:: - \mathrm{argmin}_\mathbf{d} \; (1/2) \sum_c \sum_k \left\| \sum_m - \mathbf{d}_{c,m} * \mathbf{x}_{k,m} - \mathbf{s}_{c,k} - \right\|_2^2 \quad \text{such that} \quad \mathbf{d}_{c,m} \in C - - with multi-channel dictionary filters :math:`\mathbf{d}_{c,m}` and - single-channel coefficient maps :math:`\mathbf{x}_{k,m}`. In this - latter case, normalisation of filters :math:`\mathbf{d}_{c,m}` is - performed jointly over index :math:`c` for each filter :math:`m`. - - After termination of the :meth:`solve` method, attribute :attr:`itstat` - is a list of tuples representing statistics of each iteration. The - fields of the named tuple ``IterationStats`` are: - - ``Iter`` : Iteration number - - ``DFid`` : Value of data fidelity term :math:`(1/2) \sum_k \| - \sum_m \mathbf{d}_m * \mathbf{x}_{k,m} - \mathbf{s}_k \|_2^2` - - ``Cnstr`` : Constraint violation measure - - ``Rsdl`` : Residual - - ``L`` : Inverse of gradient step parameter - - ``Time`` : Cumulative run time - """ - - - - class Options(fista.FISTADFT.Options): - r"""ConvCnstrMOD algorithm options - - Options include all of those defined in - :class:`.fista.FISTADFT.Options`, together with - additional options: - - ``ZeroMean`` : Flag indicating whether the solution - dictionary :math:`\{\mathbf{d}_m\}` should have zero-mean - components. - """ - - defaults = copy.deepcopy(fista.FISTADFT.Options.defaults) - defaults.update({'ZeroMean': False}) - - - def __init__(self, opt=None): - """ - Parameters - ---------- - opt : dict or None, optional (default None) - ConvCnstrMOD algorithm options - """ - - if opt is None: - opt = {} - fista.FISTADFT.Options.__init__(self, opt) - - - def __setitem__(self, key, value): - """Set options.""" - - fista.FISTADFT.Options.__setitem__(self, key, value) - - - itstat_fields_objfn = ('DFid', 'Cnstr') - hdrtxt_objfn = ('DFid', 'Cnstr') - hdrval_objfun = {'DFid': 'DFid', 'Cnstr': 'Cnstr'} - - - - def __init__(self, Z, S, dsz, opt=None, dimK=1, dimN=2): - """ - This class supports an arbitrary number of spatial dimensions, - `dimN`, with a default of 2. The input coefficient map array `Z` - (usually labelled X, but renamed here to avoid confusion with - the X and Y variables in the FISTA base class) is expected to - be in standard form as computed by the GenericConvBPDN class. - - The input signal set `S` is either `dimN` dimensional (no - channels, only one signal), `dimN` +1 dimensional (either - multiple channels or multiple signals), or `dimN` +2 dimensional - (multiple channels and multiple signals). Parameter `dimK`, with - a default value of 1, indicates the number of multiple-signal - dimensions in `S`: - - :: - - Default dimK = 1, i.e. assume input S is of form - S(N0, N1, C, K) or S(N0, N1, K) - If dimK = 0 then input S is of form - S(N0, N1, C, K) or S(N0, N1, C) - - The internal data layout for S, D (X here), and X (Z here) is: - :: - - dim<0> - dim : Spatial dimensions, product of N0,N1,... is N - dim : C number of channels in S and D - dim : K number of signals in S - dim : M number of filters in D - - sptl. chn sig flt - S(N0, N1, C, K, 1) - D(N0, N1, C, 1, M) (X here) - X(N0, N1, 1, K, M) (Z here) - - The `dsz` parameter indicates the desired filter supports in the - output dictionary, since this cannot be inferred from the - input variables. The format is the same as the `dsz` parameter - of :func:`.cnvrep.bcrop`. - - - Parameters - ---------- - Z : array_like - Coefficient map array - S : array_like - Signal array - dsz : tuple - Filter support size(s) - opt : ccmod.Options object - Algorithm options - dimK : int, optional (default 1) - Number of dimensions for multiple signals in input S - dimN : int, optional (default 2) - Number of spatial dimensions - """ - - # Set default options if none specified - if opt is None: - opt = ConvCnstrMOD.Options() - - # Infer problem dimensions and set relevant attributes of self - self.cri = CDU_ConvRepIndexing(dsz, S, dimK=dimK, dimN=dimN) - - # Call parent class __init__ - xshape = self.cri.shpD - super(ConvCnstrMOD, self).__init__(xshape, S.dtype, opt) - - # Set gradient step parameter - self.set_attr('L', opt['L'], dval=self.cri.K * 14.0, dtype=self.dtype) - - # Reshape S to standard layout (Z, i.e. X in cbpdn, is assumed - # to be taken from cbpdn, and therefore already in standard - # form). If the dictionary has a single channel but the input - # (and therefore also the coefficient map array) has multiple - # channels, the channel index and multiple image index have - # the same behaviour in the dictionary update equation: the - # simplest way to handle this is to just reshape so that the - # channels also appear on the multiple image index. - if self.cri.Cd == 1 and self.cri.C > 1: - self.S = S.reshape(self.cri.Nv + (1,) + - (self.cri.C * self.cri.K,) + (1,)) - else: - self.S = S.reshape(self.cri.shpS) - self.S = np.asarray(self.S, dtype=self.dtype) - - # Compute signal S in DFT domain - self.Sf = rfftn(self.S, None, self.cri.axisN) - - # Create constraint set projection function - self.Pcn = getPcn(dsz, self.cri.Nv, self.cri.dimN, self.cri.dimCd, - zm=opt['ZeroMean']) - - # Create byte aligned arrays for FFT calls - self.Y = self.X - self.X = empty_aligned(self.Y.shape, dtype=self.dtype) - self.X[:] = self.Y - - # Initialise auxiliary variable Vf: Create byte aligned arrays - # for FFT calls - self.Vf = rfftn_empty_aligned(self.X.shape, self.cri.axisN, - self.dtype) - - self.Xf = rfftn(self.X, None, self.cri.axisN) - self.Yf = self.Xf - self.store_prev() - self.Yfprv = self.Yf.copy() + 1e5 - - # Initialization needed for back tracking (if selected) - self.postinitialization_backtracking_DFT() - - if Z is not None: - self.setcoef(Z) - - - - def setcoef(self, Z): - """Set coefficient array.""" - - # If the dictionary has a single channel but the input (and - # therefore also the coefficient map array) has multiple - # channels, the channel index and multiple image index have - # the same behaviour in the dictionary update equation: the - # simplest way to handle this is to just reshape so that the - # channels also appear on the multiple image index. - if self.cri.Cd == 1 and self.cri.C > 1: - Z = Z.reshape(self.cri.Nv + (1,) + (self.cri.Cx * self.cri.K,) + - (self.cri.M,)) - self.Z = np.asarray(Z, dtype=self.dtype) - - self.Zf = rfftn(self.Z, self.cri.Nv, self.cri.axisN) - - - - def getdict(self, crop=True): - """Get final dictionary. If ``crop`` is ``True``, apply - :func:`.cnvrep.bcrop` to returned array. - """ - - D = self.X - if crop: - D = bcrop(D, self.cri.dsz, self.cri.dimN) - return D - - - - def eval_grad(self): - """Compute gradient in Fourier domain.""" - - # Compute X D - S - Ryf = self.eval_Rf(self.Yf) - - gradf = inner(np.conj(self.Zf), Ryf, axis=self.cri.axisK) - - # Multiple channel signal, single channel dictionary - if self.cri.C > 1 and self.cri.Cd == 1: - gradf = np.sum(gradf, axis=self.cri.axisC, keepdims=True) - - return gradf - - - - def eval_Rf(self, Vf): - """Evaluate smooth term in Vf.""" - - return inner(self.Zf, Vf, axis=self.cri.axisM) - self.Sf - - - - def eval_proxop(self, V): - """Compute proximal operator of :math:`g`.""" - - return self.Pcn(V) - - - - def rsdl(self): - """Compute fixed point residual in Fourier domain.""" - - diff = self.Xf - self.Yfprv - return rfl2norm2(diff, self.X.shape, axis=self.cri.axisN) - - - - def eval_objfn(self): - """Compute components of objective function as well as total - contribution to objective function. - """ - - dfd = self.obfn_dfd() - cns = self.obfn_cns() - return (dfd, cns) - - - - def obfn_dfd(self): - r"""Compute data fidelity term :math:`(1/2) \| \sum_m - \mathbf{d}_m * \mathbf{x}_m - \mathbf{s} \|_2^2`. - """ - - Ef = self.eval_Rf(self.Xf) - return rfl2norm2(Ef, self.S.shape, axis=self.cri.axisN) / 2.0 - - - - def obfn_cns(self): - r"""Compute constraint violation measure :math:`\| - P(\mathbf{y}) - \mathbf{y}\|_2`. - """ - - return np.linalg.norm((self.Pcn(self.X) - self.X)) - - - - def obfn_f(self, Xf=None): - r"""Compute data fidelity term :math:`(1/2) \| \sum_m - \mathbf{d}_m * \mathbf{x}_m - \mathbf{s} \|_2^2`. - This is used for backtracking. Since the backtracking is - computed in the DFT, it is important to preserve the - DFT scaling. - """ - if Xf is None: - Xf = self.Xf - - Rf = self.eval_Rf(Xf) - return 0.5 * np.linalg.norm(Rf.flatten(), 2)**2 - - - - def reconstruct(self, D=None): - """Reconstruct representation.""" - - if D is None: - Df = self.Xf - else: - Df = rfftn(D, None, self.cri.axisN) - - Sf = np.sum(self.Zf * Df, axis=self.cri.axisM) - return irfftn(Sf, self.cri.Nv, self.cri.axisN) - - - - -class ConvCnstrMODMask(ConvCnstrMOD): - r""" - FISTA algorithm for Convolutional Constrained MOD problem - with a spatial mask :cite:`garcia-2018-convolutional1`. - - | - - .. inheritance-diagram:: ConvCnstrMODMask - :parts: 2 - - | - - Solve the optimisation problem - - .. math:: - \mathrm{argmin}_\mathbf{d} \; - (1/2) \left\| W \left(\sum_m \mathbf{d}_m * \mathbf{x}_m - - \mathbf{s}\right) \right\|_2^2 \quad \text{such that} \quad - \mathbf{d}_m \in C \;\; \forall m - - where :math:`C` is the feasible set consisting of filters with unit - norm and constrained support, and :math:`W` is a mask array, via the - FISTA problem - - .. math:: - \mathrm{argmin}_{\mathbf{d}} \; (1/2) \left\| W \left(X - \mathbf{d} - \mathbf{s}\right) \right\|_2^2 + - \iota_C(\mathbf{d}_m) \;\;, - - where :math:`\iota_C(\cdot)` is the indicator function of feasible - set :math:`C`, and :math:`X \mathbf{d} = \sum_m \mathbf{x}_m * - \mathbf{d}_m`. - - See :class:`ConvCnstrMOD` for interface details. - """ - - - class Options(ConvCnstrMOD.Options): - """ConvCnstrMODMask algorithm options - - Options include all of those defined in - :class:`.fista.FISTA.Options`. - """ - - defaults = copy.deepcopy(ConvCnstrMOD.Options.defaults) - - - def __init__(self, opt=None): - """ - Parameters - ---------- - opt : dict or None, optional (default None) - ConvCnstrMODMask algorithm options - """ - - if opt is None: - opt = {} - ConvCnstrMOD.Options.__init__(self, opt) - - - def __init__(self, Z, S, W, dsz, opt=None, dimK=None, dimN=2): - """ - Parameters - ---------- - Z : array_like - Coefficient map array - S : array_like - Signal array - W : array_like - Mask array. The array shape must be such that the array is - compatible for multiplication with the *internal* shape of - input array S (see :class:`.cnvrep.CDU_ConvRepIndexing` for a - discussion of the distinction between *external* and - *internal* data layouts). - dsz : tuple - Filter support size(s) - opt : :class:`ConvCnstrMODMask.Options` object - Algorithm options - dimK : 0, 1, or None, optional (default None) - Number of dimensions in input signal corresponding to multiple - independent signals - dimN : int, optional (default 2) - Number of spatial dimensions - """ - - # Set default options if none specified - if opt is None: - opt = ConvCnstrMODMask.Options() - - # Infer problem dimensions and set relevant attributes of self - self.cri = CDU_ConvRepIndexing(dsz, S, dimK=dimK, dimN=dimN) - - # Append singleton dimensions to W if necessary - if hasattr(W, 'ndim'): - W = atleast_nd(self.cri.dimN + 3, W) - - # Reshape W if necessary (see discussion of reshape of S in - # ccmod base class) - if self.cri.Cd == 1 and self.cri.C > 1 and hasattr(W, 'ndim'): - # In most cases broadcasting rules make it possible for W - # to have a singleton dimension corresponding to a - # non-singleton dimension in S. However, when S is - # reshaped to interleave axisC and axisK on the same axis, - # broadcasting is no longer sufficient unless axisC and - # axisK of W are either both singleton or both of the same - # size as the corresponding axes of S. If neither of these - # cases holds, it is necessary to replicate the axis of W - # (axisC or axisK) that does not have the same size as the - # corresponding axis of S. - shpw = list(W.shape) - swck = shpw[self.cri.axisC] * shpw[self.cri.axisK] - if swck > 1 and swck < self.cri.C * self.cri.K: - if W.shape[self.cri.axisK] == 1 and self.cri.K > 1: - shpw[self.cri.axisK] = self.cri.K - else: - shpw[self.cri.axisC] = self.cri.C - W = np.broadcast_to(W, shpw) - self.W = W.reshape( - W.shape[0:self.cri.dimN] + - (1, W.shape[self.cri.axisC] * W.shape[self.cri.axisK], 1)) - else: - self.W = W - - super(ConvCnstrMODMask, self).__init__(Z, S, dsz, opt, dimK, dimN) - - # Create byte aligned arrays for FFT calls - self.WRy = empty_aligned(self.S.shape, dtype=self.dtype) - self.Ryf = rfftn_empty_aligned(self.S.shape, self.cri.axisN, - self.dtype) - - - - def eval_grad(self): - """Compute gradient in Fourier domain.""" - - # Compute X D - S - self.Ryf[:] = self.eval_Rf(self.Yf) - - # Map to spatial domain to multiply by mask - Ry = irfftn(self.Ryf, self.cri.Nv, self.cri.axisN) - # Multiply by mask - self.WRy[:] = (self.W**2) * Ry - # Map back to frequency domain - WRyf = rfftn(self.WRy, self.cri.Nv, self.cri.axisN) - - gradf = inner(np.conj(self.Zf), WRyf, axis=self.cri.axisK) - - # Multiple channel signal, single channel dictionary - if self.cri.C > 1 and self.cri.Cd == 1: - gradf = np.sum(gradf, axis=self.cri.axisC, keepdims=True) - - return gradf - - - - def obfn_dfd(self): - r"""Compute data fidelity term :math:`(1/2) \sum_k \| W (\sum_m - \mathbf{d}_m * \mathbf{x}_{k,m} - \mathbf{s}_k) \|_2^2` - """ - - Ef = self.eval_Rf(self.Xf) - E = irfftn(Ef, self.cri.Nv, self.cri.axisN) - - return (np.linalg.norm(self.W * E)**2) / 2.0 - - - - def obfn_f(self, Xf=None): - r"""Compute data fidelity term :math:`(1/2) \sum_k \| W (\sum_m - \mathbf{d}_m * \mathbf{x}_{k,m} - \mathbf{s}_k) \|_2^2`. - This is used for backtracking. Since the backtracking is - computed in the DFT, it is important to preserve the - DFT scaling. - """ - - if Xf is None: - Xf = self.Xf - - Rf = self.eval_Rf(Xf) - R = irfftn(Rf, self.cri.Nv, self.cri.axisN) - WRf = rfftn(self.W * R, self.cri.Nv, self.cri.axisN) - - return 0.5 * np.linalg.norm(WRf.flatten(), 2)**2 diff --git a/sporco/fista/fista.py b/sporco/fista/fista.py deleted file mode 100644 index 511bf2f8..00000000 --- a/sporco/fista/fista.py +++ /dev/null @@ -1,936 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright (C) 2016-2019 by Cristina Garcia-Cardona -# Brendt Wohlberg -# All rights reserved. BSD 3-clause License. -# This file is part of the SPORCO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -"""Base classes for FISTA algorithms""" - -from __future__ import division, print_function -from builtins import range - -import copy -import numpy as np - -from sporco.cdict import ConstrainedDict -from sporco.common import IterativeSolver, solve_status_str -from sporco.fft import rfftn, irfftn -from sporco.array import transpose_ntpl_list -from sporco.util import Timer - - -__author__ = """Cristina Garcia-Cardona """ - - - -class FISTA(IterativeSolver): - r"""Base class for Fast Iterative Shrinkage/Thresholding algorithm - (FISTA) algorithms :cite:`beck-2009-fast`. A robust variant - :cite:`florea-2017-robust` is also supported. - - Solve optimisation problems of the form - - .. math:: - \mathrm{argmin}_{\mathbf{x}} \; f(\mathbf{x}) + g(\mathbf{x}) \;\;, - - where :math:`f, g` are convex functions and :math:`f` is smooth. - - This class is intended to be a base class of other classes that - specialise to specific optimisation problems. - - After termination of the :meth:`solve` method, attribute - :attr:`itstat` is a list of tuples representing statistics of each - iteration. The default fields of the named tuple - ``IterationStats`` are: - - ``Iter`` : Iteration number - - ``ObjFun`` : Objective function value - - ``FVal`` : Value of smooth objective function component :math:`f` - - ``GVal`` : Value of objective function component :math:`g` - - ``F_Btrack`` : Value of objective function :math:`f + g` - (see Sec. 2.2 of :cite:`beck-2009-fast`) - - ``Q_Btrack`` : Value of Quadratic approximation :math:`Q_L` - (see Sec. 2.3 of :cite:`beck-2009-fast`) - - ``IterBtrack`` : Number of iterations in backtracking - - ``Rsdl`` : Residual - - ``L`` : Inverse of gradient step parameter - - ``Time`` : Cumulative run time - """ - - class Options(ConstrainedDict): - r"""ADMM algorithm options. - - Options: - - ``FastSolve`` : Flag determining whether non-essential - computation is skipped. When ``FastSolve`` is ``True`` and - ``Verbose`` is ``False``, the functional value and related - iteration statistics are not computed. If ``FastSolve`` is - ``True`` residuals are also not calculated, in which case the - residual-based stopping method is also disabled, with the - number of iterations determined only by ``MaxMainIter``. - - ``Verbose`` : Flag determining whether iteration status is - displayed. - - ``StatusHeader`` : Flag determining whether status header and - separator are displayed. - - ``DataType`` : Specify data type for solution variables, - e.g. ``np.float32``. - - ``X0`` : Initial value for X variable. - - ``Callback`` : Callback function to be called at the end of - every iteration. - - ``MaxMainIter`` : Maximum main iterations. - - ``IterTimer`` : Label of the timer to use for iteration times. - - ``RelStopTol`` : Relative convergence tolerance for fixed point - residual (see Sec. 4.3 of :cite:`liu-2018-first`). - - ``L`` : Inverse of gradient step parameter :math:`L`. - - ``AutoStop`` : Options for adaptive stoping strategy (fixed - point residual, see Sec. 4.3 of :cite:`liu-2018-first`). - - ``Enabled`` : Flag determining whether the adaptive stopping - relative parameter strategy is enabled. - - ``Tau0`` : numerator in adaptive criterion - (:math:`\tau_0` in :cite:`liu-2018-first`). - - ``BackTrack`` : Options for adaptive L strategy (backtracking, - see Sec. 4 of :cite:`beck-2009-fast` or Robust Fista - in :cite:`florea-2017-robust`). - - ``Enabled`` : Flag determining whether adaptive inverse step - size parameter strategy is enabled. When true, backtracking - in Sec. 4 of :cite:`beck-2009-fast` is used. In combination with - the ``Robust`` flag it enables the backtracking strategy in - :cite:`florea-2017-robust`. - - ``Robust`` : Flag determining if the robust FISTA update is to be - applied as in :cite:`florea-2017-robust`. - - ``gamma_d`` : Multiplier applied to decrease L when backtracking in - robust FISTA (:math:`\gamma_d` in :cite:`florea-2017-robust`). - - ``gamma_u`` : Multiplier applied to increase L when backtracking in - standard FISTA (corresponding to :math:`\eta` in - :cite:`beck-2009-fast`) or in robust FISTA (corresponding Total - :math:`\gamma_u` in :cite:`florea-2017-robust`). - - ``MaxIter`` : Maximum iterations of updating L when - backtracking. - """ - - defaults = {'FastSolve': False, 'Verbose': False, - 'StatusHeader': True, 'DataType': None, - 'X0': None, 'Callback': None, - 'MaxMainIter': 1000, 'IterTimer': 'solve', - 'RelStopTol': 1e-3, 'L': None, - 'BackTrack': - {'Enabled': False, 'Robust': False, - 'gamma_d': 0.9, 'gamma_u': 1.2, 'MaxIter': 100}, - 'AutoStop': {'Enabled': False, 'Tau0': 1e-2}} - - def __init__(self, opt=None): - """ - Parameters - ---------- - opt : dict or None, optional (default None) - FISTA algorithm options - """ - - if opt is None: - opt = {} - ConstrainedDict.__init__(self, opt) - - - - fwiter = 4 - """Field width for iteration count display column""" - fpothr = 2 - """Field precision for other display columns""" - - itstat_fields_objfn = ('ObjFun', 'FVal', 'GVal') - """Fields in IterationStats associated with the objective function; - see :meth:`eval_objfun`""" - itstat_fields_alg = ('Rsdl', 'F_Btrack', 'Q_Btrack', 'IterBTrack', 'L') - """Fields in IterationStats associated with the specific solver - algorithm""" - itstat_fields_extra = () - """Non-standard fields in IterationStats; see :meth:`itstat_extra`""" - - hdrtxt_objfn = ('Fnc', 'f', 'g') - """Display column headers associated with the objective function; - see :meth:`eval_objfun`""" - hdrval_objfun = {'Fnc': 'ObjFun', 'f': 'FVal', 'g': 'GVal'} - """Dictionary mapping display column headers in :attr:`hdrtxt_objfn` - to IterationStats entries""" - - - - def __new__(cls, *args, **kwargs): - """Create a FISTA object and start its initialisation timer.""" - - instance = super(FISTA, cls).__new__(cls) - instance.timer = Timer(['init', 'solve', 'solve_wo_func', - 'solve_wo_rsdl', 'solve_wo_btrack']) - instance.timer.start('init') - return instance - - - - def __init__(self, xshape, dtype, opt=None): - r""" - Parameters - ---------- - xshape : tuple of ints - Shape of working variable X - dtype : data-type - Data type for working variables (overridden by 'DataType' option) - opt : :class:`FISTA.Options` object - Algorithm options - """ - - if opt is None: - opt = FISTA.Options() - if not isinstance(opt, FISTA.Options): - raise TypeError("Parameter opt must be an instance of " - "FISTA.Options") - - self.opt = opt - - # DataType option overrides data type inferred from __init__ - # parameters of derived class - self.set_dtype(opt, dtype) - - # Initialise attributes representing step parameter and other - # parameters - self.set_attr('L', opt['L'], dval=1.0, dtype=self.dtype) - dval_gamma_u = 1.2 - if self.opt['BackTrack', 'Robust']: - dval_gamma_u = 2. - self.set_attr('L_gamma_u', opt['BackTrack', 'gamma_u'], - dval=dval_gamma_u, dtype=self.dtype) - self.set_attr('L_gamma_d', opt['BackTrack', 'gamma_d'], dval=0.9, - dtype=self.dtype) - self.set_attr('L_maxiter', opt['BackTrack', 'MaxIter'], dval=1.0, - dtype=self.dtype) - - # If using adaptative stop criterion, set tau0 parameter - if self.opt['AutoStop', 'Enabled']: - self.tau0 = self.opt['AutoStop', 'Tau0'] - - # Initialise working variable X - if self.opt['X0'] is None: - self.X = self.xinit(xshape) - else: - self.X = self.opt['X0'].astype(self.dtype, copy=True) - - # Default values for variables created only if BackTrack is enabled - if self.opt['BackTrack', 'Enabled']: - self.F = 0. - self.Q = 0. - self.iterBTrack = 0 - # Determine type of backtracking - if self.opt['BackTrack', 'Robust']: - self.Tk = 0. - self.zzinit() - self.backtracking = self.robust_backtrack - else: - self.t = 1. - self.backtracking = self.standard_backtrack - else: - self.F = None - self.Q = None - self.iterBTrack = None - self.t = 1 - - self.Y = None - - self.itstat = [] - self.k = 0 - - - - def xinit(self, xshape): - """Return initialiser for working variable X.""" - - return np.zeros(xshape, dtype=self.dtype) - - - - def zzinit(self): - """Return initialiser for working variable ZZ (required for - robust FISTA). - """ - - self.ZZ = self.X.copy() - - - - def solve(self): - """Start (or re-start) optimisation. This method implements the - framework for the iterations of a FISTA algorithm. There is - sufficient flexibility in overriding the component methods that - it calls that it is usually not necessary to override this method - in derived clases. - - If option ``Verbose`` is ``True``, the progress of the - optimisation is displayed at every iteration. At termination - of this method, attribute :attr:`itstat` is a list of tuples - representing statistics of each iteration, unless option - ``FastSolve`` is ``True`` and option ``Verbose`` is ``False``. - - Attribute :attr:`timer` is an instance of :class:`.util.Timer` - that provides the following labelled timers: - - ``init``: Time taken for object initialisation by - :meth:`__init__` - - ``solve``: Total time taken by call(s) to :meth:`solve` - - ``solve_wo_func``: Total time taken by call(s) to - :meth:`solve`, excluding time taken to compute functional - value and related iteration statistics - - ``solve_wo_rsdl`` : Total time taken by call(s) to - :meth:`solve`, excluding time taken to compute functional - value and related iteration statistics as well as time take - to compute residuals - - ``solve_wo_btrack`` : Total time taken by call(s) to - :meth:`solve`, excluding time taken to compute functional - value and related iteration statistics as well as time take - to compute residuals and implemented ``BackTrack`` mechanism - """ - - # Open status display - fmtstr, nsep = self.display_start() - - # Start solve timer - self.timer.start(['solve', 'solve_wo_func', 'solve_wo_rsdl', - 'solve_wo_btrack']) - - # Main optimisation iterations - for self.k in range(self.k, self.k + self.opt['MaxMainIter']): - - # Update record of X from previous iteration - self.store_prev() - - # Compute backtracking - if self.opt['BackTrack', 'Enabled'] and self.k >= 0: - self.timer.stop('solve_wo_btrack') - # Compute backtracking - self.backtracking() - self.timer.start('solve_wo_btrack') - else: - # Compute just proximal step - self.proximal_step() - # Update by combining previous iterates - self.combination_step() - - # Compute residuals and stopping thresholds - self.timer.stop(['solve_wo_rsdl', 'solve_wo_btrack']) - if not self.opt['FastSolve']: - frcxd, adapt_tol = self.compute_residuals() - self.timer.start('solve_wo_rsdl') - - # Compute and record other iteration statistics and - # display iteration stats if Verbose option enabled - self.timer.stop(['solve_wo_func', 'solve_wo_rsdl', - 'solve_wo_btrack']) - if not self.opt['FastSolve']: - itst = self.iteration_stats(self.k, frcxd) - self.itstat.append(itst) - self.display_status(fmtstr, itst) - self.timer.start(['solve_wo_func', 'solve_wo_rsdl', - 'solve_wo_btrack']) - - # Call callback function if defined - if self.opt['Callback'] is not None: - if self.opt['Callback'](self): - break - - # Stop if residual-based stopping tolerances reached - if not self.opt['FastSolve']: - if frcxd < adapt_tol: - break - - # Increment iteration count - self.k += 1 - - # Record solve time - self.timer.stop(['solve', 'solve_wo_func', 'solve_wo_rsdl', - 'solve_wo_btrack']) - - # Print final separator string if Verbose option enabled - self.display_end(nsep) - - return self.getmin() - - - - def getmin(self): - """Get minimiser after optimisation.""" - - return self.X - - - - def proximal_step(self, grad=None): - """Compute proximal update (gradient descent + regularization).""" - - if grad is None: - grad = self.eval_grad() - - V = self.Y - (1. / self.L) * grad - - self.X = self.eval_proxop(V) - - return grad - - - - def combination_step(self): - """Build next update by a smart combination of previous updates - (standard FISTA :cite:`beck-2009-fast`). - """ - - # Update t step - tprv = self.t - self.t = 0.5 * float(1. + np.sqrt(1. + 4. * tprv**2)) - - # Update Y - if not self.opt['FastSolve']: - self.Yprv = self.Y.copy() - self.Y = self.X + ((tprv - 1.) / self.t) * (self.X - self.Xprv) - - - - def standard_backtrack(self): - """Estimate step size L by computing a linesearch that - guarantees that F <= Q according to the standard FISTA - backtracking strategy in :cite:`beck-2009-fast`. - This also updates variable Y. - """ - - gradY = self.eval_grad() # Given Y(f), this updates computes gradY(f) - - maxiter = self.L_maxiter - - iterBTrack = 0 - linesearch = 1 - while linesearch and iterBTrack < maxiter: - - self.proximal_step(gradY) # Given gradY(f), L, this updates X(f) - - f = self.obfn_f(self.var_x()) - Dxy = self.eval_Dxy() - Q = self.obfn_f(self.var_y()) + \ - self.eval_linear_approx(Dxy, gradY) + \ - (self.L / 2.) * np.linalg.norm(Dxy.flatten(), 2)**2 - - if f <= Q: - linesearch = 0 - else: - self.L *= self.L_gamma_u - - iterBTrack += 1 - - self.F = f - self.Q = Q - self.iterBTrack = iterBTrack - # Update auxiliary sequence - self.combination_step() - - - - def robust_backtrack(self): - """Estimate step size L by computing a linesearch that - guarantees that F <= Q according to the robust FISTA - backtracking strategy in :cite:`florea-2017-robust`. - This also updates all the supporting variables. - """ - - self.L *= self.L_gamma_d - maxiter = self.L_maxiter - - iterBTrack = 0 - linesearch = 1 - - self.store_Yprev() - while linesearch and iterBTrack < maxiter: - - t = float(1. + np.sqrt(1. + 4. * self.L * self.Tk)) / (2. * self.L) - T = self.Tk + t - y = (self.Tk * self.var_xprv() + t * self.ZZ) / T - self.update_var_y(y) - - gradY = self.proximal_step() # Given Y(f), L, this updates X(f) - - f = self.obfn_f(self.var_x()) - Dxy = self.eval_Dxy() - Q = self.obfn_f(self.var_y()) + \ - self.eval_linear_approx(Dxy, gradY) + \ - (self.L / 2.) * np.linalg.norm(Dxy.flatten(), 2)**2 - - if f <= Q: - linesearch = 0 - else: - self.L *= self.L_gamma_u - - iterBTrack += 1 - - self.Tk = T - self.ZZ += (t * self.L * (self.var_x() - self.var_y())) - - self.F = f - self.Q = Q - self.iterBTrack = iterBTrack - - - - def eval_linear_approx(self, Dxy, gradY): - r"""Compute term - :math:`\langle \nabla f(\mathbf{y}), \mathbf{x} - \mathbf{y} \rangle` - that is part of the quadratic function :math:`Q_L` used - for backtracking. - """ - - return np.sum(Dxy * gradY) - - - - def eval_grad(self): - """Compute gradient. - - Overriding this method is required. - """ - - raise NotImplementedError() - - - - def eval_proxop(self, V): - """Compute proximal operator of :math:`g`. - - Overriding this method is required. - """ - - raise NotImplementedError() - - - - def store_prev(self): - """Store previous X state.""" - - self.Xprv = self.X.copy() - - - - def store_Yprev(self): - """Store previous Y state.""" - - self.Yprv = self.Y.copy() - - - - def eval_Dxy(self): - """Evaluate difference of state and auxiliary state updates.""" - - return self.X - self.Y - - - - def compute_residuals(self): - """Compute residuals and stopping thresholds.""" - - r = self.rsdl() - adapt_tol = self.opt['RelStopTol'] - - if self.opt['AutoStop', 'Enabled']: - adapt_tol = self.tau0 / (1. + self.k) - - return r, adapt_tol - - - - @classmethod - def hdrtxt(cls): - """Construct tuple of status display column title.""" - - return ('Itn',) + cls.hdrtxt_objfn + ('Rsdl', 'F', 'Q', 'It_Bt', 'L') - - - - @classmethod - def hdrval(cls): - """Construct dictionary mapping display column title to - IterationStats entries. - """ - - dict = {'Itn': 'Iter'} - dict.update(cls.hdrval_objfun) - dict.update({'Rsdl': 'Rsdl', 'F': 'F_Btrack', 'Q': 'Q_Btrack', - 'It_Bt': 'IterBTrack', 'L': 'L'}) - - return dict - - - - def iteration_stats(self, k, frcxd): - """Construct iteration stats record tuple.""" - - tk = self.timer.elapsed(self.opt['IterTimer']) - tpl = (k,) + self.eval_objfn() + \ - (frcxd, self.F, self.Q, self.iterBTrack, self.L) + \ - self.itstat_extra() + (tk,) - return type(self).IterationStats(*tpl) - - - - def eval_objfn(self): - """Compute components of objective function as well as total - contribution to objective function. - """ - - fval = self.obfn_f(self.X) - gval = self.obfn_g(self.X) - obj = fval + gval - return (obj, fval, gval) - - - - def itstat_extra(self): - """Non-standard entries for the iteration stats record tuple.""" - - return () - - - - def getitstat(self): - """Get iteration stats as named tuple of arrays instead of - array of named tuples. - """ - - return transpose_ntpl_list(self.itstat) - - - - def display_start(self): - """Set up status display if option selected. NB: this method - assumes that the first entry is the iteration count and the - last is the L value. - """ - - if self.opt['Verbose']: - # If backtracking option enabled F, Q, itBT, L are - # included in iteration status - if self.opt['BackTrack', 'Enabled']: - hdrtxt = type(self).hdrtxt() - else: - hdrtxt = type(self).hdrtxt()[0:-4] - # Call utility function to construct status display formatting - hdrstr, fmtstr, nsep = solve_status_str( - hdrtxt, fmtmap={'It_Bt': '%5d'}, fwdth0=type(self).fwiter, - fprec=type(self).fpothr) - # Print header and separator strings - if self.opt['StatusHeader']: - print(hdrstr) - print("-" * nsep) - else: - fmtstr, nsep = '', 0 - - return fmtstr, nsep - - - - def display_status(self, fmtstr, itst): - """Display current iteration status as selection of fields from - iteration stats tuple. - """ - - if self.opt['Verbose']: - hdrtxt = type(self).hdrtxt() - hdrval = type(self).hdrval() - itdsp = tuple([getattr(itst, hdrval[col]) for col in hdrtxt]) - if not self.opt['BackTrack', 'Enabled']: - itdsp = itdsp[0:-4] - - print(fmtstr % itdsp) - - - - def display_end(self, nsep): - """Terminate status display if option selected.""" - - if self.opt['Verbose'] and self.opt['StatusHeader']: - print("-" * nsep) - - - - def var_x(self): - r"""Get :math:`\mathbf{x}` variable.""" - - return self.X - - - - def var_y(self): - r"""Get :math:`\mathbf{y}` variable.""" - - return self.Y - - - - def var_xprv(self): - r"""Get :math:`\mathbf{x}` variable of previous iteration.""" - - return self.Xprv - - - - def update_var_y(self, y): - r"""Update :math:`\mathbf{y}` variable.""" - - self.Y = y - - - - def obfn_f(self, X): - r"""Compute :math:`f(\mathbf{x})` component of FISTA objective - function. - - Overriding this method is required (even if :meth:`eval_objfun` - is overriden, since this method is required for backtracking). - """ - - raise NotImplementedError() - - - - def obfn_g(self, X): - r"""Compute :math:`g(\mathbf{x})` component of FISTA objective - function. - - Overriding this method is required if :meth:`eval_objfun` - is not overridden. - """ - - raise NotImplementedError() - - - - def rsdl(self): - """Compute fixed point residual (see Sec. 4.3 of - :cite:`liu-2018-first`).""" - - return np.linalg.norm((self.X - self.Yprv).ravel()) - - - - - -class FISTADFT(FISTA): - r""" - Base class for FISTA algorithms with gradients and updates computed - in the frequency domain. - - | - - .. inheritance-diagram:: FISTADFT - :parts: 2 - - | - - Solve optimisation problems of the form - - .. math:: - \mathrm{argmin}_{\mathbf{x}} \; f(\mathbf{x}) + g(\mathbf{x}) - \;\;, - - where :math:`f, g` are convex functions and :math:`f` is smooth. - - This class specialises class FISTA, but remains a base class for - other classes that specialise to specific optimisation problems. - """ - - - class Options(FISTA.Options): - """FISTADFT algorithm options. - - Options include all of those defined in :class:`FISTA.Options`. - """ - - defaults = copy.deepcopy(FISTA.Options.defaults) - - def __init__(self, opt=None): - """ - Parameters - ---------- - opt : dict or None, optional (default None) - FISTADFT algorithm options - """ - - if opt is None: - opt = {} - FISTA.Options.__init__(self, opt) - - - - - def __init__(self, xshape, dtype, opt=None): - """ - Parameters - ---------- - xshape : tuple of ints - Shape of working variable X (the primary variable) - dtype : data-type - Data type for working variables - opt : :class:`FISTADFT.Options` object - Algorithm options - """ - - if opt is None: - opt = FISTADFT.Options() - super(FISTADFT, self).__init__(xshape, dtype, opt) - - - - def postinitialization_backtracking_DFT(self): - r""" - Computes variables needed for backtracking when the updates - are made in the DFT. (This requires the variables in DFT to - have been initialized). - """ - - if self.opt['BackTrack', 'Enabled']: - if self.opt['BackTrack', 'Robust']: - self.zzfinit() - - - - def zzfinit(self): - """Return initialiser for working variable ZZ in frequency - domain (required for robust FISTA :cite:`florea-2017-robust`). - """ - - self.ZZ = self.Xf.copy() - - - - def proximal_step(self, gradf=None): - """Compute proximal update (gradient descent + constraint). - Variables are mapped back and forth between input and - frequency domains. - """ - - if gradf is None: - gradf = self.eval_grad() - - self.Vf[:] = self.Yf - (1. / self.L) * gradf - V = irfftn(self.Vf, self.cri.Nv, self.cri.axisN) - - self.X[:] = self.eval_proxop(V) - self.Xf = rfftn(self.X, None, self.cri.axisN) - - return gradf - - - - def combination_step(self): - """Update auxiliary state by a smart combination of previous - updates in the frequency domain (standard FISTA - :cite:`beck-2009-fast`). - """ - - # Update t step - tprv = self.t - self.t = 0.5 * float(1. + np.sqrt(1. + 4. * tprv**2)) - - # Update Y - if not self.opt['FastSolve']: - self.Yfprv = self.Yf.copy() - self.Yf = self.Xf + ((tprv - 1.) / self.t) * (self.Xf - self.Xfprv) - - - - def store_prev(self): - """Store previous X in frequency domain.""" - - self.Xfprv = self.Xf.copy() - - - - def store_Yprev(self): - """Store previous Y state in frequency domain.""" - - self.Yfprv = self.Yf.copy() - - - - def eval_Dxy(self): - """Evaluate difference of state and auxiliary state in - frequency domain. - """ - - return self.Xf - self.Yf - - - - def var_x(self): - r"""Get :math:`\mathbf{x}` variable in frequency domain.""" - - return self.Xf - - - - def var_y(self): - r"""Get :math:`\mathbf{y}` variable in frequency domain.""" - - return self.Yf - - - - def var_xprv(self): - r"""Get :math:`\mathbf{x}` variable of previous iteration in - frequency domain. - """ - - return self.Xfprv - - - - def update_var_y(self, y): - r"""Update :math:`\mathbf{y}` variable in frequency domain.""" - - self.Yf = y - - - - def eval_linear_approx(self, Dxy, gradY): - r"""Compute term :math:`\langle \nabla f(\mathbf{y}), - \mathbf{x} - \mathbf{y} \rangle` (in frequency domain) that is - part of the quadratic function :math:`Q_L` used for - backtracking. Since this class computes the backtracking in - the DFT, it is important to preserve the DFT scaling. - """ - - return np.sum(np.real(np.conj(Dxy) * gradY)) diff --git a/sporco/fista/ppp.py b/sporco/fista/ppp.py deleted file mode 100644 index 0776ee65..00000000 --- a/sporco/fista/ppp.py +++ /dev/null @@ -1,153 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright (C) 2019-2020 by Brendt Wohlberg -# Ulugbek Kamilov -# All rights reserved. BSD 3-clause License. -# This file is part of the SPORCO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -"""Classes for FISTA variant of the Plug and Play Priors (PPP) algorithm.""" - -from __future__ import division, absolute_import, print_function - -import numpy as np - -from sporco.fista import fista - - -__author__ = """\n""".join(['Brendt Wohlberg ', - 'Ulugbek Kamilov ']) - - - -class GenericPPP(fista.FISTA): - """Base class for Plug and Play Priors (PPP) FISTA solvers - :cite:`kamilov-2017-plugandplay`.""" - - def __init__(self, xshape, opt=None): - """ - Parameters - ---------- - xshape : tuple of ints - Shape of working variable X - opt : :class:`GenericPPP.Options` object - Algorithm options - """ - - if opt is None: - opt = GenericPPP.Options() - - # Set dtype attribute, default is np.float32 - self.set_dtype(opt, np.dtype(np.float32)) - - super(GenericPPP, self).__init__(xshape, self.dtype, opt) - - self.store_prev() - self.Y = self.X.copy() - self.Yprv = np.zeros(self.Y.shape) - - - - itstat_fields_objfn = ('FVal',) - hdrtxt_objfn = ('FVal',) - hdrval_objfun = {'FVal': 'FVal'} - - - - def eval_grad(self): - """Compute the gradient of :math:`f`.""" - - return self.gradf(self.Y) - - - - def eval_proxop(self, V): - """Compute proximal operator of :math:`g`.""" - - return self.proxg(V, self.L) - - - - def rsdl(self): - """Compute fixed point residual.""" - - return np.linalg.norm((self.X - self.Yprv).ravel()) - - - - def eval_objfn(self): - r"""Compute components of objective function. - - In this case the regularisation term is implicit so we can only - evaluate the data fidelity term represented by the - :math:`f(\cdot)` component of the functional to be minimised. - """ - - return (self.f(self.X),) - - - - def gradf(self, X): - r"""Compute the gradient of :math:`f(\cdot)`. - - Overriding this method is required. - """ - - raise NotImplementedError() - - - - def proxg(self, X, L): - r"""Compute the proximal operator of :math:`L^{-1} g(\cdot)`. - - Overriding this method is required. Note that this method - should compute the proximal operator of - :math:`L^{-1} g(\cdot)`, *not* the proximal operator - of :math:`L g(\cdot)`. - """ - - raise NotImplementedError() - - - - def f(self, X): - r"""Evauate the data fidelity term :math:`f(\mathbf{x})`. - - Overriding this method is required. - """ - - raise NotImplementedError() - - - - - -class PPP(GenericPPP): - """Plug and Play Priors (PPP) solver :cite:`kamilov-2017-plugandplay` - that can be used without the need to derive a new class.""" - - def __init__(self, xshape, f, gradf, proxg, opt=None): - """ - Parameters - ---------- - xshape : tuple of ints - Shape of working variable X - f : function - Function evaluating the data fidelity term - gradf : function - Function computing the gradient of the data fidelity term - proxg : function - Function computing the proximal operator of the regularisation - term - opt : :class:`PPP.Options` object - Algorithm options - """ - - if opt is None: - opt = PPP.Options() - - super(PPP, self).__init__(xshape, opt) - - self.f = f - self.gradf = gradf - self.proxg = proxg diff --git a/sporco/interp.py b/sporco/interp.py index be2bc97c..d47b990f 100644 --- a/sporco/interp.py +++ b/sporco/interp.py @@ -16,7 +16,6 @@ import scipy.optimize as sco from scipy.interpolate import interp2d, griddata -from sporco._util import renamed_function __author__ = """Brendt Wohlberg """ @@ -102,7 +101,6 @@ def _linprog(*args, **kwargs): -@renamed_function(depname='lstabsdev', depmod='sporco.util') def lstabsdev(A, b): r"""Least absolute deviations (LAD) linear regression. @@ -149,7 +147,6 @@ def lstabsdev(A, b): -@renamed_function(depname='lstmaxdev', depmod='sporco.util') def lstmaxdev(A, b): r"""Least maximum deviation (least maximum error) linear regression. @@ -196,7 +193,6 @@ def lstmaxdev(A, b): -@renamed_function(depname='lanczos_kernel', depmod='sporco.util') def lanczos_kernel(x, a=3): r"""Lanczos interpolation kernel. @@ -227,7 +223,6 @@ def lanczos_kernel(x, a=3): -@renamed_function(depname='interpolation_points', depmod='sporco.util') def interpolation_points(N, include_zero=True): """Evenly spaced interpolation points. @@ -255,7 +250,6 @@ def interpolation_points(N, include_zero=True): -@renamed_function(depname='lanczos_filters', depmod='sporco.util') def lanczos_filters(sz, a=3, collapse_axes=True): """Multi-dimensional Lanczos interpolation filters. diff --git a/sporco/linalg.py b/sporco/linalg.py index 4cc7008a..5436fedb 100644 --- a/sporco/linalg.py +++ b/sporco/linalg.py @@ -22,7 +22,6 @@ else: have_numexpr = True -from sporco._util import renamed_function from sporco.array import zdivide, subsample_array @@ -202,7 +201,6 @@ def valid_adjoint(A, AT, Ashape, ATshape, eps=1e-7): -@renamed_function(depname='blockcirculant') def block_circulant(A): """Construct a block circulant matrix from a tuple of arrays. @@ -925,7 +923,6 @@ def rrs(ax, b): -@renamed_function(depname='pca', depmod='sporco.util') def pca(U, centre=False): """Compute the PCA basis for columns of input array `U`. @@ -961,7 +958,6 @@ def pca(U, centre=False): -@renamed_function(depname='nkp', depmod='sporco.util') def nkp(A, bshape, cshape): r"""Solve the Nearest Kronecker Product problem. @@ -1003,7 +999,6 @@ def nkp(A, bshape, cshape): -@renamed_function(depname='kpsvd', depmod='sporco.util') def kpsvd(A, bshape, cshape): r"""Compute the Kronecker Product SVD. diff --git a/sporco/prox/_l21.py b/sporco/prox/_l21.py index 33bd5afb..c14896db 100644 --- a/sporco/prox/_l21.py +++ b/sporco/prox/_l21.py @@ -12,7 +12,6 @@ import numpy as np from ._lp import prox_l1, prox_l2, norm_l2 -from sporco._util import renamed_function __author__ = """Brendt Wohlberg """ @@ -49,7 +48,6 @@ def norm_l21(x, axis=-1): -@renamed_function(depname='prox_l1l2') def prox_sl1l2(v, alpha, beta, axis=None): r"""Compute the proximal operator of the sum of :math:`\ell_1` and :math:`\ell_2` norms (compound shrinkage/soft thresholding) diff --git a/sporco/signal.py b/sporco/signal.py index dded35f0..802b0913 100644 --- a/sporco/signal.py +++ b/sporco/signal.py @@ -12,7 +12,6 @@ import numpy as np -from sporco._util import renamed_function from sporco.fft import is_complex_dtype, fftn, ifftn, rfftn, irfftn, fftconv from sporco import array @@ -21,7 +20,6 @@ -@renamed_function(depname='complex_randn', depmod='sporco.util') def complex_randn(*args): """Return a complex array of samples drawn from a standard normal distribution. @@ -41,7 +39,6 @@ def complex_randn(*args): -@renamed_function(depname='spnoise', depmod='sporco.util') def spnoise(s, frc, smn=0.0, smx=1.0): """Return image with salt & pepper noise imposed on it. @@ -70,7 +67,6 @@ def spnoise(s, frc, smn=0.0, smx=1.0): -@renamed_function(depname='rndmask', depmod='sporco.util') def rndmask(shp, frc, dtype=None): r"""Return random mask image with values in :math:`\{0,1\}`. @@ -96,7 +92,6 @@ def rndmask(shp, frc, dtype=None): -@renamed_function(depname='rgb2gray', depmod='sporco.util') def rgb2gray(rgb): """Convert an RGB image (or images) to grayscale. @@ -117,7 +112,6 @@ def rgb2gray(rgb): -@renamed_function(depname='grad', depmod='sporco.linalg') def grad(x, axis, zero_pad=False): r"""Compute gradient of `x` along axis `axis`. @@ -175,7 +169,6 @@ def grad(x, axis, zero_pad=False): -@renamed_function(depname='gradT', depmod='sporco.linalg') def gradT(x, axis, zero_pad=False): """Compute transpose of gradient of `x` along axis `axis`. @@ -208,7 +201,6 @@ def gradT(x, axis, zero_pad=False): -@renamed_function(depname='gradient_filters', depmod='sporco.linalg') def gradient_filters(ndim, axes, axshp, dtype=None): r"""Construct a set of filters for computing gradients in the frequency domain. @@ -249,7 +241,6 @@ def gradient_filters(ndim, axes, axshp, dtype=None): -@renamed_function(depname='tikhonov_filter', depmod='sporco.util') def tikhonov_filter(s, lmbda, npd=16): r"""Lowpass filter based on Tikhonov regularization. @@ -311,7 +302,6 @@ def tikhonov_filter(s, lmbda, npd=16): -@renamed_function(depname='gaussian', depmod='sporco.util') def gaussian(shape, sd=1.0): """Sample a multivariate Gaussian pdf, normalised to have unit sum. @@ -342,7 +332,6 @@ def gaussian(shape, sd=1.0): -@renamed_function(depname='local_contrast_normalise', depmod='sporco.util') def local_contrast_normalise(s, n=7, c=None): """Local contrast normalisation of an image.