## Copyright (C) 2006 Muthiah Annamalai <muthiah.annamalai@uta.edu> and Max Görner
##
## 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},@var{L}] =} editdistance (@var{string1}, @var{string2}, @var{weights}, @var{modus})
## Compute the Levenshtein edit distance between the strings @var{string1} and
## @var{string2}. This operation is symmetrical.
##
## 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.
##
## 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 are 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) + s) , where n and m denote the
## length of both strings.
## 
## The optional argument @var{modus} specifies the algorithm to be used. For
## @var{modus} = 0, Ukkonen's algorithm will be used whenever possible. For
## @var{modus} = 1, Ukkonens's algorithm will not be used. If @var{L} is omitted,
## and Ukkonens algorithm disabled, an algorithm using only an linear amount of
## memory with respect to the input length but O(m*n) runtime will be used,
## where, again, 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.
##
## @example
## @group
##          editdistance('marry','marie') 
##          ##returns value +2 for the distance.
## @end group
## @end example
##
## @end deftypefn

function [dist, L] = editdistance (str1, str2, weights, modus = 0)
  if (nargin < 2 || (nargin >= 3 && length (weights) < 3) )
    print_usage ();
  endif
  
  if (modus == 0 &&
      nargout < 2 && 
      (nargin >= 2 || weights(1) == weights(3) && weights(1) == 1)
     )
    dist = berghel_roach (str1,str2);
    return;
  endif

  saveMemory = nargout < 2;

  L1 = length (str1) + 1;
  L2 = length (str2) + 1;

  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
  
  if (nargin < 3)
    g = 1;%insertion
    m = 0;%match
    d = 1;%deletion
  else
    g = weights(1);
    m = weights(2);
    d = weights(3);
  endif

  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,min (m2,m3));
    endfor
    if (saveMemory)
      L(1, :) = L(2, :);
    endif
  endfor

  if (saveMemory)
    x = 2;
  else
    x = L1;
  endif

  dist = L(x,L2);
endfunction

function dist = berghel_roach (a,b)
  #Variable named according to Berghel and Roach 1996 (except s, which is called dist here)

  if (length (a) > length (b))
    ans = a;
    a = b;
    b = ans;
  endif
  m = length (a); n = length (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
