siderust 0.9.0

// SPDX-License-Identifier: AGPL-3.0-or-later
// Copyright (C) 2026 Vallés Puig, Ramon

//! # Bracketing and Seeding Policies
//!
//! ## Scientific scope
//!
//! Provides strategies for identifying candidate sub-intervals where a scalar
//! function may cross a threshold or contain an extremum.  Bracketing is the
//! prerequisite step before root-finding or golden-section search can be
//! applied.  Three policies cover progressively more complex detection
//! scenarios, from uniform probing to extremum-guided refinement.
//!
//! ## Technical scope
//!
//! All routines operate on `Interval<ModifiedJulianDate>` time windows and
//! closures `Fn(ModifiedJulianDate) → Quantity<V>`.
//!
//! | Policy | Description |
//! |--------|-------------|
//! | [`fixed_step_brackets`] | Uniform step across the window |
//! | [`adaptive_step_brackets`] | Narrows step near suspected events |
//! | [`extrema_based_brackets`] | Find extrema first, bracket crossings around each |
//!
//! ## References
//!
//! - Brent, R. P. (1973). *Algorithms for Minimization without Derivatives*.
//!   Prentice-Hall.
//! - Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P.
//!   (2007). *Numerical Recipes in C++*, 3rd ed. Cambridge University Press.

use crate::qtty::*;
use crate::time::{Interval, ModifiedJulianDate};

use super::extrema::{self, ExtremumKind};

type Mjd = ModifiedJulianDate;

#[inline]
fn opposite_sign<V: Unit>(a: Quantity<V>, b: Quantity<V>) -> bool {
    a.signum() * b.signum() < 0.0
}

// ---------------------------------------------------------------------------
// Fixed‑step seeding
// ---------------------------------------------------------------------------

/// Generate uniform brackets at a fixed step size.
///
/// Returns sign‑change brackets for `f(t) − threshold`.
pub fn fixed_step_brackets<V, F>(
    period: Interval<ModifiedJulianDate>,
    step: Days,
    f: &F,
    threshold: Quantity<V>,
) -> Vec<Interval<ModifiedJulianDate>>
where
    V: Unit,
    F: Fn(ModifiedJulianDate) -> Quantity<V>,
{
    let g = |t: Mjd| -> Quantity<V> { f(t) - threshold };

    let mut brackets = Vec::new();
    let mut t = period.start;
    let mut prev = g(t);

    while t < period.end {
        let next_t = {
            let t_next = crate::time::ModifiedJulianDate::new((t.raw() + step).value());
            if t_next.raw() <= period.end.raw() {
                t_next
            } else {
                period.end
            }
        };
        let next_v = g(next_t);

        if opposite_sign(prev, next_v) {
            brackets.push(Interval::new(t, next_t));
        }

        t = next_t;
        prev = next_v;
    }

    brackets
}

// ---------------------------------------------------------------------------
// Adaptive‑step seeding
// ---------------------------------------------------------------------------

/// Generate brackets with adaptive step: starts at `initial_step`, subdivides
/// where the function changes rapidly (derivative exceeds `rate_threshold`).
///
/// `min_step` prevents infinite subdivision.
pub fn adaptive_step_brackets<V, F>(
    period: Interval<ModifiedJulianDate>,
    initial_step: Days,
    min_step: Days,
    f: &F,
    threshold: Quantity<V>,
) -> Vec<Interval<ModifiedJulianDate>>
where
    V: Unit,
    F: Fn(ModifiedJulianDate) -> Quantity<V>,
{
    let g = |t: Mjd| -> Quantity<V> { f(t) - threshold };

    let mut brackets = Vec::new();

    struct Frame<V: Unit> {
        period: Interval<ModifiedJulianDate>,
        g_lo: Quantity<V>,
        g_hi: Quantity<V>,
    }

    let mut stack: Vec<Frame<V>> = Vec::new();

    // Initial pass at coarse step
    let mut t = period.start;
    let mut prev = g(t);
    while t < period.end {
        let next_t = {
            let t_next = crate::time::ModifiedJulianDate::new((t.raw() + initial_step).value());
            if t_next.raw() <= period.end.raw() {
                t_next
            } else {
                period.end
            }
        };
        let next_v = g(next_t);
        stack.push(Frame {
            period: Interval::new(t, next_t),
            g_lo: prev,
            g_hi: next_v,
        });
        t = next_t;
        prev = next_v;
    }

    // Process stack: subdivide large steps where sign change detected
    while let Some(frame) = stack.pop() {
        if opposite_sign(frame.g_lo, frame.g_hi) {
            let width = frame.period.end.raw() - frame.period.start.raw();
            if width <= min_step + min_step {
                brackets.push(frame.period);
            } else {
                // Subdivide to find tighter bracket
                let mid = crate::time::ModifiedJulianDate::new(
                    (frame.period.start.raw() + width * 0.5).value(),
                );
                let g_mid = g(mid);
                // Push both halves (will be processed)
                stack.push(Frame {
                    period: Interval::new(frame.period.start, mid),
                    g_lo: frame.g_lo,
                    g_hi: g_mid,
                });
                stack.push(Frame {
                    period: Interval::new(mid, frame.period.end),
                    g_lo: g_mid,
                    g_hi: frame.g_hi,
                });
            }
        }
    }

    brackets.sort_by(|a, b| a.start.partial_cmp(&b.start).unwrap());
    brackets
}

