Power-spectrum estimation
=========================

Building an estimator
---------------------

:class:`~ps_eor.pspec.PowerSpectraBuilder` constructs an estimator whose
geometry matches a Cartesian or spherical cube. Its configuration controls
the frequency transform, UV selection and binning, primary beam, weighting,
and optional exclusion of wedge or low-``k_parallel`` modes:

.. code-block:: python

   import numpy as np

   from ps_eor import pspec

   config = pspec.PowerSpectraConfig(
       window_fct="blackmanharris",
       ft_method="lssa",
   )
   config.umin = 50
   config.umax = 250
   config.du = 10

   ps_gen = pspec.PowerSpectraBuilder(config).get(cube)
   kbins = np.logspace(np.log10(ps_gen.kmin), np.log10(0.5), 10)

The estimator can use the full frequency band, a named
:class:`~ps_eor.pspec.EorBin`, an explicit MHz range, or a band centred on a
redshift. The selected bandwidth determines the line-of-sight resolution and
cosmological normalization.

.. important::

   The default ``primary_beam`` is ``"lofar_hba"`` and the beam enters the
   power-spectrum normalization. Set the appropriate instrument beam for
   non-LOFAR data, or use ``"no_pb"`` when a unit beam is intended.

Common configuration parameters
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. list-table::
   :header-rows: 1
   :widths: 22 40 22 16

   * - Parameter
     - What it controls
     - Units or values
     - Default
   * - ``window_fct``
     - Spectral window applied before the frequency-to-delay transform.
     - Any window accepted by ``scipy.signal.get_window``
     - ``"hann"``
   * - ``ft_method``
     - Frequency transform used to form the cylindrical spectrum.
     - ``"nudft"`` or ``"lssa"``
     - ``"nudft"``
   * - ``rmean_freqs``
     - Whether to remove the frequency mean before the transform.
     - Boolean
     - ``False``
   * - ``ps2d_pos_only``
     - Whether the cylindrical result retains only non-negative
       :math:`k_\parallel`.
     - Boolean
     - ``True``
   * - ``umin``, ``umax``
     - Baseline range included in the estimate.
     - Wavelengths
     - ``50``, ``250``
   * - ``du``
     - Width of transverse UV bins.
     - Wavelengths
     - ``10``
   * - ``uniform_u_bins``
     - Selects uniformly spaced UV-bin centres instead of the explicit
       ``el`` grid.
     - Boolean
     - ``False``
   * - ``weights_by_default``
     - Whether ``weighted="default"`` uses the cube's UV weights.
     - Boolean
     - ``True``
   * - ``empirical_weighting``
     - Enables additional weighting derived from measured noise.
     - Boolean
     - ``False``
   * - ``filter_kpar_min``
     - Excludes lower :math:`|k_\parallel|` modes during spherical averaging.
     - :math:`\mathrm{h\,cMpc^{-1}}` or ``None``
     - ``None``
   * - ``filter_wedge_theta``
     - Excludes modes below the foreground-wedge line during spherical
       averaging.
     - Degrees
     - ``0``
   * - ``kbins_kmax``
     - Upper limit of automatically generated spherical k bins.
     - :math:`\mathrm{h\,cMpc^{-1}}`
     - ``0.6``
   * - ``kbins_n``
     - Number of logarithmic spherical k-bin edges.
     - Integer
     - ``6``
   * - ``primary_beam``
     - Primary-beam model used in the power-spectrum normalization.
     - Registered beam name
     - ``"lofar_hba"``

These values can be stored in a
:class:`~ps_eor.pspec.PowerSpectraConfig` file or overridden for one estimator
through :meth:`~ps_eor.pspec.PowerSpectraBuilder.get`. The
:doc:`Power-spectrum API <../api/pspec>` lists the remaining settings.

Available products
------------------

Request an individual product when that is all the analysis needs:

.. code-block:: python

   ps = ps_gen.get_ps(cube)
   ps2d = ps_gen.get_ps2d(cube)
   ps3d = ps_gen.get_ps3d(kbins, cube)
   variance = ps_gen.get_variance(cube)

Use ``get_all`` to calculate all standard products together:

.. code-block:: python

   products = ps_gen.get_all(kbins, cube)

   products.ps.plot()        # spatial P(k_perp) versus frequency
   products.ps2d.plot()      # cylindrical P(k_perp, k_parallel)
   products.ps3d.plot()      # spherical Delta^2(k)
   products.variance.plot()  # variance versus frequency

``get_all`` is a convenient shorthand when all four products are wanted;
calling the individual methods is also fast and is often clearer when only a
subset is needed.

Results are stored in Kelvin squared. Their ``get`` and ``plot`` methods can
present values in mK squared, or in mK for field amplitudes. Wavenumbers use
``h cMpc^-1``.

Compatible power-spectrum products also support addition, subtraction, and
scalar multiplication. For example, ``0.5 * (ps_a + ps_b)`` averages two
products and ``ps_data - ps_noise`` forms a noise-debiased estimate. The
coordinates and product types must match; uncertainties and weights are
propagated by the result class.

A four-panel simulated example
------------------------------

:class:`~ps_eor.pspec.FourPanelPsResults` compares several cubes in four
representations: variance versus frequency, spherical
:math:`\Delta^2(k)`, cylindrical power averaged over
:math:`k_\perp`, and cylindrical power averaged over
:math:`k_\parallel`.

The figure below is generated from three independent components drawn with
:class:`~ps_eor.simu.SimuMultiGPCube`: a smooth foreground covariance, a
shorter-coherence 21-cm-like covariance, and white noise.

.. image:: /_static/power_spectra_four_panel.png
   :alt: Four panels comparing simulated foreground, 21-cm-like signal, and noise power spectra.
   :width: 100%

The complete figure-generation script is included in the documentation. It
uses the covariance-kernel simulator and therefore requires the ``ml-gpr``
installation extra.

.. raw:: html

   <details>
   <summary><strong>Show the figure-generation script</strong></summary>

.. literalinclude:: ../examples/power_spectra_four_panel.py
   :language: python

.. raw:: html

   </details>

Auto, cross, and noise-debiased spectra
---------------------------------------

An auto spectrum contains the noise bias of its input. If a matching noise
realization or expectation is available,
``get_ps3d_with_noise`` subtracts its estimated power:

.. code-block:: python

   debiased = ps_gen.get_ps3d_with_noise(kbins, data_cube, noise_cube)

For statistically independent data splits, cross power avoids the mean
uncorrelated-noise bias:

.. code-block:: python

   cross = ps_gen.get_cross_ps3d(kbins, split_a, split_b)

The two cubes must represent the same sky signal on compatible frequency and
spatial grids. Cross power does not remove correlated systematics; a negative
cross term can also cancel positive sky power in the total measured
cross-spectrum. Inspect cylindrical auto, cross, and coherence products before
interpreting a spherical average.

Weighting and excluded regions
------------------------------

By default, estimators use cube weights when they are present. Wedge and
minimum-``k_parallel`` selections act through the estimator weights during
spherical averaging; they do not erase the corresponding cells from the
cylindrical result. Keep these choices in the saved analysis configuration,
because changing them changes both the estimate and its effective number of
modes.