function params = getLDPCEncoderParameters(parityMat)
%
%   Warning: This undocumented function may be removed in a future release.
%
% getLDPCEncoderParameters Get LDPC encoder parameters
%   PARAMS = getLDPCEncoderParameters(PARITYMAT) returns LDPC encoder
%   parameters for the parity check matrix, PARITYMAT.  The parameters are
%   used by the LDPC encoder kernel.

% Copyright 2010-2016 The MathWorks, Inc.

validateLDPCParityCheckMatrix(parityMat);

N = size(parityMat,2);
K = N - size(parityMat,1);

PB = parityMat(:,K+1:end); % Extract the last (N-K) columns of the parity-check matrix.

% Check if PB is triangular
switch isfulldiagtriangular(PB)
  case 1
    Alg = 'Forward Substitution';   % PB is lower triangular and has a full diagonal.
    rowOrder = -1;                  % Don't need to reverse the order of rows in PB
  case -1
    Alg = 'Backward Substitution';  % PB is upper triangular and has a full diagonal.
    rowOrder = -1;                  % Don't need to reverse the order of rows in PB
  case 0
    % Reverse the order of rows in PB, but keep PB, since if PB is not
    % triangular, we need to factorize it in GF(2).
    Reversed_PB = PB(end:-1:1,:);
    switch isfulldiagtriangular(Reversed_PB) % Check if Reversed_PB is triangular.
      case 1
        Alg = 'Forward Substitution';
        rowOrder = (N-K):-1:1;
        PB = Reversed_PB;
      case -1
        Alg = 'Backward Substitution';
        rowOrder = (N-K):-1:1;
        PB = Reversed_PB;
      case 0
        Alg = 'Matrix Inverse';
    end
end

% The following properties won't be used for forward/backward substitution.
% But the C-MEX function expects these inputs. So just provide dummy
% values by default.
params.MatrixL_RowIndices  = int32(0);
params.MatrixL_ColumnSum   = int32(0);
params.MatrixL_RowStartLoc = int32(0);

if strcmp(Alg, 'Forward Substitution')
  params.EncodingMethod = int8(1);
  P = tril(PB, -1);  % Remove the diagonal of PB and put it in P.
elseif strcmp(Alg, 'Backward Substitution')
  params.EncodingMethod = int8(-1);
  P = triu(PB, 1);   % Remove the diagonal of PB and put it in P.
else
  % If we reach this line, then Alg = 'Matrix Inverse',
  % i.e., the last (N-K) columns of the parity-check matrix (i.e. PB) is not triangular.
  % Let's try to factorize it in GF(2).
  % See the comments in gf2factorize for details about what it does.
  [PL, PB, rowOrder, invertible] = gf2factorize(PB);
  if ~invertible
    error(message('comm:getLDPCEncoderParameters:NonInvertibleParityCheckMatrix'));
  else
    % PB is invertible in GF(2), and has been modified by gf2factorize.
    params.EncodingMethod = int8(0);
    [params.MatrixL_RowIndices, ...
     params.MatrixL_RowStartLoc, ...
     params.MatrixL_ColumnSum] = ConvertMatrixFormat(tril(PL, -1));
    PB = PB(rowOrder, :); % Need to do this before the next line.
    P = triu(PB, 1);      % Now PB is upper triangular. Remove the diagonal and put it in P.
  end
end

% Check for P to be non-zero (=> need more rows in parityMat)
if ( isequal(size(P), [1 1]) && isequal(full(P), zeros(1,1)) )
    error(message('comm:getLDPCEncoderParameters:InvParityCheckMatrix'));
end

% Update all internal data structures for the encoding operation.
params.RowOrder = int32(rowOrder-1);
[params.MatrixA_RowIndices, ...
 params.MatrixA_RowStartLoc, ...
 params.MatrixA_ColumnSum] = ConvertMatrixFormat(parityMat(:,1:K));
[params.MatrixB_RowIndices, ...
 params.MatrixB_RowStartLoc, ...
 params.MatrixB_ColumnSum] = ConvertMatrixFormat(P);

% Update all external properties.
params.NumInfoBits = K;
params.NumParityBits = N - K;
params.BlockLength = N;
params.EncodingAlgorithm = Alg;
end

