Ask Your Question
2

Find the center of a maximum inscribed sphere of a convex polytope

asked 2020-04-23 04:57:31 -0600

Sébastien gravatar image

updated 2020-04-23 06:31:17 -0600

Is there a way in SageMath to find the center of the (or of some, if not unique) sphere of maximum radius which is inside of a convex polytope? I ask the question in general dimension, but I am currently interested in a solution in dimension 2.

edit retag flag offensive close merge delete

1 answer

Sort by » oldest newest most voted
2

answered 2020-04-23 15:52:44 -0600

Sébastien gravatar image

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
edit flag offensive delete link more

Comments

I don't know if the above can be changed to find the center of the John ellipsoid instead? https://en.wikipedia.org/wiki/John_el...

Sébastien gravatar imageSébastien ( 2020-04-29 02:47:25 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2020-04-23 04:57:31 -0600

Seen: 52 times

Last updated: Apr 23