1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
//! This is an alternative pipeline for generating a Scenario, starting from origin-destination
//! data (also called desire lines), which gives a count of commuters between two zones, breaking
//! down by mode.
//!
//! Maybe someday, we'll merge the two approaches, and make the first generate DesireLines as an
//! intermediate step.

use std::collections::HashMap;

use rand::seq::SliceRandom;
use rand::Rng;
use rand_xorshift::XorShiftRng;

use abstutil::Timer;
use geom::{Duration, Polygon, Time};
use map_model::{BuildingID, BuildingType, Map};
use sim::{IndividTrip, MapBorders, PersonSpec, TripEndpoint, TripMode, TripPurpose};

/// This describes some number of commuters living in some named zone, working in another (or the
/// same zone), and commuting using some mode.
#[derive(Debug)]
pub struct DesireLine {
    pub home_zone: String,
    pub work_zone: String,
    pub mode: TripMode,
    pub number_commuters: usize,
}

// TODO Percentage of taking a lunch trip, when to do it, how far to venture out, what mode to
// use...
pub struct Options {
    /// When should somebody depart from home to work?
    pub departure_time: NormalDistribution,
    /// How long should somebody work before returning home?
    pub work_duration: NormalDistribution,
}

impl Options {
    pub fn default() -> Options {
        Options {
            departure_time: NormalDistribution::new(
                Duration::hours(8) + Duration::minutes(30),
                Duration::minutes(30),
            ),
            work_duration: NormalDistribution::new(Duration::hours(9), Duration::hours(1)),
        }
    }
}

/// Generates a scenario from aggregated origin/destination data (DesireLines). The input describes
/// an exact number of people, who live in one zone and work in another (possibly the same) and
/// commute using some mode. For each of them, we just need to pick a specific home and workplace
/// from the zones, and use the Options to pick departure times. We'll wind up creating people who
/// just take two trips daily: home -> work -> home.
///
/// The home and workplace may be a specific building, or they're snapped to a map border,
/// resulting in trips that begin and/or end off-map. The amount of the zone that overlaps with the
/// map boundary determines this. If the zone and map boundary overlap 50% by area, then half of
/// the people to/from this zone will pick buildings, and half will pick borders.
///
/// Currently, zones with no overlap with the map at all are totally filtered out. We could adjust
/// this in the future to create more highway through-traffic.
pub fn disaggregate(
    map: &Map,
    zones: HashMap<String, Polygon>,
    desire_lines: Vec<DesireLine>,
    opts: Options,
    rng: &mut XorShiftRng,
    timer: &mut Timer,
) -> Vec<PersonSpec> {
    // First decide which zones are relevant for our map. Match homes, shops, and border
    // intersections to each zone.
    timer.start("match zones");
    let zones = create_zones(map, zones);
    timer.stop("match zones");

    let mut people = Vec::new();
    timer.start("create people");
    for desire in desire_lines {
        // Skip if we filtered out either zone.
        if !zones.contains_key(&desire.home_zone) || !zones.contains_key(&desire.work_zone) {
            continue;
        }
        let home_zone = &zones[&desire.home_zone];
        let work_zone = &zones[&desire.work_zone];

        for _ in 0..desire.number_commuters {
            // Pick a specific home and workplace. It might be off-map, depending on how much the
            // zone overlaps the map.
            if let (Some((leave_home, goto_home)), Some((_, goto_work))) = (
                home_zone.pick_home(desire.mode, map, rng),
                work_zone.pick_workplace(desire.mode, map, rng),
            ) {
                // Create their schedule
                let goto_work_time = Time::START_OF_DAY + opts.departure_time.sample(rng);
                let return_home_time = goto_work_time + opts.work_duration.sample(rng);
                people.push(PersonSpec {
                    orig_id: None,
                    origin: leave_home,
                    trips: vec![
                        IndividTrip::new(goto_work_time, TripPurpose::Work, goto_work, desire.mode),
                        IndividTrip::new(
                            return_home_time,
                            TripPurpose::Home,
                            goto_home,
                            desire.mode,
                        ),
                    ],
                });
            }
        }
    }
    timer.stop("create people");
    people
}

struct Zone {
    polygon: Polygon,
    pct_overlap: f64,
    // For each building, have a value describing how many people live or work there. The exact
    // value doesn't matter; it's just a relative weighting. This way, we can use a weighted sample
    // and match more people to larger homes/stores.
    homes: Vec<(BuildingID, usize)>,
    workplaces: Vec<(BuildingID, usize)>,
    borders: MapBorders,
}

