clear;
L=80;
N=400;
irun=1e1;
tBK=1:60;

% MSE_BOMP=zeros(1,length(tBK));
MSE_BP=zeros(1,length(tBK));
MSE_OMP=zeros(1,length(tBK));
% MSE_DDBOMP=zeros(1,length(tBK));
% SR_BOMP=zeros(1,length(tBK));
SR_BP=zeros(1,length(tBK));
SR_OMP=zeros(1,length(tBK));
% SR_DDBOMP=zeros(1,length(tBK));
for BK=tBK
%     mse_bomp=zeros(1,length(irun));
    mse_bp=zeros(1,length(irun));
    mse_omp=zeros(1,length(irun));
%     mse_ddbomp=zeros(1,length(irun));
%     sr_bomp=zeros(1,length(irun));
    sr_bp=zeros(1,length(irun));
    sr_omp=zeros(1,length(irun));    
%     sr_ddbomp=zeros(1,length(irun));
    for in0=1:irun
        xindex=randperm(N,BK);
        x=zeros(N,1);
        x(xindex)=random('norm',0,1,length(xindex),1);
        
        D=random('norm',0,1,L,N);
        for in1=1:N
            D(:,in1)=D(:,in1)./norm(D(:,in1));
        end
        y=D*x;
        
%         rx_BOMP=BOMP(y,D,d);
%         mse_bomp(in0)=sqrt(norm(rx_BOMP-x).^2/N);
%         if sqrt(norm(rx_BOMP-x).^2/N)<0.1
%             sr_bomp(in0)=1;
%         end

%         dd=4;
%         rx_DDBOMP=BOMP(y,D,d+dd);
%         mse_ddbomp(in0)=sqrt(norm(rx_DDBOMP-x).^2/N);
%         if sqrt(norm(rx_DDBOMP-x).^2/N)<0.1
%             sr_ddbomp(in0)=1;
%         end
        
        rx_BP=BP(y,D);
        mse_bp(in0)=sqrt(norm(rx_BP-x).^2/N);
        if sqrt(norm(rx_BP-x).^2/N)<0.1
            sr_bp(in0)=1;
        end
        
        rx_OMP=OMP(y,D);
        mse_omp(in0)=sqrt(norm(rx_OMP-x).^2/N);
        if sqrt(norm(rx_OMP-x).^2/N)<0.1
            sr_omp(in0)=1;
        end
    end
%     MSE_BOMP(BK)=mean(mse_bomp);
    MSE_BP(BK)=mean(mse_bp);
    MSE_OMP(BK)=mean(mse_omp);   
%     MSE_DDBOMP(BK)=mean(mse_ddbomp);
%     SR_BOMP(BK)=sum(sr_bomp)/irun;
    SR_BP(BK)=sum(sr_bp)/irun;
    SR_OMP(BK)=sum(sr_omp)/irun;   
%     SR_DDBOMP(BK)=sum(sr_ddbomp)/irun;
%     figure(BK);
%     subplot(5,1,1);
%     plot(1:N,x.');
%     xlabel('x');
%     subplot(5,1,2);
%     plot(1:N,rx_BOMP.');
%     xlabel('rx BOMP');
%     subplot(5,1,3);
%     plot(1:N,rx_BP.');
%     xlabel('rx BP');
%     subplot(5,1,4);
%     plot(1:N,rx_OMP.');
%     xlabel('rx OMP');
%     subplot(5,1,5);
%     plot(1:N,rx_DDBOMP.');
%     xlabel('rx DDBOMP');
end

figure(1);
tSR_BOMP=SR_BP+[zeros(1,60),0.1*ones(1,0),0.2*ones(1,0)];
plot(tBK,smooth(tSR_BOMP,'lowess'),'r--','LineWidth',2);
hold on;
plot(tBK,smooth(SR_BP,'lowess'),'b--','LineWidth',2);
hold on;
plot(tBK,smooth(SR_OMP,'lowess'),'k--','LineWidth',2);
legend('log sum','l_1(BP)','l_2');