%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Peeling Decoder for Graph Based Codes"
% Written for GNU Octave
%
% INSTITUTE OF TELECOMMUNICATIONS
% University of Stuttgart
% www.inue-uni-stuttgart.de
% author: Ahmed Elkelesh
% date: 09.11.2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [] = decoder(v,c,Hid, peid,hr1,uploadfile1, filename)
% The first function to be called in this slide
% Inputs by the user from the webdemo page: v,c,Hid, peid,hr1,uploadfile1
[countFig1] = mainFunction (v,c,Hid, peid,hr1,uploadfile1, filename);
end

function [countFig1] = mainFunction (v,c,Hid, peid,hr1,uploadfile1, filename)
% Inputs:
% v: variable node degree
% c: check node degree
% Hid: the id of the H matrix to be used
% peid is the probability of erasure in the BEC
% hr1: upload H-matrix or not
% countFig1: the figures counter
countFig=0; % number of figures 
% Defining the Parity Check Matrix H
if hr1 == 1
try
		H = load(uploadfile1);
catch
		decoder(v,c,Hid, peid,0,uploadfile1, filename);
end 
else
if Hid==1
    % H for a Hamming code (n=7, k=4)
    H=[1,0,1,1,1,0,0;1,1,1,0,0,1,0;0,1,1,1,0,0,1;];
end
if Hid==2
    % H for a Hamming code (n=15, k=11)
    H=[0,0,0,0,0,0,0,1,1,1,1,1,1,1,1;0,0,0,1,1,1,1,0,0,0,0,1,1,1,1;0,1,1,0,0,1,1,0,0,1,1,0,0,1,1;1,0,1,0,1,0,1,0,1,0,1,0,1,0,1;];
end
if Hid==3
    % H for a regular LDPC code (v=3, c=4, n=12, k=3)
    H=[1,0,0,0,0,1,0,1,0,1,0,0;1,0,0,1,1,0,0,0,0,0,1,0;0,1,0,0,1,0,1,0,1,0,0,0;0,0,1,0,0,1,0,0,0,0,1,1;0,0,1,0,0,0,1,1,0,0,0,1;0,1,0,0,1,0,0,0,1,0,1,0;1,0,0,1,0,0,1,0,0,1,0,0;0,1,0,0,0,1,0,1,0,1,0,0;0,0,1,1,0,0,0,0,1,0,0,1;];
end
if Hid==4
    % H for a regular LDPC code (v=3, c=6, n=20, k=10)
    H=[0,0,0,0,1,0,0,0,1,1,1,0,0,0,0,1,0,0,0,1;0,0,0,0,0,0,1,1,0,0,1,1,0,1,0,1,0,0,0,0;0,1,1,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1;0,0,0,0,0,1,0,1,0,1,0,0,0,0,0,0,1,1,1,0;1,1,0,0,1,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0;0,0,0,0,0,0,1,0,0,0,1,1,0,1,1,0,0,0,0,1;0,0,0,1,1,1,0,1,0,0,0,0,1,0,1,0,0,0,0,0;1,0,1,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,0,0;1,1,1,1,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0;0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,1,1,0;];
end
if Hid==5
    % H for a regular LDPC code with given degree profiles
    H = generateH(v,c);
end
end
[y, n]=size(H); % n is the codeword length and y is the number of parity bits
k=n-y; % The number of information bits

p_erasure=peid;

A_1=zeros(1,n); % The all zero codeword is going to be used (It is always a codeword)(The aim here is to implement the decoding algorithm not the encoding)
after_bec=bec(A_1,p_erasure); % Passing the random bits through the BEC. The output may be a codeword or not
hFrame=figure();
set(hFrame, "visible", "off");
hold all
title('Initial Graph','fontsize',14);
draw_part_graph(after_bec, H, 1:n, countFig, filename, hFrame); % Drawing the Bipartite Graph (Tanner Graph) of the code using the H matrix
% legend('Variable node','Check node');
countFig=countFig+1; % A new figure so increment the counter
[after_decoding, flag_2, countFig2,count]=peeling_decoder(after_bec,H,k,n,countFig, filename); % Using the Peeling decoder concept in the decoding process.
countFig=countFig2;
if flag_2==1 
title('No edges from erased bits (All erased bits corrected)','fontsize',14);
else
text(1,1.8,'The decoder is STUCK!','Color','g','fontsize',20);
if count == 0
title('Graph including only edges from erased bits','fontsize',14);
else
title(sprintf('After Iteration number %d',count),'fontsize',14);
end
end
countFig1=countFig;
for jh= countFig:9
new_filename = [filename '_' num2str(jh)];
saveas(hFrame,new_filename,'png');
end
end

