Clarification on Behavior of intersect and & Operations in the Polytope Package

[mvnayagam]

Dear Developers and Users,

I have noticed different behaviors between the intersect and & operations when applied to a pc.Region and a pc.Polytope. Let’s say R is a pc.Region contains 11 polytopes, and P is a pc.Polytope.

When I call R & P, it yields some incorrect solutions. However, when I use the expression [Ri.intersect(P) for Ri in R], I obtain the solutions I expect.

Could you clarify the difference between R & P and [Ri.intersect(P) for Ri in R] in terms of how they are implemented in the Polytope package? Understanding this distinction would greatly enhance my comprehension of the package's code structure and implementation style.

Any suggestions and clarifications would be greatly appreciated.

Thank you!

---

Best regards
Muthu Vallinayagam
Research, Institute of Experimental Physics
TU Freiberg
Germany

[mvnayagam]

adding to the above abstract question, I enclose the best example to test & and intersect functions.
I have two regions R1 and R2. In R1 there are two polytopes.
these are

# polytope no 0 in R1
a0=np.array([[-0.57735 -0.57735 -0.57735]
 [ 0.57735  0.57735  0.57735]
 [-1.      -0.      -0.     ]
 [ 0.57735 -0.57735 -0.57735]
 [-0.      -0.      -1.     ]
 [-0.      -0.70711  0.70711]])
b0=np.array([-0.11719  0.11845 -0.2      0.11682  0.       0.     ])
r1p0=pc.Polytope(a0, b0)

# polytope no: 1 in R1
a1=np.array([[ 0.57735 -0.57735 -0.57735]
 [-0.57735  0.57735  0.57735]
 [ 1.       0.       0.     ]
 [-0.57735 -0.57735 -0.57735]
 [-0.      -0.      -1.     ]
 [-0.      -0.70711  0.70711]])
b1=np.array([ 0.11375 -0.11249  0.2     -0.11687  0.       0.     ])
r1p1=pc.Polytope(a1, b1)

---

In R2, there are 5 polytopes

# polytope no: 0 in R2
a0=np.array([
    [-0.57735 -0.57735 -0.57735]
    [ 0.57735  0.57735  0.57735]
    [-1.      -0.      -0.     ]
    [-0.      -1.      -0.     ]
    [-0.      -0.      -1.     ]
    [ 1.       0.       0.     ]
    [ 0.       1.       0.     ]
    [ 0.       0.       1.     ]])
b0=np.array([-0.11693  0.12912 -0.14286 -0.      -0.       0.21429  0.07143  0.07143])
r2p0=pc.Polytope(a0, b0)

# polytope no: 1 in R2
a1=np.array([
    [ 0.57735 -0.57735 -0.57735]
    [-0.57735  0.57735  0.57735]
    [-1.      -0.      -0.     ]
    [-0.      -1.      -0.     ]
    [-0.      -0.      -1.     ]
    [ 1.       0.       0.     ]
    [ 0.       1.       0.     ]
    [ 0.       0.       1.     ]])
b1=np.array([ 0.04803 -0.03584 -0.07143 -0.      -0.       0.14286  0.07143  0.07143])
r2p1=pc.Polytope(a1, b1)

# polytope no: 2 in R2
a2=np.array([
    [ 0.57735 -0.57735 -0.57735]
    [-0.57735  0.57735  0.57735]
    [-1.      -0.      -0.     ]
    [-0.      -1.      -0.     ]
    [-0.      -0.      -1.     ]
    [ 1.       0.       0.     ]
    [ 0.       1.       0.     ]
    [ 0.       0.       1.     ]])
b2=np.array([ 0.13051 -0.11832 -0.21429 -0.      -0.       0.28571  0.07143  0.07143])
r2p2=pc.Polytope(a2, b2)

# polytope no: 3 in R2
a3=np.array([
    [ 0.57735  0.57735 -0.57735]
    [-0.57735 -0.57735  0.57735]
    [-1.      -0.      -0.     ]
    [-0.      -1.      -0.     ]
    [-0.      -0.      -1.     ]
    [ 1.       0.       0.     ]
    [ 0.       1.       0.     ]
    [ 0.       0.       1.     ]])
b3=np.array([ 0.13051 -0.11832 -0.07143 -0.07143 -0.       0.14286  0.14286  0.07143])
r2p3=pc.Polytope(a3, b3)

# polytope no: 4 in R2
a4=np.array([
    [-0.57735  0.57735 -0.57735]
    [ 0.57735 -0.57735  0.57735]
    [-1.      -0.      -0.     ]
    [-0.      -1.      -0.     ]
    [-0.      -0.      -1.     ]
    [ 1.       0.       0.     ]
    [ 0.       1.       0.     ]
    [ 0.       0.       1.     ]])
b4=np.array([-0.03445  0.04664 -0.14286 -0.07143 -0.       0.21429  0.14286  0.07143])
r2p4=pc.Polytope(a4, b4)

