function distout = VecDist(mvec1,mvec2,inorm) 
% VecDist, VECtor DISTances, L_p version
%     Finds the distances either:
%           Between all pairs of the column vectors of mvec1
%                             or:
%           Between the column vectors of mvec1, 
%           and the corresponding columns of mvec2,
%       using the L_p norm.
%   
%   Can use first 1, 2 or 3 arguments.
% Inputs:
%     mvec1  -  d x n (or 1) matrix of d-dim column vectors
%     mvec2  -  d x n (or 1) matrix of d-dim column vectors
%                     (when not specified, returns matrix of
%                      pairwise distances between cols of mvec1)
%     inorm  -  indicator of norm type:
%                           1 - L_1 norm
%                           2 - L_2 norm    (default when not specified)
%                         inf - sup norm
%                      (can be in 2nd argument, when mvec2 not given)
%
% Output:
%     distout     -  row vector (or matrix when only mvec1 is given)
%                        of distances between vectors
%
% Assumes path can find personal function:
%    vec2mat

%    Copyright (c) J. S. Marron 1999





if nargin == 1 ;    

  calctype = 1 ;
          %  flag for type of calculation: 
          %      all pairwise dist's of mvec1

  iinorm = 2 ;

end ;



if nargin == 2 ;

  if size(mvec2,1) == 1  &  size(mvec2,2) == 1 ;
                         %  then second entry is norm type

    calctype = 1 ;
          %  flag for type of calculation: 
          %      all pairwise dist's of mvec1

    iinorm = mvec2 ;

  else ;

    calctype = 2 ;
          %  flag for type of calculation: 
          %      dist's between mvec1 and mvec2

    iinorm = 2 ;

  end ;

end ;



if nargin == 3 ;

  calctype = 2 ;
          %  flag for type of calculation: 
          %      dist's between mvec1 and mvec2

  iinorm = inorm ;

end ;




if calctype == 1 ;    %  then do all pairwise dist's of mvec1

  d = size(mvec1,1) ;
  n = size(mvec1,2) ;


  if n == 1 ;   %  then no work to do


    distout = 0 ;


  else ;    %  then do calculations


    %  Construct 3d array of differences
    %
    a3diff = [] ;
    for j = 1:d ;
      mdiff = vec2mat(mvec1(j,:),n) - vec2mat(mvec1(j,:)',n) ;
          %  matrix of pairwise differences, for dim j
      a3diff = cat(3,a3diff,mdiff) ;
    end ;
          %  a3diff is n x n x d


    if iinorm == 1 ;    %  then do L_1 norm

      a3diff = abs(a3diff) ;
      distout = sum(a3diff,3) ;

    elseif iinorm == 2 ;    %  then do L_2 (Euclidean) norm

      a3diff = a3diff.^2 ;
      distout = sqrt(sum(a3diff,3)) ;

    elseif iinorm == inf ;    %  then do L_inf (sup) norm

      a3diff = abs(a3diff) ;
      distout = max(a3diff,[],3) ;
          %  take max along third dimension

    end ;


  end ;  



elseif calctype == 2 ;    %  then do dist's between mvec1 and mvec2

  d = size(mvec1,1) ;
  n1 = size(mvec1,2) ;
  n2 = size(mvec2,2) ;

  if size(mvec2,1) == d ;    %  then OK, go ahead

    if n1 == 1  & n2 > 1 ;    %  then expand mvec1 to n2
      mvec1 = vec2mat(mvec1,n2) ;
      n1 = n2 ;
    end ;

    if n2 == 1  & n1 > 1 ;    %  then expand mvec2 to n1
      mvec2 = vec2mat(mvec2,n1) ;
      n2 = n1 ;
    end ;


    if n1 == n2 ;     %  then OK, go ahead

      if iinorm == 1 ;    %  then do L_1 norm

        distout = abs(mvec1 - mvec2) ;
        distout = sum(distout,1) ;
          %  take sum along dimension 1

      elseif iinorm == 2 ;    %  then do L_2 (Euclidean) norm

        distout = (mvec1 - mvec2).^2 ;
        distout = sqrt(sum(distout,1)) ;
          %  take sum along dimension 1

      elseif iinorm == inf ;    %  then do L_inf (sup) norm

        distout = abs(mvec1 - mvec2) ;
        distout = max(distout,[],1) ;
          %  take max along dimension 1
          %    (ensures right answer even when d = 1)

      end ;


    else ;   %  No good, give error message

      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Error from VecDist:              !!!') ;
      disp('!!!   Inputs need same number of cols  !!!') ;
      disp('!!!   Giving Empty Return              !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

      distout = [] ;

    end ;


  else ;   %  then give error message

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Error from VecDist:              !!!') ;
    disp('!!!   Inputs need same number of rows  !!!') ;
    disp('!!!   Giving Empty Return              !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

    distout = [] ;

  end ;

end ;
