"""Invariant-domain reconstruction toy implementation for 2-D ideal-gas Euler states."""
import numpy as np

GAMMA = 1.4

def pressure(w, gamma=GAMMA):
    w = np.asarray(w, dtype=float)
    rho = w[..., 0]
    kinetic = 0.5 * np.sum(w[..., 1:3] ** 2, axis=-1) / rho
    return (gamma - 1.0) * (w[..., 3] - kinetic)

def admissible(w, rho_floor=1e-10, p_floor=1e-10):
    w = np.asarray(w)
    return bool(np.all(w[..., 0] >= rho_floor) and np.all(pressure(w) >= p_floor))

def density_theta_bound(w0, wr, rho_floor=1e-10):
    """Largest theta satisfying affine density floor, assuming w0 is admissible."""
    r0, rr = float(w0[0]), float(wr[0])
    if rr >= rho_floor:
        return 1.0
    return float(np.clip((r0-rho_floor)/(r0-rr), 0.0, 1.0))

def limited_state(w0, wr, rho_floor=1e-10, p_floor=1e-10, iterations=60):
    """Return state and largest tested theta on segment w0 + theta*(wr-w0)."""
    w0, wr = np.asarray(w0, float), np.asarray(wr, float)
    if not admissible(w0, rho_floor, p_floor):
        raise ValueError("cell average must be admissible")
    hi = density_theta_bound(w0, wr, rho_floor)
    # Pressure positivity is checked by bisection; hi is already density-safe.
    if admissible(w0 + hi*(wr-w0), rho_floor, p_floor):
        return w0 + hi*(wr-w0), hi
    lo = 0.0
    for _ in range(iterations):
        mid = 0.5*(lo+hi)
        if admissible(w0 + mid*(wr-w0), rho_floor, p_floor):
            lo = mid
        else:
            hi = mid
    return w0 + lo*(wr-w0), lo

def reconstruct_face(cell, raw, floor=1e-8):
    return limited_state(cell, raw, floor, floor)

def minmod(a, b):
    return np.sign(a)*np.minimum(np.abs(a), np.abs(b)) if a*b > 0 else 0.0

def scalar_reconstruction(cells, raw_faces):
    """Toy scalar analogue: classical minmod vs hard invariant reconstruction."""
    classical=[]; safe=[]; thetas=[]
    for c, r in zip(cells, raw_faces):
        # use density-like scalar component for a transparent baseline
        slope=minmod(r[0]-c[0], c[0]-max(c[0]-0.2, 1e-5))
        classical.append(c[0] + slope)
        s,t=limited_state(c,r,1e-8,1e-8)
        safe.append(s[0]); thetas.append(t)
    return np.array(classical), np.array(safe), np.array(thetas)