src/Myers.bend

bend-edit-distance@0.2.2.0 · 722 lines · 83 definitions

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