%DEMO_DEQUANTIZATION Dequantization demo
%   This demo shows how a quantized signal, sparse in the DCT domain, can be dequantized
%   solving a convex problem using Douglas-Rachford algorithm
%
%   Suppose signal y has been quantized. In this demo we use quantization levels 
%   that are uniformly spread between the min. and max. value of the
%   signal. The resulting signal is y_Q.
%
%   The problem can be expressed as 
%
%        argmin_x  || x ||_1     s.t.   ||Dx - y_Q||_infty <= alpha/2 
%
%   where D is the synthesis dictionary (DCT in our case) and alpha is the
%   distance between quantization levels. The constraint basically
%   represents the fact that the reconstructed signal samples must stay
%   within the corresponding quantization stripes.
%
%   After sparse coordinates are found, the dequantized signal
%   is obtained simply by synthesis with the dictionary.   
%
%   The program is solved using Douglas-Rachford algorithm. We set 
%
%    f_1(x)=||x||_{1}. Its respective prox is the soft thresholding operator.
%
%    f_2(x)=i_C is the indicator function of the set C, defined as
%
%        C = { x | ||Dx - y_Q||_infty <= alpha/2 } 
%
%   Its prox is the orthogonal projection onto that set, which is realized
%   by entry-wise 1D projections onto the quantization stripes. This is
%   realized for all the entries at once by function proj_box.
%
%   As an alternative, setting algorithm = 'LP' switches to computing the
%   result via linear programming (requires Matlab optimization toolbox).
%
%   Results
%   -------
%
%   Figure 1: Original, quantized and dequantized signals
%
%       
%
%   Figure 2: Quantization error and error of reconstruction (i.e. original - reconstr.)
%
%      
%
%   Figure 3: Coefficients of original and reconstructed signals
%
%      
%
%   References:
%     P. Combettes and J. Pesquet. A douglas--rachford splitting approach to
%     nonsmooth convex variational signal recovery. Selected Topics in Signal
%     Processing, IEEE Journal of, 1(4):564--574, 2007.
%     
%
%   Url: https://epfl-lts2.github.io/unlocbox-html/doc/demos/demo_dequantization.php

% Copyright (C) 2012-2016 Nathanael Perraudin.
% This file is part of UNLOCBOX version 1.7.4
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Authors: Pavel Rajmic, Pavel Záviška
% Date: August 2015



clc
clear
close all

%% Parameters

N = 64; %length of the signal
k = 5; %sparsity of the signal in the dictionary
d = 10; %number of quantization levels

% 'DR' ... Douglas-Rachford algorithm
% 'LP' ... linear programming (requires Matlab optimization toolbox)
algorithm =  'DR' %'LP' %
                    


%% Dictionary and the original sparse signal
% A = dctmtx(N)';  % dictionary in which the signal will be sparse = inverse DCT
A = @(x) idct(x);  % dictionary in which the signal will be sparse = inverse DCT
At = @(x) dct(x); %its inverse (= transpose)

if k>N
    error('Sparsity cannot be greater then signal length/number of atoms');
end

% generate random support
support = randperm(N);
support = support(1:k);

% Determine the coefficient values
x_support =  1 + 3*rand([k 1]);
x_support =  x_support .* (((rand([k 1])>.5)*2)-1); %randomizing signs
x = zeros(N,1);
x(support) = x_support; %complete vector including zero coefs

% synthesize signal
% y = A(:, support)*x_support;
y = A(x);

%or just load a fixed dataset
% clear
% load('dequantization_dataset_01')

% show coefficients
fig_coef = figure;
h_coefs_orig = bar(x);
hold on;
title('Original coefficients of sparse signal');



%% Quantize signal
%Determine the max and min value of the signal. The quantization
%levels are spread uniformly between them
min_y = min(min(y));
max_y = max(max(y));
range = max_y - min_y;

