The vertex-blend of two cubes as an explicit H-polytope

import Mathlib
import Definitions.Def_Hirsch_model

/-!
# The vertex-blend of two cubes as an explicit H-polytope

For `d ≥ 2` the polytope `blendPoly d ⊆ ℝ^d` is cut out by the `3d` integer inequalities
```
  d x_i - s ≤ d,        -d x_i + 3 s ≤ d^2,        -d x_i - s ≤ d^2        (1 ≤ i ≤ d)
```
where `s = ∑ x_i`.  Its graph is two `d`-cubes with one vertex removed each, joined by a
perfect matching of the `d` neighbours of the removed vertices.  It has `2^(d+1) - 2`
vertices, diameter `2d - 1`, and edge expansion `d / (2^d - 1)`, so the graphs of simple
bounded H-polytopes with `O(d)` facets need not have inverse-polynomial expansion in
`n + d`.  (Campaign notes, `deep_astra_geom.md`, §2–4; the graph is the vertex blend of
two cubes, cf. Santos, Annals 2012, Lemma 6.2.)

Rows are indexed by `Fin (3 * d)` through the equivalence with `Fin 3 × Fin d`
(`finProdFinEquiv`): row `(0, i)` is `d e_i - 𝟙`, row `(1, i)` is `-d e_i + 3·𝟙`,
row `(2, i)` is `-d e_i - 𝟙`.
-/

namespace Hirsch

/-- The all-ones vector of `ℝ^d`. -/
noncomputable def ones (d : ℕ) : EuclideanSpace ℝ (Fin d) := WithLp.toLp 2 (fun _ => (1 : ℝ))

/-- The normal of the row with type `t ∈ Fin 3` and coordinate `i ∈ Fin d`. -/
noncomputable def blendRow (d : ℕ) (t : Fin 3) (i : Fin d) : EuclideanSpace ℝ (Fin d) :=
  if t = 0 then (d : ℝ) • EuclideanSpace.single i (1 : ℝ) - ones d
  else if t = 1 then -((d : ℝ) • EuclideanSpace.single i (1 : ℝ)) + (3 : ℝ) • ones d
  else -((d : ℝ) • EuclideanSpace.single i (1 : ℝ)) - ones d

/-- The `3d` normals of the cube-blend polytope, indexed by `Fin (3 * d)`. -/
noncomputable def blendA (d : ℕ) : Fin (3 * d) → EuclideanSpace ℝ (Fin d) :=
  fun j => blendRow d (finProdFinEquiv.symm j).1 (finProdFinEquiv.symm j).2

/-- The right-hand sides: `d` for type-0 rows, `d^2` for the others. -/
noncomputable def blendB (d : ℕ) : Fin (3 * d) → ℝ :=
  fun j => if (finProdFinEquiv.symm j).1 = 0 then (d : ℝ) else (d : ℝ) ^ 2

/-- The cube-blend polytope. -/
noncomputable def blendPoly (d : ℕ) : Set (EuclideanSpace ℝ (Fin d)) := Hpoly (blendA d) (blendB d)

end Hirsch