JSKenyon
commented
1 year ago I have managed to reproduce this (inconsistently) outside QC (still calling uncontroversial QC internals though):
from quartical.data_handling.angles import _make_parangles
from concurrent.futures import ThreadPoolExecutor, wait, ProcessPoolExecutor
from daskms import xds_from_storage_table, xds_from_storage_ms
import numpy as np
if __name__ == "__main__":
ms_path = "~/reductions/3C147/msdir/C147_unflagged.MS"
data_xds_list = xds_from_storage_ms(
ms_path,
group_cols=("FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER")
)
anttab = xds_from_storage_table(ms_path + "::ANTENNA")[0]
feedtab = xds_from_storage_table(ms_path + "::FEED")[0]
fieldtab = xds_from_storage_table(ms_path + "::FIELD")[0]
# We do this eagerly to make life easier.
feeds = feedtab.POLARIZATION_TYPE.values
unique_feeds = np.unique(feeds)
if np.all([feed in "XxYy" for feed in unique_feeds]):
feed_type = "linear"
elif np.all([feed in "LlRr" for feed in unique_feeds]):
feed_type = "circular"
else:
raise ValueError("Unsupported feed type/configuration.")
phase_dirs = fieldtab.PHASE_DIR.values
field_centres = [phase_dirs[0, 0] for xds in data_xds_list]
# TODO: This could be more complicated for arrays with multiple feeds.
receptor_angles = feedtab.RECEPTOR_ANGLE.values
ant_names = anttab.NAME.values
ant_positions_ecef = anttab.POSITION.values # ECEF coordinates.
epoch = "J2000" # TODO: Should this be configurable?
futures = []
with ThreadPoolExecutor(max_workers=6) as tp:
for xds, field_centre in zip(data_xds_list[:10], field_centres[:10]):
futures.append(
tp.submit(
_make_parangles,
xds.TIME.values,
ant_names,
ant_positions_ecef,
receptor_angles,
field_centre,
epoch
)
)
wait(futures)
results1 = [f.result() for f in futures]
futures = []
with ThreadPoolExecutor(max_workers=6) as tp:
for xds, field_centre in zip(data_xds_list[:10], field_centres[:10]):
futures.append(
tp.submit(
_make_parangles,
xds.TIME.values,
ant_names,
ant_positions_ecef,
receptor_angles,
field_centre,
epoch
)
)
wait(futures)
results2 = [f.result() for f in futures]
for i in range(10):
print(np.all(results1[i] == results2[i]))
print(np.max(results1[i] - results2[i]))Note that the above only gives incorrect results some of the time. Will dig a little further.
sjperkins
bennahugo
Describe the bug In a multi-threaded regime, the parallactic angles may be slightly wrong.
To Reproduce Run QuartiCal with
input_model.apply_p_jones=1and several datasets. The model will not be the same on every run. The bug is intermittent and inconsistent, similar to a race or out of bounds error. However, it is absent when no parallactic angles are applied. Note that the error will be very, very small but seems to occasionally have knock on effects.Expected behavior The parallactic angle computation should be deterministic - the same for given inputs.
Version
main