1 | initial version |
It's ok, I found the solution using MILP:
def center_inscribed_sphere_polytope(polytope, solver=None):
r"""
Return the center and radius of maximal inscribed sphere
INPUT:
- ``polytope`` -- polytope
OUTPUT:
a 2-tuple (center, radius)
EXAMPLES::
sage: P = polytopes.associahedron(['A',3])
sage: center_inscribed_sphere_polytope(P)
([0.03553390593273766, 0.5355339059327378, 0.03553390593273766],
1.4644660940672622)
sage: center_inscribed_sphere_polytope(P + vector((10,10,10)))
([10.035533905932738, 10.535533905932738, 10.035533905932738],
1.4644660940672622)
::
sage: P = Polyhedron([(0,0), (1,0), (0,10)])
sage: center_inscribed_sphere_polytope(P)
([0.47506218943955486, 0.47506218943955486], 0.47506218943955486)
"""
from math import sqrt
from sage.numerical.mip import MixedIntegerLinearProgram
p = MixedIntegerLinearProgram(solver=solver)
x = p.new_variable()
d = p.new_variable()
for ineq in polytope.inequalities_list():
constant = ineq[0]
coeffs = ineq[1:]
norm = sqrt(sum(a**2 for a in coeffs))
S = p.sum(a/norm*x[i] for (i,a) in enumerate(coeffs))
p.add_constraint(constant/norm + S >= d[0])
p.set_objective(d[0])
radius = p.solve()
soln = p.get_values(x)
center = [soln[k] for k in sorted(x.keys())]
return center, radius