Commit a81d624c authored by Gaetano Scandariato's avatar Gaetano Scandariato
Browse files

prepared plan for better plotting of nontime trends

parent 23f47885
Loading
Loading
Loading
Loading
Loading
+27 −0
Changes for .kilo/panel_alignment_explanation.md: 27 added lines, 0 removed lines.
Original line number Diff line number Diff line
# Alignment between Panel 0 and Panel 1

## Affected Modules
This logic is primarily implemented and managed in:
- `PyAETNA/reporting/plot_data_prep.py`: Handles the pre-computation of detrended data (`overlay_y`) and the preparation of the `PlotContext`.
- `PyAETNA/reporting/panel_raw.py`: Renders Panel 0, showing the raw data and the non-time correction line.
- `PyAETNA/reporting/panel_overlay.py`: Renders Panel 1, showing the detrended data and model components.

## Overview
There is a conceptual and mathematical distinction between the data presented in Panel 0 (Raw Data) and Panel 1 (Model Overlay) regarding how non-time corrections (e.g., polynomial or GP trends based on non-BJD variables) are handled.

## Panel 0: Raw Data and Corrections
Panel 0 is designed to show the unmodified observations alongside the corrections being applied. 
- **The Correction Line**: The line plotted in Panel 0 represents the total correction including the normalization baseline. 
    - For **Radial Velocity (RV)**: $\text{Line} = \text{norm} + \text{correction}$
    - For **Light Curves (LC)**: $\text{Line} = \text{norm} \times (1 + \text{correction})$
- **Purpose**: This allows the user to visually verify that the correction line follows the raw data trends and that the normalization baseline is correctly positioned.

## Panel 1: Model Overlay and Detrending
Panel 1 is aimed at showing the data, models, and confidence bands **after** non-time trends have been removed.
- **Detrended Data**: The data points in Panel 1 are computed by subtracting the non-time correction from the raw data.
    - For **RV**: $y_{\text{overlay}} = y_{\text{raw}} - (\text{norm} + \text{correction})$
    - For **LC**: $y_{\text{overlay}} = y_{\text{raw}} - (\text{norm} \times (1 + \text{correction}))$
- **Confidence Bands**: It is important to note that the confidence bands for the models shown in Panel 1 **still take into account the uncertainties from the non-time trends**. Even though the trends themselves have been removed from the plotted positions of the data and models, their contribution to the overall uncertainty budget is preserved in the bands. This need to be fixed.

## Summary of Misalignment
The corrected data in Panel 1 seems good, but the correction represented in Panel 0 shows a systematic offset wrt the data. The goal is to make Panel 0 shows the nontime correction applied to the data that are then shown in Panel 1. Panel 1, moreover, should show all the remaining components and shouldn't contain any reference to nontime trends.

PyAETNA/.kilo/step1.md

0 → 100644
+18 −0
Changes for PyAETNA/.kilo/step1.md: 18 added lines, 0 removed lines.
Original line number Diff line number Diff line
# Step 1: Enhance PlotContext and Data Preparation

## Context
`PlotContext` currently aggregates non-time trends into a single entity. To enable detailed partial residual analysis, we need to track individual non-time components separately.

## Purpose
Modify the data preparation pipeline to:
1. Identify each individual non-time component.
2. Retrieve the corresponding observable for each component from `plData.dBasis`.
3. Compute the Maximum A Posteriori (MAP) estimate and the 16-84% confidence bands for each component.

## Files to Modify
- `PyAETNA/reporting/plot_data_prep.py`

## Verification & Tests
- **Component Tracking**: Verify that `PlotContext` now contains a list or dictionary of individual non-time components rather than a single aggregate.
- **Data Integrity**: Check that the array lengths of the computed MAP and confidence bands match the length of the observables.
- **Summation Check**: Verify that the sum of all individual non-time components equals the total non-time trend previously calculated.

PyAETNA/.kilo/step2.md

0 → 100644
+17 −0
Changes for PyAETNA/.kilo/step2.md: 17 added lines, 0 removed lines.
Original line number Diff line number Diff line
# Step 2: Create the Partial Residuals Renderer

## Context
The necessary component-level data is now available in `PlotContext`. We now need a visualization layer to render this data into the final report.

