//! Trail decimation and encoding.
//!
//! A day of walking at one fix per 25 s is ~3500 points. Sent as JSON floats that
//! is ~100 kB; as a Ramer–Douglas–Peucker-simplified encoded polyline it is a few
//! kilobytes, and the browser has proportionally less to draw. Both halves matter:
//! simplification removes points that would render on top of each other anyway,
//! and the encoding removes the per-number punctuation.
//!
//! Coordinates arrive as degrees × 1e7 (the protocol and database representation)
//! and are encoded at 1e5, which is the polyline format's fixed precision — about
//! 1.1 m, well under the accuracy of any consumer GPS fix.

/// Points are `(lat_e7, lon_e7)`.
type Pt = (i64, i64);

/// Simplify to at most `max` points using Ramer–Douglas–Peucker.
///
/// The tolerance is searched rather than fixed, because the useful question is
/// "how much detail fits in the budget?" and the answer depends on how far the
/// person actually travelled. A fixed ε either keeps too much on a road trip or
/// destroys a walk around a park.
pub fn simplify(points: &[Pt], max: usize) -> Vec<Pt> {
    if points.len() <= max {
        return points.to_vec();
    }

    // Ramer–Douglas–Peucker is O(n log n) on a typical track but O(n²) in the
    // worst case — a stationary phone reporting a heartbeat all week is exactly
    // that case — and the tolerance search below runs it several times. Uniformly
    // subsampling first bounds the work at a size where that no longer matters.
    // At 8× the budget the subsampling itself removes no visible detail: anything
    // it drops would have been RDP's next victim anyway.
    let owned;
    let points = if points.len() > max.saturating_mul(PRE_DECIMATE_FACTOR) {
        owned = uniform(points, max * PRE_DECIMATE_FACTOR);
        owned.as_slice()
    } else {
        points
    };

    // ~15 m expressed in units of 1e-7 degrees.
    let base: f64 = 15.0 * 1e7 / 111_320.0;

    // Grow until something fits, to bracket the answer. Plain doubling would stop
    // at the *first* tolerance that fits, and on a track with real large-scale
    // shape that first hit overshoots: a meander can collapse to its two
    // endpoints, which is a straight line where there was a walk. So this phase
    // only establishes an upper bound.
    let mut lo = 0.0f64; // known to yield more than `max`
    let mut hi = base;
    let mut best: Option<Vec<Pt>> = None;
    for _ in 0..24 {
        let out = rdp(points, hi);
        if out.len() <= max {
            best = Some(out);
            break;
        }
        lo = hi;
        hi *= 4.0;
    }

    // Then bisect, keeping the largest result that still fits: the most detail the
    // budget allows, rather than the first thing under it.
    if let Some(mut candidate) = best {
        for _ in 0..8 {
            let mid = (lo + hi) / 2.0;
            let out = rdp(points, mid);
            if out.len() <= max {
                if out.len() >= candidate.len() {
                    candidate = out;
                }
                hi = mid;
            } else {
                lo = mid;
            }
        }
        return candidate;
    }

    // No tolerance in 24 quadruplings fit the budget. Unreachable for real data;
    // uniform sampling keeps the shape far better than an infinite tolerance.
    uniform(points, max)
}

/// How much larger than the budget the RDP input may be.
const PRE_DECIMATE_FACTOR: usize = 8;

/// Evenly spaced subsample of at most `max` points, always keeping the last one so
/// a trail still ends where the person is.
fn uniform(points: &[Pt], max: usize) -> Vec<Pt> {
    if points.len() <= max {
        return points.to_vec();
    }
    let step = points.len().div_ceil(max);
    let mut out: Vec<Pt> = points.iter().copied().step_by(step).collect();
    if out.last() != points.last() {
        out.push(*points.last().expect("non-empty"));
    }
    out
}

