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)
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()
from findsym import count_rotate,count_2disks
import numpy as np
def gene_psi(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_psi_rotate():
S=count_rotate()
m=2
gene_psi(S,m)
def gene_psi_2disks():
S=count_2disks()
m=1
gene_psi(S,m)
print("\n Generation of the function psi for the constraint two_disks")
gene_psi_2disks()
print("\n Generation of the function psi for the constraint Rotate")
gene_psi_rotate()
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
psi = { #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,psi):
sym=[]
for I in product(*sgn(Q)): sym.append(psi[(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,psi)
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(theta0.lb(),theta0.ub(),0.1):
draw_polygon_rot(P,-t,'black')
draw_polygon_rot(P,-theta0.ub(),'blue')
draw_polygon_rot(P,0,'#111111[#3333FFAA]')
vibes.endDrawing()
theta0,Y10,Y20=Interval(3,4),Interval(-10,10),Interval(5,12)
Q=[cos(theta0),sin(theta0),Y10,Y20]
print(Q)
siviaP(Q)
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(theta0.lb(),theta0.ub(),0.1):
draw_polygon_rot(P,-t,'black')
draw_polygon_rot(P,-theta0.ub(),'blue')
draw_polygon_rot(P,0,'#FF1111[#FF11DDAA]')
vibes.endDrawing()
theta0=Interval(0.5,1.5)
Y10a,Y20a=Interval(5,10),Interval(-8,8)
Y10b,Y20b=Interval(-4,11),Interval(8,10)
Qa=[cos(theta0),sin(theta0),Y10a,Y20a]
Qb=[cos(theta0),sin(theta0),Y10b,Y20b]
siviaP(Qa,Qb)
from vibes import vibes
from roblib import *
from pyibex import *
from quotient import sep_angle,RunsiviaP,draw_polygon_rot,rot
def siviaP(theta0,Y10,Y20):
S=[]
for i in range(N):
Q=[cos(theta0[i]),-sin(theta0[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,ntheta=1,0.5
thetatruth=[1.2,2.4,3.1,0.1,-1.3,-2.1]
N=len(thetatruth)
thetam=[randround(thetatruth[i],ntheta) for i in range(N)]
print("thetam",thetam)
theta0=[thetam[i]+Interval(-ntheta,ntheta) for i in range(N)]
print(theta0)
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],-thetatruth[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(theta0,Y10,Y20)
a=180/3.14
for i in range(N):
vibes.drawAUV(Pos[i][0], Pos[i][1],a*thetatruth[i], 2, 'black[yellow]')
vibes.drawPoint(Vearth[0],Vearth[1],2,'red[red]')
vibes.endDrawing()