improved reporting

# Plan: Autodetect Data Clusters for Time-Series Plotting

## Objective

Modify the plotting logic in `plot_models.py` to automatically detect data clusters and create a "zoomed-in" view of these clusters. This applies **only** for per-instrument plots. This will eliminate uninformative gaps in the time-series plots while maintaining the overall figure structure and the 6:1 ratio between the left (time-series) and right (diagnostics/KDE) columns.

## Requirements

1. **Cluster Detection**: Use `detect_data_clusters` from `model_math.py` to identify contiguous blocks of data.

2. **Layout Preservation**: 

   - Keep the 6:1 ratio between the left and right columns.

   - The right column (KDE/diagnostics) remains a single panel per row.

   - Each row of the left column must be split into subpanels, one per detected cluster.

3. **Proportional Widths**: The width of each subpanel must be proportional to the timespan of the corresponding data cluster.

4. **Zoom-in Effect**: 

   - Iterate the current plotting procedure over these subpanels.

   - Set `xlim` for each subpanel to the timespan of its corresponding cluster, adding a buffer of 5% of the cluster's timespan to the left and right.

5. **Shared Y-Axis**: Subpanels within the same row must share the same y-axis.

6. **No Logic Modification**: The actual plotting commands (errorbars, lines, fills) must remain unchanged; only the axes they are plotted on and their limits are modified.

## Implementation Plan

At the end of each step, run `/dati/pyaetnaEnv/bin/python /dati/software/pyLib/PyAETNA/runPyAETNA.py` and do not proceed until all bugs are fixed.

### 1. Integration of Cluster Detection

#### Status: COMPLETE

- Import `detect_data_clusters` into `plot_models.py`.

- In `plot_models()`, after sorting the data by time (`t`), call `detect_data_clusters(t)` to get the list of cluster boundaries.

- Calculate the total timespan of all clusters to determine the proportional widths.

- *Note: While clusters are detected for any selection, the zoomed-in layout will be applied only when plotting per-instrument data (where gaps are most prevalent).*

### 2. Modification of Figure Layout

#### Status: COMPLETE

- The current `create_figure_layout` likely creates a grid of axes. This needs to be updated or wrapped.

- Instead of a single `axes[i]` for the left column, we will need a list of sub-axes for each row `i`.

- Use `matplotlib.gridspec.GridSpec` or `subplot_mosaic` to create a layout where:

    - The right column is a fixed width.

    - The left column is divided into $N$ sub-columns (where $N$ is the number of clusters), with widths proportional to the cluster durations.

### 3. Refactoring the Plotting Loop

#### Status: COMPLETE

- Wrap the existing plotting logic for each panel (Panel 0, Panel 1, Panel 2) in a loop over the detected clusters.

- This looping logic must be applied regardless of the number of clusters (including when only one cluster is detected) to ensure consistent code paths.

- For each cluster $C_j$:

    - Select the sub-axis $Ax_{i,j}$ corresponding to row $i$ and cluster $j$.

    - Execute the existing plotting code on $Ax_{i,j}$.

    - Set `Ax_{i,j}.set_xlim(C_j.start, C_j.end)`.

    - Ensure `sharey` is enabled across $Ax_{i,j}$ for all $j$ in row $i$.

### 4. Handling Specific Panels

#### Status: COMPLETE

- **Panel 0 (Raw Data)**: Apply the same cluster-splitting logic if it is visible.

- **Panel 1 (Detrended)**: Split into subpanels per cluster.

- **Panel 2 (Residuals)**: Split into subpanels per cluster.

- **KDE/Right Column**: Remains as a single panel per row, as it does not depend on the time-axis.

### 5. Axis Formatting

- Update `format_plot_axes` to handle the new multi-axis structure.

- Ensure x-axis labels and tickmarks are only shown on the bottom-most row of subpanels.

- Ensure y-axis labels and tickmarks are only shown on the left-most column of subpanels.

- Ensure the the y ranges are shared only per-row, not all over the subpanels

- Remove redundant y-axis labels for internal subpanels.

- Promote the title of the [0,0] panel to the title of the whole page

## Technical Details

