Implement chamber counting algorithm for toric line bundles

src/AlgebraicGeometry/ToricVarieties/JToric.jl

@@ -57,6 +57,7 @@ include("ToricLineBundles/constructors.jl")
 include("ToricLineBundles/properties.jl")
 include("ToricLineBundles/attributes.jl")
 include("ToricLineBundles/standard_constructions.jl")
+include("ToricLineBundles/cech-cohomology.jl")
 
 include("CohomologyClasses/constructors.jl")
 include("CohomologyClasses/properties.jl")

src/AlgebraicGeometry/ToricVarieties/ToricLineBundles/cech-cohomology.jl

@@ -0,0 +1,128 @@
+#The following function is called, when algorithm="chambers" is specified for all_cohomologies
+function _all_cohomologies_via_cech(tl::ToricLineBundle)
+  our_points, our_maps = _toric_cech_complex(tl)
+  our_maps = transpose.(our_maps)
+  X = toric_variety(tl)
+
+  #If dim(X)<=2, then the intersection of dim(X)+1 maximal cones will be trivial and simplifications can be applied 
+  if dim(X) > 2
+    alternating_sum = abs(sum(sum(length.(values(our_points[k]))) * (-1)^k for k in dim(X)+2:n_maximal_cones(X); init = 0)) 
+  else
+    alternating_sum = abs(sum(binomial(n_maximal_cones(X), k) * (-1)^k for k in dim(X)+2:n_maximal_cones(X); init = 0))
+    alternating_sum *= Int(ncols(our_maps[end]) // binomial(n_maximal_cones(X), dim(X)+1))
+  end
+  @req ncols(our_maps[end]) >= alternating_sum "Inconsistency encountered"
+  rks_maps = push!([rref(M)[1] for M in our_maps], alternating_sum)
+  return ZZRingElem.([nrows(our_maps[1]) - rks_maps[1]; [nrows(our_maps[k]) - rks_maps[k] - rks_maps[k-1] for k in 2:dim(X)]; ncols(our_maps[dim(X)]) - rks_maps[dim(X) + 1] - rks_maps[dim(X)]])
+end
+
+function _toric_cech_complex(tl::ToricLineBundle)
+  # Extract essential information
+  X = toric_variety(tl)
+  @req is_complete(X) "Toric variety must be complete!"
+
+  # Compute the hyperplane arrangement
+  a_plane = matrix(ZZ, n_rays(X), 1, coefficients(toric_divisor(tl)))
+  sc = cone_from_inequalities(matrix(ZZ, [[-1; zeros(Int, dim(X))]]))
+  H = Polymake.fan.HyperplaneArrangement( HYPERPLANES = [a_plane matrix(ZZ, rays(X))], SUPPORT=sc.pm_cone)
+
+  # Classify the regions obtained from the hyperlane arrangement
+  pff = Polymake.fan.PolyhedralComplex(POINTS=H.CHAMBER_DECOMPOSITION.RAYS, INPUT_CONES=H.CHAMBER_DECOMPOSITION.MAXIMAL_CONES)
+  max_polys = maximal_polyhedra(polyhedral_complex(pff))
+  other_polys = filter(is_bounded, vcat([polyhedra_of_dim(polyhedral_complex(pff), i) for i in 0:dim(polyhedral_complex(pff))-1]...))
+  chamber_signs = 2*matrix(ZZ, H.CHAMBER_SIGNATURES).-1
+  sign_of_chamber = Dict(chamber_signs[i,:] => interior_lattice_points(p) for (i, p) in enumerate(max_polys) if is_bounded(p))
+
+  # For the other regions, we don't get their signs for free and have to calculate them
+  for p in other_polys
+    if dim(p) == 0
+      sign_list = ZZ.(sign.(matrix(ZZ, rays(X)) * vertices(p)[1] + a_plane)[:,1])
+    else
+      sign_list = ZZ.(sign.(matrix(ZZ, rays(X))  * sum(vertices(p)) * 1//n_vertices(p) + a_plane)[:,1])
+    end
+    sign_of_chamber[sign_list] = interior_lattice_points(p)
+  end
+
+  # Prepare information, which we use to iterate over the Cech complex and identify the relevant lattice points
+  RI = ray_indices(maximal_cones(X))
+  ray_index_list = map(row -> findall(!iszero, collect(row)), eachrow(RI))
+
+  # Now iterate over the Cech complex
+  cech_complex_points = Vector{Dict{Vector{Int64}, Vector{PointVector{ZZRingElem}}}}(undef, n_maximal_cones(X))
+  cech_complex_maps = Vector{Any}(undef, dim(X))
+  comb_dict = Dict{Combination{Int64}, Dict{PointVector{ZZRingElem}, Int64}}()
+  d_k = 0
+  previous_number_of_generators = 0
+  for k in 0:dim(X)
+
+    # Find the contributing lattice points
+    polyhedron_dict = Dict{Vector{Int64}, Vector{PointVector{ZZRingElem}}}()
+    combs = collect(combinations(n_maximal_cones(X), k+1))
+    for i in 1:length(combs)
+      list_of_lattice_points = PointVector{ZZRingElem}[]
+      generating_ray_indices = reduce(intersect, ray_index_list[combs[i]])
+      for (signs, pts) in sign_of_chamber
+        all(x -> x >= 0, signs[generating_ray_indices]) && append!(list_of_lattice_points, pts)
+      end
+      polyhedron_dict[combs[i]] = list_of_lattice_points
+    end
+
+    # Initialize comb_dict before using it
+    offset = 0
+    for i in 1:length(combs)
+      this_comb = combs[i]
+      pts = polyhedron_dict[this_comb]
+      inner_dict = Dict(pt => j + offset for (j, pt) in enumerate(pts))
+      comb_dict[this_comb] = inner_dict
+      offset += length(pts)
+    end
+
+    # Compute Cech differential maps
+    if k > 0
+      our_sparse_matrix = sparse_matrix(QQ, 0, previous_number_of_generators)
+      col_idx = 1
+      for comb in combs
+        pts = polyhedron_dict[comb]
+        sub_combs = collect(combinations(comb, k))  # returns lex-sorted k-subsets
+        for pt in pts
+          position_vector = Vector{Int}()
+          value_vector = Vector{QQFieldElem}()
+          for j in 1:length(sub_combs)
+            sub_comb = sub_combs[j]
+            if haskey(comb_dict, sub_comb) && haskey(comb_dict[sub_comb], pt)
+              row_idx = comb_dict[sub_comb][pt]
+              sign = (-1)^(j + 1)  # Čech sign convention (1-based indexing)
+              push!(position_vector, row_idx)
+              push!(value_vector, sign)
+            end
+          end
+          push!(our_sparse_matrix, sparse_row(QQ, position_vector, value_vector))
+          col_idx += 1
+        end
+      end
+      cech_complex_maps[k] = our_sparse_matrix
+    end
+    cech_complex_points[k+1] = polyhedron_dict
+    previous_number_of_generators = sum(length.(values(polyhedron_dict)))
+  end
+
+  #If dim(X)<=2, then the intersection of dim(X)+1 maximal cones will be trivial and simplifications can be applied 
+  if dim(X) > 2
+    for k in dim(X)+1:n_maximal_cones(X)-1
+
+      # Find the contributing lattice points
+      polyhedron_dict = Dict{Vector{Int64}, Vector{PointVector{ZZRingElem}}}()
+      combs = collect(combinations(n_maximal_cones(X), k+1))
+      for i in 1:length(combs)
+        list_of_lattice_points = PointVector{ZZRingElem}[]
+        generating_ray_indices = reduce(intersect, ray_index_list[combs[i]])
+        for (signs, pts) in sign_of_chamber
+          all(x -> x >= 0, signs[generating_ray_indices]) && append!(list_of_lattice_points, pts)
+        end
+        polyhedron_dict[combs[i]] = list_of_lattice_points
+      end
+      cech_complex_points[k+1] = polyhedron_dict
+    end
+  end
+  return cech_complex_points, cech_complex_maps
+end

src/AlgebraicGeometry/ToricVarieties/cohomCalg/special_attributes.jl

@@ -71,9 +71,10 @@ end
     all_cohomologies(l::ToricLineBundle)
 
 Compute the dimension of all sheaf cohomologies of the 
-toric line bundle `l` by use of the cohomCalg algorithm 
+toric line bundle `l`. The default algorithm is the cohomCalg algorithm 
 [BJRR10](@cite), [BJRR10*1](@cite) (see also [RR10](@cite),
-[Jow11](@cite) and [BJRR12](@cite)).
+[Jow11](@cite) and [BJRR12](@cite)). It is also possible to specify algorithm = "chamber counting"
+in which case the chamber counting algorithm will be used [CLS11](@cite) p.398 .
 
 # Examples
 ```jldoctest
@@ -85,9 +86,28 @@ julia> all_cohomologies(toric_line_bundle(dP3, [1, 2, 3, 4]))
  0
  16
  0
+
+julia> all_cohomologies(toric_line_bundle(dP3, [1, 2, 3, 4]); algorithm = "chamber counting")
+3-element Vector{ZZRingElem}:
+ 0
+ 16
+ 0
+
+julia> all_cohomologies(toric_line_bundle(dP3, [-3,-2,-2,-2]))
+3-element Vector{ZZRingElem}:
+ 0
+ 2
+ 0
+
+julia> all_cohomologies(toric_line_bundle(dP3, [-3,-2,-2,-2]); algorithm = "chamber counting")
+3-element Vector{ZZRingElem}:
+ 0
+ 2
+ 0
 ```
 """
-@attr Vector{ZZRingElem} function all_cohomologies(l::ToricLineBundle)
+@attr Vector{ZZRingElem} function all_cohomologies(l::ToricLineBundle; algorithm::String = "cohomCalg")
+  if occursin("cohomcalg", lowercase(algorithm))
     # check if we can apply cohomCalg
     v = toric_variety(l)
     if !((is_smooth(v) && is_complete(v)) || (is_simplicial(v) && is_projective(v)))
@@ -190,9 +210,11 @@ julia> all_cohomologies(toric_line_bundle(dP3, [1, 2, 3, 4]))
 
     # return result
     return result
+  elseif occursin("chamber", lowercase(algorithm))
+    return _all_cohomologies_via_cech(l)
+  end
 end
 
-
 @doc raw"""
     cohomology(l::ToricLineBundle, i::Int)
 
@@ -210,13 +232,17 @@ julia> cohomology(toric_line_bundle(dP3, [4, 1, 1, 1]), 0)
 12
 ```
 """
-function cohomology(l::ToricLineBundle, i::Int)
-    v = toric_variety(l)
-    if has_attribute(v, :vanishing_sets)
-        tvs = vanishing_sets(v)[i+1]
-        if contains(tvs, l)
-            return 0
-        end
+function cohomology(l::ToricLineBundle, i::Int; algorithm::String = "cohomCalg")
+  v = toric_variety(l)
+  if has_attribute(v, :vanishing_sets)
+    tvs = vanishing_sets(v)[i+1]
+    if contains(tvs, l)
+      return 0
     end
-    return all_cohomologies(l)[i+1]
+  end
+  if occursin("cohomcalg", lowercase(algorithm))
+    return all_cohomologies(l; algorithm = "cohomCalg")[i+1]
+  elseif occursin("chamber", lowercase(algorithm))
+    return _all_cohomologies_via_cech(l)[i+1]
+  end
 end

src/PolyhedralGeometry/PolyhedralComplex/properties.jl

@@ -370,10 +370,21 @@ _vector_matrix(::Val{_ray_polyhedral_complex}, PC::PolyhedralComplex; homogenize
     _ray_indices_polyhedral_complex(pm_object(PC)), (homogenized ? 1 : 2):end
   ]
 
-_maximal_polyhedron(
+function _maximal_polyhedron(
   ::Type{Polyhedron{T}}, PC::PolyhedralComplex{T}, i::Base.Integer
-) where {T<:scalar_types} =
-  Polyhedron{T}(Polymake.fan.polytope(pm_object(PC), i - 1), coefficient_field(PC))
+) where {T<:scalar_types}
+  pmp = Polymake.fan.polytope(pm_object(PC), i - 1)
+  if Polymake.exists(pmp, "VERTICES") && Polymake.exists(pmp, "VERTICES_IN_FACETS")
+    # workaround wrong column indexing in polymake until next polymake release
+    if nrows(pmp.VERTICES) < ncols(pmp.VERTICES_IN_FACETS)
+      mci = Polymake.row(maximal_polyhedra(IncidenceMatrix, PC), i)
+      pmp = Polymake.polytope.Polytope{_scalar_type_to_polymake(T)}(;
+        VERTICES=pmp.VERTICES, VERTICES_IN_FACETS=pmp.VERTICES_IN_FACETS[:, collect(mci)]
+      )
+    end
+  end
+  Polyhedron{T}(pmp, coefficient_field(PC))
+end
 
 _vertex_indices(::Val{_maximal_polyhedron}, PC::PolyhedralComplex) =
   pm_object(PC).MAXIMAL_POLYTOPES[:, _vertex_indices(pm_object(PC))]

test/AlgebraicGeometry/ToricVarieties/line_bundle_cohomologies.jl

@@ -2,6 +2,8 @@
 
   P2 = projective_space(NormalToricVariety, 2)
   dP3 = del_pezzo_surface(NormalToricVariety, 3)
+  X = dP3 * projective_space(NormalToricVariety, 1)
+  F0 = hirzebruch_surface(NormalToricVariety, 0)
   F5 = hirzebruch_surface(NormalToricVariety, 5)
 
   ray_generators = [[1, 0, 0,-2,-3], [0, 1, 0,-2,-3], [0, 0, 1,-2,-3], [-1,-1,-1,-2,-3], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1], [0, 0, 0,-2,-3]]
@@ -12,7 +14,9 @@
   l2 = toric_line_bundle(divisor_of_character(F5, [1, 2]))
   l3 = toric_line_bundle(P2, [1])
   l4 = anticanonical_bundle(weierstrass_over_p3)
-
+  l5 = toric_line_bundle(F0, [0,-3])
+  l6 = anticanonical_bundle(X)
+
   vs = vanishing_sets(dP3)
   vs2 = vanishing_sets(P2)
   R,_ = polynomial_ring(QQ, 3)
@@ -27,6 +31,11 @@
     @test all_cohomologies(l2) == [1, 0, 0]
   end
 
+  @testset "Cohomology with chamber counting" begin
+    @test all_cohomologies(l5; algorithm = "chamber counting") == [0, 2, 0]
+    @test all_cohomologies(l6; algorithm = "chamber counting") == [21, 0, 0, 0]
+  end
+
   @testset "Toric vanishing sets of dP3" begin
     @test is_projective_space(toric_variety(vs[1])) == false
     @test length(polyhedra(vs[1])) == 1