fn create_zones(map: &Map, input: HashMap<String, Polygon>) -> HashMap<String, Zone> {
    let all_borders = MapBorders::new(map);
    let mut zones = HashMap::new();
    for (name, polygon) in input {
        let mut overlapping_area = 0.0;
        for p in polygon.intersection(map.get_boundary_polygon()) {
            overlapping_area += p.area();
        }
        // Sometimes this is slightly over 100%, because funky things happen with the polygon
        // intersection.
        let pct_overlap = (overlapping_area / polygon.area()).min(1.0);

        // If the zone doesn't intersect our map at all, totally skip it.
        if pct_overlap == 0.0 {
            continue;
        }
        zones.insert(
            name,
            Zone {
                polygon,
                pct_overlap,
                homes: Vec::new(),
                workplaces: Vec::new(),
                borders: all_borders.clone(),
            },
        );
    }

    // Match all buildings to a zone.
    for b in map.all_buildings() {
        let center = b.polygon.center();
        // We're assuming zones don't overlap each other, so just look for the first match.
        if let Some((_, zone)) = zones
            .iter_mut()
            .find(|(_, z)| z.polygon.contains_pt(center))
        {
            match b.bldg_type {
                // The current heuristics for num_residents sometimes assign 0 people to a
                // building. We never want that, so just scale them all up.
                BuildingType::Residential { num_residents, .. } => {
                    zone.homes.push((b.id, num_residents + 1));
                }
                BuildingType::ResidentialCommercial(num_residents, _) => {
                    zone.homes.push((b.id, num_residents + 1));
                    // We know how many different stores are located in each building, according to
                    // OSM. A big mall might have 10 amenities, while standalone
                    // shops just have 1.
                    zone.workplaces.push((b.id, b.amenities.len()));
                }
                BuildingType::Commercial(_) => {
                    zone.workplaces.push((b.id, b.amenities.len()));
                }
                BuildingType::Empty => {}
            }
        }
    }

    // Match border intersections to a zone.
    for zone in zones.values_mut() {
        let polygon = zone.polygon.clone();
        for list in vec![
            &mut zone.borders.incoming_walking,
            &mut zone.borders.incoming_driving,
            &mut zone.borders.incoming_biking,
            &mut zone.borders.outgoing_walking,
            &mut zone.borders.outgoing_driving,
            &mut zone.borders.outgoing_biking,
        ] {
            list.retain(|(i, _)| polygon.contains_pt(map.get_i(*i).polygon.center()));
        }
    }

    zones
}

impl Zone {
    /// Returns endpoints to (leave home, goto home). These're usually the same, except in some
    /// cases of border trips using divided one-ways.
    fn pick_home(
        &self,
        mode: TripMode,
        map: &Map,
        rng: &mut XorShiftRng,
    ) -> Option<(TripEndpoint, TripEndpoint)> {
        if rng.gen_bool(self.pct_overlap) && !self.homes.is_empty() {
            let b = self.homes.choose_weighted(rng, |(_, n)| *n).unwrap().0;
            return Some((TripEndpoint::Bldg(b), TripEndpoint::Bldg(b)));
        }
        self.pick_borders(mode, map, rng)
    }

    /// Returns endpoints to (leave work, goto work). These're usually the same, except in some
    /// cases of border trips using divided one-ways.
    fn pick_workplace(
        &self,
        mode: TripMode,
        map: &Map,
        rng: &mut XorShiftRng,
    ) -> Option<(TripEndpoint, TripEndpoint)> {
        if rng.gen_bool(self.pct_overlap) && !self.workplaces.is_empty() {
            let b = self.workplaces.choose_weighted(rng, |(_, n)| *n).unwrap().0;
            return Some((TripEndpoint::Bldg(b), TripEndpoint::Bldg(b)));
        }
        self.pick_borders(mode, map, rng)
    }

    fn pick_borders(
        &self,
        mode: TripMode,
        map: &Map,
        rng: &mut XorShiftRng,
    ) -> Option<(TripEndpoint, TripEndpoint)> {
        let (incoming, outgoing) = self.borders.for_mode(mode);
        let leave_i = incoming.choose(rng)?.0;
        // If we can use the same border on the way back, prefer that.
        if outgoing.iter().any(|(i, _)| *i == leave_i) {
            return Some((TripEndpoint::Border(leave_i), TripEndpoint::Border(leave_i)));
        }
        // Otherwise, we might have to use a separate border to re-enter. Prefer the one closest to
        // the first, to have a better chance of matching up divided one-ways.
        let leave_pt = map.get_i(leave_i).polygon.center();
        let goto_i = outgoing
            .iter()
            .min_by_key(|(i, _)| map.get_i(*i).polygon.center().dist_to(leave_pt))?
            .0;
        Some((TripEndpoint::Border(leave_i), TripEndpoint::Border(goto_i)))
    }
}

/// A normal distribution of Durations.
pub struct NormalDistribution {
    pub mean: Duration,
    pub std_deviation: Duration,
}

impl NormalDistribution {
    pub fn new(mean: Duration, std_deviation: Duration) -> NormalDistribution {
        NormalDistribution {
            mean,
            std_deviation,
        }
    }

    pub fn sample(&self, rng: &mut XorShiftRng) -> Duration {
        use rand_distr::{Distribution, Normal};

        Duration::seconds(
            Normal::new(
                self.mean.inner_seconds(),
                self.std_deviation.inner_seconds(),
            )
            .unwrap()
            .sample(rng),
        )
    }
}