% Copyright (C) 2016 by Deanna Needell, Rayan Saab, and Tina Woolf

% This file is part of ''Weighted L1-Minimization for Sparse Recovery under
%                        Arbitrary Prior Information''

% 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/>.

% compareSufficientRIPBound.m

%% Cleanup and Setup
clearvars;
close all;
clc;
a        = 3;    % RIP constant parameter
saveFigs = 0;    % flag to save figures
w_large  = 0.5;  % largest weight
w_mid    = 0.4;  % middle weight
w_small  = 0.25; % smallest weight

%% Single Support Estimate
alpha = 1; % support accuracy
rho   = 1; % support size = rho*s

% RIP bound when using largest weight
b_large     = w_large + (1-w_large)*sqrt(1+rho-2*alpha*rho);
Delta_large = (a-b_large^2)/(a+b_large^2);
% RIP bound when using middle weight
b_mid       = w_mid + (1-w_mid)*sqrt(1+rho-2*alpha*rho);
Delta_mid   = (a-b_mid^2)/(a+b_mid^2);
% RIP bound when using smallest weight
b_small     = w_small + (1-w_small)*sqrt(1+rho-2*alpha*rho);
Delta_small = (a-b_small^2)/(a+b_small^2);

%% Two Support Estimates
w1     = w_large; % this weight is used with rho1
w2     = w_small; % this weight is used with rho2
alpha1 = alpha;   % accuracy of set T1
alpha2 = alpha;   % accuracy of set T2
rho1   = linspace(0,rho,100);

% Determine RIP upper bound as rho1 and rho2 vary
% Assuming rho1+rho2 = rho, so knowning rho1 determines rho2
for pp = 1:length(rho1)
    
    rho2(pp) = rho-rho1(pp);

    % Result of Needell, Saab, and Woolf
    K = w2 + (1-w1)*sqrt(1+rho1(pp)-2*alpha1*rho1(pp)+rho2(pp)-2*alpha2*rho2(pp)) ...
       + (w1-w2)*sqrt(1+rho2(pp)-2*alpha2*rho2(pp)); 

    Deltas2(pp) = (a-K^2)/(a+K^2);
    
    % Result of Mansour and Yilmaz, "Weighted-L1 minimization with multiple
    % weighting sets," Proc. SPIE, Wavelets and Sparsity, 2011.
    K_M = w1+w2-1 + (1-w1)*sqrt(1+rho1(pp)-2*alpha1*rho1(pp)) + (1-w2)*sqrt(1+rho2(pp)-2*alpha2*rho2(pp));
    
    Deltas2_M(pp) = (a-K_M^2)/(a+K_M^2);
    
end
figure; hold all;
plot(rho1,Deltas2,'-.','linewidth',2);
plot(rho1,Deltas2_M,'-r','linewidth',4);
plot(rho1,ones(1,length(rho1))*Delta_large,'-k','linewidth',2);
plot(rho1,ones(1,length(rho1))*Delta_small,'--k','linewidth',2);
hl = legend('Double Support Estimate (proposed)','Double Support Estimate (prior work)',sprintf('Single Support Estimate w=%1.2f',w_large),sprintf('Single Support Estimate w=%1.2f',w_small),1);
xlabel('\rho_1','fontsize',20);
ylabel('Sufficient Restricted Isometry Constant','fontsize',20);
%title({['Restricted Isometry Constant Bound using with a=' num2str(a) ', rho=' num2str(rho) ', alpha=' num2str(alpha)];['w1=' num2str(w1) ', alpha1= ' num2str(alpha1) ', w2=' num2str(w2) ', alpha2=' num2str(alpha2)]},'fontsize',14);
set(gca,'fontsize',20);
set(hl,'fontsize',14);
set(gcf,'PaperPositionMode','Auto'); 
set(gcf,'WindowStyle','normal')
grid on;
if (saveFigs)
    filename = sprintf('compareRIPconstants_twoWeights_%da_%1.1frho_%1.1falpha_%1.1falpha1_%1.1falpha2_%1.2fw1_%1.2fw2',a,rho,alpha,alpha1,alpha2,w1,w2);
    print('-depsc','-r0',[filename '.eps']);
    saveas(gcf,[filename '.png'],'png');
    saveas(gcf,[filename '.fig'],'fig');
