function [lcscount, lcs_path, lcs_str, lcstable] = lcs(a, b)
%LCS Longest (maximum) common subsequence
% Usage:
% [count, lcs_path, lcs_str, lcstable] = lcsm(a, b)
% a: input string 1
% b: input string 2
% count: count of LCS
% lcs_path: optimal path of dynamical programming through the lcs table
% lcs_str: LCS string
% lcstable: LCS table for applying dynamic programming
%
% Type "lcsm" for a self-demo.
% Roger Jang, 981226
% Roger Jang, 990409
if nargin == 0, selfdemo; return; end
a = a(:).';
b = b(:).';
m = length(a);
n = length(b);
lcstable = zeros(m+1, n+1);
prevx = zeros(m+1, n+1);
prevy = zeros(m+1, n+1);
% Find LCS using dynamic programming
for i=1:m,
for j = 1:n,
if a(i)==b(j),
lcstable(i+1,j+1) = lcstable(i,j)+1;
prevx(i+1,j+1) = i;
prevy(i+1,j+1) = j;
elseif lcstable(i,j+1) > lcstable(i+1,j),
lcstable(i+1,j+1) = lcstable(i,j+1);
prevx(i+1,j+1) = i;
prevy(i+1,j+1) = j+1;
else
lcstable(i+1,j+1) = lcstable(i+1,j);
prevx(i+1,j+1) = i+1;
prevy(i+1,j+1) = j;
end
end
end
% Get rid of initial conditions
lcstable = lcstable(2:end, 2:end);
prevx = prevx(2:end, 2:end)-1;
prevy = prevy(2:end, 2:end)-1;
% ====== Return length of LCS string
lcscount = lcstable(m, n);
% ====== Return the optimal path of the dynamical programming
if nargout > 1,
now = [m, n];
prev = [prevx(now(1), now(2)), prevy(now(1), now(2))];
lcs_path = now;
while all(prev>0),
now = prev;
prev = [prevx(now(1), now(2)), prevy(now(1), now(2))];
lcs_path = [lcs_path; now];
end
lcs_path = flipud(lcs_path);
end
% ====== Return the LCS string
if nargout > 2, % return LCS string
temp = lcstable((lcs_path(:,2)-1)*m+lcs_path(:,1)); % LCS count along the path
temp = [0; temp];
index = find(diff(temp));
lcs_str = a(lcs_path(index,1));
end
% ====== Self demo
function selfdemo
str1 = 'abxjhbscdrasngssfd';
str2 = 'sdkjtanhsfljskge';
m = length(str1);
n = length(str2);
figure;
%[xx, yy] = meshgrid(1:m, 1:n);
%plot(xx(:), yy(:), '.');
axis([0 m+1 0 n+1]);
box on;
set(gca, 'xtick', 1:m);
set(gca, 'ytick', 1:n);
set(gca, 'xticklabel', char(double(str1)'));
set(gca, 'yticklabel', char(double(str2)'));
% ====== invoke LCS
[count, lcs_path, lcs_str, lcstable] = feval(mfilename, str1, str2);
xlabel(['String1 = ', str1]);
ylabel(['String2 = ', str2]);
title(['LCS table and LCS path; with LCS = ', lcs_str]);
% ====== Plot LCS table
for i = 1:m,
for j = 1:n,
text(i, j, int2str(lcstable(i,j)), 'hori', 'center');
end
end
% ====== Plot LCS path
for i = 1:size(lcs_path,1)-1,
line(lcs_path(i:i+1, 1), lcs_path(i:i+1, 2));
end
% ====== Circle matched elements
temp = lcstable((lcs_path(:,2)-1)*m+lcs_path(:,1)); % LCS count along the path
temp = [0; temp];
index = find(diff(temp));
match_point = lcs_path(index, :);
line(match_point(:,1), match_point(:, 2), ...
'marker', 'o', 'markersize', 15, 'color', 'r', 'linestyle', 'none');
