Multiwavelength (2D) Analysis
==============================

This guide explains how to fit GP models to light curves observed simultaneously in
multiple photometric bands or across a range of wavelengths.

.. contents:: On this page
   :local:
   :depth: 2

Overview
--------

When a source is observed in multiple bands, ``pgmuvi`` can fit a **2D Gaussian
process** whose inputs are ``(time, wavelength)``.  The model captures:

* the temporal variability structure (periods, noise) shared across bands,
* the wavelength dependence of variability amplitude and coherence,
* band-specific noise levels.

This is more informative than fitting each band independently because the model
borrows statistical strength across bands, particularly when some bands have sparse
coverage.

Data Format for 2D Models
--------------------------

For 2D / multiband fitting, ``xdata`` must have shape ``(N, 2)`` where:

* column 0 is the **observation time** (same unit across all observations),
* column 1 is the **wavelength or band identifier** and **must be numeric**.
  Use the effective wavelength (e.g., in microns or nanometres) for
  wavelength-dependent (separable 2D) kernels, or integer band codes
  (0, 1, 2 …) for generic categorical band identification.
  A physically meaningful numeric wavelength is strongly recommended for
  separable 2D kernels so that the kernel can learn the correct wavelength
  scaling.

  .. note::

     When loading from a CSV file via :meth:`~pgmuvi.lightcurve.Lightcurve.from_csv`,
     string band labels are mapped to integer codes only if the wavelength input
     column itself is string-typed, for example with ``wavelcol="band"`` or
     ``xcol=["time", "band"]``. In that case, the encoded values are used in
     ``xdata[:, 1]`` and the original labels (e.g. ``"V"``, ``"I"``) are stored in
     :attr:`~pgmuvi.lightcurve.Lightcurve.band`. If the CSV already provides a
     numeric wavelength column, that numeric column is used directly in
     ``xdata[:, 1]`` and any separate string band-ID column is only stored in
     :attr:`~pgmuvi.lightcurve.Lightcurve.band` without remapping. The
     :class:`~pgmuvi.lightcurve.Lightcurve` constructor itself still requires
     numeric ``xdata``; passing string arrays directly will raise an error.

Human-readable band labels (e.g. ``"V"``, ``"R"``) can be stored separately in
the :attr:`~pgmuvi.lightcurve.Lightcurve.band` attribute for labelling plots,
but they play no role in the GP computation.

All bands are stacked into a single array::

    import numpy as np
    import pgmuvi

    # Two bands: times_b0/b1, fluxes_b0/b1, errors_b0/b1
    times_all  = np.concatenate([times_b0,  times_b1])
    fluxes_all = np.concatenate([fluxes_b0, fluxes_b1])
    errors_all = np.concatenate([errors_b0, errors_b1])
    wavelengths = np.concatenate([
        np.full_like(times_b0, fill_value=0.55),   # band 0: 0.55 μm
        np.full_like(times_b1, fill_value=2.20),   # band 1: 2.20 μm
    ])

    xdata = np.column_stack([times_all, wavelengths])
    lc = pgmuvi.lightcurve.Lightcurve(xdata, fluxes_all, errors_all)

``pgmuvi`` detects that this is 2D input and sets up the light-curve object
accordingly, but you must still choose a 2D model explicitly when fitting
(for example, with ``fit(model="2D")`` or ``set_model("2D"); fit()``).

Fitting a 2D Model
-------------------

The fitting workflow is the same as in 1D::

    lc.fit(model="2D")

For heterogeneous band sampling (e.g., one band has far more observations than
others), consider using the best-band initialisation option::

    lc.fit(model="2D", use_best_band_init=True)

This uses a 1-D Lomb–Scargle on the most densely sampled band instead of the
multiband periodogram to seed the spectral-mixture kernel frequencies.  This
can substantially improve frequency initialisation when one band has far more
observations than the others.

Assessing Sampling Quality per Band
-------------------------------------

Because each band may have a different cadence, it is important to assess data quality
per band before fitting::

    lc.assess_sampling_quality_per_band()

You can then filter out bands that have insufficient coverage::

    lc = lc.filter_well_sampled_bands(min_points=20)

Visualising 2D Results
-----------------------

Use :meth:`~pgmuvi.lightcurve.Lightcurve.plot` to visualise multiband fits. For
2D data, it shows each band's observations together with the GP predictive mean and
credible interval.

The :meth:`~pgmuvi.lightcurve.Lightcurve.plot_psd` method does not currently
support 2D models; plotting PSDs in more than one dimension is not implemented.

Kernel Choices for 2D Models
------------------------------

The choice of kernel depends on which 2D model family you are using:

* **``model="2D"`` (non-separable):** Uses a 2D spectral-mixture kernel that
  treats time and wavelength jointly.  There is no separate "wavelength kernel"
  to configure; the model learns the joint covariance structure directly.

* **Separable model family** (``"2DSeparable"``, ``"2DAchromatic"``,
  ``"2DWavelengthDependent"``, etc.): Uses an explicit product kernel where the
  temporal spectral-mixture kernel and the wavelength kernel are independent
  factors.  Common choices for the wavelength kernel include:

  * **RBF (Squared-Exponential):** Smooth, continuous wavelength dependence.
  * **Matérn ν = 3/2:** Slightly rougher wavelength dependence, more robust to
    outliers.

  See the ``separable_kernels_2d.py`` and ``dust_mean_spectral_mixture_2d.py``
  example scripts for illustrations of different kernel configurations.

Separable vs Non-Separable Models
-----------------------------------

``pgmuvi`` supports more than one kind of 2D kernel.  The explicit
``2DSeparable`` model family uses **product (separable) kernels**, where the time
and wavelength covariance are modeled as separate factors.  In these models, the
separability assumption means that the temporal correlation structure is the same at
every wavelength (up to a scaling factor).  This is a reasonable approximation for
many sources but may break down when:

* the period changes significantly with wavelength (e.g., in accretion disk reverberation),
* the variability mechanism differs qualitatively between bands.

By contrast, ``model="2D"`` uses a 2D spectral-mixture kernel rather than an
explicit ``k_time * k_wavelength`` ProductKernel construction, so the separability
assumption above should only be interpreted as applying to the ``2DSeparable``
family.

In such cases, consider fitting each band independently and comparing the inferred
periods, or contact the ``pgmuvi`` developers to discuss extensions.

2D Tutorial
-----------

A dedicated notebook tutorial for 2D multiwavelength analysis is provided in the
User Guide:

.. toctree::
   :maxdepth: 1

   ../notebooks/pgmuvi_tutorial_2d