%%%% Success Probability w.r.t band
clear all
close all
display('-------------------------------')

%% settings
N = 100; % the number of rows of sensing matrix A
F = 20; % refinement factor
M = 200 * F; % the number of columns of sensing matrix A
sparsity = 10; % sparsity of x

i = sqrt(-1);
ell = 3 * F; % 3 * rayleigh length
band = 2 * F - 1; % radius of excluded band = 2 * rayleigh length
neighbor = F - 1; % radius of local optimization = 1 * rayleigh length

noise = 10; % noise level; if noise = 5, 5% gaussian noise is added
range = 10; % dynamic range of x
search = 1; % alway use search = 1

%% genearete signal x whose support is separated by at least ell(3*rayleigh length)
%% and two indices are separated exactly by ell
supp = sepsignalclose(M,sparsity,ell-1);
if size(supp,1) == 1 % mask supp be a column vector
supp = supp.';
end
x = zeros(M,1);
x(supp) = linspace(1,range,sparsity);

%% construct the sensing matrix A
%% A are restricted Fourier measurement
A = zeros(N,M);
for k = 1:N
A(k,:) = exp(-i*2*pi*k*(0:M-1)/M)/sqrt(N);
end

%% compute the coherence of the A
%% coherence pattern is abs(A'*A)
% CoMatrix = A' * A;
% coherence = max(max(abs(CoMatrix - eye(M))));
% figure
% imshow(CoMatrix)
% colorbar
% title('pairwise coherence pattern')
% xlabel(['100*4000 matrix with F = 20 & coherence = ' num2str(coherence)])

%% compute measurement
b = A * x;
%%
sigma = (norm(b)/sqrt(N))*(noise/100);
e = normrnd(0,sigma,N,1);
if noise > 0 % ||e|| <= ep
ep = sqrt(N) * sigma;
else
ep = 10^(-12);
end
b = b + e;

% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% OMP
[xrec , supprec , r ] = OMP( A , b , sparsity , ep);
% xrec: recovered signal
% supprec: recovered support
% r = b - A * xrec
if size(supprec,1)~=size(supp,1) % make supprec be a column vector
supprec = supprec.';
end
dis = bottleneck(supp,supprec)/F;
% bottleneck distance between the exact and the recovered support in terms of rayleigh length
err = norm(x(sort(supp,'ascend'))-xrec(sort(supprec,'ascend')))/norm(x(supp));
% relative error
res = norm(r) / norm(b);
% relative residual
figure
stem((1:M)/F,x,'r*')
hold on
stem((1:M)/F,xrec,'bp')
legend('exact','recovered')
xlabel(['dis = ' num2str(dis) ' relative residual = ' num2str(res) ' relative error = ' num2str(err)])
title('OMP')

%% BLOOMP
[xrec , supprec , r ] = BLOOMP(A , b , sparsity , band , neighbor , ep , search);
% xrec: recovered signal
% supprec: recovered support
% r = b - A * xrec
if size(supprec,1)~=size(supp,1) % make supprec be a column vector
supprec = supprec.';
end
dis = bottleneck(supp,supprec)/F;
% bottleneck distance between the exact and the recovered support in terms of rayleigh length
err = norm(x(sort(supp,'ascend'))-xrec(sort(supprec,'ascend')))/norm(x(supp));
% relative error
res = norm(r) / norm(b);
% relative residual
figure
stem((1:M)/F,x,'r*')
hold on
stem((1:M)/F,xrec,'bp')
legend('exact','recovered')
xlabel(['dis = ' num2str(dis) ' relative residual = ' num2str(res) ' relative error = ' num2str(err)])
title('BLOOMP')

%% BP
opts = [];
opts.tol = 10^(-3);
opts.delta = ep;
xrec = yall1(A,b,opts);
% relative residual
res = norm(A*xrec-b)/norm(b);
% relative error
err = norm(x-xrec)/norm(x(supp));
figure
stem((1:M)/F,x,'r*')
hold on
stem((1:M)/F,xrec,'bp')
legend('exact','recovered')
xlabel(['relative residual = ' num2str(res) ' relative error = ' num2str(err)])
title('BP')

%% BP-BLOT
supprec = BT(xrec,sparsity,band);
[xrec , supprec , r] = LO( A , b , neighbor , sparsity , supprec , search , ep );
if size(supprec,1) ~= size(supp,1)
supprec = supprec.';
end
% bottleneck distance between the exact and the recovered support in terms of rayleigh length
dis = bottleneck(supp,supprec)/F;
% relative residual
res = norm(r)/norm(b);
% relative error
err = norm(x(sort(supp,'ascend'))-xrec(sort(supprec,'ascend')))/norm(x(supp));
figure
stem((1:M)/F,x,'r*')
hold on
stem((1:M)/F,xrec,'bp')
legend('exact','recovered')
xlabel(['dis = ' num2str(dis) ' relative residual = ' num2str(res) ' relative error = ' num2str(err)])
title('BP-BLOT')