In [14]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
%matplotlib inline
In [15]:
def plotellipse(ax, bmaj, bmin, bpa, **kwargs):
testphi = np.linspace(0, 2 * np.pi, 1001)
x = bmaj * np.cos(testphi)
y = bmin * np.sin(testphi)
xr = x * np.cos(bpa) - y * np.sin(bpa)
yr = x * np.sin(bpa) + y * np.cos(bpa)
ax.plot(xr, yr, **kwargs)
In [16]:
def boundingcircle(bmaj, bmin, bpa):
thisone = np.argmax(bmaj)
return(bmaj[thisone], bmaj[thisone], bpa[thisone])
In [17]:
def PtoA(bmaj, bmin, bpa):
A = np.zeros((2,2))
A[0, 0] = np.cos(bpa)**2 / bmaj**2 + np.sin(bpa)**2 / bmin**2
A[1, 0] = np.cos(bpa) * np.sin(bpa) * (1 / bmaj**2 - 1 / bmin**2)
A[0, 1] = A[1, 0]
A[1, 1] = np.sin(bpa)**2 / bmaj**2 + np.cos(bpa)**2 / bmin**2
return A
def BinsideA(B, A):
try:
CD = np.linalg.cholesky(B - A) #np.dot(A, A))
return True
except:
return False
def myobjective_regularized(p, bmajvec, bminvec, bpavec):
A = PtoA(*p)
test = np.zeros_like(bmajvec)
for idx, (bmx, bmn, bp) in enumerate(zip(bmaj, bmin, bpa)):
test[idx] = BinsideA(PtoA(bmx, bmn, bp), A)
obj = 1 / np.linalg.det(A)
if np.all(test):
return obj
else:
return obj * 1e30
In [19]:
# Generate some random ellipses
nEllipse=10
bmaj = np.abs(np.random.randn(nEllipse))+2
bmin = np.abs(np.random.randn(nEllipse))+1
bpa = np.random.rand(nEllipse) * 2 * np.pi
bmajmax, bminmax, bpamax = boundingcircle(bmaj, bmin, bpa)
# Initialize and optimize one way
p0 = (1.1* bmajmax, bminmax, np.pi/3.)
optdict = {'maxiter': 5000, 'ftol': 1e-11, 'maxfev':5000}
result = opt.minimize(myobjective_regularized, p0, method='Nelder-Mead',
args = (bmaj, bmin, bpa), options=optdict, tol=1e-9)
print result.viewitems()
# Initialize and optimize another
p0 = (bmajmax, bminmax, -np.pi/4)
result2 = opt.minimize(myobjective_regularized, p0, method='Nelder-Mead',
args = (bmaj, bmin, bpa), options=optdict, tol=1e-9)
print result2.viewitems()
# see what we hath wrought.
fig, ax = plt.subplots(1,1)
bmajmax, bminmax, bpamax = boundingcircle(bmaj, bmin, bpa)
for bmx, bmn, bp in zip(bmaj, bmin, bpa):
plotellipse(ax, bmx, bmn, bp, color='blue')
plotellipse(ax, result.x[0], result.x[1], result.x[2], color='red')
plotellipse(ax, result2.x[0], result2.x[1], result2.x[2], color='orange', lw=2, alpha=0.5)
ax.set_aspect('equal')
dict_items([('status', 0), ('success', True), ('final_simplex', (array([[ 3.13165165, 3.45492934, 1.0574624 ],
[ 3.13165165, 3.45492934, 1.0574624 ],
[ 3.13165165, 3.45492934, 1.0574624 ],
[ 3.13165165, 3.45492934, 1.0574624 ]]), array([ 117.06450469, 117.06450469, 117.06450469, 117.06450469]))), ('nfev', 3515), ('fun', 117.06450469018093), ('x', array([ 3.13165165, 3.45492934, 1.0574624 ])), ('message', 'Optimization terminated successfully.'), ('nit', 1777)])
dict_items([('status', 0), ('success', True), ('final_simplex', (array([[ 3.39275737, 3.18972543, -0.91833718],
[ 3.39275737, 3.18972543, -0.91833718],
[ 3.39275737, 3.18972543, -0.91833718],
[ 3.39275737, 3.18972543, -0.91833718]]), array([ 117.11491471, 117.11491471, 117.11491471, 117.11491471]))), ('nfev', 395), ('fun', 117.11491470672625), ('x', array([ 3.39275737, 3.18972543, -0.91833718])), ('message', 'Optimization terminated successfully.'), ('nit', 208)])
