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

probl_engine/report/
results.rs
