main.bend
bend-edit-distance@0.2.0.0 · 735 lines · 0 definitions
1import Base
2import bend-edit-distance@0.2.0.0/Cli.bend as CLI
3import bend-edit-distance@0.2.0.0/EdlibTypes.bend as E
4import bend-edit-distance@0.2.0.0/Alignment.bend as A
5import bend-edit-distance@0.2.0.0/AlignmentCli.bend as AC
6
7# The distance paths from vendor/edlib/edlib/src/edlib.cpp. Bend's native
8# integer is U32, so a Myers word covers 32 query characters instead of 64.
9# The scalar helper is NW, unit costs, unlimited distance (k = -1). The
10# configured API extends it with the typed modes/tasks in EdlibTypes.bend.
11
12type Peq is Data:
13 Mask{bits: U32}
14 Branch{low: Peq, high: Peq}
15
16type Block is Data:
17 Block{positive: U32, negative: U32, score: U32}
18
19# Horizontal differences are -1, 0, or +1. Keep both signs explicit rather
20# than relying on signed shifts or two's-complement arithmetic.
21type Carry is Data:
22 Carry{positive: U32, negative: U32}
23
24def split_query.prepend(h: Char, parts: String & String) -> String & String:
25 match parts:
26 case (prefix, suffix):
27 (SCon{h, prefix}, suffix)
28
29def split_query(n: Nat, query: String) -> String & String:
30 match n query:
31 case 0n _:
32 ("", query)
33 case 1n+p SNil{}:
34 ("", "")
35 case 1n+p SCon{h, t}:
36 split_query.prepend(h, split_query(p, t))
37
38# Bit i matches query[i] literally or through an explicit symmetric pair.
39# Unused high bits are wildcards regardless of the pairs, as in buildPeq.
40def build_peq_word(query: String, +symbol: Char, +bit: U32,
41 +pairs: List<&2, E.EqualityPair>) -> U32:
42 match query:
43 case SNil{}:
44 U32.not(U32.sub(bit, 1))
45 case SCon{h, t}:
46 (U32.mul(Bool.to_u32(E.equal(pairs, h, symbol)), bit) .|.
47 build_peq_word(t, symbol, U32.shl(bit), pairs) : U32)
48
49# A byte-indexed immutable table replaces the C++ transformed alphabet.
50# Low bits select branches, so lookup takes eight steps for every byte.
51# Both halves of the table have equal, independent work and run in parallel.
52def build_peq_table(+depth: Nat, +query: String, +symbol: U32, +step: U32,
53 +pairs: List<&2, E.EqualityPair>) -> Peq:
54 match depth:
55 case 0n:
56 Mask{build_peq_word(query, Chr{symbol}, 1, pairs)}
57 case 1n+p:
58 low high = build_peq_table(p, query, symbol, U32.shl(step), pairs) build_peq_table(p, query, U32.add(symbol, step), U32.shl(step), pairs)
59 Branch{low, high}
60
61def peq_get(table: Peq, +code: U32, bit: Bool) -> U32:
62 match table bit:
63 case Mask{bits} _:
64 bits
65 case Branch{low, high} True{}:
66 +next = U32.shr(code)
67 peq_get(high, next, U32.is_ne(U32.and(next, 1), 0))
68 case Branch{low, high} False{}:
69 +next = U32.shr(code)
70 peq_get(low, next, U32.is_ne(U32.and(next, 1), 0))
71
72def build_peq.step(parts: String & String, pairs: List<&2, E.EqualityPair>,
73 tail: String -> List<&2, Peq>) -> List<&2, Peq>:
74 match parts:
75 case (word, rest):
76 build_peq_table(8n, word, 0, 1, pairs) <> tail(rest)
77
78def build_peq(fuel: Nat, query: String, +pairs: List<&2, E.EqualityPair>) -> List<&2, Peq>:
79 match fuel:
80 case 0n:
81 Nil{}
82 case 1n+p:
83 build_peq.step(split_query(32n, query), pairs, rest => build_peq(p, rest, pairs))
84
85def ceil_div(+x: Nat, +y: Nat) -> Nat:
86 Nat.div(Nat.add(x, Nat.sub(y, 1n)), y)
87
88def init_blocks(n: Nat, +bottom: U32) -> List<Block>:
89 match n:
90 case 0n:
91 Nil{}
92 case 1n+p:
93 Block{4294967295, 0, bottom} <> init_blocks(p, U32.add(bottom, 32))
94
95# calculateBlock / Advance_Block from Myers. Addition intentionally wraps
96# at 32 bits. Xv uses the original Eq, before the negative carry is added.
97def calculate_block(block: Block, +eq: U32, hin: Carry) -> Block & Carry:
98 Block{+pv, +mv, score} = block
99 Carry{hin_positive, +hin_negative} = hin
100 +xv = U32.or(eq, mv)
101 +eq_in = U32.or(eq, hin_negative)
102 +xh = (U32.xor(U32.add(U32.and(eq_in, pv), pv), pv) .|. eq_in : U32)
103 +ph = U32.or(mv, U32.not(U32.or(xh, pv)))
104 +mh = U32.and(pv, xh)
105 +hout_positive = U32.shrn(ph, 31n)
106 +hout_negative = U32.shrn(mh, 31n)
107 +ph_shifted = U32.or(U32.shl(ph), hin_positive)
108 mh_shifted = U32.or(U32.shl(mh), hin_negative)
109 next_pv = U32.or(mh_shifted, U32.not(U32.or(xv, ph_shifted)))
110 next_mv = U32.and(ph_shifted, xv)
111 (Block{next_pv, next_mv, U32.sub(U32.add(score, hout_positive), hout_negative)},
112 Carry{hout_positive, hout_negative})
113
114# A column and consecutive target columns depend on their predecessor;
115# parallelizing those calls would violate the Myers carry dependency.
116def advance_column.step(result: Block & Carry,
117 tail: Carry -> List<Block>) -> List<Block>:
118 match result:
119 case (next, hout):
120 next <> tail(hout)
121
122def advance_column(blocks: List<Block>, profiles: List<&2, Peq>, +symbol: U32,
123 hin: Carry) -> List<Block>:
124 match blocks profiles:
125 case Nil{} _:
126 Nil{}
127 case _ Nil{}:
128 Nil{}
129 case block <> rest profile <> tail:
130 advance_column.step(
131 calculate_block(block, peq_get(profile, symbol, U32.is_ne(U32.and(symbol, 1), 0)), hin),
132 hout => advance_column(rest, tail, symbol, hout))
133
134# getBlockCellValues(block)[W], without constructing the other 31 cells.
135# Walking up a row subtracts its positive vertical difference and adds its
136# negative difference. W is always in 0..31, never a whole extra block.
137def get_block_cell_value(padding: Nat, +pv: U32, +mv: U32, score: U32) -> U32:
138 match padding:
139 case 0n:
140 score
141 case 1n+p:
142 get_block_cell_value(p, U32.shl(pv), U32.shl(mv),
143 U32.add(U32.sub(score, U32.shrn(pv, 31n)), U32.shrn(mv, 31n)))
144
145def last_block_score(blocks: List<Block>, +padding: Nat) -> Nat:
146 match blocks:
147 case Nil{}:
148 0n
149 case Block{pv, mv, score} <> Nil{}:
150 U32.to_nat(get_block_cell_value(padding, pv, mv, score))
151 case _ <> rest:
152 last_block_score(rest, padding)
153
154# An active entry owns its current block and shares its immutable profile.
155# Dropped leading entries are never visited again. Trailing profiles remain
156# available for expansion; retries share the full profile but reset all blocks.
157type BandEntry is Data:
158 BandEntry{profile: Peq, block: Block}
159
160type Band is Type:
161 Band{first: Nat, entries: List<&2, BandEntry>, rest: List<&2, Peq>}
162
163type UpdatedBand is Type:
164 UpdatedBand{entries: List<&2, BandEntry>, score: U32, carry: Carry, bottom: Nat}
165
166type BandStep is Type:
167 BandStep{band: Band, maximum: Nat, work: Nat}
168
169# Work counts actual calculate_block calls, consumed columns and fresh attempts.
170# The probe uses these production counters to verify pruning and retry behavior.
171type Search is Data:
172 Search{distance: E.EditDistance, blocks: Nat, columns: Nat, attempts: Nat}
173
174def init_band_entries.step(parts: List<&2, BandEntry> & List<&2, Peq>,
175 profile: Peq, bottom: U32) -> List<&2, BandEntry> & List<&2, Peq>:
176 match parts:
177 case (entries, rest):
178 (BandEntry{profile, Block{4294967295, 0, bottom}} <> entries, rest)
179
180def init_band_entries(fuel: Nat, profiles: List<&2, Peq>,
181 +bottom: U32) -> List<&2, BandEntry> & List<&2, Peq>:
182 match fuel profiles:
183 case 0n _:
184 (Nil{}, profiles)
185 case _ Nil{}:
186 (Nil{}, Nil{})
187 case 1n+p profile <> tail:
188 init_band_entries.step(init_band_entries(p, tail, U32.add(bottom, 32)), profile, bottom)
189
190def init_band.step(parts: List<&2, BandEntry> & List<&2, Peq>) -> Band:
191 match parts:
192 case (entries, rest):
193 Band{0n, entries, rest}
194
195def init_band(profiles: List<&2, Peq>, +length: Nat, +target_length: Nat,
196 maximum: Nat) -> Band:
197 # Called only after k >= target_length - length, with length <= target_length.
198 count = Nat.min(ceil_div(length, 32n),
199 ceil_div(1n+Nat.div(Nat.sub(maximum, Nat.sub(target_length, length)), 2n), 32n))
200 init_band.step(init_band_entries(count, profiles, 32))
201
202def advance_band.join(entry: BandEntry, updated: UpdatedBand) -> UpdatedBand:
203 match updated:
204 case UpdatedBand{entries, score, carry, bottom}:
205 UpdatedBand{entry <> entries, score, carry, bottom}
206
207def advance_band.last(result: Block & Carry, profile: Peq, bottom: Nat) -> UpdatedBand:
208 match result:
209 case (Block{+pv, +mv, +score}, carry):
210 UpdatedBand{[BandEntry{profile, Block{pv, mv, score}}], score, carry, bottom}
211
212def advance_band.next(result: Block & Carry, profile: Peq,
213 tail: Carry -> UpdatedBand) -> UpdatedBand:
214 match result:
215 case (next, carry):
216 advance_band.join(BandEntry{profile, next}, tail(carry))
217
218def advance_band(entries: List<&2, BandEntry>, +symbol: U32,
219 hin: Carry, +bottom: Nat) -> UpdatedBand:
220 match entries:
221 case Nil{}:
222 UpdatedBand{Nil{}, 0, Carry{1, 0}, bottom}
223 case BandEntry{+profile, block} <> Nil{}:
224 advance_band.last(calculate_block(block,
225 peq_get(profile, symbol, U32.is_ne(U32.and(symbol, 1), 0)), hin), profile, bottom)
226 case BandEntry{+profile, block} <> tail:
227 advance_band.next(calculate_block(block,
228 peq_get(profile, symbol, U32.is_ne(U32.and(symbol, 1), 0)), hin), profile,
229 carry => advance_band(tail, symbol, carry, Nat.add(bottom, 32n)))
230
231# Edlib's conservative Ukkonen block bounds, rearranged so every term is
232# nonnegative Nat arithmetic. Rows/columns here count consumed characters.
233# The 62/63 constants adapt 2*WORD_SIZE-2/-1 from the pinned 64-bit code.
234def band_keep_last(+score: Nat, bottom: Nat, +maximum: Nat,
235 length: Nat, target_length: Nat, column: Nat) -> Bool:
236 Bool.and(Nat.is_lt(score, Nat.add(maximum, 32n)),
237 Nat.is_le(Nat.add(Nat.add(bottom, score), target_length),
238 Nat.add(Nat.add(Nat.add(maximum, length), column), 63n)))
239
240def band_keep_first(+score: Nat, bottom: Nat, +maximum: Nat,
241 length: Nat, target_length: Nat, column: Nat) -> Bool:
242 Bool.and(Nat.is_lt(score, Nat.add(maximum, 32n)),
243 Nat.is_ge(Nat.add(Nat.add(bottom, maximum), target_length),
244 Nat.add(Nat.add(score, length), column)))
245
246def trim_last.choose(entry: BandEntry, entries: List<&2, BandEntry>, rest: List<&2, Peq>,
247 keep: Bool) -> List<&2, BandEntry> & List<&2, Peq>:
248 match entry entries keep:
249 case BandEntry{profile, _} Nil{} False{}:
250 (Nil{}, profile <> rest)
251 case _ _ _:
252 (entry <> entries, rest)
253
254def trim_last.join(entry: BandEntry, tail: List<&2, BandEntry> & List<&2, Peq>,
255 keep: Bool) -> List<&2, BandEntry> & List<&2, Peq>:
256 match tail:
257 case (entries, rest):
258 trim_last.choose(entry, entries, rest, keep)
259
260def trim_last(entries: List<&2, BandEntry>, rest: List<&2, Peq>, +bottom: Nat,
261 +maximum: Nat, +length: Nat, +target_length: Nat,
262 +column: Nat) -> List<&2, BandEntry> & List<&2, Peq>:
263 match entries:
264 case Nil{}:
265 (Nil{}, rest)
266 case BandEntry{profile, Block{pv, mv, +score}} <> tail:
267 trim_last.join(BandEntry{profile, Block{pv, mv, score}},
268 trim_last(tail, rest, Nat.add(bottom, 32n), maximum, length, target_length, column),
269 band_keep_last(U32.to_nat(score), bottom, maximum, length, target_length, column))
270
271def trim_first.step(entry: BandEntry, tail: List<&2, BandEntry>, rest: List<&2, Peq>,
272 first: Nat, keep: Bool, next: Unit -> Band) -> Band:
273 match keep:
274 case True{}:
275 Band{first, entry <> tail, rest}
276 case False{}:
277 next(Unit{})
278
279def trim_first(entries: List<&2, BandEntry>, +rest: List<&2, Peq>, +first: Nat,
280 +maximum: Nat, +length: Nat, +target_length: Nat,
281 +column: Nat) -> Band:
282 match entries:
283 case Nil{}:
284 Band{first, Nil{}, rest}
285 case BandEntry{profile, Block{pv, mv, +score}} <> +tail:
286 trim_first.step(BandEntry{profile, Block{pv, mv, score}}, tail, rest,
287 first,
288 band_keep_first(U32.to_nat(score), Nat.mul(1n+first, 32n),
289 maximum, length, target_length, column),
290 _ => trim_first(tail, rest, 1n+first, maximum, length, target_length, column))
291
292def reduce_band.step(parts: List<&2, BandEntry> & List<&2, Peq>, first: Nat,
293 +maximum: Nat, length: Nat, target_length: Nat,
294 column: Nat, work: Nat) -> BandStep:
295 match parts:
296 case (entries, rest):
297 BandStep{trim_first(entries, rest, first, maximum, length, target_length, column), maximum, work}
298
299def reduce_band(entries: List<&2, BandEntry>, rest: List<&2, Peq>, +first: Nat,
300 +maximum: Nat, +length: Nat, +target_length: Nat,
301 +column: Nat, work: Nat) -> BandStep:
302 reduce_band.step(trim_last(entries, rest, Nat.mul(1n+first, 32n),
303 maximum, length, target_length, column), first, maximum, length, target_length, column, work)
304
305def expand_band.step(result: Block & Carry, entries: List<&2, BandEntry>, profile: Peq,
306 rest: List<&2, Peq>, first: Nat, maximum: Nat, length: Nat,
307 target_length: Nat, column: Nat, work: Nat) -> BandStep:
308 match result:
309 case (block, _):
310 reduce_band(List.append(&2, BandEntry, entries, [BandEntry{profile, block}]),
311 rest, first, maximum, length, target_length, column, 1n+work)
312
313def expand_band(entries: List<&2, BandEntry>, rest: List<&2, Peq>, +score: U32,
314 carry: Carry, +symbol: U32, first: Nat, +maximum: Nat,
315 length: Nat, target_length: Nat, column: Nat,
316 work: Nat, expand: Bool) -> BandStep:
317 match rest carry expand:
318 case +profile <> tail Carry{+positive, +negative} True{}:
319 previous = U32.add(U32.sub(U32.add(score, negative), positive), 32)
320 expand_band.step(calculate_block(Block{4294967295, 0, previous},
321 peq_get(profile, symbol, U32.is_ne(U32.and(symbol, 1), 0)), Carry{positive, negative}),
322 entries, profile, tail, first, maximum, length, target_length, column, work)
323 case _ _ _:
324 reduce_band(entries, rest, first, maximum, length, target_length, column, work)
325
326def band_step.updated(updated: UpdatedBand, rest: List<&2, Peq>, +first: Nat,
327 symbol: U32, +length: Nat, +target_length: Nat,
328 +column: Nat, maximum: Nat, padding: Nat) -> BandStep:
329 match updated:
330 case UpdatedBand{entries, +score, carry, +bottom}:
331 # Completing the suffix by edits bounds the optimum. Padding can lower
332 # the stored bottom score, so add it back before tightening k.
333 +k = Nat.min(maximum, Nat.add(U32.to_nat(score),
334 Nat.add(Nat.max(Nat.sub(target_length, column), Nat.sub(length, bottom)),
335 Bool.pick(Nat, Nat.is_ge(bottom, length), padding, 0n))))
336 work = Nat.sub(Nat.div(bottom, 32n), first)
337 expand_band(entries, rest, score, carry, symbol, first, k, length, target_length,
338 column, work, Nat.is_le(Nat.add(Nat.add(bottom, U32.to_nat(score)), target_length),
339 Nat.add(Nat.add(Nat.add(k, length), column), 62n)))
340
341def band_step(band: Band, +symbol: U32, length: Nat, target_length: Nat,
342 column: Nat, maximum: Nat, padding: Nat) -> BandStep:
343 match band:
344 case Band{+first, entries, rest}:
345 band_step.updated(advance_band(entries, symbol, Carry{1, 0}, Nat.mul(1n+first, 32n)),
346 rest, first, symbol, length, target_length, column, maximum, padding)
347
348def band_last_score(entries: List<&2, BandEntry>, +padding: Nat) -> Nat:
349 match entries:
350 case Nil{}:
351 0n
352 case BandEntry{_, Block{pv, mv, score}} <> Nil{}:
353 U32.to_nat(get_block_cell_value(padding, pv, mv, score))
354 case _ <> tail:
355 band_last_score(tail, padding)
356
357def band_finish(band: Band, maximum: Nat, padding: Nat,
358 blocks: Nat, columns: Nat) -> Search:
359 match band:
360 case Band{_, Nil{}, _}:
361 Search{E.NotFound{}, blocks, columns, 1n}
362 case Band{_, entries, Nil{}}:
363 +score = band_last_score(entries, padding)
364 Search{Bool.pick(E.EditDistance, Nat.is_le(score, maximum), E.Found{score}, E.NotFound{}),
365 blocks, columns, 1n}
366 case _:
367 Search{E.NotFound{}, blocks, columns, 1n}
368
369def myers_calc_edit_distance_nw.step(step: BandStep, next: Band -> Nat -> Nat -> Search) -> Search:
370 match step:
371 case BandStep{band, maximum, work}:
372 next(band, maximum, work)
373
374def myers_calc_edit_distance_nw(target: String, band: Band, +length: Nat, +target_length: Nat,
375 maximum: Nat, +padding: Nat, +column: Nat, blocks: Nat) -> Search:
376 match target band:
377 case _ Band{_, Nil{}, _}:
378 Search{E.NotFound{}, blocks, column, 1n}
379 case SNil{} _:
380 band_finish(band, maximum, padding, blocks, column)
381 case SCon{h, t} _:
382 myers_calc_edit_distance_nw.step(band_step(band, Char.to_u32(h), length, target_length,
383 1n+column, maximum, padding), next => k => work =>
384 myers_calc_edit_distance_nw(t, next, length, target_length, k, padding, 1n+column, Nat.add(blocks, work)))
385
386def band_attempt(target: String, profiles: List<&2, Peq>, +length: Nat,
387 +target_length: Nat, +maximum: Nat, padding: Nat) -> Search:
388 myers_calc_edit_distance_nw(target, init_band(profiles, length, target_length, maximum),
389 length, target_length, maximum, padding, 0n, 0n)
390
391def adaptive_search.step(distance: E.EditDistance, blocks: Nat, columns: Nat,
392 attempts: Nat, retry: Unit -> Search) -> Search:
393 match distance:
394 case E.Found{value}:
395 Search{E.Found{value}, blocks, columns, attempts}
396 case E.NotFound{}:
397 retry(Unit{})
398
399def adaptive_search.result(search: Search, blocks: Nat, columns: Nat, +attempts: Nat,
400 retry: Nat -> Nat -> Search) -> Search:
401 match search:
402 case Search{distance, work, consumed, _}:
403 +next_blocks = Nat.add(blocks, work)
404 +next_columns = Nat.add(columns, consumed)
405 adaptive_search.step(distance, next_blocks, next_columns, 1n+attempts,
406 _ => retry(next_blocks, next_columns))
407
408def adaptive_search(fuel: Nat, +target: String, +profiles: List<&2, Peq>,
409 +length: Nat, +target_length: Nat, +maximum: Nat,
410 +padding: Nat, +blocks: Nat, +columns: Nat,
411 +attempts: Nat) -> Search:
412 match fuel:
413 case 0n:
414 Search{E.NotFound{}, blocks, columns, attempts}
415 case 1n+p:
416 adaptive_search.result(band_attempt(target, profiles, length, target_length, maximum, padding),
417 blocks, columns, attempts, bs => cs =>
418 adaptive_search(p, target, profiles, length, target_length,
419 Nat.min(target_length, Nat.mul(maximum, 2n)), padding, bs, cs, 1n+attempts))
420
421def search_profiles(target: String, profiles: List<&2, Peq>, +length: Nat,
422 +target_length: Nat, threshold: E.Threshold) -> Search:
423 match threshold:
424 case E.MaxDistance{maximum}:
425 band_attempt(target, profiles, length, target_length, Nat.min(maximum, target_length),
426 Nat.sub(Nat.mul(ceil_div(length, 32n), 32n), length))
427 case E.Unlimited{}:
428 adaptive_search(1n+target_length, target, profiles, length, target_length,
429 Nat.min(target_length, Nat.max(32n, Nat.sub(target_length, length))),
430 Nat.sub(Nat.mul(ceil_div(length, 32n), 32n), length), 0n, 0n, 0n)
431
432def search_query(+query: String, target: String, +length: Nat, target_length: Nat,
433 threshold: E.Threshold, pairs: List<&2, E.EqualityPair>,
434 rejected: Bool) -> Search:
435 match rejected:
436 case True{}:
437 Search{E.NotFound{}, 0n, 0n, 0n}
438 case False{}:
439 search_profiles(target, build_peq(ceil_div(length, 32n), query, pairs), length, target_length, threshold)
440
441def search_shorter(query: String, target: String, +query_length: Nat,
442 +target_length: Nat, +threshold: E.Threshold,
443 query_is_shorter: Bool, pairs: List<&2, E.EqualityPair>) -> Search:
444 match query_is_shorter:
445 case True{}:
446 search_query(query, target, query_length, target_length, threshold, pairs,
447 Bool.not(E.within_threshold(threshold, Nat.sub(target_length, query_length))))
448 case False{}:
449 search_query(target, query, target_length, query_length, threshold, pairs,
450 Bool.not(E.within_threshold(threshold, Nat.sub(query_length, target_length))))
451
452# Directional modes never swap query and target. The initial vertical scores
453# are i in both modes, while the top horizontal difference is +1 for SHW and
454# 0 for HW. NW's length/remaining-target band bounds do not apply here.
455def semiglobal_carry(mode: E.Mode) -> Carry:
456 match mode:
457 case E.HW{}:
458 Carry{0, 0}
459 case _:
460 Carry{1, 0}
461
462def init_semiglobal_band(profiles: List<&2, Peq>, length: Nat, +maximum: Nat) -> Band:
463 init_band.step(init_band_entries(
464 Nat.min(ceil_div(length, 32n), ceil_div(1n+maximum, 32n)), profiles, 32))
465
466# If the bottom score is >= k+32 every cell in that block is > k. Only
467# trailing blocks are removed. HW always retains block zero: its free top
468# boundary can seed a new alignment after an arbitrarily long target prefix.
469def trim_semiglobal(entries: List<&2, BandEntry>, rest: List<&2, Peq>,
470 +mode: E.Mode, +index: Nat, +maximum: Nat) -> List<&2, BandEntry> & List<&2, Peq>:
471 match entries mode:
472 case Nil{} _:
473 (Nil{}, rest)
474 case BandEntry{profile, Block{pv, mv, +score}} <> tail E.HW{}:
475 trim_last.join(BandEntry{profile, Block{pv, mv, score}},
476 trim_semiglobal(tail, rest, mode, 1n+index, maximum),
477 Bool.or(Nat.is_eq(index, 0n), Nat.is_lt(U32.to_nat(score), Nat.add(maximum, 32n))))
478 case BandEntry{profile, Block{pv, mv, +score}} <> tail _:
479 trim_last.join(BandEntry{profile, Block{pv, mv, score}},
480 trim_semiglobal(tail, rest, mode, 1n+index, maximum),
481 Nat.is_lt(U32.to_nat(score), Nat.add(maximum, 32n)))
482
483def reduce_semiglobal(parts: List<&2, BandEntry> & List<&2, Peq>) -> Band:
484 match parts:
485 case (entries, rest):
486 Band{0n, entries, rest}
487
488def expand_semiglobal.result(result: Block & Carry, entries: List<&2, BandEntry>,
489 profile: Peq, rest: List<&2, Peq>, mode: E.Mode,
490 maximum: Nat) -> Band:
491 match result:
492 case (block, _):
493 reduce_semiglobal(trim_semiglobal(
494 List.append(&2, BandEntry, entries, [BandEntry{profile, block}]), rest, mode, 0n, maximum))
495
496def expand_semiglobal.choose(entries: List<&2, BandEntry>, profile: Peq,
497 rest: List<&2, Peq>, previous: U32, eq: U32,
498 carry: Carry, mode: E.Mode, maximum: Nat, expand: Bool) -> Band:
499 match expand:
500 case True{}:
501 expand_semiglobal.result(calculate_block(Block{4294967295, 0, U32.add(previous, 32)}, eq, carry),
502 entries, profile, rest, mode, maximum)
503 case False{}:
504 reduce_semiglobal(trim_semiglobal(entries, profile <> rest, mode, 0n, maximum))
505
506def expand_semiglobal(updated: UpdatedBand, rest: List<&2, Peq>, +symbol: U32,
507 mode: E.Mode, +maximum: Nat) -> Band:
508 match updated rest:
509 case UpdatedBand{entries, +score, Carry{+positive, +negative}, _} +profile <> tail:
510 +previous = U32.sub(U32.add(score, negative), positive)
511 +eq = peq_get(profile, symbol, U32.is_ne(U32.and(symbol, 1), 0))
512 expand_semiglobal.choose(entries, profile, tail, previous, eq, Carry{positive, negative}, mode, maximum,
513 Bool.and(Nat.is_le(U32.to_nat(previous), maximum),
514 Bool.or(U32.is_ne(U32.and(eq, 1), 0), U32.is_ne(negative, 0))))
515 case UpdatedBand{entries, _, _, _} Nil{}:
516 reduce_semiglobal(trim_semiglobal(entries, Nil{}, mode, 0n, maximum))
517
518def semiglobal_step(band: Band, +symbol: U32, +mode: E.Mode, maximum: Nat) -> Band:
519 match band:
520 case Band{_, entries, rest}:
521 expand_semiglobal(advance_band(entries, symbol, semiglobal_carry(mode), 32n),
522 rest, symbol, mode, maximum)
523
524def semiglobal_endpoint(entries: List<&2, BandEntry>, rest: List<&2, Peq>, padding: Nat,
525 column: Nat, +best: A.Best) -> A.Best:
526 match entries rest:
527 case Nil{} _:
528 best
529 case _ Nil{}:
530 A.Best{old, _} = best
531 +score = U32.from_nat(band_last_score(entries, padding))
532 A.update_best(best, score, A.Location{E.Position{column}, 0n}, U32.cmp(score, old))
533 case _ _:
534 best
535
536def myers_calc_edit_distance_semiglobal.step(band: Band, padding: Nat, column: Nat,
537 best: A.Best, next: Band -> A.Best -> A.Best) -> A.Best:
538 match band:
539 case Band{first, +entries, +rest}:
540 next(Band{first, entries, rest}, semiglobal_endpoint(entries, rest, padding, column, best))
541
542# Decode the real last query cell at every column, without Edlib's delayed
543# padded endpoints. Tighten k after a better score, but continue for every tie.
544def myers_calc_edit_distance_semiglobal(target: String, band: Band, +mode: E.Mode,
545 +padding: Nat, +column: Nat, +maximum: Nat,
546 +best: A.Best) -> A.Best:
547 match target band:
548 case _ Band{_, Nil{}, _}:
549 A.ordered(best)
550 case SNil{} _:
551 A.ordered(best)
552 case SCon{h, t} _:
553 A.Best{score, _} = best
554 +k = Nat.min(maximum, U32.to_nat(score))
555 myers_calc_edit_distance_semiglobal.step(semiglobal_step(band, Char.to_u32(h), mode, k),
556 padding, column, best, next => found =>
557 myers_calc_edit_distance_semiglobal(t, next, mode, padding, 1n+column, k, found))
558
559def semiglobal_initial_best(+length: Nat, include_before: Bool) -> A.Best:
560 match include_before:
561 case True{}:
562 A.Best{U32.from_nat(length), [A.Location{E.BeforeStart{}, 0n}]}
563 case False{}:
564 A.Best{4294967295, Nil{}}
565
566def semiglobal_maximum(threshold: E.Threshold, +length: Nat) -> Nat:
567 match threshold:
568 case E.Unlimited{}:
569 length
570 case E.MaxDistance{maximum}:
571 Nat.min(maximum, length)
572
573def semiglobal_best(target: String, profiles: List<&2, Peq>, +length: Nat,
574 +maximum: Nat, +mode: E.Mode) -> A.Best:
575 +count = ceil_div(length, 32n)
576 # Preserve pinned Edlib's 64-bit column-zero convention independently of
577 # this implementation's 32-bit padding (especially lengths 32 and 96).
578 best = semiglobal_initial_best(length, Nat.is_ne(Nat.mod(length, 64n), 0n))
579 myers_calc_edit_distance_semiglobal(target,
580 init_semiglobal_band(profiles, length, maximum), mode,
581 Nat.sub(Nat.mul(count, 32n), length), 0n, maximum, best)
582
583# Edlib selects the LAST optimal reverse-SHW end, hence the earliest start
584# for each forward end. Preserve the special forward end -1 with start 0.
585def hw_start(+end: Nat, reverse_location: A.Location) -> A.Location:
586 match reverse_location:
587 case A.Location{E.BeforeStart{}, _}:
588 A.Location{E.Position{end}, 1n+end}
589 case A.Location{E.Position{reverse_end}, _}:
590 A.Location{E.Position{end}, Nat.sub(end, reverse_end)}
591
592def hw_start_best(+end: Nat, best: A.Best) -> A.Location:
593 match best:
594 case A.Best{_, locations}:
595 hw_start(end, Maybe.default(&2, A.Location,
596 List.last(&2, A.Location, locations), A.Location{E.BeforeStart{}, 0n}))
597
598def hw_location(location: A.Location, reverse_target: String, target_length: Nat,
599 profiles: List<&2, Peq>, query_length: Nat, maximum: Nat) -> A.Location:
600 match location:
601 case A.Location{E.BeforeStart{}, _}:
602 A.Location{E.BeforeStart{}, 0n}
603 case A.Location{E.Position{+end}, _}:
604 hw_start_best(end, semiglobal_best(
605 String.drop(reverse_target, Nat.sub(target_length, 1n+end)),
606 profiles, query_length, maximum, E.SHW{}))
607
608# Independent starts share one reversed query profile. Split their searches
609# into balanced parallel branches, then concatenate in the original end order.
610def hw_locations(fuel: Nat, locations: List<&2, A.Location>,
611 +reverse_target: String, +target_length: Nat,
612 +profiles: List<&2, Peq>, +query_length: Nat,
613 +maximum: Nat) -> List<&2, A.Location>:
614 match fuel locations:
615 case 0n _:
616 locations
617 case 1n+p Nil{}:
618 Nil{}
619 case 1n+p location <> Nil{}:
620 [hw_location(location, reverse_target, target_length, profiles, query_length, maximum)]
621 case Succ{+p} +values:
622 +half = Nat.div(List.length(&2, A.Location, values), 2n)
623 left right =
624 hw_locations(p, List.take(&2, A.Location, values, half),
625 reverse_target, target_length, profiles, query_length, maximum)
626 hw_locations(p, List.drop(&2, A.Location, values, half),
627 reverse_target, target_length, profiles, query_length, maximum)
628 List.append(&2, A.Location, left, right)
629
630def configured_found(task: E.Task, mode: E.Mode, +query: String, +target: String,
631 +pairs: List<&2, E.EqualityPair>, format: E.CigarFormat,
632 best: A.Best) -> E.AlignResult:
633 match task mode best:
634 case E.Distance{} _ _:
635 A.result(query, target, task, pairs, format, best, True{})
636 case _ E.HW{} A.Best{+score, +locations}:
637 +qn +tn = String.length(query) String.length(target)
638 rq rt = String.reverse(query) String.reverse(target)
639 resolved = hw_locations(List.length(&2, A.Location, locations), locations,
640 rt, tn, build_peq(ceil_div(qn, 32n), rq, pairs), qn, U32.to_nat(score))
641 A.result(query, target, task, pairs, format, A.Best{score, resolved}, True{})
642 case _ _ _:
643 # NW/SHW always start at zero; only PATH needs reconstruction.
644 A.result(query, target, task, pairs, format, best, True{})
645
646def configured_result(query: String, target: String, mode: E.Mode, task: E.Task,
647 pairs: List<&2, E.EqualityPair>, format: E.CigarFormat,
648 best: A.Best, accepted: Bool) -> E.AlignResult:
649 match accepted:
650 case False{}:
651 E.not_found()
652 case True{}:
653 configured_found(task, mode, query, target, pairs, format, best)
654
655def configured_finish(query: String, target: String, mode: E.Mode, task: E.Task,
656 threshold: E.Threshold, pairs: List<&2, E.EqualityPair>,
657 format: E.CigarFormat, +best: A.Best) -> E.AlignResult:
658 A.Best{score, _} = best
659 configured_result(query, target, mode, task, pairs, format, best,
660 E.within_threshold(threshold, U32.to_nat(score)))
661
662def semiglobal_align(+query: String, +target: String, +mode: E.Mode, task: E.Task,
663 +threshold: E.Threshold, +pairs: List<&2, E.EqualityPair>,
664 format: E.CigarFormat) -> E.AlignResult:
665 +length = String.length(query)
666 configured_finish(query, target, mode, task, threshold, pairs, format,
667 semiglobal_best(target, build_peq(ceil_div(length, 32n), query, pairs),
668 length, semiglobal_maximum(threshold, length), mode))
669
670def search_distance(search: Search) -> E.EditDistance:
671 match search:
672 case Search{distance, _, _, _}:
673 distance
674
675def found_value(distance: E.EditDistance) -> Nat:
676 match distance:
677 case E.Found{value}:
678 value
679 case E.NotFound{}:
680 0n
681
682def align_shorter(query: String, target: String, query_length: Nat,
683 target_length: Nat, query_is_shorter: Bool) -> Nat:
684 found_value(search_distance(search_shorter(query, target, query_length,
685 target_length, E.Unlimited{}, query_is_shorter, Nil{})))
686
687def distance(query: String, target: String) -> Nat:
688 match query target:
689 case SNil{} _:
690 String.length(target)
691 case SCon{h, t} SNil{}:
692 1n+String.length(t)
693 case SCon{qh, qt} SCon{th, tt}:
694 +q = {SCon{qh, qt} : String}
695 +t = {SCon{th, tt} : String}
696 +qn +tn = String.length(q) String.length(t)
697 align_shorter(q, t, qn, tn, Nat.is_le(qn, tn))
698
699def nw_result(query: String, target: String, task: E.Task,
700 pairs: List<&2, E.EqualityPair>, format: E.CigarFormat,
701 end: E.Coordinate, distance: E.EditDistance) -> E.AlignResult:
702 match distance:
703 case E.NotFound{}:
704 E.not_found()
705 case E.Found{score}:
706 configured_result(query, target, E.NW{}, task, pairs, format,
707 A.Best{U32.from_nat(score), [A.Location{end, 0n}]}, True{})
708
709def edlib_align(query: String, target: String, config: E.AlignConfig) -> E.AlignResult:
710 match query target config:
711 case SNil{} _ _:
712 A.align("", target, config)
713 case SCon{h, t} SNil{} _:
714 A.align(SCon{h, t}, "", config)
715 case +q +t E.AlignConfig{E.NW{}, task, threshold, +pairs, format}:
716 +qn +tn = String.length(q) String.length(t)
717 nw_result(q, t, task, pairs, format, E.end_coordinate(tn),
718 search_distance(search_shorter(q, t, qn, tn, threshold, Nat.is_le(qn, tn), pairs)))
719 case +q +t E.AlignConfig{mode, task, threshold, pairs, format}:
720 semiglobal_align(q, t, mode, task, threshold, pairs, format)
721
722def response(query: String, target: String) -> String:
723 Nat.show(distance(query, target))
724
725def run(args: List<String>) -> IO(Unit):
726 match args:
727 case _ <> query <> target <> Nil{}:
728 CLI.run(~distance, ["bend-edit-distance", query, target])
729 case _:
730 AC.run(~edlib_align, args)
731
732def main() -> IO(Unit):
733 do IO<Unit>:
734 args : List<String> <- IO.args()
735 run(args)