Invalid region returned when there are obstacles that exceed the bounds.

[enriquefernandez]

Hello,

I'm not sure if this is a bug, or if obstacles are not allowed to exceed the bounds and this is expected behavior.

Basically, IRIS is returning a degenerate region that is invalid when there are obstacles that intersect the environment bounds.

The left figure shows, in green, the region that IRIS returns when using the obstacles in blue, the seed point as the red cross and when the boundary is the large box with black edges.

On the right, the same result when the boundary used is the yellow area instead. As you can see the region returned consists of three vertices where two of them are the same point. Moreover, one of those vertices is outside the boundary.

Just leaving this here in case this is something that should be fixed.
By the way, thanks for making this available, Robin. It works fantastically well and it's really easy to use!

The vertices of the regions are:

IRIS region when using big bounds:
 [[  3.86205232   2.89598565]
 [  4.12493869   5.        ]
 [ 15.           1.83012696]
 [ 15.           5.        ]]
IRIS region when using small bounds:
 [[  3.75000000e+00   2.00000000e+00]
 [  3.75000000e+00   2.00002585e+00]
 [  4.20400665e-12   5.03153075e-13]]

The code I used to generate this problem is the following (Jupyter notebook):

%matplotlib inline
import irispy
import matplotlib.pyplot as plt
import numpy as np
import scipy

def plot_polygon(vertices, ax=None, draw_vertices=False, **kwargs):
        if ax is None:
            ax = plt.gca()
        hull = scipy.spatial.ConvexHull(vertices)
        kwargs.setdefault("facecolor", "none")
        
        if draw_vertices:
            xx, yy = list(zip(*vertices))
            plt.plot(xx, yy, 'o',color='red')
        
        return [ax.add_patch(plt.Polygon(xy=vertices[hull.vertices],**kwargs))]
def plot_bounds(bounds, ax=None, **kwargs):
    (minx, miny), (maxx, maxy) = bounds
    vertices = np.array([[minx,miny], [minx,maxy], [maxx, maxy], [maxx, miny]])
    plot_polygon(vertices, ax, **kwargs)

def plot_scenario():
    # Plot big boundary
    plt.figure(figsize=(14,7))
    for i in range(2):
        plt.subplot(1,2,i+1)
        plot_bounds(big_bounds, alpha=0.3, facecolor='none', edgecolor='black')

        plt.xlim([minx-1, maxx+1])
        plt.ylim([miny-1, maxy+1])

        for obs in obstacles:
            vertices = np.array(list(zip(*obs)))
            plot_polygon(vertices,facecolor='blue', alpha=0.5)

        # Plot seed
        plt.plot(seed[0], seed[1], 'xr')
    
    
# Case with no errors
big_bounds = (-5, -4), (15, 5)
small_bounds = (3.75, 2.0), (8, 3.88)
(minx, miny), (maxx, maxy) = big_bounds
obstacles = [np.array([[  3. ,  10. ,  10. ,   3. ,   3. ],
        [  1.5,   1.5,   2. ,   2. ,   1.5]]),
         np.array([[ 2.,  3.,  4.,  3.,  2.],
        [-4.,  4.,  4., -4., -4.]]),
         np.array([[ 7. ,  7.5,  8. ,  7.5,  7. ],
        [-2. , -2. ,  2.5,  2.5, -2. ]])]
seed = np.array([ 6.25318857,  2.48619646])

plot_scenario()

plt.subplot(1,2,1)
# with big bounds
iris_bounds = irispy.Polyhedron.fromBounds(*big_bounds)
ir = irispy.inflate_region(obstacles, seed, iris_bounds)    
reg_vertices1 = ir.getPolyhedron().getDrawingVertices()
plot_polygon(reg_vertices1, draw_vertices=True, facecolor='green', alpha=0.2)


plt.subplot(1,2,2)
# with small bounds
plot_bounds(small_bounds, alpha=0.3, facecolor='yellow', edgecolor='black')
iris_bounds = irispy.Polyhedron.fromBounds(*small_bounds)
ir = irispy.inflate_region(obstacles, seed, iris_bounds)    
reg_vertices2 = ir.getPolyhedron().getDrawingVertices()
plot_polygon(reg_vertices2, draw_vertices=True, facecolor='green', alpha=0.2)

print("IRIS region when using big bounds:\n", reg_vertices1)
print("IRIS region when using small bounds:\n", reg_vertices2)

[rdeits]

Hey, thanks for the detailed report! (And sorry it took me two months to reply...). This is definitely a bug, and I'll see if I can figure out what's going on.

[rdeits]

Thanks again for making this so easy to reproduce!

