from itertools import product

import numpy as np
from numpy.testing import assert_allclose
import pandas as pd
import pytest

from arch.covariance.kernel import (
    Andrews,
    Bartlett,
    CovarianceEstimate,
    CovarianceEstimator,
    Gallant,
    NeweyWest,
    Parzen,
    ParzenCauchy,
    ParzenGeometric,
    ParzenRiesz,
    QuadraticSpectral,
    TukeyHamming,
    TukeyHanning,
    TukeyParzen,
)
from arch.typing import ArrayLike

ESTIMATORS = [
    Andrews,
    Bartlett,
    Gallant,
    NeweyWest,
    Parzen,
    ParzenCauchy,
    ParzenGeometric,
    ParzenRiesz,
    QuadraticSpectral,
    TukeyHamming,
    TukeyHanning,
    TukeyParzen,
]

KERNEL_PARAMS = {
    "Andrews": (1.3221, 2, 2 / 25),
    "Bartlett": (1.1447, 1, 2 / 9),
    "Gallant": (2.6614, 2, 4 / 25),
    "NeweyWest": (1.1447, 1, 2 / 9),
    "Parzen": (2.6614, 2, 4 / 25),
    "ParzenCauchy": (1.0924, 2, 4 / 25),
    "ParzenGeometric": (1.0000, 1, 2 / 9),
    "ParzenRiesz": (1.1340, 2, 4 / 25),
    "QuadraticSpectral": (1.3221, 2, 2 / 25),
    "TukeyHamming": (1.6694, 2, 4 / 25),
    "TukeyHanning": (1.7462, 2, 4 / 25),
    "TukeyParzen": (1.8576, 2, 4 / 25),
}


DATA_PARAMS = list(product([1, 2], [True, False]))
DATA_IDS = [f"{ndim}-{pandas}" for ndim, pandas in DATA_PARAMS]


@pytest.fixture(scope="module", params=DATA_PARAMS, ids=DATA_IDS)
def data(request) -> ArrayLike:
    rs = np.random.RandomState([3894830, 432841, 323297, 8927821])
    ndim, use_pandas = request.param
    burn = 100
    size: tuple[int, ...] = (500 + burn,)
    if ndim == 2:
        size += (3,)
    e = rs.standard_normal(size)
    if ndim == 1:
        phi = rs.uniform(0, 0.9)
        for i in range(1, size[0]):
            e[i] += e[i - 1] * phi
    else:
        phi_arr = np.diag(rs.uniform(0, 0.9, size[1]))
        for i in range(1, size[0]):
            e[i] += e[i - 1] @ phi_arr
    e = e[burn:]
    if use_pandas:
        if ndim == 1:
            return pd.Series(e, name="x")
        else:
            return pd.DataFrame(e, columns=[f"x{i}" for i in range(e.shape[1])])
    return e


@pytest.fixture(params=ESTIMATORS)
def estimator(request) -> CovarianceEstimator:
    return request.param


def test_covariance_smoke(data: ArrayLike, estimator: type[CovarianceEstimator]):
    cov = estimator(data)
    est_cov = cov.cov
    ndim = data.shape[1] if data.ndim > 1 else 1
    assert isinstance(est_cov, CovarianceEstimate)
    assert est_cov.long_run.shape == (ndim, ndim)
    assert est_cov.short_run.shape == (ndim, ndim)
    assert est_cov.one_sided_strict.shape == (ndim, ndim)
    assert est_cov.one_sided.shape == (ndim, ndim)
    assert isinstance(str(cov), str)
    assert isinstance(repr(cov), str)


def test_covariance_errors(data: ArrayLike, estimator: type[CovarianceEstimator]):
    with pytest.raises(ValueError, match="Degrees of freedom is <= 0"):
        estimator(data, df_adjust=data.shape[0] + 1)
    with pytest.raises(ValueError, match="df_adjust must be a non-negative"):
        estimator(data, df_adjust=-2)
    with pytest.raises(ValueError, match="df_adjust must be a non-negative"):
        # Type ignored due to invalid type used in test
        estimator(data, df_adjust=np.ones(2))  # type: ignore
    with pytest.raises(ValueError, match="bandwidth must be"):
        estimator(data, bandwidth=-3)
    with pytest.raises(ValueError, match="weights must be"):
        estimator(data, weights=np.ones(7))


