editdistance.m

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

function [minDist, edPath, edTable] = editDistance(a, b, substituteCost, plotOpt)
%editDistance: Edit distance (ED) via dynamic programming
%	Usage: [minDist, edPath, edTable] = editDistance(a, b, substituteCost, plotOpt)
%		a: input string 1
%		b: input string 2
%		plotOpt: plot option
%		minDist: minimum edit distance
%		edPath: optimal path of dynamical programming through the ED table
%		edTable: ED table for applying dynamic programming
%
%	Type "editDistance" for a self-demo.

%	Roger Jang, 20081008

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

% ====== Return length of ED 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

% ====== Plot
if plotOpt
	dpPathPlot4strMatch(a, b, edPath, edTable, prevx, prevy);
	title(sprintf('Min. edit distance = %d', edTable(end, end)));
end

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