nx_arousal_lib.nx source
↩ module page · 330 lines · 14817 B
1// nx_arousal_lib.nx -- pure physiological-arousal math core. NO main(): this is the shared
2// substrate for nx_arousal (simulate) and nx_arousal_fit (fit). Extracted so the two CLIs cannot
3// drift into duplicate twins -- the exact class nx_dupfunc measures (1838 product-only copies).
4// Integer/fixed-point ONLY: time=ms, drive/perfusion/amp=permil, temp=millidegC, vel=cm/s*100.
5//
6// PROVENANCE (auditable, not trusted):
7// Kukkonen 2007 thermography 32.1C -> 33.89C | Chivers 2010 concordance r=.66 M / .26 W
8// Masters&Johnson systolic peaks at PLATEAU onset | Bohlen 1982 contraction chirp +0.1s
9// Loken 2009 C-tactile log-tuned, peak 3 cm/s | Bancroft&Janssen dual control SES/SIS1/SIS2
10// license_tier: ORIGINAL
11import "nx_syscalls.nx"
12const A_MAGIC_999999: i64 = 999999
13const A_MAGIC_1131: i64 = 1131
14const A_MAGIC_2000: i64 = 2000
15const A_MAGIC_6000000: i64 = 6000000
16const A_MAGIC_24000000000: i64 = 24000000000
17const A_MAGIC_120000000000000: i64 = 120000000000000
18const A_MAGIC_1103515245: i64 = 1103515245
19const A_MAGIC_12345: i64 = 12345
20
21const A_DT_MS: i64 = 50
22const A_MAXSTEP: i64 = 60000
23const A_RESMAX: i64 = 1000000000
24
25func a_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
26func a_putn(v: i64) -> i64 {
27 if v == 0 { sys_write(1, "0" as *u8, 1); return 0 }
28 var m: i64 = v
29 if m < 0 { sys_write(1, "-" as *u8, 1); m = 0 - m }
30 let d: *u8 = sys_mmap(32); var k: i64 = 0
31 while m > 0 { d[k] = (48 + (m % 10)) as u8; m = m / 10; k = k + 1 }
32 let o: *u8 = sys_mmap(32); var w: i64 = 0
33 while w < k { o[w] = d[k-1-w]; w = w + 1 }
34 sys_write(1, o, k)
35 sys_munmap(d, 32); sys_munmap(o, 32)
36 return 0
37}
38func a_slen(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } return n }
39func a_atoi(s: *u8) -> i64 {
40 var i: i64 = 0; var v: i64 = 0; var neg: i64 = 0
41 if s[0] == (45 as u8) { neg = 1; i = 1 }
42 while s[i] >= (48 as u8) { if s[i] > (57 as u8) { return 0 - A_MAGIC_999999 } v = v * 10 + ((s[i] as i64) - 48); i = i + 1 }
43 if neg == 1 { return 0 - v }
44 return v
45}
46func a_conf(buf: *u8, n: i64, key: *u8, dflt: i64) -> i64 {
47 if (buf as i64) == 0 { return dflt }
48 let kl: i64 = a_slen(key)
49 if kl == 0 { return dflt }
50 var i: i64 = 0
51 while i + kl < n {
52 var k: i64 = 0; var ok: i64 = 1
53 while k < kl { if buf[i+k] != key[k] { ok = 0; k = kl } k = k + 1 }
54 if ok == 1 { if buf[i+kl] == (61 as u8) { return a_atoi((buf as i64 + i + kl + 1) as *u8) } }
55 i = i + 1
56 }
57 return dflt
58}
59func a_abs(x: i64) -> i64 { if x < 0 { return 0 - x } return x }
60
61// Dual Control (Bancroft & Janssen): excitation minus TWO independent inhibitions.
62// SIS1 (performance-failure threat) scales WITH engagement; SIS2 (consequence threat) is a context load.
63func a_drive(stim: i64, ses: i64, sis1: i64, sis2: i64) -> i64 {
64 let ex: i64 = (ses * stim) / 1000
65 let i1: i64 = (sis1 * stim) / 1000
66 var d: i64 = ex - i1 - sis2
67 if d < 0 { d = 0 }
68 if d > 1000 { d = 1000 }
69 return d
70}
71// One first-order lag step, state in PPM. The ppm scale is LOAD-BEARING: at permil,
72// (1000-cur)*50/185000 floors to 0 every step -> the integrator never leaves zero and every
73// downstream tooth passes against a dead model.
74func a_lag(cur: i64, target: i64, tau_ms: i64) -> i64 {
75 if tau_ms <= 0 { return target }
76 let dl: i64 = target - cur
77 return cur + (dl * A_DT_MS) / tau_ms
78}
79// Integrate ONE compartment to t_ms. Fast neurogenic and slow vascular are this same function at
80// different tau -- do not fork it into twins. drive permil in, permil out, integration at ppm.
81func a_perf(t_ms: i64, drive: i64, tau_ms: i64) -> i64 {
82 var cur: i64 = 0
83 let target: i64 = drive * 1000
84 var steps: i64 = t_ms / A_DT_MS
85 if steps > A_MAXSTEP { steps = A_MAXSTEP }
86 var i: i64 = 0
87 while i < steps { cur = a_lag(cur, target, tau_ms); i = i + 1 }
88 return cur / 1000
89}
90func a_temp(pslow: i64, base_mc: i64, span_mc: i64) -> i64 { return base_mc + (span_mc * pslow) / 1000 }
91func a_hr(drive: i64, base: i64, span: i64) -> i64 { return base + (span * drive) / 1000 }
92func a_rr(drive: i64, base: i64, span: i64) -> i64 { return base + (span * drive) / 1000 }
93// NON-MONOTONIC BY DESIGN: systolic peaks at plateau onset and holds; it does NOT peak at orgasm.
94func a_bp(drive: i64, base: i64, span: i64, plateau: i64) -> i64 {
95 var d: i64 = drive
96 if d > plateau { d = plateau }
97 if plateau <= 0 { return base }
98 return base + (span * d) / plateau
99}
100// Subjective arousal is a SEPARATE channel, weakly coupled in women (Chivers r=.26 vs .66).
101func a_subj(pslow: i64, coup: i64) -> i64 { return (pslow * coup) / 1000 }
102func a_flush(pslow: i64, thresh: i64, propensity: i64) -> i64 {
103 if pslow < thresh { return 0 }
104 if thresh >= 1000 { return 0 }
105 return ((pslow - thresh) * propensity) / (1000 - thresh)
106}
107// --- MALE GENITAL: TUMESCENCE AND RIGIDITY ARE TWO DIFFERENT THINGS -----------------------------
108// Volume fills FIRST (fast, tau~15s -- erection develops in seconds, not the ~664s genital-thermal
109// time-to-peak, which measures skin temperature and not erection). Intracavernosal PRESSURE, and so
110// rigidity, only rises once the tunica is stretched -- so rigidity LAGS tumescence and is non-linear
111// in it. A model with one "erection" scalar is wrong the same way a female twin with one "arousal"
112// scalar is wrong.
113//
114// RigiScan measures RADIAL rigidity as percent of normal maximum at BASE and TIP. Base runs higher
115// and correlates well with axial rigidity; TIP correlates POORLY with buckling pressure. The
116// clinically load-bearing consequence: the conventional 70% rigidity threshold OVERESTIMATES organic
117// erectile dysfunction, while >=60% tip for >=10 min is associated with potency. Those two criteria
118// disagree over a real band of men, and a twin must be able to land inside that band.
119func a_rigidity(tume: i64, thresh: i64) -> i64 {
120 if tume <= thresh { return 0 }
121 if thresh >= 1000 { return 0 }
122 return ((tume - thresh) * 1000) / (1000 - thresh)
123}
124func a_rigidity_tip(rig: i64, tipfac: i64) -> i64 { return (rig * tipfac) / 1000 }
125// Post-orgasmic refractory suppression, permil (1000 = fully suppressed, decaying to 0).
126// ⚠PHENOMENOLOGICAL ONLY. Prolactin as the mechanism is DISPROVEN (Comm Biol 2021: acute
127// manipulation in either direction fails to change refractory duration) and is still repeated
128// everywhere. tau_refr is fitted, not derived -- do not attach a mechanism to it.
129func a_refractory(t_since_ms: i64, tau_refr_ms: i64) -> i64 {
130 if t_since_ms < 0 { return 1000 }
131 if tau_refr_ms <= 0 { return 0 }
132 return a_exp_neg((t_since_ms * 1000) / tau_refr_ms)
133}
134
135// --- MORPHOLOGY: discrete shape change, NOT the vasocongestion gradient -------------------------
136// Nipple erection is MYOTONIC -- areolar smooth-muscle contraction on a ~6s time constant -- not
137// vasocongestive on a ~185s one. Two different mechanisms, two different clocks; collapsing them
138// into one "arousal" gradient is what makes a twin read as a fade rather than a body reacting.
139func a_nipple_erect(t_ms: i64, drive: i64, tau_nip: i64) -> i64 {
140 if drive < 150 { return 0 }
141 return a_perf(t_ms, drive, tau_nip)
142}
143// Areolar engorgement is vasocongestive and therefore SLOW -- it tracks the perfusion envelope.
144func a_areola_engorge(pslow: i64) -> i64 { return pslow }
145// ★APPARENT projection is NON-MONOTONIC (Masters & Johnson): the nipple erects fast and early, then
146// the surrounding areola engorges and MASKS it, so visible projection PEAKS EARLY AND THEN DECLINES
147// while arousal is still climbing. A model that ramps this monotonically to a maximum is wrong in a
148// way that is invisible in numbers and obvious on sight.
149func a_nipple_apparent(erect: i64, engorge: i64, mask: i64) -> i64 {
150 let v: i64 = erect - (mask * engorge) / 1000
151 if v < 0 { return 0 }
152 if v > 1000 { return 1000 }
153 return v
154}
155// Bohlen 1982: interval grows LINEARLY through the train -- a chirp, not a fixed frequency.
156func a_interval(n: i64, c0: i64, slope: i64) -> i64 { return c0 + slope * n }
157func a_amp(n: i64, total: i64) -> i64 {
158 if total <= 0 { return 0 }
159 let half: i64 = total / 2
160 if half <= 0 { return 1000 }
161 if n <= half { return (1000 * n) / half }
162 let rem: i64 = total - half
163 if rem <= 0 { return 0 }
164 return (1000 * (total - n)) / rem
165}
166// C-tactile: inverted-U against LOG velocity, peak 3 cm/s, band 1-10 (Loken 2009).
167func a_ctoct(v: i64, peak: i64) -> i64 {
168 if v <= 0 { return 999 }
169 var r: i64 = 0
170 if v >= peak { r = (v * 100) / peak } else { r = (peak * 100) / v }
171 if r < 141 { return 0 }
172 if r < 200 { return 50 }
173 if r < 283 { return 100 }
174 if r < 400 { return 150 }
175 if r < 566 { return 200 }
176 if r < 800 { return 250 }
177 if r < A_MAGIC_1131 { return 300 }
178 return 350
179}
180func a_ct(v: i64, peak: i64) -> i64 {
181 let d: i64 = a_ctoct(v, peak)
182 var r: i64 = 1000 - (250 * d) / 100
183 if r < 0 { r = 0 }
184 return r
185}
186// --- series parsing + dense-trajectory objective ------------------------------------------------
187// parse a non-negative int at buf[p[0]]; advances p[0]. returns -1 if no digits were consumed.
188func a_pint(buf: *u8, n: i64, p: *i64) -> i64 {
189 var i: i64 = p[0]
190 var v: i64 = 0
191 var got: i64 = 0
192 var d: i64 = 0
193 while d == 0 {
194 if i >= n { d = 1 } else {
195 if buf[i] < (48 as u8) { d = 1 } else {
196 if buf[i] > (57 as u8) { d = 1 } else { v = v * 10 + ((buf[i] as i64) - 48); i = i + 1; got = 1 }
197 }
198 }
199 }
200 p[0] = i
201 if got == 0 { return 0 - 1 }
202 return v
203}
204// advance p[0] to the next digit byte (or to n)
205func a_skipnd(buf: *u8, n: i64, p: *i64) -> i64 {
206 var i: i64 = p[0]
207 var d: i64 = 0
208 while d == 0 {
209 if i >= n { d = 1 } else {
210 if buf[i] >= (48 as u8) { if buf[i] <= (57 as u8) { d = 1 } else { i = i + 1 } } else { i = i + 1 }
211 }
212 }
213 p[0] = i
214 return i
215}
216// --- integer exp(-x), permil in / permil out --------------------------------------------------
217// General primitive (lives here, not in a consumer, so it cannot be re-derived per organ).
218// exp(-x) for x in [0,1): 6-term Taylor. Scaling, with x = r/1000 and output in permil:
219// 1000*exp(-x) = 1000 - r + r^2/2000 - r^3/6e6 + r^4/2.4e10 - r^5/1.2e14
220// At r=1000 this gives 367 against a true 368 (<0.3% error). r^5 peaks near 1e15, well inside i64.
221func a_exp_frac(r: i64) -> i64 {
222 if r <= 0 { return 1000 }
223 let r2: i64 = r * r
224 let r3: i64 = r2 * r
225 let r4: i64 = r3 * r
226 let r5: i64 = r4 * r
227 var v: i64 = 1000 - r + (r2 / A_MAGIC_2000) - (r3 / A_MAGIC_6000000) + (r4 / A_MAGIC_24000000000) - (r5 / A_MAGIC_120000000000000)
228 if v < 0 { v = 0 }
229 if v > 1000 { v = 1000 }
230 return v
231}
232// exp(-a/1000) for absorbance a in permil. Splits a into whole units and a remainder:
233// exp(-a) = exp(-1)^q * exp(-rem), with exp(-1) = 368 permil.
234// ★THE POINT: an exponential NEVER REACHES ZERO. Linear subtraction clamps a channel to 0 and
235// destroys it; this preserves signal at every physically reachable absorbance.
236func a_exp_neg(a: i64) -> i64 {
237 if a <= 0 { return 1000 }
238 let q: i64 = a / 1000
239 let rem: i64 = a % 1000
240 var v: i64 = a_exp_frac(rem)
241 var i: i64 = 0
242 var stop: i64 = 0
243 while i < q {
244 if stop == 0 {
245 v = (v * 368) / 1000
246 if v <= 0 { v = 1; stop = 1 }
247 }
248 i = i + 1
249 }
250 return v
251}
252// Deterministic LCG. Sensor noise must be REPRODUCIBLE or a noise-robustness tooth cannot be a tooth:
253// same seed must give the same series on every run, on every host.
254func a_lcg(s: i64) -> i64 { return ((s * A_MAGIC_1103515245) + A_MAGIC_12345) & 0x7FFFFFFF }
255// count "t obs" pairs in a series -- needed to scale the confidence slack per-sample
256func a_series_count(buf: *u8, n: i64) -> i64 {
257 if (buf as i64) == 0 { return 0 }
258 let p: *i64 = sys_mmap(16) as *i64
259 p[0] = 0
260 var pts: i64 = 0
261 var done: i64 = 0
262 while done == 0 {
263 a_skipnd(buf, n, p)
264 if p[0] >= n { done = 1 } else {
265 let t: i64 = a_pint(buf, n, p)
266 a_skipnd(buf, n, p)
267 let o: i64 = a_pint(buf, n, p)
268 if t < 0 { done = 1 } else { if o < 0 { done = 1 } else { pts = pts + 1 } }
269 }
270 }
271 sys_munmap(p as *u8, 16)
272 return pts
273}
274// Dense-trajectory objective over a "t_ms obs_permil" series (one pair per line).
275// SINGLE forward integration that samples as it passes each timestamp: O(max_t/dt), NOT O(sum t).
276// Calling a_perf per sample would re-integrate from zero every point and make a 2D sweep intractable.
277// REFUSES a non-ascending series: an out-of-order trajectory is malformed input, not data to fit.
278func a_series_residual(buf: *u8, n: i64, tau: i64, drive: i64) -> i64 {
279 if (buf as i64) == 0 { return A_RESMAX }
280 if tau <= 0 { return A_RESMAX }
281 if drive <= 0 { return A_RESMAX }
282 let p: *i64 = sys_mmap(16) as *i64
283 p[0] = 0
284 let target: i64 = drive * 1000
285 var cur: i64 = 0
286 var cur_t: i64 = 0
287 var total: i64 = 0
288 var pts: i64 = 0
289 var bad: i64 = 0
290 var done: i64 = 0
291 while done == 0 {
292 a_skipnd(buf, n, p)
293 if p[0] >= n { done = 1 } else {
294 let t: i64 = a_pint(buf, n, p)
295 a_skipnd(buf, n, p)
296 let o: i64 = a_pint(buf, n, p)
297 if t < 0 { done = 1 } else {
298 if o < 0 { done = 1 } else {
299 if t < cur_t { bad = 1; done = 1 } else {
300 if o > 1000 { bad = 1; done = 1 } else {
301 var steps: i64 = (t - cur_t) / A_DT_MS
302 if steps > A_MAXSTEP { steps = A_MAXSTEP }
303 var k: i64 = 0
304 while k < steps { cur = a_lag(cur, target, tau); k = k + 1 }
305 cur_t = cur_t + steps * A_DT_MS
306 total = total + a_abs((cur / 1000) - o)
307 pts = pts + 1
308 }
309 }
310 }
311 }
312 }
313 }
314 sys_munmap(p as *u8, 16)
315 if bad == 1 { return A_RESMAX }
316 if pts == 0 { return A_RESMAX }
317 if total < 0 { return A_RESMAX }
318 return total
319}
320
321// Two-point fit objective. A residual is a DISTANCE: it can never be negative.
322// Returns A_RESMAX as an explicit refusal sentinel on a corrupt/degenerate input.
323func a_residual(tau: i64, drive: i64, t1: i64, o1: i64, t2: i64, o2: i64) -> i64 {
324 if tau <= 0 { return A_RESMAX }
325 let m1: i64 = a_perf(t1, drive, tau)
326 let m2: i64 = a_perf(t2, drive, tau)
327 let r: i64 = a_abs(m1 - o1) + a_abs(m2 - o2)
328 if r < 0 { return A_RESMAX }
329 return r
330}