scripts.step_9_generate_report¶

def add_heading(story, text, styles, lvl=1):
    """
    Add a heading of level `lvl` to the story.
    1 = main section, 2 = subsection, etc.
    """
    size = 18 if lvl == 1 else 14
    style = ParagraphStyle(
        name=f"Heading{lvl}",
        parent=styles["Heading1"],
        fontSize=size,
        leading=size + 2,
        spaceAfter=6,
    )
    story.append(Paragraph(text, style))
    story.append(Spacer(1, 0.3 * cm))

def add_paragraph(story, text, styles):
    """
    Add a paragraph (or multiple paragraphs) to the story.
    `text` can be a string or a list/tuple of strings.
    """
    if isinstance(text, (list, tuple)):
        text = " ".join(str(t) for t in text if t is not None)
    else:
        text = str(text)

    for para in text.split("\n\n"):
        para = para.strip()
        if not para:
            continue
        story.append(Paragraph(para, styles["BodyText"]))
        story.append(Spacer(1, 0.2 * cm))

def add_plot_if_available(story, relative_plot: str, caption: str, styles, width_cm: float = 14.0):
    """
    Embed a plot into the report if present in report/plots.
    Prefer PNG (created via ImageMagick). If missing, try embedding the PDF directly.
    """
    base = REPORT_PLOTS / relative_plot
    png_path = base.with_suffix(".png")
    pdf_path = base.with_suffix(".pdf")

    img_path = None
    if png_path.exists():
        img_path = png_path
    elif pdf_path.exists():
        # reportlab can embed PDFs as images only in limited cases;
        # but often it works when Ghostscript is installed.
        img_path = pdf_path
    else:
        return

    img = Image(str(img_path), width=width_cm * cm, height=(width_cm * 0.65) * cm)
    story.append(img)
    story.append(Spacer(1, 0.15 * cm))
    story.append(Paragraph(caption, styles["Italic"]))
    story.append(Spacer(1, 0.4 * cm))

def build_knobs_paragraph() -> str:
    """
    Describe gating + model knobs by importing step scripts.
    """
    sys.path.insert(0, str(PROJECT_ROOT / "scripts"))
    try:
        import step_1a_fit_upregulation as s1a  # type: ignore
    except Exception:
        return textwrap.dedent(
            """
            All analyses share a consistent gating strategy and modeling
            backbone. Forward- and side-scatter channels define a coarse
            rectangle gate to remove debris, followed by an FSC-H/FSC-A
            singlet gate to exclude doublets. TagBFP (BV421-A) serves as the
            repressor proxy, and mCherry (PE-A) serves as the reporter.

            Kinetic parameters are obtained from exponential rise/decay fits
            to normalized BFP trajectories, while steady-state dose–response
            curves are captured by Hill functions. These parameters then feed
            into a two-variable ODE system (repressor R, reporter Y) used to
            estimate delays and to explore a broader design space.
            """
        ).strip()

    # If import worked, use the actual constants
    try:
        fsc_min = float(s1a.BOUND_MIN[s1a.CH_FSC_A])
        fsc_max = float(s1a.BOUND_MAX[s1a.CH_FSC_A])
        ssc_min = float(s1a.BOUND_MIN[s1a.CH_SSC_A])
        ssc_max = float(s1a.BOUND_MAX[s1a.CH_SSC_A])
        ratio_lo = float(s1a.SINGLET_RATIO_LOW)
        ratio_hi = float(s1a.SINGLET_RATIO_HIGH)
    except Exception:
        fsc_min = 0.4e5
        fsc_max = 2.0e5
        ssc_min = 0.2e5
        ssc_max = 1.3e5
        ratio_lo = 0.85
        ratio_hi = 1.15

    return textwrap.dedent(
        f"""
        Across all steps, the same analysis "knobs" are used to keep the
        quantitative picture consistent:

        • Boundary gate:
            FSC-A in [{fsc_min:.1f}, {fsc_max:.1f}],
            SSC-A in [{ssc_min:.1f}, {ssc_max:.1f}].
          This removes debris and very large aggregates.

        • Singlet gate:
            FSC-H/FSC-A constrained to [{ratio_lo:.2f}, {ratio_hi:.2f}],
            enriching for single cells.

        • Channels:
            tagBFP = BV421-A (repressor proxy),
            mCherry = PE-A (reporter output).

        • Kinetic models:
            – Up-regulation and down-regulation half-times are obtained from
              single-exponential fits to min–max normalized BFP trajectories.
            – At steady state (day 4), dTAG-13 dose–response is summarized
              by a Hill function (K, n) relating normalized BFP to reporter
              fold-change.

        • ODE backbone:
            A two-variable ODE system describes repressor R(t) and reporter
            Y(t), with synthesis/decay terms parameterized by the fitted
            half-times, Hill parameters and reporter decay constant alpha.
            Explicit onset delays align the simulations with the observed
            time courses.

        These knobs define a compact but biologically interpretable model
        of degron–Cas tuners as analog input–output devices.
        """
    ).strip()

