levenshtein.c

/*
 * @author Juan Miguel Cejuela, jmcejuela
 *
 *
 * INSTALLATION
 *
 * gcc -o levenshtein.so -shared levenshtein.c -I /usr/include/mysql/
 * (mac) gcc -bundle -o levenshtein.so levenshtein.c -I/usr/local/mysql/include
 *
 * Put the shared library as described in: http://dev.mysql.com/doc/refman/5.0/en/udf-compiling.html
 *
 * Afterwards in SQL:
 *
 * CREATE FUNCTION levenshtein RETURNS INT SONAME 'levenshtein.so';
 * CREATE FUNCTION levenshtein_k RETURNS INT SONAME 'levenshtein.so';
 *
 *
 * Some Credit for simple levenshtein to: Joshua Drew, SpinWeb Net Designs
 *
 */

#ifdef STANDARD
/* STANDARD is defined, don't use any mysql functions */
#include <string.h>
#ifdef __WIN__
typedef unsigned __int64 ulonglong;	/* Microsofts 64 bit types */
typedef __int64 longlong;
#else
typedef unsigned long long ulonglong;
typedef long long longlong;
#endif /*__WIN__*/
#else
#include <my_global.h>
#include <my_sys.h>
#if defined(MYSQL_SERVER)
#include <m_string.h>
#else
/* when compiled as standalone */
#include <string.h>
#endif
#endif
#include <mysql.h>
#include <ctype.h>


#ifdef HAVE_DLOPEN


/*
 * See MySQL UDF documentation pages for details on their implementation.
 *
 */


/**
 * Levenshtein distance
 *
 * @param s string 1 to compare, length n
 * @param t string 2 to compare, length m
 * @result levenshtein distance between s and t
 *
 * @time O(nm), quadratic
 * @space O(nm)
 */

my_bool		levenshtein_init(UDF_INIT *initid, UDF_ARGS *args, char *message);
void		levenshtein_deinit(UDF_INIT *initid);
longlong 	levenshtein(UDF_INIT *initid, UDF_ARGS *args, char *is_null, char *error);



/**
 * Levenshtein distance with threshold k (maximum allowed distance)
 *
 * @param s string 1 to compare, length n
 * @param t string 2 to compare, length m
 * @param k maximum threshold
 * @result levenshtein distance between s and t or >k (not specified) if this is greater than k
 *
 * @time O(kl), linear; where l = min(n, m)
 * @space O(k), constant
 */

my_bool 	levenshtein_k_init(UDF_INIT *initid, UDF_ARGS *args, char *message);
void		levenshtein_k_deinit(UDF_INIT *initid);
longlong 	levenshtein_k(UDF_INIT *initid, UDF_ARGS *args, char *is_null, char *error);




//(Expected) maximum number of digits to return
#define LEVENSHTEIN_MAX 3




my_bool levenshtein_init(UDF_INIT *initid, UDF_ARGS *args, char *message) {
  if ((args->arg_count != 2) ||
      (args->arg_type[0] != STRING_RESULT || args->arg_type[1] != STRING_RESULT)) {
    strcpy(message, "Function requires 2 arguments, (string, string)");
    return 1;
  }

  //matrix for levenshtein calculations of size n+1 x m+1 (+1 for base values)
  int *d = (int *) malloc(sizeof(int) * (args->lengths[0] + 1) * (args->lengths[1] + 1));
  if (d == NULL) {
    strcpy(message, "Failed to allocate memory");
    return 1;
  }

  initid->ptr = (char*) d;
  initid->max_length = LEVENSHTEIN_MAX;
  initid->maybe_null = 0; //doesn't return null

  return 0;
}


void levenshtein_deinit(UDF_INIT *initid) {
  if (initid->ptr != NULL) {
    free(initid->ptr);
  }
}


inline int minimum(int a, int b, int c) {
  int min = a;
  if (b < min)
    min = b;
  if (c < min)
    min = c;
  return min;
}


