# HG changeset patch
# User Max Görner
# Date 1418594611 -3600
#      Sun Dec 14 23:03:31 2014 +0100
# Node ID c9c5a15528a7ceb6fabba79ef08c7b15a8336bf5
# Parent  e9c2e6fd23e130d79167a2dbd43604d26a35844b
editdistance: use Berghel and Roach algorithm for performance (patch #8426)

* editdistance.m: This file got almost rewritten, so that it takes some new
optional arguments, checks whether the input is suiting and which algorithm to
use and then calculates the editdistance. The algorithms used are the one by
Berghel and Roach (BaR), the one by Fisher and Wagner (FaW) and FaW but with
linear memory usage.

* editdistance_berghel_roach.m: This file was added to be a helper function of
editdistance.m but may be useful on its own. It calculates the editdistance by
using an algorithm by Berghel and Roach
(http://www.berghel.net/publications/asm/asm.php). Therefor it is possible to
calculate 'editdistance(repmat('AB',1,20000),repmat('BA',1,20000)' in less
than two seconds now.

diff -r e9c2e6fd23e1 -r c9c5a15528a7 inst/editdistance.m
--- a/inst/editdistance.m	Sun Jan 19 10:19:12 2014 +0100
+++ b/inst/editdistance.m	Sun Dec 14 23:03:31 2014 +0100
@@ -1,4 +1,5 @@
 ## Copyright (C) 2006 Muthiah Annamalai <muthiah.annamalai@uta.edu>
+## Copyright (C) 2014 Max Görner <max.goerner@tu-dresden.de>
 ##
 ## This program is free software; you can redistribute it and/or modify it under
 ## the terms of the GNU General Public License as published by the Free Software
@@ -14,20 +15,37 @@
 ## this program; if not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{dist},@var{L}] =} editdistance (@var{string1}, @var{string2}, @var{weights})
-## Compute the Levenshtein edit distance between the strings @var{string1} and
-## @var{string2}. This operation is symmetrical.
+## @deftypefn  {Function File} {[@var{dist}, @var{L}] =} editdistance (@var{str1}, @var{str2})
+## @deftypefnx {Function File} {[@var{dist}, @var{L}] =} editdistance (@var{str1}, @var{str2}, @var{weights})
+## @deftypefnx {Function File} {[@var{dist}, @var{L}] =} editdistance (@var{str1}, @var{str2}, @var{weights}, @var{modus})
 ##
-## The optional argument @var{weights} specifies weights for the
-## deletion, matched, and insertion operations; by default it is set to
-## +1, 0, +1 respectively, so that a least editdistance means a 
-## closer match between the two strings. This function implements
-## the Levenshtein edit distance as presented in Wikipedia article,
-## accessed Nov 2006. Also the levenshtein edit distance of a string
-## with an empty string is defined to be its length.
+## Compute the Levenshtein edit distance between the two strings.
+##
+## The optional argument @var{weights} specifies weights for the deletion,
+## matched, and insertion operations; by default it is set to +1, 0, +1
+## respectively, so that a least editdistance means a closer match between the
+## two strings. This function implements the Levenshtein edit distance as
+## presented in Wikipedia article, accessed Nov 2006. Also the levenshtein edit
+## distance of a string with the empty string is defined to be its length.
+##
+## For the special case that there are no weights given and the array L is not
+## requested, an algorithm of Berghel and Roach, which improves an algorithm 
+## introduced by Ukkonen in 1985, will be applied. This algorithm is
+## significantly faster most of the times. Its main strength lies in cases with
+## small edit distances, where huge speedups and memory savings are suspectible.
+## The time (and space) complexity is O(((dist^2 - (n - m)^2)/2) + dist), where
+## n and m denote the length of both strings.
 ## 
-## The default return value is @var{dist} the edit distance, and
-## the other return value  @var{L} is the distance matrix.
+## The optional argument @var{modus} specifies the algorithm to be used. For
+## @var{modus} = 0, Berghel and Roach's algorithm will be used whenever
+## possible. For @var{modus} = 1, the classic algorithm by Fisher and Wagner
+## will be used. If @var{L} is omitted, and @var{modus} = 1, a variant of Fisher
+## and Wagner's algorithm using only a linear amount of memory with respect to
+## the input length, but O(m*n) runtime, will be used. Again, n and m denote the
+## length of both strings.
+##
+## The default return value @var{dist} is the edit distance, and
+## the other return value @var{L} is the distance matrix.
 ##
 ## @example
 ## @group
@@ -38,42 +56,116 @@
 ##
 ## @end deftypefn
 
-function [dist, L] = editdistance (str1, str2, weights)
-  if(nargin < 2 || (nargin == 3 && length(weights)  < 3) )
-    print_usage();
-  end
+function [dist, L] = editdistance (str1, str2, weights = [1 0 1], modus = 0)
+  #Checking correct call of the function.
+  if (nargin < 2 || nargin > 4)
+    print_usage ();
+  endif
+  if (nargin >= 3 && numel (weights) != 3)
+    error ("editdistance: WEIGTHS must be a 3 element vector.")
+  endif 
+  if (!isvector(str1) || !isvector(str2))
+    error ("editdistance: Both strings must be a vector.")
+  endif
+
+  #Using the approach of Berghel and Roach, if possible.
+  if (modus == 0 &&
+      nargout < 2 && 
+      (weights(1) == weights(3) && weights(2) == 0)
+     )
+    dist = weights(1)*editdistance_berghel_roach (str1,str2);
+    return;
+  endif
+
+  saveMemory = nargout < 2;
+
+  L1 = numel (str1) + 1;
+  L2 = numel (str2) + 1;
+
+  #Checking whether str1 or str2 are empty strings. If so, the determination of the minimal edit distance is trivial.
+  if (L1 == 1)
+    dist = L2-1;
+    return;
+  elseif (L2 == 1)
+    dist = L1-1;
+    return;
+  endif
+
+  if (saveMemory)
+    L = zeros (2,L2);
+  else
+    L  = zeros (L1,L2);
+  endif
   
-  L1=length(str1)+1;
-  L2=length(str2)+1;
-  L=zeros(L1,L2);
-  
-  if(nargin < 3)
-    g=+1;%insertion
-    m=+0;%match
-    d=+1;%deletion
+  %Setting weights
+  g = weights(1); %insertion
+  m = weights(2); %match
+  d = weights(3); %deletion
+
+  if (not (saveMemory))
+    L(:,1) = [0:L1-1]'*g;
+  endif
+  L(1,:) = [0:L2-1]*g;
+
+  for idx = 2:L1;
+    if (saveMemory)
+      L(2, 1) = idx-1;
+    endif
+    for idy = 2:L2
+      if (str1(idx-1) == str2(idy-1))
+        score = m;
+      else
+        score = d;
+      endif
+      if (saveMemory)
+        x = 2;
+      else
+        x = idx;
+      endif
+      m1 = L(x-1,idy-1) + score;
+      m2 = L(x-1,idy) + g;
+      m3 = L(x,idy-1) + g;
+      L(x,idy) = min ([m1, m2, m3]);
+    endfor
+    if (saveMemory)
+      L(1, :) = L(2, :);
+    endif
+  endfor
+
+  if (saveMemory)
+    x = 2;
   else
-    g=weights(1);
-    m=weights(2);
-    d=weights(3);
-  end
+    x = L1;
+  endif
 
-  L(:,1)=[0:L1-1]'*g;
-  L(1,:)=[0:L2-1]*g;
-  
-  m4=0;
-  for idx=2:L1;
-    for idy=2:L2
-      if(str1(idx-1)==str2(idy-1))
-        score=m;
-      else
-        score=d;
-      end
-      m1=L(idx-1,idy-1) + score;
-      m2=L(idx-1,idy) + g;
-      m3=L(idx,idy-1) + g;
-      L(idx,idy)=min(m1,min(m2,m3));
-    end
-  end
-  
-  dist=L(L1,L2);
+  dist = L(x,L2);
 endfunction
+
+%!test
+%! l = 50;
+%! n = 20;
+%! rand('state',31513);
+%! abc = 'A':'Z';
+
+%! for it = 1:n
+%!   #Generate two Strings
+%!   #This kind of generation produces worst case examples for the new algorithm,
+%!   #so runtime comparisons won't be informative.
+%!   str1 = str2 = abc(randi([1 length(abc)],1,l));
+%!   m = randi(l/2);
+%!   str1(randi([1 l],1,m)) = abc(randi([1 length(abc)],1,m));
+%!   m = randi(l/2);
+%!   str2(randi([1 l],1,m)) = abc(randi([1 length(abc)],1,m));
+
+%!   d1 = editdistance(str1,str2); #Berghel and Roach
+%!   tempDist = editdistance(str2,str1);
+%!   assert(d1 == tempDist, "editdistance: Berghel and Roach algorithm is not symmetric.");
+%!   [d2 ~] = editdistance(str1,str2); #Fisher and Wagner
+%!   [tempDist ~] = editdistance(str2,str1);
+%!   assert(d2 == tempDist, "editdistance: Fisher and Wagner algorithm is not symmetric.");
+%!   d3 = editdistance(str1,str2,[1 0 1],1); #Fisher and Wagner with linear memory usage
+%!   tempDist = editdistance(str2,str1,[1 0 1],1);
+%!   assert(d3 == tempDist, "editdistance: Fisher and Wagner algorithm with linear memory usage is not symmetric.");
+%!   assert(d1 == d2, "editdistance: The result of the algorithm by Berghel and Roach and that of the algorithm by Fisher and Wagners differ.");
+%!   assert(d2 == d3, "editdistance: The results of the algorithm by Fisher and Wagner with and without linear memory usage differ.");
+%! endfor
diff -r e9c2e6fd23e1 -r c9c5a15528a7 inst/editdistance_berghel_roach.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/editdistance_berghel_roach.m	Sun Dec 14 23:03:31 2014 +0100
@@ -0,0 +1,113 @@
+## Copyright (C) 2014 Max Görner <max.goerner@tu-dresden.de>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{dist} =} editdistance_berghel_roach (@var{str1}, @var{str2})
+##
+## Compute the Levenshtein edit distance between the two strings using the approach of Berghel and Roach.
+##
+## This algorithm is significantly faster then the algorithm by Fisher and
+## Wagner most of the times. Its main strength lies in cases with small edit
+## distances, where huge speedups and memory savings are suspectible. The time
+## (and space) complexity is O(((dist^2 - (n - m)^2)/2) + dist), where n and m
+## denote the length of both strings.
+## 
+## The return value @var{dist} is the edit distance.
+##
+## @example
+## @group
+##          editdistance_berghel_roach('marry','marie') 
+##          ##returns value +2 for the distance.
+## @end group
+## @end example
+##
+## @end deftypefn
+
+function dist = editdistance_berghel_roach (a,b)
+  #Variables named according to Berghel and Roach 1996 (except s, which is called dist here)
+  
+  if (!isvector(a) || !isvector(b))
+    print_usage();
+  endif
+
+  if (numel (a) > numel (b))
+    ans = a;
+    a = b;
+    b = ans;
+  endif
+  m = numel (a); n = numel (b);
+  if (m == 0)
+    dist = n;
+    return;
+  endif
+
+  MIN_K = -m-1; #The diagonal with the lowest number is -m
+  MIN_P = -1; #minimum p that has to be cached is 0
+  FKP = sparse (m+n+2,n+1);
+  FKP_cached = sparse (m+n+2,n+1);
+
+  p = n-m;
+  do
+    inc = p;
+    for temp_p = 0:p-1
+      if (abs (n-m - inc) <= temp_p)
+        get_f ((n-m) - inc, temp_p);
+      endif
+      if (abs(n-m + inc) <= temp_p)
+        get_f(n-m+inc,temp_p);
+      endif
+      inc--;
+    endfor
+    get_f (n-m,p);
+    p++;
+  until (get_f (n-m,p-1) == m);
+
+  dist = p-1;
+
+  function f = get_f (k,p)
+    if (p == abs (k)-1)
+      if (k < 0)
+        f = abs (k) - 1;
+      else
+        f = -1;
+      endif
+      return;
+    endif
+    if (p < abs (k)-1)
+      f = -inf;
+      return;
+    endif
+    
+    if (FKP_cached(k-MIN_K,p-MIN_P) == 1)
+      f = FKP(k-MIN_K,p-MIN_P);
+      return;
+    endif
+
+    c1 = get_f (k,p-1) + 1;
+    c2 = get_f (k-1,p-1);
+    c3 = get_f (k+1,p-1) + 1;
+    t = max ([c1 c2 c3]);
+    #if (a([t, t+1]) == b([k+t+1, k+t]) t2 = t+1; endif #taking adjacent transpositions into account
+    while (t < m && t+k < n && a(t+1) == b(t+1+k))
+      t += 1;
+    endwhile
+    if (t > m || t+k > n)
+      f = FKP(k-MIN_K,p-MIN_P) = NaN;
+    else
+      f = FKP(k-MIN_K,p-MIN_P) = t;
+    endif
+    FKP_cached(k-MIN_K,p-MIN_P) = 1;
+  endfunction
+endfunction
