Realize relations as equations

engine-proto/engine.jl

@@ -9,34 +9,47 @@ using Groebner
 
 # --- primitve elements ---
 
-mutable struct Point{T}
+abstract type Element{T} end
+
+mutable struct Point{T} <: Element{T}
   coords::Union{Vector{MPolyRingElem{T}}, Nothing}
   vec::Union{Vector{MPolyRingElem{T}}, Nothing}
+  rel::Nothing
   
   ## [to do] constructor argument never needed?
-  Point{T}(vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing) where T = new(vec)
+  Point{T}(
+    coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
+    vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing
+  ) where T = new(coords, vec, nothing)
 end
 
 coordnames(_::Point) = [:xₚ, :yₚ, :zₚ]
 
 function buildvec(pt::Point, coordqueue)
-  pt.coords = splice!(coordqueue, 1:3)
   coordring = parent(coordqueue[1])
+  pt.coords = splice!(coordqueue, 1:3)
   pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
 end
 
-mutable struct Sphere{T}
+mutable struct Sphere{T} <: Element{T}
   coords::Union{Vector{MPolyRingElem{T}}, Nothing}
   vec::Union{Vector{MPolyRingElem{T}}, Nothing}
+  rel::Union{MPolyRingElem{T}, Nothing}
   
-  Sphere{T}(vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing) where T = new(vec)
+  Sphere{T}(
+    coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
+    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ₛ]
 
 function buildvec(sph::Sphere, coordqueue)
+  coordring = parent(coordqueue[1])
   sph.coords = splice!(coordqueue, 1:5)
   sph.vec = sph.coords
+  sph.rel = mprod(sph.coords, sph.coords) + one(coordring)
 end
 
 # --- primitive relations ---
@@ -45,52 +58,70 @@ abstract type Relation{T} end
 
 mprod(v, w) = v[1]*w[2] + w[1]*v[2] - dot(v[3:end], w[3:end])
 
+# elements: point, sphere
 struct LiesOn{T} <: Relation{T}
-  pt::Point{T}
-  sph::Sphere{T}
+  elements::Vector{Element{T}}
+  
+  LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph])
 end
 
+equation(rel::LiesOn) = dot(rel.elements[1].vec, rel.elements[2].vec)
+
+# elements: sphere, sphere
 struct AlignsWithBy{T} <: Relation{T}
-  sph_v::Sphere{T}
-  sph_w::Sphere{T}
+  elements::Vector{Element{T}}
   cos_angle::T
+  
+  LiesOn{T}(sph1::Point{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle)
 end
 
+equation(rel::AlignsWithBy) = dot(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle
+
 # --- constructions ---
 
 mutable struct Construction{T}
-  points::Vector{Point{T}}
-  spheres::Vector{Sphere{T}}
+  elements::Set{Element{T}}
+  relations::Set{Relation{T}}
   
-  Construction{T}(; points = Point{T}[], spheres = Sphere{T}[]) where T = new{T}(points, spheres)
+  function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T
+    allelements = union(elements, (rel.elements for rel in relations)...)
+    new{T}(allelements, relations)
+  end
 end
 
-function Base.push!(ctx::Construction{T}, elem::Point{T}) where T
-  push!(ctx.points, elem)
+function Base.push!(ctx::Construction{T}, elem::Element{T}) where T
+  push!(ctx.elements, elem)
 end
 
-function Base.push!(ctx::Construction{T}, elem::Sphere{T}) where T
-  push!(ctx.spheres, elem)
+function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
+  push!(ctx.relations, rel)
+  union!(ctx.elements, rel.elements)
 end
 
 function realize(ctx::Construction{T}) where T
   # collect variable names
-  allcoordnames = Symbol[]
-  elements = vcat(ctx.points, ctx.spheres)
-  for (index, elem) in enumerate(elements)
+  coordnamelist = Symbol[]
+  elemenum = enumerate(ctx.elements)
+  for (index, elem) in elemenum
     subscript = Subscripts.sub(string(index))
-    append!(allcoordnames,
+    append!(coordnamelist,
       [Symbol(name, subscript) for name in coordnames(elem)]
     )
   end
   
   # construct coordinate ring
-  coordring, coordqueue = polynomial_ring(parent_type(T)(), allcoordnames)
+  coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
   
   # construct coordinate vectors
-  for elem in elements
+  for (_, elem) in elemenum
     buildvec(elem, coordqueue)
   end
+  
+  # turn relations into equations
+  vcat(
+    equation.(ctx.relations),
+    [elem.rel for elem in ctx.elements if !isnothing(elem.rel)]
+  )
 end
 
 end
@@ -98,8 +129,14 @@ end
 # ~~~ sandbox setup ~~~
 
 a = Engine.Point{Rational{Int64}}()
-b = Engine.Point{Rational{Int64}}()
 s = Engine.Sphere{Rational{Int64}}()
-ctx = Engine.Construction{Rational{Int64}}(points = [a])
+a_on_s = Engine.LiesOn{Rational{Int64}}(a, s)
+ctx = Engine.Construction{Rational{Int64}}(elements = Set([a]), relations= Set([a_on_s]))
+eqns_a_s = Engine.realize(ctx)
+
+b = Engine.Point{Rational{Int64}}()
+b_on_s = Engine.LiesOn{Rational{Int64}}(b, s)
 Engine.push!(ctx, b)
-Engine.push!(ctx, s)
+Engine.push!(ctx, s)
+Engine.push!(ctx, b_on_s)
+eqns_ab_s = Engine.realize(ctx)