diff --git a/PyFHD/beam_setup/beam.py b/PyFHD/beam_setup/beam.py index 26c4dc1..a7ef0fb 100644 --- a/PyFHD/beam_setup/beam.py +++ b/PyFHD/beam_setup/beam.py @@ -320,6 +320,6 @@ def create_psf(obs: dict, pyfhd_config: dict, logger: Logger) -> dict | File: sys.exit(1) raise ValueError( f"Unknown beam file type {pyfhd_config['beam_file_path'].suffix}. " - "Please use a .sav, .h5, .hdf5" + "Please use a .sav, .h5, .hdf5 " "If you meant for PyFHD to do the beam forming, please set the beam_file_path to None (~ in YAML)." ) diff --git a/PyFHD/beam_setup/beam_utils.py b/PyFHD/beam_setup/beam_utils.py index f2c46a3..3b58bdb 100644 --- a/PyFHD/beam_setup/beam_utils.py +++ b/PyFHD/beam_setup/beam_utils.py @@ -156,7 +156,7 @@ def beam_image( """ psf_dim = psf["dim"] - freq_norm = psf["fnorm"] + freq_norm = psf["freq_norm"] pix_horizon = psf["pix_horizon"] group_id = psf["id"][pol_i, 0, :] if "beam_gaussian_params" in psf: @@ -170,10 +170,12 @@ def beam_image( freq_norm = freq_norm[:] pix_horizon = pix_horizon[0] dimension = elements = obs["dimension"] - xl = dimension / 2 - psf_dim / 2 + 1 - xh = dimension / 2 - psf_dim / 2 + psf_dim - yl = elements / 2 - psf_dim / 2 + 1 - yh = elements / 2 - psf_dim / 2 + psf_dim + # these should all be integers b/c dimensions are usually even numbers. + # but they have to be cast to ints to be used in slicing. + xl = int(dimension / 2 - psf_dim / 2 + 1) + xh = int(dimension / 2 - psf_dim / 2 + psf_dim) + yl = int(elements / 2 - psf_dim / 2 + 1) + yh = int(elements / 2 - psf_dim / 2 + psf_dim) group_n, _, ri_id = histogram(group_id, min=0) gi_use = np.nonzero(group_n) diff --git a/PyFHD/calibration/calibration_utils.py b/PyFHD/calibration/calibration_utils.py index cdc5c42..82a90fc 100644 --- a/PyFHD/calibration/calibration_utils.py +++ b/PyFHD/calibration/calibration_utils.py @@ -49,7 +49,7 @@ def vis_extract_autocorr( if autocorr_i.size > 0: auto_tile_i = obs["baseline_info"]["tile_a"][autocorr_i] - 1 # As auto_tile_i is used for indexing we need to make it an integer array - auto_tile_i = auto_tile_i.astype(np.integer) + auto_tile_i = auto_tile_i.astype(int) auto_tile_i_single = np.unique(auto_tile_i) # expect it as a list of 2D arrays, so there might be trouble if not pyfhd_config["cal_time_average"]: @@ -1557,6 +1557,10 @@ def vis_baseline_hist( wh_noflag = np.where(np.abs(model_vals) > 0)[0] if wh_noflag.size > 0: inds = inds[wh_noflag] + # inds changed so we need to update model_vals + # otherwise it has a different shape than vis_cal_use below + # causing errors + model_vals = (vis_model_arr[pol_i]).flatten()[inds] else: continue # if Keyword_Set(calibration_visibilities_subtract) THEN BEGIN diff --git a/PyFHD/data_setup/uvfits.py b/PyFHD/data_setup/uvfits.py index 29aef60..0a8b577 100644 --- a/PyFHD/data_setup/uvfits.py +++ b/PyFHD/data_setup/uvfits.py @@ -88,17 +88,31 @@ def extract_header( pyfhd_header["frequency_array"] = ( np.arange(pyfhd_header["n_freq"]) - freq_ref_i ) * pyfhd_header["freq_res"] + pyfhd_header["freq_ref"] + # the following is guaranteed from the uvfits memo, logic stolen from pyuvdata + if params_header["naxis"] == 7: + ra_axis = 6 + dec_axis = 7 + else: + ra_axis = 5 + dec_axis = 6 + # obsra/obsdec/ra/dec are not standard uvfits keywords try: pyfhd_header["obsra"] = params_header["obsra"] except KeyError: logger.warning("OBSRA not found in UVFITS file") - pyfhd_header["obsra"] = params_header["ra"] + if "ra" in params_header: + pyfhd_header["obsra"] = params_header["ra"] + else: + pyfhd_header["obsra"] = params_header[f"CRVAL{ra_axis}"] try: pyfhd_header["obsdec"] = params_header["obsdec"] except KeyError: logger.warning("OBSDEC not found in UVFITS file") - pyfhd_header["obsdec"] = params_header["dec"] + if "dec" in params_header: + pyfhd_header["obsdec"] = params_header["dec"] + else: + pyfhd_header["obsdec"] = params_header[f"CRVAL{dec_axis}"] # Put in locations of instrument from FITS file or from Astropy site data # If you want to see the list of current site names using EarthLocation.get_site_names() # If you want to use PyFHD with HERA in the future @@ -113,6 +127,8 @@ def extract_header( # Can also do MWA or Murchison Widefield Array location = EarthLocation("mwa") + # These are all non-standard uvfits keywords. This information is stored in + # the antenna table. See pyuvdata for the right way to do this. try: pyfhd_header["lon"] = params_header["lon"] except KeyError: diff --git a/PyFHD/gridding/gridding_utils.py b/PyFHD/gridding/gridding_utils.py index 21bac94..0771ce7 100644 --- a/PyFHD/gridding/gridding_utils.py +++ b/PyFHD/gridding/gridding_utils.py @@ -300,7 +300,7 @@ def dirty_image_generate( logger: Logger, uniform_filter_uv: NDArray[np.float64] | None = None, mask: NDArray[np.integer] | None = None, - baseline_threshold: int | float = 0, + baseline_threshold: int | float | None = None, normalization: float | NDArray[np.float64] | None = None, resize: int | None = None, width_smooth: int | float | None = None, @@ -330,7 +330,7 @@ def dirty_image_generate( mask : NDArray[np.integer] | None, optional A 2D {u,v} mask to apply before image creation, by default None baseline_threshold : int | float, optional - The maximum baseline length to include in units of pixels, by default 0 + The maximum baseline length to include in units of pixels, default None normalization : float | NDArray[np.float64] | None, optional A value by which to normalize the image by, by default None resize : int | None, optional diff --git a/PyFHD/plotting/image.py b/PyFHD/plotting/image.py index 6236dcf..2de8402 100644 --- a/PyFHD/plotting/image.py +++ b/PyFHD/plotting/image.py @@ -176,7 +176,7 @@ def quick_image( count_missing = len(wh_missing[0]) if count_missing > 0: image[wh_missing] = np.nan - missing_color = 0 + missing_color = 0 else: count_missing = 0 wh_missing = None diff --git a/PyFHD/pyfhd.py b/PyFHD/pyfhd.py index a6ea170..1bace64 100644 --- a/PyFHD/pyfhd.py +++ b/PyFHD/pyfhd.py @@ -144,6 +144,45 @@ def main(): ) pyfhd_config["checkpoint_dir"].mkdir(exist_ok=True) + obs_checkpoint_file = Path( + pyfhd_config["checkpoint_dir"], + f"{checkpoint_name}_obs_checkpoint.h5", + ) + if ( + pyfhd_config["obs_checkpoint"] is not None + and not obs_checkpoint_file.exists() + ): + logger.warning( + "obs_checkpoint is set but obs checkpoint file does not exist. Recalculating obs." + ) + pyfhd_config["obs_checkpoint"] = None + + cal_checkpoint_file = Path( + pyfhd_config["checkpoint_dir"], + f"{checkpoint_name}_calibrate_checkpoint.h5", + ) + if ( + pyfhd_config["calibrate_checkpoint"] is not None + and not cal_checkpoint_file.exists() + ): + logger.warning( + "calibrate_checkpoint is set but cal checkpoint file does not exist. Recalculating cal." + ) + pyfhd_config["calibrate_checkpoint"] = None + + grid_checkpoint_file = Path( + pyfhd_config["checkpoint_dir"], + f"{checkpoint_name}_gridding_checkpoint.h5", + ) + if ( + pyfhd_config["gridding_checkpoint"] is not None + and not grid_checkpoint_file.exists() + ): + logger.warning( + "gridding_checkpoint is set but grid checkpoint file does not exist. Recalculating grid." + ) + pyfhd_config["gridding_checkpoint"] = None + if ( pyfhd_config["obs_checkpoint"] is None and pyfhd_config["calibrate_checkpoint"] is None @@ -217,10 +256,7 @@ def main(): "vis_weights": vis_weights, } save( - Path( - pyfhd_config["checkpoint_dir"], - f"{checkpoint_name}_obs_checkpoint.h5", - ), + obs_checkpoint_file, checkpoint, "obs_checkpoint", logger=logger, @@ -231,11 +267,8 @@ def main(): ) else: # Load the checkpoint and initialize the required variables - if ( - pyfhd_config["obs_checkpoint"] - and Path(pyfhd_config["obs_checkpoint"]).exists() - ): - obs_checkpoint = load(pyfhd_config["obs_checkpoint"], logger=logger) + if pyfhd_config["obs_checkpoint"]: + obs_checkpoint = load(obs_checkpoint_file, logger=logger) obs = obs_checkpoint["obs"] params = obs_checkpoint["params"] vis_arr = obs_checkpoint["vis_arr"] @@ -249,9 +282,9 @@ def main(): # to get the observation metadata and visibility parameters if ( pyfhd_config["calibrate_checkpoint"] is not None - and Path(pyfhd_config["calibrate_checkpoint"]).exists() + and cal_checkpoint_file.exists() ): - cal_checkpoint = load(pyfhd_config["calibrate_checkpoint"], logger=logger) + cal_checkpoint = load(cal_checkpoint_file, logger=logger) obs = cal_checkpoint["obs"] params = cal_checkpoint["params"] vis_arr = cal_checkpoint["vis_arr"] @@ -272,7 +305,10 @@ def main(): ) # Check if the calibrate checkpoint has been used, if not run the calibration steps - if pyfhd_config["calibrate_checkpoint"] is None: + if ( + pyfhd_config["calibrate_checkpoint"] is None + or not cal_checkpoint_file.exists() + ): if pyfhd_config["deproject_w_term"] is not None: w_term_start = time.time() vis_arr = simple_deproject_w_term( @@ -394,10 +430,7 @@ def main(): "cal": cal, } save( - Path( - pyfhd_config["checkpoint_dir"], - f"{checkpoint_name}_calibrate_checkpoint.h5", - ), + cal_checkpoint_file, checkpoint, "calibrate_checkpoint", logger=logger, @@ -461,7 +494,7 @@ def main(): pyfhd_successful = True sys.exit(0) - if ( + if "image_info" not in psf or ( psf["image_info"]["image_power_beam_arr"] is not None and psf["image_info"]["image_power_beam_arr"].shape == 1 ): @@ -471,7 +504,8 @@ def main(): if ( pyfhd_config["recalculate_grid"] - and pyfhd_config["gridding_checkpoint"] is None + or pyfhd_config["gridding_checkpoint"] is None + or not grid_checkpoint_file.exists() ): grid_start = time.time() image_uv = np.empty( @@ -535,6 +569,8 @@ def main(): if vis_model_arr is not None: model_uv = crosspol_reformat(model_uv) if pyfhd_config["gridding_plots"]: + # TODO: move this after the checkpointing so an error in plotting + # doesn't require rerunning gridding. logger.info( f"Plotting the continuum gridding outputs into {pyfhd_config['output_dir']/'plots'/'gridding'}" ) @@ -557,10 +593,7 @@ def main(): if vis_model_arr is not None: checkpoint["model_uv"] = model_uv save( - Path( - pyfhd_config["checkpoint_dir"], - f"{checkpoint_name}_gridding_checkpoint.h5", - ), + grid_checkpoint_file, checkpoint, "gridding_checkpoint", logger=logger, @@ -573,9 +606,7 @@ def main(): _print_time_diff(grid_start, grid_end, "Visibilities gridded", logger) else: if pyfhd_config["gridding_checkpoint"]: - grid_checkpoint = load( - pyfhd_config["gridding_checkpoint"], logger=logger - ) + grid_checkpoint = load(grid_checkpoint_file, logger=logger) image_uv = grid_checkpoint["image_uv"] weights_uv = grid_checkpoint["weights_uv"] variance_uv = grid_checkpoint["variance_uv"] diff --git a/PyFHD/pyfhd_tools/pyfhd_setup.py b/PyFHD/pyfhd_tools/pyfhd_setup.py index 21e1f7a..c65fd8f 100644 --- a/PyFHD/pyfhd_tools/pyfhd_setup.py +++ b/PyFHD/pyfhd_tools/pyfhd_setup.py @@ -961,7 +961,7 @@ def _check_file_exists(config: dict, key: str) -> int: """ if config[key]: # If it doesn't exist, add error message - if not Path(config[key]).exists(): + if not Path(config[key]).expanduser().resolve().exists(): logging.error( "{} has been enabled with a path that doesn't exist, check the path.".format( key @@ -1142,6 +1142,9 @@ def pyfhd_logger(pyfhd_config: dict) -> Tuple[logging.Logger, Path]: if pyfhd_config["get_sample_data"]: output_dir = Path(pyfhd_config["output_path"]) else: + pyfhd_config["output_path"] = ( + Path(pyfhd_config["output_path"]).expanduser().resolve() + ) output_dir = Path(pyfhd_config["output_path"], dir_name) if Path.is_dir(output_dir): output_dir_exists = True @@ -1217,6 +1220,7 @@ def pyfhd_setup(options: argparse.Namespace) -> Tuple[dict, logging.Logger]: pyfhd_config["output_dir"] = output_dir pyfhd_config["top_level_dir"] = str(output_dir).split("/")[-1] # Check input_path exists and obs_id uvfits and metafits files exist (Error) + pyfhd_config["input_path"] = Path(pyfhd_config["input_path"]).expanduser().resolve() if not pyfhd_config["input_path"].exists(): logger.error( "{} doesn't exist, please check your input path".format(options.input_path) @@ -1272,14 +1276,15 @@ def pyfhd_setup(options: argparse.Namespace) -> Tuple[dict, logging.Logger]: pyfhd_config["interpolate_kernel"] = False # If the user has set a beam file, check it exists (Error) - if ( - pyfhd_config["beam_file_path"] is not None - and not Path(pyfhd_config["beam_file_path"]).exists() - ): - logger.error( - f"Beam file {pyfhd_config['beam_file_path']} does not exist, please check your input path" + if pyfhd_config["beam_file_path"] is not None: + pyfhd_config["beam_file_path"] = ( + Path(pyfhd_config["beam_file_path"]).expanduser().resolve() ) - errors += 1 + if not Path(pyfhd_config["beam_file_path"]).exists(): + logger.error( + f"Beam file {pyfhd_config['beam_file_path']} does not exist, please check your input path" + ) + errors += 1 if pyfhd_config["beam_file_path"] is None: logger.info("No beam file was set, PyFHD will calculate the beam.") @@ -1384,6 +1389,9 @@ def pyfhd_setup(options: argparse.Namespace) -> Tuple[dict, logging.Logger]: # if importing model visiblities from a uvfits file, check that file # exists if pyfhd_config["model_file_path"]: + pyfhd_config["model_file_path"] = ( + Path(pyfhd_config["model_file_path"]).expanduser().resolve() + ) errors += _check_file_exists(pyfhd_config, "model_file_path") if pyfhd_config["model_file_path"] == "sav": diff --git a/PyFHD/resources/1088285600_example/1088285600_example.yaml b/PyFHD/resources/1088285600_example/1088285600_example.yaml index 3566c78..9e61ef2 100644 --- a/PyFHD/resources/1088285600_example/1088285600_example.yaml +++ b/PyFHD/resources/1088285600_example/1088285600_example.yaml @@ -18,7 +18,7 @@ deproject-w-term : ~ # Checkpointing save-checkpoints: true obs-checkpoint: ~ -calibration-checkpoint: ~ +calibrate-checkpoint: ~ gridding-checkpoint: ~ # Instrument diff --git a/PyFHD/resources/config/pyfhd.yaml b/PyFHD/resources/config/pyfhd.yaml index da4e376..4de6191 100644 --- a/PyFHD/resources/config/pyfhd.yaml +++ b/PyFHD/resources/config/pyfhd.yaml @@ -18,7 +18,7 @@ deproject-w-term : ~ # Checkpointing save-checkpoints: true obs-checkpoint: ~ -calibration-checkpoint: ~ +calibrate-checkpoint: ~ gridding-checkpoint: ~ # Instrument diff --git a/PyFHD/source_modeling/vis_model_transfer.py b/PyFHD/source_modeling/vis_model_transfer.py index e338e2e..397cabc 100644 --- a/PyFHD/source_modeling/vis_model_transfer.py +++ b/PyFHD/source_modeling/vis_model_transfer.py @@ -105,18 +105,36 @@ def import_vis_model_from_sav( params_model : dict The parameters for said model used for flagging """ + fhd_subdirs = { + "params": "metadata", + "vis": "vis_data", + } + try: path = Path( pyfhd_config["model_file_path"], f"{pyfhd_config['obs_id']}_params.sav" ) + if not path.exists(): + path = Path( + pyfhd_config["model_file_path"], + fhd_subdirs["params"], + f"{pyfhd_config['obs_id']}_params.sav", + ) params_model = readsav(path) params_model = recarray_to_dict(params_model.params) # Read in the first polarization from pol_names pol_i = 0 + vis_path_parts = [pyfhd_config["model_file_path"]] path = Path( - pyfhd_config["model_file_path"], + *vis_path_parts, f"{pyfhd_config['obs_id']}_vis_model_{obs['pol_names'][pol_i]}.sav", ) + if not path.exists(): + vis_path_parts = [pyfhd_config["model_file_path"], fhd_subdirs["vis"]] + path = Path( + *vis_path_parts, + f"{pyfhd_config['obs_id']}_vis_model_{obs['pol_names'][pol_i]}.sav", + ) curr_vis_model = readsav(path) if isinstance(curr_vis_model, dict): curr_vis_model = curr_vis_model["vis_model_ptr"] @@ -132,7 +150,7 @@ def import_vis_model_from_sav( vis_model_arr[pol_i] = curr_vis_model for pol_i in range(1, obs["n_pol"]): path = Path( - pyfhd_config["model_file_path"], + *vis_path_parts, f"{pyfhd_config['obs_id']}_vis_model_{obs['pol_names'][pol_i]}.sav", ) curr_vis_model = readsav(path) diff --git a/docs/source/develop/contribution_guide.md b/docs/source/develop/contribution_guide.md index 7c8fcbb..b24e963 100644 --- a/docs/source/develop/contribution_guide.md +++ b/docs/source/develop/contribution_guide.md @@ -124,7 +124,7 @@ The very first question you should ask yourself before making a new function is, Take the `histogram` function for example, are you wondering why we did a histogram function? -After all, `IDL's` histogram function can be done using a combination of `np.histogram` and the reverse indices can technically be done using many indexing functions available in NumPy as well. Well it turns out the NumPy functions like `np.histogram` suffer from a fatal flaw, they're awfully slow with large arrays, and this isn't necessarily their fault either, NumPy (rightly for their particular case) prioritised compatibility over raw speed. Other than speed, a function needed to be made to create the reverse indices part of `IDL's` histogram anyway as there is no `SciPy` or `NumPy` equivalent. As such it was decided to rewrite the histogram from scratch with speed as the priority given that `histogram` in `PyFHD` get's called hundreds, maybe hundreds of times. +After all, `IDL's` histogram function can be done using a combination of `np.histogram` and the reverse indices can technically be done using many indexing functions available in NumPy as well. Well it turns out the NumPy functions like `np.histogram` suffer from a fatal flaw, they're awfully slow with large arrays, and this isn't necessarily their fault either, NumPy (rightly for their particular case) prioritised compatibility over raw speed. Other than speed, a function needed to be made to create the reverse indices part of `IDL's` histogram anyway as there is no `SciPy` or `NumPy` equivalent. As such it was decided to rewrite the histogram from scratch with speed as the priority given that `histogram` in `PyFHD` get's called hundreds, maybe thousands of times. To summarise, only make a new function if it hasn't been made before, or the ones that do exist do not meet your requirements. The requirements initally should not include anything in regards to the speed of the function unless you know before hand that the function is going to be called a lot or in a loop. @@ -570,6 +570,9 @@ In order to run these tests you need to do the following: Put the zip file in the following directory inside the repository (the directory given is relative to the root of the repository): `PyFHD/resources/test_data` 2. Unzip the zip file in the previously mentioned `PyFHD/resources/test_data` directory. +After unzipping, there will be a new folder named `pyfhd_test_data`. Move all +the contents out of that folder up into the `PyFHD/resources/test_data` folder +and delete the now empty `pyfhd_test_data` folder. 3. In a terminal, at the root of the repository, run the following command: diff --git a/docs/source/tutorial/tutorial.rst b/docs/source/tutorial/tutorial.rst index b6e1d3b..4779a16 100644 --- a/docs/source/tutorial/tutorial.rst +++ b/docs/source/tutorial/tutorial.rst @@ -507,6 +507,9 @@ likely you'll need to adjust the default configuration file to suit your needs. `pyfhd.yaml `_ + Note that to provide a ``None`` input in the yaml you must use ``~`` or ``null`` + as those are translated to ``None`` when yamls are read into python. + Some files can be discovered automatically through the ``input-path`` option of ``PyFHD`` so read through the usage help text to work out how you wish to configure your input. ``PyFHD`` is rather flexible on how you do your input as many of the files you may require can be in completely separate directories. @@ -1388,7 +1391,7 @@ are not complete (for example the model doesn't also create the params file), bu dim = H5D_READ(dim) fbin_i = H5D_OPEN(file_id, "fbin_i") fbin_i = H5D_READ(fbin_i) - fnorm = H5D_OPEN(file_id, "fnorm") + fnorm = H5D_OPEN(file_id, "freq_norm") fnorm = H5D_READ(fnorm) freq = H5D_OPEN(file_id, "freq") freq = H5D_READ(freq) diff --git a/environment.yml b/environment.yml index b79ef8b..cbb7f45 100644 --- a/environment.yml +++ b/environment.yml @@ -1,7 +1,6 @@ name: pyfhd channels: - conda-forge - - defaults dependencies: - python - astropy