def load_csv(name: str) -> Optional[pd.DataFrame]:
    """
    Load a CSV from PARAM_PATH if it exists; warn and return None otherwise.

    Args:
        name (str): Filename of the CSV to load.

    Returns:
        Optional[pd.DataFrame]: Loaded DataFrame or None if missing.
    """
    path = PARAM_PATH / name
    if not path.exists():
        print(f"[warn] Missing {name} – skipping its section.")
        return None
    df = pd.read_csv(path)
    if df.empty:
        print(f"[warn] {name} is empty.")
    return df

def mirror_plots() -> None:
    """
    Mirror all plots from `plots/` into `report/plots/` and create PNG mirrors
    for each PDF.
    """
    print("[info] Mirroring plots into report/plots ...")
    for pdf in PLOTS_PATH.rglob("*.pdf"):
        rel = pdf.relative_to(PLOTS_PATH)
        target_pdf = REPORT_PLOTS / rel
        target_pdf.parent.mkdir(parents=True, exist_ok=True)
        shutil.copy2(pdf, target_pdf)

        png_target = target_pdf.with_suffix(".png")
        if not png_target.exists():
            ok = run_convert(target_pdf, png_target)
            if ok:
                print(f"[info] Created PNG mirror: {png_target}")

def mirror_tables() -> None:
    """
    Mirror all CSVs from `parameters/` (including subfolders) into `report/tables/`.
    """
    print("[info] Mirroring tables into report/tables ...")
    for csv in PARAM_PATH.rglob("*.csv"):
        rel = csv.relative_to(PARAM_PATH)
        target_csv = REPORT_TABLES / rel
        target_csv.parent.mkdir(parents=True, exist_ok=True)
        shutil.copy2(csv, target_csv)

def run_convert(pdf: Path, png: Path) -> bool:
    """
    Convert a PDF to PNG using ImageMagick's `convert` command.
    Returns True on success, False on failure.

    Args:
        pdf (Path): Path to the input PDF file.
        png (Path): Path to the output PNG file.

    Returns:
        bool: True if conversion succeeded, False otherwise.
    """
    try:
        subprocess.run(
            ["convert", "-density", "300", str(pdf), "-quality", "95", str(png)],
            check=True,
            stdout=subprocess.DEVNULL,
            stderr=subprocess.DEVNULL,
        )
        return True
    except FileNotFoundError:
        print("[warn] ImageMagick 'convert' not found – skipping PNG creation.")
        return False
    except subprocess.CalledProcessError:
        print(f"[warn] Failed to convert {pdf} → {png}")
        return False

