editdistance - 複製.m

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

function [minDist, edPath, edStr, edTable] = editDistance(a, b, substituteCost, plotOpt)
%lcs: Longest (maximum) common subsequence
%	Usage: [count, edPath, edStr, edTable] = lcsm(a, b, plotOpt)
%		a: input string 1
%		b: input string 2
%		plotOpt: plot option
%		count: count of LCS
%		edPath: optimal path of dynamical programming through the lcs table
%		edStr: LCS string
%		edTable: 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, substituteCost=2; end
if nargin<4, plotOpt=0; end

a = a(:).'; m = length(a);
b = b(:).'; n = length(b);
edTable = zeros(m+1, n+1);
edTable(:,1)=0:m;
edTable(1,:)=(0:n)';
prevx = zeros(m+1, n+1);
prevy = zeros(m+1, n+1);
cost=zeros(3,1);
% Find LCS using dynamic programming
for i=1:m,
	for j = 1:n,
		if a(i)==b(j)
			cost=0;
		else
			cost=substituteCost;
		end
		cost(1)=edTable(i, j)+cost;
		cost(2)=edTable(i, j+1)+1;
		cost(3)=edTable(i+1, j)+1;
		[edTable(i+1,j+1), index]=min(cost);
		switch index
			case 1
				prevx(i+1,j+1) = i;
				prevy(i+1,j+1) = j;
			case 2
				prevx(i+1,j+1) = i;
				prevy(i+1,j+1) = j+1;
			case 3
				prevx(i+1,j+1) = i+1;
				prevy(i+1,j+1) = j;
		end
	end
end

% Get rid of initial conditions
edTable = edTable(2:end, 2:end);
prevx = prevx(2:end, 2:end)-1;
prevy = prevy(2:end, 2:end)-1;

% ====== Return length of LCS string
minDist = edTable(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))];
	edPath = now;
	while all(prev>0),
		now = prev;
		prev = [prevx(now(1), now(2)), prevy(now(1), now(2))];
		edPath = [edPath; now];
	end 
	edPath = flipud(edPath);
end

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

% ====== Plot
if plotOpt
	edPathPlot(a, b, edPath, edStr, edTable, prevx, prevy);
end


% ====== edPathPlot: Plot the path of ED
function edPathPlot(str1, str2, lcsPath, lcsStr, lcsTable, prevx, prevy)
m = length(str1);
n = length(str2);
%[xx, yy] = meshgrid(1:m, 1:n);
%plot(xx(:), yy(:), '.');
plot(nan, nan); axis([0 m+1 0 n+1]); box on;
set(gca, 'xtick', 1:m);
set(gca, 'ytick', 1:n);
set(gca, 'xticklabel', str1(:));
set(gca, 'yticklabel', str2(:));
xlabel(['String1 = ', str1]);
ylabel(['String2 = ', str2]);
title(['LCS table and LCS path; with LCS = ', lcsStr]);
% === Plot LCS table
for i = 1:m
	for j = 1:n
		text(i, j, int2str(lcsTable(i,j)), 'hori', 'center');
	end
end
% === Plot prevPos
for i=2:m+1
	for j=2:n+1
		now=i+j*sqrt(-1);
		next=prevx(i-1,j-1)+1+(prevy(i-1,j-1)+1)*sqrt(-1);
		arrowPlot(now, next, [1 0 0]);
	end
end
% === Plot LCS path
for i = 1:size(lcsPath,1)-1,
	line(lcsPath(i:i+1, 1), lcsPath(i:i+1, 2));
end
% === Circle matched elements
temp = lcsTable((lcsPath(:,2)-1)*m+lcsPath(:,1)); % LCS count along the path
temp = [0; temp];
index = find(diff(temp));
match_point = lcsPath(index, :);
line(match_point(:,1), match_point(:, 2), 'marker', 'o', 'markersize', 15, 'color', 'r', 'linestyle', 'none');


% ====== Self demo
function selfdemo
str1='execution';
str2='intention';
substituteCost=2;
plotOpt=1;
[minDist, edPath, edStr, edTable] = feval(mfilename, str1, str2, substituteCost, plotOpt);