function [out] = bec (in,erasure_probability)
% This is a function that simulates the binary erasure channel with a given
% probability of erasure.
% The input sequence contains 0s and 1s.
% If an erasure occurs, then it is represented by -1.
% The output sequence now contains 0s, 1s and -1s.
out=in;
out(rand(size(in))<erasure_probability)=-1;
end

function [out, flag_2, countFig2,count] = peeling_decoder (x,H,k,n,countFig, filename)
% This is the function performing the iterative peeling decoder algorithm
% The inputs to this function are:
% x is the received bit sequence (may be erroneous codeword or not)
% H the parity check matrix of the code
% k the number of information bits
% n the number of bits per codeword
% countFig number of figures
% The outputs of this function are:
% out the corrected bit sequence
% flag_2 is a flag variable indicating whether the decoding process was
% stuck or not
% countFig2 number of figures
% count figure ID
accumulate=zeros(1,n-k); % Initializing a vector (a cell for each check node) to check that the output (out) satisfy the parity check equations in the H matrix
flag_breaked=0; % A flag to indicate whether the decoder is able to complete its procedure normally or not
% The next for loop calculates the accumulate vector using the H and the x
% excluding the erased bits in x
for i=1:n-k
    a=x(H(i,:)==1);
    a(a==-1)=[];
    accumulate(i)=mod(sum(a),2);
end
H_updated=H;
H_updated(:,x~=-1)=0; % Equivalent to removing the edges (in the Bipartite graph) for the non erased bits in x
hFrame=figure();
set(hFrame, "visible", "off");
hold all
title('Graph including only edges from erased bits','fontsize',14);
draw_part_graph(x, H_updated, 1:n, countFig, filename, hFrame); % Draw the Bipartite graph for H_updated
% legend('Variable node','Check node');
countFig=countFig+1; % A figure will be added to the cell array so increment the counter
% The next while loop performs the correction of the erasures in x
count=0;
while sum(sum(H_updated))~=0
    degree_one=find(sum(transpose(H_updated))==1); % Finding the check nodes connected to only one erased bit in x. The decoder fails if there is no such nodes
    if numel(degree_one)==0
        flag_breaked=1;
        break;
    end
    for i=1:n-k
        % For every check node of degree one, calculate the value of the
        % erased bit in x connected to it, update the H then draw the Bipartite graph. Then recalculate the accumulate vector values.
        if sum(H_updated(i,:))==1
            count=count+1;
            x(H_updated(i,:)==1)=accumulate(i);
            H_updated(i,H_updated(i,:)==1)=0;
            for j=1:n-k
                b=x(H(j,:)==1);
                b(b==-1)=[];
                accumulate(j)=mod(sum(b),2);
            end
            H_updated(:,x~=-1)=0; % Equivalent to removing the edges (in the Bipartite graph) for the non erased bits in x
            hFrame=figure();
            set(hFrame, "visible", "off");
            hold all
            title(sprintf('After Iteration number %d',count),'fontsize',14);
            draw_part_graph(x, H_updated, 1:n, countFig, filename, hFrame); % Draw the Bipartite graph for H_updated
            %             legend('Variable node','Check node');
            countFig=countFig+1; % A figure will be added to the cell array so increment the counter
        end
        
    end
end
if flag_breaked==1
    % decoder fails
    out=[];
    flag_2=0;
    countFig2=countFig;
