diff --git a/src/TropicalGeometry/TropicalGeometry.jl b/src/TropicalGeometry/TropicalGeometry.jl index bff650ae63fa..c9cdc655bbea 100644 --- a/src/TropicalGeometry/TropicalGeometry.jl +++ b/src/TropicalGeometry/TropicalGeometry.jl @@ -7,15 +7,20 @@ include("semiring.jl") include("semiring_map.jl") include("matrix.jl") include("poly.jl") +include("variety_supertype.jl") include("homogenization.jl") include("groebner_basis.jl") include("initial.jl") # include("groebner_polyhedron.jl") # include("points.jl") -include("variety_supertype.jl") include("hypersurface.jl") include("curve.jl") include("linear_space.jl") include("variety.jl") +include("variety_binomial.jl") +include("variety_principal.jl") +include("variety_affine_linear.jl") +include("variety_zerodimensional.jl") +include("variety_prime.jl") include("intersection.jl") include("groebner_fan.jl") diff --git a/src/TropicalGeometry/homogenization.jl b/src/TropicalGeometry/homogenization.jl index d58d78093319..e60edd24f901 100644 --- a/src/TropicalGeometry/homogenization.jl +++ b/src/TropicalGeometry/homogenization.jl @@ -34,13 +34,14 @@ function homogenize_pre_tropicalization(I::MPolyIdeal) end -function dehomogenize_post_tropicalization(Sigma::PolyhedralComplex) - @req lineality_dim(Sigma)>0 "dehomogenizing polyhedral complex without lineality" +function dehomogenize_post_tropicalization(TropV::TropicalVarietySupertype) + + @req lineality_dim(TropV)>0 "dehomogenizing polyhedral complex without lineality" ### # Construct hyperplane {first coord = 0} ### - n = ambient_dim(Sigma) + n = ambient_dim(TropV) zerothUnitRowVector = zeros(Int,1,n) zerothUnitRowVector[1,1] = 1 dehomogenisingHyperplane = polyhedron((zeros(Int,0,n),zeros(Int,0)), (zerothUnitRowVector,[0])) @@ -52,7 +53,7 @@ function dehomogenize_post_tropicalization(Sigma::PolyhedralComplex) dehomogenizedVertices = Vector{QQFieldElem}[] incidenceMatrixRays = Vector{Int}[] dehomogenizedRays = Vector{QQFieldElem}[] - for sigma in maximal_polyhedra(Sigma) + for sigma in maximal_polyhedra(TropV) sigmaDehomogenized = intersect(sigma,dehomogenisingHyperplane) incidenceVectorVertices = Int[] V,_ = minimal_faces(sigmaDehomogenized) @@ -95,12 +96,17 @@ function dehomogenize_post_tropicalization(Sigma::PolyhedralComplex) ### # Dehomogenize lineality space ### - sigma = first(maximal_polyhedra(Sigma)) + sigma = first(maximal_polyhedra(TropV)) sigmaDehomogenized = intersect(sigma,dehomogenisingHyperplane) dehomogenizedLineality = [linealityVector[2:end] for linealityVector in lineality_space(sigmaDehomogenized)] - return polyhedral_complex(incidenceMatrixVerticesAndRays, - dehomogenizedVerticesAndRays, - collect(length(dehomogenizedVertices)+1:length(dehomogenizedVertices)+length(dehomogenizedRays)), - dehomogenizedLineality) + SigmaDehom = polyhedral_complex(incidenceMatrixVerticesAndRays, + dehomogenizedVerticesAndRays, + collect(length(dehomogenizedVertices)+1:length(dehomogenizedVertices)+length(dehomogenizedRays)), + dehomogenizedLineality) + + TropVDehom = tropical_variety(SigmaDehom,multiplicities(TropV),convention(TropV)) + copy_tropical_attributes!(TropVDehom,TropV) + return TropVDehom + end diff --git a/src/TropicalGeometry/variety.jl b/src/TropicalGeometry/variety.jl index 66f481378f70..a576d6da8718 100644 --- a/src/TropicalGeometry/variety.jl +++ b/src/TropicalGeometry/variety.jl @@ -24,10 +24,10 @@ end # ################################################################################ -function Base.show(io::IO, tv::TropicalVariety{typeof(min), true}) +function Base.show(io::IO, ::TropicalVariety{typeof(min), true}) print(io, "Min tropical variety") end -function Base.show(io::IO, tv::TropicalVariety{typeof(max), true}) +function Base.show(io::IO, ::TropicalVariety{typeof(max), true}) print(io, "Max tropical variety") end @@ -39,8 +39,34 @@ end # ################################################################################ +function copy_tropical_attributes!(TropVcopyTo::TropicalVarietySupertype, TropVcopyFrom::TropicalVarietySupertype) + inheritableAttributes = + [ + # from all + :tropical_semiring_map, + # from tropical hypersurfaces + :tropical_polynomial, + :algebraic_polynomial, + :dual_subdivision, + # from tropical linear spaces + :polymake_object, + :pluecker_indices, + :tropical_pluecker_vector, + :tropical_matrix, + :algebraic_pluecker_vector, + :algebraic_matrix, + :algebraic_ideal + # from tropical varieties + ] + for attr in inheritableAttributes + if has_attribute(TropVcopyFrom,attr) + set_attribute!(TropVcopyTo,attr,get_attribute(TropVcopyFrom,attr)) + end + end +end + @doc raw""" - tropical_variety(Sigma::PolyhedralComplex, mult, minOrMax::Union{typeof(min),typeof(max)}=min) + tropical_variety(Sigma::PolyhedralComplex, mult::Vector{ZZRingElem}, minOrMax::Union{typeof(min),typeof(max)}=min) Return the `TropicalVariety` whose polyhedral complex is `Sigma` with multiplicities `mult` and convention `minOrMax`. Here, `mult` is optional can be specified as a `Vector{ZZRingElem}` which represents a list of multiplicities on the maximal polyhedra in the order of `maximal_polyhedra(Sigma)`. If `mult` is unspecified, then all multiplicities are set to one. @@ -80,28 +106,13 @@ end function tropical_variety(TropH::TropicalHypersurface{minOrMax,true}) where {minOrMax<:Union{typeof(min),typeof(max)}} TropV = tropical_variety(polyhedral_complex(TropH),multiplicities(TropH), convention(TropH)) - - if has_attribute(TropH,:algebraic_polynomial) - set_attribute!(TropV,:algebraic_ideal,ideal([get_attribute(TropH,:algebraic_polynomial)])) - end - if has_attribute(TropH,:tropical_semiring_map) - set_attribute!(TropV,:tropical_semiring_map,get_attribute(TropH,:tropical_semiring_map)) - end - + copy_tropical_attributes!(TropV,TropH) return TropV end - function tropical_variety(TropL::TropicalLinearSpace{minOrMax,true}) where {minOrMax<:Union{typeof(min),typeof(max)}} TropV = tropical_variety(polyhedral_complex(TropL),multiplicities(TropL), convention(TropL)) - - if has_attribute(TropL,:algebraic_ideal) - set_attribute!(TropV,:algebraic_ideal,get_attribute(TropL,:algebraic_ideal)) - end - if has_attribute(TropL,:tropical_semiring_map) - set_attribute!(TropV,:tropical_semiring_map,get_attribute(TropL,:tropical_semiring_map)) - end - + copy_tropical_attributes!(TropV,TropL) return TropV end @@ -134,544 +145,71 @@ end ################################################################################ @doc raw""" - tropical_variety(I::MPolyIdeal, nu::Union{TropicalSemiringMap,Nothing}=nothing; weighted_polyhedral_complex_only::Bool=false, skip_saturation::Bool=false, skip_primary_decomposition::Bool=false) + tropical_variety(I::MPolyIdeal, nu::Union{TropicalSemiringMap,Nothing}=nothing; weighted_polyhedral_complex_only::Bool=false) -Return the tropicalization of `I` with respect to `nu` as a `Vector{TropicalVariety}`. +Return the tropicalization of `I` with respect to `nu`. If `nu==nothing`, will compute with respect to the trivial valuation and min convention. If `weighted_polyhedral_complex_only==true`, will not cache any additional information. -If `skip_saturation==true`, will not saturate `I` with respect to the product of all variables. -If `skip_primary_decomposition==true`, will not decompose `I`. !!! warning - `tropical_variety` is currently under development and only works for ideals that primary decompose into - principal, linear, and binomial ideals. + Experimental feature, only special cases supported: + - any coefficient field and any valuation: `I` principal, binomial, or affine linear + - QQ and trivial / p-adic valuation only: `I` prime # Examples ```jldoctest -julia> R,(x,y) = QQ[:x, :y]; - -julia> I = ideal([(x^2+y)*(x+y^2)*(x+y)]); - -julia> tropical_variety(I) -3-element Vector{TropicalVariety}: - Min tropical variety - Min tropical variety - Min tropical variety - -``` -""" -function tropical_variety(I::Union{MPolyIdeal,MPolyRingElem}, nu::Union{TropicalSemiringMap,Nothing}=nothing; weighted_polyhedral_complex_only::Bool=false, skip_saturation::Bool=false, skip_primary_decomposition::Bool=false) - ### - # Step 0.a: convert I to ideal if poly, - # initialize nu as the trivial valuation if not specified by user - ### - if I isa MPolyRingElem - I = ideal(parent(I),[I]) - end - if isnothing(nu) - nu = tropical_semiring_map(coefficient_ring(I)) - end - - ### - # Step 0.b: Saturate `I` and assert that `I` is not the whole ring - ### - if !skip_saturation - R = base_ring(I) - I = saturation(I,ideal([prod(gens(R))])) - end - @req !isone(I) "ideal contains a monomial, tropical varieties in OSCAR cannot be empty" - - ### - # Step 0.c: Compute primary decomposition and tropical varieties of all primary factors - ### - toTropicalize = [I] - if !skip_primary_decomposition - toTropicalize = [ P for (P,_) in primary_decomposition(I) ] - end - - tropicalVarieties = TropicalVariety[] - for P in toTropicalize - # compute a reduced GB to test whether `P` is principal, binomial, or linear - GB = groebner_basis(P,complete_reduction=true) - - if length(GB)==1 - # P is principal - TropV = tropical_variety_principal(P,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) - elseif max(length.(GB)...)==2 - # P is binomial - TropV = tropical_variety_binomial(P,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) - elseif max(total_degree.(GB)...)==1 - # P linear - TropV = tropical_variety_linear(P,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) - else - # P general - error("general tropical varieties currently unsupported") - end - push!(tropicalVarieties,TropV) - end - - return tropicalVarieties -end - - -################################################################################ -# -# Principal ideals -# -################################################################################ - -function tropical_variety_principal(I::MPolyIdeal,nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false) - ### - # Construct TropicalVariety from TropicalHypersurface - ### - g = first(gens(I)) - TropV = tropical_variety(tropical_hypersurface(g,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)) - if !weighted_polyhedral_complex_only - set_attribute!(TropV,:algebraic_ideal,I) - set_attribute!(TropV,:tropical_semiring_map,nu) - end - return TropV -end - +julia> R,(x,y,z) = QQ["x","y","z"]; +julia> nu_2 = tropical_semiring_map(QQ,2) +Map into Min tropical semiring encoding the 2-adic valuation on Rational field -################################################################################ -# -# Binomial ideals -# -################################################################################ - -function tropical_variety_binomial(I::MPolyIdeal,nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false) - ### - # Construct matrix of exponent vector differences - # and vector of coefficient valuation differences - ### - G = gens(I) - A = matrix(ZZ,[first(collect(expv))-last(collect(expv)) for expv in exponents.(G)]) - b = [ QQ(nu(last(collect(coeff))/first(collect(coeff)))) for coeff in coefficients.(G)] - - ### - # Compute tropical variety multiplicity - ### - snfAdiag = elementary_divisors(A) - weight = abs(prod([m for m in snfAdiag if !iszero(m)])) - - ### - # Constructing tropical variety set-theoretically - ### - A = QQMatrix(A) - L = transpose(kernel(A, side = :right)) - can_solve, V = can_solve_with_solution(transpose(A),matrix(QQ,[b]),side=:left) - @req can_solve "tropical variety cannot be empty" - SigmaV = polyhedral_complex(IncidenceMatrix([[1]]), V, nothing, L) - - ### - # Assemble tropical variety - ### - TropV = tropical_variety(SigmaV,[weight],convention(nu)) - if !weighted_polyhedral_complex_only - set_attribute!(TropV,:algebraic_ideal,I) - set_attribute!(TropV,:tropical_semiring_map,nu) - end - return TropV -end +julia> f1 = 8*x^2 + x*y + x*z + x + 8*y^2 + y*z + y + 8*z^2 + z + 8; +julia> f2 = x + 1; +julia> I = ideal([f1,f2]); -################################################################################ -# -# Linear ideals -# -################################################################################ - -function tropical_variety_linear(I::MPolyIdeal,nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false) - ### - # Compute reduced Groebner basis (usually already cached), - # and check whether the linear polynomials have a constant term - ### - R = base_ring(I) - G = groebner_basis(I,complete_reduction=true) - if min(total_degree.(Iterators.flatten(collect.(terms.(G))))...)==1 - # input homogneeous, construct TropicalVariety via TropicalLinearSpace - TropV = tropical_variety(tropical_linear_space(I,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)) - if !weighted_polyhedral_complex_only - set_attribute!(TropV,:algebraic_ideal,I) - set_attribute!(TropV,:tropical_semiring_map,nu) - end - return TropV - else - # input inhomogeneous, homogenise first - Ih,_,_ = homogenize_pre_tropicalization(I) - TropLh = tropical_linear_space(Ih,nu,weighted_polyhedral_complex_only=true) - Sigma = dehomogenize_post_tropicalization(polyhedral_complex(TropLh)) - - multiplicities = ones(ZZRingElem, n_maximal_polyhedra(Sigma)) - TropV = tropical_variety(Sigma,multiplicities) - if !weighted_polyhedral_complex_only - set_attribute!(TropV,:algebraic_ideal,I) - set_attribute!(TropV,:tropical_semiring_map,nu) - end - return TropV - end -end - - - -################################################################################ -# -# Zero-dimensional ideals -# -################################################################################ - -function tropical_variety_zerodimensional(I::MPolyIdeal,nu::TropicalSemiringMap{QQField,QQFieldElem,<:Union{typeof(min),typeof(max)}}) - - k,(a,_) = number_field(I) - zk = maximal_order(k) - p = uniformizer(nu) - lp = [x[1] for x = prime_decomposition(zk,p)] - ma = representation_matrix(a) - mb = representation_matrix(k(lp[1].gen_two*lp[2].gen_two^2)) - @assert iszero(ma*mb - mb*ma) - Qp = PadicField(p, 10) - TropVDict = simultaneous_diagonalization([map_entries(Qp, ma),map_entries(Qp, mb)]) - - TropVPoints = collect(values(TropVDict)) - TropVPointsUnique = unique(TropVPointsMults) - Sigma = polyhedral_complex(IncidenceMatrix([[i] for i in 1:length(TropVPointsUnique)]), TropVPointsUnique) - TropVMults = [ZZ(length(findall(isequal(p),TropVPoints))) for p in TropVPointsUnique] - TropV = tropical_variety(Sigma,TropVMults) - set_attribute!(TropV,:algebraic_points,collect(keys(TropVDict))) - return TropV -end +julia> TropI_0 = tropical_variety(I) +Min tropical variety +julia> vertices(TropI_0) +1-element SubObjectIterator{PointVector{QQFieldElem}}: + [0, 0, 0] -### -# Code by Claus: -### -function slope_eigenspace(M::MatElem{T}) where T <: Hecke.NonArchLocalFieldElem - f = charpoly(M) - lf = Hecke.slope_factorization(f) - # @req all(==(1), values(lf)) +julia> TropI_2 = tropical_variety(I,nu_2) +Min tropical variety - se = Dict{typeof(f), typeof(M)}() - k = base_ring(M) - zk = maximal_order(k) +julia> vertices(TropI_2) +2-element SubObjectIterator{PointVector{QQFieldElem}}: + [0, -3, 3] + [0, 3, -3] - for f = keys(lf) - se[f] = kernel(f(M), side = :right) #hopefully, this is in rref +``` +""" +function tropical_variety(I::MPolyIdeal, nu::Union{TropicalSemiringMap,Nothing}=nothing; weighted_polyhedral_complex_only::Bool=false) + # initialize nu as trivial valuation if not given + if isnothing(nu) + nu = tropical_semiring_map(coefficient_ring(I)) end - @assert sum(ncols(x) for x = values(se)) == nrows(M) - return se -end -function _intersect(M::MatElem{T}, N::MatElem{T}) where T <: Hecke.FieldElem - k = base_ring(M) - I = [M N] - PR = maximum(precision, I) - pr = minimum(precision, I) - if pr != PR - for i = eachindex(I) - I[i] = setprecision(I[i], pr) - end + # compute a reduced GB to test whether `I` is principal, binomial, or affine linear + GB = groebner_basis(I,complete_reduction=true) + if length(GB)==1 + # I is principal + return tropical_variety_principal(I,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) + elseif all(g->(length(g)==2), GB) + # I is binomial + return tropical_variety_binomial(I,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) + elseif all(g->(total_degree(g)<2), GB) + # I affine linear + return tropical_variety_affine_linear(I,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) end - v = kernel(I, side = :right) #precision issues... - l = M*v[1:ncols(M), 1:ncols(v)] - return transpose(rref(transpose(l))[2]) -end - -function valuation_of_roots(f::PolyRingElem{<:Hecke.NonArchLocalFieldElem}) - iszero(f) && error("polynomial must not be zero") - return (valuation(constant_coefficient(f)) - valuation(leading_coefficient(f)))//degree(f) + # I general + return tropical_variety_prime(I,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) end -function simultaneous_diagonalization(v::Vector{<:MatElem{T}}) where T <: Hecke.NonArchLocalFieldElem - - k = base_ring(v[1]) - @assert all(x->base_ring(x) == k, v) - n = nrows(v[1]) - @assert all(x->ncols(x) == nrows(x) == n, v) - - vv = map(slope_eigenspace, v) - - d = Dict(v => [valuation_of_roots(k)] for (k,v) = vv[1]) - @assert sum(ncols(x) for x = keys(d)) == n - for i=2:length(vv) - dd = typeof(d)() - for (mat, pol_vec) = d - for (p, m) = vv[i] - j = _intersect(mat, m) - if ncols(j) > 0 - dd[j] = push!(copy(pol_vec), valuation_of_roots(p)) - end - end - end - d = dd - @assert sum(ncols(x) for x = keys(d)) == n - end - return d +function tropical_variety(f::MPolyRingElem, nu::Union{TropicalSemiringMap,Nothing}=nothing; weighted_polyhedral_complex_only::Bool=false) + return tropical_variety(ideal(parent(f),[f]),nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) end - - -# #======= -# tropical variety of an ideal -# todo: proper documentation -# Example: -# import Random -# K,s = RationalFunctionField(QQ,"s"); -# Kx,(x1,x2,x3,x4) = polynomial_ring(K,4); -# val = TropicalSemiringMap(K,s); -# I = ideal([x1-s*x2+(s+1)*x3,3*x2-s^2*x3+(s^2+1)*x4]); -# Random.seed!(3847598273423); -# TropI = tropical_variety(I,val) -# =======# -# function tropical_variety(I::MPolyIdeal, val::TropicalSemiringMap, convention::Union{typeof(min),typeof(max)}=min) - -# ### -# # Part 0: Preprocessing -# # Check whether valuation is on the coefficient ring of input polynomials, -# # homogenize input ideal if not homogeneous -# ### -# if coefficient_ring(base_ring(I))!=val.valued_field -# error("input valuation not on coefficient ring of input ideal") -# end -# was_input_homogeneous = true -# for g in groebner_basis(I,complete_reduction=true) # todo: replace GB computation with interreduction -# if !sloppy_is_homogeneous(g) -# was_input_homogeneous = false -# I = homogenize(I) -# break -# end -# end - - - -# ### -# # Part 1: Tropical starting polyhedra (todo: avoid recomputation) -# # - compute and recompute starting points until they lie in the relative interior of maximal cells -# # - initialize working lists -# ### - -# # Note: a working list entry consists of a triple (w,C,G) where -# # * G is a (tropical) Groebner basis -# # * C is the Groebner polyhedron -# # * w is the sum of vertices and rays of C -# # In particular: -# # * w is a weight vector with respect to which G is a Groebner basis, -# # * w is compatible with coordinate permutations if symmetries exist, -# # * instead of comparing C or G it suffices to compare w. -# working_list_todo = [] # list of groebner polyhedra with potentially unknown neighbors -# working_list_done = [] # list of groebner polyhedra with known neighbors -# facet_points_done = [] # list of facet points whose tropical links were computed and traversed - -# compute_starting_points = true -# while compute_starting_points -# print("computing random starting points... ") -# starting_points = tropical_points(I,val) -# println("done") - -# working_list_todo = [] -# working_list_done = [] -# facet_points_done = [] -# compute_starting_points = false -# for starting_point in starting_points -# print("computing groebner_basis for starting point ",starting_point,"... ") -# G = groebner_basis(I,val,starting_point) -# println("done") -# C = groebner_polyhedron(G,val,starting_point) -# w = anchor_point(C) - -# # if C is lower-dimensional, recompute all starting points -# if dim(C)!=dim(I) -# println("starting point on lower-dimensional cell, recomputing...") -# compute_starting_points = true -# break -# end - -# # if (w,C,G) already is in working_list_todo, skip -# i = searchsortedfirst(working_list_todo,(w,C,G),by=first) -# if i<=length(working_list_todo) && working_list_todo[i][1]==w -# continue -# end -# # otherwise, add (w,C,G) to todo list -# insert!(working_list_todo, i, (w,C,G)) -# end -# end - - - -# ### -# # Part 2: Tropical traversal -# ### -# while !isempty(working_list_todo) -# print("#working_list_todo: ",length(working_list_todo)," ") -# println("#working_list_done: ",length(working_list_done)) - -# # pick a groebner polyhedron from todo list, add it to the done list, and compute its facet points -# (w,C,G) = popfirst!(working_list_todo) -# i = searchsortedfirst(working_list_done,(w,C,G),by=first) -# insert!(working_list_done, i, (w,C,G)) - -# points_to_traverse = facet_points(C) -# for point_to_traverse in points_to_traverse -# # if point was traversed before, skip -# i = searchsortedfirst(facet_points_done,point_to_traverse) -# if i<=length(facet_points_done) && facet_points_done[i]==point_to_traverse -# continue -# end -# # otherwise add point_to_traverse to facet_points_done -# insert!(facet_points_done, i, point_to_traverse) - -# directions_to_traverse = tropical_link(ideal(G),val,point_to_traverse) # todo, this output can be wrong -# for direction_to_traverse in directions_to_traverse -# # compute neighbor -# print("computing groebner_basis for ",point_to_traverse,direction_to_traverse,"... ") -# G_neighbor = groebner_flip(G,val,w,point_to_traverse,direction_to_traverse) -# println("done") -# C_neighbor = groebner_polyhedron(G_neighbor,val,point_to_traverse,perturbation=direction_to_traverse) -# w_neighbor = anchor_point(C_neighbor) - -# # if neighbor is already in done list, skip -# i = searchsortedfirst(working_list_done, -# (w_neighbor,C_neighbor,G_neighbor), -# by=first) -# if i<=length(working_list_done) && working_list_done[i][1]==w_neighbor -# continue -# end -# # if neighbor is already in todo list, skip -# i = searchsortedfirst(working_list_todo, -# (w_neighbor,C_neighbor,G_neighbor), -# by=first) -# if i<=length(working_list_todo) && working_list_todo[i][1]==w_neighbor -# continue -# end -# # otherwise, add neighbor to todo list -# insert!(working_list_todo, i, (w_neighbor,C_neighbor,G_neighbor)) -# end -# end -# end - - - -# ### -# # Part 3: Postprocessing -# ### -# # 3.0: dehomogenize data if input was homogenized -# if !was_input_homogeneous -# n = length(gens(base_ring(I))) - -# # 3.0.1: dehomogenize Groebner polyhedra -# zeroth_unit_vector_as_row_vector = zeros(Int,1,n) -# zeroth_unit_vector_as_row_vector[1,1] = 1 -# dehomogenising_hyperplane = polyhedron((zeros(Int,0,n),zeros(Int,0)), -# (zeroth_unit_vector_as_row_vector,[1])) -# for wCG in working_list_done -# wCG[2] = intersect(wCG[2],dehomogenising_hyperplane) -# end -# # 3.0.2: check that initial ideals are distinct (todo) -# end - -# # 3.1: construct PolyhedralComplex -# # 3.1.1: construct incidence_matrix, vertices_and_rays, and far_vertices -# incidence_matrix = Vector{Vector{Int}}() -# verts_rays = Vector{Vector{Polymake.Rational}}() -# far_vertices = Vector{Int}() -# for (w,C,G) in working_list_done -# incidence_vector = Vector{Int}() -# for vert in vertices(C) -# i = findfirst(isequal(vert),verts_rays) -# if i === nothing -# # if vert does not occur in verts_rays -# # add it to verts_rays -# push!(verts_rays,vert) -# push!(incidence_vector,length(verts_rays)) -# else -# push!(incidence_vector,i) -# end -# end -# for ray in rays(C) -# i = findfirst(isequal(ray),verts_rays) -# if i === nothing || !(i in far_vertices) -# # if ray does not occur in verts_rays or if it occurs but not as a ray, -# # add it to verts_rays -# push!(verts_rays,ray) -# push!(far_vertices,length(verts_rays)) -# push!(incidence_vector,length(verts_rays)) -# else -# push!(incidence_vector,i) -# end -# end -# push!(incidence_matrix,incidence_vector) - -# end -# verts_rays_matrix = permutedims(reduce(hcat, verts_rays)) # convert Vector{Vector} to Matrix - -# # 3.1.2: construct lineality space -# (w,C,G) = first(working_list_done) -# lineality_space_gens = matrix(QQ,lineality_space(C)) - - - -# # 3.2: Construct lists for weight_vectors, initial_ideals and multiplicities -# weight_vectors = [w for (w,C,G) in working_list_done] -# initial_ideals = [ideal(initial(G,val,w)) for (w,C,G) in working_list_done] -# multiplicities = [multiplicity(inI) for inI in initial_ideals] - - - -# PC = polyhedral_complex(IncidenceMatrix(incidence_matrix), -# verts_rays_matrix, -# far_vertices, -# lineality_space_gens) - -# verts_rays_perm = collect(vertices_and_rays(PC)) -# verts_rays_perm = Vector{Int64}.(verts_rays_perm) -# permutation = [findfirst(isequal(l), verts_rays_perm) for l in verts_rays] -# # inc = [findall(incidence_matrix[i, :]) for i in 1:size(incidence_matrix, 1)] -# new_incidence = [permutation[incidence] for incidence in incidence_matrix] -# mults = Dict(new_incidence[i] => multiplicities[i] for i in 1:length(working_list_done)) - -# TropI = TropicalVariety{typeof(max),true}(PC, mults) - - -# set_attribute!(TropI,:weight_vectors,weight_vectors) -# set_attribute!(TropI,:initial_ideals,initial_ideals) - -# return TropI -# end - - - -# #======= -# Example: -# P = cube(4) -# anchor_point(P) -# facet_points(P) -# =======# -# function anchor_point(P::Polyhedron) -# # compute the sum of vertices and rays in homogenized coordinates -# pt = convert(Vector{QQFieldElem},sum([vertices(P)...,rays(P)...])) -# pushfirst!(pt,n_vertices(P)) - -# # project to orthogonal complement of lineality space if necessary -# if lineality_dim(P)>0 -# pt = Polymake.Matrix{Polymake.Rational}(vcat(transpose(pt))) -# Polymake.common.project_to_orthogonal_complement(pt, P.pm_polytope.LINEALITY_SPACE) -# pt = convert(Matrix{QQFieldElem}, pt)[1,:] -# end - -# # rescale until first entry is 1 and remove it -# pt = [pt[i]//pt[1] for i in 2:length(pt)] -# return pt -# end - -# function facet_points(P::Polyhedron) -# points = [] -# for facet in faces(P,dim(P)-1) -# if length(vertices(facet))>0 # skipping facets at infinity -# push!(points,anchor_point(facet)) -# end -# end -# return points -# end diff --git a/src/TropicalGeometry/variety_affine_linear.jl b/src/TropicalGeometry/variety_affine_linear.jl new file mode 100644 index 000000000000..add0a3f48ea4 --- /dev/null +++ b/src/TropicalGeometry/variety_affine_linear.jl @@ -0,0 +1,22 @@ +################################################################################ +# +# Tropicalization of affine linear ideals +# +# WARNING: assumes without test that `gens(I)` has degree one +# +################################################################################ + +function tropical_variety_affine_linear(I::MPolyIdeal,nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false) + # compute a reduced GB to check whether I is homogeneous + G = groebner_basis(I,complete_reduction=true) + + if all(Oscar._is_homogeneous,G) + # input homogneeous, construct TropicalVariety via TropicalLinearSpace + return tropical_variety(tropical_linear_space(I,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)) + end + + # input inhomogeneous, homogenise first + Ih,_,_ = homogenize_pre_tropicalization(I) + TropLh = tropical_linear_space(Ih,nu,weighted_polyhedral_complex_only=true) + return dehomogenize_post_tropicalization(TropLh) +end diff --git a/src/TropicalGeometry/variety_binomial.jl b/src/TropicalGeometry/variety_binomial.jl new file mode 100644 index 000000000000..81b2dc997ec7 --- /dev/null +++ b/src/TropicalGeometry/variety_binomial.jl @@ -0,0 +1,41 @@ +################################################################################ +# +# Tropicalization of binomial ideals +# +# WARNING: assumes without test that `gens(I)` is binomial +# +################################################################################ + +function tropical_variety_binomial(I::MPolyIdeal,nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false) + ### + # Construct matrix of exponent vector differences + # and vector of coefficient valuation differences + ### + G = gens(I) + A = matrix(ZZ,[first(expv)-last(expv) for expv in collect.(exponents.(G))]) + b = [ QQ(nu(last(coeff)/first(coeff))) for coeff in collect.(coefficients.(G))] + + ### + # Compute tropical variety multiplicity + ### + snfAdiag = elementary_divisors(A) + weight = abs(prod([m for m in snfAdiag if !iszero(m)])) + + ### + # Constructing tropical variety set-theoretically + ### + A = matrix(QQ, A) + can_solve, V, L = can_solve_with_solution_and_kernel(transpose(A), matrix(QQ,[b]); side=:left) + @req can_solve "tropical variety cannot be empty" + SigmaV = polyhedral_complex(IncidenceMatrix([[1]]), V, nothing, L) + + ### + # Assemble tropical variety + ### + TropV = tropical_variety(SigmaV,[weight],convention(nu)) + if !weighted_polyhedral_complex_only + set_attribute!(TropV,:algebraic_ideal,I) + set_attribute!(TropV,:tropical_semiring_map,nu) + end + return TropV +end diff --git a/src/TropicalGeometry/variety_prime.jl b/src/TropicalGeometry/variety_prime.jl new file mode 100644 index 000000000000..8ac648c2f100 --- /dev/null +++ b/src/TropicalGeometry/variety_prime.jl @@ -0,0 +1,135 @@ +################################################################################ +# +# Tropicalization of prime ideals +# +# WARNING: assumes without test that `I` is prime +# +################################################################################ + +function tropical_variety_prime(I::MPolyIdeal, nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false) + # compute a reduced GB to check whether I is homogeneous + G = groebner_basis(I,complete_reduction=true) + + if all(Oscar._is_homogeneous, G) + return tropical_variety_prime_singular(I,nu; weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) + end + + Ih,_,_ = homogenize_pre_tropicalization(I) + TropIh = tropical_variety_prime_singular(Ih,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only) + return dehomogenize_post_tropicalization(TropIh) +end + +# trivial valuation +function tropical_variety_prime_singular(I::MPolyIdeal, nu::TropicalSemiringMap{QQField,Nothing,<:Union{typeof(min),typeof(max)}}; weighted_polyhedral_complex_only::Bool=false) + R = base_ring(I) + singularCommand = join(["ring r=0,("*join(string.(symbols(R)),",")*"),dp;", + "ideal I = "*join(string.(gens(I)), ",")*";", + "if (!defined(tropicalVariety)) { LIB \"tropical.lib\"; };", + "fan TropI = tropicalVariety(I);", + "string TropIString = string(TropI);"]) + Singular.call_interpreter(singularCommand) + TropIString = Singular.lookup_library_symbol("Top", "TropIString") + Sigma = gfan_fan_string_to_oscar_complex(TropIString,convention(nu)==max,false) + TropI = compute_weights_and_construct_tropical_variety(Sigma,I,nu) + if !weighted_polyhedral_complex_only + set_attribute!(TropI,:algebraic_ideal,I) + set_attribute!(TropI,:tropical_semiring_map,nu) + end + return TropI +end + +# p-adic valuation +function tropical_variety_prime_singular(I::MPolyIdeal, nu::TropicalSemiringMap{QQField,ZZRingElem,<:Union{typeof(min),typeof(max)}}; weighted_polyhedral_complex_only::Bool=false) + R = base_ring(I) + singularCommand = join(["ring r=0,("*join(string.(symbols(R)),",")*"),dp;", + "ideal I = "*join(string.(gens(I)), ",")*";", + "if (!defined(tropicalVariety)) { LIB \"tropical.lib\"; };", + "fan TropI = tropicalVariety(I,number("*string(uniformizer(nu))*"));", + "string TropIString = string(TropI);"]) + Singular.call_interpreter(singularCommand) + TropIString = Singular.lookup_library_symbol("Top", "TropIString") + Sigma = gfan_fan_string_to_oscar_complex(TropIString,convention(nu)==max,true) + TropI = compute_weights_and_construct_tropical_variety(Sigma,I,nu) + if !weighted_polyhedral_complex_only + set_attribute!(TropI,:algebraic_ideal,I) + set_attribute!(TropI,:tropical_semiring_map,nu) + end + return TropI +end + + +# Takes the string of a gfan fan (= Singular fan) and converts it to an OSCAR polyhedral complex +# if negateFan==true, return the negative of the fan (because Singular is max-convention only, and OSCAR is flexible) +# if dehomogenizeFan==true, the gfan fan is a homogenized polyhedral complex, and we need to dehomogenize it +function gfan_fan_string_to_oscar_complex(input_string::String, negateFan::Bool=false, dehomogenizeFan::Bool=false) + + # Extracting AMBIENT_DIM + ambientDim = parse(Int, match(r"AMBIENT_DIM\n(\d+)", input_string).captures[1]) + + # Extracting the RAYS, LINEALITY_SPACE and MAXIMAL_CONES sections + stringsParsed = Vector{SubString{String}}[] + for regexp in [r"RAYS\n([\s\S]*?)\nN_RAYS", r"LINEALITY_SPACE\n([\s\S]*?)\nORTH_LINEALITY_SPACE", r"MAXIMAL_CONES\n([\s\S]*)"] + sectionOfInterest = match(regexp, input_string).captures[1] + linesOfInterest = split(sectionOfInterest, "\n") + linesOfInterestFiltered = [split(line, r"\s+#")[1] for line in linesOfInterest if !isempty(line)] + push!(stringsParsed, linesOfInterestFiltered) + end + + # Convert Rays and ORTH_LINEALITY_SPACE to matrices + # and negate if necessary + rayGenerators = [parse.(Int, split(line)) for line in stringsParsed[1]] + if isempty(rayGenerators) + rayGenerators = zero_matrix(QQ,0,ambientDim) + else + rayGenerators = matrix(QQ,rayGenerators) + end + if negateFan + rayGenerators *= -1 + end + + linealityGenerators = [parse.(Int, split(line)) for line in stringsParsed[2]] + if isempty(linealityGenerators) + linealityGenerators = zero_matrix(QQ,0,ambientDim) + else + linealityGenerators = matrix(QQ,linealityGenerators) + end + + # Convert MAXIMAL_CONES to a Vector{Vector{Int}} + if length(stringsParsed[3])==1 && first(stringsParsed[3])=="{}" + coneIncidences = [Int[]] + else + coneIncidences = [parse.(Int, split(replace(line, r"[{}]" => ""), r"\s+")) .+ 1 for line in stringsParsed[3]] + end + + if dehomogenizeFan + # if the singular fan is a homogenized polyhedral complex, + # identify which rays have a non-zero first entry and actually represent vertices + # then strip the first coordinate of the rays and lineality generators + rayIndices = findall(iszero, rayGenerators[:,1]) + rayGenerators = rayGenerators[:,2:end] + linealityGenerators = linealityGenerators[:,2:end] + # in some cases, the first unit vector can be a lineality generator + # hence we need to filter out potential zero rows from linealityGenerators + linealityGenerators = linealityGenerators[findall(i->!iszero(linealityGenerators[i,:]),1:nrows(linealityGenerators)),:] + return polyhedral_complex(IncidenceMatrix(coneIncidences), rayGenerators, rayIndices, linealityGenerators) + else + # if the singular fan is an honest polyhedral fan + # we need to add the origin as a vertex and add it to all polyhedra + rayIndices = collect(1:nrows(rayGenerators)) + rayGenerators = vcat(rayGenerators, zero_matrix(QQ,1,ncols(rayGenerators))) + originIndex = nrows(rayGenerators) + coneIncidences = [ vcat(incidence,[originIndex]) for incidence in coneIncidences ] + return polyhedral_complex(IncidenceMatrix(coneIncidences), rayGenerators, rayIndices, linealityGenerators) + end + +end + + +# this function takes the tropical variety as a polyhedral complex, +# computes the weights of the tropical variety, and constructs the tropical variety +function compute_weights_and_construct_tropical_variety(Sigma::PolyhedralComplex,::MPolyIdeal,nu::TropicalSemiringMap) + # TODO: compute multiplicities + mults = ones(ZZRingElem,n_maximal_polyhedra(Sigma)) + TropI = TropicalVariety{typeof(convention(nu)),true}(Sigma,mults) + return TropI +end diff --git a/src/TropicalGeometry/variety_principal.jl b/src/TropicalGeometry/variety_principal.jl new file mode 100644 index 000000000000..451aec10ceb6 --- /dev/null +++ b/src/TropicalGeometry/variety_principal.jl @@ -0,0 +1,20 @@ +################################################################################ +# +# Tropicalization of principal ideals +# +# WARNING: assumes without test that `gens(I)` has cardinality one +# +################################################################################ + +function tropical_variety_principal(I::MPolyIdeal,nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false) + ### + # Construct TropicalVariety from TropicalHypersurface + ### + g = gen(I,1) + TropV = tropical_variety(tropical_hypersurface(g,nu,weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)) + if !weighted_polyhedral_complex_only + set_attribute!(TropV,:algebraic_ideal,I) + set_attribute!(TropV,:tropical_semiring_map,nu) + end + return TropV +end diff --git a/src/TropicalGeometry/variety_zerodimensional.jl b/src/TropicalGeometry/variety_zerodimensional.jl new file mode 100644 index 000000000000..f1c0ae3a5d75 --- /dev/null +++ b/src/TropicalGeometry/variety_zerodimensional.jl @@ -0,0 +1,98 @@ +################################################################################ +# +# Tropicalization of zero-dimensional ideals +# +# WARNING: assumes without test that `I` is zero-dimensional +# +# WARNING: currently p-adic valuation only +# +################################################################################ + +function tropical_variety_zerodimensional_eigenvalue(I::MPolyIdeal,nu::TropicalSemiringMap{QQField,QQFieldElem,<:Union{typeof(min),typeof(max)}}) + k,(a,_) = number_field(I) + zk = maximal_order(k) + p = uniformizer(nu) + lp = [x[1] for x = prime_decomposition(zk,p)] + ma = representation_matrix(a) + mb = representation_matrix(k(lp[1].gen_two*lp[2].gen_two^2)) + @assert iszero(ma*mb - mb*ma) + Qp = padic_field(p; precision=10) + TropVDict = simultaneous_diagonalization([map_entries(Qp, ma),map_entries(Qp, mb)]) + + TropVPoints = collect(values(TropVDict)) + TropVPointsUnique = unique(TropVPointsMults) + Sigma = polyhedral_complex(IncidenceMatrix([[i] for i in 1:length(TropVPointsUnique)]), TropVPointsUnique) + TropVMults = [ZZ(length(findall(isequal(p),TropVPoints))) for p in TropVPointsUnique] + TropV = tropical_variety(Sigma,TropVMults) + set_attribute!(TropV,:algebraic_points,collect(keys(TropVDict))) + return TropV +end + + +### +# Code by Claus: +### +function slope_eigenspace(M::MatElem{T}) where T <: Hecke.NonArchLocalFieldElem + f = charpoly(M) + lf = Hecke.slope_factorization(f) + # @req all(x->x==1, values(lf)) + + se = Dict{typeof(f), typeof(M)}() + k = base_ring(M) + zk = maximal_order(k) + + for f = keys(lf) + se[f] = kernel(f(M), side = :right) #hopefully, this is in rref + end + @assert sum(ncols(x) for x = values(se)) == nrows(M) + return se +end + +function _intersect(M::MatElem{T}, N::MatElem{T}) where T <: Hecke.FieldElem + k = base_ring(M) + I = [M N] + PR = maximum(precision, I) + pr = minimum(precision, I) + if pr != PR + for i = eachindex(I) + I[i] = setprecision(I[i], pr) + end + end + + v = kernel(I, side = :right) #precision issues... + l = M*v[1:ncols(M), 1:ncols(v)] + return transpose(rref(transpose(l))[2]) +end + +function valuation_of_roots(f::PolyRingElem{<:Hecke.NonArchLocalFieldElem}) + iszero(f) && error("polynomial must not be zero") + return (valuation(constant_coefficient(f)) - valuation(leading_coefficient(f)))//degree(f) +end + +function simultaneous_diagonalization(v::Vector{<:MatElem{T}}) where T <: Hecke.NonArchLocalFieldElem + + k = base_ring(v[1]) + @assert all(x->base_ring(x) == k, v) + n = nrows(v[1]) + @assert all(x->ncols(x) == nrows(x) == n, v) + + vv = map(slope_eigenspace, v) + + d = Dict(v => [valuation_of_roots(k)] for (k,v) = vv[1]) + @assert sum(ncols(x) for x = keys(d)) == n + for i=2:length(vv) + dd = typeof(d)() + for (mat, pol_vec) = d + for (p, m) = vv[i] + j = _intersect(mat, m) + if ncols(j) > 0 + dd[j] = push!(copy(pol_vec), valuation_of_roots(p)) + end + end + end + d = dd + @assert sum(ncols(x) for x = keys(d)) == n + end + + return d +end diff --git a/test/TropicalGeometry/variety.jl b/test/TropicalGeometry/variety.jl index 3b816574b357..b9be860a41cf 100644 --- a/test/TropicalGeometry/variety.jl +++ b/test/TropicalGeometry/variety.jl @@ -19,11 +19,11 @@ R,(x,y,z,w) = QQ[:x, :y, :z, :w] I = ideal(R,[x+2*y,z+4*w]) nu = tropical_semiring_map(QQ,2) - TropV = first(tropical_variety(I,nu)) + TropV = tropical_variety(I,nu) TropL = tropical_linear_space(I,nu) @test issetequal(maximal_polyhedra(TropV),maximal_polyhedra(TropL)) nu = tropical_semiring_map(QQ,2,max) - TropV = first(tropical_variety(I,nu)) + TropV = tropical_variety(I,nu) TropL = tropical_linear_space(I,nu) @test issetequal(maximal_polyhedra(TropV),maximal_polyhedra(TropL)) end @@ -33,12 +33,36 @@ f = x*y*z+2 nu = tropical_semiring_map(QQ,2) TropH = tropical_hypersurface(f,nu) - TropV = first(tropical_variety(ideal(R,f),nu)) + TropV = tropical_variety(ideal(R,f),nu) @test issetequal(maximal_polyhedra(TropH),maximal_polyhedra(TropV)) nu = tropical_semiring_map(QQ,2,max) - TropV = first(tropical_variety(ideal(R,f),nu)) + TropV = tropical_variety(ideal(R,f),nu) TropH = tropical_hypersurface(f,nu) @test issetequal(maximal_polyhedra(TropH),maximal_polyhedra(TropV)) end + # running tropical_variety and all its subroutines + @testset "testing tropical_variety" begin + + # principal ideals + R,(x,y,z) = QQ["x","y","z"] + f = x^2+y^2+z^2+1 + TropV = tropical_variety(ideal(R,f)) + @test f_vector(TropV) == [1,4,6] + + # binomial ideals + f = x^2+1 + g = y^2+1 + TropV = tropical_variety(ideal(R,[f,g])) + @test f_vector(TropV) == [0,1] + + # affine linear ideals + f = x+z+1 + g = y+z+1 + TropV = tropical_variety(ideal(R,[f,g])) + @test f_vector(TropV) == [1,3] + + # general ideals, see doctests + + end end