nx_align_match.nx source
↩ module page · 215 lines · 8201 B
1// nx_align_match.nx -- minimizer-based seed match (query vs reference).
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/li-2018-minimap2
5//
6// G1.2 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Composes the
7// minimizer-extraction primitive (G1.1) into the seeding stage of
8// read alignment.
9//
10// Pipeline:
11// reference: bases -> minimizer_extract -> (r_vals, r_pos)
12// query : bases -> minimizer_extract -> (q_vals, q_pos)
13// seeds : minimizer_match(q, r) -> {(q_pos, r_pos) : q_val == r_val}
14//
15// Each emitted (q_pos, r_pos) pair is a candidate seed: a k-mer of
16// the query that occurs at r_pos in the reference. Downstream
17// (G1.3 chaining + G1.4 SW extension) groups compatible seeds into
18// chains and extends them into full alignments.
19//
20// Naive O(q_count * r_count) impl for G1.2. For real-world use
21// the reference minimizer index is sorted-by-value or hashed so
22// match is O(q_count * log r_count) or O(q_count) -- those live in
23// nx_align_match_indexed.nx (G1.4).
24//
25// API:
26// minimizer_match(q_vals, q_pos, q_count,
27// r_vals, r_pos, r_count,
28// out_q_pos, out_r_pos, max_out) -> i64 pair count
29//
30// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
31// intended_use: "Seed match in read-to-reference alignment;
32// composes with nx_align Smith-Waterman to
33// complete the seed-and-extend pipeline"
34// sil_target: SIL2
35// asil_target: QM
36// dal_target: DAL C
37// iec_62304_class: B
38// evidence: [no_floating_point, deterministic,
39// bit_equal_reproducible,
40// composes_nx_align_minimizer_KAT,
41// nested_loop_emission_order_stable,
42// license_tier_INDEPENDENT_REDERIVE]
43// hazard_register: [bug-tape-seed-explosion-repetitive-region,
44// bug-tape-seed-strand-not-tracked,
45// bug-tape-seed-match-misses-reverse-strand,
46// bug-tape-seed-pair-capacity-overrun]
47// residual_risk: "Naive O(q*r) is fine for KAT and short reads
48// against short references; high-multiplicity
49// minimizers (telomeric / centromeric repeats)
50// produce O(n^2)-class seed counts and need
51// the multiplicity-capped indexed variant
52// (G1.4). Strand-aware seeding (canonical
53// already strand-invariant; explicit +/- tag
54// for downstream chain orientation) lives in
55// nx_align_match_stranded.nx (G1.3)."
56// verdict: NOT_YET_EVALUATED
57
58import "nx_syscalls.nx"
59
60// ============================================================
61// G1.4a -- indexed (binary-search) matcher.
62// Replaces O(q*r) naive scan with O(q*log r) per query minimizer
63// when the reference minimizer array is pre-sorted by value.
64//
65// Caller responsibility: r_sorted_vals[] is sorted ascending; the
66// r_sorted_pos[] array tracks corresponding positions in lockstep
67// (i.e., r_sorted_pos[i] is the original reference position for the
68// minimizer value r_sorted_vals[i]).
69//
70// Sort key is value-only; same-value entries may be in any order
71// as long as the two arrays stay synchronised. For a value with
72// multiple positions, emission follows the sorted-array order.
73// ============================================================
74
75// Companion to minimizer_match_indexed: in-place insertion sort of
76// paired (vals, positions) arrays by val ascending. Position-ties
77// preserved in input order (stable). Removes the "caller hand-sorts"
78// precondition from minimizer_match_indexed; together they close the
79// perf path of the seed-and-extend pipeline.
80//
81// O(n^2) reference impl per the four-pillar perf rule -- fine for
82// per-read query minimizers (n <= ~500). For reference-genome scale
83// (millions of minimizers), the radix-sort fast path (G1.5) is queued.
84func sort_minimizers_by_value(vals: *i64, positions: *i64, n: i64) -> i64 {
85 if n <= 1 { return 0 }
86
87 var i: i64 = 1
88 while i < n {
89 let key_v: i64 = vals[i]
90 let key_p: i64 = positions[i]
91 var j: i64 = i - 1
92 var keep_going: i64 = 1
93 while keep_going == 1 {
94 if j < 0 {
95 keep_going = 0
96 } else {
97 // Shift if vals[j] > key_v (strict; stable on ties).
98 if vals[j] > key_v {
99 vals[j + 1] = vals[j]
100 positions[j + 1] = positions[j]
101 j = j - 1
102 } else {
103 keep_going = 0
104 }
105 }
106 }
107 vals[j + 1] = key_v
108 positions[j + 1] = key_p
109 i = i + 1
110 }
111 return 0
112}
113
114// Internal: find lower-bound index of target in sorted_vals[0..n).
115// Returns the smallest idx such that sorted_vals[idx] >= target,
116// or n if all values are smaller.
117func nx_bsearch_lower(sorted_vals: *i64, n: i64, target: i64) -> i64 {
118 var lo: i64 = 0
119 var hi: i64 = n
120 while lo < hi {
121 let mid: i64 = (lo + hi) / 2
122 if sorted_vals[mid] < target {
123 lo = mid + 1
124 } else {
125 hi = mid
126 }
127 }
128 return lo
129}
130
131// Indexed matcher. Output (q_pos, r_pos) pairs for each query
132// minimizer that matches a sorted-reference minimizer.
133// Same return semantics as minimizer_match: positive = pair count,
134// 0 = no pairs, -1 = overflow / bad input.
135func minimizer_match_indexed(q_vals: *i64, q_pos: *i64, q_count: i64,
136 r_sorted_vals: *i64, r_sorted_pos: *i64, r_count: i64,
137 out_q_pos: *i64, out_r_pos: *i64,
138 max_out: i64) -> i64 {
139 if q_count < 0 { return -1 }
140 if r_count < 0 { return -1 }
141 if max_out <= 0 {
142 if q_count == 0 { return 0 }
143 if r_count == 0 { return 0 }
144 return -1
145 }
146
147 var emitted: i64 = 0
148 var q: i64 = 0
149 while q < q_count {
150 let qv: i64 = q_vals[q]
151 let qp: i64 = q_pos[q]
152 let start_idx: i64 = nx_bsearch_lower(r_sorted_vals, r_count, qv)
153 var r: i64 = start_idx
154 var keep_going: i64 = 1
155 while keep_going == 1 {
156 if r >= r_count {
157 keep_going = 0
158 } else {
159 if r_sorted_vals[r] != qv {
160 keep_going = 0
161 } else {
162 if emitted >= max_out { return -1 }
163 out_q_pos[emitted] = qp
164 out_r_pos[emitted] = r_sorted_pos[r]
165 emitted = emitted + 1
166 r = r + 1
167 }
168 }
169 }
170 q = q + 1
171 }
172 return emitted
173}
174
175
176// Match query minimizers against reference minimizers.
177//
178// Emission order: outer loop over query (q = 0..q_count-1), inner
179// loop over reference (r = 0..r_count-1). For each q, all matching
180// (q_pos, r_pos) pairs are emitted in reference-order. Deterministic
181// across hosts.
182//
183// Returns: count of pairs written; -1 if capacity exceeded or
184// negative count parameters.
185func minimizer_match(q_vals: *i64, q_pos: *i64, q_count: i64,
186 r_vals: *i64, r_pos: *i64, r_count: i64,
187 out_q_pos: *i64, out_r_pos: *i64,
188 max_out: i64) -> i64 {
189 if q_count < 0 { return -1 }
190 if r_count < 0 { return -1 }
191 if max_out <= 0 {
192 if q_count == 0 { return 0 }
193 if r_count == 0 { return 0 }
194 return -1
195 }
196
197 var emitted: i64 = 0
198 var q: i64 = 0
199 while q < q_count {
200 let qv: i64 = q_vals[q]
201 let qp: i64 = q_pos[q]
202 var r: i64 = 0
203 while r < r_count {
204 if r_vals[r] == qv {
205 if emitted >= max_out { return -1 }
206 out_q_pos[emitted] = qp
207 out_r_pos[emitted] = r_pos[r]
208 emitted = emitted + 1
209 }
210 r = r + 1
211 }
212 q = q + 1
213 }
214 return emitted
215}