Commit 6ce9481a authored by FETICK Romain's avatar FETICK Romain
Browse files

Add unittest

parent 0dd95031
......@@ -41,8 +41,7 @@ psf = Pmodel(param,dx=0,dy=0)
noise = ron*np.random.randn(npix,npix)
image = flux*psf + bckgd + noise
#%% Fit the PSF
#%% Fit the PSF with Psfao
guess = [0.145,2e-7,1.2,0.08,1.5]
w = np.ones_like(image)/ron**2.0
......
......@@ -15,8 +15,9 @@ from maoppy.config import _EPSILON
from maoppy.utils import binning #, getgaussnoise
#%% FITTING FUNCTION
def rmserror(model,image,weights=1.0,background=True, positive_bck=False):
def rmserror(model,image,weights=None,background=True, positive_bck=False):
"""Compute the RMS error between a PSF model and an image"""
if weights is None: weights = np.ones_like(image)
amp,bck = lsq_flux_bck(model, image, weights=weights, background=background, positive_bck=positive_bck)
diff = amp*model+bck - image
#noise_std = getgaussnoise(image)
......@@ -439,13 +440,12 @@ class Psfao(ParametricPSF):
@property
def samp(self):
return self._samp
return self._samp_over/self._k
@samp.setter
def samp(self,value):
# Manage cases of undersampling
self._samp = value
self._samp_num, self._k = oversample(value)
self._samp_over, self._k = oversample(value)
def psd(self,x0):
"""Compute the PSD model from parameters
......@@ -468,7 +468,7 @@ class Psfao(ParametricPSF):
Nx_over = self.Npix[0]*self._k
Ny_over = self.Npix[1]*self._k
pix2freq = 1.0/(self.system.D*self._samp_num)
pix2freq = 1.0/(self.system.D*self._samp_over)
f2D = np.mgrid[0:Nx_over, 0:Ny_over].astype(float)
f2D[0] -= Nx_over//2
f2D[1] -= Ny_over//2
......@@ -507,7 +507,7 @@ class Psfao(ParametricPSF):
def _otf_turbulent(self,x0):
PSD = self.psd(x0)
L = self.system.D * self._samp_num
L = self.system.D * self._samp_over
Bg = fft2(fftshift(PSD)) / L**2
Dphi = fftshift(np.real(2 * (Bg[0, 0] - Bg)))
return np.exp(-Dphi/2.)
......@@ -516,10 +516,10 @@ class Psfao(ParametricPSF):
Nx_over = self.Npix[0]*self._k
Ny_over = self.Npix[1]*self._k
NpupX = np.ceil(Nx_over/self._samp_num)
NpupY = np.ceil(Ny_over/self._samp_num)
NpupX = np.ceil(Nx_over/self._samp_over)
NpupY = np.ceil(Ny_over/self._samp_over)
tab = np.zeros((Nx_over, Ny_over), dtype=np.complex)
tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_num)
tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
return fftshift(abs(ifft2(abs(fft2(tab)) ** 2)) / np.sum(tab))
def _otf_shift(self,dx,dy):
......
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May 8 20:21:20 2020
@author: rfetick
"""
import unittest
import numpy as np
from maoppy.psfmodel import lsq_flux_bck, moffat, gauss, Psfao
from maoppy.instrument import MUSE_NFM
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,symmetric=True,samp=samp)
self.assertEqual(samp,P.samp)
self.assertGreaterEqual(P._samp_over,2.0)
self.assertEqual(P._k%1,0)
def test_symmetric(self):
npix = 100
samp = 0.3
P = Psfao((npix,npix),system=MUSE_NFM,symmetric=True,samp=samp)
self.assertEqual(len(P.bounds[0]),5)
self.assertEqual(len(P.bounds[1]),5)
P = Psfao((npix,npix),system=MUSE_NFM,symmetric=False,samp=samp)
self.assertEqual(len(P.bounds[0]),7)
self.assertEqual(len(P.bounds[1]),7)
if __name__=='__main__':
unittest.main()
\ No newline at end of file
......@@ -15,6 +15,6 @@ setup(name='MAOPPY',
author='Romain JL Fetick (LAM, Marseille, France)',
author_email='romain.fetick@lam.fr',
description='Modelization of the Adaptive Optics Psf in PYthon (MAOPPY)',
packages=find_packages(exclude=['examples','tests']),
packages=find_packages(exclude=['example','test']),
requires=['numpy','scipy','astropy'],
zip_safe=False)
Supports Markdown
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