#!/usr/bin/env python
# -*- coding: utf-8 -*-
# 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.

"""
Basis Pursuit DeNoising
=======================

This example demonstrates the use of class :class:`.admm.bpdn.BPDN` to solve the Basis Pursuit DeNoising (BPDN) problem :cite:`chen-1998-atomic`

  $$\mathrm{argmin}_\mathbf{x} \; (1/2) \| D \mathbf{x} - \mathbf{s} \|_2^2 + \lambda \| \mathbf{x} \|_1 \;,$$

where $D$ is the dictionary, $\mathbf{x}$ is the sparse representation, and $\mathbf{s}$ is the signal to be represented. In this example the BPDN problem is used to estimate the reference sparse representation that generated a signal from a noisy version of the signal.
"""


from __future__ import print_function
from builtins import input

import numpy as np

from sporco.admm import bpdn
from sporco import util
from sporco import plot


"""
Configure problem size, sparsity, and noise level.
"""

N = 512      # Signal size
M = 4*N      # Dictionary size
L = 32       # Number of non-zero coefficients in generator
sigma = 0.5  # Noise level


"""
Construct random dictionary, reference random sparse representation, and test signal consisting of the synthesis of the reference sparse representation with additive Gaussian noise.
"""

# Construct random dictionary and random sparse coefficients
np.random.seed(12345)
D = np.random.randn(N, M)
x0 = np.zeros((M, 1))
si = np.random.permutation(list(range(0, M-1)))
x0[si[0:L]] = np.random.randn(L, 1)

# Construct reference and noisy signal
s0 = D.dot(x0)
s = s0 + sigma*np.random.randn(N,1)


"""
Set BPDN solver class options.
"""

opt = bpdn.BPDN.Options({'Verbose': False, 'MaxMainIter': 500,
                    'RelStopTol': 1e-3, 'AutoRho': {'RsdlTarget': 1.0}})


"""
Select regularization parameter $\lambda$ by evaluating the error in recovering the sparse representation over a logarithmicaly spaced grid. (The reference representation is assumed to be known, which is not realistic in a real application.) A function is defined that evalues the BPDN recovery error for a specified $\lambda$, and this function is evaluated in parallel by :func:`sporco.util.grid_search`.
"""

# Function computing reconstruction error at lmbda
def evalerr(prm):
    lmbda = prm[0]
    b = bpdn.BPDN(D, s, lmbda, opt)
    x = b.solve()
    return np.sum(np.abs(x-x0))


# Parallel evalution of error function on lmbda grid
lrng = np.logspace(1, 2, 20)
sprm, sfvl, fvmx, sidx = util.grid_search(evalerr, (lrng,))
lmbda = sprm[0]

print('Minimum ℓ1 error: %5.2f at 𝜆 = %.2e' % (sfvl, lmbda))


"""
Once the best $\lambda$ has been determined, run BPDN with verbose display of ADMM iteration statistics.
"""

# Initialise and run BPDN object for best lmbda
opt['Verbose'] = True
b = bpdn.BPDN(D, s, lmbda, opt)
x = b.solve()

print("BPDN solve time: %.2fs" % b.timer.elapsed('solve'))


"""
Plot comparison of reference and recovered representations.
"""

plot.plot(np.hstack((x0, x)), title='Sparse representation',
          lgnd=['Reference', 'Reconstructed'])


"""
Plot lmbda error curve, functional value, residuals, and rho
"""

its = b.getitstat()
fig = plot.figure(figsize=(15, 10))
plot.subplot(2, 2, 1)
plot.plot(fvmx, x=lrng, ptyp='semilogx', xlbl='$\lambda$',
          ylbl='Error', fig=fig)
plot.subplot(2, 2, 2)
plot.plot(its.ObjFun, xlbl='Iterations', ylbl='Functional', fig=fig)
plot.subplot(2, 2, 3)
plot.plot(np.vstack((its.PrimalRsdl, its.DualRsdl)).T,
          ptyp='semilogy', xlbl='Iterations', ylbl='Residual',
          lgnd=['Primal', 'Dual'], fig=fig)
plot.subplot(2, 2, 4)
plot.plot(its.Rho, xlbl='Iterations', ylbl='Penalty Parameter', fig=fig)
fig.show()



# Wait for enter on keyboard
input()

