<jupyter_start><jupyter_text>Example 5: Zonotope Containment: Empirial Study of Necessity Gaps<jupyter_code>import numpy as np
import pypolycontain as pp
import time
import pydrake.solvers.mathematicalprogram as MP
# use Gurobi solver
import pydrake.solvers.gurobi as Gurobi_drake
gurobi_solver=Gurobi_drake.GurobiSolver()
np.random.seed(0)
def alpha_necessary_old(Z_i,Z_o):
    program=MP.MathematicalProgram()
    zeta=program.NewContinuousVariables(Z_o.G.shape[1],2**Z_i.G.shape[1],"zeta")
    beta=program.NewContinuousVariables(1,"beta")
    V=pp.vcube(Z_i.G.shape[1])
    for i in range(V.shape[0]):
        program.AddLinearEqualityConstraint( Aeq=Z_o.G, beq=np.dot(Z_i.G,V.T)[:,i], vars= zeta[:,i]  )
    program.AddLinearConstraint(  np.less_equal(zeta,beta*np.ones(zeta.shape),dtype='object').flatten() )
    program.AddLinearConstraint(  np.less_equal(-zeta,beta*np.ones(zeta.shape),dtype='object').flatten() )
    program.AddLinearCost(np.eye(1),np.zeros((1)),beta)
    result=gurobi_solver.Solve(program,None,None)
    if result.is_success():
        alpha=1/result.GetSolution(beta)[0]
        return alpha
    else:
        print("optimization failed")
        
def alpha_necessary_older(Z_i,Z_o):
    alpha_min=np.inf
    V=pp.vcube(Z_i.G.shape[1])
    B=np.dot(Z_i.G,V.T)
    for i in range(V.shape[0]):
        program=MP.MathematicalProgram()
        zeta=program.NewContinuousVariables(Z_o.G.shape[1],"zeta")
        beta=program.NewContinuousVariables(1,"beta")
        program.AddLinearEqualityConstraint( Aeq=Z_o.G, beq=B[:,i:i+1], vars= zeta  )
#         for j in range(zeta.shape[0]):
#             program.AddLinearConstraint(  zeta[j]-beta, -np.inf, 0 )
#             program.AddLinearConstraint(  -zeta[j]+ beta, -np.inf, 0 )
        program.AddLinearConstraint(  np.less_equal(zeta,beta*np.ones(zeta.shape),dtype='object').flatten() )
        program.AddLinearConstraint(  np.less_equal(-zeta,beta*np.ones(zeta.shape),dtype='object').flatten() )
        program.AddLinearCost(np.eye(1),np.zeros((1)),beta)
        result=gurobi_solver.Solve(program,None,None)
        if result.is_success():
            alpha=1/result.GetSolution(beta)[0]
            alpha_min=min(alpha,alpha_min)
        else:
            print("optimization failed")
    return alpha_min

def alpha_necessary(Z_i,Z_o):
    alpha_min=np.inf
    V=pp.vcube(Z_i.G.shape[1])
    B=np.dot(Z_i.G,V.T)
    for i in range(V.shape[0]):
        program=MP.MathematicalProgram()
        zeta=program.NewContinuousVariables(Z_o.G.shape[1],"zeta")
        alpha=program.NewContinuousVariables(1,"alpha")
        Aeq=np.hstack(( Z_o.G, -B[:,i:i+1]))
        program.AddLinearEqualityConstraint( Aeq=np.hstack(( Z_o.G, -B[:,i:i+1])), \
                                            beq=np.zeros((Z_o.G.shape[0])),\
                                            vars= np.hstack((zeta,alpha))  )
        program.AddBoundingBoxConstraint(-1,1,zeta)
        program.AddLinearCost(-np.eye(1),np.zeros((1)),alpha)
        result=gurobi_solver.Solve(program,None,None)
        if result.is_success():
            alpha=result.GetSolution(alpha)[0]
            alpha_min=min(alpha,alpha_min)
        else:
            print("optimization failed")
    return alpha_min
        
        
def alpha_sufficient(Z_i,Z_o):
    program=MP.MathematicalProgram()
    beta=program.NewContinuousVariables(1,"beta")
    circumbody=pp.to_AH_polytope(Z_o)
    parametric_circumbody=circumbody.copy()
    parametric_circumbody.P.h=circumbody.P.h*beta
    Theta,*_=pp.subset(program,Z_i,parametric_circumbody,k=-1)
    program.AddLinearCost(np.eye(1),np.zeros((1)),beta)
    result=gurobi_solver.Solve(program,None,None)
    if result.is_success():
        alpha=1/result.GetSolution(beta)[0]
        return alpha
    else:
        print("optimization failed")
N,n_max=100,20
D=10 # Percent to report
k,j=0,0 # Counter
d,j={},0
Q_o={ 3:[5,10,15,20], 5:[10,15,20], 10:[15,20], 15:[20] }
for n in Q_o.keys():
    print("="*20,"\n","n=",n)
    for q_o in Q_o[n]:
        d[n,q_o]=np.zeros(N)
        for i in range(N):
            j+=1
            print(n,q_o,i)
#             q_i=np.random.randint(n,q_o+1)
            z_i=pp.zonotope(G=np.eye(n))
            z_o=pp.zonotope(G=np.random.normal(size=(n,q_o)))
#             x=time.time()
            a_suf=alpha_sufficient(z_i,z_o)
#             print("sufficient:",time.time()-x)
#             x=time.time()
            a_nec=alpha_necessary(z_i,z_o)
#             print("necessary:",time.time()-x)
            d[n,q_o][i]=1-a_suf/a_nec
#             print("gap:",n,q_o,"=",d[n,q_o][i])
            j+=1
print("Done")
for key in d.keys():
    print(key[0],' & ',key[1],' & ',np.sum(d[key]<0.001), ' & ', np.sum(d[key]<0.01), ' & ',\
         np.sum(d[key]<0.05), ' & ', np.sum(d[key]<0.05), np.round(np.max(d[key]),3) , ' \\\ \\hline')<jupyter_output>3  &  5  &  97  &  97  &  100  &  0.039  \\ \hline
3  &  10  &  89  &  94  &  99  &  0.051  \\ \hline
3  &  15  &  94  &  96  &  100  &  0.041  \\ \hline
3  &  20  &  94  &  98  &  100  &  0.024  \\ \hline
5  &  10  &  82  &  86  &  96  &  0.086  \\ \hline
5  &  15  &  77  &  81  &  98  &  0.077  \\ \hline
5  &  20  &  60  &  68  &  92  &  0.104  \\ \hline
10  &  15  &  66  &  76  &  96  &  0.095  \\ \hline
10  &  20  &  29  &  41  &  70  &  0.177  \\ \hline
15  &  20  &  58  &  71  &  94  &  0.1  \\ \hline
{'fraction_non_alphanumeric': 0.16967777984727136, 'fraction_numerical': 0.04321102626187372, 'mean_word_length': 2.6505778382053027, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '292855', 'license': 'AFL-1.1', 'n_tokens_mistral': 2256, 'n_tokens_neox': 2132, 'n_words': 280, 'path': 'doc/source/Example 5.ipynb', 'repo_name': 'sadraddini/pypolycontain', 'size': '5369'}