Regularization Methods
======================

Overview
--------
Pyralysis supports a variety of regularization (penalty) functions to encode prior knowledge or promote desirable properties in reconstructed images. Regularizers can be combined with data fidelity terms in the `ObjectiveFunction` for flexible, composite objectives.

.. note::
   All regularizers support masks, multi-channel images, and have efficient gradient and proximal operator implementations for use in modern optimization algorithms.

Available Regularizers
----------------------

L1 Norm Regularizer
~~~~~~~~~~~~~~~~~~~
Promotes sparsity in the image (many pixels exactly zero).

**Mathematical Formulation:**

.. math::

   \ell_1(I) = \sum_i \|I_i\|

**Typical Use Case:** Sparse solutions, e.g., point sources, compressed sensing.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import L1Norm

   l1_term = L1Norm(penalization_factor=0.01)

   # Apply to specific channel in multi-channel image
   l1_channel = L1Norm(parameter_index=0, penalization_factor=0.01)

Entropy Regularizer
~~~~~~~~~~~~~~~~~~~
Favors similarity to a prior image (maximum entropy method).

**Mathematical Formulation:**

.. math::

   S = \sum_i I_i \log(I_i / G)

**Typical Use Case:** Biases solution toward a reference/prior, useful in MEM.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import Entropy

   entropy_term = Entropy(
       prior_image=prior_model,
       penalization_factor=0.005
   )

Tikhonov Regularizer
~~~~~~~~~~~~~~~~~~~~
Penalizes large values (L2 norm squared).

**Mathematical Formulation:**

.. math::

   f(I) = \frac{1}{2}\|I\|_2^2

**Typical Use Case:** Smooths the solution, discourages large pixel values.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import Tikhonov

   tikhonov_term = Tikhonov(penalization_factor=0.001)

Total Flux Regularizer
~~~~~~~~~~~~~~~~~~~~~~
Controls the total sum of pixel values (flux).

**Mathematical Formulation:**

.. math::

   f(I) = \sum_k I_k

**Typical Use Case:** Enforces or penalizes total flux, e.g., for normalization.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import TotalFlux

   flux_term = TotalFlux(penalization_factor=0.1)

TSV (Total Squared Variation) Regularizer
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Promotes smoothness while preserving edges (quadratic TV).

**Mathematical Formulation:**

.. math::

   TSV(I) = \sum_i [ (\partial_x I_i)^2 + (\partial_y I_i)^2 ]

**Typical Use Case:** Edge-preserving smoothing, denoising, super-resolution.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import TSV

   tsv_term = TSV(penalization_factor=0.01)

L2 Norm Regularizer
~~~~~~~~~~~~~~~~~~~
Penalizes the L2 norm (Euclidean norm) of the image.

**Mathematical Formulation:**

.. math::

   \|I\|_2 = \sqrt{\sum_i I_i^2}

**Typical Use Case:** General smoothness, different from Tikhonov (which uses L2 squared).

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import L2Norm

   l2_term = L2Norm(penalization_factor=0.001)

**Note:** This is different from Tikhonov regularization, which uses the L2 squared norm
:math:`\frac{1}{2}\|I\|_2^2`. The L2 norm is non-differentiable at zero, while Tikhonov is
differentiable everywhere.

Huber Loss Regularizer
~~~~~~~~~~~~~~~~~~~~~~
Robust loss function that combines L1 and L2 behavior, providing a smooth transition.

**Mathematical Formulation:**

.. math::

   L_\delta(x) = \begin{cases}
       \frac{1}{2}x^2 & \text{if } |x| \leq \delta \\
       \delta(|x| - \frac{\delta}{2}) & \text{if } |x| > \delta
   \end{cases}

**Typical Use Case:** Robust optimization in the presence of outliers, handling calibration errors,
RFI mitigation in radio interferometry.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import Huber

   huber_term = Huber(delta=1.0, penalization_factor=0.01)