def summarize_alpha(df: pd.DataFrame) -> str:
    """
    Summary for alphamcherry.csv table.
    """
    if "alpha" not in df.columns:
        for alt in ("alphamcherry", "alpha_mcherry", "a"):
            if alt in df.columns:
                df = df.rename(columns={alt: "alpha"})
                break
    if "alpha" not in df.columns:
        return "Could not find an 'alpha' column in alphamcherry.csv."

    vals = df["alpha"].dropna().to_numpy(float)
    if vals.size == 0:
        return "No finite alpha values found."

    n = len(vals)
    a_min, a_max = float(np.min(vals)), float(np.max(vals))
    a_med = float(np.median(vals))

    t_half = np.log(2.0) / vals
    t_min, t_max = float(np.min(t_half)), float(np.max(t_half))

    return textwrap.dedent(
        f"""
        The mCherry degradation / dilution rate alpha was fitted from Rev
        time courses for {n} instances. Alpha values span {a_min:.3g}–{a_max:.3g} h⁻¹
        (median {a_med:.3g} h⁻¹), corresponding to effective reporter half-lives
        of roughly {t_min:.2f}–{t_max:.2f} hours. This sets how quickly the
        reporter output can track changes in repressor activity.
        """
    ).strip()

def summarize_bootstrap_global(boot_all: Optional[pd.DataFrame]) -> str:
    """
    Summarize bootstrap distributions (uq_sensitivity/bootstrap_params_and_metrics_all.csv).
    """
    if boot_all is None or boot_all.empty:
        return "Bootstrap uncertainty results were not available."

    df = boot_all.copy()
    if "model" not in df.columns or "plasmid" not in df.columns:
        return "bootstrap_params_and_metrics_all.csv missing 'model'/'plasmid' columns."

    metrics = [c for c in ["dynamic_range", "t50", "t10_90", "overshoot", "auc", "y_end"] if c in df.columns]
    if not metrics:
        return "Bootstrap table did not contain expected metric columns."

    parts = []
    for model, sub in df.groupby("model"):
        lines = []
        for m in metrics:
            v = sub[m].to_numpy(float)
            v = v[np.isfinite(v)]
            if v.size < 50:
                continue
            q2, q50, q97 = np.percentile(v, [2.5, 50, 97.5])
            lines.append(f"{m}: median≈{q50:.3g}, 95%≈[{q2:.3g}, {q97:.3g}]")
        if lines:
            parts.append(f"{model}: " + "; ".join(lines))

    if not parts:
        return "Bootstrap summaries could not be extracted (too few finite values)."

    return (
        "Bootstrap uncertainty propagates measurement and fit uncertainty through the ODE model "
        "to generate distributions over predicted behaviors.\n\n"
        + "\n".join(parts)
    )

def summarize_delays(df_de: Optional[pd.DataFrame],
                     df_kd: Optional[pd.DataFrame]) -> str:
    """
    Summary for delays in derepression (df_de) and repression (df_kd).
    """
    parts: List[str] = []
    if df_de is not None and not df_de.empty:
        col = "d_rev" if "d_rev" in df_de.columns else df_de.columns[-1]
        vals = df_de[col].dropna().to_numpy(float)
        if vals.size:
            parts.append(
                f"Derepression delays (REV) range from {vals.min():.2f}–{vals.max():.2f} h."
            )
    if df_kd is not None and not df_kd.empty:
        col = "d_rev" if "d_rev" in df_kd.columns else df_kd.columns[-1]
        vals = df_kd[col].dropna().to_numpy(float)
        if vals.size:
            parts.append(
                f"Repression delays (KD) range from {vals.min():.2f}–{vals.max():.2f} h."
            )
    if not parts:
        return "No delay information was available for derepression or repression."
    return (
            "The ODE models required an explicit onset delay to align simulations with\n"
            + "the experimental time courses. "
            + " ".join(parts)
            + " These delays capture upstream processes such as degrader uptake, degron\n"
            + "processing and the time required until the effective nuclear repressor\n"
            + "concentration begins to change."
    )

