import Mathlib

namespace Gam.Terms.Basis

section BSplineFoundations

variable {numKnots : ℕ}

/-- A valid B-spline knot vector: non-decreasing with proper multiplicity. -/
structure KnotVector (m : ℕ) where
  knots : Fin m → ℝ
  sorted : ∀ i j : Fin m, i ≤ j → knots i ≤ knots j

/-- Cox-de Boor recursive definition of B-spline basis function.
    N_{i,p}(x) is the i-th basis function of degree p.
    We use a simpler formulation to avoid index bound issues. -/
noncomputable def bspline_basis_raw (t : ℕ → ℝ) : ℕ → ℕ → ℝ → ℝ
  | i, 0, x => if t i ≤ x ∧ x < t (i + 1) then 1 else 0
  | i, p + 1, x =>
    let left_denom := t (i + p + 1) - t i
    let right_denom := t (i + p + 2) - t (i + 1)
    let left := if left_denom = 0 then 0
                else (x - t i) / left_denom * bspline_basis_raw t i p x
    let right := if right_denom = 0 then 0
                 else (t (i + p + 2) - x) / right_denom * bspline_basis_raw t (i + 1) p x
    left + right

/-- Local support property: N_{i,p}(x) = 0 outside [t_i, t_{i+p+1}).

    **Geometric insight**: The support of a B-spline grows by one knot interval with
    each degree increase. N_{i,p} lives on [t_i, t_{i+p+1}). In the recursion for p+1,
    we combine N_{i,p} (starts at t_i) and N_{i+1,p} (ends at t_{i+p+2}).
    The union creates the new support [t_i, t_{i+p+2}).

    **Proof by induction on p**:
    - Base case (p=0): By definition, N_{i,0}(x) = 1 if t_i ≤ x < t_{i+1}, else 0.
    - Inductive case (p+1): Cox-de Boor recursion combines N_{i,p} and N_{i+1,p}.
      Both have zero support outside the required interval. -/
theorem bspline_local_support (t : ℕ → ℝ)
    (h_sorted : ∀ i j, i ≤ j → t i ≤ t j)
    (i p : ℕ) (x : ℝ)
    (h_outside : x < t i ∨ t (i + p + 1) ≤ x) :
    bspline_basis_raw t i p x = 0 := by
  induction p generalizing i with
  | zero =>
    simp only [bspline_basis_raw]
    split_ifs with h_in
    · obtain ⟨h_lo, h_hi⟩ := h_in
      rcases h_outside with h_lt | h_ge
      · exact absurd h_lo (not_le.mpr h_lt)
      · simp only [add_zero] at h_ge
        exact absurd h_hi (not_lt.mpr h_ge)
    · rfl
  | succ p ih =>
    simp only [bspline_basis_raw]
    rcases h_outside with h_lt | h_ge
    · -- x < t_i: both terms zero
      have h_left_zero : bspline_basis_raw t i p x = 0 := ih i (Or.inl h_lt)
      have h_i1_le : t i ≤ t (i + 1) := h_sorted i (i + 1) (Nat.le_succ i)
      have h_right_zero : bspline_basis_raw t (i + 1) p x = 0 :=
        ih (i + 1) (Or.inl (lt_of_lt_of_le h_lt h_i1_le))
      simp only [h_left_zero, h_right_zero, mul_zero, ite_self, add_zero]
    · -- x ≥ t_{i+p+2}: both terms zero
      have h_right_idx : i + 1 + p + 1 = i + p + 2 := by ring
      have h_right_zero : bspline_basis_raw t (i + 1) p x = 0 := by
        apply ih (i + 1); right; rw [h_right_idx]; exact h_ge
      have h_mono : t (i + p + 1) ≤ t (i + p + 2) := h_sorted (i + p + 1) (i + p + 2) (Nat.le_succ _)
      have h_left_zero : bspline_basis_raw t i p x = 0 := by
        apply ih i; right; exact le_trans h_mono h_ge
      simp only [h_left_zero, h_right_zero, mul_zero, ite_self, add_zero]