---

With these regions, I am trying to find a common region where R1 and R2 intersect, in other sense trying to find solution space, through two methods.

Method: 1

#  R1 = pc.Region([r1p0, r1p1])
#  R2 = pc.Region([r2p0,  r2p1,  r2p2, ,  r2p3, ,  r2p4])

f=[]
for i in R1:
    for j in R2:
        if i.intersect(j):
            f.append(i.intersect(j))

This gives two polytopes as a result, shown in Fig.1

printing them gives

Single polytope 
  [[-0.57735 -0.57735 -0.57735] |    [[-0.11719]
   [ 0.57735  0.57735  0.57735] |     [ 0.11845]
   [-1.      -0.      -0.     ] x <=  [-0.2    ]
   [ 0.57735 -0.57735 -0.57735] |     [ 0.11682]
   [-0.      -0.      -1.     ] |     [ 0.     ]
   [-0.      -0.70711  0.70711]]|     [ 0.     ]]

Single polytope 
  [[-0.57735 -0.57735 -0.57735] |    [[-0.11693]
   [ 0.57735 -0.57735 -0.57735] |     [ 0.11375]
   [-0.57735  0.57735  0.57735] x <=  [-0.11249]
   [ 1.       0.       0.     ] |     [ 0.2    ]
   [-0.      -0.      -1.     ] |     [ 0.     ]
   [-0.      -0.70711  0.70711]]|     [ 0.     ]

Method 2

I wrote a function

def R_P_intersect(r1, r2):
    u = []

    for count, i in enumerate(r2):
        v =i.intersect(r1) # Intersection between Region 'R1' and Polytope 'i' in R2
        
        # Check if the intersection result is a Polytope
        if isinstance(v, pc.Polytope):
            if not pc.is_empty(v):
                u.append(v)  # Add non-empty Polytope to the list
                
        # If the intersection result is a Region, handle each Polytope in it
        elif isinstance(v, pc.Region):
            for k in v:
                if not pc.is_empty(k):
                    u.append(k)

    return pc.Region(u)

# Calling the fn to compute the intersection
case1 = R_P_intersect (R1, R2)  # look at position of R1 and R2
case2 = R_P_intersect(R2,, R1)  # look at position of R1 and R2

in case1, I get a region that contains two polytopes and both are identical

for i in case1:
    print(i)

Single polytope 
  [[ 0.57735 -0.57735 -0.57735] |    [[ 0.11375]
   [-0.57735  0.57735  0.57735] |     [-0.11249]
   [ 1.       0.       0.     ] x <=  [ 0.2    ]
   [-0.      -0.70711  0.70711] |     [ 0.     ]
   [-0.57735 -0.57735 -0.57735] |     [-0.11693]
   [-0.      -0.      -1.     ]]|     [ 0.     ]]

Single polytope 
  [[ 0.57735 -0.57735 -0.57735] |    [[ 0.11375]
   [-0.57735  0.57735  0.57735] |     [-0.11249]
   [ 1.       0.       0.     ] x <=  [ 0.2    ]
   [-0.      -0.70711  0.70711] |     [ 0.     ]
   [-0.57735 -0.57735 -0.57735] |     [-0.11693]
   [-0.      -0.      -1.     ]]|     [ 0.     ]]

The visualization, Fig2

in case2 i get a result similar to Method 1. Also, the A and b matrices of the polytope are the same as in Method 1.

Why does the position of R1 and R2 affect the final result?

Method 3 (using & directly):

---

I know & can be used to do the intersection process., I did as

case3= R1 & R2

This results in

Single polytope 
  [[ 0.57735 -0.57735 -0.57735 ] |    [[ 0.11375]
   [-0.57735  0.57735  0.57735 ] |     [-0.11249]
   [ 1.       0.       0.                     ] x <=  [ 0.2    ]
   [-0.      -0.70711  0.70711     ] |     [ 0.     ]
   [-0.57735 -0.57735 -0.57735] |     [-0.11693]
   [-0.      -0.      -1.                   ]]|     [ 0.     ]]

Single polytope 
  [[ 0.57735 -0.57735 -0.57735] |    [[ 0.11375]
   [-0.57735  0.57735  0.57735 ] |     [-0.11249]
   [ 1.       0.       0.                     ] x <=  [ 0.2    ]
   [-0.      -0.70711  0.70711      ] |     [ 0.     ]
   [-0.57735 -0.57735 -0.57735 ] |     [-0.11693]
   [-0.      -0.      -1.                    ]]|     [ 0.     ]]

This is the same as in Fig 2. The trial case4= R2 & R1 gives the same result, identical to case3.

So the intersect and & operations behave very differently. The operation R.intersect(P) does not give any warning or message about the position of the input argument.

So the explanation about & and intersect would be handy and will help me to understand their implementation.

Thank you

---

Best regards
Muthu Vallinayagam
Research, Institute of Experimental Physics
TU Freiberg
Germany