feat: non-trivial atoms and convexify improvements

.github/workflows/main.yml

@@ -27,7 +27,7 @@ jobs:
         run: |
           set -o pipefail
           curl -sSfL https://github.com/leanprover/elan/releases/download/v1.4.5/elan-x86_64-unknown-linux-gnu.tar.gz | tar xz
-          ./elan-init -y --default-toolchain leanprover/lean4:v4.3.0-rc1
+          ./elan-init -y --default-toolchain leanprover/lean4:v4.4.0-rc1
           echo "$HOME/.elan/bin" >> $GITHUB_PATH
       - name: Install Rust and Cargo
         run: |

.gitignore

@@ -1,15 +1,15 @@
 .DS_Store
 
+/.lake
 /build
-/solver/build
-/solver/test.cbf
-/solver/test.sol
-/rewriter/*.dot 
-/rewriter/*.svg
-/ffi/build
-/ffi/lake-packages 
-/lake-packages
+/solver/*.cbf
+/solver/*.sol
+/egg-convexify/*.dot 
+/egg-convexify/*.svg
 *.olean
 
 /playground
+/evaluation
+/plots
+
 /CvxLean/Tactic/Solver/Mosek/Path.lean

CvxLean.lean

@@ -1 +1,2 @@
 import CvxLean.Command.Solve
+import CvxLean.Tactic.Convexify.Convexify

CvxLean/Command/Solve.lean

@@ -70,7 +70,7 @@ def getReducedProblemAndReduction (prob : Expr)
   trace[Meta.debug] "probOpt: {probOpt.objFun}"
   -- NOTE: We should get the value from this applied to the float sol point.
 
-  let (_, (forwardMap, backwardMap, probReduction)) ← DCP.canonizeGoalFromExpr probSol
+  let (forwardMap, backwardMap, probReduction) ← DCP.canonizeGoalFromExpr probSol
 
   let probReducedSol := match (← inferType probReduction) with
   | Expr.forallE _ r _ _ => r
@@ -123,9 +123,10 @@ unsafe def evalSolve : CommandElab := fun stx =>
       addProblemDeclaration (probName ++ `reduced) probReducedOpt false
 
       -- Call the solver on prob.reduced and get a point in E.
-      let coeffsData ← determineCoeffsFromExpr probReducedExpr
+      let (coeffsData, sections) ← determineCoeffsFromExpr probReducedExpr
       trace[Meta.debug] "coeffsData: {coeffsData}"
-      let solPointResponse ← Meta.conicSolverFromValues probReducedExpr coeffsData
+      let solPointResponse ←
+        Meta.conicSolverFromValues probReducedExpr coeffsData sections
       trace[Meta.debug] "solPointResponse: {solPointResponse}"
 
       match solPointResponse with

CvxLean/Examples/CovarianceEstimation.lean

@@ -58,8 +58,6 @@ reduction reduction₁₂/problem₂ (n : ℕ) (N : ℕ) (α : ℝ) (y : Fin N 
     intro hR
     rw [nonsing_inv_nonsing_inv R hR.isUnit_det]
 
-set_option trace.Meta.debug true
-
 #print problem₂
 
 set_option maxHeartbeats 20000000

CvxLean/Examples/OptimalVehicleSpeed.lean

@@ -1,7 +1,7 @@
 import CvxLean.Command.Equivalence
 import CvxLean.Command.Solve
 import CvxLean.Tactic.Basic.Rename
-import CvxLean.Tactic.PreDCP.ChangeOfVariables
+import CvxLean.Tactic.ChangeOfVariables.ChangeOfVariables
 
 noncomputable section OptimalVehicleSpeed
 
@@ -35,11 +35,11 @@ open FinsetInterval
 
 instance [i : Fact (0 < n)] : OfNat (Fin n) 0 := ⟨⟨0, i.out⟩⟩
 
-def optimalVehicleSpeed (_ : Fact (0 < n)) := 
+def optimalVehicleSpeed (_ : Fact (0 < n)) :=
   optimization (s : Fin n → ℝ)
-    minimize ∑ i, (d i / s i) * Φ (s i) 
-    subject to 
-      hsmin : ∀ i, smin i ≤ s i 
+    minimize ∑ i, (d i / s i) * Φ (s i)
+    subject to
+      hsmin : ∀ i, smin i ≤ s i
       hsmax : ∀ i, s i ≤ smax i
       hτmin : ∀ i, τmin i ≤ ∑ j in [[0, i]], d j / s j
       hτmax : ∀ i, ∑ j in [[0, i]], d j / s j ≤ τmax i
@@ -48,46 +48,46 @@ private lemma simp_fraction : ∀ s : Fin n → ℝ, ∀ i, d i / (d i / s i) =
   intros s i
   have h : d i ≠ 0 := by {
     have d_pos_i := d_pos i
-    linarith 
+    linarith
   }
   rw [← div_mul, div_self h, one_mul]
 }
 
-equivalence equiv/optimalVehicleSpeedConvex 
-  {n : ℕ} 
+equivalence equiv/optimalVehicleSpeedConvex
+  {n : ℕ}
   (n_pos : 0 < n)
   (d : Fin n → ℝ)
   (d_pos : ∀ i, 0 < d i)
   (τmin τmax : Fin n → ℝ)
   (smin smax : Fin n → ℝ)
   (smin_pos : ∀ i, 0 < smin i)
-  (Φ : ℝ → ℝ) : 
-  optimalVehicleSpeed d τmin τmax smin smax Φ ⟨n_pos⟩ := by 
+  (Φ : ℝ → ℝ) :
+  optimalVehicleSpeed d τmin τmax smin smax Φ ⟨n_pos⟩ := by
   unfold optimalVehicleSpeed
   apply Minimization.Equivalence.trans
   apply ChangeOfVariables.toEquivalence (fun t => d / t)
   { rintro s ⟨hsmin, _hsmax, _hτmin, _hτmax⟩ i
     split_ands
     { --suffices h : 0 < s i by exact (ne_of_lt h).symm
-      --exact lt_of_lt_of_le (smin_pos i) (hsmin i) 
+      --exact lt_of_lt_of_le (smin_pos i) (hsmin i)
       have smin_pos_i := smin_pos i
       have hsmin_i := hsmin i
       linarith
       }
-    { --exact (ne_of_lt (d_pos i)).symm 
+    { --exact (ne_of_lt (d_pos i)).symm
       have d_pos_i := d_pos i
       linarith
       } }
-  -- IDEA: This can be done by change of variables by detecting that the 
+  -- IDEA: This can be done by change of variables by detecting that the
   -- variable is a vector.
   apply Minimization.Equivalence.trans
   apply MinimizationQ.rewrite_objective
   { rintro s ⟨_hsmin, _hsmax, _hτmin, _hτmax⟩
     simp only [Pi.div_apply, simp_fraction d d_pos s]
-    rfl 
+    rfl
   }
   apply Minimization.Equivalence.trans
-  apply MinimizationQ.rewrite_constraint_3 
+  apply MinimizationQ.rewrite_constraint_3
   { rintro s _hsmin _hsmax _hτmax
     simp only [Pi.div_apply, simp_fraction d d_pos s]
     rfl
@@ -108,7 +108,7 @@ equivalence equiv/optimalVehicleSpeedConvex
 -- TODO: rename_vars in equivalence mode.
 
 -- IDEA: simplify_objective [...], simplify_constraint i [...],
--- Also rewrite_objective and rewrite_constraint. Potentially get rid of the 
+-- Also rewrite_objective and rewrite_constraint. Potentially get rid of the
 -- rewrite_constraint_x lemmas and handle that in meta.
 
 end OptimalVehicleSpeed

CvxLean/Examples/TrajectoryOptimization.lean

@@ -8,18 +8,18 @@ noncomputable section TrajectoryOptimization
 
 open Matrix
 
-def originalBezier (n d : ℕ) 
-  (K : Matrix (Fin (d + 2)) (Fin n) ℝ) 
+def originalBezier (n d : ℕ)
+  (K : Matrix (Fin (d + 2)) (Fin n) ℝ)
   (V : Matrix (Fin (d + 1)) (Fin n) ℝ)
   (A : Matrix (Fin d) (Fin n) ℝ)
   (k : Fin (d + 2) → ℝ)
   (v : Fin (d + 1) → ℝ)
   (a : Fin d → ℝ) :=
-  optimization (x : Fin n → ℝ) (T : ℝ) 
-    minimize T 
-    subject to 
-      hT : 1 ≤ T 
-      hk : K.mulVec x ≤ k 
+  optimization (x : Fin n → ℝ) (T : ℝ)
+    minimize T
+    subject to
+      hT : 1 ≤ T
+      hk : K.mulVec x ≤ k
       hv : V.mulVec x ≤ T • v
       ha : A.mulVec x ≤ T ^ 2 • a
 
@@ -42,7 +42,7 @@ def convexBezier (n d : ℕ)
 variable {R D E : Type} [Preorder R]
 variable (p : Minimization D R) (q : Minimization E R)
 
-structure Relaxation := 
+structure Relaxation :=
   (r : D → E)
   (r_injective : Function.Injective r)
   (r_feasibility : ∀ x, p.constraints x → q.constraints (r x))
@@ -58,13 +58,12 @@ def relaxationBezier (n d : ℕ)
   Relaxation (originalBezier n d K V A k v a) (convexBezier n d K V A k v a) :=
   { r := fun ⟨x, T⟩ => ⟨x, T, T ^ 2⟩,
     r_injective := fun ⟨x, T⟩ ⟨x', T'⟩ h => by rcases h with ⟨hx, hT, _⟩; rfl,
-    r_feasibility := fun ⟨x, T⟩ ⟨hT, hk, hv, ha⟩ => ⟨hT, hk, hv, ha, le_refl _⟩    
+    r_feasibility := fun ⟨x, T⟩ ⟨hT, hk, hv, ha⟩ => ⟨hT, hk, hv, ha, le_refl _⟩
     r_lower_bound := fun ⟨x, T⟩ ⟨hT, _, _, _⟩ => by {
       simp only [convexBezier, originalBezier]
       rw [sqrt_le_iff]
       have : 0 ≤ T := le_trans zero_le_one hT
       simp [this]
     } }
 
-
-end TrajectoryOptimization
+end TrajectoryOptimization

CvxLean/Lib/Cones/ExpCone.lean

@@ -2,7 +2,6 @@ import Mathlib.Data.Complex.Exponential
 
 namespace Real
 
-@[reducible]
 def expCone (x y z : Real) : Prop :=
   (0 < y ∧ y * exp (x / y) ≤ z) ∨ (y = 0 ∧ 0 ≤ z ∧ x ≤ 0)
 

CvxLean/Lib/Cones/PSDCone.lean

@@ -1,9 +1,9 @@
 import Mathlib.Data.Real.Basic
 import Mathlib.LinearAlgebra.Matrix.PosDef
 
-namespace Real 
+namespace Real
 
-def Matrix.PSDCone {m} [Fintype m] (M : Matrix m m Real) : Prop := 
+def Matrix.PSDCone {m} [Fintype m] (M : Matrix m m Real) : Prop :=
   Matrix.PosSemidef M
 
 end Real

CvxLean/Lib/Cones/SOCone.lean

@@ -20,6 +20,14 @@ def soCone (t : Real) (x : n → Real) : Prop :=
 def rotatedSoCone (v w : Real) (x : n → Real) : Prop :=
   (∑ i, x i ^ 2) ≤ (v * w) * 2 ∧ 0 ≤ v ∧ 0 ≤ w
 
+def Vec.soCone (t : m → Real) (X : Matrix m n Real) : Prop :=
+  ∀ i, Real.soCone (t i) (X i)
+
+def Vec.rotatedSoCone (v w : m → Real) (X : Matrix m n Real) : Prop :=
+  ∀ i, Real.rotatedSoCone (v i) (w i) (X i)
+
+section ConeConversion
+
 noncomputable def rotateSoCone {n : ℕ} (t : Real) (x : Fin n.succ → Real) :
   Real × Real × (Fin n → Real) :=
   ((t + x 0) / sqrt 2, (t - x 0) / sqrt 2, fun i => x i.succ)
@@ -28,29 +36,35 @@ noncomputable def unrotateSoCone {n : ℕ} (v w : Real) (x : Fin n → Real) :
    Real × (Fin n.succ → Real) :=
 ((v + w) / sqrt 2, Matrix.vecCons ((v - w) / sqrt 2) x)
 
-def Vec.soCone (t : m → Real) (X : Matrix m n Real) : Prop :=
-  ∀ i, Real.soCone (t i) (X i)
-
-def Vec.rotatedSoCone (v w : m → Real) (X : Matrix m n Real) : Prop :=
-  ∀ i, Real.rotatedSoCone (v i) (w i) (X i)
+end ConeConversion
 
 section Lemmas
 
 /-- To handle powers, a common trick is to consider for x, y ≥ 0, z ∈ ℝ,
       ((x + y), (x - y, 2 * z)^T) ∈ SO,
 which is equivalent to z ^ 2 ≤ x * y. -/
-lemma soCone_add_sub_two_mul_of_nonneg {x y : ℝ} (z : ℝ) (hx : 0 ≤ x) (hy : 0 ≤ y) :
+lemma soCone_add_sub_two_mul_of_nonneg {x y : ℝ} (z : ℝ)
+  (hx : 0 ≤ x) (hy : 0 ≤ y) :
   soCone (x + y) ![x - y, 2 * z] ↔ z ^ (2 : ℝ) ≤ x * y := by
   have hxy := add_nonneg hx hy
-  conv => lhs; simp [soCone, sqrt_le_left hxy, ←le_sub_iff_add_le']
+  conv => lhs; unfold soCone; simp [sqrt_le_left hxy, ←le_sub_iff_add_le']
   ring_nf; simp
 
+/-- Same as `soCone_add_sub_two_mul_of_nonneg` with z = 1. -/
 lemma soCone_add_sub_two_of_nonneg {x y : ℝ} (hx : 0 ≤ x) (hy : 0 ≤ y) :
   soCone (x + y) ![x - y, 2] ↔ 1 ≤ x * y := by
   have h := soCone_add_sub_two_mul_of_nonneg 1 hx hy
   rw [mul_one, one_rpow] at h
   exact h
 
+/-- Similar trick. -/
+lemma soCone_sub_add_two_mul_of_nonneg {x y : ℝ} (z : ℝ) :
+  soCone (x - y) ![x + y, 2 * z] ↔ y ≤ x ∧ z ^ (2 : ℝ) ≤ -(x * y) := by
+  conv => lhs; unfold soCone; simp [sqrt_le_iff, ←le_sub_iff_add_le']
+  apply Iff.and
+  { rfl }
+  { ring_nf!; rw [←neg_mul, ←div_le_iff (by norm_num)]; simp }
+
 end Lemmas
 
 end Real

CvxLean/Lib/Cones/ZeroCone.lean

@@ -1,11 +1,15 @@
 import Mathlib.Data.Real.Basic
+import CvxLean.Lib.Math.Data.Matrix
 
-namespace Real 
+namespace Real
 
 def zeroCone (x : Real) : Prop :=
   x = 0
 
-def Vec.zeroCone (x : Fin n → Real) : Prop := 
+def Vec.zeroCone (x : Fin n → Real) : Prop :=
   ∀ i, Real.zeroCone (x i)
 
-end Real 
+def Matrix.zeroCone (M : Matrix (Fin n) (Fin m) Real) : Prop :=
+  ∀ i j, Real.zeroCone (M i j)
+
+end Real

CvxLean/Lib/Math/CovarianceEstimation.lean

@@ -18,7 +18,6 @@ noncomputable def covarianceMatrix {N n : ℕ} (y : Fin N → Fin n → ℝ) : M
 lemma gaussianPdf_pos {n : ℕ} (R : Matrix (Fin n) (Fin n) ℝ) (y : Fin n → ℝ) (h : R.PosDef):
   0 < gaussianPdf R y := by
   refine' mul_pos (div_pos zero_lt_one (sqrt_pos.2 (mul_pos _ h.det_pos))) (exp_pos _)
-  simp [rpow_eq_pow]
   exact (pow_pos (mul_pos (by positivity) pi_pos) n)
 
 lemma prod_gaussianPdf_pos {N n : ℕ} (R : Matrix (Fin n) (Fin n) ℝ) (y : Fin N → Fin n → ℝ)
@@ -33,7 +32,6 @@ lemma log_prod_gaussianPdf {N n : ℕ} (y : Fin N → Fin n → ℝ) (R : Matrix
       i ∈ Finset.univ → gaussianPdf R (y i) ≠ 0 := λ i hi => ne_of_gt (gaussianPdf_pos _ _ hR)
     have sqrt_2_pi_n_R_det_ne_zero: sqrt ((2 * π) ^ n * R.det) ≠ 0 := by
       refine' ne_of_gt (sqrt_pos.2 (mul_pos _ hR.det_pos))
-      simp [rpow_eq_pow]
       exact (pow_pos (mul_pos (by positivity) pi_pos) _)
     rw [log_prod Finset.univ (λ i => gaussianPdf R (y i)) this]
     unfold gaussianPdf
@@ -43,15 +41,12 @@ lemma log_prod_gaussianPdf {N n : ℕ} (y : Fin N → Fin n → ℝ) (R : Matrix
     simp [rpow_eq_pow]
     exact ne_of_gt (sqrt_pos.2 (pow_pos (mul_pos (by positivity) pi_pos) _))
     exact ne_of_gt (sqrt_pos.2 hR.det_pos)
-    simp [rpow_eq_pow]
     exact pow_nonneg (mul_nonneg (by positivity) (le_of_lt pi_pos)) _
     norm_num
     exact sqrt_2_pi_n_R_det_ne_zero
     exact div_ne_zero (by norm_num) sqrt_2_pi_n_R_det_ne_zero
     exact exp_ne_zero _
 
-#check congrArg
-
 lemma sum_quadForm {n : ℕ} (R : Matrix (Fin n) (Fin n) ℝ) {m : ℕ} (y : Fin m → Fin n → ℝ) :
   (∑ i, R.quadForm (y i))
   = m * (covarianceMatrix y * Rᵀ).trace := by