% This program constructs coumn-weight two girth 12 quasi-cyclic LDPC codes given
% k( row weight, j is fixed at 2) and m (size of each sub-matrix). The
% program constructs a distance graph first then converts it into a matrix.
% For details of the algorithm see G.Malema and M. Liebelt,"Quasi-cyclic
% LDPC Codes of Column-Weight two using a Search Algorithm", EURASIP journal
% on Advances in Signal Processing
%or PhD thesis by G.Malema, LDPC codes: Construction and Implementation
%(University of Adelaide 2007).
% Author : Dr. Gabofetswe Alafang Malema,University of Botswana  
% Department of Computer science. e-mail: malemag@mopipi.ub.bw


j=2; % column-weight should be fixed at 2.
k=8; % row-weight is variable/changed.
m=70; % size of each sub-matrix can be changed.
M=j*m; % total number of rows.
shift=zeros(k,j); % stores shift values for each sub-matrix.
rows = struct('connect',0,'counter',0,'con',0); % structure used for storing connections.
emptyset = [];
for y = 1:M       % connections and number of connections for rows and columns are initialized.
    rows(y).counter = 0;
    rows(y).con = intersect(emptyset,rows(y).con);
    rows(y).connect = intersect(emptyset,rows(y).connect);
    %rows(y).mem = intersect(emptyset,rows(y).mem);
end

group1 = [1:m];  %row group one
group2 = [m+1:2*m];  % row group two,

g= 3; % determines the girth of the code, in these case 12.
found_code = 1; % true if algorithm does not fail.

% the for loop finds a rows(from groups 1 and 2) that are apart at least a distance
% of 6 (12 in matrix form). They are then connected. The rest of the groups are connected
% according to this first connection. The process is repeated k times which
% is the number of connections for each row.

for kk = 1:k
    r1 = ceil(rand*length(group1)); % randomly choose a row from group1
    i = group1(r1);
    
       % create a list of rows within a distance of g from row i. 
       mem=[]; mem1 = rows(i).con;
        mem2 = [];
        for y = 1:g
            for x = 1:length(mem1)
                x1 = mem1(x);
                mem2 = union(mem2,rows(x1).con);
            end
            mem = union(mem1,mem2);
            mem1=mem;
        end
        
    % find row not in mem in the given range.
    A = intersect(mem,group2);
    A = setxor(A,group2);     % set of rows at least a distance of g from row i.

    % the while loop below looks for rows that satisfy the distance g and
    % avoids eight-cycles.
    done=0;
    while done == 0
        if isempty(A)~=1
            r1= ceil(rand*length(A));
            row1 = A(r1);
            
            dist1 = []; dist2 = [];
            for gg = 1:rows(row1).counter
                r1 = rows(row1).connect(gg,1);
                dist1 = union(dist1,r1);
                for ff=1:rows(row1).counter
                    dist2 = union(dist2,rows(r1).connect(ff,1));
                    dist2 = union(dist2,rows(i).connect(ff,1));
                end
            end
            
            dist11=[];dist22=[];
            for xx = 1:length(dist1)
                if i < dist1(xx)
                    dist11 = union(dist11,dist1(xx)-i);
                else
                    dist11 = union(dist11,mod(dist1(xx)-i,m));
                end
            end
                
           for xx = 1:length(dist2)
               if row1 < dist2(xx)
                   dist22 = union(dist22, dist2(xx)-row1);
               else
                   dist22 = union(dist22, mod(dist2(xx)-row1,m));
               end
           end
           
           sel = intersect(dist11,dist22);        
    
           if isempty(sel)==1 
                row1= row1;
                done = 1;                 
                   
            else
                A = setxor(A,row1);
           end  
           
        else
           disp('Row one not found');
           found_code = 0;  %algorithm fails
           break;
        end % if isempty(A)~=1
    end % while done ==0
                                 
                         
    % connect all the rows by using a shift.
    u = find(group1 == i);
    y = find(group2 == row1);
        
    p0 =circshift(group1,[1 -u+1]);
    p1 =circshift(group2,[1 -y+1]);
    
    
     shift(kk,1)=u; % stores shift values of the sub-matrices.
     shift(kk,2)=y;
         
    a = kk;
    
    for b=1:m
        
        g0=p0(b);g1=p1(b);
        rows(g0).connect(a,1)=g1;               
        rows(g1).connect(a,1) = g0;
                        
        rows(g0).con = union(rows(g0).con, g1);
        rows(g1).con = union(rows(g1).con, g0);
        
        rows(g0).counter = rows(g0).counter + 1;
        rows(g1).counter = rows(g1).counter + 1;
            
    end
        
    %end % while loop. 
end % for kk =1:k

 count = zeros(1,M);
 for z = 1:M
     count(z)= rows(z).counter;
 end             
 
 %create a sparse matrix H from the row connections above.
 counter = 1; 
 H = zeros(M,5);

 if found_code == 1
 for i = 1:k
     for z = 1:m
         x1 = rows(z).connect(i,1);
         H(z,counter)=1;
         H(x1,counter)=1;
         counter = counter + 1;
     end
 end
 else
     disp('Algorithm failed to get code. Try again or increase group-size m.');
 end