[((44, 0, 44, 21), 'numpy.random.seed', 'np.random.seed', ({(44, 15, 44, 20): '(12345)'}, {}), '(12345)', True, 'import numpy as np\n'), ((45, 4, 45, 25), 'numpy.random.randn', 'np.random.randn', ({(45, 20, 45, 21): 'N', (45, 23, 45, 24): 'M'}, {}), '(N, M)', True, 'import numpy as np\n'), ((46, 5, 46, 21), 'numpy.zeros', 'np.zeros', ({(46, 14, 46, 20): '(M, 1)'}, {}), '((M, 1))', True, 'import numpy as np\n'), ((48, 14, 48, 35), 'numpy.random.randn', 'np.random.randn', ({(48, 30, 48, 31): 'L', (48, 33, 48, 34): '1'}, {}), '(L, 1)', True, 'import numpy as np\n'), ((59, 6, 60, 72), 'sporco.admm.bpdn.BPDN.Options', 'bpdn.BPDN.Options', ({(59, 24, 60, 71): "{'Verbose': False, 'MaxMainIter': 500, 'RelStopTol': 0.001, 'AutoRho': {\n    'RsdlTarget': 1.0}}"}, {}), "({'Verbose': False, 'MaxMainIter': 500, 'RelStopTol': \n    0.001, 'AutoRho': {'RsdlTarget': 1.0}})", False, 'from sporco.admm import bpdn\n'), ((76, 7, 76, 28), 'numpy.logspace', 'np.logspace', ({(76, 19, 76, 20): '1', (76, 22, 76, 23): '2', (76, 25, 76, 27): '20'}, {}), '(1, 2, 20)', True, 'import numpy as np\n'), ((77, 25, 77, 59), 'sporco.util.grid_search', 'util.grid_search', ({(77, 42, 77, 49): 'evalerr', (77, 51, 77, 58): '(lrng,)'}, {}), '(evalerr, (lrng,))', False, 'from sporco import util\n'), ((89, 4, 89, 31), 'sporco.admm.bpdn.BPDN', 'bpdn.BPDN', ({(89, 14, 89, 15): 'D', (89, 17, 89, 18): 's', (89, 20, 89, 25): 'lmbda', (89, 27, 89, 30): 'opt'}, {}), '(D, s, lmbda, opt)', False, 'from sporco.admm import bpdn\n'), ((108, 6, 108, 35), 'sporco.plot.figure', 'plot.figure', (), '', False, 'from sporco import plot\n'), ((109, 0, 109, 21), 'sporco.plot.subplot', 'plot.subplot', ({(109, 13, 109, 14): '(2)', (109, 16, 109, 17): '(2)', (109, 19, 109, 20): '(1)'}, {}), '(2, 2, 1)', False, 'from sporco import plot\n'), ((110, 0, 111, 32), 'sporco.plot.plot', 'plot.plot', (), '', False, 'from sporco import plot\n'), ((112, 0, 112, 21), 'sporco.plot.subplot', 'plot.subplot', ({(112, 13, 112, 14): '(2)', (112, 16, 112, 17): '(2)', (112, 19, 112, 20): '(2)'}, {}), '(2, 2, 2)', False, 'from sporco import plot\n'), ((113, 0, 113, 68), 'sporco.plot.plot', 'plot.plot', (), '', False, 'from sporco import plot\n'), ((114, 0, 114, 21), 'sporco.plot.subplot', 'plot.subplot', ({(114, 13, 114, 14): '(2)', (114, 16, 114, 17): '(2)', (114, 19, 114, 20): '(3)'}, {}), '(2, 2, 3)', False, 'from sporco import plot\n'), ((118, 0, 118, 21), 'sporco.plot.subplot', 'plot.subplot', ({(118, 13, 118, 14): '(2)', (118, 16, 118, 17): '(2)', (118, 19, 118, 20): '(4)'}, {}), '(2, 2, 4)', False, 'from sporco import plot\n'), ((119, 0, 119, 72), 'sporco.plot.plot', 'plot.plot', (), '', False, 'from sporco import plot\n'), ((125, 0, 125, 7), 'builtins.input', 'input', ({}, {}), '()', False, 'from builtins import input\n'), ((70, 8, 70, 35), 'sporco.admm.bpdn.BPDN', 'bpdn.BPDN', ({(70, 18, 70, 19): 'D', (70, 21, 70, 22): 's', (70, 24, 70, 29): 'lmbda', (70, 31, 70, 34): 'opt'}, {}), '(D, s, lmbda, opt)', False, 'from sporco.admm import bpdn\n'), ((99, 10, 99, 28), 'numpy.hstack', 'np.hstack', ({(99, 20, 99, 27): '(x0, x)'}, {}), '((x0, x))', True, 'import numpy as np\n'), ((52, 15, 52, 35), 'numpy.random.randn', 'np.random.randn', ({(52, 31, 52, 32): 'N', (52, 33, 52, 34): '(1)'}, {}), '(N, 1)', True, 'import numpy as np\n'), ((72, 18, 72, 30), 'numpy.abs', 'np.abs', ({(72, 25, 72, 29): '(x - x0)'}, {}), '(x - x0)', True, 'import numpy as np\n'), ((115, 10, 115, 51), 'numpy.vstack', 'np.vstack', ({(115, 20, 115, 50): '(its.PrimalRsdl, its.DualRsdl)'}, {}), '((its.PrimalRsdl, its.DualRsdl))', True, 'import numpy as np\n')]