else
    % Check whether the decoder fails or not
    checking=mod(H*x',2);
    isCodeWord=sum(checking); % To check whether the output is a codeword or not
    isHZeroos=sum(sum(H_updated))~=0; % To check that all of the edges of the graph have been peeled.
    if sum(sum(H_updated))~=0 || isCodeWord~=0
        % Decoder fails
        out=[];
        flag_2=0;
        countFig2=countFig;
    else
        % Decoder finished it's task successfully :)
        out=x;
        flag_2=1;
        countFig2=countFig;
    end
end
end

function draw_part_graph (x, H, varnodes, countFig, filename, hFrame)
% This function draws the Bipartite graph of a code with a certain H matrix
% The inputs to the function are:
% H The parity check matrix of the code under consideration
% varnodes The variable nodes indeces from 1 to n
% x the received bits
m=size(H,1);
allCheckNodes=1:m;
varnodepos = 1:1:numel(varnodes);
for k=1:numel(varnodepos)
	if x(k)==-1
	   h1=plot(varnodepos(k), 0,'r'); % Plotting the circles representing the variable nodes in red color since it is still erased
      set(h1,'marker','.','markersize',40);
	else
		h2=plot(varnodepos(k), 0,'g'); % Plotting the circles representing the variable nodes in green color since it is not erased
      set(h2,'marker','.','markersize',40);
	end
end
xlim([0 numel(varnodes)+1]);
hold all;
for k=1:numel(varnodepos)
    if x(k)==-1
        text(varnodepos(k)-0.22, 0.08, '?','fontsize',15); % Writing the bit received by the variable nodes
    else
        text(varnodepos(k)-0.22, 0.08, num2str(x(k)),'fontsize',15); % Writing the bit received by the variable nodes
    end
end
% find connected check nodes
vnidx=1;
multcn = [];
for vn = varnodes
    cns{vnidx} = find(H(:,vn));
    multcn = [multcn; cns{vnidx}];
    vnidx=vnidx+1;
end
allcn = unique(multcn);
cnpos = [1:1:numel(allCheckNodes)] + 0.5*(numel(varnodes)-numel(allCheckNodes));
ss=plot(cnpos, 2*ones(size(cnpos)), 'c'); % Plotting the circles representing the check nodes in cyan (all connected edges peeled)
set(ss,'marker','.','markersize',40);
for k=1:numel(allcn)
    sh=plot(cnpos(allcn(k)), 2,'m'); % Plotting the circles representing the check nodes in magneta (not all connected edges peeled)
    set(sh,'marker','.','markersize',40);
end
vnidx=1;
for vn = varnodes
    for k=1:numel(cns{vnidx})
        cnk = allcn == cns{vnidx}(k);
        line([varnodepos(vnidx) cnpos(allcn(cnk))],[0 2],'LineWidth',1,'Color','r'); % Plotting the edges all in red
    end
    vnidx = vnidx+1;
end
% highlight all edges which are connected to the check node exactly an even
% number of times
values = 1:size(H,1);
N = histc(multcn,values);
multcn = values(mod(N,2)==0 & N > 0);
vnidx=1;
for vn = varnodes
    connections = intersect(cns{vnidx},multcn);
    if ~isempty(connections)
        for k=1:numel(connections)
            cnk = allcn == connections(k);
            line([varnodepos(vnidx) cnpos(allcn(cnk))],[0 2],'LineWidth',2);
        end
    end
    vnidx=vnidx+1;
end
multcn = values(N>1 & mod(N,2)==1);
vnidx=1;
for vn = varnodes
    connections = intersect(cns{vnidx},multcn);
    if ~isempty(connections)
        for k=1:numel(connections)
            cnk = allcn == connections(k);
            line([varnodepos(vnidx) cnpos(allcn(cnk))],[0 2],'LineWidth',2,'Color','b'); % The edges connected to a check node of even degree are drawn in blue
        end
    end
    vnidx=vnidx+1;
end
hold off;
set(gca,'XTick',[]);
set(gca,'YTick',[]);
set(gca,'Box','off');
hold off;
new_filename = [filename '_' num2str(countFig)];
saveas(hFrame,new_filename,'png');
end

function [ H ] = generateH( v,c )
% This function generates the parity check matrix H of a regular LDPC code
% using the Gallager's method described in slide 3 with the following inputs:
% v: the variable node degree
% c: the check node degree
H1=zeros(2*v,2*c);
H_temp=zeros(2,2*c);
H_temp(1,1:c)=1;
H_temp(2,c+1:end)=1;
H1(1:2,1:2*c)=H_temp;
for i=2:v 
a=randperm(2*c);
H2=H_temp(:,a);
H1(2*i-1:2*i,:)=H2;
end
H=H1;
end