1 import Base 2 import ./Sequence.bend as S 3 import ../EdlibTypes.bend as E 4 import ./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 11 type Peq is Data: 12 Mask{bits: U32} 13 Branch{low: Peq, high: Peq} 14 15 type 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. 20 type Carry is Data: 21 Carry{positive: U32, negative: U32} 22 23 def split_query.prepend(h: Char, parts: String & String) -> String & String: 24 match parts: 25 case (prefix, suffix): 26 (SCon{h, prefix}, suffix) 27 28 def 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. 39 def 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. 51 def 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 60 def 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 71 def 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 77 def 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 84 def ceil_div(+x: Nat, +y: Nat) -> Nat: 85 Nat.div(Nat.add(x, Nat.sub(y, 1n)), y) 86 87 def init_blocks(n: Nat, +bottom: U32) -> List: 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. 96 def 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. 115 def advance_column.step(result: Block & Carry, 116 tail: Carry -> List) -> List: 117 match result: 118 case (next, hout): 119 next <> tail(hout) 120 121 def advance_column(blocks: List, profiles: List<&2, Peq>, +symbol: U32, 122 hin: Carry) -> List: 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. 136 def 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 144 def last_block_score(blocks: List, +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. 156 type BandEntry is Data: 157 BandEntry{profile: Peq, block: Block} 158 159 type Band is Type: 160 Band{first: Nat, entries: List<&2, BandEntry>, rest: List<&2, Peq>} 161 162 type UpdatedBand is Type: 163 UpdatedBand{entries: List<&2, BandEntry>, score: U32, carry: Carry, bottom: Nat} 164 165 type 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. 170 type Search is Data: 171 Search{distance: E.EditDistance, blocks: Nat, columns: Nat, attempts: Nat} 172 173 def 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 179 def 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 189 def init_band.step(parts: List<&2, BandEntry> & List<&2, Peq>) -> Band: 190 match parts: 191 case (entries, rest): 192 Band{0n, entries, rest} 193 194 def 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 201 def 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 206 def 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 211 def 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 217 def 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. 233 def 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 239 def 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 245 def 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 253 def 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 259 def 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 270 def 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 278 def 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 291 def 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 298 def 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 304 def 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 312 def 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 325 def 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 340 def 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 347 def 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 356 def 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 368 def 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 373 def 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 385 def 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 390 def 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 398 def 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 407 def 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 420 def 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 431 def 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 440 def 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. 454 def semiglobal_carry(mode: E.Mode) -> Carry: 455 match mode: 456 case E.HW{}: 457 Carry{0, 0} 458 case _: 459 Carry{1, 0} 460 461 def 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. 468 def 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 482 def reduce_semiglobal(parts: List<&2, BandEntry> & List<&2, Peq>) -> Band: 483 match parts: 484 case (entries, rest): 485 Band{0n, entries, rest} 486 487 def 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 495 def 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 505 def 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 517 def 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 523 def 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 535 def 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. 543 def 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 558 def 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 565 def 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 572 def 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. 584 def 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 591 def 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 597 def 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. 609 def 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 629 def 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 645 def 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 654 def 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 661 def 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 669 def search_distance(search: Search) -> E.EditDistance: 670 match search: 671 case Search{distance, _, _, _}: 672 distance 673 674 def found_value(distance: E.EditDistance) -> Nat: 675 match distance: 676 case E.Found{value}: 677 value 678 case E.NotFound{}: 679 0n 680 681 def 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 686 def 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 698 def 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 708 def 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 721 def response(query: String, target: String) -> String: 722 Nat.show(distance(query, target))