include("HittingSet.jl")

module Engine

include("Engine.Algebraic.jl")
include("Engine.Numerical.jl")

using .Algebraic
using .Numerical

export Construction, mprod, codimension, dimension

end

# ~~~ sandbox setup ~~~

using Random
using Distributions
using LinearAlgebra
using AbstractAlgebra
using HomotopyContinuation
using GLMakie

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 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, eqns_ab_s = Engine.realize(ctx)
##freedom = Engine.dimension(ideal_ab_s)
##println("Two points on a sphere: $freedom 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
##]
##tangencies = [
  ##Engine.LiesOn{CoeffType}(points[1], spheres[2]),
  ##Engine.LiesOn{CoeffType}(points[1], spheres[3]),
  ##Engine.LiesOn{CoeffType}(points[2], spheres[3]),
  ##Engine.LiesOn{CoeffType}(points[2], spheres[1]),
  ##Engine.LiesOn{CoeffType}(points[3], spheres[1]),
  ##Engine.LiesOn{CoeffType}(points[3], spheres[2])
##]
##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
##ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
##freedom = Engine.dimension(ideal_tan_sph)
##println("Three mutually tangent spheres: $freedom degrees of freedom")

p = Engine.Point{CoeffType}()
q = Engine.Point{CoeffType}()
a = Engine.Sphere{CoeffType}()
b = Engine.Sphere{CoeffType}()
p_on_a = Engine.LiesOn{CoeffType}(p, a)
p_on_b = Engine.LiesOn{CoeffType}(p, b)
q_on_a = Engine.LiesOn{CoeffType}(q, a)
q_on_b = Engine.LiesOn{CoeffType}(q, b)
ctx_joined = Engine.Construction{CoeffType}(
  elements = Set([p, q, a, b]),
  relations= Set([p_on_a, p_on_b, q_on_a, q_on_b])
)
ideal_joined, eqns_joined = Engine.realize(ctx_joined)
freedom = Engine.dimension(ideal_joined)
println("Two points on two spheres: $freedom degrees of freedom")

# --- test rational cut ---

coordring = base_ring(ideal_joined)
vbls = Variable.(symbols(coordring))

# test a random witness set
system = CompiledSystem(System(eqns_joined, variables = vbls))
norm2 = vec -> real(dot(conj.(vec), vec))
Random.seed!(6071)
n_planes = 3
samples = []
for _ in 1:n_planes
  real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
  for soln in real_solns
    if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
      push!(samples, soln)
    end
  end
end
println("$(length(samples)) sample solutions:")
for soln in samples
  ## display([vbls round.(soln, digits = 6)]) ## [verbose]
  k_sq = abs2(soln[1])
  if abs2(soln[end-2]) > 1e-12
    if k_sq < 1e-12
      println("  center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))")
    else
      sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq
      println("  center on z axis:   r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))")
    end
  else
    sum_sq = sum(soln[[4, 7, 10]] .^ 2)
    println("  center at origin:   r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
  end
end

# show a sample solution
function show_solution(ctx, vals)
  # evaluate elements
  real_vals = real.(vals)
  disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]
  disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres]
  
  # create scene
  scene = Scene()
  cam3d!(scene)
  scatter!(scene, disp_points, color = :green)
  for sph in disp_spheres
    mesh!(scene, sph, color = :gray)
  end
  scene
end