def summarize_design_space(design_df: Optional[pd.DataFrame],
                           top10_df: Optional[pd.DataFrame]) -> str:
    """
    Summary for design-space scan and top-10 candidate selection.
    """
    if design_df is None or design_df.empty:
        return "Design-space scan results were not found."

    txt = []

    for col in ["dynamic_range", "t50", "overshoot"]:
        if col not in design_df.columns:
            continue
        vals = design_df[col].dropna().to_numpy(float)
        if not vals.size:
            continue
        txt.append(
            f"{col} spans {vals.min():.3g}–{vals.max():.3g} "
            f"(median {np.median(vals):.3g})."
        )

    if top10_df is not None and not top10_df.empty:
        # Try to detect 'score' and 'design_id' or similar
        score_col = "Score" if "Score" in top10_df.columns else None
        if score_col is None:
            for c in top10_df.columns:
                if c.lower().startswith("score"):
                    score_col = c
                    break
        id_col = "design_id" if "design_id" in top10_df.columns else None
        if id_col is None and "plasmid" in top10_df.columns:
            id_col = "plasmid"

        if score_col is not None and id_col is not None:
            best = top10_df.sort_values(score_col, ascending=False).iloc[0]
            best_label = best[id_col]
            best_score = float(best[score_col])
            txt.append(
                f"The top-ranked candidate according to the dynamic-range/"
                f"response-time score is {best_label} (Score ≈ {best_score:.3g})."
            )

    intro = (
        "Using the measured parameter ranges as priors, a synthetic design space\n"
        "of repression ODE models was explored. For each sampled parameter set,\n"
        "the reporter trajectory Y(t) was summarized into dynamic range, response\n"
        "times (t50) and overshoot. This allows ranking hypothetical designs that\n"
        "would be difficult to realize or test exhaustively in the wet lab.\n"
    )

    if not txt:
        return intro + "Summary statistics could not be extracted from the table."
    return intro + " ".join(txt)

def summarize_half_times(df: pd.DataFrame, label: str) -> str:
    """
    Generic summary for up/down half-times table.
    """
    cols = df.columns
    if "halftime" in cols:
        col = "halftime"
    elif "t_up" in cols:
        col = "t_up"
    elif "t_down" in cols:
        col = "t_down"
    else:
        return f"No recognizable half-time column found in {label}."

    df_use = df.dropna(subset=[col])
    if df_use.empty:
        return f"No finite half-times available for {label}."

    n = len(df_use)
    v = df_use[col].to_numpy(float)
    v_min, v_max = float(np.min(v)), float(np.max(v))
    v_med = float(np.median(v))

    fastest_row = df_use.iloc[int(np.argmin(v))]
    slowest_row = df_use.iloc[int(np.argmax(v))]
    fastest_pl = fastest_row.get("plasmid", "NA")
    slowest_pl = slowest_row.get("plasmid", "NA")

    return textwrap.dedent(
        f"""
        {label} ({n} constructs) show half-times ranging from {v_min:.2f} h
        to {v_max:.2f} h (median {v_med:.2f} h).
        The fastest response is observed for {fastest_pl} (~{v_min:.2f} h),
        while the slowest tuner is {slowest_pl} (~{v_max:.2f} h).
        Together, this defines the dynamic range of how quickly degron–Cas
        systems approach their new steady state after a perturbation.
        """
    ).strip()

def summarize_hill(df: pd.DataFrame) -> str:
    """
    Summary for Hill_parameters.csv table.
    """
    df = df.copy()
    cols = {c: c.strip() for c in df.columns}
    df = df.rename(columns=cols)

    # accept K/k in any capitalization
    if "k" not in df.columns:
        for alt in ("K", "k", "midpoint", "ec50"):
            if alt in df.columns:
                df = df.rename(columns={alt: "k"})
                break

    n_rows = len(df)
    k_vals = df["k"].to_numpy(float)
    n_vals = df["n"].to_numpy(float)

    k_min, k_max = float(np.min(k_vals)), float(np.max(k_vals))
    n_min, n_max = float(np.min(n_vals)), float(np.max(n_vals))
    n_med = float(np.median(n_vals))

    steep_row = df.iloc[int(np.argmax(n_vals))]
    graded_row = df.iloc[int(np.argmin(n_vals))]

    steep_pl = steep_row.get("plasmid", "NA")
    steep_n = float(steep_row["n"])

    graded_pl = graded_row.get("plasmid", "NA")
    graded_n = float(graded_row["n"])

    return textwrap.dedent(
        f"""
        Steady-state dose–response fits (Hill curves) were obtained for {n_rows}
        constructs. Effective midpoints K span approximately {k_min:.3g}–{k_max:.3g},
        indicating how much dTAG-13 is required to tune the system into its
        responsive regime.

        Hill exponents n range from {n_min:.2f} to {n_max:.2f} (median {n_med:.2f}).
        The steepest, most switch-like response is seen in {steep_pl} (n ≈ {steep_n:.2f}),
        whereas {graded_pl} (n ≈ {graded_n:.2f}) behaves more like a graded analog tuner.
        Together, these parameters quantify how sharply degron–Cas constructs
        transition from "off" to "on" as the degrader dose increases.
        """
    ).strip()

