nx_arousal_fit.nx source
↩ module page · 202 lines · 9672 B
1// nx_arousal_fit.nx -- fits the perfusion time-constant tau to OBSERVED data by coarse-to-fine
2// sweep over the shared nx_arousal_lib objective. This is the rung that lets the twin EXCEED the
3// published literature instead of inheriting its ceiling: the orgasm/perfusion parameters rest on
4// tiny old samples (Bohlen n=11, 1982), so the honest move is to make them FITTABLE, not fixed.
5//
6// Reports seed_residual alongside the fitted one so a run that did NOT improve is visible rather
7// than silently reported as a fit. A fitter that always claims success measures nothing.
8//
9// usage: nx_arousal_fit tau <drive> <t1> <o1> <t2> <o2> [lo_ms] [hi_ms]
10// nx_arousal_fit selftest
11// license_tier: ORIGINAL expect_exit: 0
12import "nx_arousal_lib.nx"
13const F_MAGIC_4000000: i64 = 4000000
14
15const F_ROUNDS: i64 = 8
16const F_SAMPLES: i64 = 16
17// per-sample residual slack defining the confidence bracket: the set of tau that fit within an
18// average of 1 permil per sample of the best. On noisy data this bracket WIDENS -- that widening
19// is the honest signal. A point estimate on noisy input is false precision.
20const F_CI_SLACK: i64 = 1
21const F_LO_DEFAULT: i64 = 1000
22const F_HI_DEFAULT: i64 = 1000000
23
24// one coarse pass: scan [lo,hi] at `step`, return the best tau found (residual via out-param)
25func f_scan(lo: i64, hi: i64, step: i64, drive: i64, t1: i64, o1: i64, t2: i64, o2: i64, best: *i64) -> i64 {
26 var tau: i64 = lo
27 var btau: i64 = lo
28 var bres: i64 = A_RESMAX
29 while tau <= hi {
30 let r: i64 = a_residual(tau, drive, t1, o1, t2, o2)
31 if r < bres { bres = r; btau = tau }
32 tau = tau + step
33 }
34 best[0] = bres
35 return btau
36}
37
38func main(argc: i64, argv: *i64) -> i64 {
39 if argc < 2 { a_puts("REFUSED reason=no_verb\n" as *u8); return 2 }
40 let verb: *u8 = argv[1] as *u8
41
42 if verb[0] == (115 as u8) {
43 if verb[1] == (101 as u8) {
44 if verb[2] == (114 as u8) {
45 // series <path> [tlo] [thi] -- JOINT (tau, drive) fit over a dense trajectory.
46 // This is the shape real thermography emits: a sample per frame, with drive unknown.
47 if argc < 3 { a_puts("REFUSED reason=series_needs_path\n" as *u8); return 2 }
48 let sl: *i64 = sys_mmap(16) as *i64
49 sl[0] = 0
50 let sbuf: *u8 = sys_read_file(argv[2] as *u8, sl)
51 let sn: i64 = sl[0]
52 if (sbuf as i64) == 0 { a_puts("REFUSED reason=series_unreadable\n" as *u8); return 4 }
53 if sn <= 0 { a_puts("REFUSED reason=series_empty\n" as *u8); return 4 }
54 var tlo: i64 = F_LO_DEFAULT
55 var thi: i64 = F_HI_DEFAULT
56 if argc >= 4 { tlo = a_atoi(argv[3] as *u8) }
57 if argc >= 5 { thi = a_atoi(argv[4] as *u8) }
58 if tlo <= 0 { a_puts("REFUSED reason=nonpositive_lo\n" as *u8); return 3 }
59 if thi <= tlo { a_puts("REFUSED reason=empty_bracket\n" as *u8); return 3 }
60 let mt: i64 = tlo + (thi - tlo) / 2
61 let seed: i64 = a_series_residual(sbuf, sn, mt, 500)
62 // a malformed or non-ascending series is refused, never fitted
63 if seed >= A_RESMAX { a_puts("REFUSED reason=malformed_series\n" as *u8); return 4 }
64 var dlo: i64 = 1
65 var dhi: i64 = 1000
66 var btau: i64 = mt
67 var bdr: i64 = 500
68 var bres: i64 = A_RESMAX
69 var r: i64 = 0
70 while r < F_ROUNDS {
71 var ts: i64 = (thi - tlo) / F_SAMPLES
72 if ts < 1 { ts = 1 }
73 var ds: i64 = (dhi - dlo) / F_SAMPLES
74 if ds < 1 { ds = 1 }
75 var tt: i64 = tlo
76 while tt <= thi {
77 var dd: i64 = dlo
78 while dd <= dhi {
79 let rr: i64 = a_series_residual(sbuf, sn, tt, dd)
80 if rr < bres { bres = rr; btau = tt; bdr = dd }
81 dd = dd + ds
82 }
83 tt = tt + ts
84 }
85 var nl: i64 = btau - ts
86 if nl < 1 { nl = 1 }
87 tlo = nl
88 thi = btau + ts
89 var ndl: i64 = bdr - ds
90 if ndl < 1 { ndl = 1 }
91 dlo = ndl
92 dhi = bdr + ds
93 if dhi > 1000 { dhi = 1000 }
94 r = r + 1
95 }
96 if bres >= A_RESMAX { a_puts("REFUSED reason=no_admissible_fit\n" as *u8); return 5 }
97 if bres < 0 { a_puts("REFUSED reason=negative_residual\n" as *u8); return 5 }
98 // CONFIDENCE BRACKET: walk out from the optimum until residual exceeds the
99 // per-sample slack. Clean data -> tight bracket; noisy data -> wide bracket.
100 let pts: i64 = a_series_count(sbuf, sn)
101 if pts <= 0 { a_puts("REFUSED reason=no_points\n" as *u8); return 4 }
102 let thr: i64 = bres + pts * F_CI_SLACK
103 var cstep: i64 = btau / 200
104 if cstep < 1 { cstep = 1 }
105 var cilo: i64 = btau
106 var w1: i64 = 0
107 while w1 == 0 {
108 let cand: i64 = cilo - cstep
109 if cand < 1 { w1 = 1 } else {
110 if a_series_residual(sbuf, sn, cand, bdr) <= thr { cilo = cand } else { w1 = 1 }
111 }
112 }
113 var cihi: i64 = btau
114 var w2: i64 = 0
115 while w2 == 0 {
116 let cand2: i64 = cihi + cstep
117 if cand2 > F_MAGIC_4000000 { w2 = 1 } else {
118 if a_series_residual(sbuf, sn, cand2, bdr) <= thr { cihi = cand2 } else { w2 = 1 }
119 }
120 }
121 let ciw: i64 = cihi - cilo
122 if ciw < 0 { a_puts("REFUSED reason=negative_ci_width\n" as *u8); return 5 }
123 var imp: i64 = 0
124 if bres <= seed { imp = 1 }
125 a_puts("ok best_tau=" as *u8); a_putn(btau)
126 a_puts(" best_drive=" as *u8); a_putn(bdr)
127 a_puts(" residual=" as *u8); a_putn(bres)
128 a_puts(" seed_residual=" as *u8); a_putn(seed)
129 a_puts(" improved=" as *u8); a_putn(imp)
130 a_puts(" pts=" as *u8); a_putn(pts)
131 a_puts(" tau_lo=" as *u8); a_putn(cilo)
132 a_puts(" tau_hi=" as *u8); a_putn(cihi)
133 a_puts(" ci_width=" as *u8); a_putn(ciw)
134 a_puts("\n" as *u8)
135 return 0
136 }
137 }
138 a_puts("AROUSAL-FIT-SELFTEST rounds=" as *u8); a_putn(F_ROUNDS)
139 a_puts(" samples=" as *u8); a_putn(F_SAMPLES)
140 a_puts("\n" as *u8)
141 return 0
142 }
143 if verb[0] != (116 as u8) { a_puts("REFUSED reason=unknown_verb\n" as *u8); return 2 }
144 if argc < 7 { a_puts("REFUSED reason=tau_needs_5_args\n" as *u8); return 2 }
145
146 let drive: i64 = a_atoi(argv[2] as *u8)
147 let t1: i64 = a_atoi(argv[3] as *u8)
148 let o1: i64 = a_atoi(argv[4] as *u8)
149 let t2: i64 = a_atoi(argv[5] as *u8)
150 let o2: i64 = a_atoi(argv[6] as *u8)
151 var lo: i64 = F_LO_DEFAULT
152 var hi: i64 = F_HI_DEFAULT
153 if argc >= 8 { lo = a_atoi(argv[7] as *u8) }
154 if argc >= 9 { hi = a_atoi(argv[8] as *u8) }
155
156 if drive <= 0 { a_puts("REFUSED reason=nonpositive_drive\n" as *u8); return 3 }
157 if drive > 1000 { a_puts("REFUSED reason=drive_over_permil\n" as *u8); return 3 }
158 if lo <= 0 { a_puts("REFUSED reason=nonpositive_lo\n" as *u8); return 3 }
159 if hi <= lo { a_puts("REFUSED reason=empty_bracket\n" as *u8); return 3 }
160 if t1 < 0 { a_puts("REFUSED reason=negative_time\n" as *u8); return 3 }
161 if t2 < 0 { a_puts("REFUSED reason=negative_time\n" as *u8); return 3 }
162 // an observation is a permil perfusion: outside 0..1000 it is not a measurement
163 if o1 < 0 { a_puts("REFUSED reason=obs_out_of_range\n" as *u8); return 3 }
164 if o2 < 0 { a_puts("REFUSED reason=obs_out_of_range\n" as *u8); return 3 }
165 if o1 > 1000 { a_puts("REFUSED reason=obs_out_of_range\n" as *u8); return 3 }
166 if o2 > 1000 { a_puts("REFUSED reason=obs_out_of_range\n" as *u8); return 3 }
167
168 // seed = residual at the bracket midpoint, BEFORE any fitting. This is the honesty baseline.
169 let mid: i64 = lo + (hi - lo) / 2
170 let seed: i64 = a_residual(mid, drive, t1, o1, t2, o2)
171
172 let bp: *i64 = sys_mmap(16) as *i64
173 bp[0] = A_RESMAX
174 var btau: i64 = mid
175 var bres: i64 = A_RESMAX
176 var clo: i64 = lo
177 var chi: i64 = hi
178 var r: i64 = 0
179 while r < F_ROUNDS {
180 var step: i64 = (chi - clo) / F_SAMPLES
181 if step < 1 { step = 1 }
182 let t: i64 = f_scan(clo, chi, step, drive, t1, o1, t2, o2, bp)
183 if bp[0] < bres { bres = bp[0]; btau = t }
184 var nlo: i64 = t - step
185 if nlo < 1 { nlo = 1 }
186 clo = nlo
187 chi = t + step
188 r = r + 1
189 }
190
191 if bres >= A_RESMAX { a_puts("REFUSED reason=no_admissible_tau\n" as *u8); return 5 }
192 // a residual is a DISTANCE -- refuse a corrupt statistic, never report it as a fit
193 if bres < 0 { a_puts("REFUSED reason=negative_residual\n" as *u8); return 5 }
194 var improved: i64 = 0
195 if bres <= seed { improved = 1 }
196 a_puts("ok best_tau=" as *u8); a_putn(btau)
197 a_puts(" residual=" as *u8); a_putn(bres)
198 a_puts(" seed_residual=" as *u8); a_putn(seed)
199 a_puts(" improved=" as *u8); a_putn(improved)
200 a_puts("\n" as *u8)
201 return 0
202}