Refactor sampling and reporting modules for improved maintainability

# Refactoring Plan: Monolithic `sampler.py` to `sampling/` Package

## Goal

Decompose the monolithic `PyAETNA/sampler.py` into a granular, modular package structure to improve maintainability, IDE performance, and AI agent context efficiency.

## Target Architecture

```text

PyAETNA/

└── sampling/

    ├── __init__.py          # Public API (exposes runSampler_NumPyro)

    ├── precision.py         # JAX precision, caching, and dtype casting

    ├── data.py              # PlData and column validation

    ├── model/               # Sub-package for model assembly

    │   ├── __init__.py

    │   ├── core.py          # PlModel class definition

    │   ├── detrending.py    # Polynomial and Fourier model logic

    │   └── gp_models.py     # GP-specific model methods

    ├── config/              # Sub-package for configuration

    │   ├── __init__.py

    │   ├── parser.py        # YAML parsing and scoped param expansion

    │   └── priors.py        # Prior string parsing and NumPyro dist mapping

    ├── engine/              # Sub-package for the sampling loop

    │   ├── __init__.py

    │   ├── runner.py        # runSampler_NumPyro (the main loop)

    │   ├── numpyro_model.py # _build_numpyro_model and log_joint

    │   └── logging.py       # All sampling-related log formatters

    └── utils.py             # Math helpers, roll-angle, and general utilities

```

## Execution Strategy: "Incremental Migration"

To ensure stability and prevent regression, the refactoring will follow a strict **one-by-one** submodule migration process.

**Baseline Requirement:**

`runPyAETNA.py` currently runs flawlessly with the monolithic `sampler.py` implementation. It **must** continue to do so after every single submodule is created. No additional code, logic changes, or feature additions are permitted—only the movement of code, correct imports, and maintaining existing function calling signatures.

### The Migration Loop

For each submodule defined in the plan:

1. **Extract**: Move the specific functions/classes listed below from `sampler.py` to the new submodule file.

2. **Disable**: Comment out or remove the original implementation in `sampler.py`.

3. **Bridge**: Import the function/class from the new submodule back into `sampler.py` so that the rest of the monolithic file (and external callers) still works.

4. **Validate**: Run `runPyAETNA.py` to verify that the system still functions correctly.

5. **Proceed**: Only after `runPyAETNA.py` passes will we move to the next submodule.

---

## Detailed Migration Sequence

### Phase 1: Foundation & Utilities COMPLETE 

1. **`sampling/utils.py`** COMPLETE

   - Functions to move: `linearizeRollAngle`, `_safe_eval_math`, `relabel_array`, `_get_default_param_order`.

2. **`sampling/precision.py`** COMPLETE

   - Logic to move: JAX cache configuration (`_JAX_CACHE_DIR` and `jax.config.update` calls).

   - Functions to move: `_configure_jax_precision`, `_cast_warmup_state_dtype`. 

3. **`sampling/data.py`** COMPLETE

   - Classes to move: `PlData`.

   - Functions to move: `_validate_and_prepare_columns`.

### Phase 2: Model Logic (`sampling/model/`) COMPLETE

4. **`sampling/model/detrending.py`** COMPLETE

   - Functions to move: `get_polynomial_model`, `_compute_polynomial_model` (static), `get_fourier_model`, `_compute_fourier_model` (static).

5. **`sampling/model/gp_models.py`** COMPLETE

   - Functions to move: `prebuild_gp_models`, `get_gp_model_cached`, `get_gp_marginal_logprob`, `_compute_gp_marginal_logprob` (static), `get_gp_posterior_prediction`, `_compute_gp_posterior_prediction` (static), `get_gp_model`, `_compute_gp_model` (static).

6. **`sampling/model/core.py`** COMPLETE

   - Classes to move: `PlModel` (base class and `__init__`).

   - *Note: Moved last in this phase so that the methods extracted above can be easily referenced or integrated into the core class.*

### Phase 3: Configuration (`sampling/config/`) COMPLETE

7. **`sampling/config/priors.py`** COMPLETE

   - Functions to move: `_parse_prior_string`, `_get_numpyro_dist`.

8. **`sampling/config/parser.py`** COMPLETE

   - Functions to move: `parse_config`, `_parse_yaml_config`, `_expand_detrend_coeff`, `_parse_scoped_param`, `_expand_scoped_param`.

### Phase 4: Sampling Engine (`sampling/engine/`) COMPLETE

9. **`sampling/engine/logging.py`** COMPLETE

   - Functions to move: `_format_sampling_message`, `_format_rate`, `_log_warmup_energy_summary`, `_log_warmup_summary`, `_log_config_summary`.

10. **`sampling/engine/numpyro_model.py`** COMPLETE

    - Functions to move: `_build_numpyro_model`, `_compute_log_joint`.

11. **`sampling/engine/runner.py`** COMPLETE

    - Functions to move: `runSampler_NumPyro`.

### Phase 5: Finalization

12. **`sampling/__init__.py`** COMPLETE

    - Expose `runSampler_NumPyro` and other key public classes.

