-- Orthogonal projection onto a line and its least-squares characterization
-- Part of Cumulus, the Grace17 mathematics library.
-- Copyright (c) 2026 Grant Passmore.
-- Released under the Apache License 2.0. See LICENSE in the repository root.
-- Authors: Grant Passmore

module Grace.Cumulus.LinearAlgebra.Projection
import CauchySchwarz

-- All vector statements concern the first n coordinates of real sequences.
-- The inverse is total; the projection laws below require a nonzero norm.
def projectionCoefficient (u v : RealSequence, n : Nat) : Real :=
  dotProduct(u, v, n) * realInvValue(squaredNorm(v, n))

def lineProjection (u v : RealSequence, n : Nat) : RealSequence :=
  mulSequence(constantSequence(projectionCoefficient(u, v, n)), v)

def lineDistanceSquared (u v : RealSequence, C : Real, n : Nat) : Real :=
  squaredNorm(subSequence(u, mulSequence(constantSequence(C), v)), n)

theorem projectionCoefficient.real (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows Real(projectionCoefficient(u, v, n))
proof
  unfold projectionCoefficient
  typecheck

theorem lineProjection.sequence (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows RealSequence(lineProjection(u, v, n))
proof
  have constants := constantSequence.sequence at projectionCoefficient(u, v, n) with *
  unfold lineProjection
  exact mulSequence.sequence with *

theorem lineDistanceSquared.real (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows Real(lineDistanceSquared(u, v, C, n))
proof
  have constants := constantSequence.sequence with *
  have scaled := mulSequence.sequence at constantSequence(C), v with *
  have difference := subSequence.sequence at u, mulSequence(constantSequence(C), v) with *
  unfold lineDistanceSquared
  exact squaredNorm.real with *

theorem projectionCoefficient.normal_equation (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows projectionCoefficient(u, v, n) * squaredNorm(v, n) = dotProduct(u, v, n)
proof
  have inverseLaw := realInvValue.mul_right with *
  unfold projectionCoefficient
  linear_combination dotProduct(u, v, n) * inverseLaw

theorem lineProjection.orthogonal (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows Orthogonal(subSequence(u, lineProjection(u, v, n)), v, n)
proof
  have projection := lineProjection.sequence at u, v, n with *
  have normal := projectionCoefficient.normal_equation at u, v, n with *
  unfold Orthogonal
  rw dotProduct.sub_left
  unfold lineProjection
  rw dotProduct.scale_left
  unfold squaredNorm at normal
  linear_combination -normal

theorem lineDistanceSquared.expansion (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows lineDistanceSquared(u, v, C, n) =
    squaredNorm(u, n) - (2:R) * C * dotProduct(u, v, n) + C * C * squaredNorm(v, n)
proof
  have constants := constantSequence.sequence with *
  have scaled := mulSequence.sequence at constantSequence(C), v with *
  unfold lineDistanceSquared
  simp [squaredNorm.sub, squaredNorm.scale, dotProduct.scale_right]
  ring

theorem lineProjection.distance_decomposition (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows lineDistanceSquared(u, v, C, n) =
    lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) +
    (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)) * squaredNorm(v, n)
proof
  have normal := projectionCoefficient.normal_equation at u, v, n with *
  rw [lineDistanceSquared.expansion, lineDistanceSquared.expansion]
  linear_combination ((2:R) * C - (2:R) * projectionCoefficient(u, v, n)) * normal

theorem lineProjection.minimum (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) ≤
    lineDistanceSquared(u, v, C, n)
proof
  have squareNonnegative := realMulValue.square_nonnegative at C - projectionCoefficient(u, v, n) with *
  have normNonnegative := squaredNorm.nonnegative at v, n with *
  have gapNonnegative := realMulValue.nonnegative at (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)), squaredNorm(v, n) with *
  have translated := realLe.add_left at realZero, (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)) * squaredNorm(v, n), lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) with *
  have decomposition := lineProjection.distance_decomposition at u, v, C, n with *
  rw decomposition
  simpa [realAddValue.zero_right] using translated

theorem lineProjection.strict_minimum (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  requires different : C ≠ projectionCoefficient(u, v, n)
  shows lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) <
    lineDistanceSquared(u, v, C, n)
proof
  have differenceNonzero : C - projectionCoefficient(u, v, n) ≠ 0
  . intro zero
    have equal : C = projectionCoefficient(u, v, n)
    . linear_combination zero
    exact different with equal
  have squarePositive := realMulValue.square_positive at C - projectionCoefficient(u, v, n) with *
  have normNonnegative := squaredNorm.nonnegative at v, n with *
  have normPositive := Real.positive_of_nonnegative_of_ne_zero at squaredNorm(v, n) with *
  have gapPositive := realLess.mul_right_of_positive at realZero, (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)), squaredNorm(v, n) with *
  simp [realMulValue.zero_left] at gapPositive
  have translated := realLess.add_left at realZero, (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)) * squaredNorm(v, n), lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) with *
  have decomposition := lineProjection.distance_decomposition at u, v, C, n with *
  rw decomposition
  simpa [realAddValue.zero_right] using translated

