Using Pyralysis
===============

What this guide covers
----------------------
This guide provides a **practical workflow** for using Pyralysis for radio astronomy
**imaging** (also called **image synthesis** / **image formation**). You'll learn how
to load data, simulate visibilities, grid, optimize (reconstruct) an image, and save
results. Pyralysis often says *reconstruction* or *optimization* for the imaging step;
see :doc:`glossary`.

.. tip::
   Best practice: Start with small datasets to get familiar with the workflow, then scale up to larger data using Dask.
   For **pipeline** workflows (``SimulationPipeline``, ``ImagerPipeline``), optional local or cluster Dask setup is described in :doc:`pipelines_distributed`.

.. note::
   Common mistake: Forgetting to set the correct cell size or image size for your data. Double-check these parameters!

---

Basic Workflow
--------------

A typical workflow in Pyralysis involves:

1. **Loading Measurement Sets**
2. **Simulating visibilities from an image (optional)**
3. **Applying gridding techniques**
4. **Imaging / reconstruction** — run an optimization algorithm (or form a dirty image)
5. **Saving the final image**

.. note::
   Want the fastest way to generate a dataset? See :ref:`creating-dataset` for the minimal example.

---

Loading Measurement Sets
------------------------

Pyralysis processes radio interferometric data stored in **Measurement Set (MS) format**.
Under the hood, it leverages ``dask-ms`` for efficient handling of large datasets.

**Example: Loading a Measurement Set**

.. code-block:: python

   from pyralysis.io import DaskMS

   # Define the path to your Measurement Set file
   ms_file = "/path/to/data.ms"

   # Load the dataset using Dask-MS (lazy loading enabled)
   dataset = DaskMS(input_name=ms_file).read()

   print("Loaded dataset structure:", dataset)

---

Simulating a Synthetic Dataset
-----------------------------

You can use Pyralysis to generate synthetic radio interferometric datasets for testing and benchmarking. This is useful for developing and validating imaging algorithms.

**Example: Simulate an in-memory dataset (and optionally write a Measurement Set)**

.. note::
   The simulator returns a :class:`~pyralysis.base.dataset.Dataset` in memory. To
   choose the output ``.ms`` path, pass ``output_name`` to
   :class:`~pyralysis.io.daskms.DaskMS` and call ``write`` (not an argument on
   ``Simulator``). See :ref:`creating-dataset` and :doc:`io_operations`.

.. code-block:: python

   from pyralysis.io.antenna_config_io import AntennaConfigurationIo
   from pyralysis.simulation import Simulator
   from pyralysis.models.sky import PointSource
   from pyralysis.io import DaskMS

   # Load array configuration
   interferometer = AntennaConfigurationIo(input_name="/path/to/array.cfg").read()
   interferometer.configure_observation(
       min_frequency_hz=1e9, max_frequency_hz=1.1e9, frequency_step_hz=1e7,
       right_ascension="12:00:00", declination="45:00:00", integration_time=10, observation_time="1h"
   )

   # Define a simple sky model (point source)
   source = PointSource(reference_intensity=1.0, sky_position="12:00:00 45:00:00")

   # Run the simulation (Dataset stays in memory until you write)
   sim = Simulator(interferometer=interferometer, sources=source)
   dataset = sim.simulate(create_dataset=True)

   # Write to disk — this is where you set the Measurement Set name/path
   ms_io = DaskMS(output_name="/path/to/simulated.ms")
   ms_io.write(dataset=dataset, tables="ALL")

.. tip::
   See :doc:`simulation` for more advanced simulation options, including Gaussian sources, composite models, and noise injection.

---

Simulating Non-Parametric Sources
---------------------------------

Pyralysis allows simulating **non-parametric sources**, meaning we compute visibilities
from an image without assuming any predefined model.

**Methods for estimating visibilities from an image:**

- **Nearest Neighbor Interpolation** – Fastest but least accurate.
- **Bilinear Interpolation** – Uses linear interpolation for smoother results.
- **Degridding (Convolutional Gridding with a Half-Kernel)** – Most accurate method.

**Example: Simulating Model Visibilities**

.. code-block:: python

   from pyralysis.estimators import NearestNeighbor, BilinearInterpolation, Degridding
   from pyralysis.convolution import PSWF1

   # Use Degridding (Convolutional Gridding)
   kernel = PSWF1(size=3, cellsize=0.003, oversampling_factor=3)
   mv_degridding = Degridding(input_data=dataset, image=image, cellsize=0.003, ckernel_object=kernel)

   # Compute visibilities
   mv_degridding.transform()