**Parameters:**

   - ``delta``: Threshold parameter determining L1/L2 transition (default: 1.0)
   - Smaller values: More L1-like (robust to outliers)
   - Larger values: More L2-like (smooth)

Isotropic Total Variation (TV) Regularizer
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
True Total Variation using Euclidean norm of gradients (rotationally invariant).

**Mathematical Formulation:**

.. math::

   TV_{iso}(I) = \sum_{i} \sqrt{ \left(\frac{\partial I}{\partial x}\right)_i^2 +
   \left(\frac{\partial I}{\partial y}\right)_i^2 }

**Typical Use Case:** Edge-preserving smoothing, image denoising, superior to TSV for preserving
sharp features.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import TVIsotropic

   tv_iso_term = TVIsotropic(penalization_factor=0.01, epsilon=1e-8)

**Parameters:**

   - ``epsilon``: Regularization parameter to avoid division by zero in gradient (default: 1e-8)

**Note:** This is the true non-smooth TV, different from TSV (smooth approximation). Requires
specialized optimization methods (FISTA, SDMM) due to non-differentiability.

Anisotropic Total Variation (TV) Regularizer
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Total Variation using L1 norm of gradients (separable, faster than isotropic).

**Mathematical Formulation:**

.. math::

   TV_{aniso}(I) = \sum_{i} \left| \frac{\partial I}{\partial x} \right|_i +
   \left| \frac{\partial I}{\partial y} \right|_i

**Typical Use Case:** Edge-preserving smoothing, faster than isotropic TV, good for images with
axis-aligned features.

**Implementation:**

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import TVAnisotropic

   tv_aniso_term = TVAnisotropic(penalization_factor=0.01)

**Note:** Anisotropic TV is separable and typically faster to compute than isotropic TV, but is
not rotationally invariant. Both TV variants are non-smooth and require proximal optimization
methods (FISTA, SDMM).

Regularizer Compatibility with Optimizers
-----------------------------------------

**Gradient-Based Optimizers (CG, L-BFGS):**
These optimizers work with smooth (differentiable) regularizers:

- **Entropy**, **Tikhonov**, **L2Norm**, **Huber**, **TotalFlux**, **TSV**

For gradient-based optimizers and line search, see :doc:`optimizers` and
:doc:`linesearchers` (overview: :doc:`optimization`).

**Proximal Optimizers (FISTA, SDMM):**
These optimizers work with regularizers that have proximal operators (both smooth and non-smooth):

- **All regularizers** implement proximal operators and can be used with FISTA/SDMM
- **Non-differentiable regularizers** (L1Norm, TVIsotropic, TVAnisotropic, L2Norm) require proximal methods
- **Smooth regularizers** can also be used with proximal methods for flexibility

For details on proximal optimization, see :doc:`compressed_sensing`.

Key Features
------------
- **Multi-channel Support**: All regularizers support multi-channel images with optional channel selection via `parameter_index`
- **Masking**: Full support for boolean masks to restrict operations to specific regions
- **Efficient Gradients**: Optimized gradient computations using Dask arrays for scalability
- **Proximal Operators**: All regularizers implement proximal operators for FISTA/SDMM compatibility
- **Flexible Penalization**: Configurable penalization factors for each regularizer
- **Optimized Operations**: All mask-based operations use efficient arithmetic instead of conditional branching

TSV Regularizer Details
-----------------------
The TSV (Total Squared Variation) regularizer is a smooth approximation of Total Variation that promotes edge-preserving smoothness. It computes finite differences along spatial dimensions and uses FFT-based proximal operators for efficient optimization.

**Advantages over Tikhonov:**
- Preserves sharp features while smoothing noise
- Better edge preservation
- More suitable for images with complex structures

**Mathematical Details:**
TSV computes finite differences along spatial dimensions and uses FFT-based proximal operators:

.. math::

   \nabla TSV(I) = D^T D\, I

