lcs.m

一个关于数据聚类和模式识别的程序,在生物化学,化学中因该都可以用到.希望对大家有用,谢谢支持

function [lcscount, lcsPath, lcsStr, lcsTable] = lcs(a, b, plotOpt)
%lcs: Longest (maximum) common subsequence
%	Usage: [count, lcsPath, lcsStr, lcsTable] = lcsm(a, b, plotOpt)
%		a: input string 1
%		b: input string 2
%		plotOpt: plot option
%		count: count of LCS
%		lcsPath: optimal path of dynamical programming through the lcs table
%		lcsStr: LCS string
%		lcsTable: LCS table for applying dynamic programming
%
%	Type "lcs" for a self-demo.

%	Roger Jang, 981226, 990409, 20060528

if nargin<1, selfdemo; return; end
if nargin<3, plotOpt=0; end

a = a(:).'; m = length(a);
b = b(:).'; 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 | plotOpt
	now = [m, n];
	prev = [prevx(now(1), now(2)), prevy(now(1), now(2))];
	lcsPath = now;
	while all(prev>0),
		now = prev;
		prev = [prevx(now(1), now(2)), prevy(now(1), now(2))];
		lcsPath = [lcsPath; now];
	end 
	lcsPath = flipud(lcsPath);
end

% ====== Return the LCS string
if nargout>2 | plotOpt		% return LCS string
	temp = lcsTable((lcsPath(:,2)-1)*m+lcsPath(:,1));	% LCS count along the path
	temp = [0; temp];
	index = find(diff(temp));
	lcsStr = a(lcsPath(index,1));
end

% ====== Plot
if plotOpt
	dpPathPlot4strMatch(a, b, lcsPath, lcsTable, prevx, prevy);
end

% ====== Self demo
function selfdemo
str1='intention';
str2='execution';
plotOpt=1;
[lcsCount, lcsPath, lcsStr, lcsTable] = feval(mfilename, str1, str2, plotOpt);