def test_bartlett_auto(data: ArrayLike):
    nw = Bartlett(data, force_int=True)
    if data.ndim == 1:
        expected_bw = 11
    else:
        expected_bw = 12
    assert int(nw.bandwidth) == expected_bw
    expected = 1.0 - np.arange(nw.bandwidth + 1) / (nw.bandwidth + 1)
    assert_allclose(nw.kernel_weights, expected)
    resid = data - np.asarray(data.mean(0))
    resid = np.asarray(resid)
    nobs = resid.shape[0]
    expected_cov = resid.T @ resid / nobs
    assert_allclose(expected_cov, np.squeeze(nw.cov.short_run))
    expected_oss = np.zeros_like(expected_cov)
    for i in range(1, int(nw.bandwidth) + 1):
        gamma = (resid[i:].T @ resid[:-i]) / nobs
        expected_oss += expected[i] * gamma
        expected_cov += expected[i] * (gamma + gamma.T)
    assert_allclose(expected_oss, np.squeeze(nw.cov.one_sided_strict))
    assert_allclose(expected_cov, np.squeeze(nw.cov.long_run))
    ce = CovarianceEstimate(
        short_run=np.asarray(nw.cov.short_run),
        one_sided_strict=np.asarray(nw.cov.one_sided_strict),
        long_run=np.asarray(nw.cov.long_run),
        one_sided=np.asarray(nw.cov.one_sided),
    )
    assert_allclose(ce.short_run, nw.cov.short_run)
    assert_allclose(ce.one_sided_strict, nw.cov.one_sided_strict)
    assert_allclose(ce.long_run, nw.cov.long_run)
    assert_allclose(ce.one_sided, nw.cov.one_sided)


def test_parzen_auto(data: ArrayLike):
    pz = Parzen(data, force_int=True)

    if data.ndim == 1:
        # This test is noisy
        expected_bw: tuple[int, ...] = (18, 19)
        expected_weights = [
            1.0000e00,
            9.8575e-01,
            9.4600e-01,
            8.8525e-01,
            8.0800e-01,
            7.1875e-01,
            6.2200e-01,
            5.2225e-01,
            4.2400e-01,
            3.3175e-01,
            2.5000e-01,
            1.8225e-01,
            1.2800e-01,
            8.5750e-02,
            5.4000e-02,
            3.1250e-02,
            1.6000e-02,
            6.7500e-03,
            2.0000e-03,
            2.5000e-04,
        ]
    else:
        expected_bw = (17,)
        expected_weights = [
            1.00000000e00,
            9.82510288e-01,
            9.34156379e-01,
            8.61111111e-01,
            7.69547325e-01,
            6.65637860e-01,
            5.55555556e-01,
            4.45473251e-01,
            3.41563786e-01,
            2.50000000e-01,
            1.75582990e-01,
            1.17626886e-01,
            7.40740741e-02,
            4.28669410e-02,
            2.19478738e-02,
            9.25925926e-03,
            2.74348422e-03,
            3.42935528e-04,
        ]
    assert int(pz.bandwidth) in expected_bw
    if data.ndim > 1 or (pz.bandwidth + 1) == len(expected_weights):
        assert_allclose(pz.kernel_weights, np.array(expected_weights), rtol=1e-4)


def test_qs_auto(data: ArrayLike):
    qs = QuadraticSpectral(data, force_int=True)
    if data.ndim == 1:
        expected_bw = 8
        expected_weights = [
            1.0,
            0.97796879,
            0.91394558,
            0.81388998,
            0.68693073,
            0.54427698,
            0.3979104,
            0.25923469,
            0.13786058,
            0.04067932,
        ]
    else:
        expected_bw = 9
        expected_weights = [
            1.0,
            0.98256363,
            0.93155267,
            0.85073648,
            0.74599423,
            0.62475681,
            0.49531303,
            0.36605493,
            0.24474206,
            0.13786058,
        ]
    assert int(qs.bandwidth) == expected_bw
    assert_allclose(qs.kernel_weights[:10], np.array(expected_weights))


def test_force_int(data: ArrayLike, estimator: type[CovarianceEstimator]):
    bw = estimator(data, force_int=False).bandwidth
    bw_int = estimator(data, force_int=True).bandwidth
    assert bw_int >= bw
    assert bw_int == int(bw_int)


def test_first_weights(data: ArrayLike, estimator: type[CovarianceEstimator]):
    w = estimator(data).kernel_weights
    assert_allclose(w[0], 1.0)


def test_constants(data: ArrayLike, estimator: type[CovarianceEstimator]):
    cov_est = estimator(data)
    kc, bs, rate = KERNEL_PARAMS[estimator.__name__]
    assert_allclose(cov_est.kernel_const, kc)
    assert cov_est.bandwidth_scale == bs
    assert_allclose(cov_est.rate, rate)


def test_weight_len(data: ArrayLike, estimator: type[CovarianceEstimator]):
    cov_est = estimator(data, force_int=True)
    name = estimator.__name__
    is_qs = name in ("QuadraticSpectral", "Andrews")
    if is_qs:
        exp_len = data.shape[0]
    else:
        exp_len = int(cov_est.bandwidth) + 1
    assert cov_est.kernel_weights.shape[0] == exp_len


def test_kernel_weights(data: ArrayLike, estimator: type[CovarianceEstimator]):
    if data.ndim == 1:
        return
    weights = np.arange(1, data.shape[1] + 1)
    weights = weights / weights.mean()
    wcov = estimator(data, weights=weights, force_int=False)
    cov = estimator(data, force_int=False)
    assert wcov.bandwidth != cov.bandwidth


def test_center(data: ArrayLike, estimator: type[CovarianceEstimator]):
    centered_cov = estimator(data, center=False, force_int=False)
    cov = estimator(data, force_int=False)
    assert centered_cov.bandwidth != cov.bandwidth
