Fixes based on review comments.

.github/workflows/unit-tests.yml

@@ -18,11 +18,6 @@ jobs:
             python-version: 3.9
             toxenv: build_docs
 
-          - name: Python 3.9 with developer version of astropy and numpy
-            os: ubuntu-latest
-            python-version: 3.9
-            toxenv: py39-test-devdeps
-
           - name: Python 3.8 with minimal dependencies
             os: ubuntu-latest
             python-version: '3.8'

banzai/bias.py

@@ -17,6 +17,8 @@ def calibration_type(self):
         return 'BIAS'
 
     def make_master_calibration_frame(self, images):
+        for image in images:
+            image /= image.n_sub_exposures
         master_image = super(BiasMaker, self).make_master_calibration_frame(images)
         master_image.bias_level = np.mean([image.bias_level for image in images if image.bias_level is not None])
         return master_image
@@ -65,6 +67,8 @@ def __init__(self, runtime_context):
     def do_stage(self, image):
         bias_level = stats.sigma_clipped_mean(image.data, 3.5, mask=image.mask)
         image -= bias_level
+        # This is only run for bias frames and n_sub_exposures should always be 1
+        # but we divide it here for consistency to cover pathological uses cases
         image.meta['BIASLVL'] = bias_level / image.n_sub_exposures, 'Bias level that was removed after overscan'
         return image
 

banzai/data.py

@@ -132,7 +132,8 @@ def __init__(self, data: Union[np.array, Table], meta: fits.Header,
                  mask: np.array = None, name: str = '', uncertainty: np.array = None, memmap=True):
         super().__init__(data=data, meta=meta, mask=mask, name=name, memmap=memmap)
         if uncertainty is None:
-            uncertainty = self.read_noise * self.n_sub_exposures * np.ones(data.shape, dtype=data.dtype) / self.gain
+            uncertainty = self.read_noise * np.sqrt(self.n_sub_exposures) * np.ones(data.shape, dtype=data.dtype)
+            uncertainty /= self.gain
         self.uncertainty = self._init_array(uncertainty)
         self._detector_section = Section.parse_region_keyword(self.meta.get('DETSEC'))
         self._data_section = Section.parse_region_keyword(self.meta.get('DATASEC'))
@@ -185,7 +186,7 @@ def to_fits(self, context):
 
     def __del__(self):
         super().__del__()
-        del self.uncertainty
+        del self._uncertainty
 
     def __isub__(self, value):
         if isinstance(value, CCDData):
@@ -201,9 +202,14 @@ def __sub__(self, other):
         return type(self)(data=self.data - other.data, meta=self.meta, mask=self.mask | other.mask,
                           uncertainty=uncertainty)
 
-    def add_uncertainty(self, readnoise: np.array):
-        self._validate_array(readnoise)
-        self.uncertainty = self._init_array(readnoise)
+    @property
+    def uncertainty(self):
+        return self._uncertainty
+
+    @uncertainty.setter
+    def uncertainty(self, value: np.array):
+        self._validate_array(value)
+        self._uncertainty = self._init_array(value)
 
     def signal_to_noise(self):
         return np.abs(self.data) / self.uncertainty
@@ -318,6 +324,10 @@ def n_sub_exposures(self):
             n_exposures = 1
         return n_exposures
 
+    @n_sub_exposures.setter
+    def n_sub_exposures(self, value):
+        self.meta['NSUBREAD'] = value
+
     def rebin(self, binning):
         # TODO: Implement me
         return self

banzai/frames.py

@@ -82,6 +82,10 @@ def background(self, value):
     def n_sub_exposures(self):
         return self.primary_hdu.n_sub_exposures
 
+    @n_sub_exposures.setter
+    def n_sub_exposures(self, value):
+        self.primary_hdu.n_sub_exposures = value
+
     @abc.abstractmethod
     def save_processing_metadata(self, context):
         pass

banzai/readnoise.py

@@ -9,7 +9,7 @@ class ReadNoiseLoader(CalibrationUser):
     def apply_master_calibration(self, image, master_calibration_image):
         try:
             for image_extension, readnoise_extension in zip(image.ccd_hdus, master_calibration_image.ccd_hdus):
-                image_extension.add_uncertainty(readnoise_extension.data * np.sqrt(image_extension.n_sub_exposures))
+                image_extension.uncertainty = readnoise_extension.data * np.sqrt(image_extension.n_sub_exposures)
         except:
             logger.error(f"Can't add READNOISE to image, stopping reduction: {format_exception()}", image=image)
             return None