where :math:`D` denotes the finite difference operator along the spatial dimensions.

Using Regularizers in Objective Functions
-----------------------------------------
Combine multiple regularizers with data fidelity terms:

.. code-block:: python

   from pyralysis.optimization import ObjectiveFunction
   from pyralysis.optimization.terms import Chi2, L1Norm, TSV, Entropy

   # Create objective function with multiple terms
   terms = [
       Chi2(model_visibility=mv, normalize=True),           # Data fidelity
       L1Norm(penalization_factor=0.01),                    # Sparsity
       TSV(penalization_factor=0.005),                      # Edge-preserving smoothness
       Entropy(prior_image=prior_model, penalization_factor=0.001)  # Prior bias
   ]

   obj_func = ObjectiveFunction(terms=terms, parameter=image)

Advanced Usage with Masks
-------------------------
Apply regularization only to specific regions:

.. code-block:: python

   from pyralysis.reconstruction.mask import Mask

   # Create mask for specific region
   mask = Mask(data=region_of_interest)

   # Apply regularization with mask
   l1_term = L1Norm(penalization_factor=0.01)
   l1_term.function(mask=mask)  # Only considers masked pixels

Multi-channel Support
---------------------
Handle spectral and polarization data efficiently:

.. code-block:: python

   # Apply to specific channel in multi-channel image
   l1_channel = L1Norm(parameter_index=0, penalization_factor=0.01)  # Channel 0 only

   # Apply to all channels
   l1_all = L1Norm(parameter_index=-1, penalization_factor=0.01)     # All channels

Performance Optimization
--------------------------
- **Gradient Persistence**: Enable ``persist_gradient=True`` in ObjectiveFunction for caching gradients
- **Proximal Persistence**: Enable ``persist_proximal=True`` in regularizer terms for caching proximal operator results (useful for FISTA/SDMM)
- **Chunking**: Configure Dask chunking for optimal memory usage
- **Masking**: Use masks to limit computation to relevant regions
- **Channel Selection**: Use `parameter_index` to avoid unnecessary computations

**Proximal Operator Persistence:**

For proximal optimization methods (FISTA, SDMM), enable persistence to maintain consistent iteration times:

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import L1Norm, Tikhonov

   # Enable proximal persistence
   l1_term = L1Norm(penalization_factor=0.01, persist_proximal=True)
   tikh_term = Tikhonov(penalization_factor=0.001, persist_proximal=True)

   # Proximal results are cached in memory, similar to persist_gradient

Best Practices
--------------
1. **Start Simple**: Begin with basic regularizers (L1, Tikhonov) before adding complexity
2. **Balance Terms**: Ensure data fidelity and regularization terms are properly balanced
3. **Use Masks**: Apply regularization only to regions where it's meaningful
4. **Monitor Convergence**: Watch for convergence issues that may indicate regularization conflicts
5. **Validate Results**: Always validate that regularization improves image quality

**Example: Progressive Regularization**

.. code-block:: python

   # Start with basic regularization
   basic_terms = [Chi2(mv), L1Norm(0.01)]

   # Add complexity gradually
   advanced_terms = [Chi2(mv), L1Norm(0.01), TSV(0.005)]

   # Full regularization suite
   full_terms = [Chi2(mv), L1Norm(0.01), TSV(0.005), Entropy(prior, 0.001)]

See Also
--------
- :doc:`optimizers` - Conjugate gradient and L-BFGS (smooth regularizers)
- :doc:`linesearchers` - Step acceptance along the search direction
- :doc:`step_size_seeders` - Barzilai–Borwein and interpolation initial steps
- :doc:`optimization` - Workflow overview
- :doc:`compressed_sensing` - Proximal methods (non-smooth regularizers)
- :doc:`objective_functions` - Objective function framework

----

:doc:`optimizers` | :doc:`compressed_sensing` | :doc:`objective_functions`