 codex/symmetric-gaussian-overlap-log-ratio.py

"""Exact symmetric Gaussian-overlap diagnostic and its inverse."""
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def main():
    q, common = 1.5, .75
    h = np.log(q)
    u = np.linspace(-2, 2, 501)
    plus, minus = common/np.cosh(u+h)**2, common/np.cosh(u-h)**2
    ratio = np.log(plus/minus)
    analytic = -4*np.arctanh(np.tanh(u)*np.tanh(h))
    inverted = np.arctanh(-np.tanh(ratio/4)/np.tanh(h))
    assert np.allclose(ratio, analytic)
    assert np.allclose(inverted, u)
    assert np.all(np.diff(ratio) < 0)
    assert np.all(np.abs(ratio) < 4*h)
    assert np.allclose(ratio, -ratio[::-1])
    point = np.log(.8)
    vp, vm = common/np.cosh(point+h)**2, common/np.cosh(point-h)**2
    d = np.log(vp/vm)
    assert np.allclose([vp, vm, d], [.7256113947863477, .5172351863603165, .3385169256384606])
    s = -np.tanh(d/4)/np.tanh(h)
    derivative = -(1-np.tanh(d/4)**2)/(4*np.tanh(h)*(1-s*s))
    assert np.isclose(derivative, -.6780384615384614)
    fig, (raw, log) = plt.subplots(1, 2, figsize=(9.4, 4.2), dpi=100,
                                  layout="constrained")
    fig.set_facecolor("white")
    raw.plot(u, plus, color="#1558a0", lw=2, label="V₊: radius multiplied by 1.5")
    raw.plot(u, minus, "--", color="#b94b16", lw=2, label="V₋: radius divided by 1.5")
    raw.scatter([point, point], [vp, vm], color=["#1558a0", "#b94b16"], zorder=4)
    raw.axvline(0, color="#aaaaaa", lw=.8)
    raw.set(xlabel="Log radius ratio u = log ρ", ylabel="Scaled efficiency",
            title="Two readings with common scale A = 0.75",
            xlim=(-2, 2), ylim=(0, .82))
    raw.legend(fontsize=8, loc="lower center")
    log.plot(u, ratio, color="#1558a0", lw=2)
    for value in [-4*h, 4*h]:
        log.axhline(value, color="#b94b16", ls="--", lw=1)
    log.axhline(0, color="#aaaaaa", lw=.8)
    log.axvline(0, color="#aaaaaa", lw=.8)
    log.scatter([point], [d], color="#b94b16", zorder=4)
    log.annotate("ρ = 0.8\nD = 0.33852", (point, d), (point+.35, d+.4),
                 arrowprops=dict(arrowstyle="->", color="#555555"), fontsize=9)
    log.text(1.9, 4*h+.08, "+4h", ha="right", color="#b94b16", fontsize=9)
    log.text(-1.9, -4*h-.17, "−4h", color="#b94b16", fontsize=9)
    log.set(xlabel="Log radius ratio u = log ρ", ylabel="D = log(V₊ / V₋)",
            title="The log-ratio removes the common scale",
            xlim=(-2, 2), ylim=(-1.95, 1.95))
    for ax in [raw, log]:
        ax.grid(alpha=.18)
    fig.savefig(Path(__file__).with_suffix(".png").name, facecolor="white")
    plt.close(fig)


if __name__ == "__main__":
    main()