nx_align_chain.nx source
↩ module page · 153 lines · 6063 B
1// nx_align_chain.nx -- co-linear seed-chain dynamic programming.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/li-2018-minimap2-section-2.2
5//
6// G1.3 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Bridges the seed
7// pairs from nx_align_match (G1.2) to the alignment extension input
8// of Smith-Waterman (G1.0).
9//
10// A "chain" is a sequence of seeds (q_i, r_i) that are strictly
11// co-linear under the forward-strand model: both q-positions and
12// r-positions strictly increase along the chain. Co-linearity is
13// what makes a chain a coherent alignment hypothesis -- it forbids
14// the seeds from belonging to two different alignments or to
15// different strands of an inversion.
16//
17// Algorithm (O(n^2) reference DP):
18// - Precondition: caller has pre-sorted seeds by r_pos ascending
19// (and within-tie by q_pos ascending). Ties in r_pos are rare
20// for typical w >= 5 minimizers but possible; the strictly-less
21// check on r still works because the inner DP scans j < i.
22// - chain_len[i] = 1 + max(chain_len[j]) over j < i with q[j] < q[i] AND r[j] < r[i]
23// = 1 otherwise
24// - parent[i] = the j that achieved the max (or -1 if no extension)
25// - The best chain ends at argmax_i chain_len[i] (leftmost-tie wins)
26// - Backtrack from argmax through parent[] to recover the chain
27// - Output indices are in q-increasing (== r-increasing) order
28//
29// Why strict-less on both axes:
30// - Equality in q would mean the same query position contributes
31// two seeds -- a repeated minimizer landing on the same query
32// base, which would not extend a single alignment chain.
33// - Equality in r similarly maps two seeds to the same reference
34// base, which is geometrically incoherent for a single chain.
35//
36// What this G1.3 does NOT do (deferred):
37// - Gap-cost-aware scoring (minimap2 alpha + beta) -- G1.3b
38// - Striped / sparse-dp acceleration to O(n log n) -- G1.5
39// - Strand-aware chaining (forward vs reverse-complement) -- G1.3c
40// - Anchor-band pruning + secondary-chain reporting -- G1.4
41// - Multiplicity-capped seed input -- G1.4
42//
43// API:
44// seed_chain(seeds_q, seeds_r, n,
45// out_indices, max_out) -> i64
46//
47// Returns the length of the longest co-linear chain (>= 0).
48// Writes that chain's seed indices to out_indices in
49// q-increasing order. Returns -1 if max_out cannot hold the
50// chain (or if n < 0).
51//
52// Memory:
53// Internal sys_mmap of 2 * n * 8 bytes (chain_len + parent arrays).
54// For typical short-read alignment (n ~ 50-500 seeds per read)
55// this is 800-8000 bytes; for long-read alignment (n ~ 1000s) the
56// O(n log n) variant in G1.5 reduces both time and space.
57//
58// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
59// intended_use: "Seed-pair chaining -- the inner stage of
60// seed-and-extend that converts raw seed-pair
61// sets into ordered chains for SW extension"
62// sil_target: SIL2
63// asil_target: QM
64// dal_target: DAL C
65// iec_62304_class: B
66// evidence: [no_floating_point, deterministic,
67// bit_equal_reproducible,
68// co_linear_dp_textbook_recurrence,
69// leftmost_tie_break_stable,
70// composes_nx_align_match_KAT,
71// license_tier_INDEPENDENT_REDERIVE]
72// hazard_register: [bug-tape-chain-non-strict-q-double-counts-seed,
73// bug-tape-chain-backtrace-reversed-output,
74// bug-tape-chain-tie-break-rightmost-drift,
75// bug-tape-chain-input-not-r-sorted-silent,
76// bug-tape-chain-n-zero-returns-junk]
77// residual_risk: "Caller pre-sort is a contract not enforced
78// by the DP -- mis-sorted input may produce
79// a shorter-than-optimal chain. Add a debug
80// sortedness-check primitive in G1.3b."
81// verdict: NOT_YET_EVALUATED
82
83import "nx_syscalls.nx"
84
85// Longest co-linear chain via O(n^2) DP.
86// Returns chain length; writes chain indices to out_indices in
87// q-increasing order. Returns -1 on bad input.
88func seed_chain(seeds_q: *i64, seeds_r: *i64, n: i64,
89 out_indices: *i64, max_out: i64) -> i64 {
90 if n < 0 { return -1 }
91 if n == 0 { return 0 }
92 if max_out <= 0 { return -1 }
93
94 let chain_len: *i64 = sys_mmap(n * 8) as *i64
95 let parent: *i64 = sys_mmap(n * 8) as *i64
96
97 var best_len: i64 = 0
98 var best_end: i64 = 0
99
100 var i: i64 = 0
101 while i < n {
102 chain_len[i] = 1
103 parent[i] = -1
104
105 let qi: i64 = seeds_q[i]
106 let ri: i64 = seeds_r[i]
107
108 var j: i64 = 0
109 while j < i {
110 let qj: i64 = seeds_q[j]
111 let rj: i64 = seeds_r[j]
112 if qj < qi {
113 if rj < ri {
114 let candidate: i64 = chain_len[j] + 1
115 if candidate > chain_len[i] {
116 chain_len[i] = candidate
117 parent[i] = j
118 }
119 }
120 }
121 j = j + 1
122 }
123
124 // Track best end (leftmost-tie on chain length).
125 if chain_len[i] > best_len {
126 best_len = chain_len[i]
127 best_end = i
128 }
129 i = i + 1
130 }
131
132 if best_len > max_out { return -1 }
133
134 // Backtrace. Walk parent chain to fill a reverse buffer, then
135 // reverse into out_indices to deliver q-increasing order.
136 let rev: *i64 = sys_mmap(best_len * 8) as *i64
137 var k: i64 = 0
138 var cur: i64 = best_end
139 while cur >= 0 {
140 rev[k] = cur
141 k = k + 1
142 cur = parent[cur]
143 }
144
145 // k == best_len at this point. Reverse into out_indices.
146 var w: i64 = 0
147 while w < best_len {
148 out_indices[w] = rev[best_len - 1 - w]
149 w = w + 1
150 }
151
152 return best_len
153}