13. **Cleanup** COMPLETE

    - Remove the now-empty `sampler.py` (or replace it with a simple redirect).

    - Update imports in `runPyAETNA.py` and other scripts to point directly to `PyAETNA.sampling`. Consider centralizing repeated imports across submodules.

14. **New subpackaging** 

    - Consider the new layering of submodules below runner.py

    - Consider the new layering of submodules below parser/__init__.py

15. Update the repository structure in README.md paying attention to the new subpackaging of the sampling and reporting stages (sampling/ and reporting/ subfolders respectively)

16. Upgrade versioning by mentioning the new submodular structure for the sampling and reporting stages

## Success Criteria

- `runPyAETNA.py` without any other arguments executes without errors.

- No change in sampling results or performance.

- All files in `PyAETNA/sampling/` are ideally under 600 lines.

## Version 3.6.0 (unreleased)

### Architecture Refactor — Modular Sampling and Reporting Packages

This release finalizes the structural refactor of the runtime pipeline into

dedicated subpackages for both sampling and reporting, improving

maintainability, IDE navigation, and module-level testability without changing

the scientific model semantics.

- **Sampling package layout:** the former monolithic sampler implementation is

  now organized under `PyAETNA/sampling/` with explicit layers for

  `model/`, `config/` (including `config/parser/`), and `engine/`

  (including `engine/utils/`).

- **Reporting package layout:** report generation and plotting helpers are now

  organized under `PyAETNA/reporting/` as focused modules (confidence bands,

  grid building, diagnostics, periodograms, residual plotting, PDF assembly,

  and related utilities).

- **Versioning note:** this is an architectural modularization release for the

  sampling and reporting stages; interfaces are preserved where possible and no

  intentional change to inference behavior is introduced by the package split.

### Transit Timing Variations (TTV)

Added full TTV support for photometric transit fitting:

