function [A,R,Asic,Rsic] = IFsic(H,snr)

%inputs:
%H - real channel matrix
%snr - signal-to-noise ratio in linear scale

%outputs
%A - (approximated) optimal target integer-valued matrix for standart IF
%R - Computation rates (in bits) that correspond to each row of A using standard IF
%Asic - (approximated) optimal target integer-valued matrix for successive IF
%Rsic - Computation rates (in bits) that correspond to each row of Asic using successive IF



[N,M]=size(H);
Kee=inv(eye(M)+snr*H'*H);
Lt=chol(Kee);

[almostHt T]=fullLLL(Lt,1);
Atmp=T';
noiseVarTmp=diag(Atmp*Kee*Atmp');
[noiseVar idx]=sort(noiseVarTmp);
A=Atmp(idx,:);
R=-0.5*log2(noiseVar);
idx=find(R<0);
R(idx)=0;


Asic=A;
Al=A;
Kll=Kee;
for l=2:M
    L=chol(Al*Kll*Al');
    Ll=L(2:end,2:end);
    Kll=Ll'*Ll;
    [almostHt T]=fullLLL(Ll,1);
    Al=T';
    Atmp=diag([ones(1,l-1) zeros(1,M+1-l)]);
    Atmp(l:M,l:M)=Al;
    Asic=Atmp*Asic;
end
Lt=chol(Asic*Kee*Asic');
noiseVarSIC=diag(Lt).^2;
Rsic=-0.5*log2(noiseVarSIC);

end

function [Hreduced T]=fullLLL(H,delta)

[r c]=size(H);
[Qs Rs]=qr(H);

Q=Qs;
R=Rs;
T=eye(r);

ctr=1;
k=2;
while k<=c
    for l=[k-1:-1:1]
        
        u=round(R(l,k)/R(l,l));
       
        if u~=0
            R(1:l,k)=R(1:l,k)-u*R(1:l,l);
            T(:,k)=T(:,k)-u*T(:,l);
        end
    end
    
    if (delta*(R(k-1,k-1)^2)>((R(k,k)^2)+(R(k-1,k)^2)))
        tmpR=R(:,k-1);        
        R(:,k-1)=R(:,k);
        R(:,k)=tmpR;
        
        tmpT=T(:,k-1);
        T(:,k-1)=T(:,k);
        T(:,k)=tmpT;
        
        alpha=R(k-1,k-1)/sqrt(R(k-1:k,k-1)'*R(k-1:k,k-1));
        beta=R(k,k-1)/sqrt(R(k-1:k,k-1)'*R(k-1:k,k-1));
        theta=[alpha beta;-beta alpha];
        
        R(k-1:k,k-1:c)=theta*R(k-1:k,k-1:c);
        Q(:,k-1:k)=Q(:,k-1:k)*theta';
        
        k=max(k-1,2);
    else
        k=k+1;
    end
    ctr=ctr+1;
    if ctr>1000*c
        break
    end
end
Hreduced=H*T(1:c,1:c);
end