From 86dbd9ea45ce9908c6895355c936a29ffe00629c Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <fenyes@ihes.fr>
Date: Sat, 27 Jan 2024 14:21:03 -0500
Subject: [PATCH] Order variables by coordinate and then element
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

In other words, order coordinates like
  (rₛ₁, rₛ₂, sₛ₁, sₛ₂, xₛ₁, xₛ₂, xₚ₃, yₛ₁, yₛ₂, yₚ₃, zₛ₁, zₛ₂, zₚ₃)
instead of like
  (rₛ₁, sₛ₁, xₛ₁, yₛ₁, zₛ₁, rₛ₂, sₛ₂, xₛ₂, yₛ₂, zₛ₂, xₚ₃, yₚ₃, zₚ₃).

In the test cases, this really cuts down the size of the Gröbner basis.
---
 engine-proto/engine.jl | 76 +++++++++++++++++++++++++++++++-----------
 1 file changed, 57 insertions(+), 19 deletions(-)

diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl
index f672745..7e68fe4 100644
--- a/engine-proto/engine.jl
+++ b/engine-proto/engine.jl
@@ -12,46 +12,68 @@ using Groebner
 abstract type Element{T} end
 
 mutable struct Point{T} <: Element{T}
-  coords::Union{Vector{MPolyRingElem{T}}, Nothing}
+  coords::Vector{MPolyRingElem{T}}
   vec::Union{Vector{MPolyRingElem{T}}, Nothing}
   rel::Nothing
   
   ## [to do] constructor argument never needed?
   Point{T}(
-    coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
+    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
     vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing
   ) where T = new(coords, vec, nothing)
 end
 
-coordnames(_::Point) = [:xₚ, :yₚ, :zₚ]
+##coordnames(_::Point) = [:xₚ, :yₚ, :zₚ]
 
-function buildvec(pt::Point, coordqueue)
-  coordring = parent(coordqueue[1])
-  pt.coords = splice!(coordqueue, 1:3)
+function buildvec!(pt::Point)
+  coordring = parent(pt.coords[1])
   pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
 end
 
 mutable struct Sphere{T} <: Element{T}
-  coords::Union{Vector{MPolyRingElem{T}}, Nothing}
+  coords::Vector{MPolyRingElem{T}}
   vec::Union{Vector{MPolyRingElem{T}}, Nothing}
   rel::Union{MPolyRingElem{T}, Nothing}
   
+  ## [to do] constructor argument never needed?
   Sphere{T}(
-    coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
+    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
     vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
     rel::Union{MPolyRingElem{T}, Nothing} = nothing
   ) where T = new(coords, vec, rel)
 end
 
-coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
+##coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
 
-function buildvec(sph::Sphere, coordqueue)
-  coordring = parent(coordqueue[1])
-  sph.coords = splice!(coordqueue, 1:5)
+function buildvec!(sph::Sphere)
+  coordring = parent(sph.coords[1])
   sph.vec = sph.coords
   sph.rel = mprod(sph.coords, sph.coords) + one(coordring)
 end
 
+const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}(
+  nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ],
+  nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
+)
+
+coordname(elem::Element, index) = coordnames[nameof(typeof(elem))][index]
+
+function pushcoordname!(coordnamelist, indexed_elem::Tuple{Any, Element}, coordindex)
+  elemindex, elem = indexed_elem
+  name = coordname(elem, coordindex)
+  if !isnothing(name)
+    subscript = Subscripts.sub(string(elemindex))
+    push!(coordnamelist, Symbol(name, subscript))
+  end
+end
+
+function takecoord!(coordlist, indexed_elem::Tuple{Any, Element}, coordindex)
+  elem = indexed_elem[2]
+  if !isnothing(coordname(elem, coordindex))
+    push!(elem.coords, popfirst!(coordlist))
+  end
+end
+
 # --- primitive relations ---
 
 abstract type Relation{T} end
@@ -99,22 +121,38 @@ function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
 end
 
 function realize(ctx::Construction{T}) where T
-  # collect variable names
+  # collect coordinate names
   coordnamelist = Symbol[]
   elemenum = enumerate(ctx.elements)
-  for (index, elem) in elemenum
-    subscript = Subscripts.sub(string(index))
-    append!(coordnamelist,
-      [Symbol(name, subscript) for name in coordnames(elem)]
-    )
+  for coordindex in 1:5
+    for indexed_elem in elemenum
+      pushcoordname!(coordnamelist, indexed_elem, coordindex)
+    end
   end
   
+  display(collect(elemenum))
+  display(coordnamelist)
+  println()
+  
   # construct coordinate ring
   coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
   
+  # retrieve coordinates
+  for (_, elem) in elemenum
+    empty!(elem.coords)
+  end
+  for coordindex in 1:5
+    for indexed_elem in elemenum
+      takecoord!(coordqueue, indexed_elem, coordindex)
+    end
+  end
+  
   # construct coordinate vectors
   for (_, elem) in elemenum
-    buildvec(elem, coordqueue)
+    buildvec!(elem)
+    display(elem.coords)
+    display(elem.vec)
+    println()
   end
   
   # turn relations into equations