# (c) 2026 William J. Dwyer. MIT License. Cite Zenodo DOI 10.5281/zenodo.20709963.
# Figure for "The Log-Domain Distributional Atlas" — two panels:
#  (A) the log-domain Cullen-Frey diagram: (log-skewness, log-excess-kurtosis) loci of the families.
#  (B) the CV -> Var(ln X) bridges for the same families.
import numpy as np
from scipy.special import polygamma, gammaln
from scipy import optimize
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt

plt.rcParams.update({"axes.labelsize":11,"xtick.labelsize":9.5,"ytick.labelsize":9.5,
                     "legend.fontsize":8.5,"font.size":10.5})
LN="#1f6f8b"; GM="#b45f06"; IG="#7a3f8c"; WB="#1a7a4b"; LL="#b3261e"; PA="#444a55"
GUMBEL=-12*np.sqrt(6)*1.2020569031595942/np.pi**3

fig,(axA,axB)=plt.subplots(1,2,figsize=(9.6,4.4))

# ---------- Panel A: log-domain Cullen-Frey ----------
k=np.linspace(1.0,60,400)                              # gamma: from exponential (k=1) to normal limit
g_sk=polygamma(2,k)/polygamma(1,k)**1.5; g_ku=polygamma(3,k)/polygamma(1,k)**2
a=np.linspace(2.2,60,400)                              # inverse-gamma: positive-skew mirror
ig_sk=-polygamma(2,a)/polygamma(1,a)**1.5; ig_ku=polygamma(3,a)/polygamma(1,a)**2
axA.plot(g_sk,g_ku,color=GM,lw=2.0,label="gamma (curve in k)")
axA.plot(ig_sk,ig_ku,color=IG,lw=2.0,ls=(0,(4,2)),label="inverse-gamma (curve)")
axA.scatter([0],[0],color=LN,s=70,zorder=5,label="lognormal  (0, 0)")
axA.scatter([0],[1.2],color=LL,s=70,marker="s",zorder=5,label="log-logistic  (0, 1.2)")
axA.scatter([GUMBEL],[2.4],color=WB,s=80,marker="^",zorder=5,label="Weibull  (−1.14, 2.4)")
axA.scatter([2],[6],color=PA,s=80,marker="D",zorder=5,label="Pareto tail  (+2, 6)")
axA.annotate("exponential\n(gamma k=1 = Weibull)", xy=(GUMBEL,2.4), xytext=(-1.05,4.0),
             fontsize=8, color=WB, ha="center", arrowprops=dict(arrowstyle="->",color=WB))
axA.axvline(0,color="#cfd6da",lw=1,zorder=0)
axA.set_xlabel("log-skewness  $\\gamma_1(\\ln X)$"); axA.set_ylabel("log-excess-kurtosis  $\\gamma_2(\\ln X)$")
axA.set_title("A. Log-domain Cullen–Frey atlas", fontweight="bold", fontsize=11)
axA.set_xlim(-1.6,2.5); axA.set_ylim(-0.6,6.6)
axA.legend(loc="upper center", frameon=False, ncol=1, fontsize=8)
axA.grid(alpha=.22)

# ---------- Panel B: CV -> Var(ln X) bridges ----------
cv=np.linspace(0.05,1.2,300)
def k_weibull(c): return optimize.brentq(lambda kk: np.exp(gammaln(1+2/kk)-2*gammaln(1+1/kk))-1-c**2,1e-2,1e4)
def beta_ll(c):
    f=lambda b:( (2*np.pi/b)/np.sin(2*np.pi/b))/((np.pi/b)/np.sin(np.pi/b))**2-1-c**2
    return optimize.brentq(f,2.0001,1e4)
ln_v=np.log1p(cv**2)
gm_v=polygamma(1,1/cv**2)
ig_v=polygamma(1,2+1/cv**2)
ap=1+np.sqrt(1+1/cv**2); pa_v=1/ap**2
wb_v=np.array([np.pi**2/(6*k_weibull(c)**2) for c in cv])
ll_v=np.array([np.pi**2/(3*beta_ll(c)**2) for c in cv])
axB.plot(cv,cv**2,color="#9aa3aa",lw=1.2,ls=":",label="$CV^2$ (reference)")
axB.plot(cv,gm_v,color=GM,lw=2.0,label="gamma  $\\psi'(1/CV^2)$")
axB.plot(cv,ig_v,color=IG,lw=1.8,ls=(0,(4,2)),label="inverse-gamma  $\\psi'(2+1/CV^2)$")
axB.plot(cv,ln_v,color=LN,lw=2.0,label="lognormal  $\\ln(1+CV^2)$")
axB.plot(cv,wb_v,color=WB,lw=1.8,label="Weibull  $\\pi^2/6k^2$")
axB.plot(cv,ll_v,color=LL,lw=1.8,label="log-logistic  $\\pi^2/3\\beta^2$")
axB.plot(cv,pa_v,color=PA,lw=1.8,label="Pareto  $1/\\alpha^2$")
axB.set_xlabel("coefficient of variation  CV"); axB.set_ylabel("log-scale variance  Var(ln X)")
axB.set_title("B. CV → Var(ln X) bridges", fontweight="bold", fontsize=11)
axB.legend(loc="upper left", frameon=False, fontsize=8)
axB.grid(alpha=.22)

fig.suptitle("The log-domain distributional atlas — triaging positive families by CV and log-cumulants",
             fontsize=11.5, fontweight="bold", y=1.01)
fig.tight_layout()
fig.savefig("fig_distribution_atlas.png", dpi=200, bbox_inches="tight",
            metadata={"Author":"William J. Dwyer (ORCID 0009-0004-0855-7222)","Title":"fig_distribution_atlas",
                      "Subject":"Companion figure for the Log-Domain Distributional Atlas; MIT; cite Zenodo 10.5281/zenodo.20709963"})
print("saved fig_distribution_atlas.png")