%------------------------------------------------------------------
% Local functions
function shape = isfulldiagtriangular(X)
% X must be a square logical matrix.
% shape = 1  if X is lower triangular and has a full diagonal.
% shape = -1 if X is upper triangular and has a full diagonal.
% Otherwise, shape = 0.
% If X is a diagonal matrix, X is considered upper triangular.

N = size(X,1);
NumNonzeros = nnz(X);
if ~all(diag(X)) % Full diagonal?
  shape = 0;   % Must have a full diagonal.
else
  NumNonzerosInLowerPart = nnz(tril(X));
  if NumNonzerosInLowerPart == NumNonzeros
    shape = 1;  % X is lower triangular.
  elseif NumNonzerosInLowerPart == N
    shape = -1; % X is upper triangular.
  else
    shape = 0;  % X is not triangular.
  end
end
end

function [RowIndices, RowStartLoc, ColumnSum] = ConvertMatrixFormat(X)
% Create an alternative representation of a zero-one matrix.

% Find nonzero elements.
[i,~] = find(X);

% Generate zero-based row indices and use int32 for C-MEX function.
if ~isempty(i)
    RowIndices = int32(i-1);
else
    % In this case, RowIndices' value doesn't matter. But we still return
    % a nonempty matrix to avoid error from Simulink block mask.
    RowIndices = int32(0);
end

% Find the number of nonzero elements in each column.
% Use int32 for C-MEX function.
ColumnSum  = int32(full(sum(X)));

% For each column, find where the corresponding row indices start in RowIndices.
% Generate zero-based indices and use int32 for C-MEX function.
CumulativeSum = cumsum(double(ColumnSum));
RowStartLoc = int32([0 CumulativeSum(1:end-1)]);
end

function [A, B, chosen_pivots, invertible] = gf2factorize(X)
%GF2FACTORIZE  Factorize a square matrix in GF(2).
%   [A B chosen_pivots invertible] = gf2factorize(X) factorizes a square matrix X
%   in GF(2) using Gaussian elimination.
%
%   X = A * B    (in GF(2), i.e., using modulo-2 arithmetic)
%
%   X may be a sparse matrix. Nonzero elements are treated as ones.
%   A and B are sparse logical matrices.
%
%   A is always lower triangular.
%   If X is invertible in GF(2), then B(chosen_pivots,:) is upper triangular and
%   invertible = true.
%   If X is non-invertible in GF(2), then X = A * B still holds and invertible = false.
%
%   To evaluate A * B in GF(2), use mod(A*double(B),2).

n = size(X,1);
if n ~= size(X,2)
  error(message('comm:getLDPCEncoderParameters:NonSquareInputForFactorization'));
end

Y1 = logical(eye(n,n));
Y2 = full(X~=0);
chosen_pivots = zeros(n,1); % Haven't chosen any pivots yet.
invertible = true; % Assume that X is invertible in GF(2).

for col = 1:n
  candidate_rows = Y2(:,col);
  candidate_rows(chosen_pivots(1:(col-1))) = 0; % Never use a chosen pivot.
  candidate_rows = find(candidate_rows); % Convert into row indices.
  
  if isempty(candidate_rows)
    invertible = false; % X is not invertible in GF(2).
    break;
  else
    pivot = candidate_rows(1);      % Choose the first candidate row as the pivot row.
    chosen_pivots(col) = pivot;     % Record this pivot.
    
    % Find all nonzero elements in the pivot row.
    % They will be xor-ed with the corresponding elements in other
    % candidate rows.
    columnind = find(Y2(pivot,:));
    
    % Subtract the pivot row from all other candidate_rows.
    % Exploit the fact that we are working in GF(2).
    % Just use logical NOT.
    % As we continue in this loop, Y2 will become "psychologically"
    % upper triangular.
    Y2(candidate_rows(2:end), columnind) = not(Y2(candidate_rows(2:end), columnind));
    
    % Update the lower triangular matrix Y1.
    Y1(candidate_rows(2:end), pivot) = 1;
  end
end

% Output sparse matrices to save memory.
A = sparse(Y1);
B = sparse(Y2);

if ~invertible
  % Output the chosen pivots even if X is not invertible in GF(2).
  chosen_pivots = [chosen_pivots(1:(col-1)); setdiff((1:n)', chosen_pivots(1:(col-1)))];
end
end
