"""
Figure 1 of Paper IX: the optimised weights lambda_d (d <= 240, lam_L0174.npy) against Polymath's Euler mollifier
prod_{p <= 13} (1 - b_p p^{-s*}), whose coefficients are lambda_d = prod_{p | d} (-b_p) on the divisors d of 30030.
Linear axes.  Run:  python weights_fig.py   (writes weights.png next to the paper).  P. Caceres, October 2026.
"""
import itertools
import math
import os

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
T = 0.154


def b(n):
    return math.exp(T / 4 * math.log(n) ** 2)


def main():
    lam = np.load(os.path.join(HERE, "lam_L0174.npy"))
    d = np.arange(1, len(lam) + 1)
    euler = {}
    primes = [2, 3, 5, 7, 11, 13]
    for r in range(len(primes) + 1):
        for subset in itertools.combinations(primes, r):
            m = math.prod(subset)
            if m <= len(lam):
                euler[m] = math.prod(-b(p) for p in subset)
    fig, ax = plt.subplots(figsize=(9, 3.6))
    ax.bar(d, lam, width=0.8, color="#3060a0", label=r"optimised $\lambda_d$ ($d\leq240$)")
    ks = sorted(euler)
    ax.plot(ks, [euler[k] for k in ks], "o", mfc="none", mec="#c03030", ms=5,
            label=r"Euler mollifier $\prod_{p\mid d}(-b_p)$, $p\leq13$")
    ax.axhline(0, color="k", lw=0.5)
    ax.set_xlabel("$d$")
    ax.set_ylabel("weight")
    ax.set_title(r"Mollifier weights, $t_0=0.154$, $y_0=0.2$, $x=6\cdot10^{12}$")
    ax.legend(loc="lower left", fontsize=8)
    fig.tight_layout()
    out = os.path.join(os.path.dirname(HERE), "weights.png")
    fig.savefig(out, dpi=200)
    print("wrote", out)


if __name__ == "__main__":
    main()