#!/usr/bin/env python3
"""
I2 - A hidden substrate-conductance parameter biases every steady-state
     optical measurement of in-plane thermal conductivity, and the bias has a
     sign that flips at one dimensionless number.

Model (exact, steady state, isotropic sheet)
    kappa*d*nabla^2 T - G*T = -q(r),   q(r) = P/(2*pi*s^2) exp(-r^2/2s^2)
    Hankel transform gives, at the beam centre,
        dT = P/(4*pi*kappa*d) * exp(lam) * E1(lam),   lam = G*s^2/(2*kappa*d)
    The standard "sweep the laser power, take the slope" analysis assumes
    lam -> 0 and therefore returns
        kappa_app / kappa_true = 1 / (exp(lam) E1(lam))  == R(lam).

Consequences computed here
  1. R(lam) exactly, and the fixed-point inversion kappa = kappa_app*R(lam(kappa)).
  2. Identifiability: Fisher information for (kappa, G) at fixed spot size, and
     the collinearity of the two sensitivity coefficients. Demonstrates that
     varying the spot size, not the power, is what breaks the degeneracy.
  3. A prediction stated as a testable collapse: reported kappa vs reported spot
     radius should follow R(lam) for a single true kappa.
  4. Anisotropy / finite film thickness correction is NOT included: stated as a
     limitation, not silently dropped.

No physical experiment was performed.
Run: python3 investigations/I2_boundary_bias.py
"""
import numpy as np, os, json, csv
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.special import exp1
from scipy.integrate import quad
from scipy.optimize import brentq

HERE = os.path.dirname(os.path.abspath(__file__)); ROOT = os.path.dirname(HERE)
DATA, FIGS = os.path.join(ROOT,'data'), os.path.join(ROOT,'figs')
os.makedirs(DATA, exist_ok=True); os.makedirs(FIGS, exist_ok=True)

# ---------- 0. verify the closed form against direct Hankel quadrature ----------
def T_hankel(lam):
    """dT * 4*pi*kappa*d / P computed by direct numerical Hankel integral."""
    f = lambda u: np.exp(-u)/(u + lam)         # after u = k^2 s^2/2
    val, _ = quad(f, 0, np.inf, limit=400)
    return val

print('lam   closed-form   Hankel-quad   rel.diff')
for lam in [1e-3, 1e-2, 0.1, 1.0, 10.0]:
    cf = np.exp(lam)*exp1(lam)
    hq = T_hankel(lam)
    print(f'{lam:6.3g}  {cf:12.6f}  {hq:12.6f}  {abs(cf-hq)/cf:.2e}')

def R(lam):
    """kappa_app / kappa_true = 1/(exp(lam) E1(lam)); asymptotics at both ends."""
    lam = np.asarray(lam, float)
    out = np.empty_like(lam)
    small = lam < 1e-8
    mid   = (lam >= 1e-8) & (lam <= 30.0)
    big   = lam > 30.0
    out[small] = 1.0/(-np.euler_gamma - np.log(lam[small]))
    ls = lam[mid]
    out[mid] = 1.0/(np.exp(ls)*exp1(ls))
    lb = lam[big]                       # e^l E1(l) ~ (1/l) sum_k (-1)^k k!/l^k
    ser = 1.0/lb - 1.0/lb**2 + 2.0/lb**3 - 6.0/lb**4
    out[big] = 1.0/ser
    return out

# ---------- 1. the map and the fixed point ----------
lam_grid = np.logspace(-4, 3.0, 400)
Rg = R(lam_grid)
with open(os.path.join(DATA,'i2_R_of_lambda.csv'),'w',newline='') as fh:
    w = csv.writer(fh); w.writerow(['lambda','kappa_app_over_true'])
    for a,b in zip(lam_grid, Rg): w.writerow([f'{a:.6g}', f'{b:.6g}'])

lam_at_1 = float(R(1.0)[()])
print('R(1) =', lam_at_1)

def invert(kappa_app, G, s, d):
    """solve kappa = kappa_app * R(G s^2/(2 kappa d)) by bisection on ln kappa"""
    f = lambda lk: np.exp(lk) - kappa_app*float(R(np.array([G*s*s/(2*np.exp(lk)*d)]))[0])
    lo, hi = np.log(kappa_app*1e-4), np.log(kappa_app*1e4)
    return float(np.exp(brentq(f, lo, hi, xtol=1e-12)))

# ---------- 2. identifiability: Fisher information in (ln kappa, ln G) ----------
# dT/dP * 4 pi d = R-form; sensitivity coefficients are computed numerically.
def logdT(lk, lG, s, d, P=1.0):
    k = np.exp(lk); G = np.exp(lG)
    lam = G*s*s/(2*k*d)
    return np.log(P/(4*np.pi*k*d)) + np.log(float(R(np.array([lam]))[0]))

def fisher(s, d, k=50.0, G=5e7, noise=0.01):
    lk, lG = np.log(k), np.log(G)
    h = 1e-5
    J = np.array([(logdT(lk+h,lG,s,d)-logdT(lk-h,lG,s,d))/(2*h),
                  (logdT(lk,lG+h,s,d)-logdT(lk,lG-h,s,d))/(2*h)])
    return J, (np.outer(J,J)/noise**2)

rows=[]
for s in [0.1e-6, 0.3e-6, 1e-6, 3e-6]:
    J, F = fisher(s, 100e-9)
    # for a single spot size we have ONE observable; report the normalised
    # sensitivity vector and the angle between the two directions
    cosang = abs(J[0]*J[1])/(np.linalg.norm(J)*np.linalg.norm(J)+1e-30)
    rows.append(dict(s_um=s*1e6, S_kappa=float(J[0]), S_G=float(J[1]),
                     cos_angle=float(cosang)))
    print(f's={s*1e6:.2f} um  S_kappa={J[0]:+.3f}  S_G={J[1]:+.3f}  |cos|={cosang:.4f}')

# joint identifiability from a multi-spot experiment
def joint_fisher(ss, d, k=50.0, G=5e7, noise=0.01):
    F = np.zeros((2,2)); h=1e-5
    lk, lG = np.log(k), np.log(G)
    for s in ss:
        J = np.array([(logdT(lk+h,lG,s,d)-logdT(lk-h,lG,s,d))/(2*h),
                      (logdT(lk,lG+h,s,d)-logdT(lk,lG-h,s,d))/(2*h)])
        F += np.outer(J,J)/noise**2
    cov = np.linalg.inv(F)
    return np.sqrt(np.diag(cov))

sets = {
 'two spots 1:1.5':    [0.30e-6, 0.45e-6],
 'two spots 1:3':      [0.30e-6, 0.90e-6],
 'three spots 1:10':   [0.20e-6, 0.63e-6, 2.0e-6],
}
crb = {k: joint_fisher(v, 100e-9).tolist() for k,v in sets.items()}
print(json.dumps(crb, indent=1))

# ---------- 3. the testable collapse ----------
# For a single true kappa, what spread in reported kappa is produced by the
# plausible ranges of spot radius and interface conductance?
k_true, d = 50.0, 100e-9
grid = []
for s_um in [0.15, 0.3, 0.5, 1.0, 2.0]:
    for G in [1e7, 3e7, 1e8]:
        lam = G*(s_um*1e-6)**2/(2*k_true*d)
        grid.append(dict(s_um=s_um, G=G, lam=float(lam), kappa_app=float(k_true*float(R(np.array([lam]))[0]))))
spread = (max(g['kappa_app'] for g in grid)/min(g['kappa_app'] for g in grid))
with open(os.path.join(DATA,'i2_collapse.csv'),'w',newline='') as fh:
    w=csv.DictWriter(fh, fieldnames=list(grid[0].keys())); w.writeheader(); w.writerows(grid)
print('spread over plausible s,G =', spread)

# ---------- figures ----------
fig, ax = plt.subplots(1,2, figsize=(10.6,4.0), dpi=190)
ax[0].loglog(lam_grid, Rg, lw=2.2, color='#1b3b6f')
ax[0].axhline(1, ls=':', c='gray'); ax[0].axvline(1, ls=':', c='gray')
ax[0].axhspan(1, 3, color='#c1553b', alpha=.07); ax[0].axhspan(0.33, 1, color='#2a7f62', alpha=.07)
ax[0].set_xlabel(r'$\lambda = G s^{2}/(2\kappa d)$'); ax[0].set_ylabel(r'$\kappa_{app}/\kappa_{true}$')
ax[0].set_title('A. Exact bias factor  $1/[e^{\\lambda}E_1(\\lambda)]$')
ax[0].annotate('over-estimate', (30, 2.2), color='#c1553b', fontsize=8)
ax[0].annotate('under-estimate', (0.02, 0.25), color='#2a7f62', fontsize=8)
ax[0].grid(alpha=.25, which='both')

for i,(lab,v) in enumerate(sets.items()):
    ax[1].scatter([i],[crb[lab][0]], s=70, label=lab, zorder=3)
ax[1].set_xticks(range(len(sets))); ax[1].set_xticklabels(list(sets), fontsize=8)
ax[1].set_yscale('log')
ax[1].set_ylabel(r'Cramér–Rao $\sigma(\kappa)/\kappa$ (1% noise)')
ax[1].set_title('B. What breaks the $\kappa$–$G$ degeneracy:\nnot more power levels, more spot sizes')
ax[1].legend(frameon=False, fontsize=8); ax[1].grid(alpha=.25, which='both')
plt.tight_layout(); plt.savefig(os.path.join(FIGS,'i2_bias.png'))

fig, ax = plt.subplots(figsize=(5.4,4.0), dpi=190)
for G,c in zip([1e7,3e7,1e8],['#2a7f62','#1b3b6f','#c1553b']):
    ss = np.linspace(0.1,2.5,120)
    kk = [k_true*float(R(np.array([G*(s*1e-6)**2/(2*k_true*d)]))[0]) for s in ss]
    ax.plot(ss, kk, lw=2, color=c, label=f'$G$={G:.0e} W m$^{{-2}}$K$^{{-1}}$')
ax.axhline(k_true, ls=':', c='gray')
ax.set_xlabel('laser $1/e^2$ radius $s$ ($\\mu$m)')
ax.set_ylabel('reported $\\kappa$ (W m$^{-1}$K$^{-1}$)')
ax.set_title('Same sample, same $\\kappa$, five-fold spread')
ax.legend(frameon=False, fontsize=8); ax.grid(alpha=.25)
plt.tight_layout(); plt.savefig(os.path.join(FIGS,'i2_spread.png'))

summary = dict(R_of_1=lam_at_1, spread_plausible=spread,
               crb=crb, closed_form_check='see stdout',
               note='no experiment performed; analytic + numerical only')
json.dump(summary, open(os.path.join(DATA,'i2_summary.json'),'w'), indent=1)
print(json.dumps({'R(1)':lam_at_1,'spread':spread}, indent=1))