- **Cluster Widths**: If clusters are $[(t_{0,start}, t_{0,end}), \dots, (t_{N,start}, t_{N,end})]$, the width ratio for `GridSpec` will be $[(t_{0,end}-t_{0,start}), \dots, (t_{N,end}-t_{N,start})]$.

- **Y-Axis Sharing**: Use `sharey=axes[i, 0]` when creating sub-axes for row $i$.

# Direct layout creator used by plot_models

# ──────────────────────────────────────────────────────────────────────────────

def create_figure_layout(figsize, show_panel0):

def create_figure_layout(figsize, show_panel0, clusters=None):

    """Create the figure and panel axes for plot_models.

    Standard layout: 3 left-column panels (overlay / detrended / residuals)

    ``show_panel0`` is True, a 4th "raw data" panel is prepended, giving a

    taller figure.

    ``axes[0..2]`` always refer to overlay / detrended / residuals so all

    downstream rendering code is independent of whether Panel 0 is shown.

    ``ax_raw_panel`` is ``None`` in the 3-panel case.

    If ``clusters`` is provided (list of (start, end) tuples), the left column

    is split into sub-panels, one per cluster, with widths proportional to

    the cluster duration.

    Parameters

    ----------

        Base figure size.

    show_panel0 : bool

        Whether to add the extra raw-data panel at the top.

    clusters : list of tuples, optional

        List of (start, end) for each data cluster.

    Returns

    -------

    fig : matplotlib.figure.Figure

    axes : list[Axes]

        ``[ax_overlay, ax_detrended, ax_residuals]``

    axes : list[list[Axes]]

        Nested list: ``[[ax_overlay_c1, ...], [ax_detrended_c1, ...], [ax_residuals_c1, ...]]``

    ax_kde : Axes

        Right-inset KDE axis (shares y-axis with ``axes[2]``).

    ax_raw_panel : Axes or None

        Extra top panel (raw data + non-time correction), ``None`` if

        ``show_panel0`` is False.

    """

        Right-inset KDE axis (shares y-axis with the last residual sub-panel).

    ax_raw_panel : list[Axes] or Axes or None

        If clusters: list of sub-panels for raw data.

        If no clusters: single Axes or None.

    """,

    from matplotlib.gridspec import GridSpec as _GS

    # ── Handle Cluster-based Layout ──────────────────────────────────────

    if clusters is not None and len(clusters) > 1:

        n_clusters = len(clusters)

        # Widths proportional to cluster durations

        widths = [c[1] - c[0] for c in clusters]

        # Total width for left column is 6, right is 1.

        # We use a nested GridSpec or a complex width_ratios.

        # Simplest: create a GridSpec with n_clusters + 1 columns.

        # The first n_clusters columns are the zoomed views, the last is KDE.

        # Normalize widths to sum to 6

        total_w = sum(widths)

        norm_widths = [w * 6.0 / total_w for w in widths]

        full_width_ratios = norm_widths + [1.0]

        if show_panel0:

            _n_rows = 4

            fig = plt.figure(figsize=(figsize[0], figsize[1] * 1.25))

            _gs = _GS(_n_rows, n_clusters + 1, figure=fig, width_ratios=full_width_ratios,

                      hspace=0.08, wspace=0.05)

            ax_raw_panel = []

            for col in range(n_clusters):

                ax_raw_panel.append(fig.add_subplot(_gs[0, col]))

            # Each column's rows share x only within the same column so that

            # set_xlim on one cluster does not propagate to other clusters.

            ax0 = []

            for col in range(n_clusters):

                ax0.append(fig.add_subplot(_gs[1, col], sharex=ax_raw_panel[col]))

            ax1 = []

            for col in range(n_clusters):

                ax1.append(fig.add_subplot(_gs[2, col], sharex=ax_raw_panel[col],

                                           sharey=ax1[0] if col > 0 else None))

            ax2 = []

            for col in range(n_clusters):

                ax2.append(fig.add_subplot(_gs[3, col], sharex=ax_raw_panel[col],

                                           sharey=ax2[0] if col > 0 else None))

            ax_kde = fig.add_subplot(_gs[3, n_clusters], sharey=ax2[0])

        else:

            _n_rows = 3

            fig = plt.figure(figsize=figsize)

            _gs = _GS(_n_rows, n_clusters + 1, figure=fig, width_ratios=full_width_ratios,

                      hspace=0.08, wspace=0.05)

            ax_raw_panel = None

            ax0 = []

            for col in range(n_clusters):

                ax0.append(fig.add_subplot(_gs[0, col]))

            # Each column's rows share x only within the same column.

            ax1 = []

            for col in range(n_clusters):

                ax1.append(fig.add_subplot(_gs[1, col], sharex=ax0[col],

                                           sharey=ax1[0] if col > 0 else None))

            ax2 = []

            for col in range(n_clusters):

                ax2.append(fig.add_subplot(_gs[2, col], sharex=ax0[col],

                                           sharey=ax2[0] if col > 0 else None))

            ax_kde = fig.add_subplot(_gs[2, n_clusters], sharey=ax2[0])

        axes = [ax0, ax1, ax2]

        return fig, axes, ax_kde, ax_raw_panel

    # ── Standard Single-Column Layout (Fallback) ──────────────────────────

    if show_panel0:

        _n_rows = 4

        fig = plt.figure(figsize=(figsize[0], figsize[1] * 1.25))

    time_center = float(getattr(plData, 'timeCenter', 0.0))

    if time_center != 0.0:

        fmt = FuncFormatter(lambda x, pos: f"{(x + time_center):.5f}")

        _all_axes = [ax0, ax1, ax2]

        # Flatten all axes into a single list to apply the formatter

        _all_axes_flat = []

        if ax_raw_panel is not None:

            _all_axes.insert(0, ax_raw_panel)

        for ax in _all_axes:

            _all_axes_flat.extend(ax_raw_panel if isinstance(ax_raw_panel, list) else [ax_raw_panel])

        for ax_group in [ax0, ax1, ax2]:

            _all_axes_flat.extend(ax_group if isinstance(ax_group, list) else [ax_group])

        for ax in _all_axes_flat:

            ax.xaxis.set_major_formatter(fmt)

    time_zero_point = float(getattr(plData, 'timeZeroPoint', 0.0))

    # Apply x-label to the first sub-panel of the residuals row

    _ax2_first = ax2[0] if isinstance(ax2, list) else ax2

    if time_zero_point != 0.0:

        ax2.set_xlabel(f'BJD-{time_zero_point:g}')

        _ax2_first.set_xlabel(f'BJD-{time_zero_point:g}')

    else:

        ax2.set_xlabel('BJD')

        _ax2_first.set_xlabel('BJD')

    title = build_figure_title(data_type_in_mask, visit_to_plot,

                                instrument_to_plot, plModel)

    _title_ax = ax_raw_panel if ax_raw_panel is not None else ax0

    _title_ax.set_title(title)

    # Set title on the figure instead of a specific panel

    _any_ax = ax0[0] if isinstance(ax0, list) else ax0

    _any_ax.get_figure().suptitle(title, fontsize=12)

    # Hide redundant y-axis labels and ticks for internal subpanels

    # and ensure y-axis labels are only on the leftmost column.

    if isinstance(ax0, list) and len(ax0) > 1:

        for i in range(1, len(ax0)):

            # Hide redundant y-axis labels and ticks for internal subpanels

            for ax_group in [ax0, ax1, ax2]:

                ax = ax_group[i]

                plt.setp(ax.get_yticklabels(), visible=False)

                plt.setp(ax.get_yaxis().get_label(), visible=False)

            if ax_raw_panel is not None:

                ax_raw = ax_raw_panel[i] if isinstance(ax_raw_panel, list) else None

                if ax_raw:

                    plt.setp(ax_raw.get_yticklabels(), visible=False)

                    plt.setp(ax_raw.get_yaxis().get_label(), visible=False)

        # Hide x-ticklabels for all rows except the bottom-most (ax2)

        # and for the raw panel if it exists.

        for ax_group in [ax0, ax1]:

            for ax in ax_group:

                plt.setp(ax.get_xticklabels(), visible=False)

        if ax_raw_panel is not None:

            _raw_list = ax_raw_panel if isinstance(ax_raw_panel, list) else [ax_raw_panel]

            for ax in _raw_list:

                plt.setp(ax.get_xticklabels(), visible=False)

    import warnings

    with warnings.catch_warnings():

        warnings.filterwarnings('ignore', message='.*not compatible with tight_layout.*')

        warnings.filterwarnings('ignore', message='.*tight_layout.*not compatible.*')

        ax0.get_figure().tight_layout()

        # Use any axis to get the figure

        _any_ax.get_figure().tight_layout()

        # Adjust top to make room for suptitle

        _any_ax.get_figure().tight_layout(rect=[0, 0, 1, 0.96])

            num_chains = saved_data.get('num_chains', None)

            if num_chains is not None:

                if saved_sampler_type in ('warmup', 'warmup_complete'):

                    log_func(f"Checkpoint type: warmup ({existing_warmup_done} iters, "

                    log_func(logFile, f"Checkpoint type: warmup ({existing_warmup_done} iters, "

                             f"{format_time(existing_warmup_elapsed)}), "

                             f"num_chains={num_chains}")

                else:

                    log_func(f"Checkpoint verified: num_chains={num_chains} (from checkpoint)")

                    log_func(logFile, f"Checkpoint verified: num_chains={num_chains} (from checkpoint)")

            if existing_samples is not None:

                scalar_flat = min(len(v) for v in existing_samples.values())

                    first_walker_param = next(iter(existing_mcmc_state.z.values()), None)

                    if first_walker_param is not None:

                        actual_chain_count = len(first_walker_param)

                        log_func(f"Inferred {actual_chain_count} chains from MCMC state z shape")

                        log_func(logFile, f"Inferred {actual_chain_count} chains from MCMC state z shape")

                if actual_chain_count is None:

                    candidate_counts = [4, 16, 32, 34]

                    for candidate in candidate_counts:

                        if scalar_flat % candidate == 0:

                            actual_chain_count = candidate

                            log_func(f"Inferred {actual_chain_count} chains from scalar param divisibility")

                            log_func(logFile, f"Inferred {actual_chain_count} chains from scalar param divisibility")

                            break

                if actual_chain_count is None:

                    actual_chain_count = num_chains

                    log_func(f"Using saved num_chains={actual_chain_count} (may be incorrect if production used different sampler)")

                    log_func(logFile, f"Using saved num_chains={actual_chain_count} (may be incorrect if production used different sampler)")

                num_chains = actual_chain_count

                existing_num_samples = scalar_flat // num_chains

                log_func(f"Loaded checkpoint: {len(plData.time)} data points, "

                log_func(logFile, f"Loaded checkpoint: {len(plData.time)} data points, "

                         f"{existing_num_samples} samples/chain × {num_chains} chains "

                         f"= {scalar_flat} scalar-param samples ({saved_sampler_type})")

                k: v for k, v in mcmc_config_peek.items()

                if k == 'use_float64'

            })

            use_float64 = _configure_jax_precision(param_config, mcmc_config_restored, log_func=log_func)

            use_float64 = _configure_jax_precision(param_config, mcmc_config_restored, log_func=lambda *args: log_func(logFile, args[-1]))            

            if use_float64:

                for attr in ('time', 'y', 'dy', 'exptime'):

                    arr = getattr(plData, attr, None)

    # PATH A: Extending existing chains (reset=False, pickle found)

    # ==================================================================

    if existing_samples is not None:

        # Ensure we have the correct number of samples per chain from the checkpoint

        # (existing_num_samples is returned by load_checkpoint)

        samples = existing_samples.copy()

        # Reshape each parameter from flat storage to (num_chains, n_samples[, n_dim]).

        for name, vals in samples.items():

            n_total = len(vals)

            n_scalar = num_chains * existing_num_samples

            if n_total == n_scalar:

            if n_scalar > 0 and n_total == n_scalar:

                samples_grouped[name] = vals.reshape(num_chains, existing_num_samples)

            elif n_total > n_scalar and n_total % n_scalar == 0:

            elif n_scalar > 0 and n_total > n_scalar and n_total % n_scalar == 0:

                n_extra = n_total // n_scalar

                samples_grouped[name] = vals.reshape(num_chains, existing_num_samples, n_extra)

            else: