This microsite is associated to the paper
L. Jaulin (2023).
Inner and outer characterization of the projection of polynomial equations using symmetries, quotients and intervals,
International Journal of Approximate Reasoning, Volume 159, August 2023.
Using symmetries, set quotient and interval analysis we can compute efficiently an inner and an outer approximation for the set P
which corresponds to the projection of another set X defined by polynomial equations.
from pyibex import *
from vibes import *
f = Function('x', 'y', '(x-2)^2 + (3*y+x-1)^2')
S1 = SepFwdBwd(f, sqr(Interval(0,4)))
s = Function('x', 'y', '(y, x)')
S2 = SepTransform(S1,s,s)
vibes.beginDrawing()
X0 = IntervalVector(2, [-3,7])
vibes.newFigure('X')
pySIVIA(X0, S1|S2, 0.01)
Code for computing sep(Bn(X)))
Click below to see/run/modify the online program :
from numpy import *
from math import factorial
import numpy as np
import sympy as sp
import itertools
def separable_only(Bn, m):
S = [s for s in Bn if (np.abs(np.prod(s[0:m])) == factorial(m))]
return (S)
def Bn(n): # all_OctaSym
N = [i for i in range(1, n + 1)]
LX = [I for I in itertools.permutations(N)]
LY = [I for I in list(itertools.product([-1, 1], repeat=n))]
L = []
for X in LX:
for Y in LY:
s = []
for x, y in zip(X, Y):
s.append(x * y)
L.append(s)
return L
def octasym(I):
return lambda L: [
np.sign(I[i]) * L[np.abs(I[i]) - 1] for i in range(0, len(I))
]
def diff(L1, L2):
L = [x for x in L1 if not (x in L2)]
return L
def vect2matsym(s):
n = len(s)
M = np.zeros([n, n])
for i in range(n):
e = np.sign(s[i])
j = np.abs(s[i]) - 1
M[j][i] = e
return M
def mat2vectsym(M):
n = len(M)
v = []
for j in range(n):
for i in range(n):
if not (M[i, j] == 0): v.append(np.int(np.sign(M[i, j]) * (i + 1)))
return (v)
def invertsym(v):
M = vect2matsym(v)
N = np.linalg.inv(M)
return mat2vectsym(N)
def mult2vectsym(v1, v2):
return [np.sign(v2[i]) * v1[np.abs(v2[i]) - 1] for i in range(len(v1))]
def Test_solution(Eq, X, Lsymb):
m = len(Eq)
bb = True
for j in range(m):
E = Eq[j]
for i in range(len(X)):
E = E.subs(Lsymb[i], X[i])
bb = bb and (sp.simplify(E) == 0)
return bb
def Sep_BnX(Bn, m, x0, S, Eq,
Lsymb): # x0 : known solution, S : known symmetries
U = []
i = 0
B = separable_only(Bn, m)
for I in B:
i = i + 1
s = octasym(I)
if Test_solution(Eq, s(x0), Lsymb):
U.append(I)
print("Sep_BnX in U, i=", i, " ", I)
U = diff(U, S)
bb = True
while bb == True:
bb = False
for s1 in S:
Z = []
for s2 in U:
s3 = mult2vectsym(s1, s2)
if s3 in S:
bb = True
Z.append(s2)
S.append(s2)
U = diff(U, Z)
if (len(U) > 0): print("Please, add solutions or symmetries")
return S, U
def count_2disks(
): #Mutiplication : ((x1-2)**2+(x2-3)**2)*((x1+2)**2+(x2+3)**2)
n, m = 2, 1
B2 = Bn(n)
print("Number of all Hyperoctahedral symmetries found : ", len(B2))
print("2**n*factorial(n) : ", 2**n * factorial(n))
x0 = [2 + 1 / sp.sqrt(4),
3 + sp.sqrt(3) / 2] # valid solutions to be given by the user
S = [[-1, -2], [1, 2]] # valid symmetries to be given by the user
x1, x2 = sp.symbols("x1 x2")
Lsymb = [x1, x2]
eq1 = sp.sympify("((x1-2)**2+(x2-3)**2-1)*((x1+2)**2+(x2+3)**2-1)")
Eq = [eq1]
S, U = Sep_BnX(B2, 1, x0, S, Eq, Lsymb)
return S
def count_rotate(): #rotate : X1*x2=x3
# x3*x1-x4*x2=x5
# x4*x1+x3*x2=x6
# x3^2+x4^2=1
# one solution (2,4,1/2,sqrt(3)/2,1-2*sqrt(3),2+sqrt(3)
n, m = 6, 2
B6 = Bn(n)
print("len(B6)", len(B6))
print(2**n * factorial(n))
x0 = [
2, 4, 1 / sp.sqrt(4),
sp.sqrt(3) / 2, 1 - 2 * sp.sqrt(3),
sp.sqrt(3) + 2
]
Lsymb = list(sp.symbols("x1 x2 x3 x4 x5 x6"))
eq1 = sp.sympify("x3*x1-x4*x2-x5")
eq2 = sp.sympify("x4*x1+x3*x2-x6")
eq3 = sp.sympify("x3^2+x4^2-1")
Eq = [eq1, eq2, eq3]
print(Eq)
S = [[2, -1, -4, 3, 5, 6], [1, -2, 3, -4, 5, -6],
[-1, -2, -4, 3, 6, -5]] # valid symmetries to be given by the user
S, U = Sep_BnX(B6, m, x0, S, Eq, Lsymb)
return S
print("Generate the separable symmetries for two-disks")
S = count_2disks()
print("Generate the separable symmetries for rotate")
S = count_rotate()
Code for generating the function psi for the rotate constraint and for the two-disks constraint
Click below to see/run/modify the online program :
from findsym import count_rotate,count_2disks
import numpy as np
def gene_ψ(S,m):
def φ(s): return [np.sign(s[i]) for i in range(m,l)]
A=[]
for s in S:
l=len(s)
a=φ(s)
if not (a in A):
A.append(a)
Str=""
for i in range(l-m-1): Str=Str+str(a[i])+ ","
Str="(" + Str+ str(a[l-m-1]) + ") : " + str(s)+ " ,"
print(Str)
def gene_ψ_rotate():
S=count_rotate()
m=2
gene_ψ(S,m)
def gene_ψ_2disks():
S=count_2disks()
m=1
gene_ψ(S,m)
print("\n Generation of the function ψ for the constraint two_disks")
gene_ψ_2disks()
print("\n Generation of the function ψ for the constraint Rotate")
gene_ψ_rotate()
Code for the rotating rectangle
Click below to see/run/modify the online program :
from vibes import vibes
from roblib import *
from pyibex import *
from rot import sponge0, octasym
from findsym import invertsym
from itertools import product
class CtcSponge0(Ctc): #contract X1,X2 on the positive quarter
def __init__(C,Q,outer):
Ctc.__init__(C, 2)
C.outer=outer
C.Q=Q
def contract(C,P):
if not(P.is_empty()):
P[0],P[1] = sponge0(P[0],P[1],*C.Q,C.outer)
return P
class CtcSym(Ctc):
def __init__(C,C0,s):
C.m=len(s)-len(C0.Q)
Ctc.__init__(C, C.m)
C.outer=C0.outer
C.sp=s[0:C.m]
C.sq=[np.sign(s[i])*(np.abs(s[i])-C.m) for i in range(C.m,len(s))]
C.C0=C0
C.Q=C0.Q.copy()
def contract(C,P):
if not(P.is_empty()):
C.C0.Q=octasym(C.sq)(C.Q)
P1=octasym(C.sp)(P)
P2=C.C0.contract(IntervalVector(P1))
P=IntervalVector(octasym(invertsym(C.sp))(P2))
return P
def sgn(Q): #interval extension of the vector sign
def sgn(Q): # return all possible signs an interval.
L=[]
if (Q.lb()<0): L.append(-1)
if (Q.ub()>0): L.append(1)
return L
l=len(Q)
L=[sgn(Q[i]) for i in range(l)]
return L
ψ = { #obtained via the program gene_psi
( -1 , 1 , 1 , 1 ) : [2, -1, -4, 3, 5, 6] ,
( 1 , -1 , 1 , -1 ) : [1, -2, 3, -4, 5, -6] ,
( -1 , -1 , 1 , 1 ) : [-1, -2, -3, -4, 5, 6] ,
( 1 , -1 , -1 , 1 ) : [-1, 2, 3, -4, -5, 6] ,
( -1 , 1 , 1 , -1 ) : [-1, -2, -4, 3, 6, -5] ,
( 1 , 1 , 1 , 1 ) : [1, 2, 3, 4, 5, 6] ,
( 1 , 1 , 1 , -1 ) : [2, -1, 3, 4, 6, -5] ,
( 1 , -1 , 1 , 1 ) : [-2, 1, 4, -3, 5, 6] ,
( 1 , 1 , -1 , 1 ) : [2, 1, 4, 3, -5, 6] ,
( 1 , 1 , -1 , -1 ) : [-1, -2, 3, 4, -5, -6] ,
( 1 , -1 , -1 , -1 ) : [-2, -1, 3, -4, -6, -5] ,
( -1 , 1 , -1 , -1 ) : [-2, 1, -4, 3, -5, -6] ,
( -1 , -1 , 1 , -1 ) : [2, 1, -4, -3, 5, -6] ,
( -1 , 1 , -1 , 1 ) : [1, -2, -3, 4, -5, 6] ,
( -1 , -1 , -1 , -1 ) : [1, 2, -3, -4, -5, -6] ,
( -1 , -1 , -1 , 1 ) : [-2, -1, -4, -3, -5, 6]
}
def sep_sym(CtcC0,Q,ψ):
sym=[]
for I in product(*sgn(Q)): sym.append(ψ[(I)])
_C0=CtcC0(Q,True)
C0=CtcC0(Q,False)
Seps=[]
for s in sym:
_s=invertsym(s)
C1=CtcSym(C0,_s)
_C1=CtcSym(_C0,_s)
S= SepCtcPair(_C1,C1)
Seps.append(S)
return SepUnion(Seps)
def sep_angle(Q):
return sep_sym(CtcSponge0,Q,ψ)
def rot(x,y,t):
x1=np.cos(t)*x-np.sin(t)*y
y1=np.sin(t)*x+np.cos(t)*y
return [x1,y1]
def draw_polygon_rot(P,t,col='darkblue'):
R=[rot(A[0],A[1],t) for A in P]
vibes.drawPolygon(R, col)
def RunsiviaP(S):
p1min,p1max,p2min,p2max=-20,20,-20,20
vibes.beginDrawing()
vibes.newFigure('Solution')
params = { "x": 0,"y": 0,"width": 800,"height": 800}
vibes.setFigureProperties(params)
pySIVIA([[p1min,p1max],[p2min,p2max]], S, 1.2)
for i in range(-21,21,1):
vibes.drawLine([[-0.3,i],[0.3,i]],'black')
vibes.drawLine([[i,-0.3],[i,0.3]],'black')
vibes.drawLine([[p1min,0],[p1max,0]],'black')
vibes.drawLine([[0,p2min],[0,p2max]],'black')
def siviaP(Q):
S=sep_angle(Q)
RunsiviaP(S)
A1=[Q[2].lb(),Q[3].lb()]
A2=[Q[2].lb(),Q[3].ub()]
A3=[Q[2].ub(),Q[3].ub()]
A4=[Q[2].ub(),Q[3].lb()]
P=[A1,A2,A3,A4]
if True: #if we want to draw the polygon
for t in arange(θ0.lb(),θ0.ub(),0.1):
draw_polygon_rot(P,-t,'black')
draw_polygon_rot(P,-θ0.ub(),'blue')
draw_polygon_rot(P,0,'#111111[#3333FFAA]')
vibes.endDrawing()
θ0,Y10,Y20=Interval(3,4),Interval(-10,10),Interval(5,12)
Q=[cos(θ0),sin(θ0),Y10,Y20]
print(Q)
siviaP(Q)
Code for the Workspace characterization
Click below to see/run/modify the online program :
from vibes import vibes
from roblib import *
from pyibex import *
from quotient import sep_angle,RunsiviaP,draw_polygon_rot
def siviaP(Q1,Q2):
S1=sep_angle(Q1)
S2=sep_angle(Q2)
RunsiviaP(S1|S2)
A1=[Q1[2].lb(),Q1[3].lb()]
A2=[Q1[2].lb(),Q1[3].ub()]
A3=[Q2[2].lb(),Q2[3].lb()]
A4=[Q2[2].lb(),Q2[3].ub()]
A5=[Q2[2].ub(),Q2[3].ub()]
A6=[Q2[2].ub(),Q2[3].lb()]
A7=[Q1[2].ub(),Q1[3].ub()]
A8=[Q1[2].ub(),Q1[3].lb()]
P=[A1,A2,A3,A4,A5,A6,A7,A8]
print(P)
if True:
for t in arange(θ0.lb(),θ0.ub(),0.1):
draw_polygon_rot(P,-t,'black')
draw_polygon_rot(P,-θ0.ub(),'blue')
draw_polygon_rot(P,0,'#FF1111[#FF11DDAA]')
vibes.endDrawing()
θ0=Interval(0.5,1.5)
Y10a,Y20a=Interval(5,10),Interval(-8,8)
Y10b,Y20b=Interval(-4,11),Interval(8,10)
Qa=[cos(θ0),sin(θ0),Y10a,Y20a]
Qb=[cos(θ0),sin(θ0),Y10b,Y20b]
siviaP(Qa,Qb)
Code for the speed estimator
Click below to see/run/modify the online program :
from vibes import vibes
from roblib import *
from pyibex import *
from quotient import sep_angle,RunsiviaP,draw_polygon_rot,rot
def siviaP(θ0,Y10,Y20):
S=[]
for i in range(N):
Q=[cos(θ0[i]),-sin(θ0[i]),Y10[i],Y20[i]]
S1=sep_angle(Q)
S.append(S1)
Sq=SepQInter(S)
Sq.q=3
RunsiviaP(Sq)
def randround(x,e):
return 2*e*round(x/(2*e))-e
ny,nθ=1,0.5
θtruth=[1.2,2.4,3.1,0.1,-1.3,-2.1]
N=len(θtruth)
θm=[randround(θtruth[i],nθ) for i in range(N)]
print("θm",θm)
θ0=[θm[i]+Interval(-nθ,nθ) for i in range(N)]
print(θ0)
Pos=[[-3,5],[-13,-5],[5,-10],[-7,10],[-10,-10],[10,-5]]
Vearth=[10,10]
Y10,Y20=[],[]
Y1m,Y2m=[],[]
for i in range(N):
Vi=rot(Vearth[0],Vearth[1],-θtruth[i])
Y1m.append(randround(Vi[0],ny))
Y2m.append(randround(Vi[1],ny))
Y10=[Y1m[i]+Interval(-ny,ny) for i in range(N)]
Y20=[Y2m[i]+Interval(-ny,ny) for i in range(N)]
print("\nY1m",Y1m)
print("Y10",Y10)
print("\nY2m",Y2m)
print("Y20",Y20)
siviaP(θ0,Y10,Y20)
a=180/3.14
for i in range(N):
vibes.drawAUV(Pos[i][0], Pos[i][1],a*θtruth[i], 2, 'black[yellow]')
vibes.drawPoint(Vearth[0],Vearth[1],2,'red[red]')
vibes.endDrawing()