Reorganize projection stuff

PolyhedralCombinatorics.lean

@@ -1,8 +1,16 @@
 -- This module serves as the root of the `PolyhedralCombinatorics` library.
 -- Import modules here that should be built as part of the library.
 import Utils.IsEmpty
+import Utils.Finset
 
 import PolyhedralCombinatorics.LinearSystem.Basic
 import PolyhedralCombinatorics.LinearSystem.LinearConstraints
+
 import PolyhedralCombinatorics.Polyhedron.Basic
-import PolyhedralCombinatorics.Polyhedron.Projection
+import PolyhedralCombinatorics.Polyhedron.Duality
+import PolyhedralCombinatorics.Polyhedron.Notation
+
+import PolyhedralCombinatorics.Projection.Basic
+import PolyhedralCombinatorics.Projection.SemiSpace
+import PolyhedralCombinatorics.Projection.Computable
+import PolyhedralCombinatorics.Projection.FourierMotzkin

PolyhedralCombinatorics/LinearSystem/Basic.lean

@@ -48,8 +48,10 @@ variable {𝔽 m n}
 @[match_pattern] abbrev of (A : Matrix (Fin m) (Fin n) 𝔽) (b : Fin m → 𝔽) : LinearSystem 𝔽 n :=
   ⟨m, A, b⟩
 
+/-- The empty linear system. -/
 @[simp] def empty : LinearSystem 𝔽 n := of vecEmpty vecEmpty
 
+/-- A linear system is empty if and only if it has no rows. -/
 theorem eq_empty_iff : p = empty ↔ p.m = 0 := by
   constructor <;> intro h
   . rw [h]
@@ -132,6 +134,8 @@ instance : HasEquiv (LinearSystem 𝔽 n) := ⟨(·.solutions = ·.solutions)⟩
 instance isSetoid (𝔽 n) [LinearOrderedField 𝔽] : Setoid (LinearSystem 𝔽 n) :=
   ⟨instHasEquiv.Equiv, fun _ ↦ rfl, Eq.symm, Eq.trans⟩
 
+@[simp] theorem equiv_def : p ≈ q ↔ p.solutions = q.solutions := Iff.rfl
+
 @[simp] lemma mem_solutions : x ∈ p.solutions ↔ p.mat *ᵥ x ≤ p.vec := Set.mem_setOf
 
 @[simp] lemma mem_solutions_of {x : Fin n → 𝔽} : x ∈ (of A b).solutions ↔ A *ᵥ x ≤ b := Set.mem_setOf

PolyhedralCombinatorics/Polyhedron/Duality.lean

@@ -1,4 +1,4 @@
-import PolyhedralCombinatorics.Polyhedron.Projection
+import PolyhedralCombinatorics.Projection.FourierMotzkin
 import Mathlib.LinearAlgebra.Matrix.DotProduct
 
 namespace LinearSystem.Duality
@@ -16,5 +16,5 @@ theorem farkas_lemma {A : Matrix (Fin m) (Fin n) 𝔽} {b : Fin m → 𝔽}
     rw [dotProduct_mulVec, hy, zero_dotProduct] at this
     assumption
   . by_contra hc
-    simp_rw [not_exists, ←eq_empty_iff, ←recElimDimensions_eq_empty_iff P(A, b) (le_refl n)] at hc
+    simp_rw [not_exists, ←Polyhedron.eq_empty_iff] at hc
     sorry

PolyhedralCombinatorics/Projection/Basic.lean

