nx_crater_quality.nx source
↩ module page · 263 lines · 10121 B
1// nx_crater_quality.nx -- per-kind grader for nx_crater_field output.
2//
3// Reads a crater array (4 i64 per crater: cx, cy, r, depth) and grades
4// on 5 axes:
5//
6// 0. POWER_LAW_FIT: Neukum 1983 N(>D) ~ D^-2. We check that the
7// count above the median radius is roughly 25% of the total
8// (consistent with -2 exponent).
9// 1. SPATIAL_UNIFORMITY: nearest-neighbour distance distribution
10// suggests Poisson (vs. clustered). Sample N pairs.
11// 2. SIZE_RANGE: max_r / min_r should be >= 4 (real cratering covers
12// multiple decades of size).
13// 3. DEPTH_RATIO: mean(depth/r) should match Pike 1977 value 0.20
14// (fresh) or be lower (eroded); tolerance band.
15// 4. NON_OVERLAP: fraction of craters whose centres are NOT inside
16// another crater (catches generators that emit duplicates).
17//
18// EMITS LAYER_VERDICT (kind = NX_LAYER_KIND_CRATER).
19
20// nx_safety_envelope:
21// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
22// sil_target: SIL1
23// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
24// verdict: NOT_YET_EVALUATED
25
26import "nx_syscalls.nx"
27import "nx_tier.nx"
28import "nx_layer_verdict.nx"
29
30const NX_CR_Q: nx_int = 16384
31
32const NX_CR_AXIS_POWER_LAW: nx_int = 0
33const NX_CR_AXIS_SPATIAL_UNIF: nx_int = 1
34const NX_CR_AXIS_SIZE_RANGE: nx_int = 2
35const NX_CR_AXIS_DEPTH_RATIO: nx_int = 3
36const NX_CR_AXIS_NON_OVERLAP: nx_int = 4
37const NX_CR_AXIS_COUNT: nx_int = 5
38
39const NX_CR_CRATER_STRIDE: nx_int = 4
40const NX_CR_OFF_X: nx_int = 0
41const NX_CR_OFF_Y: nx_int = 1
42const NX_CR_OFF_R: nx_int = 2
43const NX_CR_OFF_DEPTH: nx_int = 3
44
45func _cr_band_score(val: nx_int, lo_q: nx_int, hi_q: nx_int) -> nx_int {
46 let q: nx_int = NX_CR_Q
47 if val < 0 { return 0 }
48 if val < lo_q {
49 if lo_q > 0 { return (val * q) / lo_q }
50 return 0
51 }
52 if val <= hi_q { return q }
53 let extra: nx_int = val - hi_q
54 let span: nx_int = hi_q - lo_q
55 if span <= 0 { return 0 }
56 var s: nx_int = q - (extra * q) / span
57 if s < 0 { s = 0 }
58 return s
59}
60
61// ===== Power-law axis ==============================================
62// For alpha=2 power law N(>D) ~ D^-2, the count above median D should
63// be 25% of total (since F(D_med) = 0.5). Practical tolerance ~ +/- 5%.
64func _cr_power_law_q14(craters: *i64, n: nx_int) -> nx_int {
65 if n < 4 { return 0 }
66 let q: nx_int = NX_CR_Q
67 // Find median radius via mean as proxy (avoids sort).
68 var sum_r: nx_int = 0
69 var i: nx_int = 0
70 while i < n {
71 sum_r = sum_r + craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_R]
72 i = i + 1
73 }
74 let mean_r: nx_int = sum_r / n
75 var n_above: nx_int = 0
76 var j: nx_int = 0
77 while j < n {
78 if craters[j * NX_CR_CRATER_STRIDE + NX_CR_OFF_R] > mean_r {
79 n_above = n_above + 1
80 }
81 j = j + 1
82 }
83 let frac_q: nx_int = (n_above * q) / n
84 // Target [20%, 30%] for alpha=2 (median is biased toward smaller
85 // craters in power-law, so we look around 25%).
86 return _cr_band_score(frac_q, q / 5, q * 3 / 10)
87}
88
89// ===== Spatial uniformity ==========================================
90// Sample N consecutive crater-pair distances; compare mean to expected
91// Poisson nearest-neighbour distance. Simplified: a tight CV of pair
92// distances indicates uniform distribution.
93func _cr_spatial_unif_q14(craters: *i64, n: nx_int) -> nx_int {
94 if n < 4 { return 0 }
95 let q: nx_int = NX_CR_Q
96 var sum: nx_int = 0
97 var min_d_sq: nx_int = 0 - 1
98 var max_d_sq: nx_int = 0
99 var count: nx_int = 0
100 var i: nx_int = 0
101 while i < n - 1 {
102 let dx: nx_int = craters[(i + 1) * NX_CR_CRATER_STRIDE + NX_CR_OFF_X] -
103 craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_X]
104 let dy: nx_int = craters[(i + 1) * NX_CR_CRATER_STRIDE + NX_CR_OFF_Y] -
105 craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_Y]
106 let d_sq: nx_int = dx * dx + dy * dy
107 sum = sum + d_sq
108 count = count + 1
109 if d_sq > max_d_sq { max_d_sq = d_sq }
110 if min_d_sq < 0 { min_d_sq = d_sq }
111 if d_sq < min_d_sq { min_d_sq = d_sq }
112 i = i + 1
113 }
114 if count == 0 { return 0 }
115 let mean_d_sq: nx_int = sum / count
116 if mean_d_sq <= 0 { return 0 }
117 let cv: nx_int = ((max_d_sq - min_d_sq) * q) / mean_d_sq
118 // Lower cv = more uniform. Score = q - cv / 4 (saturate).
119 var score: nx_int = q - cv / 4
120 if score < 0 { score = 0 }
121 if score > q { score = q }
122 return score
123}
124
125// ===== Size range axis =============================================
126func _cr_size_range_q14(craters: *i64, n: nx_int) -> nx_int {
127 if n < 2 { return 0 }
128 let q: nx_int = NX_CR_Q
129 var rmin: nx_int = craters[NX_CR_OFF_R]
130 var rmax: nx_int = craters[NX_CR_OFF_R]
131 var i: nx_int = 0
132 while i < n {
133 let r: nx_int = craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_R]
134 if r < rmin { rmin = r }
135 if r > rmax { rmax = r }
136 i = i + 1
137 }
138 if rmin <= 0 { return 0 }
139 let ratio: nx_int = rmax / rmin
140 // Target ratio >= 4. Score saturates at ratio=8.
141 if ratio >= 8 { return q }
142 if ratio <= 1 { return 0 }
143 return ((ratio - 1) * q) / 7
144}
145
146// ===== Depth ratio axis ============================================
147// Mean depth / mean radius should be ~ 0.2 (Pike 1977 fresh) or lower
148// (aged). We accept [0.05, 0.25] target band.
149func _cr_depth_ratio_q14(craters: *i64, n: nx_int) -> nx_int {
150 if n <= 0 { return 0 }
151 let q: nx_int = NX_CR_Q
152 var sum_r: nx_int = 0
153 var sum_d: nx_int = 0
154 var i: nx_int = 0
155 while i < n {
156 sum_r = sum_r + craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_R]
157 sum_d = sum_d + craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_DEPTH]
158 i = i + 1
159 }
160 if sum_r <= 0 { return 0 }
161 let ratio_q: nx_int = (sum_d * q) / sum_r
162 return _cr_band_score(ratio_q, q / 20, q / 4) // [0.05, 0.25]
163}
164
165// ===== Non-overlap axis ============================================
166// Fraction of craters whose centre is NOT inside another's bowl.
167func _cr_non_overlap_q14(craters: *i64, n: nx_int) -> nx_int {
168 if n < 2 { return NX_CR_Q }
169 let q: nx_int = NX_CR_Q
170 var n_clear: nx_int = 0
171 var i: nx_int = 0
172 while i < n {
173 let xi: nx_int = craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_X]
174 let yi: nx_int = craters[i * NX_CR_CRATER_STRIDE + NX_CR_OFF_Y]
175 var covered: nx_int = 0
176 var j: nx_int = 0
177 while j < n {
178 if j != i {
179 let xj: nx_int = craters[j * NX_CR_CRATER_STRIDE + NX_CR_OFF_X]
180 let yj: nx_int = craters[j * NX_CR_CRATER_STRIDE + NX_CR_OFF_Y]
181 let rj: nx_int = craters[j * NX_CR_CRATER_STRIDE + NX_CR_OFF_R]
182 let dx: nx_int = xi - xj
183 let dy: nx_int = yi - yj
184 if dx * dx + dy * dy < rj * rj { covered = 1; j = n }
185 }
186 j = j + 1
187 }
188 if covered == 0 { n_clear = n_clear + 1 }
189 i = i + 1
190 }
191 return (n_clear * q) / n
192}
193
194// ===== Public: grader ==============================================
195func nx_crater_quality_grade(
196 craters: *i64, n: nx_int, out_verdict: *i64
197) {
198 let pl: nx_int = _cr_power_law_q14(craters, n)
199 let su: nx_int = _cr_spatial_unif_q14(craters, n)
200 let sr: nx_int = _cr_size_range_q14(craters, n)
201 let dr: nx_int = _cr_depth_ratio_q14(craters, n)
202 let no: nx_int = _cr_non_overlap_q14(craters, n)
203 nx_layer_verdict_init(out_verdict, NX_LAYER_KIND_CRATER,
204 NX_CR_AXIS_COUNT, NX_LAYER_REFINE_MORE_FEATURES)
205 out_verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_POWER_LAW] = pl
206 out_verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_SPATIAL_UNIF] = su
207 out_verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_SIZE_RANGE] = sr
208 out_verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_DEPTH_RATIO] = dr
209 out_verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_NON_OVERLAP] = no
210 nx_layer_verdict_finalize(out_verdict)
211}
212
213// ===== Self-test ====================================================
214func main() -> i64 {
215 let q: nx_int = NX_CR_Q
216 let verdict: *i64 = (sys_mmap(NX_LV_STRIDE * NX_SIZEOF_NX_INT)) as *i64
217
218 // T1: 16 craters with varied radii (1..16), depth = 0.2*r,
219 // positions on a grid. Should pass most axes.
220 let n: nx_int = 16
221 let craters: *i64 = (sys_mmap(n * NX_CR_CRATER_STRIDE * NX_SIZEOF_NX_INT)) as *i64
222 var i: nx_int = 0
223 while i < n {
224 let base: nx_int = i * NX_CR_CRATER_STRIDE
225 // Radius alternates small/large for power-law shape.
226 var r: nx_int = 100 + (i % 4) * 50
227 if i >= n / 2 { r = 100 + (i % 4) * 400 } // some big craters
228 craters[base + NX_CR_OFF_X] = (i % 4) * 200
229 craters[base + NX_CR_OFF_Y] = (i / 4) * 200
230 craters[base + NX_CR_OFF_R] = r
231 craters[base + NX_CR_OFF_DEPTH] = (r * 2) / 10
232 i = i + 1
233 }
234 nx_crater_quality_grade(craters, n, verdict)
235 // Size range axis: ratio of max(1700)/min(100) = 17 -> Q.
236 if verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_SIZE_RANGE] != q {
237 return __syscall(93, 1, 0, 0, 0, 0, 0)
238 }
239 // Depth ratio: 0.2 exact -> Q (in band).
240 if verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_DEPTH_RATIO] != q {
241 return __syscall(93, 2, 0, 0, 0, 0, 0)
242 }
243 // Non-overlap: positions are 200 apart, max radius 1700 -- some
244 // overlap; but caller pre-filters so we just check it's > 0 and
245 // valid range.
246 let no_v: nx_int = verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_NON_OVERLAP]
247 if no_v < 0 { return __syscall(93, 3, 0, 0, 0, 0, 0) }
248 if no_v > q { return __syscall(93, 4, 0, 0, 0, 0, 0) }
249
250 // T2: All-same-radius craters -> size range axis = 0.
251 var j: nx_int = 0
252 while j < n {
253 craters[j * NX_CR_CRATER_STRIDE + NX_CR_OFF_R] = 100
254 craters[j * NX_CR_CRATER_STRIDE + NX_CR_OFF_DEPTH] = 20
255 j = j + 1
256 }
257 nx_crater_quality_grade(craters, n, verdict)
258 if verdict[NX_LV_OFF_AXIS_0 + NX_CR_AXIS_SIZE_RANGE] != 0 {
259 return __syscall(93, 10, 0, 0, 0, 0, 0)
260 }
261
262 return 0
263}