add support for serializing linear programs over non-rational fields (#5246)

Co-authored-by: Antony Della Vecchia <vecchia@math.tu-berlin.de>
Co-authored-by: antonydellavecchia <antonydellavecchia@gmail.com>

Author: benlorenz

Project.toml

@@ -40,7 +40,7 @@ LazyArtifacts = "1.6"
 Markdown = "1.6"
 Nemo = "0.51.1"
 Pkg = "1.6"
-Polymake = "0.13.0"
+Polymake = "0.13.1"
 ProgressMeter = "1.10.2"
 Random = "1.6"
 RandomExtensions = "0.4.3"

src/PolyhedralGeometry/Optimization/linear_program.jl

@@ -30,6 +30,7 @@ function linear_program(
     throw(ArgumentError("convention must be set to :min or :max."))
   end
   ambDim = ambient_dim(P)
+  objective = coefficient_field(P).(objective)
   size(objective, 1) == ambDim || error("objective has wrong dimension.")
   lp = Polymake.polytope.LinearProgram{_scalar_type_to_polymake(T)}(;
     LINEAR_OBJECTIVE=homogenize(objective, k)

src/PolyhedralGeometry/Optimization/mixed_integer_linear_program.jl

@@ -37,6 +37,7 @@ function mixed_integer_linear_program(
     integer_variables = 1:ambDim
   end
   size(objective, 1) == ambDim || error("objective has wrong dimension.")
+  objective = coefficient_field(P).(objective)
   milp = Polymake.polytope.MixedIntegerLinearProgram{_scalar_type_to_polymake(T)}(;
     LINEAR_OBJECTIVE=homogenize(objective, k), INTEGER_VARIABLES=Vector(integer_variables)
   )
@@ -128,8 +129,9 @@ function objective_function(
   milp::MixedIntegerLinearProgram{T}; as::Symbol=:pair
 ) where {T<:scalar_types}
   if as == :pair
-    return Vector{T}(milp.polymake_milp.LINEAR_OBJECTIVE[2:end]),
-    convert(T, milp.polymake_milp.LINEAR_OBJECTIVE[1])
+    cf = coefficient_field(milp)
+    return T[cf(x) for x in milp.polymake_milp.LINEAR_OBJECTIVE[2:end]],
+    cf.(milp.polymake_milp.LINEAR_OBJECTIVE[1])
   elseif as == :function
     (c, k) = objective_function(milp; as=:pair)
     return x -> sum(x .* c) + k

src/PolyhedralGeometry/Polyhedron/constructors.jl

@@ -175,6 +175,7 @@ Get the underlying polymake `Polytope`.
 pm_object(P::Polyhedron) = P.pm_polytope
 
 function ==(P0::Polyhedron{T}, P1::Polyhedron{T}) where {T<:scalar_types}
+  @req coefficient_field(P0) == coefficient_field(P1) "Cannot compare polyhedra over different coefficient fields."
   Polymake.polytope.equal_polyhedra(pm_object(P0), pm_object(P1))
 end
 

src/Serialization/PolyhedralGeometry.jl

@@ -34,6 +34,9 @@ end
 
 type_params(obj::T) where {S <: Union{QQFieldElem, Float64}, T <: PolyhedralObject{S}} = TypeParams(T, coefficient_field(obj))
 
+type_params(obj::T) where {S <: Union{QQFieldElem, Float64}, T <: LinearProgram{S}} = TypeParams(T, coefficient_field(obj))
+type_params(obj::T) where {S <: Union{QQFieldElem, Float64}, T <: MixedIntegerLinearProgram{S}} = TypeParams(T, coefficient_field(obj))
+
 function type_params(obj::T) where {S, T <: PolyhedralObject{S}}
   p_dict = _polyhedral_object_as_dict(obj)
   field = p_dict[:_coeff]
@@ -44,15 +47,19 @@ function type_params(obj::T) where {S, T <: PolyhedralObject{S}}
     :pm_params => type_params(p_dict))
 end
 
+function type_params(obj::T) where {S, T <: Union{LinearProgram{S}, MixedIntegerLinearProgram{S}}}
+  par = params(type_params(feasible_region(obj)))
+  return TypeParams(
+    T,
+    par...)
+end
+
+
 function save_object(s::SerializerState, obj::PolyhedralObject{S}) where S <: Union{QQFieldElem, Float64}
   save_object(s, pm_object(obj))
 end
 
 function save_object(s::SerializerState, obj::PolyhedralObject{<:FieldElem})
-  if typeof(obj) <: Union{MixedIntegerLinearProgram, LinearProgram}
-    T = typeof(obj)
-    error("Unsupported type $T for serialization")
-  end
   p_dict = _polyhedral_object_as_dict(obj)
   delete!(p_dict, :_coeff)
   save_data_dict(s) do
@@ -105,6 +112,15 @@ function save_object(s::SerializerState, lp::LinearProgram{QQFieldElem})
   end
 end
 
+function save_object(s::SerializerState, lp::LinearProgram{<:FieldElem})
+  lpcoeffs = lp.polymake_lp.LINEAR_OBJECTIVE
+  save_data_dict(s) do
+    save_object(s, lp.feasible_region, :feasible_region)
+    save_object(s, lp.convention, :convention)
+    save_object(s, _pmdata_for_oscar(lpcoeffs, coefficient_field(lp)), :lpcoeffs)
+  end
+end
+
 function save_object(s::SerializerState{<: LPSerializer}, lp::LinearProgram{QQFieldElem})
   lp_filename = basepath(s.serializer) * "-$(objectid(lp)).lp"
   save_lp(lp_filename, lp)
@@ -134,6 +150,28 @@ function load_object(s::DeserializerState, ::Type{<:LinearProgram}, field::QQFie
   return LinearProgram{coeff_type}(fr, lp, Symbol(conv))
 end
 
+function load_object(s::DeserializerState, ::Type{<:LinearProgram}, params::Dict)
+  if s.obj isa String
+    error("Loading this file requires using the LPSerializer")
+  end
+  field = params[:field]
+  coeff_type = elem_type(field)
+  fr = load_object(s, Polyhedron, params, :feasible_region)
+  conv = load_object(s, String, :convention)
+  lpcoeffs = load_object(s, Vector{coeff_type}, field, :lpcoeffs)
+  all = Polymake._lookup_multi(pm_object(fr), "LP")
+  lp = nothing
+  for i in 1:length(all)
+    lo = _pmdata_for_oscar(all[i].LINEAR_OBJECTIVE, field)
+    if lpcoeffs == lo
+      lp = all[i]
+      break
+    end
+  end
+  @req lp !== nothing "could not identify LP subobject"
+  return LinearProgram{coeff_type}(fr, lp, Symbol(conv), field)
+end
+
 function load_object(s::DeserializerState{LPSerializer},
                      ::Type{<:LinearProgram}, field::QQField)
   load_node(s) do _
@@ -158,7 +196,19 @@ function save_object(s::SerializerState, milp::MixedIntegerLinearProgram{QQField
     save_object(s, milp.feasible_region, :feasible_region)
     save_object(s, milp.convention, :convention)
     save_json(s, coeffs_jsonstr, :milp_coeffs)
-    save_json(s, int_vars_jsonstr, :int_vars)
+    # we can probably changes this line to just store a vector{Int} but will need upgrade
+    save_json(s, int_vars_jsonstr, :int_vars) 
+  end
+end
+
+function save_object(s::SerializerState, milp::MixedIntegerLinearProgram{<:FieldElem})
+  milp_coeffs = milp.polymake_milp.LINEAR_OBJECTIVE
+  int_vars = milp.polymake_milp.INTEGER_VARIABLES
+  save_data_dict(s) do
+    save_object(s, milp.feasible_region, :feasible_region)
+    save_object(s, milp.convention, :convention)
+    save_object(s, _pmdata_for_oscar(milp_coeffs, coefficient_field(milp)), :milp_coeffs)
+    save_object(s, _pmdata_for_oscar(int_vars, QQ), :int_vars)
   end
 end
 
@@ -193,6 +243,26 @@ function load_object(s::DeserializerState, ::Type{<: MixedIntegerLinearProgram},
   return MixedIntegerLinearProgram{T}(fr, lp, Symbol(conv), field)
 end
 
+function load_object(s::DeserializerState, ::Type{<: MixedIntegerLinearProgram}, params::Dict)
+  conv = load_object(s, String, :convention)
+  field = params[:field]
+  coeff_type = elem_type(field)
+  fr = load_object(s, Polyhedron, params, :feasible_region)
+  milp_coeffs = load_object(s, Vector{coeff_type}, field, :milp_coeffs)
+  int_vars = load_object(s, Vector{Int}, :int_vars)
+  
+  all = Polymake._lookup_multi(pm_object(fr), "MILP")
+  index = 0
+  for i in 1:length(all)
+    if _pmdata_for_oscar(all[i].LINEAR_OBJECTIVE, field) == milp_coeffs && all[i].INTEGER_VARIABLES == Set(int_vars)
+      index = i
+      break
+    end
+  end
+  milp = Polymake._lookup_multi(pm_object(fr), "MILP", index-1)
+  return MixedIntegerLinearProgram{coeff_type}(fr, milp, Symbol(conv), field)
+end
+
 # use generic serialization for the other types:
 @register_serialization_type Cone
 @register_serialization_type PolyhedralComplex

src/Serialization/polymake.jl

@@ -116,14 +116,22 @@ _pmdata_for_oscar(s::Polymake.Set, coeff::Field) = Set(_pmdata_for_oscar(e, coef
 
 function _bigobject_to_dict(bo::Polymake.BigObject, coeff::Field, parent_key::String="")
   data = Dict{Symbol,Any}()
+  bot = Polymake.bigobject_type(bo)
   for pname in Polymake.list_properties(bo)
     p = Polymake.give(bo, pname)
     key_str = parent_key == "" ? pname : parent_key * "." * pname
     if p isa Polymake.PropertyValue
       @debug "missing c++ mapping: skipping $pname of type $(Polymake.typeinfo_string(p, true))"
     elseif p isa Polymake.BigObject
-       obj = _bigobject_to_dict(p, coeff, key_str)
-       merge!(data, obj)
+      if Polymake.bigobject_prop_is_multiple(bot, pname)
+        arr = Polymake._lookup_multi(bo, pname)
+        # this must be without parent key (since will be stored below a parent key)
+        d = Oscar.Serialization._bigobject_to_dict.(arr, Ref(coeff))
+        data[Symbol(pname)] = Tuple(d)
+      else
+        obj = _bigobject_to_dict(p, coeff, key_str)
+        merge!(data, obj)
+      end
     else
       try
         obj = _pmdata_for_oscar(p, coeff)
@@ -155,13 +163,20 @@ end
 
 function _load_bigobject_from_dict!(obj::Polymake.BigObject, dict::Dict, parent_key::String="")
   delay_loading = Tuple{Symbol,Any}[]
+  bot = Polymake.bigobject_type(obj)
   for (k, v) in dict
     # keys of dict are symbols
     k = string(k)
     key_str = parent_key == "" ? k : parent_key * "." * k
     first(k) == '_' && continue
-
-    if v isa Dict
+    if v isa Tuple && Polymake.bigobject_prop_is_multiple(bot, key_str)
+      subobjtype = Polymake.bigobject_prop_type(bot, key_str)
+      for e in v
+        subobj = Polymake.BigObject(subobjtype)
+        Polymake.add(obj, key_str, subobj)
+        _load_bigobject_from_dict!(subobj, e, "") # the parent string is empty because we assign to the subobject
+      end
+    elseif v isa Dict
       _load_bigobject_from_dict!(obj, v, key_str)
     else
       pmv = convert(Polymake.PolymakeType, v)

test/Serialization/PolyhedralGeometry.jl

@@ -137,6 +137,21 @@ using Oscar: _integer_variables
         @test objective_function(LP) == objective_function(loaded)
         @test feasible_region(LP) == feasible_region(loaded)
       end
+
+      P = dodecahedron()
+      F = coefficient_field(P)
+      a = gen(number_field(F))
+      LP1 = linear_program(P, F.([3, -2, 4]);k=2,convention = :min)
+      LP2 = linear_program(P, F.([-1, a - 2, a + 5]);k=2,convention = :min)
+      ov1 = optimal_value(LP1)
+      ov2 = optimal_value(LP2)
+      test_save_load_roundtrip(path, [LP1, LP2]) do (loaded1, loaded2)
+        @test objective_function(LP1) == objective_function(loaded1)
+        @test feasible_region(LP1) == feasible_region(loaded1)
+        @test ov1 == optimal_value(loaded1)
+        @test ov2 == optimal_value(loaded2)
+        @test feasible_region(LP1) == feasible_region(LP2)
+      end
     end
 
     @testset "MixedIntegerLinearProgram" begin
@@ -153,6 +168,34 @@ using Oscar: _integer_variables
         @test feasible_region(MILP) == feasible_region(loaded)
         @test Oscar._integer_variables(MILP) == Oscar._integer_variables(loaded)
       end
+
+      P = dodecahedron()
+      F = coefficient_field(P)
+      a = gen(number_field(F))
+
+      MILP1 = mixed_integer_linear_program(
+        P,
+        [3,-2, a];
+        k=2,
+        convention = :min,
+        integer_variables=[1, 2]
+      )
+
+      MILP2 = mixed_integer_linear_program(
+        P,
+        [-3 * a,-2, 3];
+        k=2,
+        convention = :max,
+        integer_variables=[1, 2]
+      )
+
+      test_save_load_roundtrip(path, [MILP1, MILP2]) do (loaded1, loaded2)
+        @test objective_function(MILP1) == objective_function(loaded1)
+        @test feasible_region(MILP1) == feasible_region(loaded1)
+        @test feasible_region(MILP1) == feasible_region(loaded2)
+        @test Oscar._integer_variables(MILP1) == Oscar._integer_variables(loaded1)
+        @test Oscar._integer_variables(MILP2) == Oscar._integer_variables(loaded2)
+      end
     end
 
     @testset "SubdivisionOfPoints" begin