vsrg 0.3.0

use crate::math::easing::eval::easing_value;
use crate::math::easing::integration::integrate_simple;
use crate::math::{Easable, Easing};

pub fn solve_integration_evaluation<T, const N: usize>(
    ease: Easing,
    from: T,
    to: T,
    val: T,
) -> Option<T>
where
    T: Easable<N>,
{
    let from = from.components();
    let to = to.components();
    let val = val.components();

    for i in 0..N {
        // Monotonic if there's no sign flip.
        if from[i] * to[i] < 0.0 {
            return None;
        }
        if !from[i].is_finite() || !to[i].is_finite() || !val[i].is_finite() {
            return None;
        }
    }

    Some(T::compose(std::array::from_fn(|i| {
        let from = from[i];
        let to = to[i];
        let val = val[i];
        if from < 0.0 || to < 0.0 {
            solve_integration(ease, -from, -to, -val)
        } else {
            solve_integration(ease, from, to, val)
        }
    })))
}

fn solve_integration(ease: Easing, from: f64, to: f64, val: f64) -> f64 {
    let delta = to - from;

    // Clamp to the achievable range [F(0), F(1)] = [0, from*(1-I(1)) + to*I(1)].
    let f0 = 0.0;
    let i1 = integrate_simple(ease, 1.0);
    let f1 = from * (1.0 - i1) + to * i1;
    let (lo, hi) = if f0 <= f1 { (f0, f1) } else { (f1, f0) };
    if val <= lo {
        return 0.0;
    }
    if val >= hi {
        return 1.0;
    }

    match ease {
        Easing::Linear => {
            // F(t) = from*(t - t^2/2) + to*t^2/2 = from*t + delta*t^2/2.
            // Solve delta*t^2/2 + from*t - val = 0.
            // t = (-from + sqrt(from^2 + 2*delta*val)) / delta  (delta != 0).
            if delta.abs() <= 1e-12 {
                if from.abs() <= 1e-12 {
                    0.0
                } else {
                    (val / from).clamp(0.0, 1.0)
                }
            } else {
                let disc = from * from + 2.0 * delta * val;
                let t = (-from + f64::sqrt(disc.max(0.0))) / delta;
                t.clamp(0.0, 1.0)
            }
        }

        Easing::InConst => {
            if from.abs() <= 1e-12 {
                0.0
            } else {
                (val / from).clamp(0.0, 1.0)
            }
        }
        Easing::OutConst => {
            if to.abs() <= 1e-12 {
                0.0
            } else {
                (val / to).clamp(0.0, 1.0)
            }
        }
        Easing::InOutConst => {
            let f_at_half = from * 0.5;
            if val <= f_at_half {
                if from.abs() <= 1e-12 {
                    0.0
                } else {
                    (val / from).clamp(0.0, 0.5)
                }
            } else if to.abs() <= 1e-12 {
                0.5
            } else {
                0.5 + ((val - f_at_half) / to).clamp(0.0, 0.5)
            }
        }

        // All remaining easings: bisection on the monotonic forward F(t).
        _ => newton_bisect_solve(ease, from, to, val),
    }
}

/// Newton + bisection hybrid on the monotonic forward F(t) = from*(t-I(t)) + to*I(t).
///
/// F'(t) = from*(1-e(t)) + to*e(t) (the eased value, by FTC) is cheap
/// and already analytical, so Newton steps converge quadratically away from
/// flat regions. The bracket [lo, hi] is maintained with the bisection invariant
/// F(lo) < target <= F(hi); any Newton step that leaves the bracket (or when
/// F'(t) ~ 0 on a flat region) falls back to a bisection step. The `>=`
/// comparison ensures the lower bound is returned on flat regions.
fn newton_bisect_solve(ease: Easing, from: f64, to: f64, target: f64) -> f64 {
    let f = |t: f64| {
        let i = integrate_simple(ease, t);
        from * (t - i) + to * i
    };
    let fp = |t: f64| {
        let e = easing_value(ease.as_easing(), t);
        from * (1.0 - e) + to * e
    };

    let mut lo = 0.0_f64;
    let mut hi = 1.0_f64;

    for _ in 0..24 {
        if hi - lo <= f64::EPSILON * 4.0 {
            break;
        }
        let mid = 0.5 * (lo + hi);
        let fm = f(mid);
        let d = fp(mid);

        // Newton step from mid; fall back to mid if degenerate or out of bracket.
        let probe = if d.abs() > 1e-10 {
            let t_n = mid - (fm - target) / d;
            if t_n > lo && t_n < hi { t_n } else { mid }
        } else {
            mid
        };

        let pv = if probe == mid { fm } else { f(probe) };

        if pv >= target {
            hi = probe;
        } else {
            lo = probe;
        }
    }
    // hi is the lower bound of the t-range producing `target`.
    hi
}

#[cfg(test)]
mod tests {
    use crate::math::easing::integration;

    use super::*;

    const FROM: f64 = 0.0;
    const TO: f64 = 100.0;
    const SAMPLE_TS: &[f64] = &[0.1, 0.25, 0.5, 0.75, 0.9];

    const IN_SCOPE_EASINGS: &[Easing] = &[
        Easing::Linear,
        Easing::InSine,
        Easing::OutSine,
        Easing::InOutSine,
        Easing::InQuad,
        Easing::OutQuad,
        Easing::InOutQuad,
        Easing::InCubic,
        Easing::OutCubic,
        Easing::InOutCubic,
        Easing::InQuart,
        Easing::OutQuart,
        Easing::InOutQuart,
        Easing::InQuint,
        Easing::OutQuint,
        Easing::InOutQuint,
        Easing::InExpo,
        Easing::OutExpo,
        Easing::InOutExpo,
    ];

    fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() <= tol
    }

    #[test]
    fn round_trip_cumulate_then_solve() {
        // solve(ease, cumulate(ease, t)) == t
        let mut failures: Vec<(Easing, f64, f64)> = Vec::new();
        for &ease in IN_SCOPE_EASINGS {
            for &t in SAMPLE_TS {
                let val = integration::evaluate_simple(ease, FROM, TO, t);
                let solved = solve_integration_evaluation(ease, FROM, TO, val).unwrap();
                if !approx_eq(solved, t, 1e-3) {
                    failures.push((ease, t, solved));
                }
            }
        }
        assert!(
            failures.is_empty(),
            "solve(cumulate(t)) != t: {}",
            failures
                .iter()
                .map(|(e, t, s)| format!("{e:?}: t={t}, solved={s}"))
                .collect::<Vec<_>>()
                .join("; ")
        );
    }

    #[test]
    fn round_trip_solve_then_cumulate() {
        // cumulate(ease, solve(ease, V)) == V
        let mut failures: Vec<(Easing, f64, f64)> = Vec::new();
        for &ease in IN_SCOPE_EASINGS {
            for &t in SAMPLE_TS {
                let v = integration::evaluate_simple(ease, FROM, TO, t);
                let back = integration::evaluate_simple(
                    ease,
                    FROM,
                    TO,
                    solve_integration_evaluation(ease, FROM, TO, v).unwrap(),
                );
                if !approx_eq(back, v, 1e-2) {
                    failures.push((ease, v, back));
                }
            }
        }
        assert!(
            failures.is_empty(),
            "cumulate(solve(V)) != V: {}",
            failures
                .iter()
                .map(|(e, v, b)| format!("{e:?}: V={v}, back={b}"))
                .collect::<Vec<_>>()
                .join("; ")
        );
    }

    #[test]
    fn solve_integration_endpoints() {
        for &ease in IN_SCOPE_EASINGS {
            assert!(
                approx_eq(
                    solve_integration_evaluation(ease, FROM, TO, 0.0).unwrap(),
                    0.0,
                    1e-6
                ),
                "{ease:?}: solve(V=0) should be 0"
            );
            let v1 = integration::evaluate_simple(ease, FROM, TO, 1.0);
            assert!(
                approx_eq(
                    solve_integration_evaluation(ease, FROM, TO, v1).unwrap(),
                    1.0,
                    1e-3
                ),
                "{ease:?}: solve(V=F(1)) should be 1"
            );
        }
    }

    #[test]
    fn solve_integration_clamps_out_of_range() {
        for &ease in IN_SCOPE_EASINGS {
            assert!(
                approx_eq(
                    solve_integration_evaluation(ease, FROM, TO, -10.0).unwrap(),
                    0.0,
                    1e-6
                ),
                "{ease:?}: solve(V<0) should clamp to 0"
            );
            let v1 = integration::evaluate_simple(ease, FROM, TO, 1.0);
            assert!(
                approx_eq(
                    solve_integration_evaluation(ease, FROM, TO, v1 + 10.0).unwrap(),
                    1.0,
                    1e-6
                ),
                "{ease:?}: solve(V>F(1)) should clamp to 1"
            );
        }
    }

    #[test]
    fn sign_flip_endpoints_return_none() {
        // from and to with opposite signs -> non-monotonic -> None
        for &ease in IN_SCOPE_EASINGS {
            assert!(
                solve_integration_evaluation(ease, -1.0, 10.0, 5.0).is_none(),
                "{ease:?}: from<0, to>0 should return None"
            );
            assert!(
                solve_integration_evaluation(ease, 10.0, -1.0, 5.0).is_none(),
                "{ease:?}: from>0, to<0 should return None"
            );
        }
    }

    #[test]
    fn negative_same_sign_endpoints_solve() {
        // Both endpoints negative (same sign) is still monotonic; solve by
        // negating. With from=-100, to=0, Linear: F(t) = -100*t + 100*t^2/2.
        // F(0.5) = -50 + 12.5 = -37.5; solving for -37.5 should give t=0.5.
        let r = solve_integration_evaluation(Easing::Linear, -100.0, 0.0, -37.5).unwrap();
        assert!(
            approx_eq(r, 0.5, 1e-4),
            "Linear with negative endpoints: expected t=0.5, got {r}"
        );
    }

    #[test]
    fn from_equals_to_returns_lower_bound() {
        // F(t) = c*t. With from==to==c, solve any V in range -> V/c.
        let c: f64 = 50.0;
        for &ease in IN_SCOPE_EASINGS {
            let r = solve_integration_evaluation(ease, c, c, 25.0).unwrap();
            assert!(
                approx_eq(r, 0.5, 1e-6),
                "{ease:?}: from==to=c, V=c/2 should give t=0.5, got {r}"
            );
        }
    }

    #[test]
    fn linear_matches_closed_form() {
        // F(t) = from*t + delta*t^2/2. With from=0, to=100: F(t)=50*t^2.
        // Range [0, 50].  t = sqrt(V/50).
        for &v in &[2.0_f64, 12.5, 25.0, 40.0, 50.0] {
            let t = solve_integration_evaluation(Easing::Linear, 0.0, 100.0, v).unwrap();
            let expected = f64::sqrt(v / 50.0);
            assert!(
                approx_eq(t, expected, 1e-5),
                "Linear: V={v}, expected t={expected}, got {t}"
            );
        }
    }
}