1use std::collections::{BTreeMap, HashMap};
77
78use serde::{Deserialize, Serialize};
79
80use super::summary::{Apply, Check as Logged, Gate, Label, Summary};
81use super::{ARMS, Route, Unit};
82use crate::conscience::{self, Challenge};
83use crate::memory::Rng;
84
85pub const ALPHA: f64 = 0.10;
90pub const DELTA: f64 = 0.05;
93pub const DRAWS: usize = conscience::P_DRAWS;
96pub const SOLVE_MAX: usize = 16;
101
102pub const PRIOR_SD: f64 = 2.5;
106pub const ASK_PRIOR: (f64, f64) = (0.0, 0.35);
109pub const NOT_PRIOR: (f64, f64) = (-0.3, 0.25);
112pub const DILUTION_PRIOR: (f64, f64) = (0.0, 0.1);
116pub const ARD_START: f64 = 1.0;
119pub const LEVEL_START: f64 = 0.5;
121pub const SD_FLOOR: f64 = 0.01;
124pub const HALF_LIVES: [Option<u32>; 5] = [None, Some(64), Some(32), Some(16), Some(8)];
127pub const FIT_SESSIONS: usize = 256;
130pub const MAX_IPW: f64 = 10.0;
133pub const PSI_FLOOR: f64 = 1e-3;
137pub const POOL_STEPS: [i32; 5] = [-2, -1, 0, 1, 2];
140pub const LATENCY_MIN: usize = 20;
143
144const TOL: f64 = 1e-8;
146const MAX_NEWTON: usize = 50;
148const MAX_HALVINGS: usize = 30;
151const MAX_EB: usize = 50;
153const EB_TOL: f64 = 1e-4;
155const PREQ_SESSIONS: usize = 16;
158const CS_SESSIONS: f64 = 30.0;
161const CS_MIN: usize = 10;
164const PRICE_STEPS: usize = 60;
166
167const B0: usize = 0;
169const G2: usize = 1;
170const G3: usize = 2;
171const RANK: usize = 3;
172const STOP: usize = 4;
173const THREE: usize = 5;
174const ASK: usize = 6;
175const NOT: usize = 7;
176const ASK_M: usize = 8;
177const FREE_M: usize = 9;
178const ASK_STOP: usize = 10;
179const ASK_G3: usize = 11;
180const NOT_G: usize = 12;
181const POOL: usize = 13;
182const THREE_STOP: usize = 14;
183pub const P: usize = 15;
186pub const NAMES: [&str; P] = [
188 "intercept",
189 "grade2",
190 "grade3",
191 "rank",
192 "stop",
193 "three",
194 "ask",
195 "not",
196 "ask_dilution",
197 "free_dilution",
198 "ask_stop",
199 "ask_grade3",
200 "not_grade",
201 "pool",
202 "three_stop",
203];
204const SESSION: usize = 0;
206const LESSON: usize = 1;
207const JUDGE: usize = 2;
208const FAMILIES: usize = 3;
209
210#[derive(Debug, Clone, Copy)]
212struct Coef {
213 mean: f64,
214 sd: f64,
215 ard: bool,
216 nonpos: bool,
217}
218
219fn coefs() -> [Coef; P] {
222 let fixed = |(mean, sd): (f64, f64), nonpos| Coef { mean, sd, ard: false, nonpos };
223 let ard = |nonpos| Coef { mean: 0.0, sd: ARD_START, ard: true, nonpos };
224 let mut c = [fixed((0.0, PRIOR_SD), false); P];
225 c[RANK] = ard(false);
226 c[ASK] = fixed(ASK_PRIOR, false);
227 c[NOT] = fixed(NOT_PRIOR, false);
228 c[ASK_M] = fixed(DILUTION_PRIOR, true);
229 c[FREE_M] = fixed(DILUTION_PRIOR, true);
230 for i in [ASK_STOP, ASK_G3, NOT_G] {
231 c[i] = ard(false);
232 }
233 c[POOL] = ard(true);
234 c
235}
236
237const PRUNED_KEEPS: [usize; 7] = [B0, G2, G3, ASK, NOT, ASK_M, FREE_M];
240
241fn pruned() -> [Coef; P] {
244 let mut c = coefs();
245 for (k, x) in c.iter_mut().enumerate() {
246 if !PRUNED_KEEPS.contains(&k) {
247 *x = Coef { mean: 0.0, sd: SD_FLOOR, ard: false, nonpos: false };
248 }
249 }
250 c
251}
252
253#[derive(Debug, Clone, Copy, PartialEq)]
255pub struct Feat {
256 pub grade: u8,
257 pub rank: usize,
258 pub stop: bool,
260 pub three: bool,
262 pub pool: f64,
264}
265
266impl Feat {
267 pub fn of(units: &[Unit], stop: bool, tokens: u32) -> Vec<Feat> {
269 let three = units.iter().any(|u| u.grade == 3);
270 let pool = pool_x(tokens);
271 units.iter().map(|u| Feat { grade: u.grade, rank: u.rank, stop, three, pool }).collect()
272 }
273}
274
275fn pool_x(tokens: u32) -> f64 {
277 if tokens == 0 { 0.0 } else { (f64::from(tokens) / crate::whole::TOKENS as f64).ln() }
278}
279
280fn design(f: &Feat, route: Route, m: usize) -> [f64; P] {
282 let mut x = [0.0; P];
283 let g = f.grade;
284 x[B0] = 1.0;
285 x[G2] = f64::from(u8::from(g == 2));
286 x[G3] = f64::from(u8::from(g == 3));
287 x[RANK] = (1.0 + f.rank as f64).ln();
288 x[STOP] = f64::from(u8::from(f.stop));
289 x[THREE] = f64::from(u8::from(f.three));
290 x[THREE_STOP] = f64::from(u8::from(f.stop && f.three));
291 x[POOL] = f.pool;
292 let ln_m = (m.max(1) as f64).ln();
293 match route {
294 Route::Ask => {
295 x[ASK] = 1.0;
296 x[ASK_M] = ln_m;
297 x[ASK_STOP] = x[STOP];
298 x[ASK_G3] = x[G3];
299 }
300 Route::Free => x[FREE_M] = ln_m,
301 Route::Placebo | Route::Not => {
302 x[NOT] = 1.0;
303 x[NOT_G] = f64::from(g) - 2.0;
304 }
305 }
306 x
307}
308
309fn dot<const K: usize>(a: &[f64; K], b: &[f64; K]) -> f64 {
310 a.iter().zip(b).map(|(x, y)| x * y).sum()
311}
312
313fn sigmoid(x: f64) -> f64 {
314 if x >= 0.0 {
315 1.0 / (1.0 + (-x).exp())
316 } else {
317 let e = x.exp();
318 e / (1.0 + e)
319 }
320}
321
322fn log1pexp(x: f64) -> f64 {
324 if x > 0.0 { x + (-x).exp().ln_1p() } else { x.exp().ln_1p() }
325}
326
327fn cholesky(a: &mut [f64], n: usize) -> bool {
330 for j in 0..n {
331 let d = a[j * n + j] - (0..j).map(|k| a[j * n + k] * a[j * n + k]).sum::<f64>();
332 if d.is_nan() || d <= 0.0 || !d.is_finite() {
333 return false;
334 }
335 let d = d.sqrt();
336 a[j * n + j] = d;
337 for i in j + 1..n {
338 let s = a[i * n + j] - (0..j).map(|k| a[i * n + k] * a[j * n + k]).sum::<f64>();
339 a[i * n + j] = s / d;
340 }
341 for i in 0..j {
342 a[i * n + j] = 0.0;
343 }
344 }
345 true
346}
347
348fn chol_solve(l: &[f64], n: usize, b: &mut [f64]) {
350 for i in 0..n {
351 let s = b[i] - (0..i).map(|k| l[i * n + k] * b[k]).sum::<f64>();
352 b[i] = s / l[i * n + i];
353 }
354 for i in (0..n).rev() {
355 let s = b[i] - (i + 1..n).map(|k| l[k * n + i] * b[k]).sum::<f64>();
356 b[i] = s / l[i * n + i];
357 }
358}
359
360fn chol_inverse(l: &[f64], n: usize) -> Vec<f64> {
362 let mut inv = vec![0.0; n * n];
363 let mut e = vec![0.0; n];
364 for j in 0..n {
365 e.iter_mut().for_each(|x| *x = 0.0);
366 e[j] = 1.0;
367 chol_solve(l, n, &mut e);
368 for (i, x) in e.iter().enumerate() {
369 inv[i * n + j] = *x;
370 }
371 }
372 inv
373}
374
375fn normal(rng: &mut Rng) -> f64 {
377 let (u, v) = (rng.unit(), rng.unit());
378 (-2.0 * u.ln()).sqrt() * (std::f64::consts::TAU * v).cos()
379}
380
381fn gamma(shape: f64, rng: &mut Rng) -> f64 {
383 if shape < 1.0 {
384 return gamma(shape + 1.0, rng) * rng.unit().powf(1.0 / shape);
385 }
386 let d = shape - 1.0 / 3.0;
387 let c = 1.0 / (9.0 * d).sqrt();
388 for _ in 0..64 {
391 let z = normal(rng);
392 let v = (1.0 + c * z).powi(3);
393 if v <= 0.0 {
394 continue;
395 }
396 if rng.unit().ln() < 0.5 * z * z + d - d * v + d * v.ln() {
397 return d * v;
398 }
399 }
400 d
401}
402
403fn beta_draw(a: f64, b: f64, rng: &mut Rng) -> f64 {
405 let (x, y) = (gamma(a, rng), gamma(b, rng));
406 if x + y > 0.0 { x / (x + y) } else { a / (a + b) }
407}
408
409fn hash(parts: &[&[u8]]) -> u64 {
411 let mut h = blake3::Hasher::new();
412 h.update(b"borg-select-policy-v1\0");
413 for p in parts {
414 h.update(p);
415 h.update(b"\0");
416 }
417 let mut first = [0u8; 8];
418 first.copy_from_slice(&h.finalize().as_bytes()[..8]);
419 u64::from_le_bytes(first)
420}
421
422fn mix(seed: u64, i: u64) -> u64 {
424 Rng::new(seed ^ i.wrapping_add(1).wrapping_mul(0x9e37_79b9_7f4a_7c15)).next_u64()
425}
426
427#[derive(Debug, Clone)]
429struct Obs<const K: usize, const F: usize> {
430 x: [f64; K],
431 y: f64,
432 w: f64,
433 lv: [Option<usize>; F],
434}
435
436#[derive(Debug, Clone)]
438struct Est<const K: usize, const F: usize> {
439 beta: [f64; K],
440 u: [Vec<f64>; F],
441 sd: [f64; K],
443 spread: [f64; F],
444 cov: [[f64; K]; K],
446 var_u: [Vec<f64>; F],
447 converged: bool,
448}
449
450impl<const K: usize, const F: usize> Est<K, F> {
451 fn eta(&self, o: &Obs<K, F>) -> f64 {
452 dot(&self.beta, &o.x) + (0..F).filter_map(|f| o.lv[f].map(|l| self.u[f][l])).sum::<f64>()
453 }
454}
455
456fn map<const K: usize, const F: usize>(
461 obs: &[Obs<K, F>],
462 levels: [usize; F],
463 spec: &[Coef; K],
464 sd: &[f64; K],
465 spread: &[f64; F],
466 warm: Option<&Est<K, F>>,
467) -> Option<Est<K, F>> {
468 let elim = (0..F).filter(|&f| levels[f] > 0).max_by_key(|&f| (levels[f], std::cmp::Reverse(f)));
469 let mut off = [0usize; F];
470 let mut dim = K;
471 for f in 0..F {
472 if Some(f) != elim {
473 off[f] = dim;
474 dim += levels[f];
475 }
476 }
477 let members: Vec<Vec<usize>> = match elim {
478 Some(e) => {
479 let mut m = vec![Vec::new(); levels[e]];
480 for (i, o) in obs.iter().enumerate() {
481 if let Some(l) = o.lv[e] {
482 m[l].push(i);
483 }
484 }
485 m
486 }
487 None => Vec::new(),
488 };
489 let mut beta: [f64; K] = std::array::from_fn(|k| warm.map_or(spec[k].mean, |w| w.beta[k]));
490 for k in 0..K {
491 if spec[k].nonpos {
492 beta[k] = beta[k].min(0.0);
493 }
494 }
495 let mut u: [Vec<f64>; F] = std::array::from_fn(|f| match warm {
496 Some(w) if w.u[f].len() == levels[f] => w.u[f].clone(),
497 _ => vec![0.0; levels[f]],
498 });
499 let prec: [f64; K] = std::array::from_fn(|k| 1.0 / (sd[k] * sd[k]));
500 let uprec: [f64; F] = std::array::from_fn(|f| 1.0 / (spread[f] * spread[f]));
501 let objective = |beta: &[f64; K], u: &[Vec<f64>; F]| -> f64 {
502 let mut ll = 0.0;
503 for o in obs {
504 let eta = dot(beta, &o.x) + (0..F).filter_map(|f| o.lv[f].map(|l| u[f][l])).sum::<f64>();
505 ll += o.w * (o.y * eta - log1pexp(eta));
506 }
507 let pb: f64 = (0..K).map(|k| (beta[k] - spec[k].mean).powi(2) * prec[k]).sum();
508 let pu: f64 = (0..F).map(|f| u[f].iter().map(|x| x * x).sum::<f64>() * uprec[f]).sum();
509 ll - 0.5 * (pb + pu)
510 };
511 let mut converged = false;
512 let mut s = vec![0.0; dim * dim];
513 let mut d_e: Vec<f64> = Vec::new();
514 let mut b_e: Vec<Vec<(usize, f64)>> = Vec::new();
515 let mut current = objective(&beta, &u);
516 for _ in 0..MAX_NEWTON {
517 let mut g = vec![0.0; dim];
519 s.iter_mut().for_each(|x| *x = 0.0);
520 let mut r = vec![0.0; obs.len()];
521 let mut wt = vec![0.0; obs.len()];
522 let mut dense: Vec<(usize, f64)> = Vec::with_capacity(K + F);
523 for (i, o) in obs.iter().enumerate() {
524 let eta = dot(&beta, &o.x) + (0..F).filter_map(|f| o.lv[f].map(|l| u[f][l])).sum::<f64>();
525 let mu = sigmoid(eta);
526 r[i] = o.w * (o.y - mu);
527 wt[i] = o.w * mu * (1.0 - mu);
528 dense.clear();
529 dense.extend((0..K).filter(|&k| o.x[k] != 0.0).map(|k| (k, o.x[k])));
530 dense.extend((0..F).filter(|&f| Some(f) != elim).filter_map(|f| o.lv[f].map(|l| (off[f] + l, 1.0))));
531 for &(a, xa) in &dense {
532 g[a] += r[i] * xa;
533 for &(b, xb) in &dense {
534 s[a * dim + b] += wt[i] * xa * xb;
535 }
536 }
537 }
538 for k in 0..K {
539 g[k] -= (beta[k] - spec[k].mean) * prec[k];
540 s[k * dim + k] += prec[k];
541 }
542 for f in (0..F).filter(|&f| Some(f) != elim) {
543 for l in 0..levels[f] {
544 g[off[f] + l] -= u[f][l] * uprec[f];
545 s[(off[f] + l) * dim + off[f] + l] += uprec[f];
546 }
547 }
548 let mut g_e: Vec<f64> = Vec::new();
550 d_e.clear();
551 b_e.clear();
552 if let Some(e) = elim {
553 for (l, rows) in members.iter().enumerate() {
554 let mut d = uprec[e];
555 let mut gl = -u[e][l] * uprec[e];
556 let mut b: Vec<(usize, f64)> = Vec::new();
557 for &i in rows {
558 d += wt[i];
559 gl += r[i];
560 let o = &obs[i];
561 for k in 0..K {
562 add_sparse(&mut b, k, wt[i] * o.x[k]);
563 }
564 for f in (0..F).filter(|&f| f != e) {
565 if let Some(lv) = o.lv[f] {
566 add_sparse(&mut b, off[f] + lv, wt[i]);
567 }
568 }
569 }
570 for &(a, va) in &b {
571 for &(c, vc) in &b {
572 s[a * dim + c] -= va * vc / d;
573 }
574 g[a] -= va * gl / d;
575 }
576 d_e.push(d);
577 g_e.push(gl);
578 b_e.push(b);
579 }
580 }
581 let mut held: Vec<usize> = (0..K).filter(|&k| spec[k].nonpos && beta[k] >= 0.0 && g[k] > 0.0).collect();
584 let step = loop {
585 let mut delta = vec![0.0; dim];
586 for &k in &held {
587 delta[k] = -beta[k];
588 }
589 let free: Vec<usize> = (0..dim).filter(|i| !held.contains(i)).collect();
590 let n = free.len();
591 let mut a = vec![0.0; n * n];
592 let mut rhs = vec![0.0; n];
593 for (fi, &i) in free.iter().enumerate() {
594 rhs[fi] = g[i] - held.iter().map(|&h| s[i * dim + h] * delta[h]).sum::<f64>();
595 for (fj, &j) in free.iter().enumerate() {
596 a[fi * n + fj] = s[i * dim + j];
597 }
598 }
599 if !cholesky(&mut a, n) {
600 return None;
601 }
602 chol_solve(&a, n, &mut rhs);
603 for (fi, &i) in free.iter().enumerate() {
604 delta[i] = rhs[fi];
605 }
606 let out: Vec<usize> = (0..K).filter(|&k| spec[k].nonpos && !held.contains(&k) && beta[k] + delta[k] > 0.0).collect();
607 if out.is_empty() {
608 break delta;
609 }
610 held.extend(out);
611 };
612 let delta_e: Vec<f64> = (0..d_e.len()).map(|l| (g_e[l] - b_e[l].iter().map(|&(a, v)| v * step[a]).sum::<f64>()) / d_e[l]).collect();
613 let largest = step.iter().chain(&delta_e).fold(0.0f64, |m, x| m.max(x.abs()));
614 let mut t = 1.0;
616 let mut accepted = false;
617 for _ in 0..=MAX_HALVINGS {
618 let nb: [f64; K] = std::array::from_fn(|k| {
619 let x = beta[k] + t * step[k];
620 if spec[k].nonpos { x.min(0.0) } else { x }
621 });
622 let mut nu = u.clone();
623 for f in 0..F {
624 for l in 0..levels[f] {
625 nu[f][l] += t * if Some(f) == elim { delta_e[l] } else { step[off[f] + l] };
626 }
627 }
628 let next = objective(&nb, &nu);
629 if next >= current - 1e-12 * current.abs().max(1.0) {
631 (beta, u, current) = (nb, nu, next);
632 accepted = true;
633 break;
634 }
635 t *= 0.5;
636 }
637 if !accepted || largest * t < TOL {
638 converged = accepted || largest < TOL;
639 break;
640 }
641 }
642 if !cholesky(&mut s, dim) {
644 return None;
645 }
646 let inv = chol_inverse(&s, dim);
647 let cov: [[f64; K]; K] = std::array::from_fn(|a| std::array::from_fn(|b| inv[a * dim + b]));
648 let var_u: [Vec<f64>; F] = std::array::from_fn(|f| {
649 if Some(f) == elim {
650 (0..levels[f])
651 .map(|l| {
652 let b = &b_e[l];
653 let q: f64 = b.iter().map(|&(a, va)| b.iter().map(|&(c, vc)| va * inv[a * dim + c] * vc).sum::<f64>()).sum();
654 1.0 / d_e[l] + q / (d_e[l] * d_e[l])
655 })
656 .collect()
657 } else {
658 (0..levels[f]).map(|l| inv[(off[f] + l) * dim + off[f] + l]).collect()
659 }
660 });
661 Some(Est { beta, u, sd: *sd, spread: *spread, cov, var_u, converged })
662}
663
664fn add_sparse(b: &mut Vec<(usize, f64)>, i: usize, v: f64) {
666 if v == 0.0 {
667 return;
668 }
669 match b.iter_mut().find(|(j, _)| *j == i) {
670 Some(e) => e.1 += v,
671 None => b.push((i, v)),
672 }
673}
674
675fn eb<const K: usize, const F: usize>(obs: &[Obs<K, F>], levels: [usize; F], spec: &[Coef; K], warm: Option<&Est<K, F>>) -> Option<Est<K, F>> {
679 let mut sd: [f64; K] = std::array::from_fn(|k| warm.filter(|_| spec[k].ard).map_or(spec[k].sd, |w| w.sd[k]));
680 let mut spread: [f64; F] = std::array::from_fn(|f| warm.map_or(LEVEL_START, |w| w.spread[f]));
681 let mut est = map(obs, levels, spec, &sd, &spread, warm)?;
682 for _ in 0..MAX_EB {
683 let nsd: [f64; K] = std::array::from_fn(|k| match spec[k].ard {
684 true => ((est.beta[k] - spec[k].mean).powi(2) + est.cov[k][k]).sqrt().clamp(SD_FLOOR, PRIOR_SD),
685 false => sd[k],
686 });
687 let nspread: [f64; F] = std::array::from_fn(|f| match levels[f] {
688 0 => spread[f],
689 n => (est.u[f].iter().zip(&est.var_u[f]).map(|(x, v)| x * x + v).sum::<f64>() / n as f64).sqrt().clamp(SD_FLOOR, PRIOR_SD),
690 });
691 let moved = sd.iter().zip(&nsd).chain(spread.iter().zip(&nspread)).any(|(a, b)| (a - b).abs() > EB_TOL * a.abs().max(SD_FLOOR));
692 (sd, spread) = (nsd, nspread);
693 est = map(obs, levels, spec, &sd, &spread, Some(&est))?;
694 if !moved {
695 break;
696 }
697 }
698 Some(est)
699}
700
701#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
704pub struct Rates {
705 pub free_stop: f64,
707 pub free_live: f64,
709 pub later: f64,
712 pub wait: f64,
714 pub ask_live: f64,
720}
721
722#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
724pub struct PoolFit {
725 #[serde(default, skip_serializing_if = "Vec::is_empty")]
728 pub yields: Vec<(f64, f64)>,
729 #[serde(default, skip_serializing_if = "is_default")]
731 pub per_lesson: f64,
732 #[serde(default, skip_serializing_if = "is_default")]
735 pub latency: (f64, f64, f64),
736 #[serde(default, skip_serializing_if = "is_default")]
737 pub latency_n: u32,
738 #[serde(default, skip_serializing_if = "is_default")]
741 pub eta: f64,
742 #[serde(default, skip_serializing_if = "is_default")]
743 pub slope: (f64, f64),
744 #[serde(default, skip_serializing_if = "is_default")]
745 pub not: (f64, f64),
746}
747
748#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
750pub struct Model {
751 pub beta: Vec<f64>,
754 pub chol: Vec<f64>,
755 pub sd: Vec<f64>,
757 pub spread: [f64; FAMILIES],
759 pub sessions: BTreeMap<String, (f64, f64)>,
761 pub lessons: BTreeMap<String, (f64, f64)>,
762 pub judges: BTreeMap<String, (f64, f64)>,
763 #[serde(default, skip_serializing_if = "String::is_empty")]
765 pub judge: String,
766 #[serde(default, skip_serializing_if = "std::ops::Not::not")]
767 pub converged: bool,
768}
769
770impl Model {
771 #[must_use]
775 pub fn current(&self) -> bool {
776 self.beta.len() == P && self.chol.len() == P * P && self.sd.len() == P
777 }
778}
779
780#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
783pub struct Fit {
784 #[serde(default, skip_serializing_if = "String::is_empty")]
787 pub key: String,
788 #[serde(default, skip_serializing_if = "Option::is_none")]
790 pub model: Option<Model>,
791 #[serde(default, skip_serializing_if = "Option::is_none")]
793 pub half_life: Option<u32>,
794 #[serde(default, skip_serializing_if = "is_default")]
795 pub rates: Rates,
796 #[serde(default, skip_serializing_if = "is_default")]
797 pub pool: PoolFit,
798 #[serde(default, skip_serializing_if = "Option::is_none")]
801 pub budget: Option<f64>,
802 #[serde(default, skip_serializing_if = "is_default")]
803 pub lambda: f64,
804 #[serde(default, skip_serializing_if = "std::ops::Not::not")]
807 pub sensitive: bool,
808 #[serde(default, skip_serializing_if = "is_default")]
810 pub labels: u32,
811 #[serde(default, skip_serializing_if = "is_default")]
812 pub dropped: u32,
813}
814
815fn is_default<T: Default + PartialEq>(v: &T) -> bool {
816 *v == T::default()
817}
818
819pub(super) fn decay(age: u32, h: Option<u32>) -> f64 {
821 h.map_or(1.0, |h| 0.5f64.powf(f64::from(age) / f64::from(h)))
822}
823
824struct Data {
827 obs: Vec<Obs<P, FAMILIES>>,
828 tick: Vec<u32>,
831 base_w: Vec<f64>,
832 unclear: Vec<bool>,
833 levels: [usize; FAMILIES],
834 names: [Vec<String>; FAMILIES],
835 lv_of: HashMap<usize, [Option<usize>; FAMILIES]>,
837 judged_p: HashMap<usize, f64>,
838 dropped: u32,
839}
840
841fn level_of(names: &mut [Vec<String>; FAMILIES], index: &mut [HashMap<String, usize>; FAMILIES], tagged: &(usize, String)) -> usize {
843 let (f, name) = tagged;
844 *index[*f].entry(name.clone()).or_insert_with(|| {
845 names[*f].push(name.clone());
846 names[*f].len() - 1
847 })
848}
849
850fn model_route(a: &Apply) -> Option<Route> {
853 match a.route {
854 Route::Ask | Route::Free if a.delivered != super::summary::Delivered::Dropped => Some(a.route),
855 Route::Placebo => Some(Route::Not),
856 _ => None,
857 }
858}
859
860fn applies_of<'a>(s: &'a Summary, c: &Logged) -> &'a [Apply] {
862 s.applies.get(c.first as usize..c.first as usize + usize::from(c.n)).unwrap_or(&[])
863}
864
865fn logged_had(s: &Summary, c: &Logged) -> Vec<String> {
868 applies_of(s, c).iter().filter(|a| a.nofree).filter_map(|a| s.lessons.get(a.lesson as usize).cloned()).collect()
869}
870
871fn check_view<'a>(s: &'a Summary, c: &Logged) -> (Vec<Unit<'a>>, Vec<Feat>, Vec<Route>) {
873 let applies = applies_of(s, c);
874 let units: Vec<Unit> =
875 applies.iter().map(|a| Unit { lesson: s.lessons.get(a.lesson as usize).map_or("", String::as_str), grade: a.grade, rank: a.rank as usize }).collect();
876 let feats = Feat::of(&units, c.stop, c.tokens);
877 let routes = applies.iter().map(|a| if a.route == Route::Placebo { Route::Not } else { a.route }).collect();
878 (units, feats, routes)
879}
880
881impl Data {
882 fn of(s: &Summary, zero: bool) -> Data {
885 let now = s.sessions.len() as u32;
886 let from = now.saturating_sub(FIT_SESSIONS as u32);
887 let mut names: [Vec<String>; FAMILIES] = Default::default();
888 let mut index: [HashMap<String, usize>; FAMILIES] = Default::default();
889 let judged_p = judged(s);
890 let mut d = Data {
891 obs: Vec::new(),
892 tick: Vec::new(),
893 base_w: Vec::new(),
894 unclear: Vec::new(),
895 levels: [0; FAMILIES],
896 names: Default::default(),
897 lv_of: HashMap::new(),
898 judged_p: HashMap::new(),
899 dropped: 0,
900 };
901 let names_of = |c: &Logged, a: &Apply| -> [Option<(usize, String)>; FAMILIES] {
902 [
903 Some((SESSION, s.sessions.get(c.tick as usize).map_or("", |x| x.id.as_str()).to_string())),
904 s.lessons.get(a.lesson as usize).map(|n| (LESSON, n.clone())),
905 a.judge.and_then(|j| s.judges.get(j as usize).cloned()).map(|n| (JUDGE, n)),
906 ]
907 };
908 for c in &s.checks {
909 let (_, feats, _) = check_view(s, c);
910 let (m_q, m_f) = (usize::from(c.m_q), usize::from(c.m_f));
911 for (j, f) in feats.iter().enumerate() {
912 let i = c.first as usize + j;
913 let a = &s.applies[i];
914 let Some(route) = model_route(a) else { continue };
915 let y = match a.label {
916 Label::Changed => 1.0,
917 Label::Ignored => 0.0,
918 Label::Unclear if zero => 0.0,
919 _ => continue,
920 };
921 if c.tick < from {
922 d.dropped += 1;
923 continue;
924 }
925 let m = match route {
926 Route::Ask => m_q,
927 Route::Free => m_f,
928 _ => 1,
929 };
930 let p = judged_p.get(&i).copied().unwrap_or(1.0);
931 let w = 1.0 / p.max(1.0 / MAX_IPW);
932 let lv = names_of(c, a).map(|n| n.map(|n| level_of(&mut names, &mut index, &n)));
933 d.obs.push(Obs { x: design(f, route, m), y, w, lv });
934 d.tick.push(c.tick);
935 d.base_w.push(w);
936 d.unclear.push(a.label == Label::Unclear);
937 }
938 }
939 for c in &s.checks {
941 for i in c.first as usize..(c.first as usize + c.n as usize).min(s.applies.len()) {
942 let n = names_of(c, &s.applies[i]);
943 d.lv_of.insert(i, std::array::from_fn(|f| n[f].as_ref().and_then(|(_, n)| index[f].get(n).copied())));
944 }
945 }
946 d.levels = std::array::from_fn(|f| names[f].len());
947 d.names = names;
948 d.judged_p = judged_p;
949 d
950 }
951
952 fn without_levels(mut self) -> Data {
954 for o in &mut self.obs {
955 o.lv = [None; FAMILIES];
956 }
957 for lv in self.lv_of.values_mut() {
958 *lv = [None; FAMILIES];
959 }
960 self.levels = [0; FAMILIES];
961 self.names = Default::default();
962 self
963 }
964
965 fn weighted(&self, now: u32, h: Option<u32>) -> Vec<Obs<P, FAMILIES>> {
967 self.obs.iter().zip(&self.tick).zip(&self.base_w).map(|((o, t), w)| Obs { w: w * decay(now.saturating_sub(*t), h), ..o.clone() }).collect()
968 }
969}
970
971pub(crate) fn judged(s: &Summary) -> HashMap<usize, f64> {
976 const K: usize = 6;
978 let mut last: HashMap<u32, usize> = HashMap::new();
979 for (i, c) in s.checks.iter().enumerate() {
980 last.insert(c.tick, i);
981 }
982 let mut rows: Vec<(usize, Obs<K, 0>)> = Vec::new();
983 for (ci, c) in s.checks.iter().enumerate() {
984 for i in c.first as usize..(c.first as usize + c.n as usize).min(s.applies.len()) {
985 let a = &s.applies[i];
986 if model_route(a).is_none() || a.label == Label::Open {
987 continue;
988 }
989 let y = f64::from(u8::from(matches!(a.label, Label::Changed | Label::Ignored)));
990 let x = [
991 1.0,
992 f64::from(u8::from(a.route == Route::Ask)),
993 f64::from(u8::from(a.route == Route::Free)),
994 f64::from(u8::from(c.stop)),
995 (1.0 + f64::from(a.pending)).ln(),
996 f64::from(u8::from(last.get(&c.tick) == Some(&ci))),
997 ];
998 rows.push((i, Obs { x, y, w: 1.0, lv: [] }));
999 }
1000 }
1001 if rows.is_empty() {
1002 return HashMap::new();
1003 }
1004 let spec = [Coef { mean: 0.0, sd: PRIOR_SD, ard: false, nonpos: false }; K];
1005 let obs: Vec<Obs<K, 0>> = rows.iter().map(|(_, o)| o.clone()).collect();
1006 let Some(est) = map(&obs, [], &spec, &[PRIOR_SD; K], &[], None) else { return HashMap::new() };
1007 rows.iter().map(|(i, o)| (*i, sigmoid(est.eta(o)))).collect()
1008}
1009
1010fn loglik(est: &Est<P, FAMILIES>, obs: &[Obs<P, FAMILIES>]) -> f64 {
1012 obs.iter()
1013 .map(|o| {
1014 let p = sigmoid(est.eta(o)).clamp(1e-12, 1.0 - 1e-12);
1016 o.w * (o.y * p.ln() + (1.0 - o.y) * (1.0 - p).ln())
1017 })
1018 .sum()
1019}
1020
1021fn half_life(d: &Data, base: &Est<P, FAMILIES>, spec: &[Coef; P]) -> Option<u32> {
1025 let mut ticks: Vec<u32> = d.tick.clone();
1026 ticks.sort_unstable();
1027 ticks.dedup();
1028 let eval: Vec<u32> = ticks.iter().skip(1).rev().take(PREQ_SESSIONS).rev().copied().collect();
1029 if eval.is_empty() {
1030 return None;
1031 }
1032 let mut best: (f64, Option<u32>) = (f64::NEG_INFINITY, None);
1033 for h in HALF_LIVES {
1034 let mut score = 0.0;
1035 let mut warm: Option<Est<P, FAMILIES>> = None;
1036 for &t in &eval {
1037 let train: Vec<Obs<P, FAMILIES>> = d
1038 .obs
1039 .iter()
1040 .zip(&d.tick)
1041 .zip(&d.base_w)
1042 .filter(|((_, tk), _)| **tk < t)
1043 .map(|((o, tk), w)| Obs { w: w * decay(t - tk, h), ..o.clone() })
1044 .collect();
1045 let test: Vec<Obs<P, FAMILIES>> = d.obs.iter().zip(&d.tick).filter(|(_, tk)| **tk == t).map(|(o, _)| o.clone()).collect();
1046 let Some(est) = map(&train, d.levels, spec, &base.sd, &base.spread, warm.as_ref()) else { continue };
1047 score += loglik(&est, &test);
1048 warm = Some(est);
1049 }
1050 if score > best.0 + 1e-9 * best.0.abs().max(1.0) || best.0 == f64::NEG_INFINITY {
1052 best = (score, h);
1053 }
1054 }
1055 best.1
1056}
1057
1058fn model_of(est: &Est<P, FAMILIES>, d: &Data, judge: String) -> Option<Model> {
1061 let mut chol: Vec<f64> = est.cov.iter().flatten().copied().collect();
1062 let finite = est.beta.iter().chain(&est.sd).chain(&est.spread).chain(est.u.iter().flatten()).chain(est.var_u.iter().flatten()).all(|x| x.is_finite());
1063 if !finite || !cholesky(&mut chol, P) {
1064 return None;
1065 }
1066 let levels = |f: usize| -> BTreeMap<String, (f64, f64)> {
1067 d.names[f].iter().enumerate().map(|(l, n)| (n.clone(), (est.u[f][l], est.var_u[f][l].max(0.0).sqrt()))).collect()
1068 };
1069 Some(Model {
1070 beta: est.beta.to_vec(),
1071 chol,
1072 sd: est.sd.to_vec(),
1073 spread: est.spread,
1074 sessions: levels(SESSION),
1075 lessons: levels(LESSON),
1076 judges: levels(JUDGE),
1077 judge,
1078 converged: est.converged,
1079 })
1080}
1081
1082#[must_use]
1085pub fn fit(s: &Summary, budget: Option<f64>) -> Fit {
1086 fit_with(s, budget, &coefs(), true)
1087}
1088
1089#[must_use]
1095pub fn fit_pruned(s: &Summary, budget: Option<f64>) -> Fit {
1096 fit_with(s, budget, &pruned(), false)
1097}
1098
1099fn fit_with(s: &Summary, budget: Option<f64>, spec: &[Coef; P], levels: bool) -> Fit {
1101 let data = |zero: bool| if levels { Data::of(s, zero) } else { Data::of(s, zero).without_levels() };
1102 let d = data(false);
1103 let now = s.sessions.len() as u32;
1104 let mut fit = Fit { key: s.key(), labels: d.obs.len() as u32, dropped: d.dropped, budget: budget.filter(|b| *b > 0.0), ..Fit::default() };
1105 let base = if d.obs.is_empty() { None } else { eb(&d.obs, d.levels, spec, None) };
1106 if let Some(base) = &base {
1107 fit.half_life = half_life(&d, base, spec);
1108 let est = eb(&d.weighted(now, fit.half_life), d.levels, spec, Some(base));
1109 let judge = s.checks.iter().rev().flat_map(|c| applies_of(s, c)).find_map(|a| a.judge);
1110 let judge = judge.filter(|_| levels).and_then(|j| s.judges.get(j as usize)).cloned().unwrap_or_default();
1111 if let Some(est) = est {
1112 let z = data(true);
1113 if z.unclear.iter().any(|u| *u)
1114 && let Some(alt) = eb(&z.weighted(now, fit.half_life), z.levels, spec, None)
1115 {
1116 fit.sensitive = [ASK, NOT].iter().any(|&k| est.beta[k].signum() != alt.beta[k].signum());
1117 }
1118 let eta = d.obs.iter().map(|o| est.eta(o)).sum::<f64>() / d.obs.len() as f64;
1119 fit.pool.eta = eta;
1120 fit.pool.slope = (est.beta[POOL], est.cov[POOL][POOL].max(0.0).sqrt());
1121 fit.pool.not = (est.beta[NOT], est.cov[NOT][NOT].max(0.0).sqrt());
1122 fit.model = model_of(&est, &d, judge).map(|m| if levels { m } else { Model { spread: [0.0; FAMILIES], ..m } });
1124 }
1125 }
1126 fit.rates = rates(s, fit.half_life);
1127 pool_state(s, &mut fit.pool);
1128 if let Some(b) = fit.budget {
1129 fit.lambda = price(s, &fit, b);
1130 }
1131 fit
1132}
1133
1134#[must_use]
1139pub fn prior() -> Fit {
1140 prior_with(&coefs(), true)
1141}
1142
1143#[must_use]
1147pub fn prior_pruned() -> Fit {
1148 prior_with(&pruned(), false)
1149}
1150
1151fn prior_with(spec: &[Coef; P], levels: bool) -> Fit {
1154 let s = Summary::default();
1155 let d = Data::of(&s, false);
1156 let sd: [f64; P] = std::array::from_fn(|k| spec[k].sd);
1157 let model = map(&d.obs, d.levels, spec, &sd, &[LEVEL_START; FAMILIES], None)
1158 .and_then(|est| model_of(&est, &d, String::new()))
1159 .map(|m| if levels { m } else { Model { spread: [0.0; FAMILIES], ..m } });
1160 Fit { model, rates: rates(&s, None), ..Fit::default() }
1161}
1162
1163const PROBIT: f64 = std::f64::consts::PI / 8.0;
1166
1167#[must_use]
1173pub fn predict(fit: &Fit, f: &Feat, route: Route, m: usize, session: &str, lesson: &str) -> Option<f64> {
1174 let model = fit.model.as_ref().filter(|m| m.current())?;
1175 let x = design(f, route, m);
1176 let lx = (0..P).map(|k| (k..P).map(|i| model.chol[i * P + k] * x[i]).sum::<f64>());
1178 let mut var: f64 = lx.map(|v| v * v).sum();
1179 let mut mean: f64 = model.beta.iter().zip(&x).map(|(b, x)| b * x).sum();
1180 for (map, name, spread) in [
1181 (&model.sessions, session, model.spread[SESSION]),
1182 (&model.lessons, lesson, model.spread[LESSON]),
1183 (&model.judges, model.judge.as_str(), model.spread[JUDGE]),
1184 ] {
1185 let (mu, sd) = map.get(name).copied().unwrap_or((0.0, spread));
1186 mean += mu;
1187 var += sd * sd;
1188 }
1189 Some(sigmoid(mean / (1.0 + PROBIT * var).sqrt()))
1190}
1191
1192fn rates(s: &Summary, h: Option<u32>) -> Rates {
1195 let mean = |(k, miss): (f64, f64)| (1.0 + k) / (2.0 + k + miss);
1196 let [free_stop, free_live, later, wait, ask] = s.tallies(h);
1197 let ask_live = if ask.0 + ask.1 > 0.0 { mean(ask) } else { 1.0 };
1199 Rates { free_stop: mean(free_stop), free_live: mean(free_live), later: mean(later), wait: mean(wait), ask_live }
1200}
1201
1202fn pool_state(s: &Summary, pool: &mut PoolFit) {
1205 let bins: Vec<(f64, f64)> = s.pool.bins.iter().take_while(|(n, _)| *n > 0).map(|&(n, k)| (1.0 + k as f64, 1.0 + n.saturating_sub(k) as f64)).collect();
1209 pool.yields = pava(&bins);
1210 pool.per_lesson = if s.pool.shown > 0 { s.pool.tokens as f64 / s.pool.shown as f64 } else { 0.0 };
1211 let obs = &s.pool.latency;
1212 pool.latency_n = obs.len() as u32;
1213 if obs.len() >= 2 {
1214 let n = obs.len() as f64;
1215 let mx = obs.iter().map(|o| f64::from(o.0)).sum::<f64>() / n;
1216 let my = obs.iter().map(|o| f64::from(o.1)).sum::<f64>() / n;
1217 let sxx: f64 = obs.iter().map(|o| (f64::from(o.0) - mx).powi(2)).sum();
1218 let sxy: f64 = obs.iter().map(|o| (f64::from(o.0) - mx) * (f64::from(o.1) - my)).sum();
1219 let b = if sxx > 0.0 { (sxy / sxx).max(0.0) } else { 0.0 };
1220 let a = my - b * mx;
1221 let mut res: Vec<f64> = obs.iter().map(|o| f64::from(o.1) - a - b * f64::from(o.0)).collect();
1222 res.sort_by(f64::total_cmp);
1223 let q = res[((1.0 - DELTA) * (res.len() - 1) as f64).ceil() as usize];
1224 pool.latency = (a, b, q.max(0.0));
1225 }
1226}
1227
1228fn pava(bins: &[(f64, f64)]) -> Vec<(f64, f64)> {
1231 let mut blocks: Vec<(f64, f64, usize)> = Vec::new();
1232 for &(a, b) in bins {
1233 blocks.push((a, b, 1));
1234 while blocks.len() >= 2 {
1235 let (x, y) = (blocks[blocks.len() - 2], blocks[blocks.len() - 1]);
1236 if x.0 / (x.0 + x.1) >= y.0 / (y.0 + y.1) {
1237 break;
1238 }
1239 blocks.pop();
1240 if let Some(last) = blocks.last_mut() {
1241 *last = (x.0 + y.0, x.1 + y.1, x.2 + y.2);
1242 }
1243 }
1244 }
1245 blocks.iter().flat_map(|&(a, b, n)| std::iter::repeat_n((a, b), n)).collect()
1246}
1247
1248#[derive(Debug, Clone)]
1251struct Theta {
1252 b: [f64; P],
1253 base: Vec<f64>,
1254}
1255
1256struct Levels<'a> {
1259 name: &'a str,
1260 session: u64,
1261 judge: u64,
1262 lessons: Vec<u64>,
1263}
1264
1265impl<'a> Levels<'a> {
1266 fn of(session: &'a str, judge: &str, units: &[Unit]) -> Levels<'a> {
1267 Levels {
1268 name: session,
1269 session: hash(&[b"session", session.as_bytes()]),
1270 judge: hash(&[b"judge", judge.as_bytes()]),
1271 lessons: units.iter().map(|u| hash(&[b"lesson", u.lesson.as_bytes()])).collect(),
1272 }
1273 }
1274}
1275
1276fn draw_beta(m: &Model, seed: u64) -> [f64; P] {
1283 let mut rng = Rng::new(seed);
1284 let mut z: [f64; P] = std::array::from_fn(|_| normal(&mut rng));
1285 let spec = coefs();
1286 let l = |i: usize, k: usize| m.chol.get(i * P + k).copied().unwrap_or(0.0);
1287 let mut x = [0.0; P];
1288 for i in 0..P {
1289 let given = m.beta.get(i).copied().unwrap_or(0.0) + (0..i).map(|k| l(i, k) * z[k]).sum::<f64>();
1290 x[i] = given + l(i, i) * z[i];
1291 if spec[i].nonpos {
1292 x[i] = below_zero(given, l(i, i), z[i]);
1293 if l(i, i) > 0.0 {
1295 z[i] = (x[i] - given) / l(i, i);
1296 }
1297 }
1298 }
1299 x
1300}
1301
1302fn below_zero(mean: f64, sd: f64, z: f64) -> f64 {
1306 if sd <= 0.0 || !sd.is_finite() {
1307 return mean.min(0.0);
1308 }
1309 let bound = -mean / sd;
1310 let ln_lower = ln_phi(z) + ln_phi(bound);
1312 let upper = phi(-z) + phi(z) * phi(-bound);
1313 (mean + sd * probit(ln_lower, upper)).min(0.0)
1314}
1315
1316fn phi(x: f64) -> f64 {
1318 0.5 * erfc(-x / std::f64::consts::SQRT_2)
1319}
1320
1321fn ln_phi(x: f64) -> f64 {
1323 if x < 0.0 { 0.5f64.ln() - 0.5 * x * x + erfcx(-x / std::f64::consts::SQRT_2).ln() } else { (-0.5 * erfc(x / std::f64::consts::SQRT_2)).ln_1p() }
1324}
1325
1326fn erfc(x: f64) -> f64 {
1328 let a = x.abs();
1329 let tail = (-a * a).exp() * erfcx(a);
1330 if x >= 0.0 { tail } else { 2.0 - tail }
1331}
1332
1333fn erfcx(x: f64) -> f64 {
1336 const C: [f64; 10] =
1337 [-1.265_512_23, 1.000_023_68, 0.374_091_96, 0.096_784_18, -0.186_288_06, 0.278_868_07, -1.135_203_98, 1.488_515_87, -0.822_152_23, 0.170_872_77];
1338 let t = 1.0 / (1.0 + 0.5 * x);
1339 t * C.iter().rev().fold(0.0, |acc, c| acc * t + c).exp()
1340}
1341
1342fn probit(ln_lower: f64, upper: f64) -> f64 {
1346 const A: [f64; 6] =
1347 [-3.969_683_028_665_376e1, 2.209_460_984_245_205e2, -2.759_285_104_469_687e2, 1.383_577_518_672_69e2, -3.066_479_806_614_716e1, 2.506_628_277_459_239];
1348 const B: [f64; 5] = [-5.447_609_879_822_406e1, 1.615_858_368_580_409e2, -1.556_989_798_598_866e2, 6.680_131_188_771_972e1, -1.328_068_155_288_572e1];
1349 const C: [f64; 6] =
1350 [-7.784_894_002_430_293e-3, -3.223_964_580_411_365e-1, -2.400_758_277_161_838, -2.549_732_539_343_734, 4.374_664_141_464_968, 2.938_163_982_698_783];
1351 const D: [f64; 4] = [7.784_695_709_041_462e-3, 3.224_671_290_700_398e-1, 2.445_134_137_142_996, 3.754_408_661_907_416];
1352 const LOW: f64 = 0.02425;
1353 let poly = |c: &[f64], x: f64| c.iter().fold(0.0, |acc, k| acc * x + k);
1354 let tail = |ln_p: f64| {
1356 let q = (-2.0 * ln_p).sqrt();
1357 poly(&C, q) / (poly(&D, q) * q + 1.0)
1358 };
1359 let lower = ln_lower.exp();
1360 if lower < LOW {
1361 tail(ln_lower)
1362 } else if upper < LOW {
1363 -tail(upper.ln())
1364 } else {
1365 let q = lower - 0.5;
1366 let r = q * q;
1367 poly(&A, r) * q / (poly(&B, r) * r + 1.0)
1368 }
1369}
1370
1371fn theta_with(m: &Model, b: [f64; P], feats: &[Feat], units: &[Unit], lv: &Levels, seed: u64) -> Theta {
1376 theta_by(m, b, feats, units, lv, |(mean, sd), h| mean + sd * normal(&mut Rng::new(mix(seed, h))))
1377}
1378
1379fn theta_mean(m: &Model, feats: &[Feat], units: &[Unit], lv: &Levels) -> Theta {
1385 let spec = coefs();
1386 let b = std::array::from_fn(|i| {
1387 let x = m.beta.get(i).copied().unwrap_or(0.0);
1388 if spec[i].nonpos { x.min(0.0) } else { x }
1389 });
1390 theta_by(m, b, feats, units, lv, |(mean, _), _| mean)
1391}
1392
1393fn theta_by(m: &Model, b: [f64; P], feats: &[Feat], units: &[Unit], lv: &Levels, level: impl Fn((f64, f64), u64) -> f64) -> Theta {
1396 let effect = |map: &BTreeMap<String, (f64, f64)>, name: &str, spread: f64, h: u64| level(map.get(name).copied().unwrap_or((0.0, spread)), h);
1397 let shared = effect(&m.sessions, lv.name, m.spread[SESSION], lv.session) + effect(&m.judges, &m.judge, m.spread[JUDGE], lv.judge);
1398 let base = feats
1399 .iter()
1400 .zip(units)
1401 .zip(&lv.lessons)
1402 .map(|((f, u), h)| dot(&b, &design(f, Route::Free, 1)) + shared + effect(&m.lessons, u.lesson, m.spread[LESSON], *h))
1403 .collect();
1404 Theta { b, base }
1405}
1406
1407fn theta(m: &Model, feats: &[Feat], units: &[Unit], lv: &Levels, seed: u64) -> Theta {
1409 theta_with(m, draw_beta(m, seed), feats, units, lv, seed)
1410}
1411
1412struct Valuer<'a> {
1414 t: &'a Theta,
1415 feats: &'a [Feat],
1416 rates: &'a Rates,
1417}
1418
1419impl Valuer<'_> {
1420 fn n(&self) -> usize {
1421 self.feats.len()
1422 }
1423
1424 fn stop(&self) -> bool {
1425 self.feats.first().is_some_and(|f| f.stop)
1426 }
1427
1428 fn p_not(&self, j: usize) -> f64 {
1430 let b = &self.t.b;
1431 sigmoid(self.t.base[j] + b[NOT] + b[NOT_G] * (f64::from(self.feats[j].grade) - 2.0))
1432 }
1433
1434 fn ask(&self, j: usize, m: usize) -> f64 {
1435 let (b, f) = (&self.t.b, &self.feats[j]);
1436 let eta = self.t.base[j]
1437 + b[ASK]
1438 + b[ASK_M] * (m.max(1) as f64).ln()
1439 + b[ASK_STOP] * f64::from(u8::from(f.stop))
1440 + b[ASK_G3] * f64::from(u8::from(f.grade == 3));
1441 let rate = match (self.stop(), f.three) {
1446 (true, _) => 1.0,
1447 (false, true) => self.rates.ask_live,
1448 (false, false) => self.rates.free_live,
1449 };
1450 rate * (sigmoid(eta) - self.p_not(j))
1451 }
1452
1453 fn free(&self, j: usize, m: usize) -> f64 {
1454 let rate = if self.stop() { self.rates.free_stop } else { self.rates.free_live };
1455 rate * (sigmoid(self.t.base[j] + self.t.b[FREE_M] * (m.max(1) as f64).ln()) - self.p_not(j))
1456 }
1457
1458 fn value(&self, routes: &[Route]) -> f64 {
1464 let m_q = routes.iter().filter(|r| **r == Route::Ask).count();
1465 let m_f = routes.iter().filter(|r| **r == Route::Free).count();
1466 routes
1467 .iter()
1468 .enumerate()
1469 .map(|(j, r)| match r {
1470 Route::Ask => self.ask(j, m_q),
1471 Route::Free => self.free(j, m_f),
1472 Route::Placebo | Route::Not => 0.0,
1473 })
1474 .sum()
1475 }
1476}
1477
1478fn forced(f: &Feat) -> bool {
1481 f.stop && f.grade == 3
1482}
1483
1484struct Table {
1487 n: usize,
1488 ask: Vec<f64>,
1492 free: Vec<f64>,
1493 forced: Vec<bool>,
1494 nofree: Vec<bool>,
1495}
1496
1497impl Table {
1498 fn of(v: &Valuer, nofree: &[bool]) -> Table {
1501 let n = v.n();
1502 let mut ask = vec![0.0; (n + 1) * n];
1503 let mut free = vec![0.0; (n + 1) * n];
1504 for m in 1..=n {
1505 for j in 0..n {
1506 ask[m * n + j] = v.ask(j, m);
1507 free[m * n + j] = v.free(j, m);
1508 }
1509 }
1510 Table { n, ask, free, forced: v.feats.iter().map(forced).collect(), nofree: (0..n).map(|j| nofree.get(j).copied().unwrap_or(false)).collect() }
1511 }
1512
1513 #[cfg(test)]
1516 fn value(&self, routes: &[Route]) -> f64 {
1517 let m_q = routes.iter().filter(|r| **r == Route::Ask).count();
1518 let m_f = routes.iter().filter(|r| **r == Route::Free).count();
1519 let n = self.n;
1520 routes
1521 .iter()
1522 .enumerate()
1523 .map(|(j, r)| match r {
1524 Route::Ask => self.ask[m_q * n + j],
1525 Route::Free => self.free[m_f * n + j],
1526 Route::Placebo | Route::Not => 0.0,
1527 })
1528 .sum()
1529 }
1530}
1531
1532#[derive(Debug, Clone, Copy, PartialEq)]
1534struct Total(f64);
1535
1536impl Eq for Total {}
1537
1538impl PartialOrd for Total {
1539 fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
1540 Some(self.cmp(other))
1541 }
1542}
1543
1544impl Ord for Total {
1545 fn cmp(&self, other: &Self) -> std::cmp::Ordering {
1546 self.0.total_cmp(&other.0)
1547 }
1548}
1549
1550fn top_sums(values: impl Iterator<Item = f64>, k: usize, out: &mut Vec<f64>, heap: &mut std::collections::BinaryHeap<std::cmp::Reverse<Total>>) {
1553 out.clear();
1554 heap.clear();
1555 out.push(if k == 0 { 0.0 } else { f64::NEG_INFINITY });
1556 let mut sum = 0.0;
1557 for x in values {
1558 if k > 0 {
1559 if heap.len() < k {
1560 heap.push(std::cmp::Reverse(Total(x)));
1561 sum += x;
1562 } else if let Some(min) = heap.peek().map(|r| r.0.0)
1563 && x > min
1564 {
1565 heap.pop();
1566 heap.push(std::cmp::Reverse(Total(x)));
1567 sum += x - min;
1568 }
1569 }
1570 out.push(if heap.len() < k { f64::NEG_INFINITY } else { sum });
1571 }
1572}
1573
1574fn solve(t: &Table, rail: usize, cap: usize) -> Vec<Route> {
1582 let n = t.n;
1583 let nf_forced = t.forced.iter().filter(|x| **x).count();
1584 let open: Vec<usize> = (0..n).filter(|&j| !t.forced[j]).collect();
1585 let held = open.iter().filter(|&&j| t.nofree[j]).count();
1588 let q_max = rail.max(nf_forced).min(n);
1589 let f_max = cap.min(open.len() - held);
1590 let ask = |m: usize, j: usize| if m == 0 { 0.0 } else { t.ask[m * n + j] };
1592 let free = |m: usize, j: usize| if m == 0 { 0.0 } else { t.free[m * n + j] };
1593 let mut order = open.clone();
1594 let mut keys = vec![0.0; n];
1595 let (mut pre, mut suf) = (Vec::with_capacity(n + 1), Vec::with_capacity(n + 1));
1596 let mut heap = std::collections::BinaryHeap::with_capacity(n + 1);
1597 let sort = |order: &mut Vec<usize>, keys: &mut Vec<f64>, mq: usize, mf: usize| {
1598 for &j in &open {
1599 keys[j] = if t.nofree[j] { f64::INFINITY } else { ask(mq, j) - free(mf, j) };
1600 }
1601 order.clear();
1602 order.extend_from_slice(&open);
1603 order.sort_by(|&a, &b| keys[b].total_cmp(&keys[a]).then(a.cmp(&b)));
1604 };
1605 let mut best: Option<(f64, usize, usize, usize)> = None;
1606 for mq in nf_forced..=q_max {
1607 let k = mq - nf_forced;
1608 if k > open.len() {
1609 break;
1610 }
1611 let fixed: f64 = (0..n).filter(|&j| t.forced[j]).map(|j| ask(mq, j)).sum();
1612 for mf in 0..=f_max.min(open.len() - k) {
1613 sort(&mut order, &mut keys, mq, mf);
1614 top_sums(order.iter().map(|&j| ask(mq, j)), k, &mut pre, &mut heap);
1615 top_sums(order.iter().rev().map(|&j| free(mf, j)), mf, &mut suf, &mut heap);
1616 for split in k.max(held)..=open.len() - mf {
1617 let total = fixed + pre[split] + suf[open.len() - split];
1618 if best.is_none_or(|b| total > b.0) {
1619 best = Some((total, mq, mf, split));
1620 }
1621 }
1622 }
1623 }
1624 let mut routes = vec![Route::Not; n];
1625 for j in (0..n).filter(|&j| t.forced[j]) {
1626 routes[j] = Route::Ask;
1627 }
1628 let Some((_, mq, mf, split)) = best else { return routes };
1629 sort(&mut order, &mut keys, mq, mf);
1630 let pick = |within: &[usize], gain: &dyn Fn(usize) -> f64, k: usize| -> Vec<usize> {
1631 let mut c = within.to_vec();
1632 c.sort_by(|&a, &b| gain(b).total_cmp(&gain(a)).then(a.cmp(&b)));
1633 c.truncate(k);
1634 c
1635 };
1636 for j in pick(&order[..split], &|j| ask(mq, j), mq - nf_forced) {
1637 routes[j] = Route::Ask;
1638 }
1639 for j in pick(&order[split..], &|j| free(mf, j), mf) {
1640 routes[j] = Route::Free;
1641 }
1642 routes
1643}
1644
1645#[must_use]
1648pub fn fixed_routes(units: &[Unit], had: &[String]) -> Vec<Route> {
1649 let ls: Vec<Challenge> = units.iter().map(|u| Challenge { id: u.lesson.to_string(), ..Default::default() }).collect();
1650 let applied: Vec<(&Challenge, u8)> = ls.iter().zip(units).map(|(c, u)| (c, u.grade)).collect();
1651 let (ask, free) = conscience::graded_without(&applied, None, had);
1652 let holds = |v: &[&Challenge], c: &Challenge| v.iter().any(|x| std::ptr::eq(*x, c));
1653 ls.iter()
1654 .map(|c| match (holds(&ask, c), holds(&free, c)) {
1655 (true, _) => Route::Ask,
1656 (_, true) => Route::Free,
1657 _ => Route::Not,
1658 })
1659 .collect()
1660}
1661
1662fn nofree_of(units: &[Unit], had: &[String]) -> Vec<bool> {
1664 units.iter().map(|u| had.iter().any(|h| h == u.lesson)).collect()
1665}
1666
1667#[derive(Debug, Clone, Copy)]
1669pub struct Check<'a> {
1670 pub session: &'a str,
1671 pub stop: bool,
1672 pub tokens: u32,
1674}
1675
1676pub const FIXED: &str = "fixed";
1680pub const SIZE: &str = "size";
1682pub const GATE: &str = "gate";
1685pub const EVIDENCE: &str = "evidence";
1689pub const EXPERIMENT: &str = "experiment";
1692pub const EXPERIMENT_CONTROL: &str = "experiment_control";
1695
1696#[derive(Debug, Clone, PartialEq)]
1707pub struct Decision {
1708 pub routes: Vec<Route>,
1709 pub p: Vec<f64>,
1710 pub pf: Vec<f64>,
1711 pub source: &'static str,
1712 pub gate: f64,
1713 pub risk: f64,
1714 pub xk: Option<u32>,
1715 pub xs: Option<u32>,
1716}
1717
1718fn lower_quantile(v: &mut [f64]) -> f64 {
1720 if v.is_empty() {
1721 return 0.0;
1722 }
1723 v.sort_by(f64::total_cmp);
1724 v[((DELTA * v.len() as f64).ceil() as usize).clamp(1, v.len()) - 1]
1725}
1726
1727#[cfg(test)]
1730#[must_use]
1731pub(crate) fn remaining(gate: &Gate) -> f64 {
1732 gate.allowance + lower_quantile(&mut gate.gain.clone())
1733}
1734
1735fn allowance_at(gate: &Gate, alpha: f64) -> f64 {
1738 if alpha > 0.0 { gate.allowance * (alpha / ALPHA) } else { 0.0 }
1741}
1742
1743#[must_use]
1746pub fn remaining_at(gate: &Gate, alpha: f64) -> f64 {
1747 (allowance_at(gate, alpha) + lower_quantile(&mut gate.gain.clone())).max(0.0)
1748}
1749
1750fn gate_seed(fit: &Fit, k: usize) -> u64 {
1753 mix(hash(&[b"gate", fit.key.as_bytes()]), k as u64)
1754}
1755
1756struct Draws<'a> {
1759 thetas: Vec<Theta>,
1760 feats: Vec<Feat>,
1761 rates: &'a Rates,
1762 base: Vec<f64>,
1763}
1764
1765impl<'a> Draws<'a> {
1766 fn of(model: &Model, fit: &'a Fit, units: &[Unit], check: &Check, lv: &Levels, base: &[Route]) -> Draws<'a> {
1767 let feats = Feat::of(units, check.stop, check.tokens);
1768 let thetas: Vec<Theta> = (0..DRAWS).map(|k| theta(model, &feats, units, lv, gate_seed(fit, k))).collect();
1769 let mut d = Draws { thetas, feats, rates: &fit.rates, base: Vec::new() };
1770 d.base = (0..DRAWS).map(|k| d.value(k, base)).collect();
1771 d
1772 }
1773
1774 fn value(&self, k: usize, routes: &[Route]) -> f64 {
1776 Valuer { t: &self.thetas[k], feats: &self.feats, rates: self.rates }.value(routes)
1777 }
1778
1779 fn gains(&self, routes: &[Route]) -> Vec<f64> {
1781 self.base.iter().enumerate().map(|(k, b)| self.value(k, routes) - b).collect()
1782 }
1783}
1784
1785fn evident(model: &Model, d: &Draws, units: &[Unit], lv: &Levels, base: &[Route], nofree: &[bool]) -> Option<Vec<Route>> {
1789 let t = theta_mean(model, &d.feats, units, lv);
1790 let best = solve(&Table::of(&Valuer { t: &t, feats: &d.feats, rates: d.rates }, nofree), crate::attention::RAIL, conscience::CAP);
1791 (best != base && lower_quantile(&mut d.gains(&best)) > 0.0).then_some(best)
1792}
1793
1794pub(crate) fn tested(units: &[Unit], base: &[Route], had: &[String]) -> Vec<usize> {
1799 let mut order: Vec<usize> =
1800 (0..units.len()).filter(|&j| units[j].grade != 3 && matches!(base[j], Route::Ask | Route::Free) && !had.iter().any(|h| h == units[j].lesson)).collect();
1801 order.sort_by_key(|&j| std::cmp::Reverse(units[j].grade));
1802 order
1803}
1804
1805fn test_routes(base: &[Route], order: &[usize], s: usize) -> Vec<Route> {
1808 let mut routes = base.to_vec();
1809 for (i, &j) in order.iter().enumerate() {
1810 routes[j] = if i < s { Route::Ask } else { Route::Free };
1811 }
1812 routes
1813}
1814
1815fn draw_count(seed: u64, k: usize) -> usize {
1819 (hash(&[b"question-count", &seed.to_le_bytes()]) % (k as u64 + 1)) as usize
1820}
1821
1822#[must_use]
1838pub fn choose(units: &[Unit], check: &Check, fit: &Fit, gate: &Gate, alpha: f64, seed: u64, had: &[String]) -> Decision {
1839 let base = fixed_routes(units, had);
1840 let fixed = |source: &'static str| Decision {
1841 p: share(&[&base], Route::Ask),
1842 pf: share(&[&base], Route::Free),
1843 routes: base.clone(),
1844 source,
1845 gate: 0.0,
1846 risk: 0.0,
1847 xk: None,
1848 xs: None,
1849 };
1850 let Some(model) = fit.model.as_ref().filter(|m| m.current()) else { return fixed(FIXED) };
1851 if gate.allowance <= 0.0 || gate.gain.len() != DRAWS || units.is_empty() {
1852 return fixed(FIXED);
1853 }
1854 if units.len() > SOLVE_MAX {
1855 return fixed(SIZE);
1856 }
1857 let lv = Levels::of(check.session, &model.judge, units);
1858 let d = Draws::of(model, fit, units, check, &lv, &base);
1859 let allowance = allowance_at(gate, alpha);
1860 let held = lower_quantile(&mut gate.gain.clone());
1861 let remaining = remaining_at(gate, alpha);
1862 let admit = |routes: &[Route]| -> Option<f64> {
1866 if routes == base {
1867 return Some(0.0);
1868 }
1869 let mut after: Vec<f64> = d.gains(routes).iter().zip(&gate.gain).map(|(x, g)| g + x).collect();
1870 let q = lower_quantile(&mut after);
1871 (q >= -allowance).then_some((held - q).max(0.0))
1872 };
1873 if let Some(best) = evident(model, &d, units, &lv, &base, &nofree_of(units, had)) {
1874 return match admit(&best) {
1875 Some(risk) => Decision {
1876 p: share(&[&best], Route::Ask),
1877 pf: share(&[&best], Route::Free),
1878 routes: best,
1879 source: EVIDENCE,
1880 gate: remaining,
1881 risk,
1882 xk: None,
1883 xs: None,
1884 },
1885 None => Decision { gate: remaining, ..fixed(GATE) },
1886 };
1887 }
1888 let order = tested(units, &base, had);
1889 if check.stop || alpha <= 0.0 || order.is_empty() {
1890 return Decision { gate: remaining, ..fixed(EVIDENCE) };
1891 }
1892 let k = order.len();
1895 let arms: Vec<(Vec<Route>, Option<f64>)> = (0..=k)
1896 .map(|s| {
1897 let routes = test_routes(&base, &order, s);
1898 let risk = admit(&routes);
1899 (if risk.is_some() { routes } else { base.clone() }, risk)
1900 })
1901 .collect();
1902 let played: Vec<&Vec<Route>> = arms.iter().map(|(r, _)| r).collect();
1903 let s = draw_count(seed, k);
1904 let (routes, risk) = &arms[s];
1905 let source = match risk {
1906 None => GATE,
1907 Some(_) if *routes == base => EXPERIMENT_CONTROL,
1908 Some(_) => EXPERIMENT,
1909 };
1910 let count = |x: usize| Some(u32::try_from(x).unwrap_or(u32::MAX));
1911 Decision {
1912 p: share(&played, Route::Ask),
1913 pf: share(&played, Route::Free),
1914 routes: routes.clone(),
1915 source,
1916 gate: remaining,
1917 risk: risk.unwrap_or(0.0),
1918 xk: count(k),
1919 xs: count(s),
1920 }
1921}
1922
1923#[must_use]
1932pub fn decide_on_evidence(units: &[Unit], check: &Check, fit: &Fit, had: &[String]) -> Vec<Route> {
1933 let base = fixed_routes(units, had);
1934 let Some(model) = fit.model.as_ref().filter(|m| m.current()) else { return base };
1935 if units.is_empty() || units.len() > SOLVE_MAX {
1936 return base;
1937 }
1938 let lv = Levels::of(check.session, &model.judge, units);
1939 let d = Draws::of(model, fit, units, check, &lv, &base);
1940 evident(model, &d, units, &lv, &base, &nofree_of(units, had)).unwrap_or(base)
1941}
1942
1943fn share(draws: &[&Vec<Route>], route: Route) -> Vec<f64> {
1946 let n = draws.first().map_or(0, |d| d.len());
1947 (0..n).map(|j| crate::outcome::round_p(draws.iter().filter(|d| d[j] == route).count() as f64 / draws.len() as f64)).collect()
1948}
1949
1950pub fn ledger_add_without(fit: &Fit, gate: &mut Gate, check: &Check, units: &[Unit], routes: &[Route], had: &[String]) {
1956 let Some(model) = fit.model.as_ref().filter(|m| m.current()) else { return };
1957 let betas: Vec<[f64; P]> = (0..DRAWS).map(|k| draw_beta(model, gate_seed(fit, k))).collect();
1958 add_with(fit, &betas, gate, check, units, routes, had);
1959}
1960
1961fn add_with(fit: &Fit, betas: &[[f64; P]], gate: &mut Gate, check: &Check, units: &[Unit], routes: &[Route], had: &[String]) {
1964 let Some(model) = fit.model.as_ref().filter(|m| m.current()) else { return };
1965 if units.is_empty() || units.len() != routes.len() {
1966 return;
1967 }
1968 if gate.gain.len() != DRAWS {
1969 gate.gain = vec![0.0; DRAWS];
1970 }
1971 let base = fixed_routes(units, had);
1972 let routes: Vec<Route> = routes.iter().map(|r| if *r == Route::Placebo { Route::Not } else { *r }).collect();
1974 let feats = Feat::of(units, check.stop, check.tokens);
1975 let lv = Levels::of(check.session, &model.judge, units);
1976 let mut mean = 0.0;
1977 for (k, (g, b)) in gate.gain.iter_mut().zip(betas).enumerate() {
1978 let t = theta_with(model, *b, &feats, units, &lv, gate_seed(fit, k));
1979 let v = Valuer { t: &t, feats: &feats, rates: &fit.rates };
1980 let fixed = v.value(&base);
1981 *g += v.value(&routes) - fixed;
1982 mean += fixed / DRAWS as f64;
1983 }
1984 gate.allowance += ALPHA * mean.max(0.0);
1985 gate.checks += 1;
1986 if routes != base {
1987 gate.deviated += 1;
1988 }
1989}
1990
1991#[must_use]
1994pub fn ledger(s: &Summary) -> Gate {
1995 let mut gate = Gate::default();
1996 let Some(model) = s.fit.model.as_ref().filter(|m| m.current()) else { return gate };
1997 let betas: Vec<[f64; P]> = (0..DRAWS).map(|k| draw_beta(model, gate_seed(&s.fit, k))).collect();
1998 for c in s.checks.iter().filter(|c| c.learned) {
1999 let (units, _, routes) = check_view(s, c);
2000 let session = s.sessions.get(c.tick as usize).map_or("", |x| x.id.as_str());
2001 add_with(&s.fit, &betas, &mut gate, &Check { session, stop: c.stop, tokens: c.tokens }, &units, &routes, &logged_had(s, c));
2002 }
2003 gate
2004}
2005
2006#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
2008pub struct Monitor {
2009 #[serde(default, skip_serializing_if = "std::ops::Not::not")]
2011 pub reverted: bool,
2012 #[serde(default, skip_serializing_if = "is_default")]
2014 pub sessions: u32,
2015 #[serde(default, skip_serializing_if = "is_default")]
2016 pub units: u32,
2017 #[serde(default, skip_serializing_if = "is_default")]
2021 pub mean: f64,
2022 #[serde(default, skip_serializing_if = "Option::is_none")]
2023 pub ucb: Option<f64>,
2024 #[serde(default, skip_serializing_if = "is_default")]
2028 pub labelled: u32,
2029 #[serde(default, skip_serializing_if = "is_default")]
2030 pub ess: f64,
2031 #[serde(default, skip_serializing_if = "is_default")]
2032 pub max_weight: f64,
2033 #[serde(default, skip_serializing_if = "is_default")]
2034 pub truncated: u32,
2035}
2036
2037#[must_use]
2044pub fn monitor(s: &Summary) -> Monitor {
2045 let spec = coefs();
2046 let d = Data::of(s, false);
2047 let Some(model) = &s.fit.model else { return Monitor::default() };
2048 let sd: [f64; P] = std::array::from_fn(|k| model.sd.get(k).copied().unwrap_or(spec[k].sd));
2049 let spread = model.spread;
2050 let folds: Vec<Option<Est<P, FAMILIES>>> = (0..2u32)
2051 .map(|half| {
2052 let train: Vec<Obs<P, FAMILIES>> =
2053 d.weighted(s.sessions.len() as u32, s.fit.half_life).into_iter().zip(&d.tick).filter(|(_, t)| *t % 2 != half).map(|(o, _)| o).collect();
2054 map(&train, d.levels, &spec, &sd, &spread, None)
2055 })
2056 .collect();
2057 struct Term {
2059 tick: u32,
2060 diff: [f64; ARMS],
2061 p_hat: [f64; ARMS],
2062 taken: usize,
2063 y: Option<f64>,
2064 w: f64,
2065 vb: f64,
2066 }
2067 let mut terms: Vec<Term> = Vec::new();
2068 for c in s.checks.iter().filter(|c| c.learned) {
2069 let Some(est) = &folds[(c.tick % 2) as usize] else { continue };
2070 let (units, feats, routes) = check_view(s, c);
2071 let base = fixed_routes(&units, &logged_had(s, c));
2073 let (m_q, m_f) = (usize::from(c.m_q).max(1), usize::from(c.m_f).max(1));
2074 for (j, f) in feats.iter().enumerate() {
2075 let i = c.first as usize + j;
2076 let a = &s.applies[i];
2077 let lv = d.lv_of.get(&i).copied().unwrap_or([None; FAMILIES]);
2078 let eta = |route: Route, m: usize| -> f64 {
2080 dot(&est.beta, &design(f, route, m))
2081 + lv[1].map_or(0.0, |l| est.u[LESSON].get(l).copied().unwrap_or(0.0))
2082 + lv[2].map_or(0.0, |l| est.u[JUDGE].get(l).copied().unwrap_or(0.0))
2083 };
2084 let p_hat = [sigmoid(eta(Route::Ask, m_q)), sigmoid(eta(Route::Free, m_f)), sigmoid(eta(Route::Not, 1))];
2085 let mu = [a.p, a.pf, (1.0 - a.p - a.pf).max(0.0)];
2086 let (taken, fixed) = (routes[j].arm(), base[j].arm());
2087 let diff: [f64; ARMS] = std::array::from_fn(|x| mu[x] - f64::from(u8::from(x == fixed)));
2088 let y = match a.label {
2089 Label::Changed => Some(1.0),
2090 Label::Ignored => Some(0.0),
2091 _ => None,
2092 };
2093 let psi = if a.route == Route::Placebo { a.psi.max(PSI_FLOOR) } else { 1.0 };
2096 let judged_p = d.judged_p.get(&i).copied().unwrap_or(1.0).max(1.0 / MAX_IPW);
2097 let w = if y.is_some() { 1.0 / (mu[taken].max(1.0 / DRAWS as f64) * judged_p * psi) } else { 0.0 };
2098 let vb = p_hat[fixed] - p_hat[2];
2099 terms.push(Term { tick: c.tick, diff, p_hat, taken, y, w, vb });
2100 }
2101 }
2102 let labelled: Vec<f64> = terms.iter().filter(|t| t.y.is_some()).map(|t| t.w).collect();
2103 if terms.is_empty() || labelled.is_empty() {
2104 return Monitor::default();
2105 }
2106 let cap = (labelled.len() as f64).sqrt() * labelled.iter().sum::<f64>() / labelled.len() as f64;
2107 let mut norm = [0.0f64; ARMS];
2108 for t in &terms {
2109 norm[t.taken] += t.w.min(cap);
2110 }
2111 let count = terms.len() as f64;
2112 let norm = norm.map(|x| if x > 0.0 { x / count } else { 1.0 });
2113 let mut per: BTreeMap<u32, (f64, usize)> = BTreeMap::new();
2114 for t in &terms {
2115 let w = t.w.min(cap) / norm[t.taken];
2116 let dr: f64 = (0..ARMS).map(|a| t.diff[a] * (t.p_hat[a] + if a == t.taken { t.y.map_or(0.0, |y| w * (y - t.p_hat[a])) } else { 0.0 })).sum();
2117 let e = per.entry(t.tick).or_insert((0.0, 0));
2118 e.0 += dr + ALPHA * t.vb;
2119 e.1 += 1;
2120 }
2121 let z: Vec<f64> = per.values().map(|(x, n)| x / *n as f64).collect();
2122 let mean = z.iter().sum::<f64>() / z.len() as f64;
2123 let ucb = Some(upper_bound(&z)).filter(|u| u.is_finite());
2124 let used: Vec<f64> = terms.iter().filter(|t| t.y.is_some()).map(|t| t.w.min(cap) / norm[t.taken]).collect();
2126 let (sum, sq) = used.iter().fold((0.0, 0.0), |(s, q), w| (s + w, q + w * w));
2127 let ess = if sq > 0.0 { sum * sum / sq } else { 0.0 };
2128 let max_weight = if sum > 0.0 { used.iter().copied().fold(0.0, f64::max) * used.len() as f64 / sum } else { 0.0 };
2129 Monitor {
2130 reverted: ucb.is_some_and(|u| u < 0.0),
2131 sessions: z.len() as u32,
2132 units: terms.len() as u32,
2133 mean,
2134 ucb,
2135 labelled: used.len() as u32,
2136 ess,
2137 max_weight,
2138 truncated: labelled.iter().filter(|w| **w > cap).count() as u32,
2139 }
2140}
2141
2142fn upper_bound(x: &[f64]) -> f64 {
2151 let log = (1.0 / DELTA).ln();
2152 let shape = 2.0 * log + (1.0 + 2.0 * log).ln();
2153 let (mut sum, mut sum_sq) = (0.0, 0.0);
2154 let mut best = f64::INFINITY;
2155 for (i, &xi) in x.iter().enumerate() {
2156 let t = (i + 1) as f64;
2157 sum += xi;
2158 sum_sq += xi * xi;
2159 if i + 1 < CS_MIN.max(2) {
2160 continue;
2161 }
2162 let mean = sum / t;
2163 let var = ((sum_sq - t * mean * mean) / (t - 1.0)).max(f64::MIN_POSITIVE);
2164 let (v, rho) = (t * var, var * CS_SESSIONS / shape);
2165 best = best.min(mean + ((v + rho) * ((1.0 + v / rho).ln() + 2.0 * log)).sqrt() / t);
2166 }
2167 best
2168}
2169
2170#[must_use]
2172pub fn session_seed(session: &str) -> u64 {
2173 hash(&[b"pool", session.as_bytes()])
2174}
2175
2176#[must_use]
2178pub fn pool_grid() -> Vec<usize> {
2179 POOL_STEPS.iter().map(|&i| (crate::whole::TOKENS as f64 * 2f64.powf(f64::from(i) / 2.0)).round() as usize).collect()
2180}
2181
2182fn pool_value(pool: &PoolFit, tokens: usize, seed: u64) -> f64 {
2186 let mut rng = Rng::new(seed);
2187 let shown = tokens as f64 / pool.per_lesson.max(1.0);
2188 let mut applies = 0.0;
2189 for (b, &(a, bb)) in pool.yields.iter().enumerate() {
2190 let (lo, hi) = ((1u64 << b) - 1, (1u64 << (b + 1)) - 1);
2191 let inside = (shown.min(hi as f64) - lo as f64).max(0.0);
2192 applies += inside * if inside > 0.0 { beta_draw(a, bb, &mut rng) } else { 0.0 };
2193 }
2194 let slope = (pool.slope.0 + pool.slope.1 * normal(&mut rng)).min(0.0);
2195 let not = pool.not.0 + pool.not.1 * normal(&mut rng);
2196 let eta = pool.eta + slope * (tokens as f64 / crate::whole::TOKENS as f64).ln();
2197 applies * (sigmoid(eta) - sigmoid(eta + not)).max(0.0)
2198}
2199
2200fn pool_feasible(pool: &PoolFit, tokens: usize, limit_ms: Option<u64>, lambda: f64) -> bool {
2204 let baseline = crate::whole::TOKENS;
2205 let bins_to = |t: usize| -> usize {
2206 let shown = (t as f64 / pool.per_lesson.max(1.0)).ceil() as u64;
2207 (u64::BITS - shown.leading_zeros()) as usize
2208 };
2209 let known = pool.per_lesson > 0.0 && pool.yields.len() >= bins_to(tokens.max(baseline));
2210 if tokens > baseline {
2211 let (a, b, q) = pool.latency;
2212 let fits = limit_ms.is_some_and(|l| a + b * tokens as f64 + q <= l as f64);
2213 return known && fits && pool.latency_n as usize >= LATENCY_MIN;
2214 }
2215 tokens == baseline || known || lambda > 0.0
2216}
2217
2218struct PoolDraws {
2221 grid: Vec<usize>,
2222 shared: Vec<Vec<f64>>,
2223}
2224
2225impl PoolDraws {
2226 fn of(pool: &PoolFit) -> PoolDraws {
2227 let grid = pool_grid();
2228 let shared = grid.iter().map(|&t| (0..DRAWS).map(|k| pool_value(pool, t, mix(hash(&[b"pool-gate"]), k as u64))).collect()).collect();
2229 PoolDraws { grid, shared }
2230 }
2231
2232 fn choose(&self, play: &[f64], feasible: &[bool], lambda: f64) -> usize {
2236 let baseline = crate::whole::TOKENS;
2237 let Some(b) = self.grid.iter().position(|&t| t == baseline) else { return baseline };
2238 let admits = |i: usize| -> bool {
2239 let cost = lambda * (self.grid[i] as f64 - baseline as f64);
2240 let mut gain: Vec<f64> = self.shared[i].iter().zip(&self.shared[b]).map(|(x, y)| x - y - cost).collect();
2241 i == b || lower_quantile(&mut gain) >= 0.0
2242 };
2243 (0..self.grid.len())
2244 .filter(|&i| feasible[i] && admits(i))
2245 .map(|i| (play[i] - lambda * self.grid[i] as f64, self.grid[i]))
2246 .fold(None, |best: Option<(f64, usize)>, x| if best.is_none_or(|b| x.0 > b.0) { Some(x) } else { best })
2247 .map_or(baseline, |(_, t)| t)
2248 }
2249
2250 fn play(&self, pool: &PoolFit, seed: u64) -> Vec<f64> {
2252 self.grid.iter().map(|&t| pool_value(pool, t, seed)).collect()
2253 }
2254}
2255
2256#[must_use]
2262pub fn pool_tokens(fit: &Fit, active: bool, limit_ms: Option<u64>, seed: u64) -> usize {
2263 let baseline = crate::whole::TOKENS;
2264 if !active || fit.pool.yields.is_empty() && fit.lambda == 0.0 {
2265 return baseline;
2266 }
2267 let draws = PoolDraws::of(&fit.pool);
2268 let feasible: Vec<bool> = draws.grid.iter().map(|&t| pool_feasible(&fit.pool, t, limit_ms, fit.lambda)).collect();
2269 draws.choose(&draws.play(&fit.pool, seed), &feasible, fit.lambda)
2270}
2271
2272#[must_use]
2281pub fn price(s: &Summary, fit: &Fit, budget: f64) -> f64 {
2282 let sessions: Vec<(&str, u32)> = s.sessions.iter().rev().take(FIT_SESSIONS).filter(|x| x.checks > 0).map(|x| (x.id.as_str(), x.checks)).collect();
2283 if sessions.is_empty() || budget <= 0.0 {
2284 return 0.0;
2285 }
2286 let draws = PoolDraws::of(&fit.pool);
2287 let plays: Vec<Vec<f64>> = sessions.iter().map(|(id, _)| draws.play(&fit.pool, session_seed(id))).collect();
2288 let spend = |lambda: f64| -> f64 {
2289 let feasible: Vec<bool> = draws.grid.iter().map(|&t| pool_feasible(&fit.pool, t, Some(u64::MAX), lambda)).collect();
2291 sessions.iter().zip(&plays).map(|((_, n), play)| f64::from(*n) * draws.choose(play, &feasible, lambda) as f64).sum::<f64>() / sessions.len() as f64
2292 };
2293 if spend(0.0) <= budget {
2294 return 0.0;
2295 }
2296 let mut hi = 1.0 / crate::whole::TOKENS as f64;
2298 for _ in 0..PRICE_STEPS {
2299 if spend(hi) <= budget {
2300 break;
2301 }
2302 hi *= 2.0;
2303 }
2304 let mut lo = 0.0;
2305 for _ in 0..PRICE_STEPS {
2306 let mid = 0.5 * (lo + hi);
2307 if spend(mid) <= budget { hi = mid } else { lo = mid }
2308 }
2309 hi
2310}
2311
2312#[cfg(test)]
2313mod tests {
2314 use super::*;
2315 use crate::outcome::{self, Outcome};
2316 use crate::select::summary::{self, Summary};
2317
2318 fn unit(lesson: &str, grade: u8, rank: usize) -> Unit<'_> {
2319 Unit { lesson, grade, rank }
2320 }
2321
2322 fn tables(n: usize) -> Vec<Vec<u8>> {
2324 let mut out = vec![Vec::new()];
2325 for len in 1..=n {
2326 out.extend((0..3usize.pow(len as u32)).map(|mut x| {
2327 (0..len)
2328 .map(|_| {
2329 let g = (x % 3) as u8 + 1;
2330 x /= 3;
2331 g
2332 })
2333 .collect()
2334 }));
2335 }
2336 out
2337 }
2338
2339 #[derive(Clone, Copy)]
2342 struct Truth {
2343 b: [f64; P],
2344 lift: f64,
2345 }
2346
2347 impl Truth {
2348 fn p(&self, f: &Feat, route: Route, m: usize) -> f64 {
2349 sigmoid(dot(&self.b, &design(f, route, m)) + if f.stop && f.three { self.lift } else { 0.0 })
2350 }
2351 }
2352
2353 struct Sim {
2357 rng: Rng,
2358 lines: Vec<Outcome>,
2359 ms: i64,
2360 }
2361
2362 type Router<'a> = dyn FnMut(&[u8], bool, &mut Rng) -> (Vec<Route>, Vec<f64>, Vec<f64>, &'static str) + 'a;
2363
2364 impl Sim {
2365 fn new(seed: u64) -> Sim {
2366 Sim { rng: Rng::new(seed), lines: Vec::new(), ms: 0 }
2367 }
2368
2369 fn session(&mut self, session: &str, checks: usize, lessons: usize, truth: &Truth, route: &mut Router) {
2370 for c in 0..checks {
2371 let n = 1 + (self.rng.next_u64() % 4) as usize;
2372 let stop = c % 4 == 3;
2373 let grades: Vec<u8> = (0..n).map(|_| 1 + (self.rng.next_u64() % 3) as u8).collect();
2374 let ids: Vec<String> = (0..n).map(|_| format!("L{}", self.rng.next_u64() % lessons as u64)).collect();
2375 self.check(session, stop, &grades, &ids, truth, route);
2376 }
2377 self.ms += 1;
2378 let asked: Vec<String> =
2379 self.lines.iter().filter(|o| o.session == session && o.kind == outcome::DELIVERY).flat_map(|o| o.lessons.clone()).collect();
2380 if let Some(l) = asked.first() {
2381 self.lines.push(Outcome {
2382 kind: outcome::NOT_CORRECTED.into(),
2383 session: session.into(),
2384 lessons: vec![l.clone()],
2385 ms: self.ms,
2386 reason: outcome::UNJUDGED_NONE.into(),
2387 ..Default::default()
2388 });
2389 }
2390 }
2391
2392 fn check(&mut self, session: &str, stop: bool, grades: &[u8], ids: &[String], truth: &Truth, route: &mut Router) {
2393 self.ms += 1;
2394 let check = crate::select::check_id(session, if stop { "stop" } else { "live" }, self.ms, 1);
2395 let (routes, p, pf, policy) = route(grades, stop, &mut self.rng);
2396 let pick = |r: Route| ids.iter().zip(&routes).filter(|(_, x)| **x == r).map(|(l, _)| l.clone()).collect::<Vec<_>>();
2397 self.lines.push(Outcome {
2398 kind: outcome::GRADED.into(),
2399 session: session.into(),
2400 lessons: ids.to_vec(),
2401 ms: self.ms,
2402 channel: if stop { "stop" } else { "live" }.into(),
2403 grades: grades.to_vec(),
2404 ranks: (0..ids.len() as u16).collect(),
2405 applied: ids.len() as u16,
2406 chosen: pick(Route::Ask),
2407 free: pick(Route::Free),
2408 placebo: pick(Route::Placebo),
2409 psi: 1.0,
2410 p,
2411 pf,
2412 check: check.clone(),
2413 policy: policy.into(),
2414 ..Default::default()
2415 });
2416 let three = grades.contains(&3);
2417 let m = |r: Route| routes.iter().filter(|x| **x == r).count();
2418 for (j, (l, r)) in ids.iter().zip(&routes).enumerate() {
2419 if *r == Route::Not {
2420 continue;
2421 }
2422 self.ms += 1;
2423 let delivered = match r {
2424 Route::Ask => outcome::ROUTE_ASKED,
2425 Route::Free => outcome::ROUTE_FREE,
2426 _ => outcome::ROUTE_PLACEBO,
2427 };
2428 self.lines.push(Outcome {
2429 kind: outcome::DELIVERY.into(),
2430 session: session.into(),
2431 lessons: vec![l.clone()],
2432 ms: self.ms,
2433 check: check.clone(),
2434 route: delivered.into(),
2435 ..Default::default()
2436 });
2437 let f = Feat { grade: grades[j], rank: j, stop, three, pool: 0.0 };
2438 let as_route = if *r == Route::Placebo { Route::Not } else { *r };
2439 let y = self.rng.unit() < truth.p(&f, as_route, m(*r).max(1));
2440 let channel = match r {
2441 Route::Ask => {
2442 if stop {
2443 "stop"
2444 } else {
2445 "live"
2446 }
2447 }
2448 Route::Free => "queued",
2449 _ => outcome::PLACEBO,
2450 };
2451 self.lines.push(Outcome {
2452 kind: if y { outcome::CHANGED } else { outcome::IGNORED }.into(),
2453 session: session.into(),
2454 lessons: vec![l.clone()],
2455 ms: self.ms,
2456 channel: channel.into(),
2457 check: check.clone(),
2458 judge: 2,
2459 ..Default::default()
2460 });
2461 }
2462 }
2463
2464 fn text(&self) -> String {
2465 self.lines.iter().map(|o| format!("{}\n", serde_json::to_string(o).unwrap())).collect()
2466 }
2467
2468 fn summary(&self) -> Summary {
2469 summary::fold(&self.text(), Summary::default())
2470 }
2471 }
2472
2473 fn uniform(grades: &[u8], stop: bool, rng: &mut Rng) -> (Vec<Route>, Vec<f64>, Vec<f64>, &'static str) {
2475 let routes: Vec<Route> = grades
2476 .iter()
2477 .map(|g| if stop && *g == 3 { Route::Ask } else { [Route::Ask, Route::Free, Route::Placebo][(rng.next_u64() % 3) as usize] })
2478 .collect();
2479 let p = grades.iter().map(|g| if stop && *g == 3 { 1.0 } else { 0.3333 }).collect();
2480 let pf = grades.iter().map(|g| if stop && *g == 3 { 0.0 } else { 0.3333 }).collect();
2481 (routes, p, pf, "fixed")
2482 }
2483
2484 fn fixed_rule(grades: &[u8], _: bool, _: &mut Rng) -> (Vec<Route>, Vec<f64>, Vec<f64>, &'static str) {
2486 let ids: Vec<String> = (0..grades.len()).map(|i| format!("x{i}")).collect();
2487 let units: Vec<Unit> = ids.iter().zip(grades).enumerate().map(|(r, (l, g))| unit(l, *g, r)).collect();
2488 let routes = fixed_routes(&units, &[]);
2489 let p = routes.iter().map(|r| f64::from(u8::from(*r == Route::Ask))).collect();
2490 let pf = routes.iter().map(|r| f64::from(u8::from(*r == Route::Free))).collect();
2491 (routes, p, pf, "fixed")
2492 }
2493
2494 fn truth() -> Truth {
2495 let mut b = [0.0; P];
2496 b[B0] = -1.0;
2497 b[G3] = 1.2;
2498 b[ASK] = 0.6;
2499 b[NOT] = -0.8;
2500 Truth { b, lift: 0.0 }
2501 }
2502
2503 fn planted(sessions: usize, checks: usize, seed: u64) -> Summary {
2504 let mut sim = Sim::new(seed);
2505 for i in 0..sessions {
2506 sim.session(&format!("s{i}"), checks, 30, &truth(), &mut uniform);
2507 }
2508 sim.summary()
2509 }
2510
2511 #[test]
2512 fn the_fit_recovers_a_planted_world() {
2513 let s = planted(24, 40, 1);
2514 let fit = fit(&s, None);
2515 let m = fit.model.as_ref().expect("a model");
2516 assert!(m.converged, "{m:?}");
2517 assert!(fit.labels > 600, "{}", fit.labels);
2518 let cov = |k: usize| (0..=k).map(|j| m.chol[k * P + j].powi(2)).sum::<f64>().sqrt();
2519 for k in [B0, G3, ASK, NOT] {
2520 let (got, want, sd) = (m.beta[k], truth().b[k], cov(k));
2521 assert!((got - want).abs() < 3.0 * sd + 0.15, "{}: {got:.3} vs {want} (sd {sd:.3})", NAMES[k]);
2522 }
2523 assert!(m.beta[ASK] > 0.0 && m.beta[NOT] < 0.0, "{:?}", m.beta);
2524 assert!(m.sd[RANK] < ARD_START, "{}", m.sd[RANK]);
2526 assert!(m.lessons.len() <= 30 && m.sessions.len() == 24);
2527 assert!(m.beta.iter().chain(&m.chol).all(|x| x.is_finite()));
2528 }
2529
2530 #[test]
2534 fn a_dated_prior_keeps_its_sd_where_the_labels_say_nothing() {
2535 let mut sim = Sim::new(12);
2536 for i in 0..16 {
2537 sim.session(&format!("s{i}"), 16, 12, &truth(), &mut fixed_rule);
2538 }
2539 let f = fit(&sim.summary(), None);
2540 let m = f.model.as_ref().expect("a model");
2541 for k in [ASK_M, FREE_M] {
2542 assert_eq!(m.sd[k], DILUTION_PRIOR.1, "{}: {:?}", NAMES[k], m.sd);
2543 }
2544 assert_eq!((m.sd[ASK], m.sd[NOT]), (ASK_PRIOR.1, NOT_PRIOR.1));
2545 assert!([RANK, ASK_STOP, ASK_G3, NOT_G, POOL].iter().any(|&k| m.sd[k] != ARD_START), "EM ran on the rest: {:?}", m.sd);
2546 }
2547
2548 #[test]
2552 fn the_pruned_prior_is_the_pruned_models_before_any_label() {
2553 let (full, pruned) = (prior(), prior_pruned());
2554 let (fm, pm) = (full.model.as_ref().expect("a prior"), pruned.model.as_ref().expect("a pruned prior"));
2555 let spec = coefs();
2556 for (k, name) in NAMES.iter().enumerate() {
2557 let want = if PRUNED_KEEPS.contains(&k) { (spec[k].mean, spec[k].sd) } else { (0.0, SD_FLOOR) };
2558 assert!((pm.beta[k] - want.0).abs() < 1e-9 && pm.sd[k] == want.1, "{name}: {} sd {}", pm.beta[k], pm.sd[k]);
2559 }
2560 assert_eq!((pm.spread, fm.spread), ([0.0; FAMILIES], [LEVEL_START; FAMILIES]));
2561 let f = Feat { grade: 2, rank: 3, stop: true, three: false, pool: 0.5 };
2562 let at = |fit: &Fit| predict(fit, &f, Route::Not, 1, "unseen", "L99").expect("a prediction");
2564 assert_ne!(at(&pruned), at(&full), "the full prior's rank, pool and level variance are not the pruned model's");
2565 assert_eq!(pruned.rates, full.rates);
2566 }
2567
2568 #[test]
2573 fn the_pruned_model_keeps_grade_route_and_dilution_alone() {
2574 let mut t = truth();
2575 t.b[RANK] = -0.6;
2576 t.b[STOP] = 0.8;
2577 let mut sim = Sim::new(4);
2578 for i in 0..20 {
2579 sim.session(&format!("s{i}"), 30, 20, &t, &mut uniform);
2580 }
2581 let s = sim.summary();
2582 let (full, pruned) = (fit(&s, None), fit_pruned(&s, None));
2583 let (fm, pm) = (full.model.as_ref().expect("a full model"), pruned.model.as_ref().expect("a pruned model"));
2584 assert!(fm.beta[RANK] < -0.2 && fm.beta[STOP] > 0.2 && !fm.sessions.is_empty(), "the scenario: {:?}", fm.beta);
2585 assert!(pm.current() && pm.converged && pruned.labels == full.labels && pruned.key == full.key);
2586 for (k, name) in NAMES.iter().enumerate() {
2587 if PRUNED_KEEPS.contains(&k) {
2588 assert!(pm.sd[k] > SD_FLOOR, "{name} is kept");
2589 } else {
2590 assert_eq!(pm.sd[k], SD_FLOOR, "{name}");
2591 assert!(pm.beta[k].abs() < 3.0 * SD_FLOOR, "{name}: {}", pm.beta[k]);
2592 }
2593 }
2594 assert!(pm.beta[ASK] > 0.0 && pm.beta[NOT] < 0.0, "{:?}", pm.beta);
2595 assert!(pm.sessions.is_empty() && pm.lessons.is_empty() && pm.judges.is_empty() && pm.judge.is_empty());
2596 assert_eq!(pm.spread, [0.0; FAMILIES]);
2597 assert!(HALF_LIVES.contains(&pruned.half_life));
2598 let f = Feat { grade: 2, rank: 3, stop: true, three: false, pool: 0.0 };
2599 let at = |session: &str, lesson: &str| predict(&pruned, &f, Route::Free, 2, session, lesson);
2600 assert_eq!(at("s0", "L1"), at("unseen", "L99"), "no level effect");
2601 assert_ne!(predict(&full, &f, Route::Free, 2, "s0", "L1"), predict(&full, &f, Route::Free, 2, "unseen", "L99"), "the full model's levels");
2602 }
2603
2604 #[test]
2608 fn eliminating_a_family_equals_solving_it_densely() {
2609 let mut rng = Rng::new(9);
2610 let rows: Vec<([f64; 3], usize, f64)> = (0..400)
2611 .map(|i| {
2612 let x = [1.0, rng.unit() * 2.0 - 1.0, f64::from(u8::from(rng.unit() < 0.4))];
2613 let l = i % 4;
2614 let eta = -0.3 + 0.8 * x[1] - 0.5 * x[2] + [0.6, -0.4, 0.1, 0.0][l];
2615 (x, l, f64::from(u8::from(rng.unit() < sigmoid(eta))))
2616 })
2617 .collect();
2618 let spread = 0.7;
2619 let spec3 = [Coef { mean: 0.0, sd: 2.5, ard: false, nonpos: false }; 3];
2620 let obs: Vec<Obs<3, 1>> = rows.iter().map(|(x, l, y)| Obs { x: *x, y: *y, w: 1.0 + f64::from(*l as u8) * 0.1, lv: [Some(*l)] }).collect();
2621 let a = map(&obs, [4], &spec3, &[2.5; 3], &[spread], None).unwrap();
2622 let mut spec7 = [Coef { mean: 0.0, sd: spread, ard: false, nonpos: false }; 7];
2623 spec7[..3].copy_from_slice(&spec3);
2624 let sd7: [f64; 7] = std::array::from_fn(|k| spec7[k].sd);
2625 let dense: Vec<Obs<7, 0>> = rows
2626 .iter()
2627 .map(|(x, l, y)| {
2628 let mut z = [0.0; 7];
2629 z[..3].copy_from_slice(x);
2630 z[3 + l] = 1.0;
2631 Obs { x: z, y: *y, w: 1.0 + f64::from(*l as u8) * 0.1, lv: [] }
2632 })
2633 .collect();
2634 let b = map(&dense, [], &spec7, &sd7, &[], None).unwrap();
2635 assert!(a.converged && b.converged);
2636 for k in 0..3 {
2637 assert!((a.beta[k] - b.beta[k]).abs() < 1e-6, "{k}: {} {}", a.beta[k], b.beta[k]);
2638 for j in 0..3 {
2639 assert!((a.cov[k][j] - b.cov[k][j]).abs() < 1e-6, "cov {k} {j}");
2640 }
2641 }
2642 for l in 0..4 {
2643 assert!((a.u[0][l] - b.beta[3 + l]).abs() < 1e-6, "level {l}");
2644 assert!((a.var_u[0][l] - b.cov[3 + l][3 + l]).abs() < 1e-6, "var {l}: {} {}", a.var_u[0][l], b.cov[3 + l][3 + l]);
2645 }
2646 }
2647
2648 #[test]
2649 fn a_bounded_slope_stops_at_zero_and_a_negative_one_is_learned() {
2650 let mut rng = Rng::new(3);
2651 let obs: Vec<Obs<3, 0>> = (0..600)
2652 .map(|_| {
2653 let x = [1.0, rng.unit() * 3.0, rng.unit() * 3.0];
2654 let eta = -0.5 + 0.9 * x[1] - 0.7 * x[2];
2655 Obs { x, y: f64::from(u8::from(rng.unit() < sigmoid(eta))), w: 1.0, lv: [] }
2656 })
2657 .collect();
2658 let bounded = Coef { mean: 0.0, sd: 2.5, ard: false, nonpos: true };
2659 let spec = [Coef { nonpos: false, ..bounded }, bounded, bounded];
2660 let e = map(&obs, [], &spec, &[2.5; 3], &[], None).unwrap();
2661 assert_eq!(e.beta[1], 0.0, "a positive effect under a ≤ 0 bound sits at the bound");
2662 assert!(e.beta[2] < -0.3, "{:?}", e.beta);
2663 }
2664
2665 fn brute(v: &Valuer, nofree: &[bool], rail: usize, cap: usize) -> f64 {
2667 let n = v.n();
2668 let mut best = f64::NEG_INFINITY;
2669 for code in 0..3usize.pow(n as u32) {
2670 let mut x = code;
2671 let routes: Vec<Route> = (0..n)
2672 .map(|_| {
2673 let r = [Route::Ask, Route::Free, Route::Not][x % 3];
2674 x /= 3;
2675 r
2676 })
2677 .collect();
2678 let (m_q, m_f) = (routes.iter().filter(|r| **r == Route::Ask).count(), routes.iter().filter(|r| **r == Route::Free).count());
2679 let forced_ok = v.feats.iter().zip(&routes).all(|(f, r)| !forced(f) || *r == Route::Ask);
2680 let free_ok = nofree.iter().zip(&routes).all(|(x, r)| !x || *r != Route::Free);
2681 let f = v.feats.iter().filter(|f| forced(f)).count();
2682 if forced_ok && free_ok && m_q <= rail.max(f) && m_f <= cap {
2683 best = best.max(v.value(&routes));
2684 }
2685 }
2686 best
2687 }
2688
2689 #[test]
2690 fn the_solver_equals_brute_force() {
2691 let mut rng = Rng::new(11);
2692 for case in 0..400 {
2693 let n = case % 8;
2694 let stop = case % 3 == 0;
2695 let feats: Vec<Feat> = (0..n).map(|r| Feat { grade: 1 + (rng.next_u64() % 3) as u8, rank: r, stop, three: false, pool: 0.0 }).collect();
2696 let mut b = [0.0; P];
2697 for x in &mut b {
2698 *x = 2.0 * normal(&mut rng);
2699 }
2700 b[ASK_M] = -rng.unit() * 2.0;
2701 b[FREE_M] = -rng.unit() * 2.0;
2702 let t = Theta { b, base: (0..n).map(|_| normal(&mut rng)).collect() };
2703 let rates = Rates { free_stop: rng.unit(), free_live: rng.unit(), later: rng.unit(), wait: 0.0, ask_live: 1.0 };
2704 let v = Valuer { t: &t, feats: &feats, rates: &rates };
2705 let nofree: Vec<bool> = (0..n).map(|_| case % 2 == 1 && rng.unit() < 0.4).collect();
2707 let table = Table::of(&v, &nofree);
2708 for (rail, cap) in [(20, 10), (2, 1), (1, 3), (0, 0)] {
2709 let got = solve(&table, rail, cap);
2710 let want = brute(&v, &nofree, rail, cap);
2711 assert!((v.value(&got) - want).abs() < 1e-12, "case {case} rail {rail} cap {cap}: {} vs {want} {got:?}", v.value(&got));
2712 assert!((table.value(&got) - v.value(&got)).abs() < 1e-12, "the table values as the model does");
2713 assert!(feats.iter().zip(&got).all(|(f, r)| !forced(f) || *r == Route::Ask));
2714 assert!(nofree.iter().zip(&got).all(|(x, r)| !x || *r != Route::Free), "case {case}: {nofree:?} {got:?}");
2715 }
2716 }
2717 }
2718
2719 #[test]
2720 fn the_value_of_a_route_is_the_models_lift() {
2721 let f = Feat { grade: 3, rank: 2, stop: false, three: true, pool: 0.3 };
2722 let mut b = [0.0; P];
2723 for (i, x) in b.iter_mut().enumerate() {
2724 *x = 0.1 * i as f64 - 0.6;
2725 }
2726 let base = dot(&b, &design(&f, Route::Free, 1)) + 0.25;
2727 let t = Theta { b, base: vec![base] };
2728 let rates = Rates { free_stop: 0.5, free_live: 0.8, later: 0.4, wait: 0.0, ask_live: 0.6 };
2729 let v = Valuer { t: &t, feats: std::slice::from_ref(&f), rates: &rates };
2730 let p = |r: Route, m: usize| sigmoid(dot(&b, &design(&f, r, m)) + 0.25);
2731 assert!((v.ask(0, 3) - 0.6 * (p(Route::Ask, 3) - p(Route::Not, 1))).abs() < 1e-12);
2733 assert!((v.free(0, 2) - 0.8 * (p(Route::Free, 2) - p(Route::Not, 1))).abs() < 1e-12);
2734 for r in [Route::Not, Route::Placebo] {
2737 assert_eq!(v.value(&[r]), 0.0, "during the turn a lesson left for later counts nothing here: {r:?}");
2738 }
2739 let stop = Feat { stop: true, ..f };
2740 let t = Theta { b, base: vec![dot(&b, &design(&stop, Route::Free, 1)) + 0.25] };
2741 let v = Valuer { t: &t, feats: std::slice::from_ref(&stop), rates: &rates };
2742 let p = |r: Route, m: usize| sigmoid(dot(&b, &design(&stop, r, m)) + 0.25);
2743 assert!((v.ask(0, 3) - (p(Route::Ask, 3) - p(Route::Not, 1))).abs() < 1e-12, "a stop's question blocks: it arrives");
2744 assert!((v.free(0, 2) - 0.5 * (p(Route::Free, 2) - p(Route::Not, 1))).abs() < 1e-12);
2745 for r in [Route::Not, Route::Placebo] {
2746 assert_eq!(v.value(&[r]), 0.0, "at the stop a lesson not delivered counts nothing: {r:?}");
2747 }
2748 }
2749
2750 #[test]
2751 fn a_live_question_without_a_3_arrives_as_a_hand_over_does() {
2752 let mut b = [0.0; P];
2753 for (i, x) in b.iter_mut().enumerate() {
2754 *x = 0.1 * i as f64 - 0.6;
2755 }
2756 let t = Theta { b, base: vec![0.25] };
2757 let rates = Rates { free_stop: 0.5, free_live: 0.3, later: 0.4, wait: 0.0, ask_live: 0.9 };
2758 let ask = |three: bool| {
2759 let f = Feat { grade: 2, rank: 0, stop: false, three, pool: 0.0 };
2760 Valuer { t: &t, feats: std::slice::from_ref(&f), rates: &rates }.ask(0, 1)
2761 };
2762 let lift = sigmoid(0.25 + b[ASK]) - sigmoid(0.25 + b[NOT]);
2764 assert!(lift.abs() > 0.01);
2765 assert!((ask(true) - 0.9 * lift).abs() < 1e-12);
2767 assert!((ask(false) - 0.3 * lift).abs() < 1e-12);
2770 }
2771
2772 fn open_gate() -> Gate {
2774 Gate { gain: vec![0.0; DRAWS], allowance: 5.0, ..Gate::default() }
2775 }
2776
2777 #[test]
2778 fn cold_start_is_the_fixed_rule() {
2779 let fit_model = fit(&planted(6, 20, 2), None);
2780 assert!(fit_model.model.is_some());
2781 let wrong_len = Gate { gain: vec![0.0; 3], allowance: 1.0, ..Gate::default() };
2782 let ids: Vec<String> = (0..14).map(|i| format!("L{i}")).collect();
2783 let mut grades = tables(5);
2784 grades.extend([vec![3; 12], [vec![1; 11], vec![2; 3]].concat()]);
2785 for g in &grades {
2786 let units: Vec<Unit> = ids.iter().zip(g).enumerate().map(|(r, (l, g))| unit(l, *g, r)).collect();
2787 let want = fixed_routes(&units, &[]);
2788 for (fit, gate) in [(&Fit::default(), &open_gate()), (&fit_model, &Gate::default()), (&fit_model, &wrong_len)] {
2789 for stop in [true, false] {
2790 let d = choose(&units, &Check { session: "s", stop, tokens: 0 }, fit, gate, ALPHA, 7, &[]);
2791 assert_eq!(d.routes, want, "{g:?}");
2792 assert_eq!(d.p, want.iter().map(|r| f64::from(u8::from(*r == Route::Ask))).collect::<Vec<_>>());
2793 assert_eq!(d.pf, want.iter().map(|r| f64::from(u8::from(*r == Route::Free))).collect::<Vec<_>>());
2794 assert_eq!((d.source, d.gate), ("fixed", 0.0));
2795 }
2796 }
2797 }
2798 }
2799
2800 #[test]
2801 fn n_zero_and_n_one() {
2802 let empty = fit(&Summary::default(), None);
2803 assert!(empty.model.is_none() && empty.labels == 0);
2804 let mut sim = Sim::new(5);
2805 let mut one = |_: &[u8], _: bool, _: &mut Rng| (vec![Route::Ask], vec![1.0], vec![0.0], "fixed");
2806 sim.check("s", false, &[2], &["A".to_string()], &truth(), &mut one);
2807 let f = fit(&sim.summary(), None);
2808 let m = f.model.as_ref().expect("one label is a model");
2809 assert_eq!(f.labels, 1);
2810 assert!(m.beta.iter().chain(&m.chol).all(|x| x.is_finite()), "{m:?}");
2811 let gate = open_gate();
2812 assert!(choose(&[], &Check { session: "s", stop: true, tokens: 0 }, &f, &gate, ALPHA, 1, &[]).routes.is_empty());
2813 for g in 1..=3 {
2814 let d = choose(&[unit("A", g, 0)], &Check { session: "s", stop: true, tokens: 0 }, &f, &gate, ALPHA, 1, &[]);
2815 assert_eq!(d.routes.len(), 1);
2816 assert!((d.p[0] + d.pf[0]) <= 1.0 + 1e-9);
2817 }
2818 }
2819
2820 #[test]
2823 fn decisions_are_deterministic_and_only_the_test_varies() {
2824 let f = fit(&planted(8, 20, 3), None);
2825 let gate = open_gate();
2826 let units = [unit("L1", 2, 0), unit("L2", 1, 1), unit("L3", 3, 2)];
2827 let check = Check { session: "new", stop: false, tokens: 0 };
2828 let a = choose(&units, &check, &f, &gate, ALPHA, 42, &[]);
2829 assert_eq!(a, choose(&units, &check, &f, &gate, ALPHA, 42, &[]));
2830 let plays: std::collections::HashSet<Vec<Route>> = (0..60).map(|s| choose(&units, &check, &f, &gate, ALPHA, s, &[]).routes).collect();
2831 assert!(plays.len() > 1, "the test's count varies with the check: {plays:?}");
2832 let stop = Check { stop: true, ..check };
2833 let at_stop: std::collections::HashSet<Vec<Route>> = (0..60).map(|s| choose(&units, &stop, &f, &gate, ALPHA, s, &[]).routes).collect();
2834 assert_eq!(at_stop.len(), 1, "no Thompson play: {at_stop:?}");
2835 }
2836
2837 #[test]
2838 fn propensities_match_the_frequency_of_the_play() {
2839 let f = fit(&planted(6, 15, 4), None);
2840 let gate = open_gate();
2841 let units = [unit("L1", 2, 0), unit("L2", 1, 1), unit("L3", 1, 2)];
2842 let check = Check { session: "new", stop: false, tokens: 0 };
2843 let n = 300;
2844 let (mut asked, mut free) = ([0.0; 3], [0.0; 3]);
2845 let (mut p, mut pf) = ([0.0; 3], [0.0; 3]);
2846 for seed in 0..n {
2847 let d = choose(&units, &check, &f, &gate, ALPHA, seed, &[]);
2848 for j in 0..3 {
2849 asked[j] += f64::from(u8::from(d.routes[j] == Route::Ask)) / n as f64;
2850 free[j] += f64::from(u8::from(d.routes[j] == Route::Free)) / n as f64;
2851 p[j] += d.p[j] / n as f64;
2852 pf[j] += d.pf[j] / n as f64;
2853 }
2854 }
2855 for j in 0..3 {
2856 for (freq, prop) in [(asked[j], p[j]), (free[j], pf[j])] {
2857 let se = (prop * (1.0 - prop) / n as f64).sqrt();
2858 assert!((freq - prop).abs() <= 4.0 * se + 0.01, "unit {j}: played {freq:.3}, propensity {prop:.3}");
2859 }
2860 }
2861 assert!(asked.iter().chain(&free).any(|x| *x > 0.05 && *x < 0.95), "the play varies: {asked:?} {free:?}");
2862 }
2863
2864 fn free_beats_ask() -> Fit {
2866 let mut sim = Sim::new(8);
2867 let mut t = truth();
2868 t.b[ASK] = -3.0;
2869 for i in 0..10 {
2870 sim.session(&format!("s{i}"), 12, 10, &t, &mut uniform);
2871 }
2872 fit(&sim.summary(), None)
2873 }
2874
2875 #[test]
2878 fn a_stop_grade_three_is_never_routed_away_from_asking() {
2879 let f = free_beats_ask();
2880 let gate = open_gate();
2881 let units = [unit("L1", 3, 0), unit("L2", 1, 1)];
2882 let live = choose(&units, &Check { session: "n", stop: false, tokens: 0 }, &f, &gate, ALPHA, 1, &[]);
2883 assert_eq!(live.routes[0], Route::Free, "during the turn the model hands it over: {live:?}");
2884 for seed in 0..20 {
2885 let d = choose(&units, &Check { session: "n", stop: true, tokens: 0 }, &f, &gate, ALPHA, seed, &[]);
2886 assert_eq!((d.routes[0], d.p[0], d.pf[0]), (Route::Ask, 1.0, 0.0), "{d:?}");
2887 }
2888 }
2889
2890 #[test]
2894 fn a_lesson_the_session_had_never_goes_free() {
2895 let f = free_beats_ask();
2896 let gate = open_gate();
2897 let units = [unit("L1", 2, 0), unit("L2", 1, 1), unit("L3", 1, 2)];
2898 let had = ["L2".to_string()];
2899 let check = Check { session: "n", stop: false, tokens: 0 };
2900 assert_eq!(fixed_routes(&units, &had), [Route::Free, Route::Not, Route::Free]);
2901 let cold = choose(&units, &check, &Fit::default(), &gate, ALPHA, 1, &had);
2902 assert_eq!((cold.routes, cold.pf[1], cold.source), (fixed_routes(&units, &had), 0.0, "fixed"));
2903 let mut freed = 0;
2904 for seed in 0..20 {
2905 freed += usize::from(choose(&units, &check, &f, &gate, ALPHA, seed, &[]).routes[1] == Route::Free);
2906 let d = choose(&units, &check, &f, &gate, ALPHA, seed, &had);
2907 assert!(d.routes[1] != Route::Free && d.pf[1] == 0.0, "{d:?}");
2908 }
2909 assert!(freed > 0, "the scenario: without `had` the model hands L2 over");
2910 let block = [unit("L2", 3, 0), unit("L1", 1, 1)];
2912 let d = choose(&block, &Check { stop: true, ..check }, &f, &gate, ALPHA, 3, &had);
2913 assert_eq!((d.routes[0], d.p[0]), (Route::Ask, 1.0));
2914 }
2915
2916 #[test]
2919 fn a_check_past_the_solvers_size_takes_the_fixed_rule() {
2920 let f = free_beats_ask();
2921 let gate = open_gate();
2922 let ids: Vec<String> = (0..=SOLVE_MAX).map(|i| format!("L{i}")).collect();
2923 let units: Vec<Unit> = ids.iter().enumerate().map(|(r, l)| unit(l, if r == 0 { 3 } else { 2 }, r)).collect();
2925 let check = Check { session: "n", stop: false, tokens: 0 };
2926 let at = choose(&units[..SOLVE_MAX], &check, &f, &gate, ALPHA, 1, &[]);
2927 assert_eq!(at.source, EVIDENCE, "{at:?}");
2928 assert_ne!(at.routes, fixed_routes(&units[..SOLVE_MAX], &[]), "at the size the model hands them over: {at:?}");
2929 let past = choose(&units, &check, &f, &gate, ALPHA, 1, &[]);
2930 let base = fixed_routes(&units, &[]);
2931 assert_eq!((&past.routes, past.source, past.risk), (&base, "size", 0.0));
2932 assert_eq!(past.p, base.iter().map(|r| f64::from(u8::from(*r == Route::Ask))).collect::<Vec<_>>(), "its choice is its propensity");
2933 }
2934
2935 #[test]
2938 fn a_confounded_history_keeps_the_stop_three_asked() {
2939 let mut sim = Sim::new(12);
2940 for i in 0..12 {
2941 sim.session(&format!("s{i}"), 16, 12, &truth(), &mut fixed_rule);
2942 }
2943 let s = sim.summary();
2944 assert!(s.applies.iter().all(|a| a.grade != 3 || a.route == Route::Ask || a.route == Route::Not), "never free");
2945 let f = fit(&s, None);
2946 let gate = open_gate();
2947 for seed in 0..10 {
2948 let d = choose(&[unit("L1", 3, 0), unit("L2", 2, 1)], &Check { session: "n", stop: true, tokens: 0 }, &f, &gate, ALPHA, seed, &[]);
2949 assert_eq!((d.routes[0], d.p[0]), (Route::Ask, 1.0));
2950 }
2951 }
2952
2953 fn three_stop_world(lift: f64) -> Summary {
2956 let mut t = truth();
2957 t.b[ASK] = 0.0;
2958 t.b[NOT] = 0.0;
2959 t.lift = lift;
2960 let mut cold = |grades: &[u8], stop: bool, rng: &mut Rng| {
2961 let (routes, p, pf, _) = fixed_rule(grades, stop, rng);
2962 (routes, p, pf, "learned")
2963 };
2964 let mut sim = Sim::new(1);
2965 for i in 0..40 {
2966 sim.session(&format!("s{i}"), 20, 12, &t, &mut cold);
2967 }
2968 summary::with_fit(sim.summary(), None)
2969 }
2970
2971 #[test]
2976 fn a_lift_of_the_stops_with_a_three_is_not_taken_for_asking() {
2977 let (lifted, plain) = (three_stop_world(1.5), three_stop_world(0.0));
2978 let units = [unit("A", 2, 0), unit("B", 1, 1)];
2979 assert_eq!(fixed_routes(&units, &[]), [Route::Free; 2], "the fixed rule hands both over");
2980 let asked = |s: &Summary| -> f64 {
2981 assert!(remaining(&s.gate) > 0.0, "the gate opened on the history: {:?}", s.gate.allowance);
2982 let seeds = 20u32;
2983 (0..seeds)
2984 .map(|seed| {
2985 let d = choose(&units, &Check { session: "new", stop: true, tokens: 0 }, &s.fit, &s.gate, ALPHA, u64::from(seed), &[]);
2986 assert_ne!(d.source, "fixed", "{d:?}");
2987 d.p[0]
2988 })
2989 .sum::<f64>()
2990 / f64::from(seeds)
2991 };
2992 let (with_lift, without) = (asked(&lifted), asked(&plain));
2993 let m = lifted.fit.model.as_ref().expect("a model");
2994 let sd = |k: usize| (0..=k).map(|j| m.chol[k * P + j].powi(2)).sum::<f64>().sqrt();
2995 let (ask_stop, twin) = ((m.beta[ASK_STOP], sd(ASK_STOP)), (m.beta[THREE_STOP], sd(THREE_STOP)));
2996 let seen = format!("asked {with_lift:.3} with the lift, {without:.3} without; ask_stop {ask_stop:.3?}, three_stop {twin:.3?}");
2997 assert!(with_lift < 0.9, "a stop without a 3 is not blocked by the lift: {seen}");
2998 assert!((with_lift - without).abs() < 0.15, "{seen}");
2999 assert!(ask_stop.0.abs() < 2.0 * ask_stop.1, "asking at the stop is not identified: {seen}");
3000 assert!(twin.0 > 0.75, "the lift is the twin's: {seen}");
3001 }
3002
3003 #[test]
3008 fn the_gate_never_risks_more_than_its_alpha_share() {
3009 let spec = coefs();
3012 let mut beta: Vec<f64> = spec.iter().map(|c| c.mean).collect();
3013 beta[ASK] = -0.4;
3014 beta[B0] = -0.5;
3015 let mut chol = vec![0.0; P * P];
3016 for k in 0..P {
3017 chol[k * P + k] = if k == ASK { 1.0 } else { 0.3 };
3018 }
3019 let model = Model { beta, chol, sd: spec.iter().map(|c| c.sd).collect(), spread: [0.3; FAMILIES], ..Model::default() };
3020 let rates = Rates { free_stop: 1.0, free_live: 1.0, later: 0.5, wait: 0.0, ask_live: 1.0 };
3021 let f = Fit { key: "test".into(), model: Some(model), rates, ..Fit::default() };
3022 let ids: Vec<String> = (0..10).map(|i| format!("L{i}")).collect();
3023 let mut refused = 0;
3024 let mut run = |gated: bool, rng: &mut Rng| -> (bool, usize) {
3025 let mut gate = Gate::default();
3026 let mut broke = false;
3027 let mut deviated = 0;
3028 for c in 0..200 {
3029 let n = 1 + (rng.next_u64() % 3) as usize;
3030 let units: Vec<Unit> = (0..n).map(|r| unit(&ids[(rng.next_u64() % 10) as usize], 1 + (rng.next_u64() % 3) as u8, r)).collect();
3031 let session = format!("t{}", c / 10);
3032 let check = Check { session: &session, stop: c % 5 == 0, tokens: 0 };
3033 let routes = if c < 10 {
3034 fixed_routes(&units, &[])
3035 } else if gated {
3036 let d = choose(&units, &check, &f, &gate, ALPHA, c as u64, &[]);
3037 match d.source {
3040 GATE => {
3041 assert_eq!((&d.routes, d.risk), (&fixed_routes(&units, &[]), 0.0), "{d:?}");
3042 refused += 1;
3043 }
3044 _ if d.routes != fixed_routes(&units, &[]) => {
3045 assert_eq!(d.source, EXPERIMENT, "{d:?}");
3046 assert!(d.risk <= d.gate + 1e-12, "{d:?}");
3047 }
3048 _ => {}
3049 }
3050 d.routes
3051 } else {
3052 let open = Gate { allowance: f64::INFINITY, ..gate.clone() };
3053 choose(&units, &check, &f, &open, ALPHA, c as u64, &[]).routes
3054 };
3055 deviated += usize::from(routes != fixed_routes(&units, &[]));
3056 ledger_add_without(&f, &mut gate, &check, &units, &routes, &[]);
3057 broke |= lower_quantile(&mut gate.gain.clone()) < -gate.allowance - 1e-9;
3058 }
3059 (broke, deviated)
3060 };
3061 let (broke, deviated) = run(true, &mut Rng::new(77));
3062 assert!(!broke && deviated > 0, "gated: broke {broke}, deviated {deviated}");
3063 let (broke, _) = run(false, &mut Rng::new(77));
3064 assert!(broke, "ungated, the same run overspends: the scenario needs the gate");
3065 assert!(refused > 0, "the gate refused deviations, logged as `gate`");
3066 }
3067
3068 #[test]
3069 fn no_price_without_a_budget_and_spend_within_one() {
3070 let mut s = planted(10, 12, 6);
3071 s.pool.bins = (0..6).map(|b| (400 * (1u64 << b), 200 / (b as u64 + 1))).collect();
3073 s.pool.shown = 40 * 120;
3074 s.pool.tokens = 12_000 * 120;
3075 let free = fit(&s, None);
3076 assert_eq!((free.lambda, free.budget), (0.0, None), "no budget, no price");
3077 let checks = s.sessions.iter().map(|x| f64::from(x.checks)).sum::<f64>() / s.sessions.len() as f64;
3078 let spend = |f: &Fit| {
3079 s.sessions.iter().map(|x| f64::from(x.checks) * pool_tokens(f, true, None, session_seed(&x.id)) as f64).sum::<f64>() / s.sessions.len() as f64
3080 };
3081 assert!(spend(&free) >= checks * crate::whole::TOKENS as f64 - 1e-6, "unpriced, the pool keeps the baseline");
3082 let budget = checks * 9_000.0;
3083 let priced = fit(&s, Some(budget));
3084 assert!(priced.lambda > 0.0, "a budget below today's spend prices tokens");
3085 assert!(spend(&priced) <= budget + 1e-6, "{} > {budget}", spend(&priced));
3086 let roomy = fit(&s, Some(checks * 1e6));
3087 assert_eq!(roomy.lambda, 0.0, "a budget the baseline meets costs nothing");
3088 }
3089
3090 #[test]
3094 fn a_budget_holds_where_pools_may_grow() {
3095 let mut s = Summary::default();
3096 for i in 0..10 {
3097 s.sessions.push(summary::Session { id: format!("s{i}"), checks: 10, ..summary::Session::default() });
3098 }
3099 let pool = PoolFit {
3103 yields: vec![(50.0, 10.0); 8],
3104 per_lesson: 300.0,
3105 latency: (1000.0, 0.1, 500.0),
3106 latency_n: LATENCY_MIN as u32,
3107 not: (-1.0, 0.1),
3108 ..PoolFit::default()
3109 };
3110 let spend = |f: &Fit, limit: Option<u64>| {
3111 s.sessions.iter().map(|x| f64::from(x.checks) * pool_tokens(f, true, limit, session_seed(&x.id)) as f64).sum::<f64>() / s.sessions.len() as f64
3112 };
3113 let baseline = 10.0 * crate::whole::TOKENS as f64;
3114 let unpriced = Fit { pool: pool.clone(), ..Fit::default() };
3115 assert!(spend(&unpriced, Some(60_000)) > 1.5 * baseline, "the scenario: unpriced, the pools grow");
3116 for budget in [baseline, 1.2 * baseline, 0.9 * baseline] {
3117 let mut f = Fit { pool: pool.clone(), budget: Some(budget), ..Fit::default() };
3118 f.lambda = price(&s, &f, budget);
3119 assert!(f.lambda > 0.0, "budget {budget}: growth is priced");
3120 for limit in [None, Some(3_500), Some(60_000)] {
3121 assert!(spend(&f, limit) <= budget + 1e-6, "budget {budget}, limit {limit:?}: spent {} (λ {:e})", spend(&f, limit), f.lambda);
3122 }
3123 }
3124 }
3125
3126 #[test]
3127 fn the_pool_stays_at_the_baseline_without_evidence_or_a_limit() {
3128 let f = Fit::default();
3129 assert_eq!(pool_tokens(&f, true, Some(60_000), 1), crate::whole::TOKENS);
3130 assert_eq!(pool_tokens(&f, false, Some(60_000), 1), crate::whole::TOKENS);
3131 let mut s = planted(3, 6, 7);
3132 s.pool.bins = (0..6).map(|b| (100 * (1u64 << b), 90)).collect();
3133 s.pool.shown = 40 * 18;
3134 s.pool.tokens = 12_000 * 18;
3135 let f = fit(&s, None);
3136 assert!(pool_tokens(&f, true, Some(60_000), 1) <= crate::whole::TOKENS);
3138 assert_eq!(pool_grid(), [6000, 8485, 12000, 16971, 24000]);
3139 }
3140
3141 #[test]
3144 fn the_yields_stop_at_a_bin_no_lesson_was_shown_in() {
3145 let mut s = Summary::default();
3146 s.pool.bins = vec![(40, 20), (80, 30), (0, 0), (3, 1)];
3147 let mut pool = PoolFit::default();
3148 pool_state(&s, &mut pool);
3149 assert_eq!(pool.yields, pava(&[(21.0, 21.0), (31.0, 51.0)]), "the gap ends them");
3150 s.pool.bins[2] = (5, 1);
3151 pool_state(&s, &mut pool);
3152 assert_eq!(pool.yields.len(), 4, "filled, the tail's bin past it counts");
3153 }
3154
3155 #[test]
3156 fn pava_keeps_the_yields_from_rising() {
3157 let y = pava(&[(5.0, 5.0), (8.0, 2.0), (1.0, 9.0), (2.0, 8.0)]);
3158 let means: Vec<f64> = y.iter().map(|(a, b)| a / (a + b)).collect();
3159 assert!(means.windows(2).all(|w| w[0] >= w[1] - 1e-12), "{means:?}");
3160 assert_eq!(y[0], y[1], "the rising pair merged");
3161 assert_eq!(y[0], (13.0, 7.0));
3162 }
3163
3164 #[test]
3165 fn the_samplers_have_their_moments() {
3166 let mut rng = Rng::new(1);
3167 let n = 20_000;
3168 let z: Vec<f64> = (0..n).map(|_| normal(&mut rng)).collect();
3169 let mean = z.iter().sum::<f64>() / n as f64;
3170 let var = z.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / n as f64;
3171 assert!(mean.abs() < 0.03 && (var - 1.0).abs() < 0.05, "{mean} {var}");
3172 for shape in [0.4, 1.0, 3.5] {
3173 let g = (0..n).map(|_| gamma(shape, &mut rng)).sum::<f64>() / n as f64;
3174 assert!((g - shape).abs() < 0.05 * shape.max(1.0), "{shape}: {g}");
3175 }
3176 let b = (0..n).map(|_| beta_draw(2.0, 6.0, &mut rng)).sum::<f64>() / n as f64;
3177 assert!((b - 0.25).abs() < 0.01, "{b}");
3178 }
3179
3180 #[test]
3182 fn the_upper_bound_covers_the_mean_at_every_time() {
3183 let mut rng = Rng::new(21);
3184 let mut missed = 0;
3185 for _ in 0..400 {
3186 let x: Vec<f64> = (0..200).map(|_| normal(&mut rng) * 0.4 - 0.1).collect();
3187 missed += usize::from(upper_bound(&x) < -0.1);
3188 }
3189 assert!(missed <= 400 * 7 / 100, "{missed} of 400 below the mean at some time (δ = 5%)");
3190 assert_eq!(upper_bound(&[-5.0; CS_MIN - 1]), f64::INFINITY, "too few sessions to trip");
3191 let x: Vec<f64> = (0..60).map(|_| normal(&mut rng) * 0.4 - 0.3).collect();
3192 assert!(upper_bound(&x) < 0.0, "60 sessions of a clear harm trip it: {}", upper_bound(&x));
3193 }
3194
3195 fn monitored(free_helps: bool) -> Summary {
3199 let mut t = truth();
3200 t.b[ASK] = 1.5;
3201 if free_helps {
3202 t.b[ASK] = 0.0;
3203 }
3204 let mut sim = Sim::new(31);
3205 let mut half = |grades: &[u8], _: bool, rng: &mut Rng| -> (Vec<Route>, Vec<f64>, Vec<f64>, &'static str) {
3206 let routes = grades.iter().map(|_| if rng.unit() < 0.5 { Route::Ask } else { Route::Free }).collect();
3207 (routes, vec![0.5; grades.len()], vec![0.5; grades.len()], "learned")
3208 };
3209 let mut sim_session = |sim: &mut Sim, name: &str| {
3211 for _ in 0..6 {
3212 let grades = vec![3, 2];
3213 let ids = vec![format!("L{}", sim.rng.next_u64() % 8), format!("L{}", sim.rng.next_u64() % 8)];
3214 sim.check(name, false, &grades, &ids, &t, &mut half);
3215 }
3216 };
3217 for i in 0..60 {
3218 sim_session(&mut sim, &format!("s{i}"));
3219 }
3220 summary::with_fit(sim.summary(), None)
3221 }
3222
3223 #[test]
3224 fn the_monitor_reverts_a_harmful_history_and_keeps_a_benign_one() {
3225 let harmful = monitored(false);
3226 assert!(harmful.monitor.reverted && harmful.monitor.ucb.is_some_and(|u| u < 0.0), "{:?}", harmful.monitor);
3227 assert_eq!(harmful.monitor.sessions, 60);
3228 let m = &harmful.monitor;
3231 assert!(m.labelled > 0 && m.labelled <= m.units, "{m:?}");
3232 assert!(m.ess > 0.0 && m.ess <= f64::from(m.labelled) + 1e-9 && m.max_weight >= 1.0 && m.truncated <= m.labelled, "{m:?}");
3233 let benign = monitored(true);
3234 assert!(!benign.monitor.reverted, "{:?}", benign.monitor);
3235 assert!(benign.monitor.mean > harmful.monitor.mean);
3236 }
3237
3238 #[test]
3241 fn the_fixed_routes_are_the_fixed_rule_by_position() {
3242 assert!(fixed_routes(&[], &[]).is_empty());
3243 let block = [unit("A", 1, 0), unit("B", 3, 1), unit("C", 2, 2)];
3244 assert_eq!(fixed_routes(&block, &[]), [Route::Ask; 3], "a 3 blocks: every 3, 2 and 1 up to the cap");
3245 let soft = [unit("A", 1, 0), unit("B", 2, 1)];
3246 assert_eq!(fixed_routes(&soft, &[]), [Route::Free; 2], "no 3: handed over free");
3247 let many: Vec<String> = (0..12).map(|i| format!("L{i}")).collect();
3248 let over: Vec<Unit> = many.iter().enumerate().map(|(r, l)| unit(l, 1, r)).collect();
3249 let routes = fixed_routes(&over, &[]);
3250 assert_eq!(routes.iter().filter(|r| **r == Route::Free).count(), conscience::CAP);
3251 assert_eq!(&routes[conscience::CAP..], &[Route::Not; 2], "past the cap: not now");
3252 let twice: Vec<Unit> = (0..11).map(|r| unit(if r < 9 { many[r].as_str() } else { "Z" }, 1, r)).collect();
3254 assert_eq!(&fixed_routes(&twice, &[])[9..], &[Route::Free, Route::Not]);
3255 }
3256
3257 #[test]
3258 fn a_prior_fit_writes_no_key() {
3259 assert_eq!(serde_json::to_string(&Fit::default()).unwrap(), "{}");
3260 assert_eq!(serde_json::from_str::<Fit>("{}").unwrap(), Fit::default());
3261 }
3262
3263 #[test]
3267 fn the_evidence_rule_is_chooses_without_the_test() {
3268 let f = free_beats_ask();
3269 let open = Gate { gain: vec![0.0; DRAWS], allowance: f64::INFINITY, ..Gate::default() };
3270 let ids: Vec<String> = (0..6).map(|i| format!("L{i}")).collect();
3271 let mut left = 0;
3272 for (n, g) in [(3, [2, 1, 1]), (3, [3, 2, 1]), (2, [1, 1, 0])] {
3273 let units: Vec<Unit> = (0..n).map(|r| unit(&ids[r], g[r], r)).collect();
3274 for stop in [false, true] {
3275 let check = Check { session: "new", stop, tokens: 0 };
3276 let routes = decide_on_evidence(&units, &check, &f, &[]);
3277 for seed in 0..5 {
3278 let d = choose(&units, &check, &f, &open, 0.0, seed, &[]);
3279 assert_eq!((&routes, d.xk), (&d.routes, None), "{g:?} stop {stop} seed {seed}");
3280 }
3281 left += usize::from(routes != fixed_routes(&units, &[]));
3282 assert_eq!(decide_on_evidence(&units, &check, &Fit::default(), &[]), fixed_routes(&units, &[]), "no model: the fixed rule");
3283 }
3284 }
3285 assert!(left > 0, "the scenario leaves the fixed rule");
3286 }
3287
3288 #[test]
3293 fn under_the_prior_the_evidence_rule_is_the_fixed_rule_on_every_grade_table() {
3294 let p = prior();
3295 let open = Gate { gain: vec![0.0; DRAWS], allowance: f64::INFINITY, ..Gate::default() };
3296 let ids: Vec<String> = (0..8).map(|i| format!("L{i}")).collect();
3297 let mut grades = tables(5);
3298 grades.extend([vec![2; 8], vec![1; 8], [vec![3], vec![2; 7]].concat()]);
3299 for g in &grades {
3300 let units: Vec<Unit> = ids.iter().zip(g).enumerate().map(|(r, (l, g))| unit(l, *g, r)).collect();
3301 for had in [Vec::new(), vec!["L1".to_string()]] {
3302 let want = fixed_routes(&units, &had);
3303 for stop in [true, false] {
3304 let check = Check { session: "s", stop, tokens: 0 };
3305 assert_eq!(decide_on_evidence(&units, &check, &p, &had), want, "{g:?} stop {stop} had {had:?}");
3306 let d = choose(&units, &check, &p, &open, 0.0, 7, &had);
3307 assert_eq!((&d.routes, d.source), (&want, if units.is_empty() { FIXED } else { EVIDENCE }), "{g:?} stop {stop} had {had:?}");
3308 }
3309 }
3310 }
3311 }
3312
3313 #[test]
3319 fn a_fit_with_evidence_of_dilution_deviates() {
3320 let world = |dilution: f64, sd: f64| {
3321 let mut beta = vec![0.0; P];
3322 beta[ASK] = 0.2;
3323 beta[NOT] = -1.0;
3324 beta[ASK_M] = dilution;
3325 let mut chol = vec![0.0; P * P];
3326 for k in 0..P {
3327 chol[k * P + k] = if k == ASK_M { sd } else { 0.01 };
3328 }
3329 let model = Model { beta, chol, sd: vec![1.0; P], spread: [0.01; FAMILIES], ..Model::default() };
3330 let rates = Rates { free_stop: 1.0, free_live: 1.0, later: 0.5, wait: 0.0, ask_live: 1.0 };
3331 Fit { key: "dilution".into(), model: Some(model), rates, ..Fit::default() }
3332 };
3333 let ids: Vec<String> = (0..9).map(|i| format!("L{i}")).collect();
3334 let units: Vec<Unit> = ids.iter().enumerate().map(|(r, l)| unit(l, if r == 0 { 3 } else { 2 }, r)).collect();
3335 let base = fixed_routes(&units, &[]);
3336 assert_eq!(base, [Route::Ask; 9], "the scenario: the fixed rule asks all nine");
3337 let gate = open_gate();
3338 for stop in [false, true] {
3339 let check = Check { session: "n", stop, tokens: 0 };
3340 let d = choose(&units, &check, &world(-1.0, 0.01), &gate, ALPHA, 1, &[]);
3341 assert_eq!((d.source, d.xk, d.xs), (EVIDENCE, None, None), "{d:?}");
3342 assert_ne!(d.routes, base, "{d:?}");
3343 assert!(d.routes.contains(&Route::Free) && d.routes.contains(&Route::Ask), "asks few, hands over the rest: {d:?}");
3344 assert_eq!(d.p, d.routes.iter().map(|r| f64::from(u8::from(*r == Route::Ask))).collect::<Vec<_>>(), "exact");
3345 assert_eq!(d.pf, d.routes.iter().map(|r| f64::from(u8::from(*r == Route::Free))).collect::<Vec<_>>(), "exact");
3346 let wide = world(-0.2, 0.5);
3348 let m = wide.model.as_ref().expect("a model");
3349 let feats = Feat::of(&units, stop, 0);
3350 let t = theta_mean(m, &feats, &units, &Levels::of(check.session, &m.judge, &units));
3351 let at_mode = solve(&Table::of(&Valuer { t: &t, feats: &feats, rates: &wide.rates }, &[false; 9]), crate::attention::RAIL, conscience::CAP);
3352 assert_ne!(at_mode, base, "stop {stop}");
3353 for f in [wide, world(0.0, 0.01)] {
3354 assert_eq!(decide_on_evidence(&units, &check, &f, &[]), base, "stop {stop}");
3355 }
3356 }
3357 }
3358
3359 fn tested_at(units: &[Unit], had: &[String], seeds: u64) -> Vec<Decision> {
3362 let p = prior();
3363 let open = Gate { gain: vec![0.0; DRAWS], allowance: 1e9, ..Gate::default() };
3364 (0..seeds).map(|seed| choose(units, &Check { session: "s", stop: false, tokens: 0 }, &p, &open, ALPHA, seed, had)).collect()
3365 }
3366
3367 #[test]
3372 fn the_tests_count_is_uniform_and_its_propensities_exact() {
3373 let ids: Vec<String> = (0..5).map(|i| format!("L{i}")).collect();
3374 let units: Vec<Unit> = ids.iter().zip([1, 2, 1, 2, 1]).enumerate().map(|(r, (l, g))| unit(l, g, r)).collect();
3375 let order = [1, 3, 0, 2, 4];
3376 let n = 1200;
3377 let mut counts = [0usize; 6];
3378 for d in tested_at(&units, &[], n) {
3379 let (Some(5), Some(s)) = (d.xk, d.xs) else { panic!("{d:?}") };
3380 let s = s as usize;
3381 counts[s] += 1;
3382 for (i, &j) in order.iter().enumerate() {
3383 assert_eq!(d.routes[j], if i < s { Route::Ask } else { Route::Free }, "{d:?}");
3384 assert_eq!((d.p[j], d.pf[j]), (outcome::round_p((5 - i) as f64 / 6.0), outcome::round_p((i + 1) as f64 / 6.0)), "{d:?}");
3385 }
3386 assert_eq!(d.source, if s == 0 { EXPERIMENT_CONTROL } else { EXPERIMENT }, "{d:?}");
3387 }
3388 let q = 1.0 / 6.0;
3389 let se = (q * (1.0 - q) / n as f64).sqrt();
3390 for (s, c) in counts.iter().enumerate() {
3391 let freq = *c as f64 / n as f64;
3392 assert!((freq - q).abs() <= 4.0 * se, "s = {s}: {freq:.4} against {q:.4} ({counts:?})");
3393 }
3394 }
3395
3396 #[test]
3399 fn the_test_never_touches_the_stop_or_a_three() {
3400 let ids: Vec<String> = (0..3).map(|i| format!("L{i}")).collect();
3401 let units: Vec<Unit> = ids.iter().zip([2, 3, 1]).enumerate().map(|(r, (l, g))| unit(l, g, r)).collect();
3402 let base = fixed_routes(&units, &[]);
3403 let p = prior();
3404 let open = Gate { gain: vec![0.0; DRAWS], allowance: 1e9, ..Gate::default() };
3405 for seed in 0..40 {
3406 let d = choose(&units, &Check { session: "s", stop: true, tokens: 0 }, &p, &open, ALPHA, seed, &[]);
3407 assert_eq!((&d.routes, d.source, d.xk, d.xs), (&base, EVIDENCE, None, None), "{d:?}");
3408 }
3409 let mut seen = std::collections::BTreeSet::new();
3410 for d in tested_at(&units, &[], 60) {
3411 assert_eq!((d.routes[1], d.p[1], d.pf[1], d.xk), (Route::Ask, 1.0, 0.0, Some(2)), "{d:?}");
3412 let s = d.xs.unwrap_or(9);
3413 assert_eq!(d.source == EXPERIMENT_CONTROL, s == 2, "{d:?}");
3414 seen.insert(s);
3415 }
3416 assert_eq!(seen.len(), 3, "every count drawn: {seen:?}");
3417 }
3418
3419 #[test]
3422 fn the_test_never_hands_over_a_lesson_the_session_had() {
3423 let ids: Vec<String> = (0..4).map(|i| format!("L{i}")).collect();
3424 let had = ["L2".to_string()];
3425 let soft: Vec<Unit> = ids.iter().zip([2, 1, 2, 1]).enumerate().map(|(r, (l, g))| unit(l, g, r)).collect();
3426 let block: Vec<Unit> = ids.iter().zip([2, 3, 2, 1]).enumerate().map(|(r, (l, g))| unit(l, g, r)).collect();
3427 for (units, kept) in [(&soft, Route::Not), (&block, Route::Ask)] {
3428 let mut moved = 0;
3429 for d in tested_at(units, &had, 60) {
3430 assert_eq!((d.routes[2], d.pf[2], d.xk), (kept, 0.0, Some(if kept == Route::Not { 3 } else { 2 })), "{d:?}");
3431 moved += usize::from(d.source == EXPERIMENT);
3432 }
3433 assert!(moved > 0, "the scenario: the test ran");
3434 }
3435 }
3436
3437 #[test]
3441 fn the_gate_refuses_the_test_when_short_and_risk_zero_turns_it_off() {
3442 let ids: Vec<String> = (0..3).map(|i| format!("L{i}")).collect();
3443 let units: Vec<Unit> = ids.iter().zip([2, 1, 1]).enumerate().map(|(r, (l, g))| unit(l, g, r)).collect();
3444 let base = fixed_routes(&units, &[]);
3445 let exact = |r: Route| base.iter().map(|x| f64::from(u8::from(*x == r))).collect::<Vec<_>>();
3446 let p = prior();
3447 let short = Gate { gain: vec![0.0; DRAWS], allowance: 1e-12, ..Gate::default() };
3448 let check = Check { session: "s", stop: false, tokens: 0 };
3449 let (mut refused, mut control) = (0, 0);
3450 for seed in 0..40 {
3451 let d = choose(&units, &check, &p, &short, ALPHA, seed, &[]);
3452 assert_eq!((&d.routes, &d.p, &d.pf, d.risk, d.xk), (&base, &exact(Route::Ask), &exact(Route::Free), 0.0, Some(3)), "{d:?}");
3453 match d.xs {
3454 Some(0) => control += usize::from(d.source == EXPERIMENT_CONTROL),
3455 _ => refused += usize::from(d.source == GATE),
3456 }
3457 let off = choose(&units, &check, &p, &open_gate(), 0.0, seed, &[]);
3458 assert_eq!((&off.routes, off.source, off.xk, off.xs), (&base, EVIDENCE, None, None), "{off:?}");
3459 }
3460 assert!(refused > 0 && control > 0 && refused + control == 40, "refused {refused}, control {control}");
3461 let room = Gate { allowance: 1e9 * ALPHA, ..short };
3464 assert!((0..40).any(|seed| choose(&units, &check, &p, &room, ALPHA, seed, &[]).source == EXPERIMENT));
3465 assert!((0..40).all(|seed| choose(&units, &check, &p, &room, 1e-12, seed, &[]).source != EXPERIMENT), "α scales it down");
3466 }
3467
3468 #[test]
3473 fn risk_zero_refuses_an_evidence_deviation_the_allowance_alone_would_cover() {
3474 let f = free_beats_ask();
3475 let ids: Vec<String> = (0..4).map(|i| format!("L{i}")).collect();
3477 let units: Vec<Unit> = ids.iter().enumerate().map(|(r, l)| unit(l, if r == 0 { 3 } else { 2 }, r)).collect();
3478 let base = fixed_routes(&units, &[]);
3479 let check = Check { session: "n", stop: false, tokens: 0 };
3480 let loss = 2.0 * units.len() as f64;
3483 let gate = Gate { gain: vec![-loss; DRAWS], allowance: 2.0 * loss, ..Gate::default() };
3484 let at = choose(&units, &check, &f, &gate, ALPHA, 1, &[]);
3485 assert_eq!(at.source, EVIDENCE, "the scenario: the allowance admits the deviation at ALPHA: {at:?}");
3486 assert_ne!(at.routes, base, "{at:?}");
3487 let off = choose(&units, &check, &f, &gate, 0.0, 1, &[]);
3488 assert_eq!((&off.routes, off.source, off.risk, off.gate, off.xk), (&base, GATE, 0.0, 0.0, None), "{off:?}");
3489 }
3490
3491 #[test]
3494 fn the_dilution_prior_has_its_mode_at_zero() {
3495 let c = coefs();
3496 for k in [ASK_M, FREE_M] {
3497 assert_eq!((c[k].mean, c[k].sd, c[k].nonpos, c[k].ard), (0.0, 0.1, true, false), "{}", NAMES[k]);
3498 }
3499 let p = prior();
3500 let m = p.model.as_ref().expect("the prior is a model");
3501 let units = [unit("A", 2, 0), unit("B", 2, 1), unit("C", 2, 2)];
3502 let feats = Feat::of(&units, false, 0);
3503 let t = theta_mean(m, &feats, &units, &Levels::of("s", &m.judge, &units));
3504 let v = Valuer { t: &t, feats: &feats, rates: &p.rates };
3505 assert!(v.ask(0, 1) > 0.0 && (v.ask(0, 1) - v.ask(0, 3)).abs() < 1e-12, "{} vs {}", v.ask(0, 1), v.ask(0, 3));
3506 assert!((v.free(0, 1) - v.free(0, 3)).abs() < 1e-12);
3507 }
3508
3509 #[test]
3513 fn a_bounded_slope_is_drawn_from_its_truncated_prior() {
3514 let p = prior();
3515 let m = p.model.as_ref().expect("the prior is a model");
3516 const N: u64 = 4000;
3517 let draws: Vec<[f64; P]> = (0..N).map(|k| draw_beta(m, gate_seed(&p, k as usize))).collect();
3518 let half = -DILUTION_PRIOR.1 * (2.0 / std::f64::consts::PI).sqrt();
3519 for k in [ASK_M, FREE_M] {
3521 assert!(draws.iter().all(|b| b[k] < 0.0), "{}: a draw at 0", NAMES[k]);
3522 let mean = draws.iter().map(|b| b[k]).sum::<f64>() / N as f64;
3523 assert!((mean - half).abs() < 0.005, "{}: mean {mean} against {half}", NAMES[k]);
3524 }
3525 let ask = draws.iter().map(|b| b[ASK]).sum::<f64>() / N as f64;
3526 assert!((ask - ASK_PRIOR.0).abs() < 0.03, "ask {ask}");
3527 for z in [-2.0, -0.3, 0.0, 1.0, 2.5] {
3529 assert!((below_zero(-1.0, 0.1, z) - (-1.0 + 0.1 * z)).abs() < 1e-6, "z {z}");
3530 let x = below_zero(50.0, 0.1, z);
3531 assert!(x <= 0.0 && x > -0.01, "z {z}: {x}");
3532 }
3533 assert!(below_zero(0.0, 0.1, -1.0) < below_zero(0.0, 0.1, 1.0), "monotone in z");
3534 assert!((phi(1.959_964) - 0.975).abs() < 1e-6 && (probit(0.975f64.ln(), 0.025) - 1.959_964).abs() < 1e-5);
3535 }
3536
3537 #[test]
3540 fn the_prior_is_a_model_of_the_priors_alone() {
3541 let p = prior();
3542 let m = p.model.as_ref().expect("the prior is a model");
3543 assert!(m.current() && m.converged && p.key.is_empty() && p.labels == 0);
3544 for (k, c) in coefs().iter().enumerate() {
3545 assert_eq!(m.beta[k], if c.nonpos { c.mean.min(0.0) } else { c.mean }, "{}", NAMES[k]);
3546 assert!((m.chol[k * P + k] - c.sd).abs() < 1e-12, "{}", NAMES[k]);
3547 }
3548 assert_eq!(m.spread, [LEVEL_START; FAMILIES]);
3549 assert_eq!(p.rates, Rates { free_stop: 0.5, free_live: 0.5, later: 0.5, wait: 0.5, ask_live: 1.0 });
3550 assert!(fit(&Summary::default(), None).model.is_none());
3551 }
3552
3553 #[test]
3556 fn a_prediction_is_the_posterior_chance_of_a_change() {
3557 let f = Feat { grade: 3, rank: 1, stop: true, three: true, pool: 0.0 };
3558 let mut beta = vec![0.0; P];
3559 beta[B0] = 0.8;
3560 beta[G3] = 0.4;
3561 beta[ASK] = 0.3;
3562 let lessons = BTreeMap::from([("A".to_string(), (0.5, 0.0))]);
3563 let certain = Model { beta: beta.clone(), chol: vec![0.0; P * P], sd: vec![1.0; P], spread: [0.0; FAMILIES], lessons, ..Model::default() };
3564 let mut fit = Fit { model: Some(certain.clone()), ..Fit::default() };
3565 let eta = dot(&std::array::from_fn(|k| beta[k]), &design(&f, Route::Ask, 2)) + 0.5;
3566 assert!((predict(&fit, &f, Route::Ask, 2, "s", "A").unwrap() - sigmoid(eta)).abs() < 1e-12);
3567 let mut wide = certain.clone();
3568 for k in 0..P {
3569 wide.chol[k * P + k] = 1.0;
3570 }
3571 wide.spread = [0.5; FAMILIES];
3572 fit.model = Some(wide);
3573 let p = predict(&fit, &f, Route::Ask, 2, "s", "unseen").unwrap();
3574 assert!(p > 0.5 && p < sigmoid(eta - 0.5), "{p} against {}", sigmoid(eta - 0.5));
3575 fit.model = Some(Model { beta: vec![0.0; P - 1], ..certain });
3576 assert_eq!(predict(&fit, &f, Route::Ask, 2, "s", "A"), None);
3577 }
3578}