code wiki / (root) / nx_ansur_regress.nx
`nx_ansur_regress.nx` source

↩ module page · 348 lines · 14802 B
1// nx_ansur_regress.nx -- DOES HUMAN PROPORTION CHANGE WITH STATURE? Regress a landmark ratio on
2// stature over the ANSUR II sample and let the SLOPE answer.
3//
4// WHY. nx_skelgen stores every segment as PER-MILLE OF STATURE, and its own docstring states the
5// consequence as a feature: "so the same rules generate a child, an adult or a giant". Measured,
6// femur/stature is 0.2411 at 1750 mm and 0.2407/0.2411 at 1400/1900 -- CONSTANT BY CONSTRUCTION.
7// It is not emergent; it is the storage format. A rig that only scales is a doll at every size.
8// The fix is mechanical because ANSUR measures nx_skelgen's OWN landmarks:
9//   acromialheight/stature      <-> BP_ACROMION = 818 per-mille
10//   trochanterionheight/stature <-> BP_TROCH    = 530 per-mille
11// Regress (landmark_height / stature) ON stature. A NON-ZERO SLOPE IS THE MISSING PROPORTION
12// CHANGE, and the slope itself is the correction term.
13//
14//   nx_ansur_regress <ansur.csv> <landmark_column_name>
15//
16// COLUMNS ARE RESOLVED BY NAME FROM THE HEADER ROW AT RUNTIME, never by a hardcoded index: ANSUR's
17// columns are alphabetical and the male and female files must both parse, so an index literal would
18// be a magic number that silently reads the wrong measure if either file differs.
19//
20// THE THREE TRAPS, QUOTED FROM THE PUBLISHER, ALL OF WHICH BITE THIS COMPUTATION:
21//  1. Units are millimetres EXCEPT MASS, WHICH IS HECTOGRAMS DESPITE ITS COLUMN BEING NAMED
22//     weightkg -- a 10x error waiting. This organ touches no mass column, and refuses any column
23//     whose name contains "weight" so the trap cannot be walked into later.
24//  2. The trailing Heightin / Weightlbs columns are SELF-REPORTED and DISAGREE with the official
25//     measures. This organ uses `stature`, never Heightin, and refuses those names by name.
26//  3. The sample is military including reservists and the publisher states it "is still not an
27//     approximation of the US Civilian population". THAT BOUND TRAVELS WITH EVERY NUMBER THIS
28//     ORGAN PRINTS, so it is printed WITH them rather than left in a doc nobody reads.
29//
30// ARITHMETIC. Integer-only, and mean-centred ON PURPOSE: the raw cross-products (n*Sxy - Sx*Sy)
31// reach ~5.7e12 and multiplying by a precision scale would approach the i64 ceiling. Centring first
32// keeps every term small (dx ~ +/-100, dy ~ +/-30), so the scaled slope cannot overflow. Two passes,
33// no floats -- NishiLang has no float type and this needs none.
34// license_tier: ORIGINAL  expect_exit: 0
35import "nx_syscalls.nx"
36
37const AR_E6: i64 = 1000000
38const AR_PERMIL: i64 = 1000
39const AR_COMMA: i64 = 44
40const AR_NL: i64 = 10
41const AR_CR: i64 = 13
42const AR_D0: i64 = 48
43const AR_D9: i64 = 57
44const AR_DOT: i64 = 46
45const AR_MINUS: i64 = 45
46// plausibility band for a human stature in mm. NOT a filter on the data -- a COLUMN-SHIFT DETECTOR:
47// if the parsed stature column stops looking like statures, the header mapping is wrong and every
48// number below it would be confident nonsense.
49const AR_STAT_LO: i64 = 1000
50const AR_STAT_HI: i64 = 2200
51
52func ar_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
53func ar_num(v: i64) -> i64 {
  let b: *u8 = sys_mmap(32) as *u8
  var x: i64 = v
  var ng: i64 = 0
  if x < 0 { ng = 1; x = 0 - x }
  var i: i64 = 31
  if x == 0 { b[i] = AR_D0 as u8; i = i - 1 }
  while x > 0 { b[i] = (AR_D0 + x % 10) as u8; x = x / 10; i = i - 1 }
  if ng == 1 { b[i] = AR_MINUS as u8; i = i - 1 }
  sys_write(1, ((b as i64) + i + 1) as *u8, 31 - i)
  return 0
64}
65func ar_kv(k: *u8, v: i64) -> i64 { ar_puts(k); ar_num(v); ar_puts("\n" as *u8); return 0 }
66func ar_refuse(r: *u8) -> i64 { ar_puts("ANSUR-REGRESS REFUSED: " as *u8); ar_puts(r); ar_puts("\n" as *u8); return 0 }
67func ar_slen(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } return n }
68
69// does the buffer region [a,b) equal the C string s ?
70func ar_eq(buf: *u8, a: i64, b: i64, s: *u8) -> i64 {
  let n: i64 = ar_slen(s)
  if b - a != n { return 0 }
  var i: i64 = 0
  while i < n { if buf[a+i] != s[i] { return 0 } i = i + 1 }
  return 1
76}
77// does the buffer region contain the C string s as a substring (for the refused-name guards)
78func ar_contains(buf: *u8, a: i64, b: i64, s: *u8) -> i64 {
  let n: i64 = ar_slen(s)
  if n == 0 { return 0 }
  var i: i64 = a
  while i + n <= b {
      var k: i64 = 0
      var same: i64 = 1
      while k < n { if buf[i+k] != s[k] { same = 0 } k = k + 1 }
      if same == 1 { return 1 }
      i = i + 1
  }
  return 0
90}
91// parse an unsigned integer from [a,b); a decimal point truncates (ANSUR heights are whole mm)
92func ar_parse(buf: *u8, a: i64, b: i64) -> i64 {
  var v: i64 = 0
  var got: i64 = 0
  var i: i64 = a
  var go: i64 = 1
  while go == 1 {
      if i >= b { go = 0 }
      if go == 1 {
          let c: i64 = buf[i] as i64
          if c == AR_DOT { go = 0 } else {
              var isd: i64 = 0
              if c >= AR_D0 { if c <= AR_D9 { isd = 1 } }
              if isd == 1 { v = v*10 + (c - AR_D0); got = 1; i = i + 1 } else { go = 0 }
          }
      }
  }
  if got == 0 { return 0 - 1 }
  return v
110}
111// end of the current line starting at p
112func ar_eol(buf: *u8, n: i64, p: i64) -> i64 {
  var i: i64 = p
  var go: i64 = 1
  while go == 1 {
      if i >= n { go = 0 }
      if go == 1 { if buf[i] == (AR_NL as u8) { go = 0 } else { i = i + 1 } }
  }
  return i
120}
121// index of field `want` within the line [p,e), or -1
122func ar_field_bounds(buf: *u8, p: i64, e: i64, want: i64, sa: *i64, sb: *i64) -> i64 {
  var f: i64 = 0
  var a: i64 = p
  var i: i64 = p
  var go: i64 = 1
  while go == 1 {
      if i >= e {
          if f == want { sa[0] = a; sb[0] = i; return 1 }
          go = 0
      }
      if go == 1 {
          if buf[i] == (AR_COMMA as u8) {
              if f == want { sa[0] = a; sb[0] = i; return 1 }
              f = f + 1
              a = i + 1
          }
          i = i + 1
      }
  }
  return 0
142}
143// resolve a column index BY NAME from the header line
144func ar_col_by_name(buf: *u8, p: i64, e: i64, name: *u8) -> i64 {
  var f: i64 = 0
  var a: i64 = p
  var i: i64 = p
  var go: i64 = 1
  while go == 1 {
      if i >= e {
          var b2: i64 = i
          if b2 > a { if buf[b2-1] == (AR_CR as u8) { b2 = b2 - 1 } }
          if ar_eq(buf, a, b2, name) == 1 { return f }
          go = 0
      }
      if go == 1 {
          if buf[i] == (AR_COMMA as u8) {
              if ar_eq(buf, a, i, name) == 1 { return f }
              f = f + 1
              a = i + 1
          }
          i = i + 1
      }
  }
  return 0 - 1
166}
167
168// MEAN VERB (2026-09-02, the faceanat scale derivation): `nx_ansur_regress <csv> <column> mean` prints the
169// column's mean, n, min and max in the file's own units (mm) and exits -- a POPULATION reading for the
170// aesthetic canon, never a regression. Same header-resolved column, same trap guards, and the publisher's
171// population bound is printed WITH the number rather than left in a doc nobody reads. Rows whose value
172// does not parse are COUNTED and reported, never silently dropped into the denominator.
173func ar_mean(buf: *u8, n: i64, h_end: i64, c_land: i64, lname: *u8) -> i64 {
  let sa: *i64 = sys_mmap(16) as *i64
  let sb: *i64 = sys_mmap(16) as *i64
  var rows: i64 = 0
  var bad: i64 = 0
  var sum: i64 = 0
  var lo: i64 = 0
  var hi: i64 = 0
  var p: i64 = h_end + 1
  while p < n {
      let e: i64 = ar_eol(buf, n, p)
      if e > p {
          var lv: i64 = 0 - 1
          if ar_field_bounds(buf, p, e, c_land, sa, sb) == 1 { lv = ar_parse(buf, sa[0], sb[0]) }
          if lv > 0 {
              if rows == 0 { lo = lv; hi = lv }
              if lv < lo { lo = lv }
              if lv > hi { hi = lv }
              sum = sum + lv
              rows = rows + 1
          } else { bad = bad + 1 }
      }
      p = e + 1
  }
  ar_puts("verb=mean column=" as *u8); ar_puts(lname); ar_puts("\n" as *u8)
  ar_kv("rows_used=" as *u8, rows)
  ar_kv("rows_unparsable=" as *u8, bad)
  if rows < 1 { ar_refuse("no usable rows for that column -- no mean to publish" as *u8); return 6 }
  ar_kv("mean_mm=" as *u8, sum / rows)
  ar_kv("mean_mm_x10=" as *u8, (sum * 10) / rows)
  ar_kv("min_mm=" as *u8, lo)
  ar_kv("max_mm=" as *u8, hi)
  ar_puts("POPULATION BOUND, TRAVELS WITH EVERY NUMBER ABOVE: ANSUR II is a military sample including reservists; the publisher states it is still not an approximation of the US Civilian population.\n" as *u8)
  return 0
207}
208
209func main(argc: i64, argv: *i64) -> i64 {
  if argc < 3 {
      ar_puts("usage: nx_ansur_regress <ansur.csv> <landmark_column_name>\n" as *u8)
      return 3
  }
  let path: *u8 = argv[1] as *u8
  let lname: *u8 = argv[2] as *u8
216
  // TRAP GUARDS -- refuse the publisher's known-bad columns BY NAME so the trap cannot be walked
  // into by a future caller rather than merely documented.
  let ln: i64 = ar_slen(lname)
  if ar_contains(lname, 0, ln, "weight" as *u8) == 1 {
      ar_refuse("mass columns are HECTOGRAMS despite being named weightkg -- this organ does not touch mass" as *u8)
      return 4
  }
  if ar_contains(lname, 0, ln, "Heightin" as *u8) == 1 {
      ar_refuse("Heightin is SELF-REPORTED and disagrees with the official measures -- use stature" as *u8)
      return 4
  }
228
  let lnp: *i64 = sys_mmap(16) as *i64
  lnp[0] = 0
  let buf: *u8 = sys_read_file(path, lnp)
  if (buf as i64) == 0 { ar_refuse("cannot read the ANSUR csv" as *u8); return 4 }
  let n: i64 = lnp[0]
  ar_kv("bytes=" as *u8, n)
235
  let h_end: i64 = ar_eol(buf, n, 0)
  let c_stat: i64 = ar_col_by_name(buf, 0, h_end, "stature" as *u8)
  let c_land: i64 = ar_col_by_name(buf, 0, h_end, lname)
  ar_puts("landmark_column="); ar_puts(lname); ar_puts("\n" as *u8)
  ar_kv("col_stature=" as *u8, c_stat)
  ar_kv("col_landmark=" as *u8, c_land)
  if c_land < 0 { ar_refuse("named landmark column not found in the header" as *u8); return 5 }
  // the mean verb needs only the landmark column; it runs BEFORE the stature requirement below
  if argc >= 4 { let vb: *u8 = argv[3] as *u8; if ar_eq(vb, 0, ar_slen(vb), "mean" as *u8) == 1 { return ar_mean(buf, n, h_end, c_land, lname) } }
  if c_stat < 0 { ar_refuse("no column named stature in the header" as *u8); return 5 }
246
  // ---- pass 1: collect, and DETECT COLUMN SHIFT rather than trusting the mapping ----
  let cap: i64 = 8192
  let xs: *i64 = sys_mmap(cap*8) as *i64
  let ys: *i64 = sys_mmap(cap*8) as *i64
  let sa: *i64 = sys_mmap(16) as *i64
  let sb: *i64 = sys_mmap(16) as *i64
  var rows: i64 = 0
  var bad: i64 = 0
  var implausible: i64 = 0
  var p: i64 = h_end + 1
  var over: i64 = 0
  while p < n {
      let e: i64 = ar_eol(buf, n, p)
      if e > p {
          var okrow: i64 = 0
          var sv: i64 = 0 - 1
          var lv: i64 = 0 - 1
          if ar_field_bounds(buf, p, e, c_stat, sa, sb) == 1 { sv = ar_parse(buf, sa[0], sb[0]) }
          if ar_field_bounds(buf, p, e, c_land, sa, sb) == 1 { lv = ar_parse(buf, sa[0], sb[0]) }
          if sv > 0 { if lv > 0 { okrow = 1 } }
          if okrow == 1 {
              var plaus: i64 = 0
              if sv >= AR_STAT_LO { if sv <= AR_STAT_HI { plaus = 1 } }
              if plaus == 0 { implausible = implausible + 1 } else {
                  if rows < cap {
                      xs[rows] = sv
                      ys[rows] = (lv * AR_PERMIL) / sv
                      rows = rows + 1
                  } else { over = over + 1 }
              }
          } else { bad = bad + 1 }
      }
      p = e + 1
  }
  ar_kv("rows_used=" as *u8, rows)
  ar_kv("rows_unparsable=" as *u8, bad)
  ar_kv("rows_stature_outside_plausible_band=" as *u8, implausible)
  ar_kv("rows_dropped_over_capacity=" as *u8, over)
  if over > 0 {
      ar_refuse("row capacity exceeded -- refusing to publish a regression over a SILENTLY TRUNCATED sample" as *u8)
      return 6
  }
  if implausible > 0 {
      ar_refuse("stature column carries values outside 1000-2200mm -- the header mapping is wrong and every derived number would be confident nonsense" as *u8)
      return 6
  }
  if rows < 2 { ar_refuse("fewer than two usable rows -- a slope over one point is undefined" as *u8); return 6 }
294
  // ---- pass 2: mean-centred least squares (integer, overflow-safe by construction) ----
  var sx: i64 = 0
  var sy: i64 = 0
  var i: i64 = 0
  while i < rows { sx = sx + xs[i]; sy = sy + ys[i]; i = i + 1 }
  let mx: i64 = sx / rows
  let my: i64 = sy / rows
  var sxy: i64 = 0
  var sxx: i64 = 0
  var lo: i64 = xs[0]
  var hi: i64 = xs[0]
  i = 0
  while i < rows {
      let dx: i64 = xs[i] - mx
      let dy: i64 = ys[i] - my
      sxy = sxy + dx*dy
      sxx = sxx + dx*dx
      if xs[i] < lo { lo = xs[i] }
      if xs[i] > hi { hi = xs[i] }
      i = i + 1
  }
  ar_kv("mean_stature_mm=" as *u8, mx)
  ar_kv("mean_ratio_permil=" as *u8, my)
  ar_kv("stature_min_mm=" as *u8, lo)
  ar_kv("stature_max_mm=" as *u8, hi)
  if sxx <= 0 { ar_refuse("zero variance in stature -- slope undefined" as *u8); return 6 }
321
  // slope in MILLIONTHS of a per-mille per millimetre
  let slope_e6: i64 = (sxy * AR_E6) / sxx
  ar_kv("slope_permil_per_mm_e6=" as *u8, slope_e6)
325
  // the number that matters: how much the ratio moves across the OBSERVED stature range
  let span: i64 = hi - lo
  let d_e6: i64 = slope_e6 * span
  ar_kv("stature_span_mm=" as *u8, span)
  ar_kv("delta_ratio_permil_e6_over_span=" as *u8, d_e6)
  ar_kv("delta_ratio_permil_over_span=" as *u8, d_e6 / AR_E6)
332
  // predicted ratio at each end of the observed range -- the DOLL TEST in one line
  let r_lo_e6: i64 = my * AR_E6 + slope_e6 * (lo - mx)
  let r_hi_e6: i64 = my * AR_E6 + slope_e6 * (hi - mx)
  ar_kv("ratio_at_stature_min_permil_e6=" as *u8, r_lo_e6)
  ar_kv("ratio_at_stature_max_permil_e6=" as *u8, r_hi_e6)
338
  var verdict: i64 = 0
  if slope_e6 != 0 { verdict = 1 }
  if verdict == 1 {
      ar_puts("VERDICT=PROPORTION-CHANGES-WITH-STATURE (non-zero slope: per-mille-of-stature storage CANNOT express this)\n" as *u8)
  } else {
      ar_puts("VERDICT=PROPORTION-CONSTANT (zero slope: pure scaling would be anthropometrically correct)\n" as *u8)
  }
  ar_puts("POPULATION BOUND, TRAVELS WITH EVERY NUMBER ABOVE: ANSUR II is a military sample including reservists; the publisher states it is still not an approximation of the US Civilian population.\n" as *u8)
  return 0
348}
nx_ansur_regress.nx source

`nx_ansur_regress.nx` source
