Commit e1e66f43 authored by rfetick's avatar rfetick
Browse files

add reducedcoord grad, add pytest script

parent 838e9678
......@@ -191,16 +191,34 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
return result
#%% BASIC FUNCTIONS
def reduced_coord(XY,ax,ay,theta,cx,cy):
def reduced_coord(XY,ax,ay,theta,cx,cy,grad=False):
c = _np.cos(theta)
s = _np.sin(theta)
s2 = _np.sin(2.0 * theta)
Rxx = (c/ax)**2 + (s/ay)**2
Ryy = (c/ay)**2 + (s/ax)**2
Rxy = s2/ay**2 - s2/ax**2
Ryy = (c/ay)**2 + (s/ax)**2
uxx = (XY[0]-cx)**2
uxy = (XY[0]-cx)*(XY[1]-cy)
uyy = (XY[1]-cy)**2
u = Rxx*uxx + Rxy*uxy + Ryy*uyy
if grad:
g = _np.zeros((3,)+u.shape)
# du/dax
g[0,...] = -2/ax * (uxx*(c/ax)**2 - uxy*s2/ax**2 + uyy*(s/ax)**2)
# du/day
g[1,...] = -2/ay * (uxx*(s/ay)**2 + uxy*s2/ay**2 + uyy*(c/ay)**2)
# du/dtheta
drxx = -2*c*s/ax**2 + 2*s*c/ay**2
drxy = 2*_np.cos(2*theta)*(1/ay**2 - 1/ax**2)
dryy = -2*c*s/ay**2 + 2*s*c/ax**2
g[2,...] = drxx*uxx + drxy*uxy + dryy*uyy
return u,g
u = Rxx*(XY[0]-cx)**2 + Rxy*(XY[0]-cx)*(XY[1]-cy) + Ryy*(XY[1]-cy)**2
return u
def moffat(XY, param, norm=None, removeInside=0):
......@@ -538,7 +556,7 @@ class ParametricPSFfromPSD(ParametricPSF):
nx, ny = self._nxnyOver
NpupX = _np.ceil(nx/self._samp_over)
NpupY = _np.ceil(ny/self._samp_over)
tab = _np.zeros((nx, ny), dtype=_np.complex)
tab = _np.zeros((nx, ny), dtype=complex)
tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
return _fft.fftshift(_np.abs(_fft.ifft2(_np.abs(_fft.fft2(tab)) ** 2)) / _np.sum(tab))
......
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May 8 20:21:20 2020
Created on Tue Jun 1 16:30:27 2021
@author: rfetick
"""
import unittest
import pytest
import numpy as np
from maoppy.psfmodel import lsq_flux_bck, moffat, gauss, Psfao
from maoppy.instrument import zimpol, muse_nfm
from maoppy.psfmodel import (reduced_coord, lsq_flux_bck, moffat, gauss,
ParametricPSFfromPSD, Psfao)
from maoppy.instrument import muse_nfm, zimpol
class TestLsqFluxBck(unittest.TestCase):
def test_fluxbck(self):
bck = -5.0
flux = 8.0
model = np.zeros((10,10))
model[3,3] = 1.0
image = flux*model+bck
w = np.ones_like(image)
f,b = lsq_flux_bck(model,image,w)
self.assertAlmostEqual(flux,f,delta=1e-8)
self.assertAlmostEqual(bck,b,delta=1e-8)
def test_positivebck(self):
bck = -2.0
flux = 8.0
model = np.zeros((10,10))
model[3,3] = 1.0
image = flux*model+bck
w = np.ones_like(image)
f,b = lsq_flux_bck(model,image,w,positive_bck=True)
self.assertEqual(b,0.0)
class TestMoffat(unittest.TestCase):
def test_energy(self):
npix = 512
XY = np.mgrid[0:npix,0:npix] - npix/2.0
param = [5.0,6.0,np.pi/4,1.6,0,0]
m = moffat(XY,param,norm=np.inf)
self.assertAlmostEqual(np.sum(m),1.0,delta=1e-2)
m = moffat(XY,param,norm=None)
self.assertAlmostEqual(np.max(m),1.0,delta=1e-2)
class TestGauss(unittest.TestCase):
def test_energy(self):
npix = 512
XY = np.mgrid[0:npix,0:npix] - npix/2.0
param = [5.0,6.0,np.pi/4,0,0]
g = gauss(XY,param)
self.assertAlmostEqual(np.sum(g),1.0,delta=1e-2)
class TestPsfao(unittest.TestCase):
def test_oversampling(self):
npix = 100
samp = 0.3
P = Psfao((npix,npix),system=muse_nfm,samp=samp)
self.assertEqual(samp,P.samp)
self.assertGreaterEqual(P._samp_over,2.0)
self.assertEqual(P._k%1,0)
def test_bounds_length(self):
npix = 100
samp = 0.3
P = Psfao((npix,npix),system=muse_nfm,samp=samp)
self.assertEqual(len(P.bounds[0]),7)
self.assertEqual(len(P.bounds[1]),7)
def test_psd_integral(self):
npix = 1024
samp = 2.0
P = Psfao((npix,npix),system=zimpol,samp=samp)
fao = P.system.Nact/(2.0*P.system.D)
df = 1.0/(P.system.D*P._samp_over)
r0 = 0.15
C = 0
A = 0
alpha = 0.5
ellip = 1.0
theta = 0
beta = 1.5
# Integral of the halo
psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
int_num = np.sum(psd)*df*df
int_ana = 0.023 * 6*np.pi/5 * (fao*r0)**(-5.0/3.0)
self.assertAlmostEqual(int_num,int_ana,delta=1e-2)
# Integral of the constant
r0 = np.inf
C = 1e-2
psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
int_num = np.sum(psd)*df*df
int_ana = C*np.pi*fao**2.0 - df*df # remove central pixel
self.assertAlmostEqual(int_num,int_ana,delta=1e-2)
# Integral of the Moffat
C = 0.0
A = 1.0
psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
int_num = np.sum(psd)*df*df
int_ana = A
self.assertAlmostEqual(int_num,int_ana,delta=1e-2)
def test_otf_max(self):
npix = 1024
samp = 2.0
P = Psfao((npix,npix),system=zimpol,samp=samp)
r0 = 0.25
C = 1e-4
A = 1.0
ax = 0.05
ay = 0.08
theta = np.pi/4
beta = 1.5
otf = P.otf([r0,C,A,ax,ay,theta,beta])
mtf = np.max(np.abs(otf))
self.assertAlmostEqual(mtf,1.0,delta=1e-2)
self.assertLessEqual(mtf,1.0)
def errL1(a,b):
avg = np.sum(np.abs(a)+np.abs(b))/2
diff = np.sum(np.abs(a-b))
if avg==0: return 0 # if the norms of a and b are both null!
return diff/avg
if __name__=='__main__':
unittest.main()
\ No newline at end of file
def test_reduced_coord_grad():
ax = 3
ay = 4
theta = 1.3
npix = 256
XY = np.mgrid[0:npix,0:npix]
cx = npix//2 - 3
cy = npix//2 - 5
u,g = reduced_coord(XY, ax, ay, theta, cx, cy, grad=True)
# numerical gradient
gnum = 0*g
eps = 1e-8
u2 = reduced_coord(XY, ax+eps, ay, theta, cx, cy)
gnum[0,...] = (u2-u)/eps
u2 = reduced_coord(XY, ax, ay+eps, theta, cx, cy)
gnum[1,...] = (u2-u)/eps
u2 = reduced_coord(XY, ax, ay, theta+eps, cx, cy)
gnum[2,...] = (u2-u)/eps
assert errL1(g,gnum) == pytest.approx(0,abs=1e-7)
def test_fluxbck():
bck = -5.0
flux = 8.0
model = np.zeros((10,10))
model[3:5,3:5] = 1.0/4 # dummy PSF with unit sum
image = flux*model+bck
w = np.ones_like(image)
f,b = lsq_flux_bck(model,image,w)
assert f==pytest.approx(flux,rel=1e-8)
assert b==pytest.approx(bck,rel=1e-8)
def test_fluxbck_positive():
bck = -2.0
flux = 8.0
model = np.zeros((10,10))
model[3,3] = 1.0
image = flux*model+bck
w = np.ones_like(image)
f,b = lsq_flux_bck(model,image,w,positive_bck=True)
assert b==0
def test_moffat_energy():
npix = 512
XY = np.mgrid[0:npix,0:npix] - npix/2.0
param = [5.0,6.0,np.pi/4,1.6,0,0]
m = moffat(XY,param,norm=np.inf)
assert np.sum(m)==pytest.approx(1.0,abs=1e-2)
m = moffat(XY,param,norm=None)
assert np.max(m)==pytest.approx(1.0,abs=1e-8)
def test_gauss_energy():
npix = 512
XY = np.mgrid[0:npix,0:npix] - npix/2.0
param = [5.0,6.0,np.pi/4,0,0]
g = gauss(XY,param)
assert np.sum(g)==pytest.approx(1.0,abs=1e-2)
def test_psffrompsd_oversampling():
npix = 100
samp = 0.3
nparam = 8 # whatever
P = ParametricPSFfromPSD(nparam,(npix,npix),system=muse_nfm,samp=samp)
assert samp==P.samp
assert P._samp_over>=2.0
assert (P._k%1)==0
def test_psfao_bounds():
npix = 100
samp = 0.3
P = Psfao((npix,npix),system=muse_nfm,samp=samp)
assert len(P.bounds[0])==7
assert len(P.bounds[1])==7
def test_psfao_psd_integral():
npix = 1024
samp = 2.0
P = Psfao((npix,npix),system=zimpol,samp=samp)
fao = P.system.Nact/(2.0*P.system.D)
df = 1.0/(P.system.D*P._samp_over)
r0 = 0.15
C = 0
A = 0
alpha = 0.5
ellip = 1.0
theta = 0
beta = 1.5
# Integral of the halo
psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
int_num = np.sum(psd)*df*df
int_ana = 0.0229 * 6*np.pi/5 * (fao*r0)**(-5.0/3.0)
assert int_num==pytest.approx(int_ana,abs=1e-2)
# Integral of the constant
r0 = np.inf
C = 1e-2
psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
int_num = np.sum(psd)*df*df
int_ana = C*np.pi*fao**2.0 - df*df # remove central pixel
assert int_num==pytest.approx(int_ana,abs=1e-2)
# Integral of the Moffat
C = 0.0
A = 1.0
psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
int_num = np.sum(psd)*df*df
int_ana = A
assert int_num==pytest.approx(int_ana,abs=1e-2)
def test_otf_max():
npix = 1024
samp = 2.0
P = Psfao((npix,npix),system=zimpol,samp=samp)
r0 = 0.25
C = 1e-4
A = 1.0
ax = 0.05
ay = 0.08
theta = np.pi/4
beta = 1.5
otf = P.otf([r0,C,A,ax,ay,theta,beta])
mtf = np.max(np.abs(otf))
assert mtf==pytest.approx(1.0,abs=1e-2)
assert mtf<=1.0
if __name__=="__main__":
pytest.main()
\ No newline at end of file
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