@@ -0,0 +1,41 @@
+import PolyhedralCombinatorics.Polyhedron.Basic
+
+variable {𝔽} [lof : LinearOrderedField 𝔽] {n : ℕ}
+
+namespace LinearSystem
+open Matrix Finset
+
+def projection (P : LinearSystem 𝔽 n) (c : Fin n → 𝔽) := { x | ∃ α : 𝔽, x + α • c ∈ P.solutions }
+
+theorem projection_neg (P : LinearSystem 𝔽 n) (c : Fin n → 𝔽)
+  : P.projection (-c) = P.projection c := by
+  simp only [projection, Set.ext_iff, Set.mem_setOf]
+  intro x
+  constructor
+  all_goals (
+    intro h
+    obtain ⟨α, h⟩ := h
+    exists -α
+    simp_all
+  )
+
+theorem projection_concat_comm {P Q : LinearSystem 𝔽 n} {c}
+  : projection (P ++ Q) c = projection (Q ++ P) c := by
+  unfold projection
+  simp_rw [concat_solutions, Set.inter_comm]
+
+@[simp low] lemma mem_projection {P : LinearSystem 𝔽 n}
+  : x ∈ P.projection c ↔ ∃ α : 𝔽, x + α • c ∈ P.solutions := Set.mem_setOf
+
+end LinearSystem
+
+
+namespace Polyhedron
+open Matrix LinearSystem
+
+def projection (P : Polyhedron 𝔽 n) (c : Fin n → 𝔽) := {x | ∃ α : 𝔽, x + α • c ∈ P}
+
+@[simp] lemma mem_projection {P : Polyhedron 𝔽 n}
+  : x ∈ P.projection c ↔ ∃ α : 𝔽, x + α • c ∈ P := Set.mem_setOf
+
+end Polyhedron

PolyhedralCombinatorics/Projection/Computable.lean

