hyperpaths_rs/

hyperpath.rs

//! Implementation of the Spiess-Florian algorithm for transit assignment.
//! See the ref. at spiess_floarian.tex LaTeX file.

use std::collections::{HashMap, HashSet};
use std::sync::atomic::{AtomicBool, Ordering};

use crate::hyperpath_queue::PriorityQueue;
use crate::transit_network::Link;

/// Strategy is the optimal strategy as defined in the Spiess-Florian algorithm.
pub struct Strategy<'a> {
    /// u_{i} - expected travel time from node i to destination
    pub labels: HashMap<String, f64>,
    /// f_{i} - combined frequency of attractive links at node i.
    /// `f64::INFINITY` marks a node whose basket is a single no-wait link.
    pub freqs: HashMap<String, f64>,
    /// \overline{A} - attractive links forming the hyperpath
    pub a_set: Vec<&'a Link>,
}

/// The waiting-time constant of the Spiess-Florian expected travel time:
///   u_i = (1 + sum(f_a * (c_a + u_j))) / f_i
/// When the first attractive link arrives at a node the sum is empty and
/// the numerator starts from this constant.
pub(crate) const ALPHA: f64 = 1.0;

pub static VERBOSE: AtomicBool = AtomicBool::new(false);

pub(crate) fn verbose() -> bool {
    VERBOSE.load(Ordering::Relaxed)
}

