diff --git a/eastlake/__init__.py b/eastlake/__init__.py index 9894e88..24257f5 100644 --- a/eastlake/__init__.py +++ b/eastlake/__init__.py @@ -11,3 +11,4 @@ from .pipeline import register_pipeline_step # noqa from . import des_piff # noqa +from . import des_smoothpiff # noqa diff --git a/eastlake/des_piff.py b/eastlake/des_piff.py index 8ac1338..97e150b 100644 --- a/eastlake/des_piff.py +++ b/eastlake/des_piff.py @@ -1,29 +1,14 @@ import os import logging +import functools import galsim import galsim.config import piff -import numpy as np -import ngmix -if ngmix.__version__[0:2] == "v1": - NGMIX_V2 = False - from ngmix.fitting import LMSimple - from ngmix.admom import Admom -else: - NGMIX_V2 = True - from ngmix.fitting import Fitter - from ngmix.admom import AdmomFitter - -from scipy.interpolate import CloughTocher2DInterpolator - logger = logging.getLogger(__name__) -# pixel scale used for fitting the Piff models -PIFF_SCALE = 0.25 - PSF_KWARGS = { "g": {"GI_COLOR": 1.1}, @@ -33,166 +18,36 @@ } +@functools.lru_cache(maxsize=200) +def _read_piff(file_name): + return piff.read( + os.path.expanduser(os.path.expandvars(file_name)) + ) + + class DES_Piff(object): """A wrapper for Piff to use with Galsim. - This wrapper uses ngmix to fit smooth models to the Piff PSF images. The - parameters of these models are then interpolated across the SE image - and used to generate a smooth approximation to the PSF. - Parameters ---------- file_name : str The file with the Piff psf solution. - smooth : bool, optional - If True, then smooth the Piff PSFs. Default of False. """ _req_params = {'file_name': str} _opt_params = {} _single_params = [] _takes_rng = False - def __init__(self, file_name, smooth=False, psf_kwargs=None): + def __init__(self, file_name): self.file_name = file_name - self.psf_kwargs = psf_kwargs or {} - # Read the Piff file. This may fail if the Piff - # file is missing. We catch this and continue - # since if we're substituting in some different - # PSF model for rejectlisted piff files, we'll - # never actually use self._piff - try: - self._piff = piff.read( - os.path.expanduser(os.path.expandvars(file_name))) - - self._chipnum = list(set(self._piff.wcs.keys()))[0] - assert len(list(set(self._piff.wcs.keys()))) == 1 - except IOError as e: - print("failed to load Piff file, hopefully it's rejectlisted: %s" % repr(e)) - self._piff = None - self._did_fit = False - self.smooth = smooth - - def _fit_smooth_model(self): - dxy = 256 - ny = 4096 // dxy + 1 - nx = 2048 // dxy + 1 - - xloc = np.empty((ny, nx), dtype=np.float64) - yloc = np.empty((ny, nx), dtype=np.float64) - pars = np.empty((ny, nx, 3), dtype=np.float64) - for yi, yl in enumerate(np.linspace(1, 4096, ny)): - for xi, xl in enumerate(np.linspace(1, 2048, nx)): - rng = np.random.RandomState(seed=yi + nx * xi) - xloc[yi, xi] = xl - yloc[yi, xi] = yl - - pos = galsim.PositionD(x=xl, y=yl) - gs_img = self._draw(pos).drawImage( - nx=19, ny=19, scale=PIFF_SCALE, method='sb') - img = gs_img.array - nse = np.std( - np.concatenate([img[0, :], img[-1, :]])) - obs = ngmix.Observation( - image=img, - weight=np.ones_like(img)/nse**2, - jacobian=ngmix.jacobian.DiagonalJacobian( - x=9, y=9, scale=PIFF_SCALE)) - - _g1 = np.nan - _g2 = np.nan - _T = np.nan - - # there are some nutty PSFs - if gs_img.calculateFWHM() > 0.5: - for _ in range(5): - try: - if NGMIX_V2: - am = AdmomFitter(rng=rng) - res = am.go(obs, 0.3) - if res['flags'] != 0: - continue - - lm = Fitter(model='turb') - lm_res = lm.go(obs, res['pars']) - if lm_res['flags'] == 0: - _g1 = lm_res['pars'][2] - _g2 = lm_res['pars'][3] - _T = lm_res['pars'][4] - break - else: - am = Admom(obs, rng=rng) - am.go(0.3) - res = am.get_result() - if res['flags'] != 0: - continue - - lm = LMSimple(obs, 'turb') - lm.go(res['pars']) - lm_res = lm.get_result() - if lm_res['flags'] == 0: - _g1 = lm_res['pars'][2] - _g2 = lm_res['pars'][3] - _T = lm_res['pars'][4] - break - except ngmix.gexceptions.GMixRangeError: - pass - - try: - irr, irc, icc = ngmix.moments.g2mom(_g1, _g2, _T) - # this is a fudge factor that gets the overall PSF FWHM - # correct - # the naive correction for the pixel size is - # a bit too small - pixel_var = PIFF_SCALE * PIFF_SCALE / 12 * 1.73 - irr -= pixel_var - icc -= pixel_var - _g1, _g2, _T = ngmix.moments.mom2g(irr, irc, icc) - except Exception: - _g1 = np.nan - _g2 = np.nan - _T = np.nan - - pars[yi, xi, 0] = _g1 - pars[yi, xi, 1] = _g2 - pars[yi, xi, 2] = _T - - xloc = xloc.ravel() - yloc = yloc.ravel() - pos = np.stack([xloc, yloc], axis=1) - assert pos.shape == (xloc.shape[0], 2) - - # make interps - g1 = pars[:, :, 0].ravel() - msk = np.isfinite(g1) - if len(msk) < 10: - raise ValueError('DES Piff fitting failed too much!') - if np.any(~msk): - g1[~msk] = np.mean(g1[msk]) - self._g1int = CloughTocher2DInterpolator( - pos, g1, fill_value=np.mean(g1[msk])) - - g2 = pars[:, :, 1].ravel() - msk = np.isfinite(g2) - if len(msk) < 10: - raise ValueError('DES Piff fitting failed too much!') - if np.any(~msk): - g2[~msk] = np.mean(g2[msk]) - self._g2int = CloughTocher2DInterpolator( - pos, g2, fill_value=np.mean(g2[msk])) - - T = pars[:, :, 2].ravel() - msk = np.isfinite(T) - if len(msk) < 10: - raise ValueError('DES Piff fitting failed too much!') - if np.any(~msk): - T[~msk] = np.mean(T[msk]) - self._Tint = CloughTocher2DInterpolator( - pos, T, fill_value=np.mean(T[msk])) - - self._did_fit = True - - def _draw(self, image_pos, wcs=None, n_pix=None, - x_interpolant='lanczos15', gsparams=None): + self._piff = _read_piff(file_name) + self._chipnum = list(set(self._piff.wcs.keys()))[0] + assert len(list(set(self._piff.wcs.keys()))) == 1 + + def _draw( + self, image_pos, wcs=None, n_pix=None, + gsparams=None, **psf_kwargs + ): """Get an image of the PSF at the given location. Parameters @@ -204,55 +59,126 @@ def _draw(self, image_pos, wcs=None, n_pix=None, n_pix : int, optional The image size to use when drawing without smoothing. Defaults to 53 pixels if not given - x_interpolant : str, optional - The interpolant to use. gsparams : galsim.GSParams, optional - Ootional galsim configuration data to pass along. + Optional galsim configuration data to pass along. + psf_kwargs : extra keyword arguments + These are passed to the Piff model as needed. Returns ------- - psf : galsim.InterpolatedImage - The PSF at the image position. + psf_img : galsim.ImageD + The Galsim image of the PSF. + pixel_wcs : galsim.BaseWCS + The WCS of the image + offset : tuple of floats + The offset from the true center of the PSF centroid. """ + if n_pix is None: + n_pix = 53 + if wcs is not None: - if n_pix is not None: - n_pix = n_pix - else: - n_pix = 53 pixel_wcs = wcs.local(image_pos) else: - n_pix = 19 - pixel_wcs = galsim.PixelScale(PIFF_SCALE) + pixel_wcs = galsim.PixelScale(0.263) # nice and big image size here cause this has been a problem image = galsim.ImageD(ncol=n_pix, nrow=n_pix, wcs=pixel_wcs) - psf = self.getPiff().draw( + if "GI_COLOR" in self.getPiff().interp_property_names: + psf_kwargs.pop("IZ_COLOR", None) + if ( + "GI_COLOR" in psf_kwargs and ( + psf_kwargs["GI_COLOR"] is None or + psf_kwargs["GI_COLOR"] == "None" + ) + ): + psf_kwargs["GI_COLOR"] = 1.1 + + elif "IZ_COLOR" in self.getPiff().interp_property_names: + psf_kwargs.pop("GI_COLOR", None) + + if ( + "IZ_COLOR" in psf_kwargs and ( + psf_kwargs["IZ_COLOR"] is None or + psf_kwargs["IZ_COLOR"] == "None" + ) + ): + psf_kwargs["IZ_COLOR"] = 0.34 + + offset = ( + image_pos.x - int(image_pos.x + 0.5), + image_pos.y - int(image_pos.y + 0.5), + ) + + psf_img = self.getPiff().draw( image_pos.x, image_pos.y, image=image, - center=True, chipnum=self._chipnum, - **self.psf_kwargs, + center=True, + offset=offset, + **psf_kwargs, + ) + + return psf_img, pixel_wcs, offset + + def getPiff(self): + return self._piff + + def getPSF( + self, image_pos, wcs=None, + n_pix=None, depixelize=False, + gsparams=None, + **kwargs + ): + """Get an image of the PSF at the given location. + + Parameters + ---------- + image_pos : galsim.Position + The image position for the PSF. + wcs : galsim.BaseWCS or subclass, optional + The WCS to use to draw the PSF. If not given, the WCS in the Piff model is used. + n_pix : int, optional + The image size to use when drawing without smoothing. + depixelize : bool, optional + If True, the interpolated image will be depixelized. Default is False. + gsparams : galsim.GSParams, optional + Optional galsim configuration data to pass along. + **kwargs : extra keyword arguments + These are all passed on to the Piff model when drawing. + + Returns + ------- + psf : galsim.GSObject + The PSF at the image position. + """ + + psf_img, pixel_wcs, offset = self._draw( + image_pos, + wcs=wcs, + n_pix=n_pix, + gsparams=gsparams, + **kwargs ) psf = galsim.InterpolatedImage( - galsim.ImageD(psf.array), # make sure galsim is not keeping state + galsim.ImageD(psf_img.array), # make sure galsim is not keeping state wcs=pixel_wcs, gsparams=gsparams, - x_interpolant=x_interpolant + x_interpolant='lanczos15', + depixelize=depixelize, + offset=offset, ).withFlux( 1.0 ) - return psf - def getPiff(self): - return self._piff - - def getPSF( + def getPSFImage( self, image_pos, wcs=None, - smooth=False, n_pix=None, **kwargs + n_pix=None, + gsparams=None, + **kwargs ): """Get an image of the PSF at the given location. @@ -261,39 +187,27 @@ def getPSF( image_pos : galsim.Position The image position for the PSF. wcs : galsim.BaseWCS or subclass, optional - The WCS to use to draw the PSF. Currently used only when smoothing - is turned off. - smooth : bool, optional - If True, then smooth the Piff PSFs. Default of False. + The WCS to use to draw the PSF. If not given, the WCS in the Piff model is used. n_pix : int, optional The image size to use when drawing without smoothing. - **kargs : extra keyword arguments - These are all ignored. + gsparams : galsim.GSParams, optional + Optional galsim configuration data to pass along. + **kwargs : extra keyword arguments + These are all passed on to the Piff model when drawing. Returns ------- - psf : galsim.GSObject - The PSF at the image position. + psf : galsim.Image + The PSF image at the image position. """ - if smooth or self.smooth: - if not self._did_fit: - self._fit_smooth_model() - - arr = np.array([ - np.clip(image_pos.x, 1, 2048), - np.clip(image_pos.y, 1, 4096)]) - - _g1 = self._g1int(arr)[0] - _g2 = self._g2int(arr)[0] - _T = self._Tint(arr)[0] - if np.any(np.isnan(np.array([_g1, _g2, _T]))): - logger.debug("Piff smooth fit params are NaN: %s %s %s %s", image_pos, _g1, _g2, _T) - raise RuntimeError("NaN smooth Piff params at %s!" % image_pos) - pars = np.array([0, 0, _g1, _g2, _T, 1]) - obj = ngmix.gmix.make_gmix_model(pars, 'turb').make_galsim_object() - return obj.withFlux(1) - else: - return self._draw(image_pos, wcs=wcs, n_pix=n_pix) + psf_img, pixel_wcs, offset = self._draw( + image_pos, + wcs=wcs, + n_pix=n_pix, + gsparams=gsparams, + **kwargs + ) + return psf_img class PiffLoader(galsim.config.InputLoader): @@ -301,7 +215,8 @@ def getKwargs(self, config, base, logger): req = {'file_name': str} opt = {} kwargs, safe = galsim.config.GetAllParams( - config, base, req=req, opt=opt) + config, base, req=req, opt=opt + ) return kwargs, safe @@ -312,15 +227,23 @@ def getKwargs(self, config, base, logger): # and a builder def BuildDES_Piff(config, base, ignore, gsparams, logger): - des_piff = galsim.config.GetInputObj('des_piff', config, base, 'DES_Piff') - - opt = {'flux': float, - 'num': int, - 'image_pos': galsim.PositionD, - 'x_interpolant': str, - 'smooth': bool} + opt = { + 'flux': float, + 'image_pos': galsim.PositionD, + 'n_pix': int, + 'gi_color': float, + 'iz_color': float, + 'depixelize': bool, + 'file_name': str, + } params, safe = galsim.config.GetAllParams( - config, base, opt=opt, ignore=ignore) + config, base, opt=opt, ignore=ignore + ) + + if params.get("file_name", None) is None: + des_piff = galsim.config.GetInputObj('des_piff', config, base, 'DES_Piff') + else: + des_piff = DES_Piff(params.get("file_name", None)) if 'image_pos' in params: image_pos = params['image_pos'] @@ -342,9 +265,13 @@ def BuildDES_Piff(config, base, ignore, gsparams, logger): psf = des_piff.getPSF( image_pos, - wcs, - smooth=params.get('smooth', False), - gsparams=gsparams) + wcs=wcs, + n_pix=params.get("n_pix", None), + GI_COLOR=params.get("gi_color", None), + IZ_COLOR=params.get("iz_color", None), + depixelize=params.get("depixelize", False), + gsparams=gsparams, + ) if 'flux' in params: psf = psf.withFlux(params['flux']) @@ -354,24 +281,5 @@ def BuildDES_Piff(config, base, ignore, gsparams, logger): return psf, can_be_reused -def BuildDES_Piff_with_substitute(config, base, ignore, gsparams, logger): - # This builder usually just calls BuildDES_Piff, but can also - # be passed use_substitute = True, in which case it builds some - # other PSF. We use this for rejectlisted Piff files. - if "use_substitute" in config: - use_substitute = galsim.config.ParseValue(config, "use_substitute", - base, bool)[0] - else: - use_substitute = False - - if use_substitute: - return (galsim.config.BuildGSObject( - config, "substitute_psf", base=base, - gsparams=gsparams, logger=logger)) - else: - ignore += ["use_substitute", "substitute_psf"] - return BuildDES_Piff(config, base, ignore, gsparams, logger) - - galsim.config.RegisterObjectType( - 'DES_Piff', BuildDES_Piff_with_substitute, input_type='des_piff') + 'DES_Piff', BuildDES_Piff, input_type='des_piff') diff --git a/eastlake/des_smoothpiff.py b/eastlake/des_smoothpiff.py new file mode 100644 index 0000000..c3c4d5f --- /dev/null +++ b/eastlake/des_smoothpiff.py @@ -0,0 +1,369 @@ +import os +import logging + +import galsim +import galsim.config + +import piff + +import numpy as np +import ngmix +if ngmix.__version__[0:2] == "v1": + NGMIX_V2 = False + from ngmix.fitting import LMSimple + from ngmix.admom import Admom +else: + NGMIX_V2 = True + from ngmix.fitting import Fitter + from ngmix.admom import AdmomFitter + +from scipy.interpolate import CloughTocher2DInterpolator + +logger = logging.getLogger(__name__) + +# pixel scale used for fitting the Piff models +PIFF_SCALE = 0.25 + + +class DES_SmoothPiff(object): + """A wrapper for Piff to use with Galsim. + + This wrapper uses ngmix to fit smooth models to the Piff PSF images. The + parameters of these models are then interpolated across the SE image + and used to generate a smooth approximation to the PSF. + + Parameters + ---------- + file_name : str + The file with the Piff psf solution. + smooth : bool, optional + If True, then smooth the Piff PSFs. Default of False. + """ + _req_params = {'file_name': str} + _opt_params = {} + _single_params = [] + _takes_rng = False + + def __init__(self, file_name, smooth=False, psf_kwargs=None): + self.file_name = file_name + self.psf_kwargs = psf_kwargs or {} + # Read the Piff file. This may fail if the Piff + # file is missing. We catch this and continue + # since if we're substituting in some different + # PSF model for rejectlisted piff files, we'll + # never actually use self._piff + try: + self._piff = piff.read( + os.path.expanduser(os.path.expandvars(file_name))) + + self._chipnum = list(set(self._piff.wcs.keys()))[0] + assert len(list(set(self._piff.wcs.keys()))) == 1 + except IOError as e: + print("failed to load Piff file, hopefully it's rejectlisted: %s" % repr(e)) + self._piff = None + self._did_fit = False + self.smooth = smooth + + def _fit_smooth_model(self): + dxy = 256 + ny = 4096 // dxy + 1 + nx = 2048 // dxy + 1 + + xloc = np.empty((ny, nx), dtype=np.float64) + yloc = np.empty((ny, nx), dtype=np.float64) + pars = np.empty((ny, nx, 3), dtype=np.float64) + for yi, yl in enumerate(np.linspace(1, 4096, ny)): + for xi, xl in enumerate(np.linspace(1, 2048, nx)): + rng = np.random.RandomState(seed=yi + nx * xi) + xloc[yi, xi] = xl + yloc[yi, xi] = yl + + pos = galsim.PositionD(x=xl, y=yl) + gs_img = self._draw(pos).drawImage( + nx=19, ny=19, scale=PIFF_SCALE, method='sb') + img = gs_img.array + nse = np.std( + np.concatenate([img[0, :], img[-1, :]])) + obs = ngmix.Observation( + image=img, + weight=np.ones_like(img)/nse**2, + jacobian=ngmix.jacobian.DiagonalJacobian( + x=9, y=9, scale=PIFF_SCALE)) + + _g1 = np.nan + _g2 = np.nan + _T = np.nan + + # there are some nutty PSFs + if gs_img.calculateFWHM() > 0.5: + for _ in range(5): + try: + if NGMIX_V2: + am = AdmomFitter(rng=rng) + res = am.go(obs, 0.3) + if res['flags'] != 0: + continue + + lm = Fitter(model='turb') + lm_res = lm.go(obs, res['pars']) + if lm_res['flags'] == 0: + _g1 = lm_res['pars'][2] + _g2 = lm_res['pars'][3] + _T = lm_res['pars'][4] + break + else: + am = Admom(obs, rng=rng) + am.go(0.3) + res = am.get_result() + if res['flags'] != 0: + continue + + lm = LMSimple(obs, 'turb') + lm.go(res['pars']) + lm_res = lm.get_result() + if lm_res['flags'] == 0: + _g1 = lm_res['pars'][2] + _g2 = lm_res['pars'][3] + _T = lm_res['pars'][4] + break + except ngmix.gexceptions.GMixRangeError: + pass + + try: + irr, irc, icc = ngmix.moments.g2mom(_g1, _g2, _T) + # this is a fudge factor that gets the overall PSF FWHM + # correct + # the naive correction for the pixel size is + # a bit too small + pixel_var = PIFF_SCALE * PIFF_SCALE / 12 * 1.73 + irr -= pixel_var + icc -= pixel_var + _g1, _g2, _T = ngmix.moments.mom2g(irr, irc, icc) + except Exception: + _g1 = np.nan + _g2 = np.nan + _T = np.nan + + pars[yi, xi, 0] = _g1 + pars[yi, xi, 1] = _g2 + pars[yi, xi, 2] = _T + + xloc = xloc.ravel() + yloc = yloc.ravel() + pos = np.stack([xloc, yloc], axis=1) + assert pos.shape == (xloc.shape[0], 2) + + # make interps + g1 = pars[:, :, 0].ravel() + msk = np.isfinite(g1) + if len(msk) < 10: + raise ValueError('DES Piff fitting failed too much!') + if np.any(~msk): + g1[~msk] = np.mean(g1[msk]) + self._g1int = CloughTocher2DInterpolator( + pos, g1, fill_value=np.mean(g1[msk])) + + g2 = pars[:, :, 1].ravel() + msk = np.isfinite(g2) + if len(msk) < 10: + raise ValueError('DES Piff fitting failed too much!') + if np.any(~msk): + g2[~msk] = np.mean(g2[msk]) + self._g2int = CloughTocher2DInterpolator( + pos, g2, fill_value=np.mean(g2[msk])) + + T = pars[:, :, 2].ravel() + msk = np.isfinite(T) + if len(msk) < 10: + raise ValueError('DES Piff fitting failed too much!') + if np.any(~msk): + T[~msk] = np.mean(T[msk]) + self._Tint = CloughTocher2DInterpolator( + pos, T, fill_value=np.mean(T[msk])) + + self._did_fit = True + + def _draw(self, image_pos, wcs=None, n_pix=None, + x_interpolant='lanczos15', gsparams=None): + """Get an image of the PSF at the given location. + + Parameters + ---------- + image_pos : galsim.Position + The image position for the PSF. + wcs : galsim.BaseWCS or subclass, optional + The WCS to use to draw the PSF. + n_pix : int, optional + The image size to use when drawing without smoothing. Defaults to + 53 pixels if not given + x_interpolant : str, optional + The interpolant to use. + gsparams : galsim.GSParams, optional + Ootional galsim configuration data to pass along. + + Returns + ------- + psf : galsim.InterpolatedImage + The PSF at the image position. + """ + if wcs is not None: + if n_pix is not None: + n_pix = n_pix + else: + n_pix = 53 + pixel_wcs = wcs.local(image_pos) + else: + n_pix = 19 + pixel_wcs = galsim.PixelScale(PIFF_SCALE) + + # nice and big image size here cause this has been a problem + image = galsim.ImageD(ncol=n_pix, nrow=n_pix, wcs=pixel_wcs) + + psf = self.getPiff().draw( + image_pos.x, + image_pos.y, + image=image, + center=True, + chipnum=self._chipnum, + **self.psf_kwargs, + ) + + psf = galsim.InterpolatedImage( + galsim.ImageD(psf.array), # make sure galsim is not keeping state + wcs=pixel_wcs, + gsparams=gsparams, + x_interpolant=x_interpolant + ).withFlux( + 1.0 + ) + + return psf + + def getPiff(self): + return self._piff + + def getPSF( + self, image_pos, wcs=None, + smooth=False, n_pix=None, **kwargs + ): + """Get an image of the PSF at the given location. + + Parameters + ---------- + image_pos : galsim.Position + The image position for the PSF. + wcs : galsim.BaseWCS or subclass, optional + The WCS to use to draw the PSF. Currently used only when smoothing + is turned off. + smooth : bool, optional + If True, then smooth the Piff PSFs. Default of False. + n_pix : int, optional + The image size to use when drawing without smoothing. + **kargs : extra keyword arguments + These are all ignored. + + Returns + ------- + psf : galsim.GSObject + The PSF at the image position. + """ + if smooth or self.smooth: + if not self._did_fit: + self._fit_smooth_model() + + arr = np.array([ + np.clip(image_pos.x, 1, 2048), + np.clip(image_pos.y, 1, 4096)]) + + _g1 = self._g1int(arr)[0] + _g2 = self._g2int(arr)[0] + _T = self._Tint(arr)[0] + if np.any(np.isnan(np.array([_g1, _g2, _T]))): + logger.debug("Piff smooth fit params are NaN: %s %s %s %s", image_pos, _g1, _g2, _T) + raise RuntimeError("NaN smooth Piff params at %s!" % image_pos) + pars = np.array([0, 0, _g1, _g2, _T, 1]) + obj = ngmix.gmix.make_gmix_model(pars, 'turb').make_galsim_object() + return obj.withFlux(1) + else: + return self._draw(image_pos, wcs=wcs, n_pix=n_pix) + + +class SmoothPiffLoader(galsim.config.InputLoader): + def getKwargs(self, config, base, logger): + req = {'file_name': str} + opt = {} + kwargs, safe = galsim.config.GetAllParams( + config, base, req=req, opt=opt) + + return kwargs, safe + + +# add a config input section +galsim.config.RegisterInputType('des_smoothpiff', SmoothPiffLoader(DES_SmoothPiff)) + + +# and a builder +def BuildDES_SmoothPiff(config, base, ignore, gsparams, logger): + des_piff = galsim.config.GetInputObj('des_smoothpiff', config, base, 'DES_SmoothPiff') + + opt = {'flux': float, + 'num': int, + 'image_pos': galsim.PositionD, + 'x_interpolant': str, + 'smooth': bool} + params, safe = galsim.config.GetAllParams( + config, base, opt=opt, ignore=ignore) + + if 'image_pos' in params: + image_pos = params['image_pos'] + elif 'image_pos' in base: + image_pos = base['image_pos'] + else: + raise galsim.GalSimConfigError( + "DES_SmoothPiff requested, but no image_pos defined in base.") + + if 'wcs' not in base: + raise galsim.GalSimConfigError( + "DES_SmoothPiff requested, but no wcs defined in base.") + wcs = base['wcs'] + + if gsparams: + gsparams = galsim.GSParams(**gsparams) + else: + gsparams = None + + psf = des_piff.getPSF( + image_pos, + wcs, + smooth=params.get('smooth', False), + gsparams=gsparams) + + if 'flux' in params: + psf = psf.withFlux(params['flux']) + + # we make sure to declare the returned object as not safe for reuse + can_be_reused = False + return psf, can_be_reused + + +def BuildDES_SmoothPiff_with_substitute(config, base, ignore, gsparams, logger): + # This builder usually just calls BuildDES_SmoothPiff, but can also + # be passed use_substitute = True, in which case it builds some + # other PSF. We use this for rejectlisted Piff files. + if "use_substitute" in config: + use_substitute = galsim.config.ParseValue(config, "use_substitute", + base, bool)[0] + else: + use_substitute = False + + if use_substitute: + return (galsim.config.BuildGSObject( + config, "substitute_psf", base=base, + gsparams=gsparams, logger=logger)) + else: + ignore += ["use_substitute", "substitute_psf"] + return BuildDES_SmoothPiff(config, base, ignore, gsparams, logger) + + +galsim.config.RegisterObjectType( + 'DES_SmoothPiff', BuildDES_SmoothPiff_with_substitute, input_type='des_smoothpiff') diff --git a/eastlake/pipeline.py b/eastlake/pipeline.py index bfe563c..7d88c95 100644 --- a/eastlake/pipeline.py +++ b/eastlake/pipeline.py @@ -116,8 +116,6 @@ def __init__( if config is not None: galsim.config.process.ImportModules(config) - with open(os.path.join(self.base_dir, "config.yaml"), 'w') as f: - yaml.dump(config, f) self.record_file = record_file if record_file is None: @@ -193,6 +191,16 @@ def from_config_file(cls, config_file, base_dir, logger=None, verbosity=1, config = config[0] + # make sure base_dir is an absolute path + if base_dir is not None: + base_dir = os.path.abspath(base_dir) + safe_mkdir(base_dir) + else: + base_dir = "." + + with open(os.path.join(base_dir, "orig-config.yaml"), 'w') as f: + yaml.dump(config, f) + galsim.config.process.ImportModules(config) # Process templates. @@ -233,16 +241,12 @@ def from_config_file(cls, config_file, base_dir, logger=None, verbosity=1, step_names = config["pipeline"]["steps"] steps = [] - # make sure base_dir is an absolute path - if base_dir is not None: - base_dir = os.path.abspath(base_dir) - for step_name in step_names: if step_name == "galsim": steps.append(GalSimRunner(config, base_dir, logger=logger)) else: try: - step_config = config.pop(step_name) + step_config = config[step_name] except KeyError: logger.error( "no entry for step %s found in config file, continuing with empty step_config" % step_name) @@ -269,6 +273,9 @@ def from_config_file(cls, config_file, base_dir, logger=None, verbosity=1, steps.append(step_class(step_config, base_dir, logger=logger, verbosity=step_verbosity, name=step_name)) + with open(os.path.join(base_dir, "config.yaml"), 'w') as f: + yaml.dump(config, f) + return cls( steps, base_dir, logger=logger, verbosity=verbosity, diff --git a/eastlake/tests/test_des_piff.py b/eastlake/tests/test_des_piff.py index 43fc597..00765db 100644 --- a/eastlake/tests/test_des_piff.py +++ b/eastlake/tests/test_des_piff.py @@ -6,7 +6,8 @@ import pytest -from ..des_piff import DES_Piff, PSF_KWARGS +from ..des_piff import PSF_KWARGS, DES_Piff +from ..des_smoothpiff import DES_SmoothPiff def _get_piff_file(): @@ -21,6 +22,119 @@ def _get_piff_file(): ) +@pytest.mark.skipif( + os.environ.get('TEST_DESDATA', None) is None, + reason=( + 'DES_Piff can only be tested if ' + 'test data is at TEST_DESDATA')) +@pytest.mark.parametrize("x_offset", [-0.3, 0.0, 0.3]) +@pytest.mark.parametrize("y_offset", [-0.3, 0.0, 0.3]) +def test_des_piff_centering(x_offset, y_offset): + piff_fname = _get_piff_file() + piff = DES_Piff(piff_fname) + atol = 0.01 + + # test the image it makes + psf_im = piff.getPSFImage( + galsim.PositionD(20 + x_offset, 10 + y_offset), + **PSF_KWARGS["g"], + ) + + cen = (53-1)/2 + 1 + admom = psf_im.FindAdaptiveMom() + + assert np.allclose( + admom.moments_centroid.x, + cen + x_offset, + rtol=0, + atol=atol, + ) + + assert np.allclose( + admom.moments_centroid.y, + cen + y_offset, + rtol=0, + atol=atol, + ) + + # test how well it recenters + psf_im = piff.getPSF( + galsim.PositionD(20 + x_offset, 10 + y_offset), + **PSF_KWARGS["g"], + ).drawImage(nx=53, ny=53, scale=0.263, offset=(x_offset, y_offset)) + cen = (53-1)/2 + 1 + admom = psf_im.FindAdaptiveMom() + + assert np.allclose( + admom.moments_centroid.x, + cen + x_offset, + rtol=0, + atol=atol, + ) + + assert np.allclose( + admom.moments_centroid.y, + cen + y_offset, + rtol=0, + atol=atol, + ) + + +@pytest.mark.skipif( + os.environ.get('TEST_DESDATA', None) is None, + reason=( + 'DES_Piff can only be tested if ' + 'test data is at TEST_DESDATA')) +def test_des_piff_color(): + piff_fname = _get_piff_file() + piff = DES_Piff(piff_fname) + psf_im1 = piff.getPSF( + galsim.PositionD(20, 10), + GI_COLOR=1.3, + IZ_COLOR=0.4, + ).drawImage(nx=53, ny=53, scale=0.263).array + + psf_im2 = piff.getPSF( + galsim.PositionD(20, 10), + GI_COLOR=0.7, + IZ_COLOR=0.1, + ).drawImage(nx=53, ny=53, scale=0.263).array + + assert not np.allclose(psf_im1, psf_im2) + + psf_im2 = piff.getPSF( + galsim.PositionD(20, 10), + GI_COLOR=1.3, + IZ_COLOR=0.4, + ).drawImage(nx=53, ny=53, scale=0.263).array + + assert np.array_equal(psf_im1, psf_im2) + + +@pytest.mark.skipif( + os.environ.get('TEST_DESDATA', None) is None, + reason=( + 'DES_Piff can only be tested if ' + 'test data is at TEST_DESDATA')) +def test_des_piff_raises(): + piff_fname = _get_piff_file() + piff = DES_Piff(piff_fname) + with pytest.raises(Exception) as e: + piff.getPSF( + galsim.PositionD(20, 10), + ).drawImage(nx=53, ny=53, scale=0.263).array + + assert "_COLOR" in str(e.value) + + with pytest.raises(Exception) as e: + piff.getPSF( + galsim.PositionD(20, 10), + IZ_COLOR=0.4, + ).drawImage(nx=53, ny=53, scale=0.263).array + + assert "_COLOR" in str(e.value) + + @pytest.mark.skipif( os.environ.get('TEST_DESDATA', None) is None, reason=( @@ -28,7 +142,36 @@ def _get_piff_file(): 'test data is at TEST_DESDATA')) def test_des_piff_smoke(): piff_fname = _get_piff_file() - piff = DES_Piff(piff_fname, psf_kwargs=PSF_KWARGS["g"]) + piff = DES_Piff(piff_fname) + psf_im = piff.getPSF( + galsim.PositionD(20, 10), + **PSF_KWARGS["g"], + ).drawImage(nx=53, ny=53, scale=0.263).array + + if False: + import matplotlib.pyplot as plt + plt.figure() + plt.imshow(psf_im) + import pdb + pdb.set_trace() + + assert np.all(np.isfinite(psf_im)) + assert np.allclose(np.sum(psf_im), 1, rtol=1e-4, atol=0) + + y, x = np.unravel_index(np.argmax(psf_im), psf_im.shape) + cen = (psf_im.shape[0]-1)/2 + assert y == cen + assert x == cen + + +@pytest.mark.skipif( + os.environ.get('TEST_DESDATA', None) is None, + reason=( + 'DES_SmoothPiff can only be tested if ' + 'test data is at TEST_DESDATA')) +def test_des_smoothpiff_smoke(): + piff_fname = _get_piff_file() + piff = DES_SmoothPiff(piff_fname, psf_kwargs=PSF_KWARGS["g"]) psf_im = piff.getPSF(galsim.PositionD(20, 10)).drawImage(nx=53, ny=53, scale=0.263).array if False: @@ -50,22 +193,22 @@ def test_des_piff_smoke(): @pytest.mark.skipif( os.environ.get('TEST_DESDATA', None) is None, reason=( - 'DES_Piff can only be tested if ' + 'DES_SmoothPiff can only be tested if ' 'test data is at TEST_DESDATA')) -def test_des_piff_smooth(): +def test_des_smoothpiff_smooth(): piff_fname = _get_piff_file() - piff = DES_Piff(piff_fname, psf_kwargs=PSF_KWARGS["g"]) + piff = DES_SmoothPiff(piff_fname, psf_kwargs=PSF_KWARGS["g"]) psf_im = piff.getPSF( galsim.PositionD(20, 10), ).drawImage(nx=53, ny=53, scale=0.263).array - piff = DES_Piff(piff_fname, psf_kwargs=PSF_KWARGS["g"]) + piff = DES_SmoothPiff(piff_fname, psf_kwargs=PSF_KWARGS["g"]) psf_im1 = piff.getPSF( galsim.PositionD(20, 10), smooth=True, ).drawImage(nx=53, ny=53, scale=0.263).array - piff = DES_Piff(piff_fname, psf_kwargs=PSF_KWARGS["g"], smooth=True) + piff = DES_SmoothPiff(piff_fname, psf_kwargs=PSF_KWARGS["g"], smooth=True) psf_im2 = piff.getPSF( galsim.PositionD(20, 10), ).drawImage(nx=53, ny=53, scale=0.263).array