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

Add technical analysis document on RV-BIS independence in Bayesian joint modeling

- Introduced a comprehensive markdown document detailing the analysis of information redundancy between Radial Velocity (RV) and Bisector Inverse Slope (BIS) in the context of Bayesian joint modeling, referencing Rajpaul et al. (2015).
- Included mathematical formulations, empirical validations, and diagnostic criteria for assessing the independence of RV and BIS measurements.
parent efd9b8b0
Loading
Loading
Loading
Loading
Loading
+7 −1
Changes for PyAETNA/reporting/periodogram_plots.py: 7 added lines, 1 removed line.
Original line number Original line Diff line number Diff line
@@ -707,6 +707,12 @@ def plot_residuals_vs_dbasis(plData, residuals, plModel=None, visit_to_plot=None
    else:
    else:
        visit_filter = np.ones(len(residuals), dtype=bool)
        visit_filter = np.ones(len(residuals), dtype=bool)


    # Determine y-axis label based on instrument data_type
    # Use the data_type from plModel for the specific instrument
    # plModel.data_types is a dict mapping inst_id -> 'RV' or 'LC'
    dtype = getattr(plModel, 'data_types', {}).get(instrument_to_plot, None)
    y_label='Residual '+dtype

    visit_mask = inst_filter & visit_filter
    visit_mask = inst_filter & visit_filter
    residuals = residuals[visit_mask]
    residuals = residuals[visit_mask]
    
    
@@ -782,7 +788,7 @@ def plot_residuals_vs_dbasis(plData, residuals, plModel=None, visit_to_plot=None
                    pass
                    pass
        
        
        ax.set_xlabel(var_name)
        ax.set_xlabel(var_name)
        ax.set_ylabel('Residual Flux')
        ax.set_ylabel(y_label)
        ax.axhline(0, color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
        ax.axhline(0, color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
        ax.grid(True, alpha=0.3)
        ax.grid(True, alpha=0.3)
        
        
+1 −0
Changes for README.md: 1 added line, 0 removed lines.
Original line number Original line Diff line number Diff line
@@ -259,6 +259,7 @@ This project is licensed under the MIT License — see the [LICENSE](LICENSE) fi
- check the checkpoint_interval logic
- check the checkpoint_interval logic
- implement different limb darkening laws
- implement different limb darkening laws
- better arrangement of legends and text in plots
- better arrangement of legends and text in plots
- remove double INFO messages in the log


## Contact
## Contact


+16 −15
Changes for docs/configuration.md: 16 added lines, 15 removed lines.
Original line number Original line Diff line number Diff line
@@ -1067,6 +1067,13 @@ GP:
CHEOPS sees the sum $K_\text{GP\_time} + K_\text{GP\_roll}$.  TESS sees
CHEOPS sees the sum $K_\text{GP\_time} + K_\text{GP\_roll}$.  TESS sees
only $K_\text{GP\_time}$.
only $K_\text{GP\_time}$.


> **Note:** When several `GP_<name>` blocks list the **same instrument**, their
`x_variable`s are stacked into a multidimensional input and kernels are
summed for that instrument:
$K_{\text{total}} = K_1(x_{:,0}) + K_2(x_{:,1}) + \ldots$\
> When GP sections target **different instruments**, each instrument gets
only the kernel(s) from sections that include it.

#### Scenario 5 — Spectral Mixture GP
#### Scenario 5 — Spectral Mixture GP


The Spectral Mixture kernel can efficiently model complex variability and periodic signals with a flexible mixture of components (see https://tinygp.readthedocs.io/en/latest/tutorials/kernels.html). Unlike other kernels, each parameter $\sigma$, $l_s$, and $f$ in SM can be specified independently for multiple components, or broadcast to all components as needed.
The Spectral Mixture kernel can efficiently model complex variability and periodic signals with a flexible mixture of components (see https://tinygp.readthedocs.io/en/latest/tutorials/kernels.html). Unlike other kernels, each parameter $\sigma$, $l_s$, and $f$ in SM can be specified independently for multiple components, or broadcast to all components as needed.
@@ -1218,17 +1225,7 @@ GP:
For detailed benchmarks and numerical analysis, see
For detailed benchmarks and numerical analysis, see
[GP Per-Visit Analysis](gp_per_visit_analysis/GP_PER_VISIT_ANALYSIS.md).
[GP Per-Visit Analysis](gp_per_visit_analysis/GP_PER_VISIT_ANALYSIS.md).


#### Multiple GPs
### About Multi-GPs

When several `GP_<name>` blocks list the **same instrument**, their
`x_variable`s are stacked into a multidimensional input and kernels are
summed for that instrument:
$K_{\text{total}} = K_1(x_{:,0}) + K_2(x_{:,1}) + \ldots$

When GP sections target **different instruments**, each instrument gets
only the kernel(s) from sections that include it.

#### About Multi-GPs


PyAETNA’s GP system is fundamentally an **instrumentwise GP framework**.
PyAETNA’s GP system is fundamentally an **instrumentwise GP framework**.
Each GP section defines a kernel for one or more instruments, and the
Each GP section defines a kernel for one or more instruments, and the
@@ -1255,7 +1252,7 @@ factorised as independent GP blocks. In other words:
- PyAETNA does not currently do **full joint output covariance** between
- PyAETNA does not currently do **full joint output covariance** between
  RV and BIS.
  RV and BIS.


### Comparison with Rajpaul et al. (2015)
#### Comparison with Rajpaul et al. (2015)


Rajpaul et al. build a model from a single latent stellar activity GP
Rajpaul et al. build a model from a single latent stellar activity GP
$G(t)$ and its derivative $\dot{G}(t)$, then express RV and BIS as linear
$G(t)$ and its derivative $\dot{G}(t)$, then express RV and BIS as linear
@@ -1269,7 +1266,9 @@ analytic cross-covariance terms:
This is a stronger coupling than the “same period/lengthscale” fit, and it is
This is a stronger coupling than the “same period/lengthscale” fit, and it is
mathematically different from PyAETNA’s instrumentwise approach.
mathematically different from PyAETNA’s instrumentwise approach.


### Comparison with PyORBIT multiGP
> **Important note:** For a detailed technical analysis of this formalism, advantages and risks, see [RV-BIS Independence Analysis](rajpaul_bis_analysis.md).

#### Comparison with PyORBIT multiGP


PyORBIT does provide a true multidimensional TinyGP implementation. Its
PyORBIT does provide a true multidimensional TinyGP implementation. Its
multidimensional GP classes use a single shared kernel family for all
multidimensional GP classes use a single shared kernel family for all
@@ -1286,7 +1285,7 @@ This means:
  kernel forms for different outputs; the kernel family is shared across the
  kernel forms for different outputs; the kernel family is shared across the
  joint dataset.
  joint dataset.


### Pros and cons
#### Pros and cons


**PyAETNA pros**
**PyAETNA pros**
- Very flexible: different instruments may use different kernels.
- Very flexible: different instruments may use different kernels.
@@ -1312,7 +1311,7 @@ This means:
- The multiGP block must use a single kernel family for all datasets it
- The multiGP block must use a single kernel family for all datasets it
  couples.
  couples.


### Practical guidance
#### Practical guidance


- If you need speed and flexibility, and you only want to share physical
- If you need speed and flexibility, and you only want to share physical
  timescales like rotation period or decay lengthscale, PyAETNA’s existing
  timescales like rotation period or decay lengthscale, PyAETNA’s existing
@@ -1354,6 +1353,8 @@ In this example:
> `INSTRUMENT` ID/name or placing them in a separate dataset entry, rather
> `INSTRUMENT` ID/name or placing them in a separate dataset entry, rather
> than keeping them as an ancillary column attached to the RV rows.
> than keeping them as an ancillary column attached to the RV rows.


> **Important note:** When analyzing multiple indicators extracted from the same CCF profile (e.g., RV and BIS), be aware of the risk of information redundancy. For a detailed technical analysis of when this is a concern and how it relates to the model structure, see [RV-BIS Independence Analysis](rajpaul_bis_analysis.md).

---
---


## Jitter
## Jitter
+136 −0
Changes for docs/convert_md_to_pdf.py: 136 added lines, 0 removed lines.
Original line number Original line Diff line number Diff line
"""
Convert *.md files in the current directory to PDF.

Pipeline:
  1. pandoc  : Markdown -> self-contained HTML (with embedded CSS for readability)
  2. google-chrome --headless : HTML -> PDF

Requirements (must be in PATH):
  - pandoc  (https://pandoc.org)
  - google-chrome  (any Chrome/Chromium binary)
"""

import os
import subprocess
import sys
import tempfile
from pathlib import Path

# Get the directory where this script is located
SCRIPT_DIR = Path(__file__).parent

CHROME_CANDIDATES = ["google-chrome", "chromium-browser", "chromium", "google-chrome-stable"]

CSS = """
body {
    font-family: Georgia, 'Times New Roman', serif;
    font-size: 11pt;
    line-height: 1.6;
    max-width: 900px;
    margin: 40px auto;
    padding: 0 20px;
    color: #111;
}
h1 { font-size: 1.8em; border-bottom: 2px solid #333; padding-bottom: 6px; }
h2 { font-size: 1.3em; border-bottom: 1px solid #aaa; padding-bottom: 3px; margin-top: 1.6em; }
h3 { font-size: 1.1em; margin-top: 1.2em; }
img { max-width: 100%; height: auto; display: block; margin: 1em auto; }
table { border-collapse: collapse; width: 100%; margin: 1em 0; font-size: 10pt; }
th, td { border: 1px solid #bbb; padding: 6px 10px; text-align: left; }
th { background: #eee; font-weight: bold; }
tr:nth-child(even) td { background: #f9f9f9; }
code { background: #f0f0f0; padding: 1px 4px; border-radius: 3px; font-size: 0.92em; }
pre { background: #f0f0f0; padding: 12px; border-radius: 4px; overflow-x: auto; }
blockquote { border-left: 4px solid #aaa; margin-left: 0; padding-left: 16px; color: #555; }
a { color: #1a5a99; }
"""


def find_chrome():
    for name in CHROME_CANDIDATES:
        result = subprocess.run(["which", name], capture_output=True, text=True)
        if result.returncode == 0:
            return result.stdout.strip()
    sys.exit(
        "ERROR: No Chrome/Chromium binary found. "
        "Tried: " + ", ".join(CHROME_CANDIDATES)
    )


def md_to_html(md_path, html_path, css_path):
    cmd = [
        "pandoc",
        str(md_path),
        "--from=markdown+tex_math_single_backslash",
        "--to=html5",
        "--standalone",
        "--css=" + str(css_path),
        "--metadata", "title=" + md_path.stem,
        "--output", str(html_path),
        "--mathml",  # Use MathML for native browser rendering (better for headless PDF)
        "--embed-resources",  # Embed images and CSS directly into the HTML
        "--resource-path=" + str(md_path.parent)
    ]
    subprocess.run(cmd, check=True)


def html_to_pdf(chrome, html_path, pdf_path):
    cmd = [
        chrome,
        "--headless",
        "--disable-gpu",
        "--no-sandbox",
        "--disable-dev-shm-usage",
        "--allow-file-access-from-files",
        "--disable-web-security",
        "--run-all-compositor-stages-before-draw",
        "--virtual-time-budget=10000",
        "--print-to-pdf=" + str(pdf_path),
        "--print-to-pdf-no-header",
        str(html_path),
    ]
    subprocess.run(cmd, check=True, capture_output=True)


def main():
    chrome = find_chrome()
    print("Using Chrome:", chrome)

    # Find all markdown files in the current directory
    md_files = sorted(
        p for p in SCRIPT_DIR.iterdir()
        if p.suffix == ".md"
    )

    if not md_files:
        print("No .md files found.")
        return

    with tempfile.TemporaryDirectory() as tmpdir:
        css_path = Path(tmpdir) / "style.css"
        css_path.write_text(CSS)

        for md_path in md_files:

            if not 'rajpaul' in md_path.name:
                continue

            pdf_path = md_path.with_suffix(".pdf")
            html_path = Path(tmpdir) / (md_path.stem + ".html")

            print("  " + md_path.name + "  ->  " + pdf_path.name + " ... ", end="", flush=True)
            try:
                md_to_html(md_path, html_path, css_path)
                html_to_pdf(chrome, html_path, pdf_path)
                print("OK")
            except subprocess.CalledProcessError as exc:
                print("FAILED (return code " + str(exc.returncode) + ")")
                print("    Command: " + " ".join(str(a) for a in exc.cmd))
                if exc.stderr:
                    print("    stderr: " + exc.stderr[:400])

    print("\nDone.")


if __name__ == "__main__":
    main()
+524 −0

File added.

Preview size limit exceeded, changes collapsed.