nx_sequence_edit.nx source
↩ module page · 119 lines · 4342 B
1// nx_sequence_edit.nx -- Levenshtein edit distance for byte sequences.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/wagner-fischer-1974
5//
6// G0.2 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Classic
7// Wagner-Fischer dynamic-programming algorithm:
8//
9// d[i][j] = min(d[i-1][j] + 1, // delete a[i-1]
10// d[i][j-1] + 1, // insert b[j-1]
11// d[i-1][j-1] + (a[i-1] != b[j-1])) // substitute / match
12// d[0][0] = 0; d[i][0] = i; d[0][j] = j
13//
14// Works on raw byte sequences -- equality is byte-equal. Caller may
15// pass ASCII bases, packed 2-bit codes one per byte, amino acids, or
16// any other byte alphabet; the algorithm is alphabet-agnostic.
17//
18// Why classic O(mn) DP for the reference impl:
19// - Bit-equal reproducible across hosts (no SIMD-order nondeterminism)
20// - Trivially auditable against the textbook recurrence
21// - Reference behaviour to validate the Myers' bitvector fast path
22// (nx_sequence_edit_bv.nx, G0.5) and Edlib parity bench (G14)
23//
24// Memory:
25// Internal sys_mmap of (n+1) * (m+1) * 8 bytes for the DP table.
26// For n=m=200 -> 320 KB; for n=m=1000 -> 8 MB. Long-read alignment
27// (n,m ~ 10000) will require the rolling two-row optimisation
28// (nx_sequence_edit_lin.nx, G0.5).
29//
30// API:
31// edit_distance(a, n, b, m) -> i64
32// Returns the Levenshtein distance, or -1 if either length is
33// negative. n + m == 0 returns 0.
34//
35// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
36// intended_use: "Edit distance for read-to-reference compare,
37// guide-RNA off-target enumeration, primer
38// validation, sequence-similarity scoring"
39// sil_target: SIL2
40// asil_target: QM
41// dal_target: DAL C
42// iec_62304_class: B
43// evidence: [no_floating_point, deterministic,
44// bit_equal_reproducible,
45// wagner_fischer_textbook_recurrence,
46// published_kat_vectors_verified,
47// license_tier_INDEPENDENT_REDERIVE]
48// hazard_register: [bug-tape-edit-distance-off-by-one-boundary,
49// bug-tape-edit-distance-mn-overflow-table-size,
50// bug-tape-edit-distance-min3-tie-break-drift]
51// residual_risk: "O(mn) memory blows up for n,m > ~1500 on a
52// 4096-byte page; callers using long-read data
53// must wait for the rolling-row variant in
54// G0.5 or stage allocation themselves."
55// verdict: NOT_YET_EVALUATED
56
57import "nx_syscalls.nx"
58
59// 3-way min on i64. Pulled out for clarity.
60func nx_min3(a: i64, b: i64, c: i64) -> i64 {
61 var m: i64 = a
62 if b < m { m = b }
63 if c < m { m = c }
64 return m
65}
66
67// Wagner-Fischer Levenshtein distance.
68//
69// Flattens d[i][j] as table[i * (m+1) + j] for cache-friendly access
70// along the inner j loop.
71func edit_distance(a: *u8, n: i64, b: *u8, m: i64) -> i64 {
72 if n < 0 { return -1 }
73 if m < 0 { return -1 }
74 if n == 0 { return m }
75 if m == 0 { return n }
76
77 let cols: i64 = m + 1
78 let rows: i64 = n + 1
79 let bytes: i64 = rows * cols * 8
80
81 let table: *i64 = sys_mmap(bytes) as *i64
82
83 // Boundary row + column.
84 var j0: i64 = 0
85 while j0 <= m {
86 table[j0] = j0
87 j0 = j0 + 1
88 }
89 var i0: i64 = 1
90 while i0 <= n {
91 table[i0 * cols] = i0
92 i0 = i0 + 1
93 }
94
95 // Inner DP. Iterate i = 1..n, j = 1..m.
96 var i: i64 = 1
97 while i <= n {
98 let a_byte: i64 = a[i - 1] & 0xff
99 var j: i64 = 1
100 while j <= m {
101 let b_byte: i64 = b[j - 1] & 0xff
102 let row_off: i64 = i * cols
103 let prev_row_off: i64 = (i - 1) * cols
104
105 let del_cost: i64 = table[prev_row_off + j] + 1
106 let ins_cost: i64 = table[row_off + j - 1] + 1
107
108 var sub_penalty: i64 = 1
109 if a_byte == b_byte { sub_penalty = 0 }
110 let sub_cost: i64 = table[prev_row_off + j - 1] + sub_penalty
111
112 table[row_off + j] = nx_min3(del_cost, ins_cost, sub_cost)
113 j = j + 1
114 }
115 i = i + 1
116 }
117
118 return table[n * cols + m]
119}