Commit cbd95c63 authored by rfetick's avatar rfetick
Browse files

compute moffat derivatives

parent 5d13de5e
......@@ -258,34 +258,42 @@ def moffat(XY, param, norm=None, removeInside=0, grad=False):
ax,ay,theta,beta,cx,cy = param
if grad:
u,ug = reduced_coord(XY,ax,ay,theta,cx,cy,grad=True)
g = _np.zeros((4,)+u.shape)
u,du = reduced_coord(XY,ax,ay,theta,cx,cy,grad=True)
else:
u = reduced_coord(XY,ax,ay,theta,cx,cy,grad=False)
V = (1.0+u)**(-beta)
E = 1.0
V = (1.0+u)**(-beta) # Moffat shape
E = 1.0 # normalization factor (eventually up to infinity)
F = 1.0 # normalization factor (eventually up to a limited radius)
if grad:
dVdu = -beta*V/(1.0+u)
dV = _np.zeros((4,)+u.shape)
for i in range(3): dV[i,...] = dVdu*du[i,...]
dV[3,...] = -V*_np.log(1.0+u)
if norm is None:
if grad:
dVdu = -beta*V/(1.0+u)
g[0,...] = dVdu*ug[0,...]
g[1,...] = dVdu*ug[1,...]
g[2,...] = dVdu*ug[2,...]
g[3,...] = -V*_np.log(1.0+u)
return V,g
return V,dV
else: # norm can be float or np.inf
if (beta<=1) and (norm==_np.inf): raise ValueError("Cannot compute Moffat energy for beta<=1")
if beta==1: raise ValueError("Energy computation for beta=1.0 not implemented yet. Sorry!")
E = (beta-1) / (_np.pi*ax*ay)
Fout = 1 - (1 + (norm**2) / (ax*ay))**(1-beta)
Fin = 1 - (1 + (removeInside**2) / (ax*ay))**(1-beta)
E = E / (Fout-Fin)
if grad:
# TODO: finish gradient computation
raise ValueError("Gradient not implemented yet!")
return V*E
Fout = (1 + (norm**2)/(ax*ay))**(1-beta)
Fin = (1 + (removeInside**2)/(ax*ay))**(1-beta)
F = 1/(Fin-Fout)
if grad: # TODO: debug gradient for norm=float
dE = [-E/ax,-E/ay,0,E/(beta-1)]
k = (1-beta)*Fout/(1 + (norm**2)/(ax*ay))*(norm**2)/(ax*ay)
dFout = _np.array([-k/ax,-k/ay,0,-Fout*_np.log(1 + (norm**2)/(ax*ay))]) if norm<_np.inf else _np.zeros(4)
k = (1-beta)*Fin/(1 + (removeInside**2)/(ax*ay))*(removeInside**2)/(ax*ay)
dFin = _np.array([-k/ax,-k/ay,0,-Fin*_np.log(1 + (removeInside**2)/(ax*ay))])
dF = -(dFin-dFout)/(Fin-Fout)**2
dm = 0*dV
for i in range(4): dm[i,...] = V*E*dF[i] + V*dE[i]*F + dV[i,...]*E*F
return V*E*F, dm
return V*E*F
def gauss(XY,param):
"""
......
......@@ -82,17 +82,18 @@ def test_moffat_grad():
npix = 512
XY = np.mgrid[0:npix,0:npix]-npix//2
param = [35,25,1.3,1.5,0,0]
# for norm=None
m,g = moffat(XY, param, norm=None, grad=True)
gnum = g*0
eps = 1e-6
for i in range(gnum.shape[0]):
dp = copy.copy(param)
dp[i] += eps
m2 = moffat(XY, dp, norm=None)
gnum[i,...] = (m2-m)/eps
assert errL1(g,gnum)==pytest.approx(0,abs=1e-5)
# TODO: similar gradient test for norm=np.inf or any float
for norm in [None,np.inf,50]: # test the different normalization
for rm in [0,4]: # test different removeInside
m,g = moffat(XY, param, norm=norm, removeInside=rm, grad=True)
gnum = g*0
eps = 1e-7
for i in range(gnum.shape[0]):
dp = copy.copy(param)
dp[i] += eps
m2 = moffat(XY, dp, norm=norm, removeInside=rm)
gnum[i,...] = (m2-m)/eps
assert errL1(g,gnum)==pytest.approx(0,abs=1e-6)
#%% TEST PSFfromPSD and PSFAO MODELS
def test_psffrompsd_oversampling():
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment