engine-proto/Engine.jl

@@ -2,13 +2,12 @@ include("HittingSet.jl")
 
 module Engine
 
-export Construction, mprod, codimension, dimension
+export Construction, mprod
 
 import Subscripts
 using LinearAlgebra
 using AbstractAlgebra
 using Groebner
-using HomotopyContinuation: Variable, Expression, System
 using ..HittingSet
 
 # --- commutative algebra ---
@@ -28,34 +27,6 @@ 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
@@ -218,75 +189,39 @@ 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), eqns)
+  Generic.Ideal(coordring, 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)
+ideal_a_s = Engine.realize(ctx)
-##println("A point on a sphere: ", Engine.dimension(ideal_a_s), " degrees of freedom")
+println("A point on a sphere: ", Engine.dimension(ideal_a_s), " degrees of freeom")
 
 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, eqns_ab_s = Engine.realize(ctx)
+ideal_ab_s = Engine.realize(ctx)
-println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freedom")
+println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freeom")
 
-##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
+spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
-##tangencies = [
+tangencies = [
-##  Engine.AlignsWithBy{CoeffType}(
+  Engine.AlignsWithBy{CoeffType}(
-##    spheres[n],
+    spheres[n],
-##    spheres[mod1(n+1, length(spheres))],
+    spheres[mod1(n+1, length(spheres))],
-##    CoeffType(-1//1)
+    CoeffType(-1//1)
-##  )
+  )
-##  for n in 1:3
+  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])))
 ]
-cut_ideal_ab_s = Generic.Ideal(base_ring(ideal_ab_s), [gens(ideal_ab_s); cut])
+ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
-cut_dim = Engine.dimension(cut_ideal_ab_s)
+ideal_tan_sph = Engine.realize(ctx_tan_sph)
-println("Two points on a sphere, after cut: ", cut_dim, " degrees of freedom")
+println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freeom")
-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, 2])
-  witness = witness_set(System(eqns_ab_s, variables = vbls), cut_subspace)
-  println("witness solutions:")
-  for wtns in solutions(witness)
-    display(wtns)
-  end
-end