nx_align_match_test.nx source

1// nx_align_match_test.nx -- end-to-end seed-pipeline KAT. 2// 3// Builds the full chain: pack reference + query DNA, extract 4// minimizers from each, match query-vs-reference, assert exact 5// (q_pos, r_pos) pairs come out. Demonstrates G1.0+G1.1+G1.2 6// + nx_sequence composition as a single coherent pipeline. 7// 8// Reference: "ACGTACGTAC" (10 bases), k=4 w=3 9// minimizers: (0x1B, 0), (0x6C, 1), (0x1B, 4) 10// 11// Query: "GTACGT" (6 bases), k=4 w=3 12// k-mers / canonical at pos 0..2: 13// pos 0: GTAC = 0xB1 (palindrome) 14// pos 1: TACG = 0xC6 / rc CGTA 0x6C -> canon 0x6C 15// pos 2: ACGT = 0x1B (palindrome) 16// 1 window, leftmost min of [0xB1, 0x6C, 0x1B] = 0x1B at pos 2. 17// minimizers: (0x1B, 2) 18// 19// Match: query value 0x1B hits reference 0x1B at r_pos 0 and r_pos 4. 20// pairs: (q=2, r=0), (q=2, r=4) 21// count = 2 22// 23// expect_exit: 0 24// 25// license_tier: ORIGINAL 26 27import "nx_syscalls.nx" 28import "nx_sequence.nx" 29import "nx_align_minimizer.nx" 30import "nx_align_match.nx" 31 32func main() -> i64 { 33 34 // -------- pack reference "ACGTACGTAC" -------- 35 let r_bases: *u8 = sys_mmap(8) 36 r_bases[0] = 0x1B // ACGT 37 r_bases[1] = 0x1B // ACGT 38 r_bases[2] = 0x10 // AC + padding 39 40 let r_vals: *i64 = sys_mmap(128) as *i64 41 let r_pos: *i64 = sys_mmap(128) as *i64 42 let r_cnt: i64 = minimizer_extract(r_bases, 10, 4, 3, r_vals, r_pos, 16) 43 if r_cnt != 3 { return 1 } 44 45 // Sanity-check reference minimizers (hand-traced). 46 if r_vals[0] != 0x1B { return 2 } 47 if r_pos[0] != 0 { return 3 } 48 if r_vals[1] != 0x6C { return 4 } 49 if r_pos[1] != 1 { return 5 } 50 if r_vals[2] != 0x1B { return 6 } 51 if r_pos[2] != 4 { return 7 } 52 53 // -------- pack query "GTACGT" -------- 54 // bases: G T A C G T 55 // pack: byte 0 = 10 11 00 01 = 0xB1 (GTAC) 56 // byte 1 = 10 11 00 00 = 0xB0 (GT + padding) 57 let q_bases: *u8 = sys_mmap(4) 58 q_bases[0] = 0xB1 59 q_bases[1] = 0xB0 60 61 let q_vals: *i64 = sys_mmap(64) as *i64 62 let q_pos: *i64 = sys_mmap(64) as *i64 63 let q_cnt: i64 = minimizer_extract(q_bases, 6, 4, 3, q_vals, q_pos, 16) 64 if q_cnt != 1 { return 10 } 65 if q_vals[0] != 0x1B { return 11 } 66 if q_pos[0] != 2 { return 12 } 67 68 // -------- match -------- 69 let out_q: *i64 = sys_mmap(128) as *i64 70 let out_r: *i64 = sys_mmap(128) as *i64 71 let n_pairs: i64 = minimizer_match(q_vals, q_pos, q_cnt, 72 r_vals, r_pos, r_cnt, 73 out_q, out_r, 16) 74 if n_pairs != 2 { return 20 } 75 76 // Emission order: outer q, inner r. For q[0] = (0x1B, 2) the 77 // inner walk hits r[0] = (0x1B, 0) first then r[2] = (0x1B, 4). 78 if out_q[0] != 2 { return 21 } 79 if out_r[0] != 0 { return 22 } 80 if out_q[1] != 2 { return 23 } 81 if out_r[1] != 4 { return 24 } 82 83 // -------- no-match case: query that shares no minimizer -------- 84 // Query "TTTT" -> not enough bases for k=4 w=3 (needs n>=6). 85 // Use query "TTTTTT" instead. 86 let q2_bases: *u8 = sys_mmap(4) 87 q2_bases[0] = 0xFF // TTTT 88 q2_bases[1] = 0xF0 // TT + padding 89 90 let q2_vals: *i64 = sys_mmap(64) as *i64 91 let q2_pos: *i64 = sys_mmap(64) as *i64 92 let q2_cnt: i64 = minimizer_extract(q2_bases, 6, 4, 3, q2_vals, q2_pos, 16) 93 if q2_cnt <= 0 { return 30 } 94 // TTTT canonical = revcomp("TTTT") = "AAAA" = 0, so canon = 0. 95 if q2_vals[0] != 0 { return 31 } 96 97 let n2: i64 = minimizer_match(q2_vals, q2_pos, q2_cnt, 98 r_vals, r_pos, r_cnt, 99 out_q, out_r, 16) 100 if n2 != 0 { return 32 } // reference has 0x1B / 0x6C, query has 0, no hits 101 102 // -------- capacity overflow surfaces -1 -------- 103 let n3: i64 = minimizer_match(q_vals, q_pos, q_cnt, 104 r_vals, r_pos, r_cnt, 105 out_q, out_r, 1) 106 if n3 != -1 { return 40 } 107 108 // -------- empty inputs -------- 109 let n4: i64 = minimizer_match(q_vals, q_pos, 0, 110 r_vals, r_pos, r_cnt, 111 out_q, out_r, 16) 112 if n4 != 0 { return 50 } 113 114 let n5: i64 = minimizer_match(q_vals, q_pos, q_cnt, 115 r_vals, r_pos, 0, 116 out_q, out_r, 16) 117 if n5 != 0 { return 51 } 118 119 // ============================================================ 120 // Section G -- nx_bsearch_lower helper. 121 // sorted_vals = [10, 20, 20, 30, 40] 122 // target=10 -> lower-bound idx 0 123 // target=20 -> lower-bound idx 1 (first 20) 124 // target=25 -> lower-bound idx 3 (first > 25) 125 // target=40 -> idx 4 126 // target=99 -> idx 5 (==n, all smaller) 127 // target=5 -> idx 0 128 // ============================================================ 129 130 let sv: *i64 = sys_mmap(64) as *i64 131 sv[0]=10; sv[1]=20; sv[2]=20; sv[3]=30; sv[4]=40 132 133 if nx_bsearch_lower(sv, 5, 10) != 0 { return 60 } 134 if nx_bsearch_lower(sv, 5, 20) != 1 { return 61 } 135 if nx_bsearch_lower(sv, 5, 25) != 3 { return 62 } 136 if nx_bsearch_lower(sv, 5, 40) != 4 { return 63 } 137 if nx_bsearch_lower(sv, 5, 99) != 5 { return 64 } 138 if nx_bsearch_lower(sv, 5, 5) != 0 { return 65 } 139 140 // ============================================================ 141 // Section H -- minimizer_match_indexed. 142 // 143 // Use the same reference minimizers as Section E ("ACGTACGTAC" 144 // k=4 w=3), but pre-sorted by value: 145 // Original (in position order): 146 // (0x1B, 0), (0x6C, 1), (0x1B, 4) 147 // Sorted by value: 148 // (0x1B, 0), (0x1B, 4), (0x6C, 1) 149 // 150 // Query: same single minimizer (0x1B, 2) -- same as Section E. 151 // Expected pairs: (q=2, r=0), (q=2, r=4) -- count 2. 152 // ============================================================ 153 154 let r_sv: *i64 = sys_mmap(64) as *i64 155 let r_sp: *i64 = sys_mmap(64) as *i64 156 r_sv[0]=0x1B; r_sp[0]=0 157 r_sv[1]=0x1B; r_sp[1]=4 158 r_sv[2]=0x6C; r_sp[2]=1 159 160 let oqi: *i64 = sys_mmap(128) as *i64 161 let ori: *i64 = sys_mmap(128) as *i64 162 let n_ix: i64 = minimizer_match_indexed(q_vals, q_pos, q_cnt, 163 r_sv, r_sp, 3, 164 oqi, ori, 16) 165 if n_ix != 2 { return 70 } 166 if oqi[0] != 2 { return 71 } 167 if ori[0] != 0 { return 72 } 168 if oqi[1] != 2 { return 73 } 169 if ori[1] != 4 { return 74 } 170 171 // No match: query with value not in reference. 172 let q_miss_v: *i64 = sys_mmap(64) as *i64 173 let q_miss_p: *i64 = sys_mmap(64) as *i64 174 q_miss_v[0] = 0x99 175 q_miss_p[0] = 10 176 let n_miss: i64 = minimizer_match_indexed(q_miss_v, q_miss_p, 1, 177 r_sv, r_sp, 3, 178 oqi, ori, 16) 179 if n_miss != 0 { return 80 } 180 181 // Capacity overflow -> -1. 182 let n_ov: i64 = minimizer_match_indexed(q_vals, q_pos, q_cnt, 183 r_sv, r_sp, 3, 184 oqi, ori, 1) 185 if n_ov != -1 { return 90 } 186 187 // Empty reference -> 0. 188 let n_er: i64 = minimizer_match_indexed(q_vals, q_pos, q_cnt, 189 r_sv, r_sp, 0, 190 oqi, ori, 16) 191 if n_er != 0 { return 100 } 192 193 // ============================================================ 194 // Section I -- sort_minimizers_by_value. 195 // Input vals = [0x6C, 0x1B, 0xB1, 0x1B, 0x06] 196 // pos = [ 1, 0, 3, 4, 6] 197 // Sorted by val (position-stable on ties): 198 // vals = [0x06, 0x1B, 0x1B, 0x6C, 0xB1] 199 // pos = [ 6, 0, 4, 1, 3] 200 // ============================================================ 201 202 let sv2: *i64 = sys_mmap(64) as *i64 203 let sp2: *i64 = sys_mmap(64) as *i64 204 sv2[0]=0x6C; sp2[0]=1 205 sv2[1]=0x1B; sp2[1]=0 206 sv2[2]=0xB1; sp2[2]=3 207 sv2[3]=0x1B; sp2[3]=4 208 sv2[4]=0x06; sp2[4]=6 209 210 sort_minimizers_by_value(sv2, sp2, 5) 211 212 if sv2[0] != 0x06 { return 110 } 213 if sp2[0] != 6 { return 111 } 214 if sv2[1] != 0x1B { return 112 } 215 if sp2[1] != 0 { return 113 } // stable: 0x1B from idx 1 (pos=0) comes first 216 if sv2[2] != 0x1B { return 114 } 217 if sp2[2] != 4 { return 115 } // stable: 0x1B from idx 3 (pos=4) comes second 218 if sv2[3] != 0x6C { return 116 } 219 if sp2[3] != 1 { return 117 } 220 if sv2[4] != 0xB1 { return 118 } 221 if sp2[4] != 3 { return 119 } 222 223 // No-op on n=0 / n=1. 224 if sort_minimizers_by_value(sv2, sp2, 0) != 0 { return 120 } 225 if sort_minimizers_by_value(sv2, sp2, 1) != 0 { return 121 } 226 227 // Already-sorted preserved. 228 sort_minimizers_by_value(sv2, sp2, 5) 229 if sv2[0] != 0x06 { return 122 } 230 if sv2[4] != 0xB1 { return 123 } 231 232 // End-to-end composition: build sorted ref then run indexed match. 233 // This mirrors how a real caller would integrate the perf path. 234 let sv3: *i64 = sys_mmap(64) as *i64 235 let sp3: *i64 = sys_mmap(64) as *i64 236 sv3[0]=0x1B; sp3[0]=0 237 sv3[1]=0x6C; sp3[1]=1 238 sv3[2]=0x1B; sp3[2]=4 239 sort_minimizers_by_value(sv3, sp3, 3) 240 // After sort: vals=[0x1B, 0x1B, 0x6C], pos=[0, 4, 1]. 241 let n_e2e: i64 = minimizer_match_indexed(q_vals, q_pos, q_cnt, 242 sv3, sp3, 3, 243 oqi, ori, 16) 244 if n_e2e != 2 { return 130 } 245 if oqi[0] != 2 { return 131 } 246 if ori[0] != 0 { return 132 } 247 if oqi[1] != 2 { return 133 } 248 if ori[1] != 4 { return 134 } 249 250 return 0 251}