pub fn find_optimal_strategy<'a>(
    all_links: &'a [Link],
    all_stops: &HashSet<String>,
    destination: &str,
) -> Strategy<'a> {
    /* 1.1 Initialization */
    if verbose() {
        println!("1.1 Initialization \\\\");
    }
    let mut u: HashMap<String, f64> = HashMap::with_capacity(all_stops.len());
    let mut f: HashMap<String, f64> = HashMap::with_capacity(all_stops.len());
    for stop in all_stops {
        if verbose() {
            println!("$f_{{{}}} = 0$ \\\\ ", stop);
        }
        f.insert(stop.clone(), 0.0);
        if stop == destination {
            if verbose() {
                println!("$u_{{{}}} = 0$ \\\\ ", destination);
            }
            u.insert(stop.clone(), 0.0);
            continue;
        }
        if verbose() {
            println!("$u_{{{}}} = Infinity$ \\\\ ", stop);
        }
        u.insert(stop.clone(), f64::INFINITY);
    }

    let mut overline_a: Vec<Option<&'a Link>> = Vec::with_capacity(all_links.len() / 2);
    // Positions of each node's basket links inside overline_a, so that a
    // no-wait link can replace the whole basket. Replaced entries are set
    // to None and compacted at the end.
    let mut a_set_idx: HashMap<&'a str, Vec<usize>> = HashMap::new();

    let mut links_by_to_node: HashMap<&'a str, Vec<&'a Link>> = HashMap::new();
    for link in all_links {
        links_by_to_node
            .entry(link.to_node.as_str())
            .or_default()
            .push(link);
    }

    let mut entries: HashMap<&'a str, Vec<usize>> = HashMap::with_capacity(all_links.len());
    let mut pq = PriorityQueue::with_capacity(all_links.len());
    for link in all_links {
        let priority = u.get(&link.to_node).copied().unwrap_or(0.0) + link.travel_cost;
        let id = pq.push(link, priority);
        entries.entry(link.from_node.as_str()).or_default().push(id);
    }
    pq.init();
    if verbose() {
        pq.print();
    }
    while pq.len() > 0 {
        /* 1.2 Get next link */
        if verbose() {
            pq.print();
        }
        let entry_id = match pq.pop() {
            Some(id) => id,
            None => break,
        };
        let priority = pq.priority(entry_id);
        if priority.is_infinite() && priority > 0.0 {
            break;
        }
        let a = pq.link(entry_id);
        let i = a.from_node.as_str();
        let j = a.to_node.as_str();
        let sum_uc = u.get(j).copied().unwrap_or(0.0) + a.travel_cost;

        /* 1.3 Update node label */
        if verbose() {
            println!("Process: $a = (i, j) = ({}, {})$, \\\\ ", i, j);
        }
        // A node already served by a no-wait link is final: the no-wait
        // link absorbs all flow (its share f_a/f_i is 1 in the limit),
        // so no other link may enter the basket
        let f_i = f.get(i).copied().unwrap_or(0.0);
        if f_i.is_infinite() {
            continue;
        }
        let u_i = u.get(i).copied().unwrap_or(0.0);
        // Strict improvement test: a link is accepted only if it
        // strictly improves the label. Step 1.3 of Spiess & Florian
        // (1989) prints the nonstrict u_i >= u_j + c_a, but the two
        // rules differ only at exact equality, where the update is a
        // no-op (the combination formula returns u_i unchanged; for
        // f_a = inf the basket is replaced at the same value): labels,
        // expected travel times and every number published in the
        // paper are identical either way. The strict form is what
        // part 2 needs. Step 2.2 loads links "in reverse topological
        // order (decreasing u_j + c_a)" (p. 94) and Proposition 4
        // claims flow conservation "by construction" - both presume an
        // acyclic strategy, which the nonstrict rule does not
        // guarantee: in an expanded route graph a boarding link (cost
        // 0) into a route node whose label came from its own alighting
        // link (cost 0) has key exactly u_i, so >= admits a zero-cost
        // stop -> node -> stop cycle and the one-pass loading strands
        // the volume entering it (see
        // test_board_alight_loop_conservation). Rejecting at equality
        // keeps the strategy acyclic and stays optimal: for the
        // rejected link mu_a = 0 satisfies dual feasibility (20) as an
        // equality and complementary slackness (24) holds since
        // v_a = 0, a degenerate optimum. The prose of p. 94 ("if this
        // time is smaller than u_i, link a is included") describes
        // exactly this strict rule. All step, equation and page
        // references above are to the original paper, not to the
        // spiess_floarian.tex excerpt in this repository.
        if u_i <= sum_uc {
            continue;
        }
        if verbose() {
            println!(
                "\\quad $u_i \\leq u_j + c_a : {} \\leq {}$ - FALSE \\\\ ",
                u_i, sum_uc
            );
        }
        if a.headway <= 0.0 {
            // No-wait link (infinite frequency): the modified step 1.3
            // given by the paper on p. 96 - the exact limit of the label
            // update formula as f_a -> inf. The link replaces the whole
            // attractive basket:
            //   u_i := u_j + c_a, f_i := inf, A_i := {a}
            u.insert(i.to_string(), sum_uc);
            f.insert(i.to_string(), f64::INFINITY);
            let indices = a_set_idx.entry(i).or_default();
            for idx in indices.iter() {
                overline_a[*idx] = None;
            }
            indices.clear();
            overline_a.push(Some(a));
            indices.push(overline_a.len() - 1);
            if verbose() {
                println!(
                    "\\quad no-wait link: $u_i = u_j + c_a = {}$, $f_i = \\infty$, basket replaced by $({}, {})$ \\\\ ",
                    sum_uc, i, j
                );
            }
        } else {
            let freq = 1.0 / a.headway;
            if verbose() {
                println!("\\quad $f_a = {}$ \\\\ ", freq);
                println!("\\quad $u_j + c_a = {}$ \\\\ ", sum_uc);
                println!("\\quad $u_i = {}$ \\\\ ", u_i);
            }
            let new_u = if f_i == 0.0 {
                // First link in the basket: u_i = (1 + f_a*(u_j+c_a)) / f_a
                (ALPHA + freq * sum_uc) / freq
            } else {
                (f_i * u_i + freq * sum_uc) / (f_i + freq)
            };
            u.insert(i.to_string(), new_u);
            f.insert(i.to_string(), f_i + freq);
            overline_a.push(Some(a));
            a_set_idx.entry(i).or_default().push(overline_a.len() - 1);
            if verbose() {
                println!(
                    "\\quad$u_i = \\frac{{f_i * u_i + f_a * (u_j + c_a)}}{{f_i + f_a}} = {}$, $f_i = {}$ \\\\ ",
                    new_u,
                    f_i + freq
                );
                println!(
                    "\\quad $\\overline{{A}} = \\overline{{A}} \\cup {{({}, {})}}$ \\\\ ",
                    i, j
                );
            }
        }

        if let Some(links_to_update) = links_by_to_node.get(i) {
            for link in links_to_update {
                if let Some(i_entries) = entries.get(link.from_node.as_str()) {
                    for &eid in i_entries {
                        let entry_link = pq.link(eid);
                        if entry_link.to_node == i && entry_link.from_node == link.from_node {
                            let new_priority = u.get(i).copied().unwrap_or(0.0) + link.travel_cost;
                            pq.update(eid, new_priority);
                            break;
                        }
                    }
                }
            }
        }
        if verbose() {
            println!("Node labels: \\\\");
            for s in all_stops {
                println!(
                    "${} -> (u_i, f_i) = ({}, {})$ \\\\ ",
                    s,
                    u.get(s).copied().unwrap_or(0.0),
                    f.get(s).copied().unwrap_or(0.0)
                );
            }
        }
    }

    // Compact the attractive set: drop entries replaced by no-wait links.
    // The append order is preserved, i.e. non-decreasing u_j + c_a.
    let a_set: Vec<&'a Link> = overline_a.into_iter().flatten().collect();

    Strategy {
        labels: u,
        freqs: f,
        a_set,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_hyper_paths() {
        VERBOSE.store(true, Ordering::Relaxed);
        let all_nodes: HashSet<String> = ["A", "X", "X2", "Y", "Y3", "B"]
            .iter()
            .map(|s| s.to_string())
            .collect();
        let all_links = vec![
            Link::new("A", "B", "Line 1", 25.0, 6.0),
            Link::new("A", "X2", "Line 2", 7.0, 6.0),
            Link::new("X2", "X", "Line 2", 0.0, 0.0),
            Link::new("X", "X2", "Line 2", 0.0, 6.0),
            Link::new("X2", "Y", "Line 2", 6.0, 0.0),
            Link::new("Y3", "Y", "Line 3", 0.0, 15.0),
            Link::new("Y", "B", "Line 4", 10.0, 3.0),
            Link::new("X", "Y3", "Line 3", 4.0, 15.0),
            Link::new("Y", "Y3", "Line 3", 0.0, 15.0),
            Link::new("Y3", "B", "Line 3", 4.0, 0.0),
        ];
        let destination_node = "B";
        let ops = find_optimal_strategy(&all_links, &all_nodes, destination_node);

        const EPS: f64 = 1e-9;

        // With exact no-wait handling the labels match the paper exactly:
        // no big-M artifacts like 4.000000000000001
        let expected_labels: HashMap<&str, f64> = HashMap::from([
            ("A", 27.75),
            ("X", 19.071428571428573),
            ("X2", 17.5),
            ("Y", 11.5),
            ("Y3", 4.0),
            ("B", 0.0),
        ]);
        // +Inf marks nodes whose basket is a single no-wait link
        let expected_freqs: HashMap<&str, f64> = HashMap::from([
            ("A", 1.0 / 3.0),
            ("X", 7.0 / 30.0),
            ("X2", f64::INFINITY),
            ("Y", 0.4),
            ("Y3", f64::INFINITY),
            ("B", 0.0),
        ]);
        // Matches the paper order (Spiess & Florian 1989, p. 93-94)
        let expected_a_set: Vec<&Link> = vec![
            // Y3->B
            &all_links[9],
            // Y->Y3
            &all_links[8],
            // X->Y3
            &all_links[7],
            // Y->B
            &all_links[6],
            // X2->Y
            &all_links[4],
            // X->X2
            &all_links[3],
            // A->X2
            &all_links[1],
            // A->B
            &all_links[0],
        ];

        assert_eq!(
            ops.labels.len(),
            expected_labels.len(),
            "Incorrect number of labels"
        );
        assert_eq!(
            ops.freqs.len(),
            expected_freqs.len(),
            "Incorrect number of frequencies"
        );
        assert_eq!(
            ops.a_set.len(),
            expected_a_set.len(),
            "Incorrect number of links in attractive set"
        );

        for (k, v) in &ops.labels {
            assert!(
                expected_labels.contains_key(k.as_str()),
                "Incorrect label key {} has met",
                k
            );
            let want = expected_labels[k.as_str()];
            assert!(
                (v - want).abs() <= EPS,
                "Incorrect label value for node {}: got {}, want {}",
                k,
                v,
                want
            );
        }
        for (k, v) in &ops.freqs {
            assert!(
                expected_freqs.contains_key(k.as_str()),
                "Incorrect frequency key {} has met",
                k
            );
            let want = expected_freqs[k.as_str()];
            if want.is_infinite() {
                assert!(
                    v.is_infinite() && *v > 0.0,
                    "Frequency for node {} must be +Inf, got {}",
                    k,
                    v
                );
            } else {
                assert!(
                    (v - want).abs() <= EPS,
                    "Incorrect frequency value for node {}: got {}, want {}",
                    k,
                    v,
                    want
                );
            }
        }
        for (i, v) in ops.a_set.iter().enumerate() {
            println!("{:?} {:?}", v, expected_a_set[i]);
            assert!(
                std::ptr::eq(*v, expected_a_set[i]),
                "Incorrect link in attractive set at index {}",
                i
            );
        }
    }

    #[test]
    fn test_no_wait_replaces_basket() {
        // A boarding link enters the basket of I first (key 4), then a
        // cheaper no-wait chain I->W->D (key 5 < current u_I = 10) must
        // replace it entirely: exact label, infinite frequency, single link.
        let all_nodes: HashSet<String> = ["I", "W", "D"].iter().map(|s| s.to_string()).collect();
        let all_links = vec![
            // boarding link, key u_D + 4 = 4, accepted first: u_I = 6 + 4 = 10
            Link::new("I", "D", "Bus", 4.0, 6.0),
            // no-wait walk, key u_W + 3 = 5, replaces the basket: u_I = 5
            Link::new("I", "W", "Walk", 3.0, 0.0),
            // no-wait walk, key 2
            Link::new("W", "D", "Walk", 2.0, 0.0),
        ];
        let ops = find_optimal_strategy(&all_links, &all_nodes, "D");

        assert!((ops.labels["I"] - 5.0).abs() <= 1e-12);
        assert!((ops.labels["W"] - 2.0).abs() <= 1e-12);
        assert!(ops.freqs["I"].is_infinite());
        assert!(ops.freqs["W"].is_infinite());

        // The replaced boarding link I->D must not remain attractive
        assert_eq!(
            ops.a_set.len(),
            2,
            "basket of I must hold only the no-wait link"
        );
        for link in &ops.a_set {
            assert_eq!(
                link.headway, 0.0,
                "only no-wait links expected in the attractive set"
            );
        }
    }
}