Repository: feder-observatory/stellarphot

Issue:
Calculation of AIJ relative fluxes fails
A consequence of #421 is the calculation of `aij_relative_flux` fails because it relies on aggregation to calculate the relative flux. 

Discussion:

Patch:
diff --git a/stellarphot/differential_photometry/aij_rel_fluxes.py b/stellarphot/differential_photometry/aij_rel_fluxes.py
index 730479d..cf6a37c 100644
--- a/stellarphot/differential_photometry/aij_rel_fluxes.py
+++ b/stellarphot/differential_photometry/aij_rel_fluxes.py
@@ -1,7 +1,7 @@
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import SkyCoord
-from astropy.table import Table
+from astropy.table import QTable, Table
 
 __all__ = ["add_in_quadrature", "calc_aij_relative_flux"]
 
@@ -125,9 +125,11 @@ def calc_aij_relative_flux(
     # stars.
 
     comp_fluxes = star_data["date-obs", counts_column_name, error_column_name][good]
-    # print(np.isnan(comp_fluxes[flux_column_name]).sum(),
-    #       np.isnan(comp_fluxes[error_column_name]).sum())
-    # print(star_data[good][flux_column_name][np.isnan(comp_fluxes[flux_column_name])])
+    # Convert comp_fluxes to a regular Table, not a QTable, to work around
+    # https://github.com/astropy/astropy/issues/10944
+    # in which it was reported that QTable columns with units cannot be aggregated.
+
+    comp_fluxes = Table(comp_fluxes)
 
     # Check whether any of the columns are masked, but with no masked values,
     # and convert to regular column...eventually
@@ -138,7 +140,7 @@ def calc_aij_relative_flux(
     comp_errors = comp_fluxes.groups.aggregate(add_in_quadrature)[error_column_name]
 
     comp_total_vector = np.ones_like(star_data[counts_column_name])
-    comp_error_vector = np.ones_like(star_data[counts_column_name])
+    comp_error_vector = np.ones_like(star_data[error_column_name])
 
     if len(set(comp_num_stars)) > 1:
         raise RuntimeError("Different number of stars in comparison sets")
@@ -147,17 +149,22 @@ def calc_aij_relative_flux(
 
     # Have to remove the flux of the star if the star is a comparison
     # star.
-    is_comp = np.zeros_like(star_data[counts_column_name])
+    # Use the .value below so that we can set the array to 1 and multiply
+    # by it without affecting units of the result.
+    is_comp = np.zeros_like(star_data[counts_column_name]).value
     is_comp[good] = 1
     flux_offset = -star_data[counts_column_name] * is_comp
 
+    # Convert comp_fluxes back to a QTable and redo groups
+    comp_fluxes = QTable(comp_fluxes)
+    comp_fluxes = comp_fluxes.group_by("date-obs")
     # This seems a little hacky; there must be a better way
     for date_obs, comp_total, comp_error in zip(
         comp_fluxes.groups.keys, comp_totals, comp_errors, strict=True
     ):
         this_time = star_data["date-obs"] == date_obs[0]
         comp_total_vector[this_time] *= comp_total
-        comp_error_vector[this_time] = comp_error
+        comp_error_vector[this_time] = comp_error * comp_fluxes[error_column_name].unit
 
     relative_flux = star_data[counts_column_name] / (comp_total_vector + flux_offset)
     relative_flux = relative_flux.flatten()


Test patch:
diff --git a/stellarphot/differential_photometry/tests/test_aij_rel_fluxes.py b/stellarphot/differential_photometry/tests/test_aij_rel_fluxes.py
index 6afd173..fb8dd39 100644
--- a/stellarphot/differential_photometry/tests/test_aij_rel_fluxes.py
+++ b/stellarphot/differential_photometry/tests/test_aij_rel_fluxes.py
@@ -29,8 +29,8 @@ def _raw_photometry_table():
     # and four stars
     star_ra = 250.0 * u.degree + np.arange(n_stars) * 10 * u.arcmin
     star_dec = np.array([45.0] * n_stars) * u.degree
-    fluxes = np.array([10000.0, 20000, 30000, 40000])
-    errors = np.sqrt(fluxes) + 50
+    fluxes = np.array([10000.0, 20000, 30000, 40000]) * u.adu
+    errors = (np.sqrt(fluxes.value) + 50) * u.electron
     star_ids = np.arange(1, 5, dtype="int")
 
     # Stars 2, 3 and 4 will be the comparison stars
@@ -65,10 +65,21 @@ def _raw_photometry_table():
             "noise_electrons",
             "star_id",
         ],
+        units=[
+            None,
+            u.degree,
+            u.degree,
+            u.adu,
+            u.electron,
+            None,
+        ],
     )
 
-    _ = PhotometryData(raw_table)
-    return expected_flux_ratios, expected_flux_error, raw_table, raw_table[1:4]
+    photom = PhotometryData(raw_table)
+    # MAKE SURE to return photom, not raw_table, below to trigger the bug
+    # https://github.com/feder-observatory/stellarphot/issues/421
+    # in which, it turns out, QTable columns with units cannot be aggregated.
+    return expected_flux_ratios, expected_flux_error, photom, photom[1:4]
 
 
 @pytest.mark.parametrize("comp_ra_dec_have_units", [True, False])
@@ -77,6 +88,8 @@ def _raw_photometry_table():
 def test_relative_flux_calculation(
     in_place, star_ra_dec_have_units, comp_ra_dec_have_units
 ):
+    # In addition to checking the flux calculation values, this is also a regression
+    # test for #421.
     expected_flux, expected_error, input_table, comp_star = _raw_photometry_table()
 
     # Try doing it all at once
@@ -95,7 +108,7 @@ def test_relative_flux_calculation(
     output_table = calc_aij_relative_flux(input_table, comp_star, in_place=in_place)
     output_flux = output_table["relative_flux"]
     output_error = output_table["relative_flux_error"]
-    print(all_expected_flux - output_flux)
+
     np.testing.assert_allclose(output_flux, all_expected_flux)
     np.testing.assert_allclose(output_error, all_expected_error)
     if in_place:
@@ -122,8 +135,7 @@ def test_bad_comp_star(bad_thing):
             ra=last_one["ra"][0], dec=last_one["dec"][0], unit=u.degree
         )
         coord_bad_ra = coord_inp.ra + 3 * u.arcsecond
-        print(len(last_one), coord_inp)
-        input_table["ra"][-1] = coord_bad_ra.degree
+        input_table["ra"][-1] = coord_bad_ra
     elif bad_thing == "NaN":
         input_table["aperture_net_cnts"][-1] = np.nan