/-- B-spline basis functions are non-negative everywhere.

    **Geometry of the "Zero-Out" Property** (Key insight from user):
    The Cox-de Boor recursion uses linear weights: (x - t_i) / (t_{i+p+1} - t_i).
    These weights become NEGATIVE when x < t_i. The ONLY reason the spline remains
    non-negative is that the lower-order basis function N_{i,p}(x) "turns off"
    (becomes exactly zero) precisely when the weight becomes negative.

    Therefore, bspline_local_support is a strict prerequisite for this proof.

    **Proof by induction on p**:
    - Base case (p=0): N_{i,0}(x) is either 0 or 1, both ≥ 0.
    - Inductive case (p+1): For each term α(x) * N_{i,p}(x):
      * If x ∈ [t_i, t_{i+p+1}): α(x) ≥ 0 and N_{i,p}(x) ≥ 0 by IH
      * If x ∉ [t_i, t_{i+p+1}): N_{i,p}(x) = 0 by local_support, so product = 0 -/
theorem bspline_nonneg (t : ℕ → ℝ) (h_sorted : ∀ i j, i ≤ j → t i ≤ t j)
    (i p : ℕ) (x : ℝ) : 0 ≤ bspline_basis_raw t i p x := by
  induction p generalizing i with
  | zero =>
    simp only [bspline_basis_raw]
    split_ifs
    · exact zero_le_one
    · exact le_refl 0
  | succ p ih =>
    simp only [bspline_basis_raw]
    apply add_nonneg
    · -- Left term: (x - t_i) / (t_{i+p+1} - t_i) * N_{i,p}(x)
      split_ifs with h_denom
      · exact le_refl 0
      · by_cases h_in_support : x < t i
        · -- x < t_i: N_{i,p}(x) = 0 by local support, so product = 0
          have : bspline_basis_raw t i p x = 0 :=
            bspline_local_support t h_sorted i p x (Or.inl h_in_support)
          simp only [this, mul_zero, le_refl]
        · -- x ≥ t_i: weight (x - t_i)/denom ≥ 0, and N_{i,p}(x) ≥ 0 by IH
          push_neg at h_in_support
          have h_num_nn : 0 ≤ x - t i := sub_nonneg.mpr h_in_support
          have h_denom_pos : 0 < t (i + p + 1) - t i := by
            have h_le : t i ≤ t (i + p + 1) := h_sorted i (i + p + 1) (by omega)
            exact lt_of_le_of_ne (sub_nonneg.mpr h_le) (ne_comm.mp h_denom)
          exact mul_nonneg (div_nonneg h_num_nn (le_of_lt h_denom_pos)) (ih i)
    · -- Right term: (t_{i+p+2} - x) / (t_{i+p+2} - t_{i+1}) * N_{i+1,p}(x)
      split_ifs with h_denom
      · exact le_refl 0
      · by_cases h_in_support : t (i + p + 2) ≤ x
        · -- x ≥ t_{i+p+2}: N_{i+1,p}(x) = 0 by local support
          have h_idx : i + 1 + p + 1 = i + p + 2 := by ring
          have : bspline_basis_raw t (i + 1) p x = 0 := by
            apply bspline_local_support t h_sorted (i + 1) p x; right; rw [h_idx]; exact h_in_support
          simp only [this, mul_zero, le_refl]
        · -- x < t_{i+p+2}: weight (t_{i+p+2} - x)/denom ≥ 0, and N_{i+1,p}(x) ≥ 0 by IH
          push_neg at h_in_support
          have h_num_nn : 0 ≤ t (i + p + 2) - x := sub_nonneg.mpr (le_of_lt h_in_support)
          have h_denom_pos : 0 < t (i + p + 2) - t (i + 1) := by
            have h_le : t (i + 1) ≤ t (i + p + 2) := h_sorted (i + 1) (i + p + 2) (by omega)
            exact lt_of_le_of_ne (sub_nonneg.mpr h_le) (ne_comm.mp h_denom)
          exact mul_nonneg (div_nonneg h_num_nn (le_of_lt h_denom_pos)) (ih (i + 1))

end BSplineFoundations

end Gam.Terms.Basis
