function[MM,Zb_str]=bptmaxmatch(Zb)
% This code takes a structured matrix Zb, and produce a matrix Zb_str such
% that none of its entries present in a same row or column. Furthermore, a
% matrix H is computed, whose non-zero rows corresponds to the edges
% associated with a maximum matching of the bipartite graph. For this,
% first, the matrix Zb is represented with a bipartite graph, and then,
% represented it as a flow network graph, by adding extra two nodes:
% soiurce and sink. Then, using the "maxflow" command in MATLAB,
% a maximum matching of the underlying bipartite graph is
% obtained. Using these edges the matrix Zb_str is constructued.
% This code is written on 21.09.2022.
%% ========================================================================
%% Bipartite graph representation of matrix Zb
sz_Zb = size(Zb);
szr_Zb = sz_Zb(1);
szc_Zb = sz_Zb(2);
st_edg = []; % Starting edges of the bipartite graph
tr_edg =[]; % Terminating edges of the bipartite graph
for i = 1:1:szr_Zb
for j = 1:1:szc_Zb
if Zb(i,j) ~= 0
st_edg = [st_edg j+1]; % column vertices are started with
% vertex number 2 and so on, since the
% source node will be indexed as 1.
tr_edg = [tr_edg szc_Zb+1+i]; % row vertices are started with
% the number of columns of Zb plus
% 2.
end
end
end
%% Network Flow representation of bipartite graph
cw = 2:1:szc_Zb+1;
rw = szc_Zb+2:1:szc_Zb+szr_Zb+1 ;
st_edg_ext = [ones(1,szc_Zb) st_edg rw ];% All the starting nodes (c_i) of
% bipartite graph are connected
% from source node 1.
tr_edg_ext = [cw tr_edg (szc_Zb+szr_Zb+2)*ones(1,szr_Zb)]; % All the
% terminating nodes of the
% bipartite graph connected to the
% sink node.
Zb_bpt_edg = [st_edg_ext;tr_edg_ext]; %
%% Maxflow and maximum matching edges extraction
siz_Zb_bpt = size(Zb_bpt_edg);
szc_Zb_bpt = siz_Zb_bpt(2);
wt_mf = ones(1,szc_Zb_bpt);
Dg = digraph(st_edg_ext,tr_edg_ext,wt_mf);
% plot(Dg,'EdgeLabel',Dg.Edges.Weight,'Layout','layered')
[~,GF] = maxflow(Dg,1,szc_Zb+szr_Zb+2,'augmentpath');
H = GF.Edges.EndNodes; % This will produce the edges associated with maxflow
siz_H = size(H);
Lh = siz_H(1,1);
for i = 1:1:Lh
if ((H(i,1) == 1) || (H(i,2) == szc_Zb+szr_Zb+2))
H(i,:) = [0 0];
end
end
MM = H; % The non-zero rows of MM are the maximum
% matching edges of the bipartite graph
%% Zsb matrix computation
Zb_str = zeros(szr_Zb,szc_Zb); % This matrix is next constructed using the
% nonzero row informations of H matrix
for k = 1:1:Lh
if H(k,1) ~= 0
nz_col = H(k,1)-1; % The vertex indices are mapped to the original
% row column vertex indices
nz_row = H(k,2)-(szc_Zb+1);
Zb_str(nz_row,nz_col) = Zb(nz_row,nz_col);
end
end