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

@@ -127,7 +127,7 @@ class FiberTrace {
         ndarray::Array<double, 2, 1> const& profiles,
         ndarray::Array<bool, 2, 1> const& good,
         ndarray::Array<double, 1, 1> const& centers,
-        ndarray::Array<double, 1, 1> const& norm=ndarray::Array<double, 1, 1>()
+        ndarray::Array<Spectrum::ImageT, 1, 1> const& norm=ndarray::Array<Spectrum::ImageT, 1, 1>()
     );
 
   private:

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<lsst::afw::image::MaskPixel>;
     using ImageArray = ndarray::Array<ImageT, 1, 1>;
     using ConstImageArray = ndarray::Array<const ImageT, 1, 1>;

include/pfs/drp/stella/SpectrumSet.h

@@ -20,8 +20,8 @@ class SpectrumSet {
     using SpectrumPtr = std::shared_ptr<Spectrum>;
     using ConstSpectrumPtr = std::shared_ptr<Spectrum const>;
     using Collection = std::vector<SpectrumPtr>;
-    using ImageArray = ndarray::Array<double, 2, 1>;
-    using CovarianceArray = ndarray::Array<double, 3, 1>;
+    using ImageArray = ndarray::Array<Spectrum::ImageT, 2, 1>;
+    using CovarianceArray = ndarray::Array<Spectrum::ImageT, 3, 1>;
     using MaskArray = ndarray::Array<lsst::afw::image::MaskPixel, 2, 1>;
     using WavelengthArray = ndarray::Array<double, 2, 1>;
     using iterator = Collection::iterator;

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

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")

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)

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<ndarray::Array<double const, 1, 1> const&,
+            py::overload_cast<ndarray::Array<Spectrum::ImageT const, 1, 1> const&,
                               ndarray::Array<lsst::afw::image::MaskPixel const, 1, 1> const&,
-                              float, float, lsst::afw::image::MaskPixel, std::size_t>(&fitLine<double>),
+                              float, float, lsst::afw::image::MaskPixel,
+                              std::size_t>(&fitLine<Spectrum::ImageT>),
             "flux"_a, "mask"_a, "peakPosition"_a, "rmsSize"_a, "badBitMask"_a,
             "fittingHalfSize"_a=0);
 

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)):

src/FiberTrace.cc

@@ -82,7 +82,7 @@ FiberTrace<ImageT, MaskT, VarianceT> FiberTrace<ImageT, MaskT, VarianceT>::fromP
     ndarray::Array<double, 2, 1> const& profiles,
     ndarray::Array<bool, 2, 1> const& good,
     ndarray::Array<double, 1, 1> const& centers,
-    ndarray::Array<double, 1, 1> const& norm
+    ndarray::Array<Spectrum::ImageT, 1, 1> const& norm
 ) {
     int const width = dims.getX();
     std::size_t const height = dims.getY();

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)
 

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)

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):

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"])

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

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):