#!/usr/bin/env julia

"""
A module implementing subgradient and prox-linear algorithms for optimizing
a robust objective for the real-valued blind deconvolution problem.
Additionally, it employs a spectral implementation procedure for provably
obtaining initial points which are sufficiently close to optima.

This code implements the algorithms described in [1].

[1]: Vasileios Charisopoulos, Damek Davis, Mateo Diaz and Dmitriy Drusvyatskiy,
"Composite optimization for robust blind deconvolution". arXiv:1901.01624
"""
module BlindDeconv

	include("Utils.jl")

	using LinearAlgebra
	using Random
	using Statistics
	using Arpack  # required for eigs

	# Symbols and functions to export
	export NOISE, PROXM, MatType
	export BCProb
	export genmat, generate_measurements!
	export subgradient_method, proximal_method
	export gen_problem

	# counter of matrix-vector multiplications
	matvecCount = [0]

	"""
		@enum NOISE indep adversarial

	Distinguish between different types of noise. `indep` will result in
	independent corruptions, while `adversarial` will result in arbitrary
	corruptions that might depend on the true signal.
	"""
	@enum NOISE indep=1 adversarial=2


	"""
		@enum PROXM pogs

	Distinguish between the method used to solve proximal subproblems. `pogs`
	(default) will use a variant of ADMM.
	"""
	@enum PROXM pogs=1

	"""
		@enum MatType gaussian complex_gaussian binary randhadm dethadm ortho pdft

	Distinguish between matrix designs: `gaussian` stands for gaussian random
	matrices, `binary` will result in ``\\ell_{i,j} \\in \\{-1, 0, 1\\}``,
	chosen uniformly at random, `randhadm` will result in an ensemble of
	Hadamard matrices multiplied with random sign patterns, and `dethadm`
	will result in a single, deterministic Hadamard matrix.
	`ortho` will result in an orthonormal matrix, which is generated by taking
	the QR decomposition of a random gaussian matrix ``A = Q R``, and `pdft`
	will result in a partial DFT matrix, whose first ``d`` columns are kept.
	"""
	@enum MatType gaussian complex_gaussian binary randhadm dethadm ortho pdft


	"""
		struct BCProb

	A struct for holding the parameters of a blind deconvolution problem
	instance. A breakdown of the parameters follows:

		L      # left measurement matrix, or a callable implementing it
		R      # right measurement matrix, or a callable implementing it
		LT     # a callable implementing the adjoint operator of L
		RT     # a callable implementing the adjoint operator of R
		y :: Array{<:Number, 1}   # the array of measurements
		w :: Array{<:Number, 1}   # true signal `w`
		x :: Array{<:Number, 1}   # true signal `x`
		w0 :: Array{<:Number, 1}  # initial guess for `w`
		x0 :: Array{<:Number, 1}  # initial guess for `x`
		pfail :: Float64         # corruption level
	"""
	struct BCProb
		L
		R
		LT
		RT
		y :: Array{<:Number, 1}
		w :: Array{<:Number, 1}
		x :: Array{<:Number, 1}
		w0 :: Array{<:Number, 1}
		x0 :: Array{<:Number, 1}
		pfail :: Float64
	end


	"""
		opAx!(r, A, x; nconj=false)

	Compute the matrix-vector multiplication ``r = Ax``, when `A` is either
	given as a matrix or as an operator implementing the multiplication. The
	result is stored in `r`. If `nconj` is set to true, does not conjugate `A`
	before computing the result.
	"""
	function opAx!(r, A, x; nconj=false)
		if isa(A, Array{<:Number, 2})
			copyto!(r, (nconj ? A : conj(A)) * x)
		else
			A(r, x)
		end
		matvecCount[1] += 1
	end


	"""
		opATx!(r, AT, x, tmp; nconj=false)

	Compute the matrix-vector multiplication ``r = A^\\top x``, when `AT` is
	either given as a matrix or as an operator implementing the multiplication
	of the adjoint. The result is stored in `r`. If `nconj` is set to true,
	does not conjugate `A` before computing the result.
	"""
	function opATx!(r, AT, x, tmp; nconj=false)
		if isa(AT, Array{<:Number, 2})
			copyto!(r, (nconj ? AT : conj(AT)) * x)
		else
			AT(r, x, tmp)
		end
		matvecCount[1] += 1
	end


	"""
		residual!(r, Lw, Rx, prob, w, x)

	Computes the residual of the blind deconvolution objective, given the
	two input signals `w` and `x`. Two placeholders for temporary results, `Lw`
	and `Rx`, should also be provided. These will be overwritten as follows:
	``Lw = A_L w, \\; Rx = A_R x``. Additionally, the result is stored in `r`.

	The residual is given by
	``
	\\frac{1}{m} \\sum_{i=1}^m
			\\langle \\ell_i, w \\rangle \\langle x, r_i \\rangle - y_i.
	``

	In the case of complex vectors, the inner product is interpreted as
	``
		\\langle x, y \\rangle := x^H y.
	``
	"""
	function residual!(r, Lw, Rx, prob::BCProb, w, x)
		# Lw = L * w, Rx = R * x
		opAx!(Lw, prob.L, w)
		opAx!(Rx, prob.R, conj(x), nconj=true)
		# r := Lw .* Rx .- prob.y
		broadcast!(*, r, Lw, Rx)
		broadcast!(-, r, r, prob.y)
		rmul!(r, 1 / (length(prob.y)))
	end


	"""
		robust_loss(r, Lw, Rx, prob::BCProb, w, x)

	Computes the robust loss of the blind deconvolution objective, given the
	two input signals `w` and `x`. Two placeholders for temporary results, `Lw`
	and `Rx`, should also be provided. These will be overwritten as follows:
	``Lw = A_L w, \\; Rx = A_R x``. Additionally, the residual is stored in `r`.

	The robust loss is given by
	``
	\\frac{1}{m} \\sum_{i=1}^m \\left|
			\\langle \\ell_i, w \\rangle \\langle r_i, x \\rangle - y_i
			\\right|.
	``
	"""
	function robust_loss(r, Lw, Rx, prob::BCProb, w, x)
		residual!(r, Lw, Rx, prob, w, x)
		return norm(r, 1)
	end


	csign(x) = begin
		rnorm = abs(real(x)) / abs(x); inorm = abs(imag(x)) / abs(x)
		return complex(sign(real(x)) * rnorm, sign(imag(x)) * inorm)
	end

	"""
		subgrad!(gw, gx, r, prob::BCProb, w, x, Lw, Rx; tmpw, tmpx)

	Compute the subgradient of the robust loss in-place, given the problem
	`prob`. `r` is a vector that is assumed to have size equal to `prob.n`.
	`w` and `x` are assumed to be the input signals. Two placeholders for
	temporary results, `Lw` and `Rx`, must be provided. Two more placeholders
	`tmpw, tmpx` for the adjoint operators, ``L^\\top, R^\\top`` should be
	provided when `prob.L, prob.R` are not matrices.

	The subgradient of the robust loss is given by
	``
	\\frac{1}{m} \\cdot
	\\sum_{i=1}^m
		\\mathrm{sign}(
			\\langle \\ell_i, w \\rangle \\cdot
			\\langle r_i, x \\rangle - y_i)
		\\cdot
	\\begin{pmatrix}
		\\langle r_i, x \\rangle \\ell_i^\\top \\
		\\langle \\ell_i, w \\rangle r_i^\\top
	\\end{pmatrix}
	``
	for w and x, respectively. When the measurements are complex, the above
	becomes
	``
	\\frac{1}{2m} \\cdot
	\\sum_{i=1}^m
		\\frac{1}{
			\\left| \\langle \\ell_i, w \\rangle \\cdot
			\\langle x, r_i \\rangle - y_i \\right|
		\\cdot
	\\begin{pmatrix}
		\\left( \\langle \\ell_i, w \\rangle \\cdot
		\\langle x, r_i \\rangle - y_i \\right)
		\\langle r_i, x \\rangle \\ell_i \\\\
		\\overline{\\left( \\langle \\ell_i, w \\rangle \\cdot
		\\langle x, r_i \\rangle - y_i \\right)}
		\\langle \\ell_i, w \\rangle r_i
	\\end{pmatrix}
	``

	Sizes: `Lw`, `Rx` should be `m`. `tmpw`, `tmpx` should have sizes `d1` and
	`d2`, respectively, if provided.
	"""
	function subgrad!(gw, gx, r, prob::BCProb, w, x, Lw, Rx;
					  tmpw=nothing, tmpx=nothing)
		# compute residual and sign first
		residual!(r, Lw, Rx, prob, w, x)
		# if complex, abs() is actually norm of complex number
		if (isa(prob.L, Array{Complex{Float64}, 2}) ||
		    isa(prob.R, Array{Complex{Float64}, 2}))
			map!(csign, r, r)
	   	else
			map!(sign, r, r)  # sign(e_i)
		end
		# Lw = L * w .* r; Rx = R * x .* r => coefficients for subgrads
		broadcast!(*, Rx, conj(Rx), r); broadcast!(*, Lw, Lw, conj(r))
		if isa(prob.L, Array{<:Number, 2})
			copyto!(gw, transpose(prob.L) * Rx)
		else
			prob.LT(gw, Rx, tmpw)
		end
		if isa(prob.R, Array{<:Number, 2})
			copyto!(gx, transpose(prob.R) * Lw)
		else
			prob.RT(gx, Lw, tmpx)
		end
		matvecCount[1] += 2
	end


	"""
		function genmat(m::Int, d::Int, mtype::MatType)

	Generates a measurement matrix for the blind deconvolution problem. If
	`mtype` is `gaussian` or `binary`, it returns the appropriate randomly
	generated matrices. If `mtype` is `randhadm`, it returns two callables
	that correspond to the operators ``x \\to Ax`` and ``x \\to A^\\top x``.
	The first callable accepts two argument, the result (to be modified in
	place) and the input vector, while the second callable accepts an extra
	argument which has to be a temporary vector of length `d`.
	If `mtype` is `dethadm`, it returns two callables, each of which takes
	two arguments and stores the FWHT of the second argument on the first.
	If `mtype` is `ortho`, it generates a Gaussian matrix, computes its QR
	factorization, and returns ``\\sqrt{m} Q``.
	If `mtype` is `pdft`, it generates two callables, of which the first
	corresponds to the operation ``r = \\frac{1}{\\sqrt{m}} F R_{1:d} x``,
	where ``F`` is a DFT matrix and ``R_{1:d}`` is its restriction to the
	first ``d`` columns, and the second corresponds to the adjoint operation.
	"""
	function genmat(m::Int, d::Int, mtype::MatType)
		if mtype == gaussian
			return randn(m, d)
		elseif mtype == complex_gaussian
			return Utils.randn_complex(m, d)
		elseif mtype == binary
			return Utils.binary_matrix(m, d)
		elseif mtype == randhadm
			ispow2(d) || throw(ArgumentError("d must be a power of 2!"))
			k = trunc(Int, m / d)
			s = Utils.signpat(d, k)
			# return a pair of callables
			Aop = ((r, x) -> Utils.A_safe!(r, x, s, k, d))  # Aop => r = Ax
			# ATop => r = A^T x
			ATop = ((r, x, tmp) -> Utils.AT_safe!(r, x, s, k, d, tmp))
			return Aop, ATop
		elseif mtype == dethadm
			ispow2(m) || throw(ArgumentError("m must be a power of 2"))
			temp = fill(0., m)
			Aop = ((r, x) -> Utils.Adet_safe!(r, x))
			ATop = ((r, x, _) -> Utils.ATdet_safe!(r, x, d, temp))
			return Aop, ATop
		elseif mtype == ortho
			gens = randn(m, d); Q, _ = qr(gens)
			return sqrt(m) * Array(Q)  # Q is orthonormal
		elseif mtype == pdft
			Aop = ((r, x) -> Utils.dft_partial!(r, x))
			ATop = ((r, x, _) -> Utils.dftT_partial!(r, x, d))
			return Aop, ATop
		end
	end


	"""
		generate_measurements!(y, L, R, w, x, pfail, noise=nothing)

	Generates a set of measurements which are corrupted by a fraction `pfail` by
	noise. When the noise is `nothing`, a set of i.i.d corruptions is generated.
	`L`, `R` can either be matrices or callables that implement Hadamard
	matrix-vector multiplication. `y` is a placeholder for the result which is
	modified in place.
	"""
	function generate_measurements!(y, L, R, w, x, pfail, noise=nothing)
		opAx!(y, L, w)
		if isa(R, Array{<:Number, 2})
			broadcast!(*, y, y, R * conj(x))
		else
			tempr = zero(y); R(tempr, x)  # store on tempr
			broadcast!(*, y, y, tempr)
		end
		num_corr = trunc.(Int, pfail * length(y))
		if noise == nothing
			noise = randn(num_corr)
			y[randperm(length(y))[1:num_corr]] += noise
		else
			# replace the indices with the noise provided
			y[randperm(length(y))[1:num_corr]] = noise[1:num_corr]
		end
	end


	"""
		corrupt_measurements(prob, pfail, noise, method="additive")

	Corrupts the measurements of a blind deconvolution problem instance, with
	noise provided explicitly as a 1-D array. This noise is either of type
	"additive" or "arbitrary"; in the latter case, the noise vector
	completely replaces a subset of the measurements ``y``.
	"""
	function corrupt_measurements(prob, pfail, noise, method="additive")
		num_corr = trunc.(Int, pfail * length(prob.y))
		if method == "additive"
			prob.y[randperm(length(prob.y))[1:num_corr]] += noise[1:num_corr]
		elseif method == "arbitrary"
			prob.y[randperm(length(prob.y))[1:num_corr]] = noise[1:num_corr]
		else
			throw(ErrorException("Invalid corruption method!"))
		end
	end


	"""
		add_noise(y::Array, noise::Array, pfail::Float, replace=false)

	Add noise to a set of measurements ```y``. `noise` is an array containing
	noise, which is either added to ```y`` or replaces ``y`` at randomly chosen
	indices (when `replace == true`). `pfail` is a parameter that controls
	the level of corruption. The measurements are corrupted in place.
	"""
	function add_noise!(y::Array, noise::Array, pfail::Float64, replace=false)
		m = length(y); ncorr = trunc(Int, pfail * m)
		(length(noise) == m) || throw(ArgumentError(
			"[noise] must have same length as [y]"))
		indices = randperm(m)[1:ncorr]
		if replace == true
			y[indices] = noise[indices]
		else
			y[indices] += noise[indices]
		end
	end


	"""
		direction_init!(wev, xev, y, L, R LT, RT; kw=1, kx=1)

	Computes an initial estimate for the directions of the solution using a
	variant of Duchi's spectral initialization method. `wev, xev` are initial
	estimates of the directions, which are modified in place. `y` is the set
	of measurements, `L, R` are either the measurement matrices or callables
	implementing ``x \\to Lx, \\; x \\to Rx`` respectively, and `LT, RT` are
	the corresponding callables implementing the adjoint operators.
	`kw` and `kx` are the number of sign pattern matrices that were used, in
	the case of random Hadamard ensembles. If left unspecified, they are set
	to `1` by default.
	"""
	function direction_init!(wev, xev, y, L, R, LT=nothing, RT=nothing;
							 kw=1, kx=1)
		med = median(abs.(y))
		ind = abs.(y) .<= med
		if isa(L, Array{<:Number, 2})
			Lind = transpose(L)[:, ind]
			Linit = (1 / length(y)) * Lind * Lind'
			_, wev[:] = eigs(Linit, nev=1, which=:SM)
		else
			Utils.specinit_op!(wev, trunc.(Int, ind), kw, L, LT)
		end
		if isa(R, Array{<:Number, 2})
			Rind = transpose(R)[:, ind]
			Rinit = (1 / length(y)) * Rind * Rind'
			_, xev[:] = eigs(Rinit, nev=1, which=:SM)
		else
			Utils.specinit_op!(xev, trunc.(Int, ind), kx, R, RT)
		end
	end


	"""
		radius_init!(wev, xev, Lw, Rx, y, L, R)

	Computes an estimate for the radius of the true vectors, given two
	directional estimates. The radius estimate is incorporated by modifying
	the estimates in place. `Lw` and `Rx` are placeholder vectors which are
	used to store ``L w_{\\mathrm{ev}}`` and ``R x_{\\mathrm{ev}}``, and must
	be of length ``m``.
	This version of `radius_init` uses an approximate line search and a
	subgradient update.
	"""
	function radius_init!(wev, xev, Lw, Rx, y, L, R)
		m = length(y)
		opAx!(Lw, L, wev); opAx!(Rx, R, xev)
		broadcast!(*, Lw, Lw, Rx)
		f = (r -> (1 / m) * norm(r * Lw - y, 1))
		gradf = (r -> (1 / m) * sum(
			sign.(r * Lw - y) .* Lw))
		eta = 5
		rbest = Utils.fmin_alin(f, gradf, eta, 1.0, 1000, eps=1e-8)
		copyto!(wev, sign(rbest) * sqrt(abs(rbest)) * wev)
		copyto!(xev, sqrt(abs(rbest)) * xev)
	end


	"""
		project2ball!(x, gamma=1.0)

	Projects a vector `x` to the ``\\ell_2``-norm ball, scaled by ``\\gamma``.
	The modification is performed in place.
	"""
	function project2ball!(x, gamma=1.0)
		if norm(x) > gamma
			rmul!(x, gamma / norm(x))
		end
	end


	"""
		subgradient_method(prob, lr, qexp iters; eps, gamma, Lw, Rx, polyak_step)

	Runs the subgradient method for the robust blind deconvolution objective
	using a geometrically decaying rate, ``\\lambda q^k``, where
	``\\lambda > 0, q \\in (0, 1)``. If `gamma == nothing`, the norm constraint
	is not enforced. Additionally, if the normalized frobenius distance falls
	below `eps`, terminates. `eps` is set to 1e-12 by default. If `polyak_step`
	is set to `true`, uses the polyak step with ``\\min f = 0``.

	Returns the two approximate solutions found and an array containing the
	frobenius distance between the vector estimates and the true signals.
	"""
	function subgradient_method(prob, lr, qexp, iters;
		eps=1e-12, gamma=nothing, Lw=nothing, Rx=nothing, polyak_step=false)
		d1, d2 = length(prob.w), length(prob.x)
		wxf = norm(prob.w) * norm(prob.x)
		frob = (wv, xv) -> Utils.frob_opt(wv, xv, prob.w, prob.x) / wxf
		xk, wk = copy(prob.x0), copy(prob.w0)
		m, oups = length(prob.y), fill(0., iters)
		gw, gx, r = zero(wk), zero(xk), zero(prob.y)
		# initialize temporary placeholders
		Lw = (Lw == nothing) ? zero(prob.y) : Lw
		Rx = (Rx == nothing) ? zero(prob.y) : Rx
		tmpw = isa(prob.L, Array{<:Number, 2}) ? nothing : zero(wk)
		tmpx = isa(prob.R, Array{<:Number, 2}) ? nothing : zero(xk)
		rate = lr * qexp
		for i = 1:iters
			subgrad!(gw, gx, r, prob, wk, xk, Lw, Rx; tmpw=tmpw, tmpx=tmpx)
			norm_mag = BLAS.nrm2(d1, gw, 1)^2 + BLAS.nrm2(d2, gx, 1)^2
			if polyak_step   # use polyak step size
				loss = robust_loss(r, Lw, Rx, prob, wk, xk)
				step = loss / norm_mag
				copyto!(wk, wk - m * step * gw)
				copyto!(xk, xk - m * step * gx)
			else
				copyto!(wk, wk - rate * ( gw / sqrt(norm_mag) ))
				copyto!(xk, xk - rate * ( gx / sqrt(norm_mag) ))
			end
			oups[i] = frob(wk, xk)
			if oups[i] <= eps
				oups = oups[1:i]
				break
			end
			rate *= qexp
			if gamma != nothing
				project2ball!(wk, gamma); project2ball!(xk, gamma)
			end
		end
		return wk, xk, oups
	end


	"""
		proximal_method(prob, alpha, k, prox=pogs, gamma;
						get_iters, eps, inner_eps)

	Runs the proximal method on the robust loss for `k` iterations with a prox
	step size of `alpha`. The user can choose between using POGS and Convex.jl
	for the proximal subproblems by choosing `prox=pogs` or `prox=cvx`.
	Returns the two approximate solutions found and an array containing the
	frobenius distance of the vector estimates and the true signals.

	The subproblems solved at every step look like

	``
		\\arg \\min_x f_{x_k}(x) + \\frac{1}{2 \\alpha} \\| x - x_k \\|^2
		s.t. \\| x \\|_2 \\leq \\gamma
	``
	If `gamma` is `nothing`, then the norm constraint is not enforced.
	If `get_iters` is set to true, returns the total number of "inner"
	iterations performed for all subproblems. `inner_eps` can be either an
	absolute number or a callable with a single integer parameter (representing
	the current iteration).
	"""
	function proximal_method(prob::BCProb, alpha, k, prox=pogs, gamma=nothing;
		get_iters=false, eps=1e-12, inner_eps=(j -> min(1e-2, 2.0^(-j))))
		# stop user if Hadamard matrices are used
		(isa(prob.L, Array{Float64, 2}) && isa(prob.R, Array{Float64, 2})) ||
			throw(ArgumentError(
				"Proximal method is currently unavailable for matrices " *
				"implemented as callables."))
		prox_pogs(prob, alpha, k, gamma, get_iters=get_iters, eps=eps,
			inner_eps=inner_eps)
	end


	"""
		prox_pogs(prob::BCProb, alpha, k, gamma; get_iters, eps, inner_eps)

	Run the prox-linear method on the robust loss for `k` iterations with prox
	parameter `alpha`, using direct steps according to the POGS method for every
	subproblem.
	Return the two approximate solutions found and an array containing the
	frobenius distance of the vector estimates and the true signals, as well
	as the total count of "inner" iterations performed, if the `get_iters` flag
	is set to true.
	"""
	function prox_pogs(prob::BCProb, alpha, k, gamma;
		get_iters=false, eps=1e-11, inner_eps=1e-3)
		wxfnorm = norm(prob.w) * norm(prob.x)
		frob = (wv, xv) -> Utils.frob_opt(wv, xv, prob.w, prob.x) / wxfnorm
		xk, wk = copy(prob.x0), copy(prob.w0)
		d, n, m = size(wk)[1], size(xk)[1], size(prob.y)[1]
		dists = fill(0., k); dists[1] = frob(wk, xk)
		initers = 0
		if isa(inner_eps, Number)
			eps_fun = (_ -> inner_eps)
		else
			eps_fun = inner_eps
		end
		for i = 2:k
			wk, xk, itcount = pogs_solve(
				wk, xk, prob, alpha, gamma, eps_fun(i))
			dists[i] = frob(wk, xk)
			initers += itcount
			if dists[i] <= eps
				dists = dists[1:i]
				break  # stop if desired accuracy was reached
			end
		end
		if get_iters
			return wk, xk, dists, initers
		else
			return wk, xk, dists
		end
	end


	"""
		soft_thres(x, alpha)

	Implements the soft thresholding operator with threshold alpha.
	"""
	function soft_thres(x, alpha)
		return sign.(x) .* max.(abs.(x) .- alpha, 0)
	end


	"""
		pogs_solve(wk, xk, prob, alpha, gamma, tol=1e-3)

	Solves the following problem:

	``
	\\mbox{Minimize } \\frac{1}{m} \\left\\| Ax - c \\right\\|_1 +
	\\frac{1}{2a} \\left\\| x - d \\right\\|_2^2 \\
	\\mbox{s.t. } \\| x \\|_2 \\leq \\gamma
	``

	which corresponds to the proximal subproblems of the linearized blind
	deconvolution objective. If `gamma == nothing`, the norm constraint is
	not enforced.
	"""
	function pogs_solve(wk, xk, prob, alpha, gamma, tol=1e-3)
		m = length(prob.y); d = length(wk); n = length(xk)
		Rx = prob.R * xk; Lw = prob.L * wk
		# update counts
		matvecCount[1] += 2
		# matrix A
		Aw = Rx .* prob.L; Ax = Lw .* prob.R; A = hcat(Aw, Ax)
		# offset vector c
		c = prob.y - (Lw .* Rx)
		# variables z, y, initialized to 0
		z = fill(0., length(wk) + length(xk)); y = fill(0., m)
		zold = copy(z); yold = copy(y)
		C = inv(Matrix(1.0I, length(z), length(z)) + A' * A)  # cache inv
		# update counts (take into account matrix structure)
		matvecCount[1] += size(A')[1]
		# dual variables
		l = fill(0., length(z)); nu = fill(0., m);
		rho = (1 / m); fx = rho * alpha; fy = (rho * m)
		# residuals
		res_p = res_d = 0
		# accuracies
		stopn = sqrt(length(z))
		eps_p = (v1, v2) -> tol * (stopn + max(norm(v1), norm(v2)))
		eps_d = (l1, l2) -> tol * (stopn + max(norm(l1), norm(l2)))
		iters = 1
		while true
			zhalf = (fx / (1 + fx)) * (z - l)
			if gamma != nothing
				project2ball!(view(zhalf, 1:d), gamma)   # project w-part
				project2ball!(view(zhalf, (d+1):(d+n)), gamma)  # project x-part
			end
			yhalf = c + soft_thres(y - nu - c, 1 / fy)
			tk = zhalf + l + A' * (yhalf + nu)
			# update counts
			matvecCount[1] += 1
			copyto!(z, C * tk); copyto!(y, A * z)
			# update counts
			matvecCount[1] += 2
			l = l + zhalf - z; nu = nu + yhalf - y
			res_p = norm(vcat(z - zhalf, y - yhalf))
			res_d = rho * (norm(vcat(zold - z, yold - y)))
			if (res_p < eps_p(z, y)) && (res_d < eps_d(l, nu))
				break
			end
			copyto!(zold, z); copyto!(yold, y)
			iters += 1
		end
		return z[1:d] + wk, z[(d+1):end] + xk, iters
	end


	"""
		gen_problem(m::Int, w::Array, x::Array, Ltype::MatType, Rtype::MatType,
					pfail=0.0, adv_noise=false)

	Set up a problem where matrix `L` is of `Ltype` and `R` is of `Rtype`,
	given the true signals `w` and `x`. If `adv_noise` is `true`, impute a
	randomly generated signal into the set of measurements, at indices chosen
	uniformly at random.
	Return the corresponding problem struct, including all the problem data up
	until initialization.
	"""
	function gen_problem(m :: Int, w::Array{<:Number, 1}, x::Array{<:Number, 1},
		Ltype::MatType, Rtype::MatType, pfail=0.0, adv_noise=false)
		Lw = fill(0., m); Rx = fill(0., m)
		d1 = length(w); d2 = length(x); y = fill(0., m)
		if ((Ltype == complex_gaussian) || (Rtype == complex_gaussian) ||
			(Ltype == pdft) || (Rtype == pdft))
			return gen_problem_complex(m, complex(w), complex(x),
									   Ltype, Rtype, pfail)
		end
		L = R = LT = RT = nothing; y = fill(0., m)
		if (Ltype == dethadm) || (Ltype == randhadm)
			L, LT = genmat(m, d1, Ltype)
		else
			L = genmat(m, d1, Ltype)
		end
		if (Rtype == dethadm) || (Rtype == randhadm)
			R, RT = genmat(m, d2, Rtype)
		else
			R = genmat(m, d2, Rtype)
		end
		if adv_noise
			# generate measurements and impute a different signal
			generate_measurements!(y, L, R, w, x, 0.0)
			opAx!(Lw, L, randn(d1))
			opAx!(Rx, R, randn(d2))
			# replace indices
			add_noise!(y, (Lw .* Rx), pfail, true)
		else
			generate_measurements!(y, L, R, w, x, pfail)
		end
		# initialize via Duchi
		w0 = randn(d1); x0 = randn(d2)
		direction_init!(w0, x0, y, L, R, LT, RT, kw=(m / d1), kx=(m / d2))
		radius_init!(w0, x0, Lw, Rx, y, L, R)
		# reset counts
		matvecCount[1] = 0
		return BCProb(L, R, LT, RT, y, w, x, w0, x0, pfail)
	end


	"""
		gen_problem(m::Int, d1::Int, d2::Int, Ltype::MatType, Rtype::MatType,
					pfail::Float64)

	If only the dimensions of the input signals are given instead, generate
	them randomly.
	"""
	function gen_problem(m :: Int, d1::Int, d2::Int,
		Ltype::MatType, Rtype::MatType, pfail=0.0, adv_noise=false)
		w = randn(d1); x = randn(d2)
		return gen_problem(m, w, x, Ltype, Rtype, pfail, adv_noise)
	end


	"""
		gen_problem(w::Array, x::Array, L::Array, R::Array, pfail=0.0,
					adv_noise=false)

	Set up a problem where the user provides the measurement vectors explicitly,
	in addition to the true signals.
	"""
	function gen_problem(w::Array, x::Array, L::Array, R::Array,
						 pfail=0.0, adv_noise=false)
		m = size(L)[1]; d1 = length(w); d2 = length(x);
		if (m != size(R)[1]) || (size(L)[2] != d1) || (size(R)[2] != d2)
			throw(ArgumentError("Matrix sizes are incompatible!"));
		end
		Lw = fill(0.0, m); Rx = fill(0.0, m); y = fill(0.0, m)
		LT = RT = nothing
		if adv_noise
			# generate measurements and impute a different signal
			generate_measurements!(y, L, R, w, x, 0.0)
			opAx!(Lw, L, randn(d1))
			opAx!(Rx, R, randn(d2))
			# replace indices
			add_noise!(y, (Lw .* Rx), pfail, true)
		else
			generate_measurements!(y, L, R, w, x, pfail)
		end
		# initialize via Duchi
		w0 = randn(d1); x0 = randn(d2)
		direction_init!(w0, x0, y, L, R, LT, RT, kw=(m / d1), kx=(m / d2))
		radius_init!(w0, x0, Lw, Rx, y, L, R)
		# reset counts
		matvecCount[1] = 0
		return BCProb(L, R, LT, RT, y, w, x, w0, x0, pfail)
	end


	"""
		gen_problem_complex(m::Int, w::Array, x::Array, Ltype::MatType,
							Rtype::MatType, pfail::Float64)

	Set up a blind deconvolution problem involving complex measurements, given
	the true signals `w` and `x`, matrix types `Ltype` and `Rtype`, and the
	corruption level `pfail`. If `adv_noise` is `true`, impute a randomly
	generated signal into the set of measurements, at indices chosen
	uniformly at random.
	Return a struct containing all generated problem data up until
	initialization.
	"""
	function gen_problem_complex(m :: Int, w :: Array, x :: Array,
		Ltype::MatType, Rtype::MatType, pfail=0.0, adv_noise=false)
		d1 = length(w); d2 = length(x)
		Lw = fill(0.0 + 0.0j, m); Rx = fill(0.0 + 0.0j, m)
		L = R = LT = RT = nothing; y = fill(0.0 + 0.0j, m)
		if (Ltype == pdft)
			L, LT = genmat(m, d1, Ltype)
		else
			L = genmat(m, d1, Ltype)
		end
		if (Rtype == pdft)
			R, RT = genmat(m, d2, Rtype)
		else
			R = genmat(m, d2, Rtype)
		end
		if adv_noise
			# generate measurements and impute a different signal
			generate_measurements!(y, L, R, w, x, 0.0)
			opAx!(Lw, L, Utils.randn_complex(d1))
			opAx!(Rx, R, Utils.randn_complex(d2), nconj=true)
			# replace indices
			add_noise!(y, (Lw .* Rx), pfail, true)
		else
			generate_measurements!(y, L, R, w, x, pfail)
		end
		# initialize via Duchi
		w0 = Utils.randn_complex(d1); x0 = Utils.randn_complex(d2)
		direction_init!(w0, x0, y, L, R, LT, RT, kw=(m / d1), kx=(m / d2))
		radius_init!(w0, x0, Lw, Rx, y, L, R)
		# reset matvecCount
		matvecCount[1] = 0
		return BCProb(L, R, LT, RT, y, w, x, w0, x0, pfail)
	end


	"""
		gen_problem_complex(m::Int, d1::Int, d2::Int, Ltype::MatType,
							Rtype::MatType, pfail=0.0, adv_noise=false)

	If only the dimensions of the true signals are given, generate them at
	random by sampling complex standard normal variables.
	"""
	function gen_problem_complex(m::Int, d1::Int, d2::Int, Ltype::MatType,
								 Rtype::MatType, pfail=0.0, adv_noise=false)
		w = Utils.randn_complex(d1); x = Utils.randn_complex(d2)
		return gen_problem_complex(m, w, x, Ltype, Rtype, pfail, adv_noise)
	end


	"""
		clearMatvecCount()

	Resets the counter of matrix-vector products.
	"""
	function clearMatvecCount()
		matvecCount[1] = 0;
	end

	"""
		getMatvecCount()

	Return the counter of matrix-vector products.
	"""
	function getMatvecCount()
		return matvecCount[1]
	end
end
