diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 821d7e1b..191e9cea 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -6,6 +6,10 @@ on: - main pull_request: null +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: tests: name: tests @@ -25,7 +29,7 @@ jobs: - uses: conda-incubator/setup-miniconda@v2 with: - python-version: 3.8 + python-version: 3.11 channels: conda-forge,defaults channel-priority: strict show-channel-urls: true @@ -36,7 +40,7 @@ jobs: - name: configure conda and install code shell: bash -l {0} run: | - mamba install pytest flake8 + conda uninstall --force --yes piff python -m pip install -e . -vv - name: run cli @@ -46,13 +50,19 @@ jobs: eastlake-src-extractor -dd run-eastlake-sim --help - - name: test + - name: test eastlake shell: bash -l {0} run: | git clone https://github.com/beckermr/des-test-data export TEST_DESDATA=`pwd`/des-test-data pytest -vvs eastlake + - name: test piff + shell: bash -l {0} + run: | + cd piff_package/tests + coverage run -m pytest -v + - name: lint shell: bash -l {0} run: | diff --git a/.gitignore b/.gitignore index 5a1411b3..0afd5052 100644 --- a/.gitignore +++ b/.gitignore @@ -132,3 +132,9 @@ eastlake/astromatic/swarp eastlake/astromatic/src-extractor src/sextractor*/* src/swarp*/* + +piff_package/* +!piff_package/apodize.patch +!piff_package/README.md +!piff_package/Piff-1.3.3.tar.gz +!piff_package/LICENSE diff --git a/eastlake/__init__.py b/eastlake/__init__.py index 24257f57..9894e882 100644 --- a/eastlake/__init__.py +++ b/eastlake/__init__.py @@ -11,4 +11,3 @@ 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 97e150bb..7f890fb2 100644 --- a/eastlake/des_piff.py +++ b/eastlake/des_piff.py @@ -10,14 +10,56 @@ logger = logging.getLogger(__name__) +PSF_KWARGS_COLOR_RANGES = { + "GI_COLOR": [0.0, 3.5], + "IZ_COLOR": [0.0, 0.65], +} +PSF_COLOR_DEFAULTS = { + "GI_COLOR": 1.1, + "IZ_COLOR": 0.34, +} + PSF_KWARGS = { - "g": {"GI_COLOR": 1.1}, - "r": {"GI_COLOR": 1.1}, - "i": {"GI_COLOR": 1.1}, - "z": {"IZ_COLOR": 0.34}, + "g": {"GI_COLOR": PSF_COLOR_DEFAULTS["GI_COLOR"]}, + "r": {"GI_COLOR": PSF_COLOR_DEFAULTS["GI_COLOR"]}, + "i": {"GI_COLOR": PSF_COLOR_DEFAULTS["GI_COLOR"]}, + "z": {"IZ_COLOR": PSF_COLOR_DEFAULTS["IZ_COLOR"]}, } +def _process_color_kwargs(piff_psf, kwargs): + if "GI_COLOR" in piff_psf.interp_property_names: + kwargs.pop("IZ_COLOR", None) + + if "GI_COLOR" in kwargs: + if ( + kwargs["GI_COLOR"] is None or + kwargs["GI_COLOR"] == "None" + ): + kwargs["GI_COLOR"] = PSF_COLOR_DEFAULTS["GI_COLOR"] + + if kwargs["GI_COLOR"] < PSF_KWARGS_COLOR_RANGES["GI_COLOR"][0]: + kwargs["GI_COLOR"] = PSF_KWARGS_COLOR_RANGES["GI_COLOR"][0] + if kwargs["GI_COLOR"] > PSF_KWARGS_COLOR_RANGES["GI_COLOR"][1]: + kwargs["GI_COLOR"] = PSF_KWARGS_COLOR_RANGES["GI_COLOR"][1] + + elif "IZ_COLOR" in piff_psf.interp_property_names: + kwargs.pop("GI_COLOR", None) + if "IZ_COLOR" in kwargs: + if ( + kwargs["IZ_COLOR"] is None or + kwargs["IZ_COLOR"] == "None" + ): + kwargs["IZ_COLOR"] = PSF_COLOR_DEFAULTS["IZ_COLOR"] + + if kwargs["IZ_COLOR"] < PSF_KWARGS_COLOR_RANGES["IZ_COLOR"][0]: + kwargs["IZ_COLOR"] = PSF_KWARGS_COLOR_RANGES["IZ_COLOR"][0] + if kwargs["IZ_COLOR"] > PSF_KWARGS_COLOR_RANGES["IZ_COLOR"][1]: + kwargs["IZ_COLOR"] = PSF_KWARGS_COLOR_RANGES["IZ_COLOR"][1] + + return kwargs + + @functools.lru_cache(maxsize=200) def _read_piff(file_name): return piff.read( @@ -84,26 +126,7 @@ def _draw( # nice and big image size here cause this has been a problem image = galsim.ImageD(ncol=n_pix, nrow=n_pix, wcs=pixel_wcs) - 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 + psf_kwargs = _process_color_kwargs(self.getPiff(), psf_kwargs) offset = ( image_pos.x - int(image_pos.x + 0.5), @@ -126,8 +149,11 @@ def getPiff(self): return self._piff def getPSF( - self, image_pos, wcs=None, - n_pix=None, depixelize=False, + self, + image_pos, + wcs=None, + n_pix=None, + depixelize=False, gsparams=None, **kwargs ): @@ -175,7 +201,9 @@ def getPSF( return psf def getPSFImage( - self, image_pos, wcs=None, + self, + image_pos, + wcs=None, n_pix=None, gsparams=None, **kwargs diff --git a/eastlake/des_smoothpiff.py b/eastlake/des_smoothpiff.py deleted file mode 100644 index c3c4d5ff..00000000 --- a/eastlake/des_smoothpiff.py +++ /dev/null @@ -1,369 +0,0 @@ -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/steps/_meds.py b/eastlake/steps/_meds.py index 4191f334..400a778f 100644 --- a/eastlake/steps/_meds.py +++ b/eastlake/steps/_meds.py @@ -535,7 +535,10 @@ def _make_psf_data(self, stash, coadd_wcs, tilename, band, is_rejectlisted, img_ def _make_psf_data_piff( self, stash, coadd_wcs, tilename, band, is_rejectlisted, img_files, img_ext, head_files, ): - smooth = stash["psf_config"].get("smooth", False) + smooth = ( + stash["psf_config"].get("smooth", False) + or stash["psf_config"].get("type") == "DES_SmoothPiff" + ) if smooth: assert stash["draw_method"] == "auto" else: diff --git a/eastlake/steps/pizza_cutter.py b/eastlake/steps/pizza_cutter.py index fbda06ce..94f75e3e 100644 --- a/eastlake/steps/pizza_cutter.py +++ b/eastlake/steps/pizza_cutter.py @@ -128,7 +128,7 @@ def _prep_input_config_and_info_file(self, tilename, band, stash, odir, allbands with open(pzyml_pth, "r") as fp: pzyml = yaml.safe_load(fp.read()) - if stash["psf_config"]["type"] in ["DES_Piff"]: + if stash["psf_config"]["type"] in ["DES_Piff", "DES_SmoothPiff"]: # this is the default and let's make sure assert pzyml["single_epoch"]["psf_type"] == "piff" else: diff --git a/eastlake/tests/test_des_piff.py b/eastlake/tests/test_des_piff.py index 00765db2..57d3c075 100644 --- a/eastlake/tests/test_des_piff.py +++ b/eastlake/tests/test_des_piff.py @@ -7,7 +7,6 @@ import pytest from ..des_piff import PSF_KWARGS, DES_Piff -from ..des_smoothpiff import DES_SmoothPiff def _get_piff_file(): @@ -29,9 +28,10 @@ def _get_piff_file(): '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): +@pytest.mark.parametrize("cls", [DES_Piff]) +def test_des_piff_centering(x_offset, y_offset, cls): piff_fname = _get_piff_file() - piff = DES_Piff(piff_fname) + piff = cls(piff_fname) atol = 0.01 # test the image it makes @@ -85,9 +85,10 @@ def test_des_piff_centering(x_offset, y_offset): reason=( 'DES_Piff can only be tested if ' 'test data is at TEST_DESDATA')) -def test_des_piff_color(): +@pytest.mark.parametrize("cls", [DES_Piff]) +def test_des_piff_color(cls): piff_fname = _get_piff_file() - piff = DES_Piff(piff_fname) + piff = cls(piff_fname) psf_im1 = piff.getPSF( galsim.PositionD(20, 10), GI_COLOR=1.3, @@ -116,9 +117,10 @@ def test_des_piff_color(): reason=( 'DES_Piff can only be tested if ' 'test data is at TEST_DESDATA')) -def test_des_piff_raises(): +@pytest.mark.parametrize("cls", [DES_Piff]) +def test_des_piff_raises(cls): piff_fname = _get_piff_file() - piff = DES_Piff(piff_fname) + piff = cls(piff_fname) with pytest.raises(Exception) as e: piff.getPSF( galsim.PositionD(20, 10), @@ -140,9 +142,10 @@ def test_des_piff_raises(): reason=( 'DES_Piff can only be tested if ' 'test data is at TEST_DESDATA')) -def test_des_piff_smoke(): +@pytest.mark.parametrize("cls", [DES_Piff]) +def test_des_piff_smoke(cls): piff_fname = _get_piff_file() - piff = DES_Piff(piff_fname) + piff = cls(piff_fname) psf_im = piff.getPSF( galsim.PositionD(20, 10), **PSF_KWARGS["g"], @@ -162,65 +165,3 @@ def test_des_piff_smoke(): 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: - 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) - - 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_smooth(): - 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 - - 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_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 - - for im in [psf_im, psf_im1, psf_im2]: - assert np.all(np.isfinite(im)) - assert np.allclose(np.sum(im), 1) - y, x = np.unravel_index(np.argmax(im), im.shape) - cen = (im.shape[0]-1)/2 - assert y == cen - assert x == cen - - assert np.allclose(psf_im1, psf_im2) - assert not np.allclose(psf_im, psf_im1) - assert not np.allclose(psf_im, psf_im2) diff --git a/environment.yaml b/environment.yaml index c2723fba..0f3ac4d7 100644 --- a/environment.yaml +++ b/environment.yaml @@ -10,20 +10,33 @@ dependencies: - des-pizza-cutter-metadetect - esutil - fitsio - - galsim + - galsim<2.5 - importlib_metadata - joblib - lsstdesc-coord - meds - ngmix - numpy - - piff - pixmappy - psfex - pyyaml - scipy - setuptools + # piff stuff + - lsstdesc-coord>=1.2 + - treecorr>=4.3.1 + - fitsio>=1.0 + - matplotlib-base>=3.3 + - galsim>=2.3 + - treegp>=0.6 + - threadpoolctl>=3.1 + - nose + - coverage + - pytest<8 + - nbval + - ipykernel # host setup.py + - pip - setuptools_scm - setuptools_scm_git_archive # host build @@ -33,3 +46,7 @@ dependencies: - compilers - fftw - libtool + - gnuconfig + # eastlake testing + - flake8 + - pytest diff --git a/piff b/piff new file mode 120000 index 00000000..d335a1b4 --- /dev/null +++ b/piff @@ -0,0 +1 @@ +piff_package/piff \ No newline at end of file diff --git a/piff_package/LICENSE b/piff_package/LICENSE new file mode 100644 index 00000000..e37e374d --- /dev/null +++ b/piff_package/LICENSE @@ -0,0 +1,24 @@ +Copyright (c) 2016 by Mike Jarvis and the other collaborators on GitHub at +https://github.com/rmjarvis/Piff All rights reserved. + +Piff is free software: Redistribution and use in source and binary forms, +with or without modification, are permitted provided that the following +conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/piff_package/Piff-1.3.3.tar.gz b/piff_package/Piff-1.3.3.tar.gz new file mode 100644 index 00000000..3a1be016 Binary files /dev/null and b/piff_package/Piff-1.3.3.tar.gz differ diff --git a/piff_package/README.md b/piff_package/README.md new file mode 100644 index 00000000..1d335bf4 --- /dev/null +++ b/piff_package/README.md @@ -0,0 +1,5 @@ +# Piff for DES + +Piff has been vendored into eastlake in order to add apodization support to avoid square stamps. + +The code in this directory is unpacked and then installed as a package after the patch has been applied. diff --git a/piff_package/apodize.patch b/piff_package/apodize.patch new file mode 100644 index 00000000..17a03ca3 --- /dev/null +++ b/piff_package/apodize.patch @@ -0,0 +1,154 @@ +diff --git a/piff/_version.py b/piff/_version.py +index 4f4c623a..c1ba101f 100644 +--- a/piff/_version.py ++++ b/piff/_version.py +@@ -12,5 +12,5 @@ + # this list of conditions and the disclaimer given in the documentation + # and/or other materials provided with the distribution. + +-__version__ = '1.3.3' ++__version__ = '1.3.3.1' + __version_info__ = tuple(map(int, __version__.split('.'))) +diff --git a/piff/pixelgrid.py b/piff/pixelgrid.py +index 4d113c77..a7765afd 100644 +--- a/piff/pixelgrid.py ++++ b/piff/pixelgrid.py +@@ -25,6 +25,31 @@ + from .model import Model + from .star import Star, StarData, StarFit + ++APODIZE_PARAMS = (1.0 * 0.263, 4.25 * 0.263) ++ ++ ++def set_apodize_params(pars): ++ global APODIZE_PARAMS ++ APODIZE_PARAMS = pars ++ ++ ++def _ap_kern_kern(x, m, h): ++ # cumulative triweight kernel ++ y = (x - m) / h + 3 ++ apval = np.zeros_like(m) ++ msk = y > 3 ++ apval[msk] = 1 ++ msk = (y > -3) & (~msk) ++ apval[msk] = ( ++ -5 * y[msk] ** 7 / 69984 ++ + 7 * y[msk] ** 5 / 2592 ++ - 35 * y[msk] ** 3 / 864 ++ + 35 * y[msk] / 96 ++ + 1 / 2 ++ ) ++ return apval ++ ++ + class PixelGrid(Model): + """A PSF modeled as interpolation between a grid of points. + +@@ -445,6 +470,21 @@ def getProfile(self, params): + :returns: a galsim.GSObject instance + """ + im = galsim.Image(params.reshape(self.size,self.size), scale=self.scale) ++ ++ if APODIZE_PARAMS is not None: ++ xpix, ypix = im.get_pixel_centers() ++ # use_true_center = False below ++ dx = xpix - im.center.x ++ dy = ypix - im.center.y ++ r2 = dx**2 + dy**2 ++ ++ apwidth, aprad = APODIZE_PARAMS # in arcsec ++ _apwidth = apwidth / self.scale # convert to pixels ++ _aprad = aprad / self.scale # convert to pixels ++ ++ apim = im._array * _ap_kern_kern(_aprad, np.sqrt(r2), _apwidth) ++ im._array = apim / np.sum(apim) * np.sum(im._array) ++ + return galsim.InterpolatedImage(im, x_interpolant=self.interp, + use_true_center=False, flux=1.) + +diff --git a/requirements.txt b/requirements.txt +index d03b7b19..cfd6693b 100644 +--- a/requirements.txt ++++ b/requirements.txt +@@ -1,10 +1,10 @@ +-numpy>=1.17 ++numpy>=1.17,<2 + scipy>=1.2 + pyyaml>=5.1 + lsstdesc.coord>=1.2 + treecorr>=4.3.1 + fitsio>=1.0 + matplotlib>=3.3 +-galsim>=2.3 ++galsim>=2.3,<2.5 + treegp>=0.6 + threadpoolctl>=3.1 +diff --git a/tests/conftest.py b/tests/conftest.py +new file mode 100644 +index 00000000..cf6d85c5 +--- /dev/null ++++ b/tests/conftest.py +@@ -0,0 +1,4 @@ ++# turn off apodization ++import piff.pixelgrid ++ ++piff.pixelgrid.set_apodize_params(None) +diff --git a/tests/test_wcs.py b/tests/test_wcs.py +index 2a373671..68f54ab3 100644 +--- a/tests/test_wcs.py ++++ b/tests/test_wcs.py +@@ -747,6 +747,53 @@ def test_des_wcs(): + np.testing.assert_allclose(wcs3.toWorld(im.center).y, wcs1.toWorld(im.center).y, + rtol=0.04) + ++@timer ++def test_newdes_apodize(): ++ # This is a DES Y6 PSF file made by Robert Gruendl using python 2, so ++ # check that this also works correctly. ++ try: ++ import pixmappy ++ except ImportError: ++ print('pixmappy not installed. Skipping test_newdes()') ++ return ++ # Also make sure pixmappy is recent enough to work. ++ if 'exposure_file' not in pixmappy.GalSimWCS._opt_params: ++ print('pixmappy not recent enough version. Skipping test_newdes()') ++ return ++ ++ import piff ++ import piff.pixelgrid ++ ++ if __name__ == '__main__': ++ logger = piff.config.setup_logger(verbose=2) ++ else: ++ logger = piff.config.setup_logger(log_file='output/test_newdes.log') ++ ++ fname = os.path.join('input', 'D00232418_i_c19_r5006p01_piff-model.fits') ++ with warnings.catch_warnings(): ++ # This file was written with GalSim 2.1, and now raises a deprecation warning for 2.2. ++ warnings.simplefilter("ignore", galsim.GalSimDeprecationWarning) ++ warnings.simplefilter("ignore", DeprecationWarning) ++ psf = piff.PSF.read(fname, logger=logger) ++ ++ ims = [] ++ for appars in [None, (1.0 * 0.263, 4.25 * 0.263)]: ++ piff.pixelgrid.set_apodize_params(appars) ++ ims.append(psf.draw(x=103.3, y=592.0, logger=logger)) ++ ++ print('sum = ',ims[1].array.sum()) ++ assert not np.allclose(ims[0].array, ims[1].array) ++ assert np.allclose(ims[1].array[0, :], 0, rtol=1.e-2) ++ assert np.allclose(ims[1].array[-1, :], 0, rtol=1.e-2) ++ assert np.allclose(ims[1].array[:, 0], 0, rtol=1.e-2) ++ assert np.allclose(ims[1].array[:, -1], 0, rtol=1.e-2) ++ assert ims[1].array.sum() > 0 ++ np.testing.assert_allclose( ++ ims[0].array[23:26,22:25] / ims[0].array[23:26,22:25].sum(), ++ ims[1].array[23:26,22:25] / ims[1].array[23:26,22:25].sum(), ++ rtol=1.e-5, ++ ) ++ + if __name__ == '__main__': + #import cProfile, pstats + #pr = cProfile.Profile() diff --git a/setup.py b/setup.py index 60643a2a..30a26c3e 100644 --- a/setup.py +++ b/setup.py @@ -30,13 +30,13 @@ def _get_cc_flags_prefix(): if "CONDA_BUILD" in os.environ: print(">>> building in conda build", flush=True) cc = "${CC}" - cldflags = "\"${CFLAGS} -fcommon\"" + cldflags = "\"${CFLAGS} -fcommon -Wno-implicit-function-declaration\"" prefix_var = "PREFIX" elif all(v in os.environ for v in ["CONDA_PREFIX", "CC", "CFLAGS", "LDFLAGS"]): # assume compilers package is installed print(">>> building w/ conda-forge compilers package", flush=True) cc = os.environ["CC"] - cldflags = "\"" + os.environ["CFLAGS"] + " -fcommon\"" + cldflags = "\"" + os.environ["CFLAGS"] + " -fcommon -Wno-implicit-function-declaration\"" prefix_var = "CONDA_PREFIX" elif sys.platform == "linux": if "CONDA_PREFIX" in os.environ: @@ -47,6 +47,7 @@ def _get_cc_flags_prefix(): "-Wl,-rpath,${CONDA_PREFIX}/lib " "-Wl,-rpath-link,${CONDA_PREFIX}/lib " "-L${CONDA_PREFIX}/lib " + "-Wno-implicit-function-declaration " "-ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 " "-ffunction-sections -pipe -fcommon\"" ) @@ -66,6 +67,7 @@ def _get_cc_flags_prefix(): "-Wl,-rpath,${CONDA_PREFIX}/lib " "-L${CONDA_PREFIX}/lib " "-Wl,-pie -Wl,-headerpad_max_install_names " + "-Wno-implicit-function-declaration " "-Wl,-dead_strip_dylibs -fcommon\"" ) prefix_var = "CONDA_PREFIX" @@ -82,6 +84,23 @@ def _get_cc_flags_prefix(): return cc, cldflags, prefix_var +def _update_autotools(): + if ( + "CONDA_BUILD" in os.environ + and "BUILD_PREFIX" in os.environ + and os.path.exists(os.path.expandvars("${BUILD_PREFIX}/share/gnuconfig")) + ): + print(">>> updating config.sub and config.guess", flush=True) + _run_shell("mkdir -p autoconf && cp ${BUILD_PREFIX}/share/gnuconfig/config.* autoconf/.") + + elif ( + "CONDA_PREFIX" in os.environ + and os.path.exists(os.path.expandvars("${CONDA_PREFIX}/share/gnuconfig")) + ): + print(">>> updating config.sub and config.guess", flush=True) + _run_shell("mkdir -p autoconf && cp ${CONDA_PREFIX}/share/gnuconfig/config.* autoconf/.") + + def _build_swarp(): cc, cldflags, prefix_var = _get_cc_flags_prefix() cwd = os.path.abspath(os.getcwd()) @@ -90,14 +109,20 @@ def _build_swarp(): _run_shell("cp %s/src/swarp-2.40.1.tar.gz ." % cwd) _run_shell("tar -xzvf swarp-2.40.1.tar.gz") with pushd("swarp-2.40.1"): - if cc is not None and cldflags is not None: - _run_shell( - "CC=%s " - "CFLAGS=%s " - "./configure --prefix=%s" % (cc, cldflags, tmpdir) - ) - else: - _run_shell("./configure --prefix=%s" % tmpdir) + try: + _update_autotools() + if cc is not None and cldflags is not None: + _run_shell( + "CC=%s " + "CFLAGS=%s " + "./configure --prefix=%s" % (cc, cldflags, tmpdir) + ) + else: + _run_shell("./configure --prefix=%s" % tmpdir) + except Exception: + _run_shell("cat config.log") + raise + _run_shell("make") _run_shell("make install") _run_shell("cp bin/swarp %s/eastlake/astromatic/." % cwd) @@ -112,6 +137,7 @@ def _build_src_ext(): _run_shell("tar -xzvf src-extractor-2.24.4.tar.gz") with pushd("sextractor-2.24.4"): try: + _update_autotools() _run_shell("sh autogen.sh") if cc is not None and cldflags is not None: _run_shell( @@ -155,22 +181,64 @@ def run(self): _build_src_ext() super().run() +# unpack and patch Piff +for fname in os.listdir("piff_package"): + if fname not in ["Piff-1.3.3.tar.gz", "apodize.patch", "README.md", "LICENSE"]: + subprocess.run( + ["rm", "-rf", fname], + check=True, + cwd="piff_package", + ) + +subprocess.run( + ["tar", "--strip-components=1", "-xzvf", "Piff-1.3.3.tar.gz"], + check=True, + cwd="piff_package", +) +subprocess.run( + ["patch", "-p1", "-u", "-i", "apodize.patch"], + check=True, + cwd="piff_package", +) scripts = glob.glob('./bin/*') scripts = [os.path.basename(f) for f in scripts if f[-1] != '~'] scripts = [os.path.join('bin', s) for s in scripts] +piff_scripts = [ + os.path.join("piff_package", "scripts", s) + for s in ['piffify', 'plotify', 'meanify'] +] + +pkgs = find_packages( + where=".", + include=["eastlake", "eastlake.*", "piff", "piff.*"], + exclude=[ + "tests", "tests.*", + "docs", "docs.*", + "devel", "devel.*", + "examples", "examples.*", + "scripts", "scripts.*" + ], +) + +print("\n\nfind_packages: %r\n\n" % pkgs, flush=True) + +assert "piff" in pkgs +assert "eastlake" in pkgs + setup( name='eastlake', - packages=find_packages(), + packages=pkgs, include_package_data=True, - scripts=scripts, + scripts=scripts + piff_scripts, cmdclass={ 'build_py': build_py, 'build_ext': build_ext, }, package_data={ "eastlake": ["astromatic/*", "data/*", "config/*"], + "piff": ["piff/pupils/*"], }, use_scm_version=True, setup_requires=['setuptools_scm', 'setuptools_scm_git_archive'],