theorem lineProjection.unique_minimizer (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  requires minimizing : lineDistanceSquared(u, v, C, n) ≤
    lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n)
  shows C = projectionCoefficient(u, v, n)
proof
  by_contra different
  have strict := lineProjection.strict_minimum with *
  have candidate := lineDistanceSquared.real at u, v, C, n with *
  have cycle := realLe.transitive_less with *
  have member : lineDistanceSquared(u, v, C, n) ∈ ℝ
  . settle with Real
  have irreflexive := realLess.irreflexive_on
  exact irreflexive with lineDistanceSquared(u, v, C, n), cycle

-- One orthogonalization step: retain the component perpendicular to v.
def orthogonalResidual (u v : RealSequence, n : Nat) : RealSequence :=
  subSequence(u, lineProjection(u, v, n))

theorem orthogonalResidual.sequence (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows RealSequence(orthogonalResidual(u, v, n))
proof
  have projection := lineProjection.sequence at u, v, n with *
  unfold orthogonalResidual
  exact subSequence.sequence with *

theorem orthogonalResidual.orthogonal (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows Orthogonal(orthogonalResidual(u, v, n), v, n)
proof
  unfold orthogonalResidual
  exact lineProjection.orthogonal with *

theorem orthogonalResidual.reconstruction (u v : Set, n i : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows u[i] = (lineProjection(u, v, n)[i]:R) + (orthogonalResidual(u, v, n)[i]:R)
proof
  have projection := lineProjection.sequence at u, v, n with *
  have ui := RealSequence.value at u, i with *
  have pi := RealSequence.value at lineProjection(u, v, n), i with *
  unfold orthogonalResidual
  rw subSequence.value
  ring

theorem projectionCoefficient.on_line (v : Set, C : Real, n : Nat)
  requires sequence : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows projectionCoefficient(mulSequence(constantSequence(C), v), v, n) = C
proof
  have constants := constantSequence.sequence with *
  have scaled := mulSequence.sequence at constantSequence(C), v with *
  have inverseLaw := realInvValue.mul_right at squaredNorm(v, n) with *
  unfold projectionCoefficient
  rw dotProduct.scale_left
  unfold squaredNorm at inverseLaw
  unfold squaredNorm
  linear_combination C * inverseLaw

theorem lineProjection.idempotent (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows lineProjection(lineProjection(u, v, n), v, n) = lineProjection(u, v, n)
proof
  unfold lineProjection
  rw projectionCoefficient.on_line
  rfl

theorem projectionCoefficient.of_orthogonal (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires perpendicular : Orthogonal(u, v, n)
  shows projectionCoefficient(u, v, n) = (0:R)
proof
  unfold Orthogonal at perpendicular
  simp [projectionCoefficient]
  ring

theorem orthogonalResidual.projection_coefficient (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows projectionCoefficient(orthogonalResidual(u, v, n), v, n) = (0:R)
proof
  have residual := orthogonalResidual.sequence at u, v, n with *
  have perpendicular := orthogonalResidual.orthogonal at u, v, n with *
  exact projectionCoefficient.of_orthogonal with *

theorem orthogonalResidual.pythagorean (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires nonzero : squaredNorm(v, n) ≠ (0:R)
  shows squaredNorm(u, n) = squaredNorm(lineProjection(u, v, n), n) +
    squaredNorm(orthogonalResidual(u, v, n), n)
proof
  have projection := lineProjection.sequence at u, v, n with *
  have constants := constantSequence.sequence at projectionCoefficient(u, v, n) with *
  have scaled := mulSequence.sequence at constantSequence(projectionCoefficient(u, v, n)), v with *
  have normal := projectionCoefficient.normal_equation at u, v, n with *
  unfold orthogonalResidual
  rw squaredNorm.sub
  unfold lineProjection
  simp [squaredNorm.scale, dotProduct.scale_right]
  linear_combination -(2:R) * projectionCoefficient(u, v, n) * normal

theorem projectionCoefficient.add (u w v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceW : RealSequence(w)
  requires sequenceV : RealSequence(v)
  shows projectionCoefficient(addSequence(u, w), v, n) =
    projectionCoefficient(u, v, n) + projectionCoefficient(w, v, n)
proof
  have sum := addSequence.sequence at u, w with *
  unfold projectionCoefficient
  rw dotProduct.add_left
  ring

theorem projectionCoefficient.scale (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows projectionCoefficient(mulSequence(constantSequence(C), u), v, n) =
    C * projectionCoefficient(u, v, n)
proof
  have constants := constantSequence.sequence with *
  have scaled := mulSequence.sequence at constantSequence(C), u with *
  unfold projectionCoefficient
  rw dotProduct.scale_left
  ring