From 228ecc352d60c2e31d013d1f5fd657fa9fb142be Mon Sep 17 00:00:00 2001 From: Paul Price Date: Tue, 3 May 2022 15:49:59 -0400 Subject: [PATCH] use double precision only for wavelength Wavelength needs to be in double precision in order to preserve the precision in the measurements, but other quantities (esp. flux) don't need double precision, and datamodel I/O has been reverted to single-precision for these to save disk space. This commit adapts to these changes. --- include/pfs/drp/stella/FiberTrace.h | 2 +- include/pfs/drp/stella/Spectrum.h | 4 ++-- include/pfs/drp/stella/SpectrumSet.h | 4 ++-- python/pfs/drp/stella/fiberProfile.py | 8 +++++--- python/pfs/drp/stella/findLines.py | 2 +- python/pfs/drp/stella/fitContinuum.py | 6 +++--- python/pfs/drp/stella/fitLine.cc | 5 +++-- python/pfs/drp/stella/mergeArms.py | 2 +- src/FiberTrace.cc | 2 +- tests/test_PfsArm.py | 14 +++++++------- tests/test_Spectrum.py | 10 +++++----- tests/test_SpectrumSet.py | 8 ++++---- tests/test_extractSpectra.py | 10 +++++----- tests/test_fitLine.py | 2 +- tests/test_profileNorm.py | 2 +- 15 files changed, 42 insertions(+), 39 deletions(-) diff --git a/include/pfs/drp/stella/FiberTrace.h b/include/pfs/drp/stella/FiberTrace.h index 23649696a..707e9a55d 100644 --- a/include/pfs/drp/stella/FiberTrace.h +++ b/include/pfs/drp/stella/FiberTrace.h @@ -127,7 +127,7 @@ class FiberTrace { ndarray::Array const& profiles, ndarray::Array const& good, ndarray::Array const& centers, - ndarray::Array const& norm=ndarray::Array() + ndarray::Array const& norm=ndarray::Array() ); private: diff --git a/include/pfs/drp/stella/Spectrum.h b/include/pfs/drp/stella/Spectrum.h index e2083c21c..3e77ea7a4 100644 --- a/include/pfs/drp/stella/Spectrum.h +++ b/include/pfs/drp/stella/Spectrum.h @@ -13,8 +13,8 @@ namespace pfs { namespace drp { namespace stella { */ class Spectrum { public: - using ImageT = double; - using VarianceT = double; + using ImageT = float; + using VarianceT = float; using Mask = lsst::afw::image::Mask; using ImageArray = ndarray::Array; using ConstImageArray = ndarray::Array; diff --git a/include/pfs/drp/stella/SpectrumSet.h b/include/pfs/drp/stella/SpectrumSet.h index a3905e64a..b58221f7d 100644 --- a/include/pfs/drp/stella/SpectrumSet.h +++ b/include/pfs/drp/stella/SpectrumSet.h @@ -20,8 +20,8 @@ class SpectrumSet { using SpectrumPtr = std::shared_ptr; using ConstSpectrumPtr = std::shared_ptr; using Collection = std::vector; - using ImageArray = ndarray::Array; - using CovarianceArray = ndarray::Array; + using ImageArray = ndarray::Array; + using CovarianceArray = ndarray::Array; using MaskArray = ndarray::Array; using WavelengthArray = ndarray::Array; using iterator = Collection::iterator; diff --git a/python/pfs/drp/stella/fiberProfile.py b/python/pfs/drp/stella/fiberProfile.py index 9df34d8ef..4bf4f4d0b 100644 --- a/python/pfs/drp/stella/fiberProfile.py +++ b/python/pfs/drp/stella/fiberProfile.py @@ -27,7 +27,7 @@ class FiberProfile: profiles : array_like of `float`, shape ``(N, M)`` Profiles for each swath, each of width ``M = 2*(radius + 1)*oversample + 1``. - norm : array_like of `float`, optional + norm : array_like of `np.float32`, optional Normalisation for each spectral pixel. """ def __init__(self, radius, oversample, rows, profiles, norm=None): @@ -220,7 +220,8 @@ def makeFiberTraceFromDetectorMap(self, detectorMap, fiberId): """ return FiberTrace.fromProfile(fiberId, detectorMap.bbox.getDimensions(), self.radius, self.oversample, self.rows, self.profiles, ~self.profiles.mask, - detectorMap.getXCenter(fiberId), self.norm) + detectorMap.getXCenter(fiberId), + self.norm.astype(np.float32) if self.norm is not None else None) def makeFiberTrace(self, dimensions, centerFunc, fiberId): """Make a FiberTrace object @@ -244,7 +245,8 @@ def makeFiberTrace(self, dimensions, centerFunc, fiberId): """ rows = np.arange(dimensions.getY(), dtype=float) return FiberTrace.fromProfile(fiberId, dimensions, self.radius, self.oversample, self.rows, - self.profiles, ~self.profiles.mask, centerFunc(rows), self.norm) + self.profiles, ~self.profiles.mask, centerFunc(rows), + self.norm.astype(np.float32) if self.norm is not None else None) def extractSpectrum(self, maskedImage, detectorMap, fiberId): """Extract a single spectrum from an image diff --git a/python/pfs/drp/stella/findLines.py b/python/pfs/drp/stella/findLines.py index 9fb57cda1..1fa4bb377 100644 --- a/python/pfs/drp/stella/findLines.py +++ b/python/pfs/drp/stella/findLines.py @@ -125,7 +125,7 @@ def convolve(self, spectrum, continuum=None): size = 2*halfSize + 1 xx = np.arange(size, dtype=float) - halfSize sigma = self.config.width - kernel = np.exp(-0.5*xx**2/sigma**2)/sigma/np.sqrt(2.0*np.pi) + kernel = (np.exp(-0.5*xx**2/sigma**2)/sigma/np.sqrt(2.0*np.pi)).astype(spectrum.flux.dtype) flux = np.convolve(spectrum.normFlux if continuum is None else spectrum.normFlux - continuum, kernel, mode="same") diff --git a/python/pfs/drp/stella/fitContinuum.py b/python/pfs/drp/stella/fitContinuum.py index d8fa58565..bdf5e4948 100644 --- a/python/pfs/drp/stella/fitContinuum.py +++ b/python/pfs/drp/stella/fitContinuum.py @@ -141,13 +141,13 @@ def _fitContinuumImpl(self, values, good): fit : `numpy.ndarray`, floating-point Fit array. """ - indices = np.arange(len(values), dtype=float) + indices = np.arange(len(values), dtype=values.dtype) knots, binned = binData(indices, values, good, self.config.numKnots) use = np.isfinite(knots) & np.isfinite(binned) if not np.any(use): raise FitContinuumError("No finite knots when fitting continuum") interp = makeInterpolate(knots[use], binned[use], self.fitType) - fit = np.array(interp.interpolate(indices)) + fit = np.array(interp.interpolate(indices)).astype(values.dtype) if lsstDebug.Info(__name__).plot: import matplotlib.pyplot as plt @@ -285,4 +285,4 @@ def binData(xx, yy, good, numBins): select = good[low:high] xBinned[ii] = np.median(xx[low:high][select]) if np.any(select) else np.nan yBinned[ii] = np.median(yy[low:high][select]) if np.any(select) else np.nan - return xBinned, yBinned + return xBinned.astype(xx.dtype), yBinned.astype(yy.dtype) diff --git a/python/pfs/drp/stella/fitLine.cc b/python/pfs/drp/stella/fitLine.cc index df08a6f13..bfdb36486 100644 --- a/python/pfs/drp/stella/fitLine.cc +++ b/python/pfs/drp/stella/fitLine.cc @@ -36,9 +36,10 @@ PYBIND11_PLUGIN(fitLine) { lsst::utils::python::addOutputOp(cls, "__repr__"); mod.def("fitLine", - py::overload_cast const&, + py::overload_cast const&, ndarray::Array const&, - float, float, lsst::afw::image::MaskPixel, std::size_t>(&fitLine), + float, float, lsst::afw::image::MaskPixel, + std::size_t>(&fitLine), "flux"_a, "mask"_a, "peakPosition"_a, "rmsSize"_a, "badBitMask"_a, "fittingHalfSize"_a=0); diff --git a/python/pfs/drp/stella/mergeArms.py b/python/pfs/drp/stella/mergeArms.py index 26861dddd..67ab8ea5b 100644 --- a/python/pfs/drp/stella/mergeArms.py +++ b/python/pfs/drp/stella/mergeArms.py @@ -328,7 +328,7 @@ def normalizeSpectra(self, spectra): assert all(np.all(ss.fiberId == fiberId) for ss in spectra) # Consistent fibers # Collect normalisations from each arm, interpolated to common wavelength sampling - norm = np.zeros((fiberId.size, wavelength.size)) + norm = np.zeros((fiberId.size, wavelength.size), dtype=np.float32) for ss in spectra: badBitmask = ss.flags.get(*self.config.mask) for ii in range(len(ss)): diff --git a/src/FiberTrace.cc b/src/FiberTrace.cc index e1fd62c17..ba2fca83e 100644 --- a/src/FiberTrace.cc +++ b/src/FiberTrace.cc @@ -82,7 +82,7 @@ FiberTrace FiberTrace::fromP ndarray::Array const& profiles, ndarray::Array const& good, ndarray::Array const& centers, - ndarray::Array const& norm + ndarray::Array const& norm ) { int const width = dims.getX(); std::size_t const height = dims.getY(); diff --git a/tests/test_PfsArm.py b/tests/test_PfsArm.py index a1bf0b586..3146423aa 100644 --- a/tests/test_PfsArm.py +++ b/tests/test_PfsArm.py @@ -30,11 +30,11 @@ def setUp(self): self.identity = Identity(self.visit, self.arm, self.spectrograph, self.pfsDesignId) self.fiberId = self.synthConfig.fiberId self.wavelength = np.vstack([np.linspace(600, 900, self.synthConfig.height)]*num) - self.flux = self.rng.uniform(size=shape) - self.mask = self.rng.uniform(high=2**31, size=shape).astype(np.int32) - self.sky = self.rng.uniform(size=shape) - self.norm = self.rng.uniform(size=shape) - self.covar = self.rng.uniform(size=(num, 3, self.synthConfig.height)) + self.flux = self.rng.uniform(size=shape).astype(np.float32) + self.mask = self.rng.uniform(high=2**31, size=shape).astype(np.uint32) + self.sky = self.rng.uniform(size=shape).astype(np.float32) + self.norm = self.rng.uniform(size=shape).astype(np.float32) + self.covar = self.rng.uniform(size=(num, 3, self.synthConfig.height)).astype(np.float32) self.flags = MaskHelper() self.metadata = dict(FOO=12345, BAR=0.9876) # FITS keywords get capitalised @@ -48,9 +48,9 @@ def makeSpectra(self): def assertSpectra(self, spectra): """Check that the spectra match what's expected""" for name in ("visit", "arm", "spectrograph", "pfsDesignId"): - self.assertEqual(getattr(spectra.identity, name), getattr(self, name)) + self.assertEqual(getattr(spectra.identity, name), getattr(self, name), msg=name) for name in ("fiberId", "wavelength", "flux", "mask", "sky", "norm", "covar"): - self.assertFloatsEqual(getattr(spectra, name), getattr(self, name)) + self.assertFloatsEqual(getattr(spectra, name), getattr(self, name), msg=name) self.assertDictEqual(spectra.flags.flags, self.flags.flags) self.assertDictEqual({**self.metadata, **spectra.metadata}, spectra.metadata) diff --git a/tests/test_Spectrum.py b/tests/test_Spectrum.py index de7281b63..ed5417d26 100644 --- a/tests/test_Spectrum.py +++ b/tests/test_Spectrum.py @@ -23,12 +23,12 @@ def setUp(self): # Spectrum self.length = 123 self.fiberId = 456 - self.image = self.rng.uniform(size=self.length).astype(float) + self.image = self.rng.uniform(size=self.length).astype(np.float32) self.mask = lsst.afw.image.Mask(self.length, 1) self.mask.array[:] = self.rng.randint(0, 2**30, self.length).astype(np.int32) - self.background = self.rng.uniform(size=self.length).astype(float) - self.norm = self.rng.uniform(size=self.length).astype(float) - self.covariance = self.rng.uniform(size=(3, self.length)).astype(float) + self.background = self.rng.uniform(size=self.length).astype(np.float32) + self.norm = self.rng.uniform(size=self.length).astype(np.float32) + self.covariance = self.rng.uniform(size=(3, self.length)).astype(np.float32) self.wavelengthArray = np.arange(1, self.length + 1, dtype=float) def makeSpectrum(self): @@ -145,7 +145,7 @@ def testBasics(self): # Set the variance and ensure it changes the covariance spectrum.covariance[:] = self.covariance self.assertFloatsEqual(spectrum.covariance, self.covariance) - variance = self.rng.uniform(size=self.length) + variance = self.rng.uniform(size=self.length).astype(np.float32) spectrum.variance = variance self.covariance[0] = variance self.assertFloatsEqual(spectrum.covariance, self.covariance) diff --git a/tests/test_SpectrumSet.py b/tests/test_SpectrumSet.py index 84ecbddbe..0c2154ec5 100644 --- a/tests/test_SpectrumSet.py +++ b/tests/test_SpectrumSet.py @@ -21,12 +21,12 @@ def setUp(self): # Spectrum self.length = 123 self.fiberId = 456 - self.image = self.rng.uniform(size=self.length).astype(float) + self.image = self.rng.uniform(size=self.length).astype(np.float32) self.mask = lsst.afw.image.Mask(self.length, 1) self.mask.array[:] = self.rng.randint(0, 2**30, self.length).astype(np.int32) - self.background = self.rng.uniform(size=self.length).astype(float) - self.norm = self.rng.uniform(size=self.length).astype(float) - self.covariance = self.rng.uniform(size=(3, self.length)).astype(float) + self.background = self.rng.uniform(size=self.length).astype(np.float32) + self.norm = self.rng.uniform(size=self.length).astype(np.float32) + self.covariance = self.rng.uniform(size=(3, self.length)).astype(np.float32) self.wavelengthArray = np.arange(1, self.length + 1, dtype=float) def makeSpectrum(self): diff --git a/tests/test_extractSpectra.py b/tests/test_extractSpectra.py index a8cbc7611..5aeccf183 100644 --- a/tests/test_extractSpectra.py +++ b/tests/test_extractSpectra.py @@ -37,7 +37,7 @@ def setUp(self): lsst.afw.image.makeExposure(self.image), detectorMap=self.detMap ).profiles for fiberId in self.fiberProfiles: - self.fiberProfiles[fiberId].norm = np.full(self.synthConfig.height, self.flux, dtype=float) + self.fiberProfiles[fiberId].norm = np.full(self.synthConfig.height, self.flux, dtype=np.float32) self.fiberTraces = self.fiberProfiles.makeFiberTracesFromDetectorMap(self.detMap) self.assertEqual(len(self.fiberTraces), self.synthConfig.numFibers) @@ -78,7 +78,7 @@ def assertSpectra(self, spectra, flux=None, mask=None): expectMask = 0 bbox = ft.trace.getBBox() - expectNorm = np.zeros(self.synthConfig.height, dtype=float) + expectNorm = np.zeros(self.synthConfig.height, dtype=np.float32) expectNorm[bbox.getMinY():bbox.getMaxY() + 1] = ft.trace.image.array.sum(axis=1) self.assertEqual(ss.fiberId, ft.fiberId) @@ -126,7 +126,7 @@ def testUndersizedFiberTrace(self): newTrace = type(trace)(trace.trace.Factory(trace.trace, newBox), trace.fiberId) self.fiberTraces[index] = newTrace spectra = self.fiberTraces.extractSpectra(self.image) - expectFlux = np.zeros(self.synthConfig.height, dtype=float) + expectFlux = np.zeros(self.synthConfig.height, dtype=np.float32) expectFlux[middle:] = self.flux expectMask = np.zeros(self.synthConfig.height, dtype=np.int32) expectMask[:middle] = trace.trace.mask.getPlaneBitMask("NO_DATA") @@ -137,7 +137,7 @@ def testUndersizedFiberTrace(self): newTrace = type(trace)(trace.trace.Factory(trace.trace, newBox), trace.fiberId) self.fiberTraces[index] = newTrace spectra = self.fiberTraces.extractSpectra(self.image) - expectFlux = np.zeros(self.synthConfig.height, dtype=float) + expectFlux = np.zeros(self.synthConfig.height, dtype=np.float32) expectFlux[:middle] = self.flux expectMask = np.zeros(self.synthConfig.height, dtype=np.int32) expectMask[middle:] = trace.trace.mask.getPlaneBitMask("NO_DATA") @@ -154,7 +154,7 @@ def testSingular(self): index = self.synthConfig.numFibers//2 row = self.synthConfig.height//2 fiberId = self.synthConfig.fiberId[index] - expectFlux = np.full(self.synthConfig.height, self.flux, dtype=float) + expectFlux = np.full(self.synthConfig.height, self.flux, dtype=np.float32) expectMask = np.zeros(self.synthConfig.height, dtype=np.int32) expectFlux[row] = 0.0 expectMask[row] = self.image.mask.getPlaneBitMask(["BAD_FIBERTRACE", "NO_DATA"]) diff --git a/tests/test_fitLine.py b/tests/test_fitLine.py index 1c71f5878..8e8499aef 100644 --- a/tests/test_fitLine.py +++ b/tests/test_fitLine.py @@ -61,7 +61,7 @@ def makeSpectrum(length, center, amplitude, rmsSize, bgConst=0.0, bgSlope=0.0): return spectrum -class FindLinesTestCase(lsst.utils.tests.TestCase): +class FitLinesTestCase(lsst.utils.tests.TestCase): def setUp(self): self.length = 256 self.center = 123.45 diff --git a/tests/test_profileNorm.py b/tests/test_profileNorm.py index a3707c585..3d5ab2cba 100644 --- a/tests/test_profileNorm.py +++ b/tests/test_profileNorm.py @@ -59,7 +59,7 @@ def testProfileNorm(self): newSpectra = traces.extractSpectra(image) for new, old in zip(newSpectra, spectra): self.assertFloatsAlmostEqual(new.norm, old.normFlux, rtol=1.0e-7) - self.assertFloatsAlmostEqual(new.normFlux, 1.0, atol=1.0e-7) + self.assertFloatsAlmostEqual(new.normFlux, 1.0, atol=1.0e-6) class TestMemory(lsst.utils.tests.MemoryTestCase):