longlong levenshtein(UDF_INIT *initid, UDF_ARGS *args, char *is_null, char *error) {
  const char *s = args->args[0];
  const char *t = args->args[1];

  int n = (s == NULL) ? 0 : args->lengths[0];
  int m = (t == NULL) ? 0 : args->lengths[1];

  if (0 == n)
    return m;
  if (0 == m)
    return n;

  int *d = (int*) initid->ptr;

  /* Initialization */
  n++; m++;

  int i;
  for (i = 0; i < n; i++)
    d[i] = i;

  int j;
  for (j = 0; j < m; j++)
    d[n * j] = j;


  /* Recurrence */

  int p, h; //indices for d matrix seen as a vector, see below

  int im1 = 0; //i minus 1
  for (i = 1; i < n; i++) {
    p = i;
    int jm1 = 0; //j minus 1
    for (j = 1; j < m; j++) {
      h = p; // d[h] = d[i,j-1], h = (j * n + i - n)  = ((j - 1) * n + i)
      p += n; // d[p] = d[i,j], p = (j * n + i)

      if (s[im1] == t[jm1]) {
	d[p] = d[h-1]; //no operation required, d[i-1,j-1]
      }

      else {
	d[p] = minimum(d[p-1], 	//deletion, d[i-1, j]
		       d[h], 	//insertion, d[i, j-1]
		       d[h-1]	//substitution, d[i-1,j-1]
		       ) + 1; 	//can put +1 outside because the cost is the same
      }

      jm1 = j;
    }

    im1 = i;
  }

  return (longlong) d[p];
}





my_bool levenshtein_k_init(UDF_INIT *initid, UDF_ARGS *args, char *message) {
  if ((args->arg_count != 3) ||
      (args->arg_type[0] != STRING_RESULT || args->arg_type[1] != STRING_RESULT || args->arg_type[2] != INT_RESULT)) {
    strcpy(message, "Function requires 3 arguments, (string, string, int)");
    return 1;
  }

  initid->max_length = LEVENSHTEIN_MAX;
  initid->maybe_null = 0; //doesn't return null

  return 0;
}


//Not necessary
//
/* void levenshtein_k_deinit(UDF_INIT *initid) { */
/*   //nothing */
/*   fflush(stderr); //for debugging */
/* } */



/*
 * 1st observation: time O(kl)
 * (see Algorithms on Strings, Trees, and Sequences, Dan Gusfield, pg. 263)
 *
 * When the levenshtein distance is limited to k the alignment path cannot be >k
 * cells off the main diagonal, either from the left (insertions) or from the
 * right (deletions) This means we don't have to fill the whole recurrence
 * matrix; it suffices to fill a strip of 2k + 1 cells in a row. Also, r <= k,
 * for r = m - n, is a necessary condition for there to be any solution.
 *
 *
 * 2nd observation: exactly k + 1
 *
 * The number of _unpaired_ insertions/deletions (the absolute difference in the
 * number of these operations) is as maximum k. When n == m (r == 0) any
 * movement to the left, insertion, in the recurrence matrix is corresponded to
 * a movement to the right, deletion, or viceversa. Analogously when m > n only
 * a maximum of k insertions/deletions operations can be left unpaired. This
 * makes the required strip size to be k + 1 cells.
 *
 * Call l = min(n, m). Place the string with length l on the y axis, and the
 * other on the x axis. By doing this, the number of rows is equal or less than
 * the number of columns and so the algorithm effectively runs in O(kl) Also, by
 * this arrangement of s and t in the recurrence matrix, r indicates the forced
 * number of unpaired deletions, i.e., since r>=0, there cannot be unpaired
 * deletions. Observe that for there to be any insertion k must be equal or
 * greater than 2.
 *
 * This arrangement leaves (k-r)/2 possible insertions (left) and (k-r)/2 + r
 * possible deletions (right) --> (k-r)/2 + (k-r)/2 + r + 1 == k + 1
 *
 *
 * 3rd observation: space O(k)
 *
 * As with the original levenshtein algorithm, only the current and last rows
 * are needed (that is, when we don't do traceback because we need not the
 * alignment path). By the previous observation, we just need 2 rows of k + 1
 * cells. For this, considering the last row, the previous diagonal is actually
 * in the same column (substitution) and the above cell is actually in the next
 * column (deletion). Observe that by doing this, we're virtually creating a
 * recurrence matrix of n * (k+1) being now the original diagonal in the middle
 * column (-r) (matrix which could be used to do the traceback)
 *
 */



