import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp
import scipy

np.random.seed(1)
b = np.array([[.5, 1.]]).T
α = 0.3
A = np.array([[1., α], [α, 1.5]])
A_inv = np.linalg.inv(A)
A_sqrt = scipy.linalg.sqrtm(A)
n_pts = 10000
pts_base = (np.random.rand(2, n_pts)-.5)*5
pts = pts_base+b
ε_1 = pts[:, [(pts-b)[:,i].T@A_inv@(pts-b)[:,i]<=1 for i in range(n_pts)]]
plt.scatter(*ε_1, c="C1", alpha=.5, label=r"$\mathcal{E} = \left \{ x | (x-x_c)^\mathsf{T}P^{-1}(x-x_c) \leq 1 \right \}$")
pts = pts_base[:, np.linalg.norm(pts_base, axis=0) <= 1]
ε_2 = b+A_sqrt@pts
plt.scatter(*ε_2, c="C2", alpha=.5, label=r"$\mathcal{E} = \left \{ x_c+Au \,|\, \| u \| \leq 1 \right \}$")
plt.legend()
plt.show()