import numpy as np
from scipy.special import betaln
from scipy.optimize import minimize
w=lambda c:c/(1-c+1e-15)
ab=lambda f,c:(w(c)*f,w(c)*(1-f))
h2=lambda a1,b1,a2,b2:1-np.exp(betaln((a1+a2)/2,(b1+b2)/2)-0.5*(betaln(a1,b1)+betaln(a2,b2)))
rao=lambda a1,b1,a2,b2:2*np.arcsinh(np.sqrt(np.clip(h2(a1,b1,a2,b2),0,0.999)/(1-np.clip(h2(a1,b1,a2,b2),0,0.999)+1e-15)))
pairs=[((0.8,0.7),(0.85,0.6)),((0.9,0.9),(0.7,0.8)),((0.5,0.5),(0.6,0.5))]
print("2-agent: NAL rev vs Frechet mean")
for (f1,c1),(f2,c2) in pairs:
    a1,b1=ab(f1,c1); a2,b2=ab(f2,c2)
    w1,w2_=w(c1),w(c2)
    fr=(w1*f1+w2_*f2)/(w1+w2_); cr=(w1+w2_)/(w1+w2_+1)
    betas=[(a1,b1),(a2,b2)]
    cost=lambda x,bs=betas:sum(rao(w(x[1])*x[0],w(x[1])*(1-x[0]),a,b)**2 for a,b in bs)
    res=minimize(cost,[0.7,0.6],bounds=[(0.02,0.98),(0.02,0.98)],method="L-BFGS-B")
    print(f"  NAL f={fr:.4f} c={cr:.4f} | Frechet f={res.x[0]:.4f} c={res.x[1]:.4f} | df={abs(fr-res.x[0]):.4f} dc={abs(cr-res.x[1]):.4f}")