def _get_passband_data(instrument_name, *, passband_svo=None, cache_dir=None,

                       log_func=None):

                       log_func=None, log_target=None):

    """Load passband throughput for an instrument.

    Resolution order:

        cache_path = os.path.join(cache_dir, f'{safe_name}.passband.dat')

        if os.path.isfile(cache_path):

            raw = np.loadtxt(cache_path)

            _log(f"  Passband '{svo_id}' loaded from cache: {cache_path}")

            _log(log_target, f"  Passband '{svo_id}' loaded from cache: {cache_path}")

            return raw[:, 0], raw[:, 1]

    # --- Download from SVO via urllib (plain HTTP to avoid expired-cert issues) ---

                np.savetxt(cache_path, np.column_stack([wl_A, tr]),

                           header=f'SVO filter: {svo_id}\nWavelength[A]  Throughput',

                           fmt='%.6f')

                _log(f"  Passband '{svo_id}' downloaded from SVO and cached: {cache_path}")

                _log(log_target, f"  Passband '{svo_id}' downloaded from SVO and cached: {cache_path}")

            else:

                _log(f"  Passband '{svo_id}' downloaded from SVO (not cached)")

                _log(log_target, f"  Passband '{svo_id}' downloaded from SVO (not cached)")

            return wl_A, tr

        except Exception as exc:

            _log(f"  WARNING: SVO download for '{svo_id}' failed ({exc}); "

            _log(log_target, f"  WARNING: SVO download for '{svo_id}' failed ({exc}); "

                 f"no passband available for '{instrument_name}'.")

    return None

def _download_model_spectrum(teff, logg, metallicity, cache_dir,

                             log_func=None):

                             log_func=None, log_target=None):

    """Download a stellar model spectrum (BT-Settl preferred, PHOENIX as fallback).

    Sources are tried in order:

    z_str  = f'{z_s:+.1f}'    # e.g. '+0.0', '-0.5'

    if teff_s != teff or logg_s != logg or z_s != metallicity:

        _log(f'  Model spectrum grid snap (PHOENIX grid used for cache key): requested '

        _log(log_target, f'  Model spectrum grid snap (PHOENIX grid used for cache key): requested '

             f'(Teff={teff:.0f}, logg={logg:.2f}, [Fe/H]={metallicity:+.2f}) '

             f'→ cache key '

             f'(Teff={teff_s}, logg={logg_s:.2f}, [Fe/H]={z_s:+.1f})')

            _first_line = _fh.readline()

        raw = np.loadtxt(cache_path)

        _src = 'btsettl' if 'BT-Settl' in _first_line else 'phoenix'

        _log(f'  Stellar spectrum ({_src}) loaded from cache: {cache_path}')

        _log(log_target, f'  Stellar spectrum ({_src}) loaded from cache: {cache_path}')

        return raw[:, 0], raw[:, 1], float(teff_s), float(logg_s), float(z_s), _src

    # ── Attempt 1: BT-Settl via SVO SSA (wider Teff/logg coverage, plain HTTP) ──────

    _log(f'  Trying BT-Settl via SVO '

    _log(log_target, f'  Trying BT-Settl via SVO '

         f'(Teff={teff_s} K, logg={logg_s:.2f}, [Fe/H]={z_s:+.1f}) …')

    bts_result = _download_btsettl_spectrum_via_svo(

        teff_s, logg_s, z_s, cache_dir, log_func=_log)

        teff_s, logg_s, z_s, cache_dir, log_func=_log, log_target=log_target)

    if bts_result is not None:

        wl_A, flux, best_teff, best_logg, best_meta = bts_result

        return wl_A, flux, best_teff, best_logg, best_meta, 'btsettl'

    # ── Attempt 2: PHOENIX LoRes from Göttingen (HTTP then HTTPS) ────────────────

    _log(f'  BT-Settl unavailable; trying PHOENIX from Göttingen '

    _log(log_target, f'  BT-Settl unavailable; trying PHOENIX from Göttingen '

         f'(Teff={teff_s} K, logg={logg_s:.2f}, [Fe/H]={z_s:+.1f}) …')

    # Build URL (Husser+ 2013, Göttingen PHOENIX-ACES-AGSS-COND-2011 LoRes)

    # Try HTTP first (more reliable in networks with proxy/firewall), fall back to HTTPS

        try:

            with urllib.request.urlopen(attempt_url, timeout=120) as resp:

                raw_bytes = resp.read()

            _log(f'  Successfully downloaded PHOENIX spectrum')

            _log(log_target, f'  Successfully downloaded PHOENIX spectrum')

            break

        except Exception as exc:

            _log(f'    Attempt {attempt_url[:45]}… failed: {exc}')

            _log(log_target, f'    Attempt {attempt_url[:45]}… failed: {exc}')

    if raw_bytes is None:

        _log(f'  WARNING: Could not download any model spectrum '

        _log(log_target, f'  WARNING: Could not download any model spectrum '

             f'(BT-Settl via SVO and PHOENIX from Göttingen both unavailable). '

             f'Falling back to blackbody.')

        _log(f'    Tip: pre-cache a spectrum at: {cache_path}')

        _log(log_target, f'    Tip: pre-cache a spectrum at: {cache_path}')

        return None, None, float(teff_s), float(logg_s), float(z_s), None

    try:

                f'(CRVAL1={crval}, CDELT1={cdelt}).'

            )

    except Exception as exc:

        _log(f'  WARNING: Could not parse PHOENIX FITS data ({exc}). '

        _log(log_target, f'  WARNING: Could not parse PHOENIX FITS data ({exc}). '

             f'Falling back to blackbody.')

        return None, None, float(teff_s), float(logg_s), float(z_s), None

        ),

        fmt='%.6e',

    )

    _log(f'  PHOENIX spectrum cached: {cache_path}')

    _log(log_target, f'  PHOENIX spectrum cached: {cache_path}')

    return wl_A, flux, float(teff_s), float(logg_s), float(z_s), 'phoenix'

def _download_btsettl_spectrum_via_svo(teff, logg, metallicity, cache_dir,

                                       log_func=None):

    """Fetch a BT-Settl spectrum from the SVO Theoretical Spectra web server.

                                           log_func=None, log_target=None):

    """

    Uses the IVOA Simple Spectral Access (SSA) interface at

    ``svo2.cab.inta-csic.es`` via plain HTTP (no SSL required).  The catalogue

    metadata (~2.5 MB VOTable for all 13 722 BT-Settl models) is fetched once;

    cache_path  = os.path.join(cache_dir, cache_fname)

    if os.path.isfile(cache_path):

        raw = np.loadtxt(cache_path)

        _log(f'  BT-Settl spectrum loaded from cache: {cache_path}')

        _log(log_target, f'  BT-Settl spectrum loaded from cache: {cache_path}')

        return raw[:, 0], raw[:, 1], float(teff), float(logg), float(metallicity)

    _SVO_SSA = 'http://svo2.cab.inta-csic.es/theory/newov2/ssap.php'

    # ── Step 1: fetch the BT-Settl model catalogue (metadata only) ─────────

    cat_url = f'{_SVO_SSA}?model=bt-settl&REQUEST=queryData'

    _log(f'  Fetching BT-Settl catalogue from SVO…')

    _log(log_target, f'  Fetching BT-Settl catalogue from SVO…')

    try:

        with urllib.request.urlopen(cat_url, timeout=60) as resp:

            cat_bytes = resp.read()

    except Exception as exc:

        _log(f'  BT-Settl catalogue fetch failed: {exc}')

        _log(log_target, f'  BT-Settl catalogue fetch failed: {exc}')

        return None

    try:

    build_gp_kernel_for_plotting,

    get_gp_x_variable,

)

from PyAETNA.model_utils import build_mean_function

from PyAETNA.model_utils import build_mean_function, detrend_data

from PyAETNA.body_models import (

    compute_lightcurve as _compute_lightcurve_raw,

    compute_lightcurve_supersampled as _compute_lightcurve_ss,

import jax

from . import sampler

from . import sampling

from .reporting import generate_reports

            random_seed = int(time.time() * 1000) % (2**31)

        try:

            sampler.runSampler_NumPyro(

            sampling.runSampler_NumPyro(

                config_yaml_abs, dataFile, logFilePath,

                outDir='PyAETNAouts/', num_chains=mcmc_config.get('num_chains', None),

                random_seed=random_seed,