## Purpose
Create a new rendering module `panel_nontime_residuals.py` to generate a PDF page. For each identified non-time component:
1. **Compute Partial Residuals**: Calculate $y_{raw} - (\text{Total Model} - \text{Component})$. This isolates the effect of the specific component.
2. **Visualization**: Plot these partial residuals against the corresponding observable.
3. **Overlay**: Plot the MAP line and the 16-84% confidence bands over the scatter plot to visualize the fit quality for that specific component.

## Files to Create
- `PyAETNA/reporting/panel_nontime_residuals.py` (New File)

## Verification & Tests
- **Visual Inspection**: Generate a test PDF and verify the layout of the partial residuals page.
- **Trend Validation**: Verify that the MAP line correctly follows the trend of the partial residuals for each component.

PyAETNA/.kilo/step3.md

0 → 100644
+16 −0
Changes for PyAETNA/.kilo/step3.md: 16 added lines, 0 removed lines.
Original line number Diff line number Diff line
# Step 3: Integrate into PDF Report

## Context
The data preparation and the rendering module are complete, but the new page is not yet integrated into the main report generation workflow.

## Purpose
Modify the main PDF report logic to trigger the partial residuals renderer.
1. Detect if non-time trends are present in the current fit results.
2. If present, call the `panel_nontime_residuals.py` renderer to insert the partial residuals page into the final PDF document.

## Files to Modify
- `PyAETNA/reporting/pdf_report.py`