banzai/tests/test_bias_level_subtractor.py

@@ -52,3 +52,16 @@ def test_bias_master_level_subtraction_is_reasonable(set_random_seed):
 
     np.testing.assert_allclose(np.zeros(image.data.shape), image.data, atol=8 * read_noise)
     np.testing.assert_allclose(image.meta.get('BIASLVL'), input_bias, atol=1.0)
+
+
+def test_multiread_bias_level():
+    input_bias = 500.0
+    nx, ny = 101, 103
+    nreads = 5
+    subtractor = BiasMasterLevelSubtractor(None)
+    image = FakeLCOObservationFrame(hdu_list=[FakeCCDData(image_multiplier=input_bias*nreads, nx=nx, ny=ny)])
+    image.n_sub_exposures = nreads
+    image = subtractor.do_stage(image)
+
+    np.testing.assert_allclose(np.zeros(image.data.shape), image.data, atol=1.0)
+    np.testing.assert_allclose(image.meta.get('BIASLVL'), input_bias, atol=1.0)

banzai/tests/test_bias_maker.py

@@ -71,7 +71,7 @@ def test_makes_a_sensible_master_bias(mock_namer):
 
     images = [FakeLCOObservationFrame(hdu_list=[FakeCCDData(data=np.random.normal(0.0, size=(99, 99), scale=expected_readnoise),
                                                             bias_level=0.0,
-                                                            read_noise= expected_readnoise,
+                                                            read_noise=expected_readnoise,
                                                             meta=Header({'DATE-OBS': '2019-12-04T14:34:00',
                                                                          'DETSEC': '[1:100,1:100]',
                                                                          'DATASEC': '[1:100,1:100]',
@@ -84,3 +84,29 @@ def test_makes_a_sensible_master_bias(mock_namer):
     assert np.abs(np.mean(master_bias)) < 0.1
     actual_readnoise = np.std(master_bias)
     assert np.abs(actual_readnoise - expected_readnoise / (nimages ** 0.5)) < 0.2
+
+
+def test_multiread_bias_maker():
+    nimages = 5
+    nreads = 20
+    images = []
+    pattern_scale = 8
+    nx = 101
+    ny = 105
+    bias_level = 100.0
+    bias_pattern = np.random.normal(0.0, pattern_scale, size=(ny, nx))
+    for i in range(nimages):
+        data = nreads * bias_pattern
+        image = FakeLCOObservationFrame(hdu_list=[FakeCCDData(data=data, bias_level=bias_level,
+                                                              meta=Header({'DATE-OBS': '2019-12-04T14:34:00',
+                                                                           'DETSEC': f'[1:{nx},1:{ny}]',
+                                                                           'DATASEC': f'[1:{nx},1:{ny}]',
+                                                                           'OBSTYPE': 'BIAS', 'RA': 0.0, 'DEC': 0.0}))])
+        image.n_sub_exposures = nreads
+        images.append(image)
+    maker = BiasMaker(FakeContext())
+    stacked_image = maker.do_stage(images)
+    np.testing.assert_allclose(stacked_image.meta['BIASLVL'], bias_level, atol=0.1)
+
+    # With 20 x 5 reads, we should get down to better than 1 count at a read_noise of 9
+    np.testing.assert_allclose(stacked_image.data, bias_pattern, atol=1.0)

banzai/tests/test_bias_subtractor.py

@@ -69,10 +69,30 @@ def test_bias_subtraction_is_reasonable(mock_master_cal_name, mock_master_frame)
     input_level = 2000.0
     nx = 101
     ny = 103
-    mock_master_frame.return_value = FakeLCOObservationFrame(hdu_list=[FakeCCDData(bias_level=input_bias,
-                                                                                   data=np.random.normal(0.0, input_readnoise, size=(ny, nx)))])
+    stacked_hdu = FakeCCDData(bias_level=input_bias, data=np.random.normal(0.0, input_readnoise, size=(ny, nx)))
+    mock_master_frame.return_value = FakeLCOObservationFrame(hdu_list=[stacked_hdu])
     subtractor = BiasSubtractor(FakeContext())
     image = FakeLCOObservationFrame(hdu_list=[FakeCCDData(image_multiplier=input_level)])
     image = subtractor.do_stage(image)
     assert np.abs(image.meta.get('BIASLVL') - input_bias) < 1.0
     assert np.abs(np.mean(image.data) - (input_level - input_bias)) < 1.0
+
+
+@mock.patch('banzai.lco.LCOFrameFactory.open')
+@mock.patch('banzai.calibrations.CalibrationUser.get_calibration_file_info')
+def test_multiread_bias_subtraction(mock_super_cal_name, mock_super_frame):
+    mock_super_cal_name.return_value = {'filename': 'test.fits'}
+    input_bias = 500.0
+    nreads = 5
+    nx = 101
+    ny = 103
+    bias_pattern = np.random.normal(0.0, 15.0, size=(ny, nx))
+    stacked_hdu = FakeCCDData(bias_level=input_bias, data=bias_pattern)
+    mock_super_frame.return_value = FakeLCOObservationFrame(hdu_list=[stacked_hdu])
+    data = bias_pattern * nreads + input_bias * nreads
+    image = FakeLCOObservationFrame(hdu_list=[FakeCCDData(bias_level=input_bias, data=data)])
+    image.n_sub_exposures = nreads
+    subtractor = BiasSubtractor(FakeContext())
+    image = subtractor.do_stage(image)
+    assert np.abs(image.meta.get('BIASLVL') - input_bias) < 1.0
+    assert np.abs(np.mean(image.data)) < 1.0

banzai/tests/test_end_to_end.py

@@ -159,7 +159,7 @@ def run_check_if_stacked_calibrations_are_in_db(raw_file_pattern, calibration_ty
 
 
 def observation_portal_side_effect(*args, **kwargs):
-    # To produce the mock observation portal resonse, we need to modify the response from
+    # To produce the mock observation portal response, we need to modify the response from
     # this type of url
     # https://observe.lco.global/api/observations/?enclosure=aqwa&telescope=0m4a&priority=&state=COMPLETED&time_span=&start_after=2024-07-23&start_before=&end_after=&end_before=2024-07-25&modified_after=&created_after=&created_before=&request_id=&request_group_id=&user=&proposal=calibrate&instrument_type=0M4-SCICAM-QHY600&configuration_type=SKY_FLAT&ordering=
     site = kwargs['params']['site']

banzai/tests/test_frames.py

@@ -140,7 +140,9 @@ def test_all_datatypes_wrong():
 
 
 def test_subtract():
-    test_data = FakeCCDData(image_multiplier=4, uncertainty=3)
+    nx = 101
+    ny = 103
+    test_data = FakeCCDData(image_multiplier=4, uncertainty=3*np.ones((ny, nx)), nx=nx, ny=ny)
     test_data -= 1
 
     assert (test_data.data == 3 * np.ones(test_data.data.shape)).all()
@@ -179,7 +181,9 @@ def test_trim():
 
 def test_init_poisson_uncertainties():
     # Make sure the uncertainties add in quadrature
-    test_data = FakeCCDData(image_multiplier=16, uncertainty=3)
+    nx = 101
+    ny = 103
+    test_data = FakeCCDData(image_multiplier=16, uncertainty=3 * np.ones((ny, nx)), nx=nx, ny=ny)
     test_data.init_poisson_uncertainties()
     assert (test_data.uncertainty == 5 * np.ones(test_data.data.shape)).all()
 

banzai/tests/utils.py

@@ -39,11 +39,11 @@ def __init__(self, image_multiplier=1.0, nx=101, ny=103, name='test_image', read
             self.mask = np.zeros(self.data.shape, dtype=np.uint8)
         else:
             self.mask = mask
+        self.memmap = True
         if uncertainty is None:
             self.uncertainty = self.read_noise * np.ones(self.data.shape, dtype=self.data.dtype)
         else:
             self.uncertainty = uncertainty
-        self.memmap = True
         if 'SATURATE' not in self.meta:
             self.meta['SATURATE'] = 65535.0
         if 'GAIN' not in self.meta:
@@ -69,10 +69,7 @@ def __init__(self, hdu_list=None, file_path='/tmp/test_image.fits', instrument=N
         self._file_path = file_path
         self.is_bad = False
         self.hdu_order = ['SCI', 'CAT', 'BPM', 'ERR']
-
-        @property
-        def n_sub_exposures(self):
-            return 1
+        self.n_sub_exposures = 1
         for keyword in kwargs:
             setattr(self, keyword, kwargs[keyword])