from __future__ import annotations
import numpy as np

def jacobian_numeric(beta_fn, x, h: float = 1e-7, method: str = "central") -> np.ndarray:
    """
    Finite-difference Jacobian with scale-aware step: dx_j = h*max(1,|x_j|)
    """
    x = np.asarray(x, dtype=float)
    f0 = np.asarray(beta_fn(x), dtype=float)
    n = x.size
    m = f0.size
    J = np.zeros((m, n), dtype=float)

    for j in range(n):
        dx = np.zeros_like(x)
        dx[j] = h * max(1.0, abs(x[j]))
        if method == "central":
            fp = np.asarray(beta_fn(x + dx), dtype=float)
            fm = np.asarray(beta_fn(x - dx), dtype=float)
            J[:, j] = (fp - fm) / (2.0 * dx[j])
        elif method == "forward":
            fp = np.asarray(beta_fn(x + dx), dtype=float)
            J[:, j] = (fp - f0) / dx[j]
        else:
            raise ValueError("method must be 'central' or 'forward'")
    return J

def compare_jacobians(Ja: np.ndarray, Jn: np.ndarray, eps: float = 1e-30):
    """
    rel_err_ij = |Ja-Jn| / max(eps, |Ja|, |Jn|)
    """
    Ja = np.asarray(Ja, dtype=float)
    Jn = np.asarray(Jn, dtype=float)
    denom = np.maximum(eps, np.maximum(np.abs(Ja), np.abs(Jn)))
    rel = np.abs(Ja - Jn) / denom
    max_rel = float(np.max(rel))
    max_abs = float(np.max(np.abs(Ja - Jn)))
    idx = np.unravel_index(int(np.argmax(rel)), rel.shape)
    i, j = int(idx[0]), int(idx[1])
    worst = {
        "i": i, "j": j,
        "analytic": float(Ja[i, j]),
        "numeric":  float(Jn[i, j]),
        "rel_err":  float(rel[i, j]),
        "abs_err":  float(abs(Ja[i, j] - Jn[i, j])),
    }
    return max_rel, max_abs, worst
