probl_engine/report/
results.rs

1//! What each report says, as numbers. The renderer formats these, and the
2//! `probl` library gives them to programs that embed Probl, so the two can't
3//! disagree (docs/library.md).
4
5use super::{Acc, Format, Sink, summary_quantile};
6use crate::continuous::Mixture;
7use crate::value::Value;
8use crate::weight::Weight;
9use probl_sema::ir::{ReportKind, ReportSite};
10
11/// One report's results.
12#[derive(Clone, Debug)]
13pub struct ReportResult {
14    /// The share of the worlds' weight (or of the runs) that reached the
15    /// report. `None` when every visit counts, when nothing reached it, when
16    /// no world finished, or when worlds failed before evidence they'd have
17    /// met, so their weight can't be compared.
18    pub reach: Option<Reach>,
19    /// One group per `by` key, in key order; one keyed `()` without `by`.
20    pub groups: Vec<GroupResult>,
21}
22
23#[derive(Clone, Copy, Debug, PartialEq)]
24pub struct Reach {
25    pub share: f64,
26    /// When sampling: its standard error.
27    pub se: Option<f64>,
28}
29
30/// What one report saw for one key.
31#[derive(Clone, Debug)]
32pub struct GroupResult {
33    pub key: Value,
34    /// When every reported value was a fact: the chance that it's true.
35    pub fact: Option<Quantity>,
36    /// Unless a value is continuous: each value's probability, in the order
37    /// of the values (facts as `true` and `false`, for tables that mix them
38    /// with other values).
39    pub values: Option<Vec<(Value, Quantity)>>,
40    /// The values and their probabilities, among the resolved worlds, in the
41    /// order of the values (facts count as `true` and `false`).
42    pub distribution: Vec<(Value, f64)>,
43    /// For real numbers, and continuous values: their summaries.
44    pub numeric: Option<Numeric>,
45    /// The share of the group's weight that isn't resolved, as the renderer
46    /// counts it.
47    pub unresolved_share: f64,
48    /// When sampling: how many runs it rests on.
49    pub support: Option<Support>,
50}
51
52/// How many sampled runs a group rests on.
53#[derive(Clone, Copy, Debug, PartialEq)]
54pub struct Support {
55    /// Independent runs that reached it; several visits in one run count once.
56    pub contributing_runs: u64,
57    /// The effective sample size of its denominator.
58    pub effective: f64,
59}
60
61/// One reported quantity: a probability, a mean, a quantile.
62#[derive(Clone, Copy, Debug, PartialEq)]
63pub struct Quantity {
64    /// Over the resolved worlds (or runs).
65    pub point: Option<f64>,
66    /// No unresolved weight, or missing mass, affects it.
67    pub complete: bool,
68    /// When it isn't complete and enumerating: bounds on the quantity over
69    /// all the weight, as if the unresolved weight had gone either way.
70    pub bounds: Option<(f64, f64)>,
71    /// When sampling: its Monte Carlo error.
72    pub sampling: Option<Uncertainty>,
73}
74
75#[derive(Clone, Copy, Debug, PartialEq)]
76pub struct Uncertainty {
77    pub status: Status,
78    /// When `Estimated`, and zero when `IntegratedZero`.
79    pub se: Option<f64>,
80    /// A 95% Wilson interval, where its rules allow one.
81    pub wilson: Option<(f64, f64)>,
82    pub support: Support,
83}
84
85/// What's known about a sampled quantity's Monte Carlo error.
86#[derive(Clone, Copy, Debug, PartialEq, Eq)]
87pub enum Status {
88    /// Its standard error is estimated.
89    Estimated,
90    /// The runs can't establish a useful estimate: they all agreed, without
91    /// integrated outcomes.
92    NotEstimable,
93    /// The engine doesn't compute one for this statistic.
94    NotComputed,
95    /// No empirical error: every run integrated the same outcomes.
96    IntegratedZero,
97}
98
99/// The summaries of a group of real numbers, or of continuous values.
100#[derive(Clone, Debug)]
101pub struct Numeric {
102    pub mean: Quantity,
103    pub sd: Quantity,
104    /// Every value was a probability (shown as percentages).
105    pub percent: bool,
106    shape: Shape,
107    /// For the quantiles: the group's completeness and sampling support.
108    complete: bool,
109    support: Option<Support>,
110}
111
112#[derive(Clone, Debug)]
113enum Shape {
114    /// A continuous marginal, possibly mixed with point masses.
115    Mixture(Mixture),
116    /// Point masses, as values (integers keep their exact value for printing).
117    Points(Vec<(Value, f64)>),
118}
119
120impl Numeric {
121    /// The continuous marginal, when there's one.
122    pub fn mixture(&self) -> Option<&Mixture> {
123        match &self.shape {
124            Shape::Mixture(m) => Some(m),
125            Shape::Points(_) => None,
126        }
127    }
128
129    /// The `q` quantile, by the report convention: the median is the
130    /// midpoint of the middle values, the others select an outcome. `q` must
131    /// be in [0, 1].
132    pub fn quantile(&self, q: f64) -> Option<Quantity> {
133        let point = match &self.shape {
134            Shape::Mixture(m) => Some(if q == 0.5 { m.median() } else { m.quantile(q) }),
135            Shape::Points(points) => summary_quantile(points, q).and_then(|v| v.as_f64()),
136        }?;
137        Some(Quantity {
138            point: Some(point),
139            complete: self.complete,
140            bounds: None,
141            sampling: self.support.map(not_computed),
142        })
143    }
144
145    /// The `q` quantile of point masses as a value, for printing integers
146    /// exactly.
147    pub fn quantile_value(&self, q: f64) -> Option<Value> {
148        match &self.shape {
149            Shape::Points(points) => summary_quantile(points, q),
150            Shape::Mixture(_) => None,
151        }
152    }
153
154    /// The values and their probabilities, as `f64`.
155    pub(super) fn points(&self) -> Option<Vec<(f64, f64)>> {
156        match &self.shape {
157            Shape::Points(points) => Some(
158                points
159                    .iter()
160                    .map(|(v, p)| (v.as_f64().unwrap_or(f64::NAN), *p))
161                    .collect(),
162            ),
163            Shape::Mixture(_) => None,
164        }
165    }
166}
167
168fn not_computed(support: Support) -> Uncertainty {
169    Uncertainty {
170        status: Status::NotComputed,
171        se: None,
172        wilson: None,
173        support,
174    }
175}
176
177/// The results of every report, in source order. `unresolved` is the weight
178/// the run left unresolved, which decides what's complete; `format` is what
179/// the renderer uses.
180pub fn results(sites: &[ReportSite], sinks: &[Sink], format: Format, unresolved: Weight) -> Vec<ReportResult> {
181    sites
182        .iter()
183        .zip(sinks)
184        .map(|(site, sink)| {
185            // A simple report's facts are judged as reached once per world;
186            // a table's by how its site is reached.
187            let kind = if site.key_label.is_none() {
188                ReportKind::Once
189            } else {
190                site.kind
191            };
192            ReportResult {
193                reach: reach(site.kind, sink, format),
194                groups: sink
195                    .groups
196                    .iter()
197                    .map(|(key, acc)| group(key, acc, format, kind, unresolved))
198                    .collect(),
199            }
200        })
201        .collect()
202}
203
204fn reach(kind: ReportKind, sink: &Sink, format: Format) -> Option<Reach> {
205    if kind == ReportKind::PerVisit || sink.groups.is_empty() || format.program_total.is_zero() || !format.reach_known {
206        return None;
207    }
208    if let Some(all_squares) = format.run_squares {
209        // The share of the runs' weight that reached the report, counting
210        // each run once, and its standard error (section 14).
211        let total = format.program_total;
212        let p = sink.reached.ratio(total);
213        let squares = sink.reached_squares;
214        let spread = squares.scale((1.0 - p) * (1.0 - p)) + all_squares.saturating_sub(squares).scale(p * p);
215        let se = spread.ratio(total * total).sqrt();
216        return Some(Reach { share: p, se: Some(se) });
217    }
218    Some(Reach {
219        share: sink.reach().ratio(format.program_total),
220        se: None,
221    })
222}
223
224pub(super) fn group(key: &Value, acc: &Acc, format: Format, kind: ReportKind, unresolved: Weight) -> GroupResult {
225    let sampled = acc.sampled();
226    let support = sampled.then(|| Support {
227        contributing_runs: acc.contributing_runs(),
228        effective: acc.effective(),
229    });
230    let complete = (unresolved + acc.missing).is_zero();
231    let distribution = acc.distribution();
232    let fact = acc
233        .is_event()
234        .then(|| fact(acc, format, kind, complete, unresolved, support));
235    let continuous = distribution
236        .iter()
237        .any(|(v, _)| matches!(v, Value::Analytic(_) | Value::Continuous(_)));
238    let values = (!continuous).then(|| values(acc, &distribution, complete, unresolved, support));
239    let numeric = if acc.is_event() {
240        None
241    } else {
242        numeric(acc, &distribution, complete, support)
243    };
244    GroupResult {
245        key: key.clone(),
246        fact,
247        values,
248        numeric,
249        unresolved_share: acc.unresolved_share(format.unresolved),
250        distribution,
251        support,
252    }
253}
254
255/// The chance that a reported fact is true.
256fn fact(
257    acc: &Acc,
258    format: Format,
259    kind: ReportKind,
260    complete: bool,
261    unresolved: Weight,
262    support: Option<Support>,
263) -> Quantity {
264    let p = acc.chance();
265    let sampling = support.map(|support| {
266        let se = acc.chance_se();
267        // A Wilson interval needs ordinary, equally weighted observations, and
268        // is given when the standard error isn't useful: every run agreed, or
269        // too few of them contributed.
270        let wilson = acc
271            .chance_interval95(kind)
272            .filter(|_| !format.weighted)
273            .filter(|_| p == 0.0 || p == 1.0 || support.contributing_runs < 30);
274        let integrated = acc.moments.as_ref().is_some_and(|m| m.integrated);
275        let (status, se) = if se != 0.0 {
276            (Status::Estimated, Some(se))
277        } else if integrated {
278            (Status::IntegratedZero, Some(0.0))
279        } else {
280            (Status::NotEstimable, None)
281        };
282        Uncertainty {
283            status,
284            se,
285            wilson,
286            support,
287        }
288    });
289    Quantity {
290        point: Some(p),
291        complete,
292        bounds: (!complete && support.is_none()).then(|| acc.chance_bounds(unresolved)),
293        sampling,
294    }
295}
296
297/// Each reported value's probability, with bounds over the unresolved
298/// weight: as if it had all gone to that value, or none of it.
299fn values(
300    acc: &Acc,
301    distribution: &[(Value, f64)],
302    complete: bool,
303    unresolved: Weight,
304    support: Option<Support>,
305) -> Vec<(Value, Quantity)> {
306    let u = unresolved + acc.missing;
307    let resolved = Weight::sum(acc.values.values().copied()) + acc.facts;
308    distribution
309        .iter()
310        .map(|(value, p)| {
311            let bounds = (!complete && support.is_none()).then(|| {
312                // Facts are counted apart from the other values.
313                let w = match value {
314                    Value::Bool(true) => acc.yes,
315                    Value::Bool(false) => acc.facts.saturating_sub(acc.yes),
316                    _ => acc.values.get(value).copied().unwrap_or_default(),
317                };
318                let total = resolved + u;
319                (w.ratio(total), (w + u).ratio(total).min(1.0))
320            });
321            let sampling = support.map(|support| {
322                let se = acc.value_se(value, *p);
323                let (status, se) = if se == 0.0 {
324                    (Status::NotEstimable, None)
325                } else {
326                    (Status::Estimated, Some(se))
327                };
328                Uncertainty {
329                    status,
330                    se,
331                    wilson: None,
332                    support,
333                }
334            });
335            let quantity = Quantity {
336                point: Some(*p),
337                complete,
338                bounds,
339                sampling,
340            };
341            (value.clone(), quantity)
342        })
343        .collect()
344}
345
346/// The mean and standard deviation of real numbers, or of a continuous
347/// marginal mixed with point masses, as the renderer prints them.
348fn numeric(acc: &Acc, distribution: &[(Value, f64)], complete: bool, support: Option<Support>) -> Option<Numeric> {
349    let summary = |mean: f64, sd: f64, mean_sampling: Option<Uncertainty>, percent: bool, shape: Shape| {
350        let quantity = |x: f64, sampling| Quantity {
351            point: Some(x),
352            complete,
353            bounds: None,
354            sampling,
355        };
356        Some(Numeric {
357            mean: quantity(mean, mean_sampling),
358            sd: quantity(sd, support.map(not_computed)),
359            percent,
360            shape,
361            complete,
362            support,
363        })
364    };
365    if let Some(m) = super::analytic_mixture(distribution) {
366        let (mean, sd) = (m.mean(), m.sd());
367        return summary(mean, sd, support.map(not_computed), false, Shape::Mixture(m));
368    }
369    let nums: Vec<(f64, f64)> = distribution
370        .iter()
371        .map(|(v, p)| match v {
372            Value::Int(n) => n
373                .to_f64()
374                .filter(|x| n.cmp_f64(*x).is_some_and(|c| c.is_eq()))
375                .map(|x| (x, *p)),
376            Value::Float(x) | Value::Prob(x) => Some((*x, *p)),
377            _ => None,
378        })
379        .collect::<Option<_>>()?;
380    let percent = distribution.iter().all(|(v, _)| matches!(v, Value::Prob(_)));
381    let total: f64 = nums.iter().map(|(_, p)| p).sum();
382    let mean = nums.iter().map(|(x, p)| x * p).sum::<f64>() / total;
383    let sd = (nums.iter().map(|(x, p)| (x - mean).powi(2) * p).sum::<f64>() / total).sqrt();
384    if !mean.is_finite() || !sd.is_finite() {
385        return None;
386    }
387    let mean_sampling = support.map(|support| {
388        let se = acc.mean_se().1;
389        let (status, se) = if se > 0.0 {
390            (Status::Estimated, Some(se))
391        } else {
392            (Status::NotEstimable, None)
393        };
394        Uncertainty {
395            status,
396            se,
397            wilson: None,
398            support,
399        }
400    });
401    summary(mean, sd, mean_sampling, percent, Shape::Points(distribution.to_vec()))
402}
probl_engine/report/results.rs

probl_engine/report/
results.rs
