Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix: allow to use numpy 2 and pyuvdata 3 #959

Merged
merged 19 commits into from
Jul 16, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ jobs:
strategy:
matrix:
os: [ubuntu-latest, macos-latest]
python-version: [3.9, "3.10", "3.11", "3.12"]
python-version: ["3.10", "3.11", "3.12"]
fail-fast: false

steps:
Expand Down Expand Up @@ -53,7 +53,7 @@ jobs:
strategy:
matrix:
os: [ubuntu-latest]
python-version: ["3.11"]
python-version: ["3.12"]
fail-fast: false

steps:
Expand Down
37 changes: 23 additions & 14 deletions hera_cal/abscal.py
Original file line number Diff line number Diff line change
Expand Up @@ -363,6 +363,8 @@ def TT_phs_logcal(model, data, antpos, wgts=None, refant=None, assume_2D=True,
# angle after divide
with warnings.catch_warnings():
warnings.filterwarnings("ignore", message="invalid value encountered in divide")
warnings.filterwarnings("ignore", message="divide by zero encountered in divide")

ydata = {k: np.angle(data[k] / model[k]) for k in keys}

# make unit weights if None
Expand Down Expand Up @@ -476,6 +478,7 @@ def amp_logcal(model, data, wgts=None, verbose=True):
# difference of log-amplitudes is ydata independent variable
with warnings.catch_warnings():
warnings.filterwarnings("ignore", message="invalid value encountered in divide")
warnings.filterwarnings("ignore", message="divide by zero encountered in divide")
ydata = odict([(k, np.log(np.abs(data[k] / model[k]))) for k in keys])

# make weights if None
Expand Down Expand Up @@ -547,6 +550,7 @@ def phs_logcal(model, data, wgts=None, refant=None, verbose=True):
# angle of visibility ratio is ydata independent variable
with warnings.catch_warnings():
warnings.filterwarnings("ignore", message="invalid value encountered in divide")
warnings.filterwarnings("ignore", message="divide by zero encountered in divide")
ydata = odict([(k, np.angle(data[k] / model[k])) for k in keys])

# make weights if None
Expand Down Expand Up @@ -2623,7 +2627,7 @@ def __init__(self, model, data, refant=None, wgts=None, antpos=None, freqs=None,
#!/usr/bin/env python
uvd = pyuvdata.UVData()
uvd.read_miriad(<filename>)
antenna_pos, ants = uvd.get_ENU_antpos()
antenna_pos, ants = uvd.telescope.get_enu_antpos()
antpos = dict(zip(ants, antenna_pos))
----
This is needed only for Tip Tilt, phase slope, and delay slope calibration.
Expand Down Expand Up @@ -3554,8 +3558,11 @@ def get_all_times_and_lsts(hd, solar_horizon=90.0, unwrap=True):

# remove times when sun was too high
if solar_horizon < 90.0:
lat, lon, alt = hd.telescope_location_lat_lon_alt_degrees
solar_alts = utils.get_sun_alt(all_times, latitude=lat, longitude=lon)
solar_alts = utils.get_sun_alt(
all_times,
latitude=hd.telescope.location.lat.deg,
longitude=hd.telescope.location.lon.deg
)
solar_flagged = solar_alts > solar_horizon
return all_times[~solar_flagged], all_lsts[~solar_flagged]
else: # skip this step for speed
Expand Down Expand Up @@ -4051,8 +4058,8 @@ def post_redcal_abscal_run(data_file, redcal_file, model_files, raw_auto_file=No
warnings.warn(f"Warning: Overwriting redcal gain_scale of {hc.gain_scale} with model gain_scale of {hdm.vis_units}", RuntimeWarning)
hc.gain_scale = hdm.vis_units # set vis_units of hera_cal based on model files.
hd_autos = io.HERAData(raw_auto_file)
assert hdm.x_orientation.lower() == hd.x_orientation.lower(), 'Data x_orientation, {}, does not match model x_orientation, {}'.format(hd.x_orientation.lower(), hdm.x_orientation.lower())
assert hc.x_orientation.lower() == hd.x_orientation.lower(), 'Data x_orientation, {}, does not match redcal x_orientation, {}'.format(hd.x_orientation.lower(), hc.x_orientation.lower())
assert hdm.telescope.x_orientation.lower() == hd.telescope.x_orientation.lower(), 'Data x_orientation, {}, does not match model x_orientation, {}'.format(hd.telescope.x_orientation.lower(), hdm.telescope.x_orientation.lower())
assert hc.telescope.x_orientation.lower() == hd.telescope.x_orientation.lower(), 'Data x_orientation, {}, does not match redcal x_orientation, {}'.format(hd.telescope.x_orientation.lower(), hc.telescope.x_orientation.lower())
pol_load_list = [pol for pol in hd.pols if split_pol(pol)[0] == split_pol(pol)[1]]

# get model bls and antpos to use later in baseline matching
Expand Down Expand Up @@ -4110,7 +4117,7 @@ def post_redcal_abscal_run(data_file, redcal_file, model_files, raw_auto_file=No
model_flags.select_or_expand_times(model_times_to_load)
model_blvecs = {bl: model.antpos[bl[0]] - model.antpos[bl[1]] for bl in model.keys()}
utils.lst_rephase(model, model_blvecs, model.freqs, data.lsts - model.lsts,
lat=hdm.telescope_location_lat_lon_alt_degrees[0], inplace=True)
lat=hdm.telescope.location.lat.deg, inplace=True)

# Flag frequencies and times in the data that are entirely flagged in the model
model_flag_waterfall = np.all([f for f in model_flags.values()], axis=0)
Expand Down Expand Up @@ -4346,19 +4353,21 @@ def run_model_based_calibration(data_file, model_file, output_filename, auto_fil
refant_init = str(refant)
# initialize HERACal to store new cal solutions
hc = UVCal()
hc = hc.initialize_from_uvdata(uvdata=hdd, gain_convention='divide', cal_style='sky',
ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
metadata_only=False, cal_type='gain',
future_array_shapes=True)
hc = hc.initialize_from_uvdata(
uvdata=hdd, gain_convention='divide', cal_style='sky',
ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
metadata_only=False, cal_type='gain',
)

hc = io.to_HERACal(hc)
hc.update(flags=data_ant_flags)
# generate cal object from model to hold model flags.
hcm = UVCal()
hcm = hcm.initialize_from_uvdata(uvdata=hdm, gain_convention='divide', cal_style='sky',
ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
metadata_only=False, cal_type='gain',
future_array_shapes=True)
hcm = hcm.initialize_from_uvdata(
uvdata=hdm, gain_convention='divide', cal_style='sky',
ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
metadata_only=False, cal_type='gain',
)
hcm = io.to_HERACal(hcm)
hcm.update(flags=model_ant_flags)
# init all gains to unity.
Expand Down
2 changes: 1 addition & 1 deletion hera_cal/chunker.py
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ def chunk_files(filenames, inputfile, outputfile, chunk_size, type="data", # no
if spw_range is None:
spw_range = (0, chunked_files.Nfreqs)
# convert polarizations to jones integers.
jones = [uvutils.polstr2num(pol, x_orientation=chunked_files.x_orientation) for pol in polarizations]
jones = [uvutils.polstr2num(pol, x_orientation=chunked_files.telescope.x_orientation) for pol in polarizations]
chunked_files.select(freq_chans=np.arange(spw_range[0], spw_range[1]).astype(int), jones=jones)
# throw away fully flagged baselines.
if throw_away_flagged_ants:
Expand Down
2 changes: 1 addition & 1 deletion hera_cal/flag_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ def solar_flag(flags, times=None, flag_alt=0.0, longitude=21.42830, latitude=-30
elif isinstance(flags, UVData):
if verbose:
print("Note: using latitude and longitude in given UVData object")
latitude, longitude, altitude = flags.telescope_location_lat_lon_alt_degrees
latitude, longitude, altitude = flags.telescope._location.lat_lon_alt_degrees()
times = np.unique(flags.time_array)
dtype = 'uvd'
if dtype in ['ndarr', 'DC']:
Expand Down
26 changes: 13 additions & 13 deletions hera_cal/frf.py
Original file line number Diff line number Diff line change
Expand Up @@ -285,8 +285,9 @@ def sky_frates(uvd, keys=None, frate_standoff=0.0, frate_width_multiplier=1.0,
"""
if keys is None:
keys = uvd.get_antpairpols()
antpos, antnums = uvd.get_ENU_antpos()
sinlat = np.sin(np.abs(uvd.telescope_location_lat_lon_alt[0]))
antpos = uvd.telescope.get_enu_antpos()
antnums = uvd.telescope.antenna_numbers
sinlat = np.sin(np.abs(uvd.telescope.location.lat))
frate_centers = {}
frate_half_widths = {}

Expand Down Expand Up @@ -412,11 +413,10 @@ def build_fringe_rate_profiles(uvd, uvb, keys=None, normed=True, combine_pols=Tr
if keys is None:
keys = uvd.get_antpairpols()

antpos_trf = uvd.antenna_positions # earth centered antenna positions
antnums = uvd.antenna_numbers # antenna numbers.
antpos_trf = uvd.telescope.antenna_positions # earth centered antenna positions
antnums = uvd.telescope.antenna_numbers # antenna numbers.

lat, lon, alt = uvd.telescope_location_lat_lon_alt_degrees
location = EarthLocation(lon=lon * units.deg, lat=lat * units.deg, height=alt * units.m)
location = uvd.telescope.location

# get topocentricl AzEl Beam coordinates.
hp = HEALPix(nside=uvb.nside, order=uvb.ordering)
Expand Down Expand Up @@ -462,7 +462,7 @@ def build_fringe_rate_profiles(uvd, uvb, keys=None, normed=True, combine_pols=Tr
# for if we are going to sum over polarizations.
# get redundancies (will only compute fr-profile once for each red group).
if reds is None:
reds = _get_key_reds(dict(zip(*uvd.get_ENU_antpos()[::-1])), keys)
reds = _get_key_reds(utils.get_ENU_antpos(uvd, asdict=True), keys)

for redgrp in reds:
# only explicitly calculate fr profile for the first vis in each redgroup.
Expand Down Expand Up @@ -610,7 +610,7 @@ def get_fringe_rate_limits(uvd, uvb=None, frate_profiles=None, percentile_low=5.
dfr = 1. / (nfr * np.mean(np.diff(np.unique(uvd.time_array))) * 1e-3 * SDAY_SEC)
fr_grid = np.arange(-nfr // 2, nfr // 2) * dfr

reds = _get_key_reds(dict(zip(*uvd.get_ENU_antpos()[::-1])), keys)
reds = _get_key_reds(utils.get_ENU_antpos(uvd, asdict=True), keys)

for redgrp in reds:
bl0 = redgrp[0]
Expand Down Expand Up @@ -1500,8 +1500,8 @@ def time_avg_data_and_write(input_data_list, output_data, t_avg, baseline_list=N
times=avg_times, lsts=avg_lsts, filetype=filetype, overwrite=clobber)
if flag_output is not None:
uv_avg = UVData()
uv_avg.read(output_data_name, use_future_array_shapes=True)
uvf = UVFlag(uv_avg, mode='flag', copy_flags=True, use_future_array_shapes=True)
uv_avg.read(output_data_name)
uvf = UVFlag(uv_avg, mode='flag', copy_flags=True)
uvf.to_waterfall(keep_pol=False, method='and')
uvf.write(os.path.join(os.path.dirname(flag_output),
os.path.basename(flag_output).replace('.h5', f'.interleave_{inum}.h5')), clobber=clobber)
Expand All @@ -1512,8 +1512,8 @@ def time_avg_data_and_write(input_data_list, output_data, t_avg, baseline_list=N
nsamples=fr.avg_nsamples, times=fr.avg_times, lsts=fr.avg_lsts)
if flag_output is not None:
uv_avg = UVData()
uv_avg.read(output_data, use_future_array_shapes=True)
uvf = UVFlag(uv_avg, mode='flag', copy_flags=True, use_future_array_shapes=True)
uv_avg.read(output_data)
uvf = UVFlag(uv_avg, mode='flag', copy_flags=True)
uvf.to_waterfall(keep_pol=False, method='and')
uvf.write(flag_output, clobber=clobber)

Expand Down Expand Up @@ -1817,7 +1817,7 @@ def load_tophat_frfilter_and_write(
# Read in the beam file if provided.
if beamfitsfile is not None:
uvb = UVBeam()
uvb.read_beamfits(beamfitsfile, use_future_array_shapes=True)
uvb.read_beamfits(beamfitsfile)
else:
uvb = None

Expand Down
Loading
Loading