Some issue with the function _calc_resulting_polygons()

We have some issues regarding the function _calc_resulting_polygons, specifically with the for loop when a polygon d is inside a polygon b (and vice versa). Since if b.contains(d) returns True (at Line 914 and 934), then d.intersects(b) should also return True. I was wondering whether this seems to conflict with the logic implemented in your if..else statements, because it would not reach elif b.contains(d) in the next else since d.intersects(b) is True first?

Here is an example to illustrate my point:

from shapely.geometry import Polygon

# Create two squares
d = Polygon([(0.1, 0.1), (1.1, 0.1), (1.1, 1.1), (0.1, 1.1), (0.1, 0.1)])
b = Polygon([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)])

d.intersects(b) # return True

b.contains(d) # return True

This code has been helping us to resolve our issues so far:

def _calc_resulting_polygons(self, boundary_polygons, incl_excls):
    """
    Parameters
    ----------
    boundary_polygons : list
        list of shapely polygons as specifed or inferred from user input
    Returns
    -------
    list of merged shapely polygons. Resolves issues arrising if any are overlapping, touching or contained in each other
    """
    included_polygons = [
        boundary_polygons[i] for i, x in enumerate(incl_excls) if x == 1
    ]
    excluded_polygons = [
        boundary_polygons[i] for i, x in enumerate(incl_excls) if x == 0
    ]

    Inc_poly = unary_union(included_polygons)
    Exc_poly = unary_union(excluded_polygons)
    Remain_polys = Inc_poly.difference(Exc_poly)

    # Filter the polygons based on their area
    if isinstance(Remain_polys, type(Polygon())):
        if Remain_polys.area > 1e-3:
            return [Remain_polys]
        else:
            return []
    elif isinstance(Remain_polys, type(MultiPolygon())):
        filtered_polygons = [
            poly for poly in Remain_polys.geoms if poly.area > 1e-3
        ]
        return filtered_polygons

Thank you for your consideration.