## Verification & Tests
- **End-to-End Test**: Run a full fit using a configuration that includes non-time components.
- **PDF Validation**: Open the resulting PDF and verify that the "Partial Residuals" page is present and correctly populated with data.
+84 −77
Changes for PyAETNA/reporting/panel_overlay.py: 84 added lines, 77 removed lines.
Original line number Diff line number Diff line
@@ -39,6 +39,14 @@ def _render_rv_overlay(overlay_axes: list, ctx: PlotContext, ax_legend=None) ->
    # otherwise use the full fine-grid trend.
    tr_mf_ov = ctx.tr_time_mf if ctx.has_non_time_trends else ctx.tr_mf

    # Coarse-grid trend for band centres: time-only when nontime trends present.
    # ctx.ov_tr_time_grid is the posterior median of (tr - tr_nontime) on the grid.
    tr_grid_ov = (
        ctx.ov_tr_time_grid
        if (ctx.has_non_time_trends and ctx.ov_tr_time_grid is not None)
        else ctx.tr_grid
    )

    for j, ax in enumerate(overlay_axes):
        # 1. Baseline (Normalization γ)
        ax.plot(ctx.t_mf, ctx.norm_mf, linestyle='--', color='C0',
@@ -49,76 +57,72 @@ def _render_rv_overlay(overlay_axes: list, ctx: PlotContext, ax_legend=None) ->
                            color='C0', alpha=0.12, zorder=1.5,
                            label='_nolegend_')

        # 2. Baseline + Trend
        # 2. Baseline + Trend (time-only)
        if has_trend:
            tlbl = ctx.trend_label
            ax.plot(ctx.t_mf, ctx.norm_mf + tr_mf_ov, '-', color='C4',
                    linewidth=1.2, zorder=3, alpha=0.8,
                    label=f'{tlbl} + Norm' if j == 0 else '_nolegend_')
            if ctx.tr_time_centred_q16_grid is not None:
                # The original bands were [norm_q16, norm_q84] and [tr_q16, tr_q84]
                # We want to combine the widths of these bands in quadrature.
                # Width of norm band: (norm_q84 - norm_q16) / 2
                # Width of trend band: (tr_q84 - tr_q16) / 2
                norm_half_width = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_half_width = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0
                combined_half_width = np.sqrt(norm_half_width**2 + tr_half_width**2)
                
                # tr_time_centred_q16/q84 are already percentiles of
                # (norm + time_trend), nontime-free in both branches of
                # confidence_bands.py — use them directly as band bounds.
                ax.fill_between(ctx.t_grid,
                                (ctx.norm_grid + ctx.tr_grid) - combined_half_width,
                                (ctx.norm_grid + ctx.tr_grid) + combined_half_width,
                                ctx.tr_time_centred_q16_grid,
                                ctx.tr_time_centred_q84_grid,
                                color='C4', alpha=0.12, zorder=2.5,
                                label='_nolegend_')

        # 3. Baseline + Trend + Planet
        if ctx.has_planets:
            # We add the planet on top of the trend (if any) and norm
            planet_val = ctx.norm_mf + (tr_mf_ov if has_trend else 0) + ctx.pm_mf
            planet_label = 'Planet + Trend + Norm' if has_trend else 'Planet + Norm'
            ax.plot(ctx.t_mf, planet_val, '-', color='C1',
                    linewidth=1.5, zorder=4,
                    label=planet_label if j == 0 else '_nolegend_')
            if ctx.pn_q16_grid is not None:
                # Combine widths of previous bands in quadrature
                norm_half_width = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_half_width = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_half_width = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0
                combined_half_width = np.sqrt(norm_half_width**2 + tr_half_width**2 + pn_half_width**2)
                
            if ctx.has_non_time_trends and ctx.ov_pt_q16_grid is not None:
                # Pre-computed nontime-free overlay band — use directly.
                ax.fill_between(ctx.t_grid,
                                ctx.ov_pt_q16_grid, ctx.ov_pt_q84_grid,
                                color='C1', alpha=0.15, zorder=3.5,
                                label='_nolegend_')
            elif ctx.pn_q16_grid is not None:
                # No nontime trends: combine widths in quadrature.
                norm_hw = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_hw   = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_hw   = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0
                hw = np.sqrt(norm_hw**2 + tr_hw**2 + pn_hw**2)
                centre = ctx.norm_grid + (tr_grid_ov if has_trend else 0) + ctx.pm_grid
                ax.fill_between(ctx.t_grid,
                                (ctx.norm_grid + (ctx.tr_grid if has_trend else 0) + ctx.pm_grid) - combined_half_width,
                                (ctx.norm_grid + (ctx.tr_grid if has_trend else 0) + ctx.pm_grid) + combined_half_width,
                                centre - hw, centre + hw,
                                color='C1', alpha=0.15, zorder=3.5,
                                label='_nolegend_')

        # 4. Baseline + Trend + Planet + GP
        if has_gp:
            # We add the GP on top of everything else
            # Note: ctx.gp_grid is on coarse grid, so we use t_grid
            gp_val = ctx.norm_grid + (ctx.tr_grid if has_trend else 0) + \
            gp_val = ctx.norm_grid + (tr_grid_ov if has_trend else 0) + \
                     ctx.pm_grid + ctx.gp_grid
            
            # Build label based on what is actually present
            gp_label_parts = ['GP']
            if ctx.has_planets: gp_label_parts.append('Planet')
            if has_trend: gp_label_parts.append('Trend')
            gp_label_parts.append('Norm')
            gp_label = ' + '.join(gp_label_parts)

            ax.plot(ctx.t_grid, gp_val, '-', color='C3',
                    linewidth=1.2, zorder=5, alpha=0.8,
                    label=gp_label if j == 0 else '_nolegend_')
            if ctx.gp_centred_q16_grid is not None:
                # Combine widths of all previous bands in quadrature
                norm_half_width = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_half_width = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_half_width = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0 if ctx.has_planets and ctx.pn_q16_grid is not None else 0
                gp_half_width = (ctx.gp_centred_q84_grid - ctx.gp_centred_q16_grid) / 2.0
                combined_half_width = np.sqrt(norm_half_width**2 + tr_half_width**2 + pn_half_width**2 + gp_half_width**2)
                
                    label=' + '.join(gp_label_parts) if j == 0 else '_nolegend_')
            if ctx.has_non_time_trends and ctx.ov_tot_q16_grid is not None:
                # Pre-computed nontime-free total band — use directly.
                ax.fill_between(ctx.t_grid,
                                ctx.ov_tot_q16_grid, ctx.ov_tot_q84_grid,
                                color='C3', alpha=0.12, zorder=4.5,
                                label='_nolegend_')
            elif ctx.gp_centred_q16_grid is not None:
                norm_hw = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_hw   = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_hw   = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0 if ctx.has_planets and ctx.pn_q16_grid is not None else 0
                gp_hw   = (ctx.gp_centred_q84_grid - ctx.gp_centred_q16_grid) / 2.0
                hw = np.sqrt(norm_hw**2 + tr_hw**2 + pn_hw**2 + gp_hw**2)
                ax.fill_between(ctx.t_grid,
                                gp_val - combined_half_width,
                                gp_val + combined_half_width,
                                gp_val - hw, gp_val + hw,
                                color='C3', alpha=0.12, zorder=4.5,
                                label='_nolegend_')

@@ -132,6 +136,13 @@ def _render_lc_overlay(overlay_axes: list, ctx: PlotContext, ax_legend=None) ->
    has_gp    = ctx.show_gp and np.any(ctx.gp != 0)
    tr_mf_ov  = ctx.tr_time_mf if ctx.has_non_time_trends else ctx.tr_mf

    # Coarse-grid trend for band centres: time-only when nontime trends present.
    tr_grid_ov = (
        ctx.ov_tr_time_grid
        if (ctx.has_non_time_trends and ctx.ov_tr_time_grid is not None)
        else ctx.tr_grid
    )

    z_grid = np.zeros_like(ctx.pm_grid)
    z_mf   = np.zeros_like(ctx.pm_mf)
    mode   = ctx.normalization_mode
@@ -156,15 +167,17 @@ def _render_lc_overlay(overlay_axes: list, ctx: PlotContext, ax_legend=None) ->
            ax.plot(ctx.t_mf, trend_baseline, '-', color='C4',
                    linewidth=1.2, zorder=3, alpha=0.8,
                    label=f'{ctx.trend_label} + Baseline' if j == 0 else '_nolegend_')
            if ctx.tr_time_centred_q16_grid is not None:
                # Combine widths of Norm and Trend bands in quadrature
            if ctx.has_non_time_trends and ctx.ov_tr_q16_grid is not None:
                ax.fill_between(ctx.t_grid,
                                ctx.ov_tr_q16_grid,
                                ctx.ov_tr_q84_grid,
                                color='C4', alpha=0.12, zorder=2.5,
                                label='_nolegend_')
            elif ctx.tr_time_centred_q16_grid is not None:
                norm_half_width = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_half_width = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0
                combined_half_width = np.sqrt(norm_half_width**2 + tr_half_width**2)
                
                # For LC, we apply the combined uncertainty to the baseline + trend model
                # Note: build_partial_for_plot is used for the center line
                baseline_trend_grid = build_partial_for_plot(ctx.norm_grid, z_grid, ctx.tr_grid, mode, 'LC')
                baseline_trend_grid = build_partial_for_plot(ctx.norm_grid, z_grid, tr_grid_ov, mode, 'LC')
                ax.fill_between(ctx.t_grid,
                                baseline_trend_grid - combined_half_width,
                                baseline_trend_grid + combined_half_width,
@@ -173,7 +186,6 @@ def _render_lc_overlay(overlay_axes: list, ctx: PlotContext, ax_legend=None) ->

        # 3. Baseline + Trend + Planet
        if ctx.has_planets:
            # We add the planet on top of the trend (if any) and baseline
            trend_comp = tr_mf_ov if has_trend else z_mf
            planet_baseline = build_partial_for_plot(
                ctx.norm_mf, ctx.pm_mf, trend_comp, mode, 'LC')
@@ -182,32 +194,29 @@ def _render_lc_overlay(overlay_axes: list, ctx: PlotContext, ax_legend=None) ->
            ax.plot(ctx.t_mf, planet_baseline, '-', color='C1',
                    linewidth=1.5, zorder=4,
                    label=planet_label if j == 0 else '_nolegend_')
            if ctx.pn_q16_grid is not None:
                # Combine widths of Norm, Trend, and Planet bands in quadrature
                norm_half_width = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_half_width = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_half_width = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0
                combined_half_width = np.sqrt(norm_half_width**2 + tr_half_width**2 + pn_half_width**2)
                
                # Center line for the band
                planet_grid = build_partial_for_plot(ctx.norm_grid, ctx.pm_grid, (ctx.tr_grid if has_trend else z_grid), mode, 'LC')
            if ctx.has_non_time_trends and ctx.ov_pt_q16_grid is not None:
                ax.fill_between(ctx.t_grid,
                                ctx.ov_pt_q16_grid, ctx.ov_pt_q84_grid,
                                color='C1', alpha=0.15, zorder=3.5,
                                label='_nolegend_')
            elif ctx.pn_q16_grid is not None:
                norm_hw = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_hw   = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_hw   = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0
                hw = np.sqrt(norm_hw**2 + tr_hw**2 + pn_hw**2)
                planet_grid = build_partial_for_plot(ctx.norm_grid, ctx.pm_grid, (tr_grid_ov if has_trend else z_grid), mode, 'LC')
                ax.fill_between(ctx.t_grid, 
                                planet_grid - combined_half_width, 
                                planet_grid + combined_half_width,
                                planet_grid - hw, 
                                planet_grid + hw,
                                color='C1', alpha=0.15, zorder=3.5,
                                label='_nolegend_')

        # 4. Baseline + Trend + Planet + GP
        if has_gp:
            # We add the GP on top of everything else
            # Note: GP is on coarse grid, so we use t_grid and coarse components
            # We use build_partial_for_plot to handle normalization mode (additive/multiplicative)
            # For the GP, we need the combined planet + trend components
            combined_pt = ctx.pm_grid + (ctx.tr_grid if has_trend else 0)
            combined_pt = ctx.pm_grid + (tr_grid_ov if has_trend else 0)
            gp_baseline = build_partial_for_plot(
                ctx.norm_grid, combined_pt, ctx.gp_grid, mode, 'LC')
            
            # Build label based on what is actually present
            gp_label_parts = ['GP']
            if ctx.has_planets: gp_label_parts.append('Planet')
            if has_trend: gp_label_parts.append('Trend')
@@ -217,24 +226,22 @@ def _render_lc_overlay(overlay_axes: list, ctx: PlotContext, ax_legend=None) ->
            ax.plot(ctx.t_grid, gp_baseline, '-', color='C3',
                    linewidth=1.2, zorder=5, alpha=0.8,
                    label=gp_label if j == 0 else '_nolegend_')
            if ctx.gp_centred_q16_grid is not None:
                # Combine widths of all components including GP in quadrature
                norm_half_width = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_half_width = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_half_width = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0 if ctx.has_planets and ctx.pn_q16_grid is not None else 0
                gp_half_width = (ctx.gp_centred_q84_grid - ctx.gp_centred_q16_grid) / 2.0
                combined_half_width = np.sqrt(norm_half_width**2 + tr_half_width**2 + pn_half_width**2 + gp_half_width**2)
                
            if ctx.has_non_time_trends and ctx.ov_tot_q16_grid is not None:
                ax.fill_between(ctx.t_grid,
                                gp_baseline - combined_half_width,
                                gp_baseline + combined_half_width,
                                ctx.ov_tot_q16_grid, ctx.ov_tot_q84_grid,
                                color='C3', alpha=0.12, zorder=4.5,
                                label='_nolegend_')
            elif ctx.gp_centred_q16_grid is not None:
                norm_hw = (ctx.norm_q84_grid - ctx.norm_q16_grid) / 2.0 if ctx.norm_q16_grid is not None else 0
                tr_hw   = (ctx.tr_time_centred_q84_grid - ctx.tr_time_centred_q16_grid) / 2.0 if has_trend and ctx.tr_time_centred_q16_grid is not None else 0
                pn_hw   = (ctx.pn_q84_grid - ctx.pn_q16_grid) / 2.0 if ctx.has_planets and ctx.pn_q16_grid is not None else 0
                gp_hw   = (ctx.gp_centred_q84_grid - ctx.gp_centred_q16_grid) / 2.0
                hw = np.sqrt(norm_hw**2 + tr_hw**2 + pn_hw**2 + gp_hw**2)
                ax.fill_between(ctx.t_grid,
                                gp_baseline - hw, gp_baseline + hw,
                                color='C3', alpha=0.12, zorder=4.5,
                                label='_nolegend_')

    # Legend is now handled in render_panel_overlay to include all elements
    _set_cluster_xlims(overlay_axes, ctx)

    # Legend is now handled in render_panel_overlay to include all elements
    _set_cluster_xlims(overlay_axes, ctx)


Loading