// ---------------------------------------------------------------------------
// Extrema‑based seeding (for satellite passes / fast movers)
// ---------------------------------------------------------------------------

/// Find extrema of `f` first, then generate crossing brackets around each
/// extremum where `f(extremum) > threshold` (or `< threshold` for minima that
/// dip below).
///
/// This is ideal for satellite pass detection: find the altitude peak of each
/// pass, then bracket the rise/set crossings on either side.
pub fn extrema_based_brackets<V, F>(
    period: Interval<ModifiedJulianDate>,
    extrema_step: Days,
    f: &F,
    threshold: Quantity<V>,
) -> Vec<Interval<ModifiedJulianDate>>
where
    V: Unit,
    F: Fn(ModifiedJulianDate) -> Quantity<V>,
{
    let extrema = extrema::find_extrema(period, extrema_step, f);

    let mut brackets = Vec::new();

    for ext in &extrema {
        if ext.kind == ExtremumKind::Maximum && ext.value > threshold {
            // This maximum is above threshold → there must be a rising crossing
            // before it and a setting crossing after it.
            // Search backward from the extremum for the rising crossing
            let rise_bracket =
                search_crossing_backward(Interval::new(period.start, ext.t), f, threshold);
            if let Some(br) = rise_bracket {
                brackets.push(br);
            }

            // Search forward from the extremum for the setting crossing
            let set_bracket =
                search_crossing_forward(Interval::new(ext.t, period.end), f, threshold);
            if let Some(br) = set_bracket {
                brackets.push(br);
            }
        }
    }

    // Deduplicate overlapping brackets
    brackets.sort_by(|a, b| a.start.partial_cmp(&b.start).unwrap());
    brackets.dedup_by(|a, b| {
        const TOL: Days = Days::new(1e-8); // 0.00000001 days ~ 0.86 seconds
        let dt_start = (a.start.raw() - b.start.raw()).abs();
        let dt_end = (a.end.raw() - b.end.raw()).abs();
        dt_start < TOL && dt_end < TOL
    });
    brackets
}

/// Search backward from `search_period.end` to `search_period.start` for a sign change in `f(t) − threshold`.
fn search_crossing_backward<V, F>(
    search_period: Interval<ModifiedJulianDate>,
    f: &F,
    threshold: Quantity<V>,
) -> Option<Interval<ModifiedJulianDate>>
where
    V: Unit,
    F: Fn(ModifiedJulianDate) -> Quantity<V>,
{
    let g = |t: Mjd| f(t) - threshold;
    let range = search_period.end.raw() - search_period.start.raw();

    // Expanding search: start near the extremum and step backward
    let mut bracket = Interval::new(search_period.end, search_period.end);
    let mut g_hi = g(bracket.end);
    let mut step = range * 0.1;
    if step < Days::new(1e-10) {
        step = range * 0.5;
    }

    bracket.start = {
        let t_prev = crate::time::ModifiedJulianDate::new((bracket.end.raw() - step).value());
        if t_prev.raw() >= search_period.start.raw() {
            t_prev
        } else {
            search_period.start
        }
    };
    let mut g_lo = g(bracket.start);

    while bracket.start > search_period.start {
        if opposite_sign(g_lo, g_hi) {
            return Some(bracket);
        }
        bracket.end = bracket.start;
        g_hi = g_lo;
        bracket.start = {
            let t_prev = crate::time::ModifiedJulianDate::new((bracket.start.raw() - step).value());
            if t_prev.raw() >= search_period.start.raw() {
                t_prev
            } else {
                search_period.start
            }
        };
        g_lo = g(bracket.start);
    }

    // Check the final segment
    if opposite_sign(g_lo, g_hi) {
        Some(bracket)
    } else {
        None
    }
}