quant_step   = range / (d-1); %quantization step
dec_bounds = (min_y+quant_step/2) : quant_step : (max_y-quant_step/2);  %decision boundaries lie in the middle
quant_levels = min_y : quant_step : max_y; %quantization levels

%use Matlab function to quantize signal
[index, y_quant] = quantiz(y, dec_bounds, quant_levels);
y_quant = y_quant';

%constraints for the signal samples;
%quantized signal lies on the quant. levels and the true signal cannot be outside the boundaries
lower_dec_bounds = y_quant - (quant_step/2);
upper_dec_bounds = y_quant + (quant_step/2);

%the highest and lowest boundary coincide precisely with the quantization levels!:
min_quant_level = quant_levels(1);
max_quant_level = quant_levels(end);
upper_dec_bounds(upper_dec_bounds > max_quant_level) = max_quant_level;
lower_dec_bounds(lower_dec_bounds < min_quant_level) = min_quant_level;


%% Show original vs. quantized signal
%define colors
grey = 0.6* ones(1,3);
lightgrey = 0.8* ones(1,3);
black = [0,0,0];
blue = [0.251,0.584,0.808];
orange = [0.851,0.325,0.098];
green = [0 1 0];

%plot of the signals
fig_time = figure;
h_orig = plot(y, '.-', 'Color', blue);
hold on;
h_quant = plot(y_quant, '.-', 'Color', orange);
title('Original and quantized signal');

%plot signal constraints
h_signal_constr = plot(upper_dec_bounds, 'Color', lightgrey);
plot(lower_dec_bounds, 'Color', lightgrey);

for j=1:d %quantization levels
    yPos = quant_levels(j);
    h_q_lev = plot(ones(1,N) * yPos, 'Color', grey);
end

for j=1:(d-1) %decision boundaries
    yPos = dec_bounds(j);
    h_dec_b = plot(ones(1,N) * yPos, ':', 'Color', grey); 
end

axis tight

%legend
uistack(h_orig,'top');
uistack(h_quant,'top');
h_legend_time = legend([h_orig, h_quant, h_q_lev, h_dec_b, h_signal_constr ], 'original', 'quantized', 'quantiz. levels', 'decision bounds', 'signal constraints');

%quanization noise
fig_quant_noise = figure;
quant_noise = y - y_quant;
h_quant_noise = plot(quant_noise, 'Color', blue);
hold on;
title('Quantization noise');


%% Sparse dequantization
switch algorithm
    case 'LP'  %dequantization using linear programming (doubles the number of variables)
        f = (ones([1 2*N]));
        
        b = [-lower_dec_bounds; upper_dec_bounds];
        Amatrix = A(eye(N)); %generate explicit dictionary matrix
        % A_ = [A -A; -A A];
        A_ = [Amatrix -Amatrix; -Amatrix Amatrix];
        lb = zeros(2*N,1);  %all variables must be nonnegative
        
        w = linprog(f,A_,b,[],[],lb);  %l1-minimization via lin. program
        
        uv = reshape(w,N,2);      %split w into two
        u = uv(:, 1);
        v = uv(:, 2);
        x_reconstructed = v - u;    %sparse vector is determined by subtracting the two non-negative
        y_dequant = A(x_reconstructed);    %reconstruct signal
        
        sol = x_reconstructed;
        
    case 'DR'  % Douglas-Rachford
        
        param.lower_lim = lower_dec_bounds;
        param.upper_lim = upper_dec_bounds;
        
        indi_thr = 10e-5; % threshold for identifying the point to lie in the set
        
        f1.eval = @(x) norm(x,1);
        f1.prox = @(x,T) sign(x).*max(abs(x)-T, 0);
        
        f2.eval = @(x) 1 / (1 - ( any((A(x(:))-param.upper_lim)>indi_thr)) || any((A(x(:))-param.lower_lim) < -indi_thr) ) - 1; %zero if x is inside boundaries, Inf otherwise
        f2.prox = @(x,T) At(proj_box(A(x),[],param)); %box projection in the signal space (thanks to DCT being orthogonal)
        
        %%%%%%%% UNLOCBOX version %%%%%%%%%%%%%
        % setting different parameter for the simulation
        paramsolver.verbose = 5;  % display parameter
        paramsolver.maxit = 300;        % maximum iteration
        paramsolver.tol = 10e-7;        % tolerance to stop iterating
        paramsolver.lambda = 1;   % step for DR algorithm (default 1)
        paramsolver.gamma = 1e-2;        % here threshold for soft thresholding
        
        [sol, info] = douglas_rachford(At(y_quant), f1, f2, paramsolver);
        info
        
        sol = f2.prox(sol,[]); %final projection into the constraints
        y_dequant = A(sol);
        
