engine-proto/Engine.jl

@@ -2,12 +2,13 @@ include("HittingSet.jl")
 
 module Engine
 
-export Construction, mprod
+export Construction, mprod, codimension, dimension
 
 import Subscripts
 using LinearAlgebra
 using AbstractAlgebra
 using Groebner
+using HomotopyContinuation: Variable, Expression, System
 using ..HittingSet
 
 # --- commutative algebra ---
@@ -27,6 +28,34 @@ end
 dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} =
   length(gens(base_ring(I))) - codimension(I, maxdepth)
 
+# hat tip Sascha Timme
+#   https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl/issues/520#issuecomment-1317681521
+function Base.convert(::Type{Expression}, f::MPolyRingElem)
+  variables = Variable.(symbols(parent(f)))
+  f_data = zip(coefficients(f), exponent_vectors(f))
+  sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data)
+end
+
+# create a ModelKit.System from an ideal in a multivariate polynomial ring. the
+# variable ordering is taken from the polynomial ring
+function System(I::Generic.Ideal)
+  eqns = Expression.(gens(I))
+  variables = Variable.(symbols(base_ring(I)))
+  System(eqns, variables = variables)
+end
+
+## [to do] not needed right now
+# create a ModelKit.System from a list of elements of a multivariate polynomial
+# ring. the variable ordering is taken from the polynomial ring
+##function System(eqns::AbstractVector{MPolyRingElem})
+##  if isempty(eqns)
+##    return System([])
+##  else
+##    variables = Variable.(symbols(parent(f)))
+##    return System(Expression.(eqns), variables = variables)
+##  end
+##end
+
 # --- primitve elements ---
 
 abstract type Element{T} end
@@ -189,39 +218,75 @@ function realize(ctx::Construction{T}) where T
     append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4])
   end
   
-  Generic.Ideal(coordring, eqns)
+  (Generic.Ideal(coordring, eqns), eqns)
 end
 
 end
 
 # ~~~ sandbox setup ~~~
 
+using AbstractAlgebra
+using HomotopyContinuation
+
 CoeffType = Rational{Int64}
 
 a = Engine.Point{CoeffType}()
 s = Engine.Sphere{CoeffType}()
 a_on_s = Engine.LiesOn{CoeffType}(a, s)
 ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
-ideal_a_s = Engine.realize(ctx)
-println("A point on a sphere: ", Engine.dimension(ideal_a_s), " degrees of freeom")
+##ideal_a_s = Engine.realize(ctx)
+##println("A point on a sphere: ", Engine.dimension(ideal_a_s), " degrees of freedom")
 
 b = Engine.Point{CoeffType}()
 b_on_s = Engine.LiesOn{CoeffType}(b, s)
 Engine.push!(ctx, b)
 Engine.push!(ctx, s)
 Engine.push!(ctx, b_on_s)
-ideal_ab_s = Engine.realize(ctx)
-println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freeom")
+ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
+println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freedom")
 
-spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
-tangencies = [
-  Engine.AlignsWithBy{CoeffType}(
-    spheres[n],
-    spheres[mod1(n+1, length(spheres))],
-    CoeffType(-1//1)
-  )
-  for n in 1:3
+##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
+##tangencies = [
+##  Engine.AlignsWithBy{CoeffType}(
+##    spheres[n],
+##    spheres[mod1(n+1, length(spheres))],
+##    CoeffType(-1//1)
+##  )
+##  for n in 1:3
+##]
+##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
+##ideal_tan_sph = Engine.realize(ctx_tan_sph)
+##println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freedom")
+
+# --- test rational cut ---
+
+cut = [
+  sum(vcat(a.coords, (s.coords - [0, 0, 0, 0, 1])))
+  sum(vcat([2, 1, 1] .* a.coords, [1, 2, 1, 1, 1] .* s.coords - [0, 0, 0, 0, 1]))
+  sum(vcat([1, 2, 0] .* a.coords, [1, 1, 0, 1, 2] .* s.coords - [0, 0, 0, 0, 1]))
 ]
-ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
-ideal_tan_sph = Engine.realize(ctx_tan_sph)
-println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freeom")
+cut_ideal_ab_s = Generic.Ideal(base_ring(ideal_ab_s), [gens(ideal_ab_s); cut])
+cut_dim = Engine.dimension(cut_ideal_ab_s)
+println("Two points on a sphere, after cut: ", cut_dim, " degrees of freedom")
+if cut_dim == 0
+  vbls = Variable.(symbols(base_ring(ideal_ab_s)))
+  cut_system = System([eqns_ab_s; cut], variables = vbls)
+  cut_result = HomotopyContinuation.solve(cut_system)
+  println("non-singular solutions:")
+  for soln in solutions(cut_result)
+    display(soln)
+  end
+  println("singular solutions:")
+  for sing in singular(cut_result)
+    display(sing.solution)
+  end
+  
+  # test corresponding witness set
+  cut_matrix = [1 1 1 1 0 1 1 0 1 1 0; 1 2 1 2 0 1 1 0 1 1 0; 1 1 0 1 0 1 2 0 2 0 0]
+  cut_subspace = LinearSubspace(cut_matrix, [1, 1, 1])
+  witness = witness_set(System(eqns_ab_s, variables = vbls), cut_subspace)
+  println("witness solutions:")
+  for wtns in solutions(witness)
+    display(wtns)
+  end
+end