Branch data Line data Source code
1 : : /** @file editdistance.cc
2 : : * @brief Edit distance calculation algorithm.
3 : : *
4 : : * Based on that described in:
5 : : *
6 : : * "An extension of Ukkonen's enhanced dynamic programming ASM algorithm"
7 : : * by Hal Berghel, University of Arkansas
8 : : * and David Roach, Acxiom Corporation
9 : : *
10 : : * http://berghel.net/publications/asm/asm.php
11 : : */
12 : : /* Copyright (C) 2003 Richard Boulton
13 : : * Copyright (C) 2007,2008,2009 Olly Betts
14 : : *
15 : : * This program is free software; you can redistribute it and/or modify
16 : : * it under the terms of the GNU General Public License as published by
17 : : * the Free Software Foundation; either version 2 of the License, or
18 : : * (at your option) any later version.
19 : : *
20 : : * This program is distributed in the hope that it will be useful,
21 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : : * GNU General Public License for more details.
24 : : *
25 : : * You should have received a copy of the GNU General Public License
26 : : * along with this program; if not, write to the Free Software
27 : : * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
28 : : */
29 : :
30 : : #include <config.h>
31 : :
32 : : #include "editdistance.h"
33 : :
34 : : #include "omassert.h"
35 : :
36 : : #include <algorithm>
37 : : #include <cstdlib>
38 : :
39 : : using namespace std;
40 : :
41 : : template<class CHR>
42 : : struct edist_seq {
43 : 720 : edist_seq(const CHR * ptr_, int len_) : ptr(ptr_), len(len_) { }
44 : : const CHR * ptr;
45 : : int len;
46 : : };
47 : :
48 : : template<class CHR>
49 : : class edist_state {
50 : : /// Don't allow assignment.
51 : : void operator=(const edist_state &);
52 : :
53 : : /// Don't allow copying.
54 : : edist_state(const edist_state &);
55 : :
56 : : edist_seq<CHR> seq1;
57 : : edist_seq<CHR> seq2;
58 : :
59 : : /* Array of f(k,p) values, where f(k,p) = the largest index i such that
60 : : * d(i,j) = p and d(i,j) is on diagonal k.
61 : : * ie: f(k,p) = largest i s.t. d(i,k+i) = p
62 : : * Where: d(i,j) = edit distance between substrings of length i and j.
63 : : */
64 : : int * fkp;
65 : : int fkp_cols;
66 : :
67 : : /* Maximum possible edit distance (this is referred to as ZERO_K in
68 : : * the algorithm description by Berghel and Roach). */
69 : : int maxdist;
70 : :
71 : 21412 : int calc_index(int k, int p) const {
72 : 21412 : return (k + maxdist) * fkp_cols + p + 1;
73 : : }
74 : :
75 : : public:
76 : :
77 : : edist_state(const CHR * ptr1, int len1, const CHR * ptr2, int len2);
78 : :
79 : : ~edist_state();
80 : :
81 : 3360 : int get_f_kp(int k, int p) const {
82 : 3360 : return fkp[calc_index(k, p)];
83 : : }
84 : :
85 : 18052 : void set_f_kp(int k, int p, int val) {
86 : 18052 : fkp[calc_index(k, p)] = val;
87 : 18052 : }
88 : :
89 : 934 : bool is_transposed(int pos1, int pos2) const {
90 [ + + ][ + + ]: 934 : if (pos1 <= 0 || pos2 <= 0 || pos1 >= seq1.len || pos2 >= seq2.len) return false;
[ + + ][ - + ]
91 : : return (seq1.ptr[pos1 - 1] == seq2.ptr[pos2] &&
92 [ + + ][ + + ]: 934 : seq1.ptr[pos1] == seq2.ptr[pos2 - 1]);
93 : : }
94 : :
95 : : void edist_calc_f_kp(int k, int p);
96 : : };
97 : :
98 : : template<class CHR>
99 : 934 : void edist_state<CHR>::edist_calc_f_kp(int k, int p)
100 : : {
101 : 934 : int maxlen = get_f_kp(k, p - 1) + 1; /* dist if do substitute */
102 : 934 : int maxlen2 = get_f_kp(k - 1, p - 1); /* dist if do insert */
103 : 934 : int maxlen3 = get_f_kp(k + 1, p - 1) + 1; /* dist if delete */
104 : :
105 [ + + ]: 934 : if (is_transposed(maxlen, maxlen + k)) {
106 : : // Transposition.
107 : 60 : ++maxlen;
108 : : }
109 : :
110 [ + + ]: 934 : if (maxlen >= maxlen2) {
111 [ + + ]: 707 : if (maxlen >= maxlen3) {
112 : : // Transposition or Substitution.
113 : : } else {
114 : : // Deletion.
115 : 9 : maxlen = maxlen3;
116 : : }
117 : : } else {
118 [ + - ]: 227 : if (maxlen2 >= maxlen3) {
119 : : // Insertion.
120 : 227 : maxlen = maxlen2;
121 : : } else {
122 : : // Deletion.
123 : 0 : maxlen = maxlen3;
124 : : }
125 : : }
126 : :
127 : : /* Check for exact matches, and increase the length until we don't have
128 : : * one. */
129 [ + + ][ + - ]: 2292 : while (maxlen < seq1.len &&
[ + + ][ + + ]
130 : : maxlen + k < seq2.len &&
131 : : seq1.ptr[maxlen] == seq2.ptr[maxlen + k]) {
132 : 1358 : ++maxlen;
133 : : }
134 : 934 : set_f_kp(k, p, maxlen);
135 : 934 : }
136 : :
137 : : #define INF 1000000
138 : : template<class CHR>
139 : 360 : edist_state<CHR>::edist_state(const CHR * ptr1, int len1,
140 : : const CHR * ptr2, int len2)
141 : 360 : : seq1(ptr1, len1), seq2(ptr2, len2), maxdist(len2)
142 : : {
143 : : Assert(len2 >= len1);
144 : : /* Each row represents a value of k, from -maxdist to maxdist. */
145 : 360 : int fkp_rows = maxdist * 2 + 1;
146 : : /* Each column represents a value of p, from -1 to maxdist. */
147 : 360 : fkp_cols = maxdist + 2;
148 : : /* fkp is stored as a rectangular array, row by row. */
149 : 360 : fkp = new int[fkp_rows * fkp_cols];
150 : :
151 [ + + ]: 4528 : for (int k = -maxdist; k <= maxdist; k++) {
152 [ + + ]: 36500 : for (int p = -1; p <= maxdist; p++) {
153 [ + + ]: 32332 : if (p == abs(k) - 1) {
154 [ + + ]: 4168 : if (k < 0) {
155 : 1904 : set_f_kp(k, p, abs(k) - 1);
156 : : } else {
157 : 2264 : set_f_kp(k, p, -1);
158 : : }
159 [ + + ]: 28164 : } else if (p < abs(k)) {
160 : 12950 : set_f_kp(k, p, -INF);
161 : : }
162 : : }
163 : : }
164 : 360 : }
165 : :
166 : : template<class CHR>
167 : 360 : edist_state<CHR>::~edist_state() {
168 [ + - ]: 360 : delete [] fkp;
169 : 360 : }
170 : :
171 : : template<class CHR>
172 : : static int
173 : 360 : seqcmp_editdist(const CHR * ptr1, int len1, const CHR * ptr2, int len2,
174 : : int max_distance)
175 : : {
176 : 360 : int lendiff = len2 - len1;
177 : : /* Make sure second sequence is longer (or same length). */
178 [ + + ]: 360 : if (lendiff < 0) {
179 : 112 : lendiff = -lendiff;
180 : 112 : swap(ptr1, ptr2);
181 : 112 : swap(len1, len2);
182 : : }
183 : :
184 : : /* Special case for if one or both sequences are empty. */
185 [ - + ]: 360 : if (len1 == 0) return len2;
186 : :
187 : 360 : edist_state<CHR> state(ptr1, len1, ptr2, len2);
188 : :
189 : 360 : int p = lendiff; /* This is the minimum possible edit distance. */
190 [ + + ]: 573 : while (p <= max_distance) {
191 [ + + ]: 1117 : for (int temp_p = 0; temp_p != p; ++temp_p) {
192 : 559 : int inc = p - temp_p;
193 [ + + ]: 559 : if (abs(lendiff - inc) <= temp_p) {
194 : 361 : state.edist_calc_f_kp(lendiff - inc, temp_p);
195 : : }
196 [ + + ]: 559 : if (abs(lendiff + inc) <= temp_p) {
197 : 15 : state.edist_calc_f_kp(lendiff + inc, temp_p);
198 : : }
199 : : }
200 : 558 : state.edist_calc_f_kp(lendiff, p);
201 : :
202 [ + + ]: 558 : if (state.get_f_kp(lendiff, p) == len1) break;
203 : 213 : ++p;
204 : : }
205 : :
206 : 360 : return p;
207 : : }
208 : :
209 : : int
210 : 360 : edit_distance_unsigned(const unsigned * ptr1, int len1,
211 : : const unsigned * ptr2, int len2,
212 : : int max_distance)
213 : : {
214 : 360 : return seqcmp_editdist<unsigned>(ptr1, len1, ptr2, len2, max_distance);
215 : : }
|
