This microsite is associated to the paper
Actions of the hyperoctahedral group to compute minimal contractors.
The hyperoctahedral group Bn is the group of symmetries of a hypercube of R^n. For instance permutations, or symmetries along each canonical plane of all belong to Bn.
Now, many sets of equations we meet in our applications contain hyperoctahedral symmetries.
This is the case of the addition constraint x1+x2=x3 or the multiplication x1x2=x3.
The addition and multiplication constraints for matrices such as A+B=C or AB=C also contain many hyperoctahedral symmetries.
In robotics, many specific geometrical constraints such as for instance constraints involving distances or angles have hyperoctahedral symmetries.
In this paper, I present an algorithm which will allow us to build minimal contractors associated with constraints containing hyperoctahedral symmetries.
from sympy import *
x1,x2,x3,x4,x5,x6 = symbols('x1 x2 x3 x4 x5 x6')
G1=groebner([x3*x1-x4*x2-x5 ,x4*x1+x3*x2-x6,x3**2+x4**2-1],x1,x2,x3,x4,x5,x6, order='lex')
G2=groebner([x3*x5+x4*x6-x1,-x4*x5+x3*x6-x2,x3**2+(-x4)**2-1],x1,x2,x3,x4,x5,x6, order='lex')
print(G1)
print(G2)
print(G1==G2)
from codac import *
import itertools
def ToSign(L):
La=[]
for Z in L:
s=""
for z in Z:
if (z==Interval(0,oo)): s=s+"+"
elif (z==Interval(-oo,0)): s = s + "-"
else: s=s+"0"
La.append(s)
return La
def diff(L1, L2):
L = [x for x in L1 if not(x in L2)]
return L
def consistent(X1,X2,X3):
X3=X3&(X1*X2)
return not(X3.diam()==0)
def BuildCoverX(): # compute a cover of X with Cartans
LX= [I for I in list(itertools.product([Interval(-oo,0),Interval(0,oo)], repeat=dimx)) if consistent(*I)]
print("Number of boxes of the paving covering X: ",len(LX))
print("LX=", ToSign(LX))
return LX
def Validate_Sequence(Ls):
Ak=[[Interval(0,oo)]*dimx]
for sk in Ls: Ak=Ak+[sk(*Z) for Z in Ak]
return (len(diff(ToSign(CoverX),ToSign(Ak)))==0)
def FindSequence():
lmax=0
notfound=True
while notfound:
lmax = lmax + 1
print('Search for a sequence of length ',lmax)
for I in list(itertools.permutations(list(range(0, Nsym)), lmax)):
if Validate_Sequence([Ls[I[i]] for i in range(0, len(I))]):
notfound=False
print('found',I)
def s1(X1, X2, X3): return -X1, -X2, X3
def s2(X1, X2, X3): return X2, X1, X3 # not needed, but will be checked by the algo
def s3(X1, X2, X3): return -X1, X2, -X3
Nsym=3
dimx=3
CoverX=BuildCoverX()
Ls=[s1,s2,s3]
FindSequence()
from vibes import vibes
from codac import *
def union_tuple(L1, L2):
return tuple([x[0] | x[1] for x in tuple(zip(L1, L2))])
def Cmult(X,Y,Z):
def C0(X, Y, Z):
X,Y,Z=X&Interval(10**-10,oo),Y&Interval(10**-10,oo),Z&Interval(10**-10,oo)
return X & Interval(Z.lb()/Y.ub(),Z.ub()/Y.lb()), Y & Interval(Z.lb()/X.ub(),Z.ub()/X.lb()), Z&Interval(X.lb()*Y.lb(),X.ub()*Y.ub())
def SymCtr(X, Y, Z, C, s, _s):
return union_tuple(C(X,Y,Z),s(*C(*_s(X,Y,Z))))
def ac(C1,s,_s):
return lambda X, Y,Z : SymCtr(X,Y,Z, C1, s,_s)
def s1(X1, X2, X3): return -X1, -X2, X3
def s3(X1, X2, X3): return -X1, X2, -X3
return ac(ac(C0,s1,s1),s3,s3)(X,Y,Z)
class CtcMult(Ctc):
def __init__(C): Ctc.__init__(C, 2)
def contract(C, X): X[0],X[1],_ = Cmult(X[0],X[1],Interval(-9,-2))
vibes.beginDrawing()
pySIVIA([[-10,10],[-10,10]], CtcMult(), 0.1, **{'color_out':'blue[cyan]', 'color_maybe':'white[yellow]'})
vibes.axisAuto()
vibes.setFigureSize(500,500)
vibes.endDrawing()
from codac import *
import itertools
def ToSign(L):
La=[]
for Z in L:
s=""
for z in Z:
if (z==Interval(0,oo)): s=s+"+"
elif (z==Interval(-oo,0)): s = s + "-"
else: s=s+"0"
La.append(s)
return La
def diff(L1, L2): return [x for x in L1 if not(x in L2)]
def consistent(X1,X2,X3,X4,X5,X6):
X5=X5&(X3*X1-X4*X2)
X6=X6&(X4*X1+X3*X2)
return not((X5.diam()==0) | (X6.diam()==0))
def BuildCoverX(): # compute a cover of X with Cartans
LX= [I for I in list(itertools.product([Interval(-oo,0),Interval(0,oo)], repeat=dimx)) if consistent(*I)]
print("Number of boxes of the paving covering X: ",len(LX))
print("LX=", ToSign(LX))
return LX
def Validate_Sequence(Ls):
Ak=[[Interval(0,oo)]*dimx]
for sk in Ls: Ak=Ak+[sk(*Z) for Z in Ak]
return (len(diff(ToSign(CoverX),ToSign(Ak)))==0)
def FindSequence():
lmax=0
notfound=True
while notfound:
lmax = lmax + 1
print('Search for a sequence of length ',lmax)
for I in list(itertools.permutations(list(range(0, Nsym)), lmax)):
if Validate_Sequence([Ls[I[i]] for i in range(0,len(I))]):
notfound=False
print('found',I)
def s1(X1, X2, X3, X4, X5, X6): return [X5, X6, X3, -X4, X1, X2]
def s2(X1, X2, X3, X4, X5, X6): return [X2, -X1, -X4, X3, X5, X6]
def s3(X1, X2, X3, X4, X5, X6): return [X1, -X2, X3, -X4, X5, -X6]
def s4(X1, X2, X3, X4, X5, X6): return [-X1, -X2, -X3, -X4, X5, X6]
def s5(X1, X2, X3, X4, X5, X6): return [-X1, X2, X3, -X4, -X5, X6]
Nsym = 5
dimx = 6
CoverX=BuildCoverX()
Ls=[s1,s2,s3,s4,s5]
FindSequence()
from vibes import vibes
from roblib import *
from codac import *
x1min,x1max,x2min,x2max=-30,30,-5,30
dt = 0.1
a=[[-10,0],[5,0],[10,0]]
np.random.seed(0)
def generate_data():
ax = init_figure(x1min, x1max, x2min, x2max)
x = [10, 20]
v = [40, 20]
err_αdot = 10.0
err_d = 1 # >1 : multiplicative error
err_dd = 1
clear(ax)
D,dD,α,αdot=[],[],[],[]
x1,x2=x[0],x[1]
v1,v2=v[0],v[1]
draw_tank(array([[x1],[x2],[arctan2(v2,v1)]]),'red',0.5,1)
for i in range(0,len(a)):
a1,a2=a[i][0],a[i][1]
α_=arctan2(x2-a2,x1-a1)
draw_tank(array([[a1],[a2],[α_]]),'blue',0.4,1)
d_=np.sqrt((x1-a1)**2+(x2-a2)**2)
D.append(np.round(d_)+err_d*Interval(-1,1))
dd_=np.round(np.cos(α_)*v1+np.sin(α_)*v2)
dD.append(dd_+err_dd*Interval(-1,1))
αdot_=(-np.sin(α_) * v1 + np.cos(α_) * v2) / d_
αdot.append(Interval(-100,100))
α.append(Interval(np.floor(α_),np.ceil(α_)))
V10 = Interval(-100, 100)
V20 = Interval(- 100, 100)
print("V10=",V10); print("V20=",V20); print("D=",D); print("dD=",dD); print("α=",α); print("αdot=",αdot)
return V10,V20,D,dD,α,αdot
generate_data()
V10= [-1000, 1000]
V20= [-1000, 1000]
D= [[28, 30], [20, 22], [19, 21]]
dD= [[42, 44], [29, 31], [20, 22]]
α= [[0, 1], [1, 2], [1, 2]]
αdot= [[-1000, 100], [-1000, 1000], [-1000, 1000]]
from rot import Crot
class CtcRadar(Ctc):
def __init__(C): Ctc.__init__(C, 2)
def contract(C, X):
V1=Interval(V10[0],V10[1])
V2=Interval(V20[0],V20[1])
for i in range(0,len(a)):
Di=Interval(D[i][0], D[i][1])
dDi=Interval(dD[i][0], dD[i][1])
αi=Interval(α[i][0], α[i][1])
αdoti=Interval(αdot[i][0], αdot[i][1])
C1i,C2i=cos(αi),sin(αi)
Wi2= dDi*αdoti
Ai=IntervalVector([Interval(a[i][0], a[i][0]), Interval(a[i][1], a[i][1])])
P=X-Ai
dDi,Wi2,C1i,C2i,V1,V2 = Crot(dDi,Wi2,C1i,C2i,V1,V2)
Zi=Interval(0,0)
Di,Zi, C1i,C2i,P[0],P[1] = Crot(Di,Zi, C1i,C2i,P[0],P[1] )
X[0]=X[0]&(P[0]+Ai[0])
X[1]=X[1]&(P[1]+Ai[1])
Cr=CtcRadar()
vibes.beginDrawing()
pySIVIA([[x1min,x1max],[x2min,x2max]], Cr, 0.5, **{'color_out':'blue[cyan]', 'color_maybe':'yellow[yellow]'})
vibes.endDrawing()
from vibes import vibes
from pyibex import *
from pyibex.geometry import SepPolygon
import numpy as np
def union_tuple(L1, L2):
return tuple([x[0] | x[1] for x in tuple(zip(L1, L2))])
def Crot(X1,X2,X3,X4,X5,X6):
Cnorm = CtcFwdBwd(Function('x1', 'x2', 'a1', 'a2', 'y1', 'y2', '(x1^2+x2^2)-(y1^2+y2^2)'))
Crot = CtcFwdBwd(Function('x1', 'x2', 'a1', 'a2', 'y1', 'y2', "(a1*x1-a2*x2-y1;a2*x1+a1*x2-y2)"))
Ca = CtcFwdBwd(Function('x1', 'x2', 'a1', 'a2', 'y1', 'y2', "a1^2+a2^2-1"))
def contract_sweep_X(X1, X2, X3,X4, X5, X6):
Rp=np.array([[X3.lb(),X4.ub()],[-X4.ub(),X3.lb()]])
Rq=np.array([[X3.ub(),X4.lb()],[-X4.lb(),X3.ub()]])
P = Rp @ [[X5.lb(), X5.ub(), X5.ub()], [X6.lb(), X6.lb(), X6.ub()]]
Q = Rq @ [[X5.ub(), X5.lb(), X5.lb()], [X6.ub(), X6.ub(), X6.lb()]]
S = SepPolygon([[Q[0, 2], P[1, 0]], [P[0, 0], P[1, 0]], [P[0, 1], P[1, 1]], [P[0, 2], P[1, 2]],
[P[0, 2], Q[1, 0]], [Q[0, 0], Q[1, 0]], [Q[0, 1], Q[1, 1]], [Q[0, 2], Q[1, 2]]])
Xa, Xb = IntervalVector([X1, X2]),IntervalVector([X1, X2])
S.separate(Xa, Xb)
return Xb[0], Xb[1]
def C0(X1, X2, X3, X4, X5, X6):
Z = IntervalVector([X1, X2, X3, X4, X5, X6])&IntervalVector(6,Interval(0,oo))
Cnorm.contract(Z)
Crot.contract(Z)
Ca.contract(Z)
X1, X2, X3,X4, X5, X6 = tuple(Z)
X1, X2 = contract_sweep_X(X1, X2, X3,X4, X5, X6)
X5, X6 = contract_sweep_X(X5, X6, X3,-X4, X1, X2)
return X1, X2, X3, X4, X5, X6
def ac(C,s,_s): return lambda X1,X2,X3,X4,X5,X6 : \
union_tuple(C(X1, X2, X3,X4, X5, X6),s(*C(*_s(X1, X2, X3,X4, X5, X6))))
def s1(X1, X2, X3,X4, X5, X6): return X5,X6,X3,-X4, X1, X2
def s2(X1, X2, X3,X4, X5, X6): return X2,-X1, -X4,X3, X5, X6
def _s2(X1, X2, X3,X4, X5, X6): return -X2, X1,X4,-X3, X5, X6
def s3(X1, X2, X3, X4, X5, X6): return X1,-X2, X3, -X4, X5, -X6
def s4(X1, X2, X3, X4, X5, X6): return -X1,-X2, -X3, -X4, X5, X6
def s5(X1, X2, X3, X4, X5, X6): return -X1, X2, X3, -X4, -X5, X6
return ac(ac(ac(ac(ac(C0,s1,s1),s2,_s2),s3,s3),s4,s4),s5,s5)(X1,X2,X3,X4,X5,X6)
class CtcRot(Ctc):
def __init__(C): Ctc.__init__(C, 2)
def contract(C, X): X[0], X[1],_,_,_,_ = Crot(X[0], X[1], cos(θ0),sin(θ0), Y10, Y20)