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
%		lcsCount: 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, 20081009

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, n);
prevx=zeros(m, n);
prevy=zeros(m, n);
cost=zeros(3,1);
% Compute the first element
lcsTable(1,1)=a(1)==b(1);
prevx(1,1)=0;		% not plotted
prevy(1,1)=0;		% not plotted
% Compute the first row & column
for i=2:m
	if a(i)==b(1)
		lcsTable(i,1)=1;
		prevx(i,1)=i-1;
		prevy(i,1)=0;
	else
		lcsTable(i,1)=lcsTable(i-1,1);
		prevx(i,1)=i-1;
		prevy(i,1)=1;
	end
end
for j=2:n
	if a(1)==b(j)
		lcsTable(1,j)=1;
		prevx(1,j) = 0;
		prevy(1,j) = j-1;
	else
		lcsTable(1,j)=lcsTable(1,j-1);
		prevx(1,j) = 1;
		prevy(1,j) = j-1;
	end
	
end
% Compute all the other elements
for i=2:m,
	for j=2:n,
		if a(i)==b(j)
			lcsTable(i,j)=lcsTable(i-1,j-1)+1;
			prevx(i,j)=i-1;
			prevy(i,j)=j-1;
		else
			[lcsTable(i,j), index]=max([lcsTable(i-1, j), lcsTable(i, j-1)]);
			switch index
				case 1
					prevx(i,j)=i-1;
					prevy(i,j)=j;
				case 2
					prevx(i,j)=i;
					prevy(i,j)=j-1;
			end
		end
	end
end

% ====== Return length of ED 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);
	title(sprintf('LCS = %s', lcsStr));
end

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