nx_sequence_fm.nx source
↩ module page · 239 lines · 8519 B
1// nx_sequence_fm.nx -- FM-index construction + backward-search count.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/ferragina-manzini-2000
5//
6// G0.3 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Built bits-up
7// against the Ferragina-Manzini construction:
8//
9// 1. Suffix array SA of T (where T ends in a unique sentinel '$').
10// G0.3 uses a naive O(n^2 log n) insertion sort over suffix
11// indices, lexicographically comparing the suffixes. Correct
12// reference behaviour; SA-IS fast path is G0.5.
13// 2. Burrows-Wheeler Transform: BWT[i] = T[(SA[i] - 1) mod n].
14// 3. C array: C[c] = number of bytes in BWT lexicographically
15// smaller than c (0..255 entries; cumulative histogram).
16// 4. Backward search (Ferragina-Manzini fm_count algorithm):
17// c <- P[plen-1]
18// lo <- C[c]
19// hi <- C[c+1] - 1
20// for i = plen-2 downto 0:
21// c <- P[i]
22// lo <- C[c] + Occ(c, lo - 1) + 1
23// hi <- C[c] + Occ(c, hi)
24// if lo > hi: return 0
25// return hi - lo + 1
26// where Occ(c, i) = count of byte c in BWT[0..i] inclusive.
27// G0.3 recomputes Occ by walking BWT each query (O(plen * n));
28// the sampled-Occ wavelet-tree fast path is G0.5.
29//
30// Sentinel convention:
31// Caller passes T with a trailing $ byte (0x24) which must not
32// appear elsewhere in T. All construction + query helpers assume
33// the sentinel is in place and is the lexicographically smallest
34// byte present. fm_check_sentinel returns -1 if the convention is
35// violated -- callers should run it once before fm_build_sa.
36//
37// API:
38// fm_check_sentinel(t, n) -> i64 0 ok / -1 bad
39// fm_build_sa(t, n, sa_out) -> i64
40// fm_build_bwt(t, n, sa, bwt_out) -> i64
41// fm_build_c(bwt, n, c_out_256) -> i64
42// fm_occ(bwt, n, c, i) -> i64 count of c in bwt[0..i]
43// fm_count(bwt, n, c_arr, p, plen) -> i64 # occurrences of p
44//
45// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
46// intended_use: "Substring search index for reference genome
47// + CRISPR off-target enumeration + read seeding
48// (composes with future nx_align primitive)"
49// sil_target: SIL2
50// asil_target: QM
51// dal_target: DAL C
52// iec_62304_class: B
53// evidence: [no_floating_point, deterministic,
54// bit_equal_reproducible,
55// ferragina_manzini_textbook_construction,
56// banana_canonical_kat_verified,
57// sentinel_invariant_enforced,
58// license_tier_INDEPENDENT_REDERIVE]
59// hazard_register: [bug-tape-suffix-array-tie-on-shared-prefix,
60// bug-tape-bwt-wrap-mod-off-by-one,
61// bug-tape-c-array-uses-bwt-not-text,
62// bug-tape-occ-inclusive-vs-exclusive-bound,
63// bug-tape-backward-search-empty-range-not-caught]
64// residual_risk: "G0.3 naive SA build is O(n^2 log n) -- fine
65// for fixtures up to ~1000 bases; the human
66// reference (3 Gbp) requires SA-IS in G0.5.
67// Naive Occ is O(plen * n) per query -- the
68// sampled-Occ wavelet-tree fast path lives
69// in G0.5 as well."
70// verdict: NOT_YET_EVALUATED
71
72import "nx_syscalls.nx"
73import "nx_const.nx"
74
75// FM-index sentinel byte is NX_FM_NX_FM_SENTINEL_BYTE from nx_const.
76// NishiLang's `const` requires integer literals so we reference
77// NX_FM_NX_FM_SENTINEL_BYTE directly throughout this module rather than
78// re-aliasing. See docs/NISHI_CONSTANT_REGISTRY.md.
79
80// Verify the caller's text has the sentinel only at the end.
81// Returns 0 on success, -1 if sentinel is missing from t[n-1] or if
82// it appears anywhere in t[0..n-2].
83func fm_check_sentinel(t: *u8, n: i64) -> i64 {
84 if n <= 0 { return -1 }
85 if (t[n - 1] & 0xff) != NX_FM_SENTINEL_BYTE { return -1 }
86 var i: i64 = 0
87 while i < n - 1 {
88 if (t[i] & 0xff) == NX_FM_SENTINEL_BYTE { return -1 }
89 i = i + 1
90 }
91 return 0
92}
93
94// Lexicographic compare of two suffixes of t starting at indices a, b.
95// Returns < 0 if suffix(a) < suffix(b); > 0 if greater; 0 if equal
96// (only possible when a == b given the unique-sentinel invariant).
97func fm_suffix_compare(t: *u8, n: i64, a: i64, b: i64) -> i64 {
98 var i: i64 = a
99 var j: i64 = b
100 while i < n {
101 if j >= n { return 1 } // a-suffix longer at this point
102 let ca: i64 = t[i] & 0xff
103 let cb: i64 = t[j] & 0xff
104 if ca < cb { return -1 }
105 if ca > cb { return 1 }
106 i = i + 1
107 j = j + 1
108 }
109 if j < n { return -1 } // b-suffix is longer
110 return 0
111}
112
113// Build the suffix array by insertion sort over indices 0..n-1.
114// O(n^2 log n) worst case but n^2 in expectation for low-shared-prefix
115// alphabets like DNA. Correct reference impl; SA-IS lives in G0.5.
116// sa_out must have capacity >= n * 8 bytes.
117func fm_build_sa(t: *u8, n: i64, sa_out: *i64) -> i64 {
118 if n <= 0 { return -1 }
119
120 // Seed with identity permutation 0..n-1.
121 var k: i64 = 0
122 while k < n {
123 sa_out[k] = k
124 k = k + 1
125 }
126
127 // Insertion sort: for i = 1..n-1, slide sa_out[i] left into place.
128 var i: i64 = 1
129 while i < n {
130 let key: i64 = sa_out[i]
131 var j: i64 = i - 1
132 var keep_going: i64 = 1
133 while keep_going == 1 {
134 if j < 0 {
135 keep_going = 0
136 } else {
137 let cmp: i64 = fm_suffix_compare(t, n, sa_out[j], key)
138 if cmp > 0 {
139 sa_out[j + 1] = sa_out[j]
140 j = j - 1
141 } else {
142 keep_going = 0
143 }
144 }
145 }
146 sa_out[j + 1] = key
147 i = i + 1
148 }
149 return 0
150}
151
152// Compute BWT from text + suffix array.
153// BWT[i] = T[(SA[i] - 1) mod n]. bwt_out must have capacity >= n bytes.
154func fm_build_bwt(t: *u8, n: i64, sa: *i64, bwt_out: *u8) -> i64 {
155 if n <= 0 { return -1 }
156 var i: i64 = 0
157 while i < n {
158 let s: i64 = sa[i]
159 var src_idx: i64 = s - 1
160 if src_idx < 0 { src_idx = n - 1 }
161 bwt_out[i] = t[src_idx] & 0xff
162 i = i + 1
163 }
164 return 0
165}
166
167// Build the C array (cumulative byte histogram of BWT).
168// C[c] = number of bytes in BWT that are strictly less than c.
169// c_out must have capacity >= 256 * 8 bytes.
170func fm_build_c(bwt: *u8, n: i64, c_out: *i64) -> i64 {
171 if n <= 0 { return -1 }
172
173 // Histogram of bytes in bwt.
174 var i: i64 = 0
175 while i < 256 {
176 c_out[i] = 0
177 i = i + 1
178 }
179 var j: i64 = 0
180 while j < n {
181 let c: i64 = bwt[j] & 0xff
182 c_out[c] = c_out[c] + 1
183 j = j + 1
184 }
185
186 // Prefix-sum into a "less than c" cumulative.
187 var prev: i64 = 0
188 var k: i64 = 0
189 while k < 256 {
190 let cnt: i64 = c_out[k]
191 c_out[k] = prev
192 prev = prev + cnt
193 k = k + 1
194 }
195 return 0
196}
197
198// Occ(c, i) = count of byte c in bwt[0..i] (inclusive).
199// If i < 0 returns 0 (empty range).
200func fm_occ(bwt: *u8, n: i64, c: i64, i: i64) -> i64 {
201 if i < 0 { return 0 }
202 var count: i64 = 0
203 let cc: i64 = c & 0xff
204 var k: i64 = 0
205 while k <= i {
206 if (bwt[k] & 0xff) == cc { count = count + 1 }
207 k = k + 1
208 }
209 return count
210}
211
212// Count occurrences of pattern p (len plen) in original text via
213// backward search on the FM-index (bwt, n, c_arr).
214// Returns 0..n (occurrence count) or -1 on bad input.
215func fm_count(bwt: *u8, n: i64, c_arr: *i64, p: *u8, plen: i64) -> i64 {
216 if n <= 0 { return -1 }
217 if plen <= 0 { return -1 }
218
219 var i: i64 = plen - 1
220 let c0: i64 = p[i] & 0xff
221 var lo: i64 = c_arr[c0]
222 // Upper bound: C[c+1] - 1; if c==255 use n-1.
223 var hi: i64 = 0
224 if c0 == 255 { hi = n - 1 } else { hi = c_arr[c0 + 1] - 1 }
225 if lo > hi { return 0 }
226
227 i = plen - 2
228 while i >= 0 {
229 let c: i64 = p[i] & 0xff
230 let occ_lo_m1: i64 = fm_occ(bwt, n, c, lo - 1)
231 let occ_hi: i64 = fm_occ(bwt, n, c, hi)
232 lo = c_arr[c] + occ_lo_m1
233 hi = c_arr[c] + occ_hi - 1
234 if lo > hi { return 0 }
235 i = i - 1
236 }
237
238 return hi - lo + 1
239}