ras_commander_fast_map/
whole_simulation.rs

1//! Whole-simulation stored-map products: arrival time, duration and percent time
2//! inundated (`native-whole-simulation-horizontal/1`).
3//!
4//! Every stored output time is rendered by the unchanged horizontal method
5//! (cell-constant WSE, native float32 visibility) and accumulated per pixel in bounded memory: first wet time, last wet
6//! time, wet count and the number of separate wet periods.
7//!
8//! A pixel is wet at an output time when its rendered depth is strictly greater
9//! than `threshold` **and** its 2D cell passes RAS Mapper's thresholded cell gate
10//! `wse > cell_minimum + threshold` (float32). Physical 2D cells only; storage
11//! areas are not rendered (Commander refuses such plans). With times `t` as float32 hours and `n` output times:
12//!
13//! - percent = float32(100 * wet_count / n);
14//! - arrival = t[first wet];
15//! - duration = t[end] - t[first wet] in float32, where `end` is the first dry
16//!   output time after the last wet one, or the last wet one when that dry
17//!   time would be the final output time or there is none;
18//! - NoData (NaN) where the pixel is never wet.
19//!
20//! Qualified against native HEC-RAS 7.0.1 (threshold 0.5 ft) on Davis p02 and
21//! BaldEagleCrkMulti2D bridges p05: every valid pixel, arrival and percent value
22//! bit for bit, including initially wet, multi-period and threshold-equal
23//! samples; RAS Mapper's horizontal and sloping products are identical. Duration
24//! is unqualified (some native recession durations are one interval shorter)
25//! and is delivered only with `allow_unqualified`, as are irregular intervals.
26use crate::{MapError, horizontal, sink};
27use rayon::prelude::*;
28use serde::Deserialize;
29use serde_json::{Value, json};
30use std::{fs, path::Path, time::Instant};
31
32pub const SCHEMA: &str = "ras-commander-whole-simulation-request/1";
33pub const METHOD_ID: &str = "native-whole-simulation-horizontal/1";
34pub const PRODUCTS: [&str; 3] = ["arrival_time", "duration", "percent_inundated"];
35
36#[derive(Deserialize)]
37#[serde(deny_unknown_fields)]
38pub struct Request {
39    pub schema: String,
40    /// Horizontal request whose profiles are the stored output times, in order.
41    pub horizontal: horizontal::Request,
42    /// Output time of each profile in hours, strictly increasing.
43    pub times_hours: Vec<f64>,
44    /// Wet when the rendered depth is strictly greater than this (model units).
45    pub threshold: f32,
46    /// Little-endian float32 cell minimum elevations, in the horizontal
47    /// request's cell order (the native thresholded cell gate).
48    pub cell_minimum_file: String,
49    #[serde(default = "all_products")]
50    pub products: Vec<String>,
51    /// Deliver products even when unqualified cases are present (still reported).
52    #[serde(default)]
53    pub allow_unqualified: bool,
54}
55
56/// Default: the products qualified bit for bit. Duration is opt-in (unqualified).
57fn all_products() -> Vec<String> {
58    vec!["arrival_time".into(), "percent_inundated".into()]
59}
60
61pub fn method_record() -> Value {
62    json!({
63        "method_id": METHOD_ID,
64        "render": "horizontal-cell-constant/2 per stored output time; physical 2D cells only",
65        "wet": "rendered depth valid and depth > threshold (float32, strict)",
66        "percent_inundated": "float32(100 * wet_count / n_times)",
67        "arrival_time": "first wet stored output time (hours)",
68        "duration": "duration_end - arrival (hours)",
69        "nodata": "never wet, including unmapped pixels",
70        "memory": "32 bytes per mapped pixel plus 12 bytes per window pixel for products; output times are evaluated sequentially, pixels in parallel",
71        "cell_gate": "wse > cell_minimum + threshold (float32)",
72        "duration_end": "first dry output time after the last wet one, unless that is the final output time or none (then the last wet time)",
73        "time_arithmetic": "float32 hours",
74        "qualified": "native HEC-RAS 7.0.1 horizontal and sloping stored maps, Davis p02 (217 times) and BaldEagleCrkMulti2D bridges p05 (433 times), threshold 0.5 ft; arrival_time and percent_inundated",
75        "unqualified_refused": ["irregular_intervals", "duration"],
76        "reported": ["initially_wet_pixels", "multiple_period_pixels", "threshold_equal_samples"],
77        "known_differences": "duration: Davis 103/4,879 and bridges 11/959,043 pixels one output interval longer than native; arrival and percent identical",
78        "not_implemented": "sloping rendering; recession; subinterval interpolation"
79    })
80}
81
82#[derive(Clone, Copy)]
83struct Acc {
84    first: u32,
85    last: u32,
86    count: u32,
87    wet: bool,
88    periods: u8,
89}
90
91const NEVER: u32 = u32::MAX;
92
93/// Products and unqualified-case statistics from the accumulated pixels.
94struct Accumulated {
95    arrival: Vec<f32>,
96    duration: Vec<f32>,
97    percent: Vec<f32>,
98    wet_pixels: u64,
99    initially_wet: u64,
100    multiple_periods: u64,
101}
102
103/// Scatters compact accumulators (window pixel `index[k]`) to window products.
104fn finish(acc: &[Acc], index: &[usize], window: usize, times: &[f64]) -> Accumulated {
105    let n = times.len() as f64;
106    let mut out = Accumulated {
107        arrival: vec![f32::NAN; window],
108        duration: vec![f32::NAN; window],
109        percent: vec![f32::NAN; window],
110        wet_pixels: 0,
111        initially_wet: 0,
112        multiple_periods: 0,
113    };
114    for (a, &i) in acc.iter().zip(index) {
115        if a.count == 0 {
116            continue;
117        }
118        out.wet_pixels += 1;
119        out.initially_wet += u64::from(a.first == 0);
120        out.multiple_periods += u64::from(a.periods > 1);
121        let first = times[a.first as usize] as f32;
122        let last = a.last as usize;
123        let end = if last + 2 < times.len() {
124            last + 1
125        } else {
126            last
127        };
128        out.arrival[i] = first;
129        out.duration[i] = times[end] as f32 - first;
130        out.percent[i] = (100.0 * f64::from(a.count) / n) as f32;
131    }
132    out
133}
134
135#[inline]
136fn step(a: &mut Acc, position: u32, depth: Option<f32>, threshold: f32) -> u64 {
137    let wet = depth.is_some_and(|d| d > threshold);
138    if wet {
139        if a.first == NEVER {
140            a.first = position;
141        }
142        if !a.wet {
143            a.periods = a.periods.saturating_add(1);
144        }
145        a.last = position;
146        a.count += 1;
147    }
148    a.wet = wet;
149    u64::from(depth == Some(threshold))
150}
151
152/// Mapped pixels only (2D-owned with terrain); every other window pixel is
153/// never wet and stays NoData.
154#[derive(Default)]
155struct Accumulator {
156    threshold: f32,
157    /// Per 2D cell: minimum + threshold (float32), the native cell gate.
158    gate: Vec<f32>,
159    index: Vec<usize>,
160    source: Vec<u32>,
161    ground: Vec<f32>,
162    acc: Vec<Acc>,
163    equal: u64,
164}
165
166impl horizontal::Consumer for Accumulator {
167    fn prepare(&mut self, owner: &[u32], ground: &[f32]) -> Result<(), MapError> {
168        for (i, (&label, &g)) in owner.iter().zip(ground).enumerate() {
169            if label != 0 && g.is_finite() {
170                self.index.push(i);
171                self.source.push(label - 1);
172                self.ground.push(g);
173            }
174        }
175        self.acc = vec![
176            Acc {
177                first: NEVER,
178                last: NEVER,
179                count: 0,
180                wet: false,
181                periods: 0,
182            };
183            self.index.len()
184        ];
185        Ok(())
186    }
187
188    fn profile(&mut self, position: usize, values: &[f32]) -> Result<(), MapError> {
189        let (threshold, gate) = (self.threshold, &self.gate);
190        let position = position as u32;
191        self.equal += self
192            .acc
193            .par_iter_mut()
194            .zip(self.source.par_iter())
195            .zip(self.ground.par_iter())
196            .map(|((a, &source), &g)| -> Result<u64, MapError> {
197                let value = values[source as usize];
198                // Fails the thresholded cell gate: never wet at this time.
199                let passes = value > gate[source as usize];
200                if !passes {
201                    a.wet = false;
202                    return Ok(0);
203                }
204                Ok(step(a, position, horizontal::plot(value, g)?, threshold))
205            })
206            .try_reduce(|| 0, |x, y| Ok(x + y))?;
207        Ok(())
208    }
209}
210
211fn regular(times: &[f64]) -> bool {
212    times.len() < 3 || {
213        let step = times[1] - times[0];
214        times
215            .windows(2)
216            .all(|w| ((w[1] - w[0]) - step).abs() <= 1e-6)
217    }
218}
219
220pub fn run(req: &Request, base: &Path, out: &Path) -> Result<Value, MapError> {
221    let started = Instant::now();
222    let h = &req.horizontal;
223    let n = req.times_hours.len();
224    if req.schema != SCHEMA
225        || n == 0
226        || n != h.profiles.len()
227        || n >= u32::MAX as usize
228        || req.times_hours.iter().any(|t| !t.is_finite())
229        || req.times_hours.windows(2).any(|w| w[1] <= w[0])
230        || !req.threshold.is_finite()
231        || req.threshold < 0.0
232    {
233        return Err(MapError::Invalid(
234            "whole-simulation schema, times (finite, strictly increasing, one per profile) or threshold".into(),
235        ));
236    }
237    if h.stream_outputs || h.min_depth != 0.0 {
238        return Err(MapError::Invalid(
239            "whole-simulation horizontal requests must not stream and must not add a depth cutoff"
240                .into(),
241        ));
242    }
243    let mut products = req.products.clone();
244    products.sort();
245    products.dedup();
246    if products.is_empty()
247        || products.len() != req.products.len()
248        || products.iter().any(|p| !PRODUCTS.contains(&p.as_str()))
249    {
250        return Err(MapError::Invalid(
251            "products must be unique arrival_time, duration or percent_inundated".into(),
252        ));
253    }
254    if products.iter().any(|p| p == "duration") && !req.allow_unqualified {
255        return Err(MapError::Unsupported(
256            "duration is not qualified (native durations are one output interval shorter on some recession pixels); set allow_unqualified to deliver it".into(),
257        ));
258    }
259    let [_, _, width, height] = h.window;
260    let count = width
261        .checked_mul(height)
262        .filter(|&c| c > 0 && c <= crate::model::MAX_PROFILE_PIXELS as usize)
263        .ok_or_else(|| MapError::Budget("whole-simulation window empty or too large".into()))?;
264    // Compact accumulators (at most 32 B per window pixel) and the three
265    // float32 products (12 B) are admitted out of the caller's budget before
266    // the horizontal plan is checked.
267    let own = (count as u64)
268        .checked_mul(44)
269        .ok_or_else(|| MapError::Budget("whole-simulation plan overflow".into()))?;
270    let budget = h
271        .max_working_bytes
272        .checked_sub(own)
273        .ok_or_else(|| MapError::Budget("whole-simulation accumulators exceed budget".into()))?;
274    let horizontal_req = horizontal::Request {
275        max_working_bytes: budget,
276        ..h.clone()
277    };
278    let minimum = horizontal::read_cell_values(
279        &crate::model::resolve(base, &req.cell_minimum_file),
280        h.cell_count,
281    )?;
282    if minimum.iter().any(|m| !m.is_finite()) {
283        return Err(MapError::Invalid(
284            "cell minimum elevations must be finite".into(),
285        ));
286    }
287    let mut accumulator = Accumulator {
288        threshold: req.threshold,
289        gate: minimum.iter().map(|m| m + req.threshold).collect(),
290        ..Default::default()
291    };
292    let mut report = horizontal::run_with(&horizontal_req, base, out, Some(&mut accumulator))?;
293    report["schema"] = "ras-commander-whole-simulation-report/1".into();
294    report["whole_simulation"] = json!({"method": method_record(), "times": n,
295        "threshold": req.threshold, "products": products, "status": "not_run"});
296    if report["status"] != "complete" {
297        report["whole_simulation"]["status"] = "aborted".into();
298        sink::write_json(&out.join("report.json"), &report)?;
299        return Ok(report);
300    }
301    let tick = Instant::now();
302    let result = finish(
303        &accumulator.acc,
304        &accumulator.index,
305        count,
306        &req.times_hours,
307    );
308    let equal_samples = accumulator.equal;
309    drop(accumulator);
310    let irregular = !regular(&req.times_hours);
311    let unqualified = json!({
312        "irregular_intervals": irregular,
313        "initially_wet_pixels": result.initially_wet,
314        "multiple_period_pixels": result.multiple_periods,
315        "threshold_equal_samples": equal_samples,
316    });
317    let any = irregular;
318    report["whole_simulation"]["wet_pixels"] = result.wet_pixels.into();
319    report["whole_simulation"]["unqualified"] = unqualified;
320    if any && !req.allow_unqualified {
321        report["status"] = "aborted".into();
322        report["error"] =
323            "unsupported: unqualified whole-simulation cases present (see whole_simulation.unqualified); set allow_unqualified to deliver them".into();
324        report["whole_simulation"]["status"] = "refused".into();
325        sink::write_json(&out.join("report.json"), &report)?;
326        return Err(MapError::Unsupported(
327            "unqualified whole-simulation cases present".into(),
328        ));
329    }
330    let dir = out.join("whole_simulation");
331    fs::create_dir(&dir)?;
332    for product in &products {
333        let values = match product.as_str() {
334            "arrival_time" => &result.arrival,
335            "duration" => &result.duration,
336            _ => &result.percent,
337        };
338        horizontal::write_values(&dir.join(format!("{product}.bin")), values)?;
339    }
340    sink::write_json(
341        &dir.join("unit.json"),
342        &json!({"id":"whole_simulation","status":"complete","products":products}),
343    )?;
344    report["whole_simulation"]["status"] = "complete".into();
345    report["whole_simulation"]["qualified"] = (!any).into();
346    report["whole_simulation"]["product_seconds"] = tick.elapsed().as_secs_f64().into();
347    report["total_seconds"] = started.elapsed().as_secs_f64().into();
348    sink::write_json(&out.join("report.json"), &report)?;
349    Ok(report)
350}
351
352#[cfg(test)]
353mod tests {
354    use super::*;
355    use crate::horizontal::Consumer;
356
357    fn run_series(series: &[Vec<f32>], threshold: f32, times: &[f64]) -> (Accumulated, u64) {
358        // One cell per pixel on ground 0 (exact float32 depths); NaN = dry cell.
359        let n = series[0].len();
360        let owner: Vec<u32> = (1..=n as u32).collect();
361        let mut a = Accumulator {
362            threshold,
363            gate: vec![f32::NEG_INFINITY; n],
364            ..Default::default()
365        };
366        a.prepare(&owner, &vec![0.0; n]).unwrap();
367        for (k, values) in series.iter().enumerate() {
368            a.profile(k, values).unwrap();
369        }
370        (finish(&a.acc, &a.index, n, times), a.equal)
371    }
372
373    #[test]
374    fn strict_threshold_first_last_count_and_periods() {
375        let times = [0.0, 0.25, 0.5, 0.75, 1.0];
376        let series = vec![
377            vec![0.1, 0.6, 0.6, 0.6, 0.4, 0.6, 0.5, f32::NAN],
378            vec![0.2, 0.7, 0.4, 0.1, 0.4, 0.6, 0.5, f32::NAN],
379            vec![0.3, 0.8, 0.4, 0.7, 0.6, 0.4, 0.5, f32::NAN],
380            vec![0.4, 0.4, 0.9, 0.1, 0.7, 0.4, 0.5, f32::NAN],
381            vec![0.1, 0.1, 0.1, 0.1, 0.4, 0.4, 0.5, f32::NAN],
382        ];
383        let (r, equal) = run_series(&series, 0.5, &times);
384        assert_eq!(equal, 5); // pixel 6 equals the threshold: reported, dry
385        assert!(r.arrival[0].is_nan() && r.percent[0].is_nan() && r.duration[0].is_nan());
386        // Pixel 1 wet 0..2, dry at 3 (not the final time): ends at t[3].
387        assert_eq!(
388            (r.arrival[1], r.duration[1], r.percent[1]),
389            (0.0, 0.75, 60.0)
390        );
391        // Pixel 5 wet 0..1, dry at 2: ends at t[2].
392        assert_eq!((r.arrival[5], r.duration[5]), (0.0, 0.5));
393        // Pixels 2 and 4: last wet at 3, dry only at the final time: end at t[3].
394        assert_eq!(
395            (r.arrival[2], r.duration[2], r.percent[2]),
396            (0.0, 0.75, 40.0)
397        );
398        assert_eq!(
399            (r.arrival[4], r.duration[4], r.percent[4]),
400            (0.5, 0.25, 40.0)
401        );
402        assert!(r.arrival[6].is_nan() && r.arrival[7].is_nan());
403        assert_eq!(r.wet_pixels, 5);
404        assert_eq!(r.initially_wet, 4); // pixels 1, 2, 3 and 5
405        assert_eq!(r.multiple_periods, 2); // pixels 2 and 3
406    }
407
408    #[test]
409    fn single_wet_sample_and_plot_tolerance() {
410        let series: Vec<Vec<f32>> = [0.0_f32, 0.0, 0.0, 1.0, 0.0]
411            .iter()
412            .map(|d| vec![*d])
413            .collect();
414        let (r, _) = run_series(&series, 0.5, &[0.0, 5.0, 10.0, 15.0, 20.0]);
415        assert_eq!(
416            (r.arrival[0], r.duration[0], r.percent[0]),
417            (15.0, 0.0, 20.0)
418        );
419        assert_eq!((r.initially_wet, r.multiple_periods), (0, 0));
420        // Threshold 0: a depth under the native plot tolerance is not visible, so not wet.
421        let (r, _) = run_series(&[vec![0.0005], vec![0.002]], 0.0, &[0.0, 1.0]);
422        assert_eq!((r.arrival[0], r.percent[0]), (1.0, 50.0));
423    }
424
425    #[test]
426    fn plot_matches_the_horizontal_kernel() {
427        let ground = [1000.0_f32, 1000.0, 100000.0, 10.0, f32::NAN, 5.0];
428        let values = [
429            f32::from_bits(0x447a000f),
430            f32::from_bits(0x447a0010),
431            100000.0,
432            10.0005,
433            12.0,
434            f32::NAN,
435        ];
436        for (&v, &g) in values.iter().zip(&ground) {
437            let mut depth = [f32::NAN];
438            let mut wse = [f32::NAN];
439            horizontal::evaluate_for_test(&[1], &[g], &[v], &mut depth, &mut wse);
440            let plotted = horizontal::plot(v, g).unwrap();
441            assert_eq!(
442                plotted.map(f32::to_bits),
443                (!depth[0].is_nan()).then(|| depth[0].to_bits())
444            );
445        }
446    }
447
448    #[test]
449    fn thresholded_cell_gate_and_wet_to_the_end() {
450        // Two cells, pixel depth always 0.6 above ground 0; cell 0 minimum 0.2
451        // (gate 0.7 fails), cell 1 minimum -1 (gate -0.5 passes).
452        let mut a = Accumulator {
453            threshold: 0.5,
454            gate: vec![0.2 + 0.5, -1.0 + 0.5],
455            ..Default::default()
456        };
457        a.prepare(&[1, 2], &[0.0, 0.0]).unwrap();
458        for k in 0..3 {
459            a.profile(k, &[0.6, 0.6]).unwrap();
460        }
461        let r = finish(&a.acc, &a.index, 2, &[0.0, 1.0, 2.0]);
462        assert!(r.arrival[0].is_nan());
463        // Wet at the final time: duration ends at the last wet time.
464        assert_eq!(
465            (r.arrival[1], r.duration[1], r.percent[1]),
466            (0.0, 2.0, 100.0)
467        );
468    }
469
470    #[test]
471    fn times_are_float32_hours() {
472        let times: Vec<f64> = (0..865).map(|k| k as f64 / 12.0).collect();
473        let mut acc = vec![Acc {
474            first: 244,
475            last: 488,
476            count: 245,
477            wet: false,
478            periods: 1,
479        }];
480        let r = finish(&acc, &[0], 1, &times);
481        assert_eq!(
482            r.duration[0].to_bits(),
483            ((489.0_f64 / 12.0) as f32 - (244.0_f64 / 12.0) as f32).to_bits()
484        );
485        acc[0].last = 864;
486        let r = finish(&acc, &[0], 1, &times);
487        assert_eq!(r.duration[0], (72.0_f32) - (244.0_f64 / 12.0) as f32);
488    }
489
490    #[test]
491    fn interval_regularity() {
492        assert!(regular(&[0.0, 1.0 / 12.0, 2.0 / 12.0, 3.0 / 12.0]));
493        assert!(regular(&[0.0, 7.0]));
494        assert!(!regular(&[0.0, 1.0, 2.5]));
495    }
496}
ras_commander_fast_map/whole_simulation.rs

ras_commander_fast_map/
whole_simulation.rs
