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

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
##]
##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 ---

coordring = base_ring(ideal_ab_s)
vbls = Variable.(symbols(coordring))
##cut_system = CompiledSystem(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 a random witness set
system = CompiledSystem(System(eqns_ab_s, variables = vbls))
sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring))
println("sphere z variables: ", vbls[sph_z_ind])
trivial_soln = fill(0, length(gens(coordring)))
trivial_soln[sph_z_ind] .= 1
println("trivial solutions: $trivial_soln")
norm2 = vec -> real(dot(conj.(vec), vec))
is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring))
##max_slope = 5
##binom = Binomial(2max_slope, 1/2)
Random.seed!(6071)
n_planes = 36
for through_trivial in [false, true]
  samples = []
  for _ in 1:n_planes
    cut_matrix = transpose(hcat(
      (normalize(randn(length(gens(coordring)))) for _ in 1:freedom)...
    ))
    ##cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope
    ##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
    ##]
    ## display(cut_matrix) ## [verbose]
    if through_trivial
      cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
      ## display(cut_offset) ## [verbose]
      cut_subspace = LinearSubspace(cut_matrix, cut_offset)
    else
      cut_subspace = LinearSubspace(cut_matrix, fill(0., freedom))
    end
    wtns = witness_set(system, cut_subspace)
    real_solns = solution.(filter(isreal, results(wtns)))
    nontrivial_solns = filter(is_nontrivial, real_solns)
    println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found")
    for soln in nontrivial_solns
    ## [test] for soln in filter(is_nontrivial, solution.(filter(isreal, results(wtns))))
      if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
        push!(samples, soln)
      end
    end
  end
  if through_trivial
    println("--- planes through trivial solution ---")
  else
    println("--- planes through origin ---")
  end
  println("$(length(samples)) sample solutions, not including the trivial one:")
  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
end