weighted.c

#define MIN(a, b) ((a) < (b) ? (a) : (b))

typedef struct dt_weighter {
	double (*insert)(char);
	double (*delete)(char);
	double (*replace)(char, char);
} dt_weighter;

int
dt_weighted(dt_weighter *weighter, const char *seq1, const char *seq2,
            size_t len1, size_t len2, double *dist)
{
	size_t i, j;
	double ic, dc, rc;
	double last, old, *column;
	
	assert(weighter != NULL);
	assert(weighter->insert && weighter->delete && weighter->replace);

	/* don't swap the sequences, or this is gonna be painful */
	
	if (len1 == 0 || len2 == 0) {
		*dist = 0.0;
		while (len1)
			*dist += weighter->delete(seq1[--len1]);
		while (len2)
			*dist += weighter->insert(seq2[--len2]);
		return 0;
	}

	column = malloc((len2 + 1) * sizeof(double));
	if (!column) return -1;

	column[0] = 0.;
	for (j = 1; j <= len2; ++j)
		column[j] = column[j - 1] + weighter->insert(seq2[j - 1]);
	
	for (i = 1; i <= len1; ++i) {
		last = column[0]; /* m[i-1][0] */
		column[0] += weighter->delete(seq1[i - 1]); /* m[i][0] */
		for (j = 1; j <= len2; ++j) {
			old = column[j];
			if (seq1[i - 1] == seq2[j - 1]) {
				column[j] = last; /* m[i-1][j-1] */
			} else {
				ic = column[j - 1] + weighter->insert(seq2[j - 1]);      /* m[i][j-1] */
				dc = column[j] + weighter->delete(seq1[i - 1]);          /* m[i-1][j] */
				rc = last + weighter->replace(seq1[i - 1], seq2[j - 1]); /* m[i-1][j-1] */
				column[j] = MIN(ic, MIN(dc, rc));
			}
			last = old;
		}
	}

	*dist = column[len2];
	free(column);
	return 0;
}

#if 1

static double insert(char c) { return 1.; }
static double delete(char c) { return 0.5; }
static double replace(char c, char d) { return 0.3; }

int main(int argc, char **argv)
{

	char *seq1 = argv[1];
	char *seq2 = argv[2];
	size_t len1 = strlen(seq1);
	size_t len2 = strlen(seq2);

	dt_weighter weighter;
	weighter.insert = insert;
	weighter.delete = delete;
	weighter.replace = replace;
	
	double dist;
	
	dt_weighted(&weighter, seq1, seq2, len1, len2, &dist);
	printf("%s -> %s -> %g\n", seq1, seq2, dist);
	
	return EXIT_SUCCESS;
}

#endif