%         pause
%         
%         %%%%%%%% manual version %%%%%%%%%%%%%      
%         %starting point
%         DR_y = A(y_quant);
%         DR_x_old = DR_y;
%         
%         relat_change_coefs = 1;
%         relat_change_obj = 1;
%         cnt = 1; %iteration counter
%         obj_eval = [];
%         
% %         while relat_change_coefs > 0.00001
%         while relat_change_obj > paramsolver.tol
%             % DR: gamma = 1
%             DR_x = f2.prox(DR_y,[]);
%             obj_eval = [obj_eval, f1.eval(DR_x) + f2.eval(DR_x)]; %record values of objective function
%             DR_y = DR_y + paramsolver.lambda*(f1.prox(2*DR_x-DR_y, paramsolver.gamma)-DR_x);
%             if cnt > 1
%                 relat_change_coefs = norm(DR_x-DR_x_old) / norm(DR_x_old);
%                 relat_change_obj = norm(obj_eval(end) - obj_eval(end-1)) / norm(obj_eval(end-1));
%                 if paramsolver.verbose > 1
%                     fprintf('  relative change in coefficients: %e \n', relat_change_coefs);
%                     fprintf('  relative change in objective fun: %e \n', relat_change_obj);
%                     fprintf('\n');
%                 end
%             end
%             DR_x_old = DR_x;
%             cnt = cnt + 1;
%             
%         end
%         
%         DR_x = f2.prox(DR_y); %final projection into the constraints
%         y_dequant = A(DR_x); %dequantized signal
%         
%         disp(['Finished after ' num2str(cnt) ' iterations.'])
%         
%         
%         %compare UNLOCBOX with manual
%         figure
%         plot([y_dequant A(sol)])
%         norm(y_dequant - A(sol))
%         title('UNLOCBOX vs. manual solution')
%         
%         %behaviour of objective through iterations
%         figure
%         plot(obj_eval)
%         title('Objective function value (after projection into constraints in each iteration)')
    
end




%% Show results

figure(fig_time)
h_dequant = plot(y_dequant, '.-', 'Color', green);
uistack(h_dequant,'top');
delete(h_legend_time)
h_legend_time = legend([h_orig, h_quant, h_dequant, h_q_lev, h_dec_b, h_signal_constr ], 'original', 'quantized', 'dequantized', 'quantiz. levels', 'decision bounds', 'signal constraints');
title('Original, quantized and dequantized signals')

%quantization and reconstruction errors
figure(fig_quant_noise)
h_dequant_error = plot(y-y_dequant, 'Color', green);
title('Quantization error and error of reconstruction (i.e. original - reconstr.)');
axis tight
legend([h_quant_noise h_dequant_error], 'Quantizat. error', 'Error of reconstr.')

%coefficients of reconstructed signal
figure(fig_coef)
hold on
h_coefs_dequant = bar(sol,'FaceColor',green);
title('Coefficients of original and reconstructed signals');
legend([h_coefs_orig h_coefs_dequant], 'Coefs of orig. signal', 'Coefs of dequant. signal')
axis tight


figure(fig_time)
axis tight
