import Mathlib.Algebra.Order.BigOperators.Ring.Finset
import Mathlib.Data.Real.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.Positivity

/-! Newman (1960), Kantorovich's inequality: the finite weighted scalar case.
No diagonalization or operator functional calculus is hidden in the proof.
The endpoint equality example shows the universal constant is sharp. -/
namespace PaperToLean.Kantorovich
noncomputable section
open scoped BigOperators
variable {ι : Type*}

theorem scalar_bound (m M x : ℝ) (hm : 0 < m) (hlo : m ≤ x) (hhi : x ≤ M) :
    x + m * M / x ≤ m + M := by
  have hx : 0 < x := lt_of_lt_of_le hm hlo
  apply (le_of_sub_nonneg ?_)
  have hp := mul_nonneg (sub_nonneg.mpr hlo) (sub_nonneg.mpr hhi)
  have hi : m + M - (x + m * M / x) = (x - m) * (M - x) / x := by
    field_simp
    ring
  rw [hi]
  exact div_nonneg hp (le_of_lt hx)

theorem weighted_bound (I : Finset ι) (w x : ι → ℝ) (m M : ℝ)
    (hm : 0 < m) (hw : ∀ i ∈ I, 0 ≤ w i) (hs : ∑ i ∈ I, w i = 1)
    (hlo : ∀ i ∈ I, m ≤ x i) (hhi : ∀ i ∈ I, x i ≤ M) :
    (∑ i ∈ I, w i * x i) + m * M * (∑ i ∈ I, w i / x i) ≤ m + M := by
  have h := Finset.sum_le_sum (fun i hi =>
    mul_le_mul_of_nonneg_left (scalar_bound m M (x i) hm (hlo i hi) (hhi i hi)) (hw i hi))
  have heq : (∑ i ∈ I, w i * (x i + m * M / x i)) =
      (∑ i ∈ I, w i * x i) + m * M * (∑ i ∈ I, w i / x i) := by
    simp only [Finset.mul_sum, ← Finset.sum_add_distrib]
    apply Finset.sum_congr rfl
    intro i _
    ring
  rw [heq, ← Finset.sum_mul, hs, one_mul] at h
  exact h

theorem product_from_linear_bound (A C m M : ℝ) (hm : 0 < m) (hM : 0 < M)
    (hA : 0 ≤ A) (hC : 0 ≤ C) (h : A + m * M * C ≤ m + M) :
    A * C ≤ (m + M) ^ 2 / (4 * m * M) := by
  apply (le_div_iff₀ (by positivity : 0 < 4 * m * M)).mpr
  have hsum : 0 ≤ A + m * M * C := by positivity
  have hsq : (A + m * M * C) ^ 2 ≤ (m + M) ^ 2 := by nlinarith
  nlinarith [sq_nonneg (A - m * M * C)]

theorem inequality (I : Finset ι) (w x : ι → ℝ) (m M : ℝ)
    (hm : 0 < m) (hM : 0 < M) (hw : ∀ i ∈ I, 0 ≤ w i)
    (hs : ∑ i ∈ I, w i = 1) (hlo : ∀ i ∈ I, m ≤ x i) (hhi : ∀ i ∈ I, x i ≤ M) :
    (∑ i ∈ I, w i * x i) * (∑ i ∈ I, w i / x i) ≤ (m + M) ^ 2 / (4 * m * M) := by
  apply product_from_linear_bound _ _ m M hm hM
  · exact Finset.sum_nonneg (fun i hi => mul_nonneg (hw i hi) (le_trans (le_of_lt hm) (hlo i hi)))
  · exact Finset.sum_nonneg (fun i hi => div_nonneg (hw i hi) (le_trans (le_of_lt hm) (hlo i hi)))
  · exact weighted_bound I w x m M hm hw hs hlo hhi

theorem endpoint_sharp (m M : ℝ) (hm : m ≠ 0) (hM : M ≠ 0) :
    ((m + M) / 2) * ((1 / m + 1 / M) / 2) = (m + M) ^ 2 / (4 * m * M) := by
  field_simp
  ring

#print axioms scalar_bound
#print axioms weighted_bound
#print axioms product_from_linear_bound
#print axioms inequality
#print axioms endpoint_sharp
end
end PaperToLean.Kantorovich