使用matlab对二部图进行匹配（匈牙利算法）

使用matlab应用匈牙利算法对二部图进行匹配。
算法的代码主体来自知乎用户洪九(李戈)的专栏文章策略算法工程师之路-图优化算法(一)(二分图&最小费用最大流)

具体算法思想文章里和其它地方都有，不在这里赘述。主要通过例子来说明算法具体运行过程中所作的操作对匈牙利算法进行说明。

首先对于如下图的一个二部图

左侧有6个点，右侧有7个点，点之间的边都画出来了。

第一步：对$x_1$进行匹配，找到它的第一个可匹配项$y_1$，发现此时$y_1$还未被匹配，因此将$x_1$和$y_1$进行配对，完成对$x_1$的处理，进入下一步。

第二步：对$x_2$进行匹配，找到它的第一个可匹配项$y_2$，发现此时$y_2$还未被匹配，因此将$x_2$和$y_2$进行配对，完成对$x_2$的处理，进入下一步。

第三步：对$x_3$进行匹配，找到它的第一个可匹配项$y_1$，但是现在$y_1$已经被$x_1$给匹配了。
因此询问$x_1$你能不能腾腾地儿，把$y_1$让出来。
$x_1$一看自己还可以选$y_2$，但是$y_2$又已经被$x_2$给匹配了。
因此询问$x_2$能不能腾腾地儿，把$y_2$让出来。
$x_2$发现自己还能匹配$y_5$，并且$y_5$还未被匹配，不冲突。
于是$x_2$将自己的匹配项从$y_2$改到$y_5$，并告诉$x_1$我已经把$y_2$腾出来了。
$x_1$接到$y_2$已经被腾出来的信号后，将自己的匹配项从$y_1$改到$y_2$，并告诉$x_3$我已经把$y_1$腾出来了。
$x_3$接到$y_1$已经被腾出来的信号后，直接将$x_3$和$y_1$进行配对。完成对$x_3$的处理，进入下一步。

第四步：对$x_4$进行匹配，找到它的第一个可匹配项$y_3$，发现此时$y_3$还未被匹配，因此将$x_4$和$y_3$进行配对，完成对$x_4$的处理，进入下一步。

第五步：对$x_5$进行匹配，找到它的第一个可匹配项$y_4$，发现此时$y_4$还未被匹配，因此将$x_5$和$y_4$进行配对，完成对$x_5$的处理，进入下一步。

第六步：对$x_6$进行匹配，找到它的第一个可匹配项$y_4$，但是现在$y_4$已经被$x_5$给匹配了。
因此询问$x_5$你能不能腾腾地儿，把$y_4$让出来。
$x_5$一看自己只有$y_4$一个可选项，没法让，因此告诉$x_6$我让不出来。
$x_6$接到$y_4$已经让不出来的信号后，再检视自身是否有其它可选项，结果是没有，因此$x_6$不进行配对，完成对$x_6$的处理。

至此所有$x$点都被处理完毕，算法结束。最大匹配包含 5条匹配边。

具体代码如下，可以将M和N的值改变，并通过随机数生成M与N之间的连接边进行测试：

clc
clear
close all

global M
global N
global Map
global p
global vis
% M,N分别表示左右侧集合的元素数量
M = 6;
N = 7;
% 邻接矩阵存图
Map = [1 1 0 1 0 0 0;
    0 1 0 0 1 0 0;
    1 0 0 1 0 0 1;
    0 0 1 0 0 1 0;
    0 0 0 1 0 0 0;
    0 0 0 1 0 0 0];
% Map = rand(M, N) > 0.7;

% 记录当前右侧元素所对应的左侧元素
p = zeros(N, 1);
% 记录右侧元素是否已被访问过
vis = zeros(N, 1);

cnt = Hungarian();


function cnt = Hungarian()
global M
global N
global Map
global p
global vis

graph_num = 1;
cnt = 0;
display_graph(graph_num, '原始二部图');
for i = 1: M
    % 重置vis数组
    vis = zeros(N, 1);
    if match(i)
        cnt = cnt + 1;
    end
    graph_num = graph_num + 1;
    display_graph(graph_num, ['第' num2str(i) '步']);
end
graph_num = graph_num + 1;
display_graph(graph_num, '最终结果');
end

function result = match(i)
global M
global N
global Map
global p
global vis
for j = 1: N
    % 有边且未访问
    if Map(i, j) && ~vis(j)
        % 记录状态为访问过
        vis(j) = 1;
        % 如果暂无匹配，或者原来匹配的左侧元素可以找到新的匹配
        if p(j) == 0 || match(p(j))
            % 当前左侧元素成为当前右侧元素的新匹配
            p(j) = i;
            % 返回匹配成功
            result = true;
            return;
        end
    end
end
% 循环结束，仍未找到匹配，返回匹配失败
result = false;
return;
end

function display_graph(graph_num, title_name)
global M
global N
global Map
global p
global vis

max_num = max(M, N);
temp = max_num - 1;
figure(graph_num);
cla
set( gca, 'XTick', [], 'YTick', [] );
set( gca, 'TickLength', [0 0]);
box on
hold on
xlim([-1, 2]);
ylim([-max_num, 1]);
% 绘制边
[row,col] = find(Map);
for j = 1: length(row)
   plot([0, 1], [(row(j) - 1) * -temp/(M-1), (col(j) - 1) * -temp/(N-1)], 'k', 'LineWidth', 1); 
end
% 绘制当前匹配边
for j = 1: length(p)
    if p(j) ~= 0
        plot([0, 1], [(p(j) - 1) * -temp/(M-1), (j - 1) * -temp/(N-1)], 'r', 'LineWidth', 2); 
    end
end
% 绘制点
scatter(zeros(1, M), 0:-temp/(M-1):-temp, 20, [217/255 83/255 25/255], 'filled');
scatter(ones(1, N), 0:-temp/(N-1):-temp, 20, [0 114/255 189/255], 'filled');
for j = 1: M
    text(-0.2, (j - 1) * -temp/(M-1), ['x_' num2str(j)]);
end
for j = 1: N
    text(1.1, (j - 1) * -temp/(N-1), ['y_' num2str(j)]);
end

title(title_name);
end