/// Search forward from `search_period.start` to `search_period.end` for a sign change.
fn search_crossing_forward<V, F>(
    search_period: Interval<ModifiedJulianDate>,
    f: &F,
    threshold: Quantity<V>,
) -> Option<Interval<ModifiedJulianDate>>
where
    V: Unit,
    F: Fn(ModifiedJulianDate) -> Quantity<V>,
{
    let g = |t: Mjd| f(t) - threshold;
    let range = search_period.end.raw() - search_period.start.raw();

    let mut bracket = Interval::new(search_period.start, search_period.start);
    let mut g_lo = g(bracket.start);
    let mut step = range * 0.1;
    if step < Days::new(1e-10) {
        step = range * 0.5;
    }

    bracket.end = {
        let t_next = crate::time::ModifiedJulianDate::new((bracket.start.raw() + step).value());
        if t_next.raw() <= search_period.end.raw() {
            t_next
        } else {
            search_period.end
        }
    };
    let mut g_hi = g(bracket.end);

    while bracket.end < search_period.end {
        if opposite_sign(g_lo, g_hi) {
            return Some(bracket);
        }
        bracket.start = bracket.end;
        g_lo = g_hi;
        bracket.end = {
            let t_next = crate::time::ModifiedJulianDate::new((bracket.end.raw() + step).value());
            if t_next.raw() <= search_period.end.raw() {
                t_next
            } else {
                search_period.end
            }
        };
        g_hi = g(bracket.end);
    }

    if opposite_sign(g_lo, g_hi) {
        Some(bracket)
    } else {
        None
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    type Radians = Quantity<Radian>;

    fn mjd(v: f64) -> Mjd {
        crate::time::ModifiedJulianDate::new((Days::new(v)).value())
    }
    fn period(a: f64, b: f64) -> Interval<ModifiedJulianDate> {
        Interval::new(mjd(a), mjd(b))
    }

    fn mjd_f64(t: Mjd) -> f64 {
        t.quantity().value()
    }

    #[test]
    fn fixed_step_finds_sine_crossings() {
        // Shift to avoid exact zeros at grid points
        let f = |t: Mjd| Radians::new((2.0 * std::f64::consts::PI * (mjd_f64(t) + 0.03)).sin());
        let brackets =
            fixed_step_brackets(period(0.0, 1.0), Days::new(0.05), &f, Radians::new(0.0));
        // Shifted sin crosses 0 twice in [0,1] (near 0.47 and 0.97)
        assert!(
            brackets.len() >= 2,
            "got {} brackets: {:?}",
            brackets.len(),
            brackets
        );
    }

    #[test]
    fn adaptive_step_finds_crossings() {
        let f = |t: Mjd| Radians::new((2.0 * std::f64::consts::PI * (mjd_f64(t) + 0.03)).sin());
        let brackets = adaptive_step_brackets(
            period(0.0, 1.0),
            Days::new(0.2),
            Days::new(0.01),
            &f,
            Radians::new(0.0),
        );
        assert!(brackets.len() >= 2, "got {} brackets", brackets.len());
    }

    #[test]
    fn extrema_based_finds_satellite_pass() {
        // Simulate a satellite pass: altitude peaks at t=5 with max=30°
        let f = |t: Mjd| Radians::new(-2.0 * (mjd_f64(t) - 5.0).powi(2) + 30.0);
        let brackets =
            extrema_based_brackets(period(0.0, 10.0), Days::new(0.5), &f, Radians::new(10.0));

        // Should find rise and set brackets around the peak
        assert!(brackets.len() >= 2, "got {} brackets", brackets.len());

        // Verify brackets contain actual crossings
        for br in &brackets {
            let g_lo = f(br.start) - Radians::new(10.0);
            let g_hi = f(br.end) - Radians::new(10.0);
            assert!(
                g_lo.signum() * g_hi.signum() <= 0.0,
                "bracket [{}, {}] doesn't contain crossing: g_lo={}, g_hi={}",
                br.start,
                br.end,
                g_lo,
                g_hi
            );
        }
    }

    #[test]
    fn extrema_based_no_pass_above_threshold() {
        // Satellite never reaches threshold
        let f = |t: Mjd| Radians::new(-2.0 * (mjd_f64(t) - 5.0).powi(2) + 5.0);
        let brackets =
            extrema_based_brackets(period(0.0, 10.0), Days::new(0.5), &f, Radians::new(10.0));
        assert!(brackets.is_empty(), "no passes above threshold");
    }

    #[test]
    fn fixed_step_no_crossings() {
        let brackets = fixed_step_brackets(
            period(0.0, 10.0),
            Days::new(1.0),
            &|_: Mjd| Radians::new(5.0),
            Radians::new(0.0),
        );
        assert!(brackets.is_empty());
    }
}