fn rdp(points: &[Pt], epsilon: f64) -> Vec<Pt> {
    if points.len() < 3 {
        return points.to_vec();
    }
    let mut keep = vec![false; points.len()];
    keep[0] = true;
    keep[points.len() - 1] = true;
    // Iterative rather than recursive: a 100k-point input would otherwise be able
    // to blow the stack, and this runs on attacker-influenced data volumes.
    let mut stack = vec![(0usize, points.len() - 1)];
    while let Some((start, end)) = stack.pop() {
        if end <= start + 1 {
            continue;
        }
        let mut worst = 0.0;
        let mut worst_i = start;
        for i in (start + 1)..end {
            let d = perpendicular_distance(points[i], points[start], points[end]);
            if d > worst {
                worst = d;
                worst_i = i;
            }
        }
        if worst > epsilon {
            keep[worst_i] = true;
            stack.push((start, worst_i));
            stack.push((worst_i, end));
        }
    }
    points
        .iter()
        .zip(keep)
        .filter_map(|(p, k)| k.then_some(*p))
        .collect()
}

fn perpendicular_distance(p: Pt, a: Pt, b: Pt) -> f64 {
    let (px, py) = (p.1 as f64, p.0 as f64);
    let (ax, ay) = (a.1 as f64, a.0 as f64);
    let (bx, by) = (b.1 as f64, b.0 as f64);
    let dx = bx - ax;
    let dy = by - ay;
    let len_sq = dx * dx + dy * dy;
    if len_sq == 0.0 {
        return ((px - ax).powi(2) + (py - ay).powi(2)).sqrt();
    }
    ((dx * (ay - py) - (ax - px) * dy).abs()) / len_sq.sqrt()
}

/// Google's encoded polyline algorithm at 1e5 precision.
pub fn encode(points: &[Pt]) -> String {
    let mut out = String::with_capacity(points.len() * 6);
    let mut prev_lat = 0i64;
    let mut prev_lon = 0i64;
    for &(lat_e7, lon_e7) in points {
        // 1e7 -> 1e5, rounding rather than truncating so error stays centred.
        let lat = div_round(lat_e7, 100);
        let lon = div_round(lon_e7, 100);
        encode_value(lat - prev_lat, &mut out);
        encode_value(lon - prev_lon, &mut out);
        prev_lat = lat;
        prev_lon = lon;
    }
    out
}

fn div_round(v: i64, d: i64) -> i64 {
    if v >= 0 {
        (v + d / 2) / d
    } else {
        -((-v + d / 2) / d)
    }
}

fn encode_value(value: i64, out: &mut String) {
    let mut v = if value < 0 { !(value << 1) } else { value << 1 };
    while v >= 0x20 {
        out.push(char::from(((0x20 | (v & 0x1f)) + 63) as u8));
        v >>= 5;
    }
    out.push(char::from((v + 63) as u8));
}

/// Decode. Used by the tests to prove the encoder is reversible, and kept public
/// because a Rust consumer of `/api/users/:id/track` needs it.
#[cfg_attr(not(test), allow(dead_code))]
pub fn decode(s: &str) -> Vec<Pt> {
    let bytes = s.as_bytes();
    let mut out = Vec::new();
    let mut i = 0;
    let mut lat = 0i64;
    let mut lon = 0i64;
    while i < bytes.len() {
        let Some(dlat) = decode_value(bytes, &mut i) else {
            break;
        };
        let Some(dlon) = decode_value(bytes, &mut i) else {
            break;
        };
        lat += dlat;
        lon += dlon;
        out.push((lat * 100, lon * 100));
    }
    out
}

