[Merged by Bors] - feat(Tactic/Linarith): Simplex Algorithm oracle

Mathlib.lean

@@ -3620,6 +3620,10 @@ import Mathlib.Tactic.Linarith.Frontend
 import Mathlib.Tactic.Linarith.Lemmas
 import Mathlib.Tactic.Linarith.Parsing
 import Mathlib.Tactic.Linarith.Preprocessing
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.Datatypes
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.Gauss
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.PositiveVector
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.SimplexAlgorithm
 import Mathlib.Tactic.Linarith.Verification
 import Mathlib.Tactic.LinearCombination
 import Mathlib.Tactic.Lint

Mathlib/Analysis/Calculus/BumpFunction/FiniteDimension.lean

@@ -556,7 +556,6 @@ instance (priority := 100) {E : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E
           _ = 1 - (R - 1) / (R + 1) := by field_simp; ring
       support := fun R hR => by
         have A : 0 < (R + 1) / 2 := by linarith
-        have A' : 0 < R + 1 := by linarith
         have C : (R - 1) / (R + 1) < 1 := by apply (div_lt_one _).2 <;> linarith
         simp only [hR, if_true, support_comp_inv_smul₀ A.ne', y_support _ (IR R hR) C,
           _root_.smul_ball A.ne', Real.norm_of_nonneg A.le, smul_zero]

Mathlib/Tactic.lean

@@ -99,6 +99,10 @@ import Mathlib.Tactic.Linarith.Frontend
 import Mathlib.Tactic.Linarith.Lemmas
 import Mathlib.Tactic.Linarith.Parsing
 import Mathlib.Tactic.Linarith.Preprocessing
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.Datatypes
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.Gauss
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.PositiveVector
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.SimplexAlgorithm
 import Mathlib.Tactic.Linarith.Verification
 import Mathlib.Tactic.LinearCombination
 import Mathlib.Tactic.Lint

Mathlib/Tactic/Linarith/Datatypes.lean

@@ -307,7 +307,8 @@ If successful, it returns a map `coeff : Nat → Nat` as a certificate.
 This map represents that we can find a contradiction by taking the sum `∑ (coeff i) * hyps[i]`.
 
 The default `CertificateOracle` used by `linarith` is
-`Linarith.FourierMotzkin.produceCertificate`.
+`Linarith.SimplexAlgo.produceCertificate`.
+`Linarith.FourierMotzkin.produceCertificate` is also available (though has some bugs).
 -/
 def CertificateOracle : Type :=
   List Comp → Nat → MetaM (Std.HashMap Nat Nat)

Mathlib/Tactic/Linarith/Elimination.lean

@@ -4,6 +4,7 @@ Released under Apache 2.0 license as described in the file LICENSE.
 Authors: Robert Y. Lewis
 -/
 import Mathlib.Tactic.Linarith.Datatypes
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.PositiveVector
 
 /-!
 # The Fourier-Motzkin elimination procedure
@@ -335,4 +336,33 @@ def FourierMotzkin.produceCertificate : CertificateOracle :=
   | (Except.ok _) => failure
   | (Except.error contr) => return contr.src.flatten
 
+namespace SimplexAlgorithm
+
+/-- Preprocess the goal to pass it to `findPositiveVector`. -/
+def preprocess (hyps : List Comp) (maxVar : ℕ) : Matrix (maxVar + 1) (hyps.length) × List Nat :=
+  let mdata : Array (Array ℚ) := Array.ofFn fun i : Fin (maxVar + 1) =>
+    Array.mk <| hyps.map (·.coeffOf i)
+  let strictIndexes : List ℕ := hyps.findIdxs (·.str == Ineq.lt)
+  ⟨⟨mdata⟩, strictIndexes⟩
+
+/-- Extract the certificate from the `vec` found by `findPositiveVector`. -/
+def postprocess (vec : Array ℚ) : HashMap ℕ ℕ :=
+  let common_den : ℕ := vec.foldl (fun acc item => acc.lcm item.den) 1
+  let vecNat : Array ℕ := vec.map (fun x : ℚ => (x * common_den).floor.toNat)
+  HashMap.ofList <| vecNat.toList.enum.filter (fun ⟨_, item⟩ => item != 0)
+
+/--
+`produceCertificate hyps vars` tries to derive a contradiction from the comparisons in `hyps`
+by eliminating all variables ≤ `maxVar`.
+If successful, it returns a map `coeff : ℕ → ℕ` as a certificate.
+This map represents that we can find a contradiction by taking the sum `∑ (coeff i) * hyps[i]`.
+-/
+def produceCertificate : CertificateOracle :=
+  fun hyps maxVar => do
+    let ⟨A, strictIndexes⟩ := preprocess hyps maxVar
+    let vec := findPositiveVector A strictIndexes
+    return postprocess vec
+
+end SimplexAlgorithm
+
 end Linarith

Mathlib/Tactic/Linarith/SimplexAlgorithm/Datatypes.lean

@@ -0,0 +1,48 @@
+/-
+Copyright (c) 2024 Vasily Nesterov. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Vasily Nesterov
+-/
+import Std.Data.Rat.Basic
+
+/-!
+# Datatypes for Simplex Algorithm implementation
+-/
+
+namespace Linarith.SimplexAlgorithm
+
+/--
+Structure for matrices over ℚ.
+
+So far it is just a 2d-array carrying dimensions (that are supposed to match with the actual
+dimensions of `data`), but the plan is to add some `Prop`-data and make the structure strict and
+safe.
+
+Note: we avoid using the `Matrix` from `Mathlib.Data.Matrix` because it is far more efficient to
+store matrix as its entries than as function between `Fin`-s.
+-/
+structure Matrix (n m : Nat) where
+  /-- The content of the matrix. -/
+  data : Array (Array Rat)
+  -- hn_pos : n > 0
+  -- hm_pos : m > 0
+  -- hn : data.size = n
+  -- hm (i : Fin n) : data[i].size = m
+deriving Repr
+
+instance (n m : Nat) : GetElem (Matrix n m) Nat (Array Rat) fun _ i => i < n where
+  getElem mat i _ := mat.data[i]!
+
+/--
+`Table` is a structure Simplex Algorithm operates on. The `i`-th row of `mat` expresses the
+variable `basic[i]` as a linear combination of variables from `free`.
+-/
+structure Table where
+  /-- Array containing the basic variables' indexes -/
+  basic : Array Nat
+  /-- Array containing the free variables' indexes -/
+  free : Array Nat
+  /-- Matrix of coefficients the basic variables expressed through the free ones. -/
+  mat : Matrix basic.size free.size
+
+end Linarith.SimplexAlgorithm

Mathlib/Tactic/Linarith/SimplexAlgorithm/Gauss.lean

@@ -0,0 +1,97 @@
+/-
+Copyright (c) 2024 Vasily Nesterov. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Vasily Nesterov
+-/
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.Datatypes
+
+/-!
+# Gaussian Elimination algorithm
+
+The first step of `Linarith.SimplexAlgorithm.findPositiveVector` is finding initial feasible
+solution which is done by standard Gaussian Elimination algorithm implemented in this file.
+-/
+
+namespace Linarith.SimplexAlgorithm.Gauss
+
+/-- The monad for the Gaussian Elimination algorithm. -/
+abbrev GaussM (n m : Nat) := StateM <| Matrix n m
+
+/-- Finds the first row starting from the current row with nonzero element in current column. -/
+def findNonzeroRow (row col : Nat) {n m : Nat} : GaussM n m <| Option Nat := do
+  for i in [row:n] do
+    if (← get)[i]![col]! != 0 then
+      return i
+  return .none
+
+/-- Swaps two rows. -/
+def swapRows {n m : Nat} (i j : Nat) : GaussM n m Unit := do
+  if i != j then
+    modify fun mat =>
+      let swapped : Matrix n m := ⟨mat.data.swap! i j⟩
+      swapped
+
+/-- Subtracts `i`-th row * `coef` from `j`-th row. -/
+def subtractRow {n m : Nat} (i j : Nat) (coef : Rat) : GaussM n m Unit :=
+  modify fun mat =>
+    let newData : Array (Array Rat) := mat.data.modify j fun row =>
+      row.zipWith mat[i]! fun x y => x - coef * y
+    ⟨newData⟩
+
+/-- Divides row by `coef`. -/
+def divideRow {n m : Nat} (i : Nat) (coef : Rat) : GaussM n m Unit :=
+  modify fun mat =>
+    let newData : Array (Array Rat) := mat.data.modify i (·.map (· / coef))
+    ⟨newData⟩
+
+/-- Implementation of `getTable` in `GaussM` monad. -/
+def getTableImp {n m : Nat} : GaussM n m Table := do
+  let mut free : Array Nat := #[]
+  let mut basic : Array Nat := #[]
+
+  let mut row : Nat := 0
+  let mut col : Nat := 0
+
+  while row < n && col < m do
+    match ← findNonzeroRow row col with
+    | .none =>
+      free := free.push col
+      col := col + 1
+      continue
+    | .some rowToSwap =>
+      swapRows row rowToSwap
+
+    divideRow row (← get)[row]![col]!
+
+    for i in [:n] do
+      if i == row then
+        continue
+      let coef := (← get)[i]![col]!
+      subtractRow row i coef
+
+    basic := basic.push col
+    row := row + 1
+    col := col + 1
+
+  for i in [col:m] do
+    free := free.push i
+
+  let ansData : Array (Array Rat) := ← do
+    let mat := (← get)
+    return Array.ofFn (fun row : Fin row => free.map fun f => -mat[row]![f]!)
+
+  return {
+    free := free
+    basic := basic
+    mat := ⟨ansData⟩
+  }
+
+/--
+Given matrix `A`, solves the linear equation `A x = 0` and returns the solution as a table where
+some variables are free and others (basic) variable are expressed as linear combinations of the free
+ones.
+-/
+def getTable {n m : Nat} (A : Matrix n m) : Table := Id.run do
+  return (← getTableImp.run A).fst
+
+end Linarith.SimplexAlgorithm.Gauss

Mathlib/Tactic/Linarith/SimplexAlgorithm/PositiveVector.lean

@@ -0,0 +1,80 @@
+/-
+Copyright (c) 2024 Vasily Nesterov. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Vasily Nesterov
+-/
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.SimplexAlgorithm
+import Mathlib.Tactic.Linarith.SimplexAlgorithm.Gauss
+
+/-!
+# `linarith` certificate search a LP problem
+
+`linarith` certificate search can easily be reduced to this LP problem: given the matrix `A` and the
+list `strictIndexes`, find the non-negative vector `v` such that some of its coordinates from
+the `strictIndexes` are positive and `A v = 0`.
+
+The function `findPositiveVector` solves this problem.
+-/
+
+namespace Linarith.SimplexAlgorithm
+
+/--
+Given matrix `A` and list `strictIndexes` of strict inequalities' indexes, we want to state the
+Linear Programming problem which solution produces solution for the initial problem (see
+`findPositiveVector`).
+
+As an objective function (that we are trying to maximize) we use sum of coordinates from
+`strictIndexes`: it suffices to find the non-negative vector that makes this function positive.
+
+We introduce two auxiliary variables and one constraint:
+* The variable `y` is interpreted as "homogenized" `1`. We need it because dealing with a
+homogenized problem is easier, but having some "unit" is necessary.
+* To bound the problem we add the constraint `x₁ + ... + xₘ + z = y` introducing new variable `z`.
+
+The objective function also interpreted as an auxiliary variable with constraint
+`f = ∑ i ∈ strictIndexes, xᵢ`.
+
+The variable `f` has to always be basic while `y` has to be free. Our Gauss method implementation
+greedy collects basic variables moving from left to right. So we place `f` before `x`-s and `y`
+after them. We place `z` between `f` and `x` because in this case `z` will be basic and
+`Gauss.getTable` produce table with non-negative last column, meaning that we are starting from
+a feasible point.
+-/
+def stateLP {n m : Nat} (A : Matrix n m) (strictIndexes : List Nat) : Matrix (n + 2) (m + 3) :=
+  Id.run do
+  let mut objectiveRow : Array Rat := #[-1, 0] ++ (Array.mkArray m 0) ++ #[0]
+  for idx in strictIndexes do
+    objectiveRow := objectiveRow.set! (idx + 2) 1 -- +2 due to shifting by `f` and `z`
+
+  let constraintRow : Array Rat := #[0, 1] ++ (Array.mkArray m 1) ++ #[-1]
+
+  let data : Array (Array Rat) := #[objectiveRow, constraintRow]
+    ++ A.data.map (#[0, 0] ++ · ++ #[0])
+
+  return ⟨data⟩
+
+/-- Extracts target vector from the table, putting auxilary variables aside (see `stateLP`). -/
+def exctractSolution (table : Table) : Array Rat := Id.run do
+  let mut ans : Array Rat := Array.mkArray (table.basic.size + table.free.size - 3) 0
+  for i in [1:table.mat.data.size] do
+    ans := ans.set! (table.basic[i]! - 2) table.mat.data[i]!.back
+
+  return ans
+
+/--
+Finds nonnegative vector `v`, such that `A v = 0` and some of its coordinates from `strictCoords`
+are positive, in the case such `v` exists.
+-/
+def findPositiveVector {n m : Nat} (A : Matrix n m) (strictIndexes : List Nat) : Array Rat :=
+  /- State the linear programming problem. -/
+  let B := stateLP A strictIndexes
+
+  /- Using Gaussian elimination split variable into free and basic forming the table that will be
+  operated by Simplex Algorithm.  -/
+  let initTable := Gauss.getTable B
+
+  /- Run Simplex Algorithm and extract the solution. -/
+  let resTable := runSimplexAlgorithm initTable
+  exctractSolution resTable
+
+end Linarith.SimplexAlgorithm