end

%% Three Support Estimate
w1     = w_large;    % this weight is used with rho1
w2     = w_mid;      % this weight is used with rho2
w3     = w_small;    % this weight is used with rho3
alpha1 = alpha;      % accuracy of set T1
alpha2 = alpha;      % accuracy of set T2
alpha3 = alpha;      % accuracy of set T3
mystep = 0.003;      % stepsize in rho for plot
rho1   = 0:mystep:1;

% Determine RIP upper bound as rho1, rho2, and rho3 vary
% Assuming rho1+rho2+rho3 = 1, so knowning rho1 and rho2 determines rho3
Deltas3 = nan(length(rho1));
for pp = 1:length(rho1)
    
    % At most, rho2 can be 1-rho1
    rho2 = 0:mystep:1-rho1(pp);
    
    for qq = 1:length(rho2)
        
        rho3 = 1-rho1(pp)-rho2(qq);

        K = w3 + (1-w1)*sqrt(1+rho1(pp)-2*alpha1*rho1(pp)+rho2(qq)-2*alpha2*rho2(qq)+rho3-2*alpha3*rho3) ...
            + (w1-w2)*sqrt(1+rho2(qq)-2*alpha2*rho2(qq)+rho3-2*alpha3*rho3) ...
            + (w2-w3)*sqrt(1+rho3-2*alpha3*rho3);
        
        Deltas3(pp,qq) = (a-K^2)/(a+K^2);
        RHO1(pp,qq) = rho1(pp);
        RHO2(pp,qq) = rho2(qq);
    
    end

end

% Want to determine when the 3 weight RIP bound matches the single weight RIP bounds
mytol     = 1e-4; % tolerance for matching
idx_large = find(abs(real(Deltas3) - Delta_large) < mytol); % indices that match Delta_large
idx_mid   = find(abs(real(Deltas3) - Delta_mid)   < mytol); % indices that match Delta_mid
idx_small = find(abs(real(Deltas3) - Delta_small) < mytol); % indices that match Delta_small

figure; hold on;
plot3(RHO1(idx_large),RHO2(idx_large),1*ones(size(idx_large,1),1),'kd','linewidth',3);
plot3(RHO1(idx_mid),  RHO2(idx_mid),  1*ones(size(idx_mid,1),1), '-k.','linewidth',3);
plot3(RHO1(idx_small),RHO2(idx_small),1*ones(size(idx_small,1),1),'ko','linewidth',3);
surface(RHO1,RHO2,real(Deltas3),'edgecolor','none');
xlabel('\rho_1','fontsize',20); 
ylabel('\rho_2','fontsize',20); 
set(gca,'fontsize',20);
hc = colorbar;
set(hc,'fontsize',16);
%title({['Restricted Isometry Constant Bound with a=' num2str(a) ', rho=' num2str(rho) ', alpha=' num2str(alpha)];['w1=' num2str(w1) ', alpha1= ' num2str(alpha1) ', w2=' num2str(w2) ', alpha2= ' num2str(alpha2) ', w3=' num2str(w3) ', alpha3= ' num2str(alpha3)]},'fontsize',14);
xlim([-0.1 1.1]);
ylim([-0.1 1.1]);
hl = legend(sprintf('Single Support Estimate w=%1.2f',w_large),sprintf('Single Support Estimate w=%1.2f',w_mid),sprintf('Single Support Estimate w=%1.2f',w_small));
set(hl,'fontsize',14);
set(gcf,'PaperPositionMode','Auto'); 
set(gcf,'WindowStyle','normal')
if (saveFigs)
    filename = sprintf('compareRIPconstants_threeWeights_%da_%1.1frho_%1.1falpha_%1.1falpha1_%1.1falpha2_%1.1falpha3_%1.2fw1_%1.2fw2_%1.2fw3',a,rho,alpha,alpha1,alpha2,alpha3,w1,w2,w3);
    print('-depsc','-painters','-r600',[filename '.eps']);
    saveas(gcf,[filename '.png'],'png');
    saveas(gcf,[filename '.fig'],'fig');
end