Skip to main content

dpm_schedule/
pert.rs

1//! Classical analytic PERT approximation of the project finish along one mean critical path.
2//!
3//! PERT treats the finish as normal, with the mean and variance of a single critical path of the
4//! deterministic projection. Any near-critical path that overtakes it in some outcomes is ignored,
5//! so where paths merge the approximation usually overstates the probability of finishing by a
6//! horizon. It is reported beside the sampled probability, never in place of it.
7
8use crate::network::{EPSILON, Network, Times};
9use crate::placement::Placement;
10use crate::remaining_duration::RemainingDuration;
11use crate::{RemainingOptions, ScheduleError};
12use dpm_model::{DependencyId, Endpoint, Plan, ThreePointEstimate, WorkItemId};
13use serde::{Deserialize, Serialize};
14use std::collections::{BTreeMap, BTreeSet};
15
16/// What the approximation leaves out, stated wherever its probability is reported.
17pub const PERT_APPROXIMATION_NOTE: &str = "approximation: classical PERT normal over one mean critical path; it ignores near-critical paths and therefore usually overstates the probability of finishing by the horizon";
18
19/// Classical PERT variance `((P − O) / 6)²` of a three-point estimate, in hours squared.
20///
21/// It is not the variance of the beta-PERT distribution the simulation samples, whose variance
22/// depends on where the most-likely value lies.
23#[must_use]
24pub fn pert_variance_hours2(estimate: ThreePointEstimate) -> f64 {
25    let spread = (estimate.pessimistic_hours - estimate.optimistic_hours).max(0.0) / 6.0;
26    spread * spread
27}
28
29#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
30/// Normal approximation of the project finish along one mean critical path.
31pub struct PertApproximation {
32    /// Deterministic project finish in elapsed hours; without calendars it is the sum of the
33    /// path's PERT means and lags.
34    pub mean_hours: f64,
35    /// Sum of the classical PERT variances of the durations the path's finish depends on. On a
36    /// calendar these are working hours, not the elapsed hours of `mean_hours`, so the spread of
37    /// the elapsed finish is understated and the approximation is more confident still.
38    pub variance_hours2: f64,
39    /// Square root of `variance_hours2`.
40    pub sigma_hours: f64,
41    /// The mean critical path from its first activity to the one that sets the project finish.
42    pub path: Vec<WorkItemId>,
43    /// Classical PERT variance of every projected activity's remaining duration: zero for
44    /// finished or exactly known work, and for started work the range its elapsed time allows.
45    pub activity_variance_hours2: BTreeMap<WorkItemId, f64>,
46}
47
48impl PertApproximation {
49    /// `Φ((H − mean) / σ)`; a path without variance finishes by `H` exactly when its mean does.
50    #[must_use]
51    pub fn probability_by(&self, horizon_hours: f64) -> f64 {
52        if self.sigma_hours <= EPSILON {
53            return if horizon_hours + EPSILON >= self.mean_hours {
54                1.0
55            } else {
56                0.0
57            };
58        }
59        let z = (horizon_hours - self.mean_hours) / self.sigma_hours;
60        0.5 * (1.0 + libm::erf(z / std::f64::consts::SQRT_2))
61    }
62}
63
64/// The approximation of the full baseline, every task at its whole PERT duration.
65pub fn pert_approximation(plan: &Plan) -> Result<PertApproximation, ScheduleError> {
66    plan.validate()?;
67    let durations = plan
68        .work_items
69        .iter()
70        .map(|(id, work)| Ok((*id, RemainingDuration::of(*id, work, false, 0.0)?)))
71        .collect::<Result<Vec<_>, ScheduleError>>()?;
72    let network = Network::compile(plan)?;
73    approximate(&network, &durations, None)
74}
75
76/// The approximation of the remaining graph that
77/// [`deterministic_remaining_with`](crate::deterministic_remaining_with) projects.
78///
79/// Started work contributes its remaining mean and the variance of the range still possible;
80/// finished work contributes neither. Excluded work is on no path and has no variance.
81pub fn pert_approximation_remaining_with(
82    plan: &Plan,
83    timeline: &dpm_model::Timeline,
84    options: RemainingOptions,
85) -> Result<PertApproximation, ScheduleError> {
86    options.validate()?;
87    plan.validate()?;
88    let remaining = crate::cpm::remaining_plan_at(plan, timeline)?;
89    let durations = remaining.durations()?;
90    let network = Network::compile(&remaining.plan)?;
91    let expected: Vec<f64> = durations.iter().map(|(_, d)| d.expected()).collect();
92    let mut placement = remaining.placement(&network, timeline, &expected, options)?;
93    let mut approximation = approximate(&network, &durations, placement.as_mut())?;
94    approximation
95        .activity_variance_hours2
96        .retain(|id, _| !remaining.excluded.contains(id));
97    approximation
98        .path
99        .retain(|id| !remaining.excluded.contains(id));
100    Ok(approximation)
101}
102
103/// [`pert_approximation_remaining_with`] without adjustments.
104pub fn pert_approximation_remaining_at(
105    plan: &Plan,
106    timeline: &dpm_model::Timeline,
107) -> Result<PertApproximation, ScheduleError> {
108    pert_approximation_remaining_with(plan, timeline, RemainingOptions::default())
109}
110
111fn approximate(
112    network: &Network,
113    durations_by_id: &[(WorkItemId, RemainingDuration)],
114    mut placement: Option<&mut Placement>,
115) -> Result<PertApproximation, ScheduleError> {
116    let count = network.order().len();
117    let (mut means, mut variances) = (vec![0.0; count], vec![0.0; count]);
118    for (id, source) in durations_by_id {
119        let at = network
120            .position(id)
121            .ok_or(ScheduleError::MissingWorkItem(*id))?;
122        if let (Some(mean), Some(variance)) = (means.get_mut(at), variances.get_mut(at)) {
123            *mean = source.expected();
124            *variance = source.pert_variance();
125        }
126    }
127    let times = match placement.as_deref_mut() {
128        Some(placement) => network.times_placed(&means, placement)?,
129        None => network.times(&means)?,
130    };
131    let critical = (0..count)
132        .map(|at| Ok(network.total_float(&times, at)? <= EPSILON))
133        .collect::<Result<Vec<bool>, ScheduleError>>()?;
134    let critical_at = |at: usize| critical.get(at).copied().unwrap_or(false);
135    let critical_links: Vec<_> = network
136        .relation_slacks(&means, &times, placement.as_deref())?
137        .into_iter()
138        .filter(|(_, from, to, slack)| *slack <= EPSILON && critical_at(*from) && critical_at(*to))
139        .collect();
140    let links = Links {
141        ids: critical_links.iter().map(|(id, ..)| *id).collect(),
142        continued: critical_links.iter().map(|(_, from, ..)| *from).collect(),
143    };
144    let (path, variance) = mean_critical_path(network, &times, &means, &variances, &links)?;
145    Ok(PertApproximation {
146        mean_hours: times.finish,
147        variance_hours2: variance,
148        sigma_hours: variance.sqrt(),
149        path,
150        activity_variance_hours2: network.order().iter().copied().zip(variances).collect(),
151    })
152}
153
154/// Critical driving constraints, and the positions a critical path continues from.
155struct Links {
156    ids: BTreeSet<DependencyId>,
157    continued: BTreeSet<usize>,
158}
159
160/// Variance of an activity's start and finish events along the best path reaching it, and the
161/// predecessor that path arrives from.
162#[derive(Debug, Clone, Copy)]
163struct Reach {
164    start: f64,
165    finish: f64,
166    from: Option<usize>,
167}
168
169/// One deterministic mean critical path and the variance of its finish.
170///
171/// Each activity is reached over critical driving constraints; where several reach it, the one
172/// that gives its finish the larger variance wins, ties going to the smaller predecessor id. A
173/// constraint fixes one event of its successor from one event of its predecessor, and the other
174/// event of the successor lies one duration away, so every duration on the path counts once.
175/// The path ends at a critical activity that sets the project finish and from which no critical
176/// constraint continues, such as a closing milestone, chosen the same way.
177fn mean_critical_path(
178    network: &Network,
179    times: &Times,
180    means: &[f64],
181    variances: &[f64],
182    links: &Links,
183) -> Result<(Vec<WorkItemId>, f64), ScheduleError> {
184    let order = network.order();
185    let id_at = |at: usize| order.get(at).copied();
186    let mut reach: Vec<Reach> = Vec::with_capacity(order.len());
187    for at in 0..order.len() {
188        let own = variances.get(at).copied().unwrap_or(0.0);
189        let mut best: Option<(f64, usize, Endpoint, f64)> = None;
190        for edge in network
191            .incoming(at)
192            .iter()
193            .filter(|e| links.ids.contains(&e.id))
194        {
195            let Some(before) = reach.get(edge.other) else {
196                continue;
197            };
198            let carried = match edge.kind.predecessor_endpoint() {
199                Endpoint::Start => before.start,
200                Endpoint::Finish => before.finish,
201            };
202            // Links are compared by the variance they give this activity's finish: a link fixing
203            // its start adds its own duration, one fixing its finish does not.
204            let produced = match edge.kind.successor_endpoint() {
205                Endpoint::Start => carried + own,
206                Endpoint::Finish => carried,
207            };
208            let better = best.is_none_or(|(_, other, _, finish)| {
209                prefers((produced, id_at(edge.other)), (finish, id_at(other)))
210            });
211            if better {
212                best = Some((
213                    carried,
214                    edge.other,
215                    edge.kind.successor_endpoint(),
216                    produced,
217                ));
218            }
219        }
220        reach.push(match best {
221            None => Reach {
222                start: 0.0,
223                finish: own,
224                from: None,
225            },
226            Some((carried, from, Endpoint::Start, _)) => Reach {
227                start: carried,
228                finish: carried + own,
229                from: Some(from),
230            },
231            Some((carried, from, Endpoint::Finish, _)) => Reach {
232                start: carried + own,
233                finish: carried,
234                from: Some(from),
235            },
236        });
237    }
238    let mut end: Option<(f64, usize)> = None;
239    for (at, reached) in reach.iter().enumerate() {
240        let finish = times.finishes(network, means, at)?.0;
241        let critical = network.total_float(times, at)? <= EPSILON;
242        if critical
243            && !links.continued.contains(&at)
244            && finish + EPSILON >= times.finish
245            && end.is_none_or(|(variance, other)| {
246                prefers((reached.finish, id_at(at)), (variance, id_at(other)))
247            })
248        {
249            end = Some((reached.finish, at));
250        }
251    }
252    let Some((variance, mut at)) = end else {
253        return Ok((Vec::new(), 0.0));
254    };
255    let mut path = Vec::new();
256    // Each step moves to an earlier topological position, so the walk ends within the network.
257    while let Some(id) = id_at(at) {
258        path.push(id);
259        match reach.get(at).and_then(|r| r.from) {
260            Some(from) => at = from,
261            None => break,
262        }
263    }
264    path.reverse();
265    Ok((path, variance))
266}
267
268/// Whether a branch with `candidate`'s variance and work id wins over `current`'s: the larger
269/// variance, then the smaller id, so the choice never depends on traversal order.
270fn prefers(candidate: (f64, Option<WorkItemId>), current: (f64, Option<WorkItemId>)) -> bool {
271    match candidate.0.total_cmp(&current.0) {
272        std::cmp::Ordering::Greater => true,
273        std::cmp::Ordering::Equal => candidate.1 < current.1,
274        std::cmp::Ordering::Less => false,
275    }
276}
277
278#[cfg(test)]
279mod tests;