lccs.m

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

function [lccscount, lccs_path, lccs_str, lccstable] = lccs(a, b)
%LCCS Longest (maximum) common consecutive subsequence
%	Usage:
%	[count, lccs_path, lccs_str, lccstable] = lccs(a, b)
%	a: input string 1
%	b: input string 2
%	count: count of LCCS
%	lccs_path: optimal path of dynamical programming through the lccs table
%	lccs_str: LCCS string
%	lccstable: LCCS table for applying dynamic programming
%
%	Type "lccs" for a self-demo.

% Cosh Hsu, 20000428	

if nargin == 0, selfdemo; return; end

a = a(:).';
b = b(:).';
m = length(a);
n = length(b);
lccstable = zeros(m+1, n+1);
prevx = zeros(m+1, n+1);
prevy = zeros(m+1, n+1);
% Find LCCS using dynamic programming
for i=1:m,
	for j = 1:n,
		if a(i)==b(j),
			lccstable(i+1,j+1) = lccstable(i,j)+1;
		else
			lccstable(i+1,j+1) = lccstable(i+1,j);
			lccstable(i+1,j+1) = 0;
		end 
	end
end

% Get rid of initial conditions
lccstable = lccstable(2:end, 2:end);

% ====== Return length of LCCS string
lccscount =-1;
for i=(-(m-1)):n-1,
  	temp = diag(lccstable,i);
   tempCount = max(temp);
   if tempCount > lccscount,
     	lccsIdx = i;
     	lccscount = tempCount;         
  	end         
end

% ===== Return the optimal path of the dynamical programming
if nargout > 1,
   if lccscount>0,
   	temp = diag(lccstable,lccsIdx);
   	[tempCount,idx]=max(temp);
   	if lccsIdx<0,
         lccs_path=[((idx-lccscount+1):idx)-lccsIdx ; (idx-lccscount+1):idx]';
      elseif lccsIdx>=0,
         lccs_path=[((idx-lccscount+1):idx) ; ((idx-lccscount+1):idx)+lccsIdx]';         
      end
  	else
   	lccs_path=[];
   end
end


% ====== Return the LCCS string
if nargout > 2,	% return LCCS string
	temp = lccstable((lccs_path(:,2)-1)*m+lccs_path(:,1));  % LCCS count along the path
	temp = [0; temp];
	index = find(diff(temp));
	lccs_str = a(lccs_path(index,1));
end

% ====== Self demo
function selfdemo
str1 = 'abbscbdrasd';
str2 = 'kjabbdrassax';
m = length(str1);
n = length(str2);
figure;
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 LCCS
[count, lccs_path, lccs_str, lccstable] = feval(mfilename, str1, str2);
xlabel(['String1 = ', str1]);
ylabel(['String2 = ', str2]);
title(['LCCS table and LCCS path; with LCCS = ', lccs_str]);

% ====== Plot LCCS table
for i = 1:m,
	for j = 1:n,
		text(i, j, int2str(lccstable(i,j)), 'hori', 'center');
	end
end

% ====== Plot LCCS path
for i = 1:size(lccs_path,1)-1,
	line(lccs_path(i:i+1, 1), lccs_path(i:i+1, 2));
end

% ====== Circle matched elements
temp = lccstable((lccs_path(:,2)-1)*m+lccs_path(:,1)); % LCCS count along the path
temp = [0; temp];
index = find(diff(temp));
match_point = lccs_path(index, :);
line(match_point(:,1), match_point(:, 2), ...
	'marker', 'o', 'markersize', 15, 'color', 'r', 'linestyle', 'none');