longlong levenshtein_k(UDF_INIT *initid, UDF_ARGS *args, char *is_null, char *error) {
  char *s = args->args[0];
  char *t = args->args[1];

  int n = (s == NULL) ? 0 : args->lengths[0];
  int m = (t == NULL) ? 0 : args->lengths[1];

  //order the strings so that the first always has the minimum length l
  if (n > m) {
    int aux = n;
    n = m;
    m = aux;
    char *auxs = s;
    s = t;
    t = auxs;
  }

  const int k = *((int*) args->args[2]);
  const int ignore = k + 1; //lev dist between s and t is at least greater than k
  const int r = m - n;

  if (0 == n)
    return (m > k) ? ignore : m;
  if (0 == m)
    return (n > k) ? ignore : n;
  if (r > k)
    return ignore;

  const int lsize = (((k > m) ? m : k) - r) / 2; //left space for insertions
  const int rsize = lsize + r; //right space for deletions, rsize >= lsize (rsize == lsize iff r == 0)
  const int stripsize = lsize + rsize + 1; // + 1 for the diagonal cell
  const int stripsizem1 = stripsize - 1; //see later, not to repeat calculations

  int d[2 * stripsize]; //Current and last rows
  int currentrow;
  int lastrow;


  /* Initialization */

  //currentrow = 0
  int i;
  for (i = lsize; i < stripsize; i++) //start from diagonal cell
    d[i] = i - lsize;


  /* Recurrence */

  currentrow = stripsize;
  lastrow = 0;

  //j index for virtual recurrence matrix, jv index for rows
  //bl & br = left & right bounds for j
  int j, jv, bl, br;
  int im1 = 0, jm1;
  int a, b, c, min; //for minimum function, coded directly here for maximum speed
  for (i = 1; i <= n; i++) {

    //bl = max(i - lsize, 0), br = min(i + rsize, m)
    bl = i - lsize;
    if (bl < 0) {
      jv = abs(bl); //no space for all allowed insertions
      bl = 0;
    }
    else
      jv = 0;
    br = i + rsize;
    if (br > m)
      br = m;

    jm1 = bl - 1;
    for (j = bl; j <= br; j++) {
      if (0 == j) //postponed part of initialization
  	d[currentrow + jv] = i;
      else {
	//By observation 3, the indices change for the lastrow (always +1)
  	if (s[im1] == t[jm1]) {
  	  d[currentrow + jv] = d[lastrow + jv];
	}
  	else {
	  //get the minimum of these 3 operations
	  a = (0 == jv) ? ignore : d[currentrow + jv - 1]; //deletion
	  b = (stripsizem1 == jv) ? ignore : d[lastrow + jv + 1]; //insertion
	  c = d[lastrow + jv]; //substitution

	  min = a;
	  if (b < min)
	    min = b;
	  if (c < min)
	    min = c;

	  d[currentrow + jv] = min + 1;
  	}
      }
      jv++;
      jm1 = j;
    }

    //obsv: the cost of a following diagonal never decreases
    if (d[currentrow + lsize + r] > k)
      return ignore;

    im1 = i;

    //swap
    currentrow = currentrow ^ stripsize;
    lastrow = lastrow ^ stripsize;
  }

  //only here if levenhstein(s, t) <= k
  return (longlong) d[lastrow + lsize + r]; //d[n, m]
}



#endif /* HAVE_DLOPEN */