#[cfg_attr(not(test), allow(dead_code))]
fn decode_value(bytes: &[u8], i: &mut usize) -> Option<i64> {
    let mut shift = 0;
    let mut result = 0i64;
    loop {
        let b = *bytes.get(*i)? as i64 - 63;
        *i += 1;
        result |= (b & 0x1f) << shift;
        shift += 5;
        if b < 0x20 {
            break;
        }
        if shift > 60 {
            return None; // malformed; refuse to shift forever
        }
    }
    Some(if result & 1 != 0 {
        !(result >> 1)
    } else {
        result >> 1
    })
}

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

    #[test]
    fn a_known_vector_matches_the_reference_implementation() {
        // The example from Google's polyline documentation, in 1e7 units.
        // (38.5, -120.2), (40.7, -120.95), (43.252, -126.453)
        let points = vec![
            (385_000_000, -1_202_000_000),
            (407_000_000, -1_209_500_000),
            (432_520_000, -1_264_530_000),
        ];
        assert_eq!(encode(&points), "_p~iF~ps|U_ulLnnqC_mqNvxq`@");
    }

    #[test]
    fn encode_decode_round_trips_within_the_formats_precision() {
        let points = vec![
            (525_200_080, 134_050_000),
            (525_210_000, 134_060_000),
            (-338_688_000, -1_754_500_000),
        ];
        let back = decode(&encode(&points));
        assert_eq!(back.len(), points.len());
        for (a, b) in points.iter().zip(&back) {
            // 1e5 precision means the last two 1e7 digits are lost: ≤ 50 units,
            // about 0.5 cm. Well inside any GPS accuracy.
            assert!((a.0 - b.0).abs() <= 50, "lat drifted: {a:?} vs {b:?}");
            assert!((a.1 - b.1).abs() <= 50, "lon drifted: {a:?} vs {b:?}");
        }
    }

    #[test]
    fn an_empty_track_encodes_to_an_empty_string() {
        assert_eq!(encode(&[]), "");
        assert_eq!(decode(""), Vec::<Pt>::new());
    }

    #[test]
    fn decoding_garbage_does_not_panic_or_hang() {
        for s in ["~", "?????", "\u{1}\u{2}\u{3}", &"~".repeat(1000)] {
            let _ = decode(s);
        }
    }

    #[test]
    fn simplify_keeps_short_tracks_untouched() {
        let points: Vec<Pt> = (0..10).map(|i| (i * 1000, i * 1000)).collect();
        assert_eq!(simplify(&points, 2000), points);
    }

    #[test]
    fn simplify_respects_the_budget() {
        // A 20k-point meander: a long sinusoidal path with per-fix jitter on top.
        // Large-scale shape matters here — a straight line with noise legitimately
        // simplifies to two points, so it would not test anything.
        let points: Vec<Pt> = (0..20_000i64)
            .map(|i| {
                let jitter = if i % 2 == 0 { 300 } else { -300 };
                let wave = (5_000_000.0 * ((i as f64) / 300.0).sin()) as i64;
                (525_200_000 + i * 40 + jitter, 134_050_000 + wave - jitter)
            })
            .collect();
        let out = simplify(&points, 500);
        assert!(out.len() <= 500, "budget exceeded: {} points", out.len());
        // And it must actually *use* the budget. Returning two endpoints would
        // satisfy the limit while turning a wander into a straight line.
        assert!(
            out.len() > 250,
            "budget under-used: only {} points",
            out.len()
        );
    }

    #[test]
    fn simplify_keeps_the_endpoints() {
        let points: Vec<Pt> = (0..5_000)
            .map(|i| (525_200_000 + i * 100, 134_050_000))
            .collect();
        let out = simplify(&points, 100);
        assert_eq!(
            out.first(),
            points.first(),
            "the start of a trail must survive"
        );
        assert_eq!(out.last(), points.last(), "the end of a trail must survive");
    }

    #[test]
    fn simplify_drops_collinear_interior_points() {
        // A dead-straight line: everything between the ends is redundant.
        let points: Vec<Pt> = (0..3_000)
            .map(|i| (525_200_000 + i * 1_000, 134_050_000))
            .collect();
        let out = simplify(&points, 2_000);
        assert!(
            out.len() < 50,
            "a straight line should collapse, got {}",
            out.len()
        );
    }

    #[test]
    fn simplify_survives_a_track_that_never_moves() {
        // Every point identical: len_sq == 0 in the distance function.
        let points: Vec<Pt> = std::iter::repeat_n((525_200_000, 134_050_000), 5_000).collect();
        let out = simplify(&points, 100);
        assert!(out.len() <= 100);
    }

    #[test]
    fn negative_coordinates_encode_correctly() {
        // The zig-zag encoding of negatives is the classic place to get an
        // off-by-one, and half the planet is at a negative longitude.
        let points = vec![(-338_688_000, -1_754_500_000)];
        let back = decode(&encode(&points));
        assert_eq!(back.len(), 1);
        assert!(back[0].0 < 0 && back[0].1 < 0);
    }
}