diff --git a/tests/create_reference_file.py b/tests/create_reference_file.py new file mode 100644 index 0000000000000000000000000000000000000000..f62021e5e84d8fe736485b26745ab913ec9c28bb --- /dev/null +++ b/tests/create_reference_file.py @@ -0,0 +1,240 @@ +"""Create reference files.""" + +import logging +import os + +import numpy as np +import xarray as xr +from rich.logging import RichHandler +from ruamel.yaml import YAML + +import mvcm.spectral_tests as tst + +_LOG_FORMAT = "%(message)s" +logging.basicConfig(level="NOTSET", datefmt="[%X]", format=_LOG_FORMAT, handlers=[RichHandler()]) +logger = logging.getLogger(__name__) + +_scene_list = [ + "Land_Day", + "Land_Day_Coast", + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + "Polar_Day_Land", + "Polar_Day_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + "Polar_Day_Snow", + "Land_Night", + "Polar_Night_Land", + "Polar_Night_Snow", + "Day_Snow", + "Night_Snow", +] + + +_threshold_file = "/home/pveglio/mvcm/tests/fixtures/thresholds.mvcm.snpp.v0.0.1.yaml" +_test_path = "/ships19/hercules/pveglio/mvcm_git_tests" + + +def main() -> None: + """Create reference file.""" + with open(_threshold_file) as f: + text = f.read() + thresholds = YAML(typ="safe").load(text) + + output_dict: dict = {} + + for scene_name in _scene_list: + logger.info(f"Running tests for {scene_name} \n") + + if os.path.isfile(f"{_test_path}/ref_confidence_{scene_name}.nc"): + logger.info(f"Skipping {scene_name}, file already exists") + continue + + if not os.path.isfile(f"{_test_path}/test_scene_{scene_name}.nc"): + logger.info(f"Skipping {scene_name}, test file not present in {_test_path}") + continue + + viirs_data = xr.open_dataset(f"{_test_path}/test_scene_{scene_name}.nc") + + bits = { + "test": np.zeros(viirs_data.M15.shape, dtype=np.int8), + "qa": np.zeros(viirs_data.M15.shape, dtype=np.int8), + } + + if np.all(viirs_data[scene_name].values == 0): + logger.info("Skipping, no pixels in scene.") + continue + + my_scene = tst.CloudTests(data=viirs_data, scene_name=scene_name, thresholds=thresholds) # type: ignore + + # 11um Test + logger.info(f"Running 11um test for {scene_name}") + confidence, bits = my_scene.test_11um("M15", np.ones(viirs_data.M15.shape), bits) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "11um" + ) + + # surface_temperature_test() is not performed right now + + # SST Test + logger.info(f"Running SST test for {scene_name} \n") + confidence, bits = my_scene.sst_test("M15", "M16", np.ones(viirs_data.M15.shape), bits) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "sst" + ) + + # 8.6 - 11 um Test + logger.info(f"Running 8.6 - 11 um test for {scene_name} \n") + confidence, bits = my_scene.bt_diff_86_11um("M14-M15", np.ones(viirs_data.M15.shape), bits) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "86_11um_diff" + ) + + # 11 - 12um Test + logger.info(f"Running 11 - 12um test for {scene_name} \n") + confidence, bits = my_scene.test_11_12um_diff( + "M15-M16", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "11_12um_diff" + ) + + # 11um Variability Test + logger.info(f"Running 11um variability test for {scene_name} \n") + confidence, bits = my_scene.variability_11um_test( + "M15", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, + confidence, + bits, + output_dict, + "variability_11um", + ) + + # 11 - 4um Test Ocean + logger.info(f"Running 11 - 4um test for {scene_name} \n") + confidence, bits = my_scene.bt_difference_11_4um_test_ocean( + "M15-M12", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "11_4um_ocean" + ) + + # 11 - 4um Test Land + # confidence, bits = my_scene.bt_difference_11_4um_test_land( + # "M15-M12", np.ones(viirs_data.M15.shape), bits + # ) + # bits, output_dict = update_and_clean( + # viirs_data.M15.shape, confidence, bits, output_dict, "11_4um_land" + # ) + + # 11 - 4um Oceanic Stratus Test + logger.info(f"Running 11 - 4um Oceanic Stratus test for {scene_name} \n") + confidence, bits = my_scene.oceanic_stratus_11_4um_test( + "M15-M12", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, + confidence, + bits, + output_dict, + "oceanic_stratus_11_4um", + ) + + # NIR Reflectance Test + logger.info(f"Running nir reflectance test for {scene_name} \n") + confidence, bits = my_scene.nir_reflectance_test("M07", np.ones(viirs_data.M15.shape), bits) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, + confidence, + bits, + output_dict, + "nir_reflectance", + ) + + # Vis/NIR Ratio test + logger.info(f"Running vis/nir ratio test for {scene_name} \n") + confidence, bits = my_scene.vis_nir_ratio_test( + "M07-M05ratio", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "vis_nir_ratio" + ) + + # 1.6 - 2.1um Test + logger.info(f"Running 1.6 - 2.1um test for {scene_name} \n") + confidence, bits = my_scene.test_16_21um_reflectance( + "M10", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, + confidence, + bits, + output_dict, + "16_21um_reflectance", + ) + + # Visible Reflectance Test + logger.info(f"Running visible reflectance test for {scene_name} \n") + confidence, bits = my_scene.visible_reflectance_test( + "M05", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "vis_refl" + ) + + # GEMI Test + logger.info(f"Running GEMI test for {scene_name} \n") + confidence, bits = my_scene.gemi_test("GEMI", np.ones(viirs_data.M15.shape), bits) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "gemi" + ) + + # 1.38 High Clouds Test + logger.info(f"Running 1.38um high clouds test for {scene_name} \n") + confidence, bits = my_scene.test_1_38um_high_clouds( + "M09", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "high_clouds" + ) + + # 4 -12um Thin Cirrus Test + logger.info(f"Running 4-12um thin cirrus test for {scene_name} \n") + confidence, bits = my_scene.thin_cirrus_4_12um_BTD_test( + "M13-M16", np.ones(viirs_data.M15.shape), bits + ) + bits, output_dict = update_and_clean( + viirs_data.M15.shape, confidence, bits, output_dict, "thin_cirrus" + ) + + logger.info(f"Saving reference file for scene {scene_name} \n") + + output = xr.Dataset.from_dict(output_dict) + output.to_netcdf(f"{_test_path}/ref_confidence_{scene_name}.nc") + + +def update_and_clean( + shape: tuple, confidence: np.ndarray, bits: dict, output_dict: dict, test_name: str +) -> tuple[dict, dict]: + """Update and clean output dictionary.""" + output_dict[f"confidence_{test_name}"] = {"dims": ("x", "y"), "data": confidence} + output_dict[f"qa_bit_{test_name}"] = {"dims": ("x", "y"), "data": bits["qa"]} + output_dict[f"test_bit_{test_name}"] = {"dims": ("x", "y"), "data": bits["test"]} + + bits = { + "test": np.zeros(shape, dtype=np.int8), + "qa": np.zeros(shape, dtype=np.int8), + } + + return bits, output_dict + + +if __name__ == "__main__": + main() diff --git a/tests/create_test_file.py b/tests/create_test_file.py index 65c2ddc00d38078dc18ab31f3e2755477b75c98d..183c3c37e07e2f5ac4005500250d6be4b915ce75 100644 --- a/tests/create_test_file.py +++ b/tests/create_test_file.py @@ -1,19 +1,26 @@ """Create test files.""" -# import os.path -# from glob import glob + +import logging from os.path import exists import numpy as np import numpy.typing as npt import xarray as xr +from rich.logging import RichHandler +from ruamel.yaml import YAML import mvcm.read_data as rd +_LOG_FORMAT = "%(message)s" +logging.basicConfig(level="NOTSET", datefmt="[%X]", format=_LOG_FORMAT, handlers=[RichHandler()]) +logger = logging.getLogger(__name__) + # #################################################################### # # TEST CASES # data: _datapath = "/ships19/hercules/pveglio/mvcm_viirs_hires" _outpath = "/ships19/hercules/pveglio/mvcm_git_tests" +_thresholds_file = f"{_outpath}/thresholds.mvcm.snpp.v0.0.1.yaml" # The following sets of files try to cover all the possible scenes found # in the mvcm. the list of scenes per granule is reported on top of every @@ -26,8 +33,8 @@ _outpath = "/ships19/hercules/pveglio/mvcm_git_tests" # - Land_Day_Desert # - Land_Day # - Ocean_Day -_fname_mod02 = f"{_datapath}/VNP02MOD.A2022173.1312.001.2022174011547.uwssec_bowtie_corrected.nc" -_fname_mod03 = f"{_datapath}/VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc" +# _fname_mod02 = f"{_datapath}/VNP02MOD.A2022173.1312.001.2022174011547.uwssec_bowtie_corrected.nc" +# _fname_mod03 = f"{_datapath}/VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc" # _fname_img02 = f"{_datapath}/VNP02IMG.A2022173.1312.001.2022174011547.uwssec_bowtie_corrected.nc" # _fname_img03 = f"{_datapath}/VNP03IMG.A2022173.1312.001.2022174012746.uwssec.nc" @@ -37,8 +44,8 @@ _fname_mod03 = f"{_datapath}/VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc" # - Polar_Day_Desert_Coast # - Polar_Day_Desert # - Polar_Day_Ocean -# _fname_mod02 = f"{_datapath}/VNP02MOD.A2022173.1324.001.2022174014257.uwssec_bowtie_corrected.nc" -# _fname_mod03 = f"{_datapath}/VNP03MOD.A2022173.1324.001.2022174012743.uwssec.nc" +_fname_mod02 = f"{_datapath}/VNP02MOD.A2022173.1324.001.2022174014257.uwssec_bowtie_corrected.nc" +_fname_mod03 = f"{_datapath}/VNP03MOD.A2022173.1324.001.2022174012743.uwssec.nc" # _fname_img02 = f"{_datapath}/VNP02IMG.A2022173.1324.001.2022174014257.uwssec_bowtie_corrected.nc" # _fname_img03 = f"{_datapath}/VNP03IMG.A2022173.1324.001.2022174012743.uwssec.nc" @@ -51,14 +58,14 @@ _fname_mod03 = f"{_datapath}/VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc" _out_fname = f"{_outpath}/test_scene" # ancillary files: -_geos_atm_1 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4" -_geos_atm_2 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4" -_geos_land = "GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4" -_geos_ocean = "GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4" -_geos_constants = "GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4" -_ndvi_file = "NDVI.FM.c004.v2.0.WS.00-04.177.hdf" -_sst_file = "oisst.20220622" -_eco_file = "goge1_2_img.v1" +_geos_atm_1 = f"{_datapath}/GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4" +_geos_atm_2 = f"{_datapath}/GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4" +_geos_land = f"{_datapath}/GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4" +_geos_ocean = f"{_datapath}/GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4" +_geos_constants = f"{_datapath}/GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4" +_ndvi_file = f"{_datapath}/NDVI.FM.c004.v2.0.WS.00-04.177.hdf" +_sst_file = f"{_datapath}/oisst.20220622" +_eco_file = f"{_datapath}/goge1_2_img.v1" # #################################################################### # @@ -127,6 +134,7 @@ def write_test_file() -> None: def main( mod02: str, mod03: str, + thresholds_file: str, geos_atm_1: str, geos_atm_2: str, geos_land: str, @@ -155,12 +163,15 @@ def main( "ANC_DIR": f"{data_path}/ancillary", } - viirs_data = rd.get_data("snpp", "viirs", file_names, hires=False) + with open(thresholds_file) as f: + text = f.read() + thresholds = YAML(typ="safe").load(text) + viirs_data = rd.get_data("snpp", "viirs", file_names, thresholds, hires=False) out_ds = xr.Dataset() for scene in _scene_list: - print(f"scene: {scene}") + logger.info(f"scene: {scene}") if exists(f"{out_fname}_{scene}.nc"): continue if len(viirs_data[scene].values[viirs_data[scene].values == 1]) < _default_scene_size: @@ -189,6 +200,8 @@ def main( for var in viirs_data: sample_data[var] = viirs_data[var][scn_idx[0], scn_idx[1]] + sample_data.attrs["scene_name"] = scene + out_ds = sample_data out_ds.to_netcdf(f"{out_fname}_{scene}.nc") @@ -197,6 +210,7 @@ if __name__ == "__main__": main( mod02=_fname_mod02, mod03=_fname_mod03, + thresholds_file=_thresholds_file, geos_atm_1=_geos_atm_1, geos_atm_2=_geos_atm_2, geos_land=_geos_land, diff --git a/tests/fixtures/thresholds.mvcm.snpp.v0.0.1.yaml b/tests/fixtures/thresholds.mvcm.snpp.v0.0.1.yaml index 99cc058820454e00b280e11c51a5e159d636eaa6..a3451722ea59da972a973d75db6c23e9e51ff3fe 100644 --- a/tests/fixtures/thresholds.mvcm.snpp.v0.0.1.yaml +++ b/tests/fixtures/thresholds.mvcm.snpp.v0.0.1.yaml @@ -270,10 +270,10 @@ Polar_Day_Land: coeffs: [3.0, 1.0] cmult: 0.3 adj: 0.3 - perform: True + perform: False 11-4um_Oceanic_Stratus_Test: thr: [-16.0, -14.0, -12.0, 1.0, 1.0] - perform: True + perform: False pdlh20: [215.0, 220.0, 225.0, 1.0, 1.0] Visible_Reflectance_Test: thr: [0.207, 0.169, 0.132, 1.0, 1.0] @@ -339,10 +339,10 @@ Polar_Day_Coast: coeffs: [3.0, 1.0] cmult: 0 # I NEED TO WORK ON THIS adj: 0 # I NEED TO WORK ON THIS - perform: True + perform: False 11-4um_Oceanic_Stratus_Test: thr: [-16.0, -14.0, -12.0, 1.0, 1.0] - perform: True + perform: False Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0] Visible_Reflectance_Test: thr: [0.207, 0.169, 0.132, 1.0, 1.0] @@ -358,11 +358,11 @@ Polar_Day_Desert: coeffs: [3.5, 1.0] cmult: 0 # I NEED TO WORK ON THIS adj: 0 # I NEED TO WORK ON THIS - perform: True + perform: False 11-4um_Oceanic_Stratus_Test: # thr: [-22.0, -20.0, -18.0, -2.0, 0.0, 2.0, 1.0, 1.0] thr: [2.0, 0.0, -2.0, -18.0, -20.0, -22.0, 1.0, 1.0] - perform: True # I disable this test while I figure out if I need + perform: False # I disable this test while I figure out if I need bt_cutoff: 320.0 # the bt_cutoff parameter Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0] Visible_Reflectance_Test: @@ -384,11 +384,11 @@ Polar_Day_Desert_Coast: coeffs: [3.5, 1.0] cmult: 0 # I NEED TO WORK ON THIS adj: 0 # I NEED TO WORK ON THIS - perform: True + perform: False 11-4um_Oceanic_Stratus_Test: # thr: [-23.0, -21.0, -19.0, -2.0, 0.0, 2.0, 1.0, 1.0] thr: [2.0, 0.0, -2.0, -19.0, -21.0, -23.0, 1.0, 1.0] - perform: True # see Polar_Day_Desert + perform: False # see Polar_Day_Desert bt_cutoff: 320.0 Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0] Visible_Reflectance_Test: @@ -405,21 +405,22 @@ Polar_Day_Snow: coeffs: [3.0, 1.0] cmult: 0 # I NEED TO WORK ON THIS adj: 0 # I NEED TO WORK ON THIS - perform: True - dpsbt1: 230.0 + perform: False Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0] 1.38um_High_Cloud_Test: thr: [0.060, 0.0525, 0.045, 1.0, 1.0] tpw: 0.75 perform: True - 11-4um_BT_Difference_Test_Land: + # 11-4um_BT_Difference_Test_Land: + 11-4um_Oceanic_Stratus_Test: low: [20.0, 18.0, 16.0, 1.0] mid1: [18.0, 16.0, 2.0, 1.0] mid2: [0.0, 0.0, 0.0, 0.0] mid3: [0.0, 0.0, 0.0, 0.0] high: [18.0, 16.0, 14.0, 1.0] - bt_11_bounds: [230.0, 0.0, 0.0, 245.0] - perform: True + bt11_bounds: [230.0, 0.0, 0.0, 245.0] + bt_cutoff: 230.0 + perform: False dps11_12adj: 0.8 Polar_Night_Snow: @@ -435,7 +436,7 @@ Polar_Night_Snow: mid2: [0.00, 0.00, 0.00, 0.0] # pn_4_12m2 mid3: [0.00, 0.00, 0.00, 0.0] # pn_4_12m3 high: [2.50, 2.00, 1.50, 1.0] # pn_4_12h - bt11_bounds: [235.0, 0.0, 0.0, 265.0] # bt_11_bounds + bt11_bounds: [235.0, 0.0, 0.0, 265.0] # bt11_bounds perform: True 7.3-11um_BTD_Mid_Level_Cloud_Test: low: [-1.00, 0.00, 1.00, 1.0, 1.0] # pn_7_11l @@ -479,12 +480,12 @@ Polar_Day_Ocean: perform: True 8.6-11um_Test: thr: [-0.50, -1.00, -1.50, 1.0, 1.0] - perform: True + perform: False 11-12um_Cirrus_Test: coeffs: [3.0, 1.0] cmult: 0.3 adj: 1.25 - perform: True + perform: False NIR_Reflectance_Test: thr: [0.062, 0.043, 0.029, 1.0, 1.0] coeffs: [1.7291, 0.0715, -0.0026, 0.000025889] @@ -514,7 +515,7 @@ Polar_Day_Ocean: perform: True 11-4um_Oceanic_Stratus_Test: thr: [-11.0, -9.0, -7.0, 1.0, 1.0] - perform: True + perform: False Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0] pdo_b26pfm: 1.0 pdo_b2bias_adj: 0.96 @@ -566,7 +567,7 @@ Day_Snow: 11-12um_Cirrus_Test: coeffs: [3.0, 1.0] cmult: 0.3 - adj: 0 # I NEED TO WORK ON THIS + adj: 0.8 perform: True CO2_High_Clouds_Test: [222.0, 224.0, 226.0, 1.0, 1.0] Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0] @@ -574,13 +575,22 @@ Day_Snow: thr: [0.060, 0.0525, 0.045, 1.0, 1.0] tpw: 0.75 perform: True - 11-4um_BT_Difference_Test_Land: + # 11-4um_BT_Difference_Test_Land: + # low: [20.0, 18.0, 16.0, 1.0] + # mid1: [18.0, 16.0, 2.0, 1.0] + # mid2: [0.0, 0.0, 0.0, 0.0] + # mid3: [0.0, 0.0, 0.0, 0.0] + # high: [18.0, 16.0, 14.0, 1.0] + # bt11_bounds: [230.0, 0.0, 0.0, 245.0] + # perform: True + 11-4um_Oceanic_Stratus_Test: low: [20.0, 18.0, 16.0, 1.0] mid1: [18.0, 16.0, 2.0, 1.0] mid2: [0.0, 0.0, 0.0, 0.0] mid3: [0.0, 0.0, 0.0, 0.0] high: [18.0, 16.0, 14.0, 1.0] - bt_11_bounds: [230.0, 0.0, 0.0, 245.0] + bt11_bounds: [230.0, 0.0, 0.0, 245.0] + bt_cutoff: 230.0 perform: True ds11_12adj: 0.8 ds11_12lcmult: 0.3 @@ -605,7 +615,7 @@ Night_Snow: mid2: [0.0, 0.0, 0.0, 0.0] mid3: [0.0, 0.0, 0.0, 0.0] high: [18.0, 16.0, 14.0, 1.0] - bt_11_bounds: [230.0, 0.0, 0.0, 245.0] + bt11_bounds: [230.0, 0.0, 0.0, 245.0] bt_thr: 230.0 tpw_thr: 0.2 perform: True @@ -679,18 +689,18 @@ Night_Snow_Inversion: n65_11: [10.0, 1.0] Snow_Mask: - sm_bt11: 280.0 - sm_ndsi: 0.31 - Grnlnd_ndsi: 0.66 - sm_ref2: 0.1056 - sm85_11: 0.0 - sm85_11hel: 1.5 - sm37_11: 11.0 - sm37_11hel: 16.0 - sm_mnir: 0.196 - sm_lsfcdif: 20.0 - sm_wsfcdif: 20.0 - sm_bt1: 273.0 + bt11: 280.0 + ndsi: 0.31 + greenland_ndsi: 0.66 + ref2: 0.1056 + diff85_11: 0.0 + diff85_11hel: 1.5 + diff37_11: 11.0 + diff37_11hel: 16.0 + mnir: 0.196 + lsfcdif: 20.0 + wsfcdif: 20.0 + bt1: 273.0 prd_ndvi_const: -0.2015 Coastal_NDVI_Thresholds: diff --git a/tests/test_conf.py b/tests/test_conf.py index ad6df62b8795b53c89d5f07ac3ea23b8264e65c6..6da43a3ffb97412f67dc2969fd86d7e9e538358d 100644 --- a/tests/test_conf.py +++ b/tests/test_conf.py @@ -70,19 +70,21 @@ def test_c_single_threshold(rad, single_threshold, reference_data): ref_confidence = np.load(reference_data)["ref_sgl"] ref_confidence_flipped = np.load(reference_data)["ref_sgl_flipped"] - locut = np.array(single_threshold[0], dtype=float) - midpt = np.array(single_threshold[1], dtype=float) - hicut = np.array(single_threshold[2], dtype=float) - power = np.array(single_threshold[3], dtype=float) + locut = np.full(np.shape(rad), single_threshold[0], dtype=np.float32) + midpt = np.full(np.shape(rad), single_threshold[1], dtype=np.float32) + hicut = np.full(np.shape(rad), single_threshold[2], dtype=np.float32) + power = np.full(np.shape(rad), single_threshold[3], dtype=np.float32) c = anc.py_conf_test(np.array(rad, dtype=np.float32), locut, hicut, power, midpt) - assert np.all(c == ref_confidence) + # NOTE: I have to use allclose() because the ref_confidence has been computed with doouble + # precision and the C version returns single precision, which has some rounding differences + assert np.allclose(c, ref_confidence) single_threshold[0:-1] = single_threshold[-2::-1] - locut = np.array(single_threshold[0], dtype=float) - midpt = np.array(single_threshold[1], dtype=float) - hicut = np.array(single_threshold[2], dtype=float) - power = np.array(single_threshold[3], dtype=float) + locut = np.full(np.shape(rad), single_threshold[0], dtype=np.float32) + midpt = np.full(np.shape(rad), single_threshold[1], dtype=np.float32) + hicut = np.full(np.shape(rad), single_threshold[2], dtype=np.float32) + power = np.full(np.shape(rad), single_threshold[3], dtype=np.float32) c = anc.py_conf_test(np.array(rad, dtype=np.float32), locut, hicut, power, midpt) - assert np.all(c == ref_confidence_flipped) + assert np.allclose(c, ref_confidence_flipped)