It turns out that, as far as I can tell, this is not exactly an IRIS bug, but instead an issue with cddlib. IRIS is actually doing the right thing, but the points returned by getDrawingVertices() (internally computed by cddlib) are just total garbage. To demonstrate, try looking at the ellipsoids that IRIS produces:

plot_bounds(big_bounds, alpha=0.3, facecolor='none', edgecolor='black')
plt.xlim([minx-1, maxx+1])
plt.ylim([miny-1, maxy+1])

for obs in obstacles:
    vertices = np.array(list(zip(*obs)))
    plot_polygon(vertices,facecolor='blue', alpha=0.5)

# Plot seed
plt.plot(seed[0], seed[1], 'xr')
ax = plt.gca()
ir.getEllipsoid().draw(ax)
plot_bounds(small_bounds, alpha=0.3, facecolor='yellow', edgecolor='black')

The ellipsoid (whose drawing points are not computed by cddlib) is right where I would expect to see it, which means that the polyhedron is probably correct. Internally, IRIS only ever uses the half-space representation of the polyhedron, so its vertices are only computed when you try to draw it.

To verify that something wonky is going on, we can pull out just a subset of the half-spaces from the final polyhedron and try to draw them. Let's use all but the last half-space:

p = ir.getPolyhedron()
plot_polygon(irispy.Polyhedron(p.getA()[:-1],
                               p.getB()[:-1]).getDrawingVertices(), 
             draw_vertices=True,
             facecolor="green",
             alpha=0.2)

That works fine. The last constraint is:

>>> ir.getPolyhedron().getA()[-1], ir.getPolyhedron().getB()[-1]
(array([-0., -1.]), -2.0)

which is just the constraint that y >= 2.0. That's very clearly consistent with the polyhedron shown above. But actually including that half space gives us the garbage vertices that you're seeing:

p = ir.getPolyhedron()
plot_polygon(irispy.Polyhedron(p.getA()[:],
                               p.getB()[:]).getDrawingVertices(), 
             draw_vertices=True,
             facecolor="green",
             alpha=0.2)

To see if this was an issue with the way IRIS uses cddlib, I passed the exact same A and B to https://github.com/JuliaPolyhedra/CDDLib.jl in Julia, and I got exactly the same garbage vertices that you're seeing here:

julia> using Polyhedra, CDDLib

julia> A = [
       7.2719127679597916902309862052788957953453063964843750000000000000000000000000000000000000000000000000e-01 -6.8643488180003897625169884122442454099655151367187500000000000000000000000000000000000000000000000000e-01;
       -9.9227812777442014890993959852494299411773681640625000000000000000000000000000000000000000000000000000e-01 1.2403272608667256782233323519903933629393577575683593750000000000000000000000000000000000000000000000e-01;
       4.1480457874537551265698201961862334741226732148788869380950927734375000000000000000000000000000000000e-08 -9.9999999999999911182158029987476766109466552734375000000000000000000000000000000000000000000000000000e-01;
       1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00;
       0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00 1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00;
       -1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00 -0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00;
       -0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00 -1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00;
       ];

julia> b = [
           3.7378473714697459939770851633511483669281005859375 ,
           -3.47297834459160004172417757217772305011749267578125,
           -1.9999989615082682803404168225824832916259765625    ,
           8.                                                  ,
           3.87999999999999989341858963598497211933135986328125,
           -3.75                                                ,
           -2.
       ];

julia> h = SimpleHRepresentation(A, b);

julia> p = polyhedron(h, CDDLibrary());

julia> vrep(p)
V-representation
begin
 2 3 real
 1.0 3.750000823973488 2.0
 3.419486915845482e-14 8.557499217096555e-13 6.838973831690964e-14
end

The same exact A and B inputs, when passed to CDD in "exact" (rational) mode give correct vertices:

julia> p = polyhedron(h, CDDLibrary(:exact));


julia> vrep(p)
V-representation
begin
 5 3 rational
 1//1 160402435284077913733595294199767//40251582856530797108415980109824 4368491638549381//1125899906842624
 1//1 67032365826159611//17875293625971088 2//1
 1//1 23016623785096468//3274978363205447 2//1
 1//1 8//1 4683029438162320//1545713938944383
 1//1 8//1 4368491638549381//1125899906842624
end

julia> Float64.(SimpleVRepresentation(vrep(p)).V)
5×2 Array{Float64,2}:
 3.985    3.88
 3.75     2.0
 7.02802  2.0
 8.0      3.02969
 8.0      3.88

So, ultimately, this just boils down to a lesson I've learned many times over: never trust cddlib in floating-point mode.