@@ -0,0 +1,122 @@
+import PolyhedralCombinatorics.Projection.SemiSpace
+
+import Mathlib.Data.Finset.Sort
+import Mathlib.Data.Sum.Order
+
+import Mathlib.Tactic.LiftLets
+
+import Utils.IsEmpty
+import Utils.Finset
+
+variable {𝔽} [lof : LinearOrderedField 𝔽] {n : ℕ}
+
+namespace LinearSystem
+open Matrix Finset
+
+def computeProjection (S : LinearSystem 𝔽 n) (c : Fin n → 𝔽) : LinearSystem 𝔽 n :=
+  let N : Finset (Fin S.m) := {i | S.mat i ⬝ᵥ c < 0}
+  let Z : Finset (Fin S.m) := {i | S.mat i ⬝ᵥ c = 0}
+  let P : Finset (Fin S.m) := {i | S.mat i ⬝ᵥ c > 0}
+  let R := Z ⊕ₗ (N ×ₗ P)
+  let r := Fintype.card R
+  let p : Fin r ≃o R := Fintype.orderIsoFinOfCardEq R rfl
+  let D : Matrix (Fin r) (Fin n) 𝔽 := Matrix.of fun i ↦
+    match p i with
+    | .inl s => S.mat s
+    | .inr (s, t) => (S.mat t ⬝ᵥ c) • S.mat s - (S.mat s ⬝ᵥ c) • S.mat t
+  let d : Fin r → 𝔽 := fun i ↦
+    match p i with
+    | .inl s => S.vec s
+    | .inr (s, t) => (S.mat t ⬝ᵥ c) • S.vec s - (S.mat s ⬝ᵥ c) • S.vec t
+  of D d
+
+@[simp] theorem mem_computeProjection {S : LinearSystem 𝔽 n} {c} {x}
+  : x ∈ (computeProjection S c).solutions ↔ x ∈ S.projection c := by
+  let A := S.mat
+  let b := S.vec
+  unfold computeProjection
+  lift_lets
+  extract_lets _ _ _ _ r p D d
+  rw [projection_inter_pairs']
+  show (∀ (i : Fin r), D i ⬝ᵥ x ≤ d i) ↔ _
+  constructor
+  . intro h
+    constructor
+    . intro i hi
+      rw [mem_semiSpace]
+      have := h $ p.symm $ Sum.inl ⟨i, mem_filter_univ.mpr hi⟩
+      simp only [D, d, dotProduct, Matrix.of_apply, OrderIso.apply_symm_apply] at this
+      exact this
+    . intro i j hi hj
+      have := h $ p.symm $ Sum.inr ⟨⟨i, mem_filter_univ.mpr hi⟩, ⟨j, mem_filter_univ.mpr hj⟩⟩
+      simp only [D, d, dotProduct, Matrix.of_apply, OrderIso.apply_symm_apply] at this
+      rw [proj_concat_semiSpace_nonorthogonal_neg_pos hi hj, mem_semiSpace]
+      exact this
+  . intro ⟨hz, hnp⟩ i
+    rcases hi : p i with s | ⟨s, t⟩
+    . have hD : D i = A s := by simp only [D, funext_iff, of_apply, hi, implies_true]
+      have hd : d i = b s := by simp only [d, hi]
+      have hs := Finset.mem_filter_univ.mp s.property
+      replace := hz s hs
+      rw [mem_semiSpace] at this
+      rw [hD, hd]
+      exact this
+    . have hD : D i = ((A t ⬝ᵥ c) • A s - (A s ⬝ᵥ c) • A t) := by simp_all only [D, funext_iff, of_apply, implies_true]
+      have hd : d i = (A t ⬝ᵥ c) • b s - (A s ⬝ᵥ c) • b t := by simp_all only [d]
+      have hs := mem_filter_univ.mp s.property
+      have ht := mem_filter_univ.mp t.property
+      replace := hnp s t hs ht
+      rw [proj_concat_semiSpace_nonorthogonal_neg_pos hs ht, mem_semiSpace] at this
+      rw [hD, hd]
+      exact this
+
+@[simp] lemma solutions_computeProjection (S : LinearSystem 𝔽 n) (c : Fin n → 𝔽)
+  : (computeProjection S c).solutions = S.projection c := by
+  simp_rw [Set.ext_iff, mem_projection]
+  apply mem_computeProjection
+
+end LinearSystem
+
+namespace Polyhedron
+
+def computeProjection (P : Polyhedron 𝔽 n) (c : Fin n → 𝔽) : Polyhedron 𝔽 n :=
+  Quotient.liftOn P (fun S ↦ S.computeProjection c) fun a b h ↦ by
+    apply toSet_inj.mp
+    simp_rw [toSet_ofLinearSystem, LinearSystem.solutions_computeProjection, LinearSystem.projection]
+    rw [h]
+
+@[simp]
+theorem mem_computeProjection {P : Polyhedron 𝔽 n} {c : Fin n → 𝔽}
+  : x ∈ P.computeProjection c ↔ ∃ α : 𝔽, x + α • c ∈ P := by
+  induction P using Quotient.ind
+  unfold computeProjection
+  simp_rw [Quotient.lift_mk, mem_ofLinearSystem, LinearSystem.solutions_computeProjection]
+  rfl
+
+@[simp] theorem subset_computeProjection {P : Polyhedron 𝔽 n} {c : Fin n → 𝔽}
+  : P ⊆ P.computeProjection c := by
+  intro x
+  rw [mem_computeProjection]
+  intro x_mem_S
+  exists 0
+  rw [zero_smul, add_zero]
+  assumption
+
+theorem self_inter_computeProjection (P : Polyhedron 𝔽 n) (c) : P ∩ P.computeProjection c = P := by
+  apply subset_antisymm
+  . apply inter_subset_left
+  . apply subset_inter
+    . apply subset_refl
+    . apply subset_computeProjection
+
+@[simp] theorem computeProjection_eq_empty_iff (P : Polyhedron 𝔽 n) (c)
+  : P.computeProjection c = ∅ ↔ P = ∅ := by
+  constructor <;> intro h
+  . have := self_inter_computeProjection P c
+    simp_rw [Polyhedron.ext_iff, mem_inter, mem_computeProjection] at this
+    simp_rw [eq_empty_iff, mem_computeProjection] at h ⊢
+    simp_all
+  . simp_rw [eq_empty_iff, mem_computeProjection]
+    simp_all [not_mem_empty]
+
+end Polyhedron

PolyhedralCombinatorics/Projection/FourierMotzkin.lean

@@ -0,0 +1,73 @@
+import PolyhedralCombinatorics.Projection.Computable
+
+import PolyhedralCombinatorics.Polyhedron.Notation
+
+variable {𝔽} [lof : LinearOrderedField 𝔽] {n : ℕ}
+
+namespace LinearSystem
+open Matrix
+
+def elim0 (S : LinearSystem 𝔽 (n + 1)) : LinearSystem 𝔽 n :=
+  LinearSystem.of (vecTail ∘ S.mat) S.vec
+
+theorem mem_elim0_of_vecCons_zero_mem {S : LinearSystem 𝔽 (n + 1)} {x : Fin n → 𝔽}
+  (h : vecCons 0 x ∈ S.solutions) : x ∈ S.elim0.solutions := by
+  simp only [elim0, mem_solutions, Pi.le_def, mulVec, of_apply] at h ⊢
+  intro i
+  replace h := h i
+  rw [dotProduct_cons, mul_zero, zero_add] at h
+  exact h
+
+end LinearSystem
+
+namespace Polyhedron
+open Matrix
+
+def fourierMotzkin (P : Polyhedron 𝔽 n) (j : Fin n) :=
+  !P{x_[j] = 0} ∩ P.computeProjection x_[j]
+
+theorem mem_fourierMotzkin (P : Polyhedron 𝔽 n) (j : Fin n) :
+  x ∈ P.fourierMotzkin j ↔ x j = 0 ∧ ∃ (α : 𝔽), x + Pi.single j α ∈ P := by
+  simp_rw [
+    fourierMotzkin, mem_inter, mem_computeProjection, mem_ofConstraints,
+    List.mem_singleton, forall_eq, LinearConstraint.valid,
+    single_dotProduct, one_mul, and_congr_right_iff,
+    ←Pi.single_smul, smul_eq_mul, mul_one, implies_true
+  ]
+
+-- def elimSucc (P : Polyhedron 𝔽 (n + 1)) : Polyhedron 𝔽 n :=
+--   Quotient.recOn P
+--     (fun S ↦ LinearSystem.of (Matrix.of fun i j ↦ S.mat i j.castSucc) S.vec)
+--     fun S S' h ↦ by
+--       ext x
+--       simp only [eq_rec_constant, mem_ofLinearSystem_of]
+--       simp_rw [LinearSystem.equiv_def, Set.ext_iff, LinearSystem.mem_solutions] at h
+--       replace h := h $ fun i ↦ if h : i.val < n then x (i.castLT h) else 0
+--       simp_rw [Pi.le_def, mulVec] at h
+--       sorry
+
+-- def recElimDimensions (p : Polyhedron 𝔽 n) {k : ℕ} (h : k ≤ n) :=
+--   match k with
+--   | 0 => p
+--   | k + 1 => (p.recElimDimensions $ le_of_add_le_left h).computeProjection x_[⟨k, h⟩]
+
+-- @[simp] theorem recElimDimensions_eq_empty_iff (p : Polyhedron 𝔽 n) {k : ℕ} (hk : k ≤ n)
+--   : p.recElimDimensions hk = ∅ ↔ p = ∅ := by
+--   unfold recElimDimensions
+--   split
+--   . rfl
+--   . rw [computeProjection_eq_empty_iff]
+--     apply p.recElimDimensions_eq_empty_iff
+
+-- theorem recElimDimensions_all_empty_or_univ (p : Polyhedron 𝔽 n) {k : ℕ}
+--   : let p' := p.recElimDimensions (le_refl _); p' = ∅ ∨ p' = ⊤ := by
+--   unfold recElimDimensions
+--   split
+--   . by_cases finZeroElim ∈ p <;> simp_all [Polyhedron.ext_iff, not_mem_empty, mem_univ]
+--   . simp only [Nat.succ_eq_add_one, computeProjection_eq_empty_iff, recElimDimensions_eq_empty_iff]
+--     by_cases p = ∅
+--     . left; assumption
+--     . right
+--       simp_rw [Polyhedron.ext_iff, mem_computeProjection, mem_univ, iff_true]
+--       intro x
+--       sorry