#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""复现 Nature Communications 2026 (s41467-026-74406-0) Fig.1a
同源 E 蛋白结合抗体纵向衰减轨迹：原发（P，虚线）vs 继发（S，实线）
含分段线性拟合 + 95% bootstrap CI 阴影带。数据为依据原图轨迹模拟（明确标注）。
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "DejaVu Sans"

ISOTYPES = ["IgG", "IgA", "IgM", "IgG1", "IgG2", "IgG3", "IgG4"]
# knots: 30/90/180/540 days post symptom onset —— 依据原图读数模拟
P_MEAN = {
    "IgG":  [5.0, 5.5, 4.9, 5.1],
    "IgA":  [4.5, 2.5, 1.7, 2.2],
    "IgM":  [5.0, 0.0, -0.1, -0.2],
    "IgG1": [5.5, 6.3, 5.6, 5.6],
    "IgG2": [3.2, 5.0, 3.5, 4.0],
    "IgG3": [4.6, 3.0, 0.8, 0.7],
    "IgG4": [-1.0, -0.7, -0.9, 1.0],
}
S_MEAN = {
    "IgG":  [6.4, 6.1, 5.3, 5.2],
    "IgA":  [2.9, 2.3, 1.6, 1.5],
    "IgM":  [1.2, -0.5, -0.6, -0.5],
    "IgG1": [7.1, 6.9, 5.9, 5.7],
    "IgG2": [7.2, 6.5, 4.3, 4.0],
    "IgG3": [3.3, 1.8, 0.4, 0.2],
    "IgG4": [0.3, 0.4, 0.9, 1.6],
}
KNOTS = np.array([30, 90, 180, 540])
POS = {30: 0.0, 90: 1.0, 180: 2.0, 540: 4.0}  # 原图为分段间距坐标：90→180 短距，180→540 长距
XP = np.linspace(0, 4, 220)
rng = np.random.default_rng(20261020)

def interp(means):
    kp = [POS[k] for k in KNOTS]
    return np.interp(XP, kp, means)

def ci_band(mean):
    """模拟 bootstrap 95% CI 带：中心线±噪声包络"""
    w = 0.22 + 0.10 * np.abs(mean) / 8.0
    n = 400
    lo, hi = np.empty((n, len(XP))), np.empty((n, len(XP)))
    for i in range(n):
        jitter = rng.normal(0, w, size=len(XP))
        jitter += rng.normal(0, 0.08) * np.sin(XP / 2.0)
        lo[i] = mean + jitter - 0.18
        hi[i] = mean + jitter + 0.18
    return np.percentile(lo, 2.5, axis=0), np.percentile(hi, 97.5, axis=0)

COL_P, COL_S = "#C9A227", "#5A5A5A"   # P: 金色虚线, S: 灰色实线
BG_COL = "#E8E8E8"

fig, axes = plt.subplots(1, 7, figsize=(16.5, 4.6), sharey=True)
for ax, iso in zip(axes, ISOTYPES):
    pm, sm = interp(P_MEAN[iso]), interp(S_MEAN[iso])
    plo, phi = ci_band(pm)
    slo, shi = ci_band(sm)
    ax.axhspan(-2, 0, color=BG_COL, alpha=0.7, zorder=0)  # background 阴影区
    ax.fill_between(XP, plo, phi, color=COL_P, alpha=0.30, linewidth=0)
    ax.fill_between(XP, slo, shi, color=COL_S, alpha=0.30, linewidth=0)
    ax.plot(XP, pm, color=COL_P, linestyle=(0, (3, 2.5)), lw=1.8)
    ax.plot(XP, sm, color=COL_S, linestyle="-", lw=1.8)
    ax.set_title(iso, fontsize=14)
    ax.set_xlim(0, 4)
    ax.set_ylim(-2, 9)
    ax.set_xticks([0, 1, 4])
    ax.set_xticklabels(["30", "90", "540"])
    ax.tick_params(labelsize=11)
    ax.axhline(0, color="black", lw=0.8)
    for spine in ("top", "right"):
        ax.set_visible(False) if False else None
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

axes[0].set_ylabel("Log$_2$(MFI) - BG", fontsize=14)
fig.text(0.5, 0.02, "Days post symptoms onset", ha="center", fontsize=14)
fig.text(0.01, 0.95, "a", fontsize=24, fontweight="bold")
# 图例
from matplotlib.lines import Line2D
leg = [Line2D([0], [0], color=COL_P, lw=1.8, linestyle=(0, (3, 2.5)), label="P"),
       Line2D([0], [0], color=COL_S, lw=1.8, label="S")]
fig.legend(handles=leg, title="Infection group", loc="center left",
           bbox_to_anchor=(0.92, 0.55), fontsize=12, title_fontsize=13, frameon=False)
fig.subplots_adjust(left=0.055, right=0.90, bottom=0.16, top=0.88, wspace=0.12)
fig.savefig("/home/hatch/workspace/arbovirus-site/day-2026-10-20-figure/repro_py.png",
            dpi=150, bbox_inches="tight")
print("saved repro_py.png")