def summarize_noise(noise_ts: Optional[pd.DataFrame],
                    noise_hier: Optional[pd.DataFrame]) -> str:
    """
    Summary for single-cell noise tables.
    """
    if noise_ts is None or noise_ts.empty:
        return "Single-cell noise tables were not available."

    df = noise_ts.copy()
    text = []

    for ch, mean_col, cv_col in [
        ("BFP", "mean_BFP", "cv2_BFP"),
        ("mCherry", "mean_mCherry", "cv2_mCherry"),
    ]:
        if mean_col not in df.columns or cv_col not in df.columns:
            continue
        hi = df.groupby("plasmid")[cv_col].median().sort_values(ascending=False)
        lo = df.groupby("plasmid")[cv_col].median().sort_values(ascending=True)
        if hi.empty:
            continue

        noisy_pl = hi.index[0]
        quiet_pl = lo.index[0]
        noisy_cv = float(hi.iloc[0])
        quiet_cv = float(lo.iloc[0])

        text.append(
            f"For {ch}, median CV² across time and replicates spans roughly "
            f"{lo.iloc[0]:.3g}–{hi.iloc[-1]:.3g}. The noisiest construct "
            f"is {noisy_pl} (CV² ≈ {noisy_cv:.3g}), whereas {quiet_pl} shows "
            f"the quietest expression (CV² ≈ {quiet_cv:.3g})."
        )

    if not text:
        return "Could not extract channel-specific noise summaries."

    prefix = (
        "Single-cell flow cytometry events were propagated through the same gates\n"
        "as the kinetic analysis, and CV² was computed for both repressor (BFP)\n"
        "and reporter (mCherry). This quantifies how noisy each construct is at\n"
        "the single-cell level, beyond mean-level behavior.\n"
    )
    return prefix + "\n".join(text)

def summarize_sobol_global(sobol_all: Optional[pd.DataFrame]) -> str:
    """
    Summarize global Sobol indices (uq_sensitivity/sobol_indices_all.csv).
    Expected cols: model, metric, param, ST.
    """
    if sobol_all is None or sobol_all.empty:
        return "Global Sobol sensitivity results were not available."

    df = sobol_all.copy()
    need = {"model", "metric", "param", "ST"}
    if not need.issubset(df.columns):
        return f"sobol_indices_all.csv missing expected columns: {sorted(need)}."

    agg = (df.groupby(["model", "metric", "param"], as_index=False)["ST"]
             .mean()
             .rename(columns={"ST": "ST_mean"}))

    lines = []
    for (model, metric), sub in agg.groupby(["model", "metric"]):
        sub = sub.sort_values("ST_mean", ascending=False).head(3)
        tops = ", ".join([f"{r['param']} (ST≈{r['ST_mean']:.2f})" for _, r in sub.iterrows()])
        lines.append(f"{model} / {metric}: top drivers are {tops}.")

    return (
        "Sobol total-order sensitivity (ST) quantifies which fitted parameters most strongly "
        "control each performance metric under the ODE model.\n\n"
        + "\n".join(lines)
    )