lpdistm.m

模式识别工具箱,希望对大家有用!

%LPDISTM  An L^p (p > 0) (non)-metric distance matrix between two data sets
% 
%     D = lpdistm (A,B,p)
%       or
%     D = lpdistm (A,B)
%       or
%     D = lpdistm (A,p)
%       or
%     D = lpdistm (A)
% 
%  Computation of the distance matrix D between two sets of vectors: 
%  A and B. If A is n-by-p and B is m-by-p, then D is n-by-m matrix.
%  Distances are computed using the L^p (p > 0) (non)-metric:
%
%     d(x,y) = (sum (|x_i - y_i|.^p))^(1/p)   
%                 i
%  Important:
%  p >= 1      => d is a metric
%  p in (0,1)  => d is not a metric              
%
%  Remark:
%  p = 1,2   => d - city block distance, Euclidean distance     
%
%  DEFAULT:
%    p = 1
%

% Elzbieta Pekalska, [e-mail]
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands



function d = lpdistm (a,b,p)

isda = isa(a,'dataset');
laba = getlab(a);
a    = +a;

[ra,ca] = size(a);
bisa = 0; 
if nargin < 2,
  p = 1;
  bisa = 1; 
else
  [rb,cb] = size(b);

  if nargin < 3,
	 if max (rb,cb) == 1, 
   	p = b;
   	bisa = 1; 
	 else
   	p = 1; 
	 end;    
  end
end  


if p < 0,  
  error ('The parameter p must be positive.');
end;

if ~bisa & ca ~= cb,  
  error ('The matrices should have the same number of columns.');
end;


if bisa,
  d = zeros(ra,ra);
  [i j] = find (triu(ones(ra),1));
  d(i + ra*(j-1)) = sum ((abs(a(i,:) - a(j,:))).^p, 2);
  d(j + ra*(i-1)) = d(i + ra*(j-1));

  if isda,
    d = dataset(d, laba, laba);
  end

else

  isdb = isa(b,'dataset');
  if ~xor(isda, isdb),
	 labb = getlab(b);
	 b    = +b;
  else
	 error('One matrix is a dataset and the other not.')
  end

  d = zeros(ra,rb);
  d = sum ((abs (repmat (permute(a,[1 3 2]), [1 rb 1]) - ... 
                 repmat (permute(b,[3 1 2]), [ra 1 1]))).^p,3); 

%  a = permute(a, [ 1 3 2 ]);
%  b = permute(b, [ 3 1 2 ]);
%  d = sum ((a(:,ones(1,rb),:) - b(ones(1,ra),:,:)).^p,3);

  if isda & isdb,
    d = dataset(d, laba, labb);
  end
end
d = d.^(1/p);

d(find (d < eps)) = 0;
return