import Mathlib.LinearAlgebra.BilinearForm.Properties
import Mathlib.Data.Real.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.Positivity

/-!
Hildreth (1957), DOI 10.1002/nav.3800040113.
Modern dual-coordinate form of the method, minimizing a quadratic over the
nonnegative orthant. Proves exact coordinate minimization and finite-sweep
descent. No infinite-iteration convergence or primal recovery is claimed.
-/
namespace PaperToLean.Hildreth
noncomputable section
open LinearMap (BilinForm)
variable {n : ℕ}
abbrev Vec (n : ℕ) := Fin n → ℝ

def basis (k : Fin n) : Vec n := fun i => if i = k then 1 else 0
def move (x : Vec n) (k : Fin n) (t : ℝ) : Vec n :=
  x + (t - x k) • basis k
def objective (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ) (x : Vec n) : ℝ :=
  B x x / 2 - l x
def coordinate (u g d : ℝ) : ℝ := max 0 (u - g / d)
def step (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ) (k : Fin n) (x : Vec n) : Vec n :=
  move x k (coordinate (x k) (B x (basis k) - l (basis k)) (B (basis k) (basis k)))
def sweep (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ) : List (Fin n) → Vec n → Vec n
  | [], x => x
  | k :: ks, x => sweep B l ks (step B l k x)

theorem move_at (x : Vec n) (k : Fin n) (t : ℝ) : move x k t k = t := by
  simp [move, basis]

theorem move_other (x : Vec n) (k i : Fin n) (t : ℝ) (hi : i ≠ k) :
    move x k t i = x i := by
  simp [move, basis, hi]

theorem objective_change (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ)
    (hB : B.IsSymm) (x : Vec n) (k : Fin n) (t : ℝ) :
    objective B l (move x k t) = objective B l x +
      (B x (basis k) - l (basis k)) * (t - x k) +
      B (basis k) (basis k) / 2 * (t - x k) ^ 2 := by
  simp only [objective, move, map_add, LinearMap.add_apply,
    map_smul, LinearMap.smul_apply, smul_eq_mul]
  rw [hB.eq (basis k) x]
  ring

/-- Clipping the unconstrained minimizer at zero solves the subproblem.
The gap bound also proves uniqueness, since d is strictly positive. -/
theorem coordinate_gap (u g d t : ℝ) (hd : 0 < d) (ht : 0 ≤ t) :
    let v := coordinate u g d
    g * (v - u) + d / 2 * (v - u) ^ 2 + d / 2 * (t - v) ^ 2 ≤
      g * (t - u) + d / 2 * (t - u) ^ 2 := by
  dsimp only
  have hd0 : d ≠ 0 := ne_of_gt hd
  by_cases h : 0 ≤ u - g / d
  · rw [coordinate, max_eq_right h]
    have heq : g * (u - g / d - u) + d / 2 * (u - g / d - u) ^ 2 +
        d / 2 * (t - (u - g / d)) ^ 2 = g * (t - u) + d / 2 * (t - u) ^ 2 := by
      field_simp
      ring
    exact le_of_eq heq
  · rw [coordinate, max_eq_left (le_of_not_ge h)]
    have hg : d * u ≤ g := by
      have hh : u ≤ g / d := by linarith
      have := (le_div_iff₀ hd).mp hh
      nlinarith
    nlinarith [mul_nonneg (sub_nonneg.mpr hg) ht]

theorem step_feasible (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ)
    (k : Fin n) (x : Vec n) (hx : ∀ i, 0 ≤ x i) : ∀ i, 0 ≤ step B l k x i := by
  intro i
  by_cases hi : i = k
  · subst i
    unfold step
    rw [move_at]
    exact le_max_left _ _
  · unfold step
    rw [move_other x k i _ hi]
    exact hx i

theorem coordinate_optimal (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ)
    (hB : B.IsSymm) (x : Vec n) (k : Fin n)
    (hd : 0 < B (basis k) (basis k)) (t : ℝ) (ht : 0 ≤ t) :
    objective B l (step B l k x) ≤ objective B l (move x k t) := by
  unfold step
  rw [objective_change B l hB, objective_change B l hB]
  have hg := coordinate_gap (x k) (B x (basis k) - l (basis k))
    (B (basis k) (basis k)) t hd ht
  have hp : 0 ≤ B (basis k) (basis k) / 2 *
      (t - coordinate (x k) (B x (basis k) - l (basis k)) (B (basis k) (basis k))) ^ 2 := by
    positivity
  dsimp only at hg
  linarith

theorem step_descent (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ)
    (hB : B.IsSymm) (x : Vec n) (k : Fin n)
    (hd : 0 < B (basis k) (basis k)) (hx : 0 ≤ x k) :
    objective B l (step B l k x) ≤ objective B l x := by
  have h := coordinate_optimal B l hB x k hd (x k) hx
  simpa [move] using h

/-- Every finite schedule, including repeated coordinates, remains feasible
and cannot increase the objective. This makes no fairness assumption. -/
theorem sweep_guarantee (B : BilinForm ℝ (Vec n)) (l : Vec n →ₗ[ℝ] ℝ)
    (hB : B.IsSymm) (hd : ∀ k, 0 < B (basis k) (basis k))
    (ks : List (Fin n)) (x : Vec n) (hx : ∀ i, 0 ≤ x i) :
    (∀ i, 0 ≤ sweep B l ks x i) ∧ objective B l (sweep B l ks x) ≤ objective B l x := by
  induction ks generalizing x with
  | nil => exact ⟨hx, le_rfl⟩
  | cons k ks ih =>
    obtain ⟨hf, ho⟩ := ih (step B l k x) (step_feasible B l k x hx)
    exact ⟨hf, ho.trans (step_descent B l hB x k (hd k) (hx k))⟩

#print axioms move_at
#print axioms move_other
#print axioms objective_change
#print axioms coordinate_gap
#print axioms step_feasible
#print axioms coordinate_optimal
#print axioms step_descent
#print axioms sweep_guarantee
end
end PaperToLean.Hildreth