Allow no model

docs/source/changelog/changelog.md

@@ -27,6 +27,7 @@ in pyuvdata.UVBeam and just called directly.
 * Analytic beams (subclassed from pyuvdata.AnalyticBeam) are now supported. Use
 the new `analytic-beam-yaml` option to configure them.
 * Added handling for `~` in paths in config yamls.
+* Altered the necessary code to allow for a run to be done without a model.
 
 ### Bug Fixes
 * Fixed a couple of bugs in beam and psf setup where the FFT direction or

docs/source/installation/installation.md

@@ -1,7 +1,6 @@
 # Installation
 
-The dependencies on pyfhd to run with FHD have been removed, this makes installing
-and `pyfhd` much easier than FHD. `pyfhd` is currently supported for Python 3.11+.
+The dependencies on pyfhd to run with FHD have been removed. This makes installing `pyfhd` much easier than FHD. `pyfhd` is currently supported for Python 3.11+.
 
 ## Installing to do some development on pyfhd?
 
@@ -17,8 +16,8 @@ can't install (HPC environments). `uv` is best used for development as it suppor
 a lot of features that Python developers of packages will appreciate (oh,
 and by the way, _it's fast_, like _really fast_ at installing packages).
 `conda` is best used where you need to use more than just Python, but perhaps a
-combination of languages and tools, `conda` has the ability to manage any
-compilers, CUDA versions, system tools etc.
+combination of languages and tools. `conda` has the ability to manage many
+compilers: CUDA versions, system tools etc.
 
 ### Make a venv environment
 This assumes you have python installed.
@@ -54,7 +53,7 @@ conda env create --file environment.yml python=3.11 # You can make it 3.11 or 3.
 conda activate pyfhd
 ```
 
-You can deactivate the environment anytime by using `conda deactivate`
+You can deactivate the environment at any time by using `conda deactivate`
 
 ### Installing pyfhd
 
@@ -70,7 +69,7 @@ You can Verify the installation using the version command
 pyfhd -v
 ```
 
-It should give you output that looks like this:
+It should give you an output that looks like this:
 
 ```
     ________________________________________________________________________

src/pyfhd/beam_setup/beam_utils.py

@@ -166,7 +166,7 @@ def beam_image(
     group_id = (psf["id"][pol_i, 0, :]).astype(int)
     # IDL psf structures often have "beam_gaussian_params" defined but set to zero
     # have to test that it's present *and* nonzero
-    if "beam_gaussian_params" in psf and psf["beam_gaussian_params"] != 0:
+    if "beam_gaussian_params" in psf and (psf["beam_gaussian_params"] != 0).all():
         beam_gaussian_params = psf["beam_gaussian_params"][:]
     else:
         beam_gaussian_params = None

src/pyfhd/healpix/healpix_utils.py

@@ -14,7 +14,6 @@
 from pyfhd.pyfhd_tools.pyfhd_utils import (
     angle_difference,
     histogram,
-    meshgrid,
     region_grow,
 )
 from pyfhd.pyfhd_tools.unit_conv import radec_to_altaz, radec_to_pixel
@@ -451,8 +450,8 @@ def phase_shift_uv_image(obs: dict) -> NDArray[np.complex128]:
     dx = (x - (obs["dimension"] / 2)) * (2 * np.pi / obs["dimension"])
     dy = (y - (obs["elements"] / 2)) * (2 * np.pi / obs["dimension"])
 
-    xvals = meshgrid(obs["dimension"], obs["elements"], 1) - (obs["dimension"] / 2)
-    yvals = meshgrid(obs["dimension"], obs["elements"], 2) - (obs["elements"] / 2)
+    xvals, yvals = np.meshgrid(np.arange(obs["dimension"]) - (obs["dimension"] / 2),
+                               np.arange(obs["elements"]) - (obs["elements"] / 2), indexing="ij")
 
     phase = xvals * dx + yvals * dy
     rephase_vals = np.cos(phase) + np.sin(phase) * 1j
@@ -567,6 +566,10 @@ def vis_model_freq_split(
         ]
         if np.size(fi_use) == 0:
             continue
+        if vis_model_arr is not None:
+            vis_model = vis_model_arr[polarization]
+        else:
+            vis_model = None
         gridding_dict = visibility_grid(
             vis_arr[polarization],
             vis_weights[polarization],
@@ -576,7 +579,7 @@ def vis_model_freq_split(
             polarization,
             pyfhd_config,
             logger,
-            model=vis_model_arr[polarization],
+            model=vis_model,
             fi_use=fi_use,
             bi_use=bi_use,
             verbose_logging=False,
@@ -597,7 +600,10 @@ def vis_model_freq_split(
             variance_uv_arr[fi] = (
                 gridding_dict["variance"] * rephase_use * gridding_dict["n_vis"]
             )
-            model_uv_arr[fi] = gridding_dict["model_return"] * gridding_dict["n_vis"]
+            if vis_model is not None:
+                model_uv_arr[fi] = (
+                    gridding_dict["model_return"] * gridding_dict["n_vis"]
+                )
 
         if fft:
             # No x_range and y_range hence no check for it here
@@ -619,13 +625,14 @@ def vis_model_freq_split(
                 degpix=obs["degpix"],
             )
             variance_arr[fi] *= gridding_dict["n_vis"]
-            model_arr[fi], _, _ = dirty_image_generate(
-                gridding_dict["model_return"],
-                pyfhd_config,
-                logger,
-                degpix=obs["degpix"],
-            )
-            model_arr[fi] *= gridding_dict["n_vis"]
+            if vis_model is not None:
+                model_arr[fi], _, _ = dirty_image_generate(
+                    gridding_dict["model_return"],
+                    pyfhd_config,
+                    logger,
+                    degpix=obs["degpix"],
+                )
+                model_arr[fi] *= gridding_dict["n_vis"]
         else:
             dirty_arr[fi] = gridding_dict["image_uv"] * gridding_dict["n_vis"]
             weights_arr[fi] = (
@@ -634,7 +641,8 @@ def vis_model_freq_split(
             variance_arr[fi] = (
                 gridding_dict["variance"] * rephase_use * gridding_dict["n_vis"]
             )
-            model_arr[fi] = gridding_dict["model_return"] * gridding_dict["n_vis"]
+            if vis_model is not None:
+                model_arr[fi] = gridding_dict["model_return"] * gridding_dict["n_vis"]
     obs["n_vis"] = n_vis_use
 
     if save_uvf:
@@ -673,15 +681,16 @@ def vis_model_freq_split(
             f"{pyfhd_config['obs_id']}_{uvf_name}_{obs['pol_names'][polarization]}_variance_uv_arr_gridded_uvf.h5",
             logger=logger,
         )
-        save(
-            Path(
-                uvf_dir,
+        if vis_model is not None:
+            save(
+                Path(
+                    uvf_dir,
+                    f"{pyfhd_config['obs_id']}_{uvf_name}_{obs['pol_names'][polarization]}_model_uv_arr_gridded_uvf.h5",
+                ),
+                h5_save_dict,
                 f"{pyfhd_config['obs_id']}_{uvf_name}_{obs['pol_names'][polarization]}_model_uv_arr_gridded_uvf.h5",
-            ),
-            h5_save_dict,
-            f"{pyfhd_config['obs_id']}_{uvf_name}_{obs['pol_names'][polarization]}_model_uv_arr_gridded_uvf.h5",
-            logger=logger,
-        )
+                logger=logger,
+            )
 
     cube_split = {
         "obs": obs,

src/pyfhd/io/pyfhd_quickview.py

@@ -79,14 +79,14 @@ def quickview(
     obs: dict,
     psf: dict,
     params: dict,
-    cal: dict,
+    cal: dict | None,
     vis_arr: NDArray[np.complex128],
     vis_weights: NDArray[np.float64],
     image_uv: NDArray[np.complex128],
     weights_uv: NDArray[np.complex128],
     variance_uv: NDArray[np.float64],
     uniform_filter_uv: NDArray[np.float64],
-    model_uv: NDArray[np.complex128],
+    model_uv: NDArray[np.complex128] | None,
     pyfhd_config: dict,
     logger: Logger,
 ) -> None:
@@ -273,26 +273,28 @@ def quickview(
             beam_ptr=beam_base_out[pol_i],
         )
         filter_arr[pol_i] = filter
-        instr_model_arr[pol_i], filter, _ = dirty_image_generate(
-            model_uv[pol_i],
-            pyfhd_config,
-            logger,
-            uniform_filter_uv=uniform_filter_uv,
-            degpix=obs_out["degpix"],
-            weights=weights_uv[pol_i],
-            pad_uv_image=pyfhd_config["pad_uv_image"],
-            filter=filter,
-            not_real=complex_flag,
-            beam_ptr=beam_base_out[pol_i],
-        )
+        if model_uv is not None:
+            instr_model_arr[pol_i], filter, _ = dirty_image_generate(
+                model_uv[pol_i],
+                pyfhd_config,
+                logger,
+                uniform_filter_uv=uniform_filter_uv,
+                degpix=obs_out["degpix"],
+                weights=weights_uv[pol_i],
+                pad_uv_image=pyfhd_config["pad_uv_image"],
+                filter=filter,
+                not_real=complex_flag,
+                beam_ptr=beam_base_out[pol_i],
+            )
     renorm_factor = get_image_renormalization(
         obs_out, weights_uv, beam_base_out, filter_arr, pyfhd_config, logger
     )
     # Reshape renorm factor to multiply per polarization without loop to [obs["n_pol"], 1, 1]
     renorm_factor = np.expand_dims(renorm_factor.reshape([obs_out["n_pol"], 1]), -1)
     instr_dirty_arr *= renorm_factor
-    instr_model_arr *= renorm_factor
-    instr_residual_arr = instr_dirty_arr - instr_model_arr
+    if model_uv is not None:
+        instr_model_arr *= renorm_factor
+        instr_residual_arr = instr_dirty_arr - instr_model_arr
     # Get the pol_names
     pol_names = obs["pol_names"]
     # The cross-polarization XY and YX images are both complex, but are conjugate mirrors of each other
@@ -387,9 +389,10 @@ def quickview(
         png_output.mkdir(exist_ok=True)
     for pol_i in range(obs["n_pol"]):
         logger.info(f"Saving the FITS files for polarization {pol_names[pol_i]}")
-        instr_residual = instr_residual_arr[pol_i] * beam_correction_out[pol_i]
         instr_dirty = instr_dirty_arr[pol_i] * beam_correction_out[pol_i]
-        instr_model = instr_model_arr[pol_i] * beam_correction_out[pol_i]
+        if model_uv is not None:
+            instr_residual = instr_residual_arr[pol_i] * beam_correction_out[pol_i]
+            instr_model = instr_model_arr[pol_i] * beam_correction_out[pol_i]
         beam_use = beam_base_out[pol_i]
 
         # Write the fits files for the dirty images
@@ -400,20 +403,21 @@ def quickview(
         fits_file_apparent.writeto(
             Path(fits_output, f"{instr_dirty_name}.fits"), overwrite=True
         )
-        fits_file_apparent.data = instr_model
-        instr_model_name = (
-            f"{pyfhd_config['obs_id']}_{filter_name}_model_{pol_names[pol_i]}"
-        )
-        fits_file_apparent.writeto(
-            Path(fits_output, f"{instr_model_name}.fits"), overwrite=True
-        )
-        fits_file_apparent.data = instr_residual
-        instr_residual_name = (
-            f"{pyfhd_config['obs_id']}_{filter_name}_residual_{pol_names[pol_i]}"
-        )
-        fits_file_apparent.writeto(
-            Path(fits_output, f"{instr_residual_name}.fits"), overwrite=True
-        )
+        if model_uv is not None:
+            fits_file_apparent.data = instr_model
+            instr_model_name = (
+                f"{pyfhd_config['obs_id']}_{filter_name}_model_{pol_names[pol_i]}"
+            )
+            fits_file_apparent.writeto(
+                Path(fits_output, f"{instr_model_name}.fits"), overwrite=True
+            )
+            fits_file_apparent.data = instr_residual
+            instr_residual_name = (
+                f"{pyfhd_config['obs_id']}_{filter_name}_residual_{pol_names[pol_i]}"
+            )
+            fits_file_apparent.writeto(
+                Path(fits_output, f"{instr_residual_name}.fits"), overwrite=True
+            )
         fits_file.data = beam_use
         beam_name = f"{pyfhd_config['obs_id']}_beam_{pol_names[pol_i]}"
         fits_file.writeto(Path(fits_output, f"{beam_name}.fits"), overwrite=True)
@@ -431,18 +435,19 @@ def quickview(
                 title=f"Dirty Image {pol_names[pol_i]}",
                 logger=logger,
             )
-            plot_fits_image(
-                Path(fits_output, f"{instr_model_name}.fits"),
-                Path(png_output, f"{instr_model_name}.png"),
-                title=f"Model Image {pol_names[pol_i]}",
-                logger=logger,
-            )
-            plot_fits_image(
-                Path(fits_output, f"{instr_residual_name}.fits"),
-                Path(png_output, f"{instr_residual_name}.png"),
-                title=f"Residual Image {pol_names[pol_i]}",
-                logger=logger,
-            )
+            if model_uv is not None:
+                plot_fits_image(
+                    Path(fits_output, f"{instr_model_name}.fits"),
+                    Path(png_output, f"{instr_model_name}.png"),
+                    title=f"Model Image {pol_names[pol_i]}",
+                    logger=logger,
+                )
+                plot_fits_image(
+                    Path(fits_output, f"{instr_residual_name}.fits"),
+                    Path(png_output, f"{instr_residual_name}.png"),
+                    title=f"Residual Image {pol_names[pol_i]}",
+                    logger=logger,
+                )
             plot_fits_image(
                 Path(fits_output, f"{beam_name}.fits"),
                 Path(png_output, f"{beam_name}.png"),

src/pyfhd/pyfhd.py

@@ -275,6 +275,9 @@ def main():
                     f"Checkpoint Loaded: Uncalibrated visibility parameters, array and weights and the observation metadata dictionary loaded from {Path(pyfhd_config['output_dir'], 'obs_checkpoint.h5')}"
                 )
 
+        # Give cal and vis_model_arr a value to avoid no value error (Default = None, but can be overwritten below)
+        cal = None
+        vis_model_arr = None
         # If the calibration checkpoint exists, load it now before loading in the beam
         # to get the observation metadata and visibility parameters
         if pyfhd_config["calibrate_checkpoint"]:
@@ -338,16 +341,20 @@ def main():
                 logger,
             )
 
-            # Get the vis_model_arr from a UVFITS file or SAV files and flag any issues
-            vis_model_arr_start = time.time()
-            vis_model_arr = vis_model_transfer(pyfhd_config, obs, params, logger)
-            vis_model_arr_end = time.time()
-            _print_time_diff(
-                vis_model_arr_start,
-                vis_model_arr_end,
-                "Model Imported and Flagged From UVFITS",
-                logger,
-            )
+            if (
+                pyfhd_config["calibrate_visibilities"]
+                or pyfhd_config["model_file_path"] is not None
+            ):
+                # Get the vis_model_arr from a UVFITS file or SAV files and flag any issues
+                vis_model_arr_start = time.time()
+                vis_model_arr = vis_model_transfer(pyfhd_config, obs, params, logger)
+                vis_model_arr_end = time.time()
+                _print_time_diff(
+                    vis_model_arr_start,
+                    vis_model_arr_end,
+                    "Model Imported and Flagged From UVFITS",
+                    logger,
+                )
 
             # Skipped initializing the cal structure as it mostly just copies values from the obs, params, config and the skymodel from FHD
             # However, there is resulting cal structure for logging and output purposes to store the resulting gain and any other associated
@@ -511,6 +518,8 @@ def main():
                     (obs["n_pol"], obs["elements"], obs["dimension"]),
                     dtype=np.complex128,
                 )
+            else:
+                model_uv = None
             # Since it's done per polarization, we can do multi-processing if it's not fast enough
             for pol_i in range(obs["n_pol"]):
                 logger.info(
@@ -524,6 +533,10 @@ def main():
                     no_conjugate = True
                 else:
                     no_conjugate = False
+                if vis_model_arr is not None:
+                    vis_model = vis_model_arr[pol_i]
+                else:
+                    vis_model = None
                 gridding_dict = visibility_grid(
                     vis_arr[pol_i],
                     vis_weights[pol_i],
@@ -535,7 +548,7 @@ def main():
                     logger,
                     uniform_flag=uniform_flag,
                     no_conjugate=no_conjugate,
-                    model=vis_model_arr[pol_i],
+                    model=vis_model,
                 )
                 if len(gridding_dict.keys()) != 0:
                     image_uv[pol_i] = gridding_dict["image_uv"]
@@ -605,6 +618,8 @@ def main():
             uniform_filter_uv = grid_checkpoint["uniform_filter_uv"]
             if "model_uv" in grid_checkpoint:
                 model_uv = grid_checkpoint["model_uv"]
+            else:
+                model_uv = None
             del grid_checkpoint
             logger.info(
                 f"Checkpoint Loaded: The Gridded UV Planes loaded from {Path(pyfhd_config['output_dir'], 'gridding_checkpoint.h5')}"