"""Seller value function: existence/uniqueness, Proposition 1 and 2 comparative statics, boundary cases."""
import numpy as np
import pytest

from src.model import Pareto, Exponential
from src.solver import seller_value, solve_sellers


def test_residual_and_uniqueness(params):
    sb = solve_sellers(8.0, params.k, params)
    assert np.abs(sb.resid).max() < 1e-10
    # Phi is strictly decreasing: perturbing V changes the sign of the residual
    for i, c in enumerate(params.c):
        V = sb.V[i]
        up = -c + params.delta * params.w0 + 8.0 * params.beta * float(params.G.H(V + 1e-3 + params.k)) - (params.r + params.delta) * (V + 1e-3)
        dn = -c + params.delta * params.w0 + 8.0 * params.beta * float(params.G.H(V - 1e-3 + params.k)) - (params.r + params.delta) * (V - 1e-3)
        assert up < 0 < dn


def test_prop1_monotonicity(params):
    sb = solve_sellers(8.0, params.k, params)
    assert np.all(np.diff(sb.V) < 0)          # reservation price falls with c
    assert np.all(np.diff(sb.eta) > 0)        # hazard rises with c
    assert np.all(np.diff(sb.pbar) < 0)       # expected price falls with c
    # analytic slope  V'(c) = -1/[(r+delta) + Lam beta (1-G)]  vs finite differences
    slope_an = -1.0 / ((params.r + params.delta) + 8.0 * params.beta * sb.acc)
    slope_fd = np.gradient(sb.V, params.c)
    assert np.allclose(slope_an[1:-1], slope_fd[1:-1], rtol=0.05)
    # |pbar'| <= |V'| under IFR
    assert np.all(np.abs(np.diff(sb.pbar)) <= np.abs(np.diff(sb.V)) + 1e-12)
    assert np.all(-1.0 / (params.r + params.delta) < slope_an) and np.all(slope_an < 0)


def test_prop1_participation_threshold(params):
    Lam = 8.0
    cbar = Lam * params.beta * float(params.G.H(params.w0 + params.k)) - params.r * params.w0
    p = params.copy(c=np.array([cbar - 1e-4, cbar + 1e-4]), f=np.array([0.5, 0.5]))
    sb = solve_sellers(Lam, params.k, p)
    assert sb.V[0] >= params.w0 - 1e-9 and sb.V[1] < params.w0


def test_prop2_arrival_rate(params):
    lams = np.array([2.0, 4.0, 8.0, 16.0])
    V = np.array([solve_sellers(L, params.k, params).V for L in lams])
    eta = np.array([solve_sellers(L, params.k, params).eta for L in lams])
    pbar = np.array([solve_sellers(L, params.k, params).pbar for L in lams])
    assert np.all(np.diff(V, axis=0) > 0)
    assert np.all(np.diff(pbar, axis=0) > 0)
    assert np.all(np.diff(eta, axis=0) > 0)                 # IFR => hazard increasing in Lam_s
    gap = V[:, 0] - V[:, -1]
    assert np.all(np.diff(gap) < 0)                         # reservation gap narrows with Lam_s
    # analytic dV/dLam vs finite difference
    sb = solve_sellers(8.0, params.k, params)
    an = params.beta * sb.H / ((params.r + params.delta) + 8.0 * params.beta * sb.acc)
    fd = (solve_sellers(8.0 + 1e-4, params.k, params).V - solve_sellers(8.0 - 1e-4, params.k, params).V) / 2e-4
    assert np.allclose(an, fd, rtol=1e-4)


def test_prop2_pareto_counterexample(params):
    """With Pareto G (DFR), d eta / d Lam_s < 0 iff c > delta*w0 + (r+delta)*U, independently of Lam_s (Proposition 2(ii))."""
    lams = (2.0, 5.0, 10.0, 20.0, 40.0)
    for w0, c in ((0.0, 0.05), (0.8, 0.2)):            # condition holds -> hazard decreasing in Lam_s
        p = params.copy(G=Pareto(a=1.5, lo=0.7), c=np.array([c]), f=np.array([1.0]), delta=0.1, w0=w0)
        assert c > p.delta * w0 + (p.r + p.delta) * p.k
        eta = np.array([solve_sellers(L, p.k, p).eta[0] for L in lams])
        assert np.all(np.diff(eta) < 0)
    for w0, c in ((0.0, 0.002), (0.8, 0.05)):          # condition fails -> hazard increasing
        p = params.copy(G=Pareto(a=1.5, lo=0.7), c=np.array([c]), f=np.array([1.0]), delta=0.1, w0=w0)
        assert c < p.delta * w0 + (p.r + p.delta) * p.k
        eta = np.array([solve_sellers(L, p.k, p).eta[0] for L in lams])
        assert np.all(np.diff(eta) > 0)


def test_prop6_exponential_elasticity(params):
    """epsilon(c) = (r+delta)/(r+delta+beta*eta) for exponential G."""
    p = params.copy(G=Exponential(rate=8.0, lo=0.7))
    L = 8.0
    sb = solve_sellers(L, p.k, p)
    fd = (np.log(solve_sellers(L * 1.001, p.k, p).eta) - np.log(solve_sellers(L / 1.001, p.k, p).eta)) / (2 * np.log(1.001))
    an = (p.r + p.delta) / ((p.r + p.delta) + p.beta * sb.eta)
    assert np.allclose(fd, an, rtol=2e-3)
    assert np.all(np.diff(an) < 0)


def test_boundary_cases(params):
    # zero urgency, no withdrawal: V solves r V = Lam beta H(V+k)
    p = params.copy(c=np.array([0.0]), f=np.array([1.0]), delta=0.0)
    sb = solve_sellers(8.0, p.k, p)
    assert abs(p.r * sb.V[0] - 8.0 * p.beta * float(p.G.H(sb.V[0] + p.k))) < 1e-10
    # frozen market: Lam -> 0 gives V -> (-c + delta w0)/(r+delta)
    sbf = solve_sellers(1e-9, params.k, params)
    assert np.allclose(sbf.V, (-params.c + params.delta * params.w0) / (params.r + params.delta), atol=1e-7)
    # very liquid market: v* -> w_H, prices -> w_H - k
    sbl = solve_sellers(1e7, params.k, params)
    assert np.all(sbl.vstar > params.G.hi - 1e-2)
    assert np.all(np.abs(sbl.pbar - (params.G.hi - params.k)) < 2e-2)
    # survival and hazards
    sb = solve_sellers(8.0, params.k, params)
    assert np.all(sb.eta >= 0) and np.all(sb.acc >= 0) and np.all(sb.acc <= 1)
    tau = np.linspace(0, 3, 50)
    surv = np.exp(-np.outer(tau, sb.eta + params.delta))
    assert np.all(surv <= 1) and np.all(np.diff(surv, axis=0) <= 0)


def test_lemma3_spurious_duration_dependence(params):
    sb = solve_sellers(8.0, params.k, params)
    mu = sb.eta + params.delta
    tau = np.linspace(0, 2, 400)
    S = np.exp(-np.outer(tau, mu)) @ params.f
    hazard = -np.gradient(np.log(S), tau)
    assert np.all(np.diff(hazard)[5:-5] < 0)