---

Gridding Visibilities
----------------------

Gridding is essential for transforming irregularly sampled visibilities into a **uniform Fourier grid**.
Pyralysis provides **DirtyMapper** to compute a **dirty image** and **dirty beam**.

**Example: Applying Gridding**

.. code-block:: python

   from pyralysis.transformers import DirtyMapper
   from pyralysis.convolution import PSWF1

   # Define a convolution kernel
   kernel = PSWF1(size=3, cellsize=0.003, oversampling_factor=3)

   # Create a DirtyMapper object for gridding
   dirty_mapper = DirtyMapper(input_data=dataset, imsize=512, cellsize=0.003, padding_factor=1.2, ckernel_object=kernel)

   # Compute dirty image and dirty beam
   dirty_images, dirty_beam = dirty_mapper.transform()

Both ``dirty_images`` and ``dirty_beam`` are returned as :class:`~pyralysis.reconstruction.image.Image` objects.

🔹 **For more details on gridding, see** :doc:`gridding`. For image processing and I/O operations, see :doc:`image_processing` and :doc:`io_operations`.

---

Performing Image Reconstruction
-------------------------------

To reconstruct an image, Pyralysis uses **gradient-based optimization techniques**.

**Example: Running an Optimizer**

.. code-block:: python

   from pyralysis.optimization.optimizer import HagerZhang
   from pyralysis.optimization.terms import ChiSquared, L1Norm
   from pyralysis.optimization import ObjectiveFunction

   # Define the objective function
   term_list = [ChiSquared(model_visibility=mv, normalize=True), L1Norm(penalization_factor=0.005)]
   of = ObjectiveFunction(term_list=term_list, image=image, persist_gradient=True)

   # Run the optimizer
   optim = HagerZhang(image=image, objective_function=of, max_iter=100)
   reconstructed_image = optim.optimize()

🔹 **For details on optimizers, line searchers, and seeders, see** :doc:`optimizers`,
:doc:`linesearchers`, and :doc:`step_size_seeders` (overview: :doc:`optimization`).

---

Phase Center Shifting
--------------------

Pyralysis provides the **Shifter** transformer for changing the effective phase center of observations. This is essential for:

- **Multi-field imaging** where different fields have different phase centers
- **Phase center corrections** for proper source positioning
- **Wide-field imaging** with non-coplanar baselines
- **Source positioning** in visibility space

**Example: Shifting Phase Center**

.. code-block:: python

   from pyralysis.transformers import Shifter
   from astropy.coordinates import SkyCoord
   import astropy.units as u

   # Shift to a new phase center
   new_center = SkyCoord(ra=12.5*u.deg, dec=-45.2*u.deg)
   shifter = Shifter(phase_center=new_center)
   shifter.transform()

   # Shift relative to a specific reference
   ref_center = SkyCoord(ra=12.0*u.deg, dec=-45.0*u.deg)
   shifter = Shifter(phase_center=new_center, reference_phase_center=ref_center)
   shifter.transform()

The phase shift is applied using the full UVW coordinate system:

.. math::
   V'(u,v,w) = V(u,v,w) \cdot e^{-2\pi i (u l_0 + v m_0 + w (n_0 - 1))}

where :math:`(l_0, m_0, n_0)` are the direction cosines of the phase center offset.

.. tip::
   The Shifter automatically handles visibility-like columns (``DATA``, ``MODEL_DATA``,
   ``CORRECTED_DATA``, ``RESIDUAL_DATA``, and other ``*_DATA`` columns) and updates
   both the visibility data and field phase center information.

---

Saving Reconstructed Images
---------------------------

Pyralysis provides **I/O handlers** to save reconstructed images.

**Example: Saving Images in FITS or Zarr Format**

.. code-block:: python

   from pyralysis.io import FITS, ZarrArray

   ioh_fits = FITS()  # FITS format
   ioh_zarr = ZarrArray()  # Zarr format

   # Use the chosen I/O handler when running the optimizer
   optim = HagerZhang(image=image, objective_function=of, max_iter=100, io_handler=ioh_fits)

🔹 **For details on I/O handlers, see** :py:mod:`pyralysis`.

---

Next Steps
----------

Now that you've seen the basic usage of Pyralysis, explore the advanced topics:

- **Gridding techniques:** :doc:`gridding`
- **Optimization methods:** :doc:`optimization`
- **API reference:** :py:mod:`pyralysis`

See also
--------
- :doc:`quickstart`
- :doc:`gridding`
- :doc:`optimization`

----

:doc:`quickstart` | :doc:`optimization`