# we cannot import test_poisson_2indep directly, pytest treats that as test
from statsmodels.compat.python import PYTHON_IMPL_WASM

import warnings

import numpy as np
from numpy import arange
from numpy.testing import assert_allclose, assert_equal
import pytest
from scipy import stats

import statsmodels.stats.rates as smr
from statsmodels.stats.rates import (  # test_poisson, # cannot import functions that start with test
    confint_poisson,
    confint_poisson_2indep,
    confint_quantile_poisson,
    etest_poisson_2indep,
    method_names_poisson_1samp,
    method_names_poisson_2indep,
    nonequivalence_poisson_2indep,
    power_equivalence_neginb_2indep,
    power_equivalence_poisson_2indep,
    power_negbin_ratio_2indep,
    power_poisson_diff_2indep,
    power_poisson_ratio_2indep,
    tolerance_int_poisson,
)

methods = ["wald", "score", "exact-c", "waldccv", "sqrt-a", "sqrt-v", "midp-c",
           "sqrt",
           ]


@pytest.mark.parametrize("method", methods)
def test_rate_poisson_consistency(method):
    # check consistency between test and confint for one poisson rate
    count, nobs = 15, 400
    ci = confint_poisson(count, nobs, method=method)
    pv1 = smr.test_poisson(count, nobs, value=ci[0], method=method).pvalue
    pv2 = smr.test_poisson(count, nobs, value=ci[1], method=method).pvalue

    rtol = 1e-10
    if method in ["midp-c"]:
        # numerical root finding, lower precision
        rtol = 1e-6
    assert_allclose(pv1, 0.05, rtol=rtol)
    assert_allclose(pv2, 0.05, rtol=rtol)

    # check one-sided, note all confint are central
    pv1 = smr.test_poisson(count, nobs, value=ci[0], method=method,
                           alternative="larger").pvalue
    pv2 = smr.test_poisson(count, nobs, value=ci[1], method=method,
                           alternative="smaller").pvalue

    assert_allclose(pv1, 0.025, rtol=rtol)
    assert_allclose(pv2, 0.025, rtol=rtol)


def test_rate_poisson_r():
    count, nobs = 15, 400

    # DescTools
    # PoissonCI(x=15, n=400, method = c("exact","score", "wald", "byar"))
    # exact 0.0375 0.0209884653319583 0.0618505471787146
    # score 0.0375 0.0227264749053794 0.0618771721463559
    # wald  0.0375 0.0185227303217751 0.0564772696782249
    # byar  0.0375 0.0219084369245707 0.0602933510943048

    # > rr = poisson.exact(15, 400, r=0.05, tsmethod="central")
    # > rr$p.value
    # > rr$conf.int
    pv2 = 0.313026269279486
    ci2 = (0.0209884653319583, 0.0618505471787146)
    rt = smr.test_poisson(count, nobs, value=0.05, method="exact-c")
    ci = confint_poisson(count, nobs, method="exact-c")
    assert_allclose(rt.pvalue, pv2, rtol=1e-12)
    assert_allclose(ci, ci2, rtol=1e-12)

    # R package ratesci
    # sr = scoreci(15, 400, contrast="p", distrib = "poi", skew=FALSE,
    #      theta0=0.05)
    # pval computed from 2 * min(pval_left, pval_right)
    pv2 = 0.263552477282973
    # ci2 = (0.022726, 0.061877)  # no additional digits available
    ci2 = (0.0227264749053794, 0.0618771721463559)  # from DescTools
    rt = smr.test_poisson(count, nobs, value=0.05, method="score")
    ci = confint_poisson(count, nobs, method="score")
    assert_allclose(rt.pvalue, pv2, rtol=1e-12)
    assert_allclose(ci, ci2, rtol=1e-12)

    # > jeffreysci(15, 400, distrib = "poi")
    ci2 = (0.0219234232268444, 0.0602898619930649)
    ci = confint_poisson(count, nobs, method="jeff")
    assert_allclose(ci, ci2, rtol=1e-12)

    # from DescTools
    ci2 = (0.0185227303217751, 0.0564772696782249)
    ci = confint_poisson(count, nobs, method="wald")
    assert_allclose(ci, ci2, rtol=1e-12)

    # from ratesci
    # rateci(15, 400, distrib = "poi")
    # I think their lower ci is wrong, it also doesn't match for exact_c
    ci2 = (0.0243357599260795, 0.0604627555786095)
    ci = confint_poisson(count, nobs, method="midp-c")
    assert_allclose(ci[1], ci2[1], rtol=1e-5)


cases_tolint = [
    ("wald", 15, 1, 1, (3, 32), (3, np.inf), (0, 31)),
    ("score", 15, 1, 1, (4, 35), (4, np.inf), (0, 33)),
    ("wald", 150, 100, 100, (104, 200), (108, np.inf), (0, 196)),
    ("score", 150, 100, 100, (106, 202), (109, np.inf), (0, 198)),
    ("exact-c", 150, 100, 100, (105, 202), (109, np.inf), (0, 198)),
    ]


@pytest.mark.parametrize("case", cases_tolint)
def test_tol_int(case):
    # results values are from R package `tolerance, e.g.
    # poistol.int(15, 1, m = 1, alpha = 0.05, P = 0.975, side = 1, method="LS")
    # poistol.int(150, 100, m = 100, alpha = 0.05, P = 0.95, side = 2,
    #     method="SC")
    # exact-c is method="TAB"

    prob = 0.95
    prob_one = 0.975
    meth, count, exposure, exposure_new, r2, rs, rl = case
    ti = tolerance_int_poisson(
        count, exposure, prob, exposure_new=exposure_new,
        method=meth, alpha=0.05, alternative="two-sided")
    assert_equal(ti, r2)
    ti = tolerance_int_poisson(
        count, exposure, prob_one, exposure_new=exposure_new,
        method=meth, alpha=0.05, alternative="larger")
    assert_equal(ti, rl)
    ti = tolerance_int_poisson(
        count, exposure, prob_one, exposure_new=exposure_new,
        method=meth, alpha=0.05, alternative="smaller")
    assert_equal(ti, rs)

    # check without parameter uncertainty, confint at alpha=1
    if meth not in ["exact-c"]:
        # confint for "exact-c" does not collapse to point if alpha=1, check
        ti = tolerance_int_poisson(
            count, exposure, prob, exposure_new=exposure_new,
            method=meth, alpha=0.99999, alternative="two-sided")
        ci = stats.poisson.interval(prob, count / exposure * exposure_new)
        assert_equal(ti, ci)

    # check confint_quantile_poisson,
    # should same as upper tol_int
    ciq = confint_quantile_poisson(
        count, exposure, prob_one, exposure_new=exposure_new,
        method=meth, alpha=0.05, alternative="two-sided")

    assert_equal(ciq[1], r2[1])

    ciq = confint_quantile_poisson(
        count, exposure, prob_one, exposure_new=exposure_new,
        method=meth, alpha=0.05, alternative="larger")

    assert_equal(ciq[1], rl[1])

    # should same as lower tol_int
    prob_low = 0.025
    ciq = confint_quantile_poisson(
        count, exposure, prob_low, exposure_new=exposure_new,
        method=meth, alpha=0.05, alternative="two-sided")

    assert_equal(ciq[0], r2[0])

    ciq = confint_quantile_poisson(
        count, exposure, prob_low, exposure_new=exposure_new,
        method=meth, alpha=0.05, alternative="smaller")

    assert_equal(ciq[0], rs[0])


class TestMethodsCompar1samp():
    # check that all methods for test and confint produce similar results

    @pytest.mark.parametrize("meth", method_names_poisson_1samp["test"])
    def test_test(self, meth):
        count1, n1 = 60, 514.775
        tst = smr.test_poisson(count1, n1, method=meth, value=0.1,
                               alternative="two-sided")

        assert_allclose(tst.pvalue, 0.25, rtol=0.1)

    @pytest.mark.parametrize("meth", method_names_poisson_1samp["confint"])
    def test_confint(self, meth):
        count1, n1 = 60, 514.775
        ci = confint_poisson(count1, n1, method=meth, alpha=0.05)
        assert_allclose(ci, [0.089, 0.158], rtol=0.1)


# ######## below are 2indep tests


methods_diff = ["wald", "score", "waldccv",
                ]


@pytest.mark.parametrize("method", methods_diff)
def test_rate_poisson_diff_consistency(method):
    # check consistency between test and confint for diff of poisson rates
    count1, n1, count2, n2 = 30, 400 / 10, 7, 300 / 10   # avoid low=0
    ci = confint_poisson_2indep(count1, n1, count2, n2, method=method,
                                compare="diff")
    pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[0],
                                  method=method, compare="diff").pvalue
    pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[1],
                                  method=method, compare="diff").pvalue

    rtol = 1e-10
    if method in ["score"]:
        # numerical root finding, lower precision
        rtol = 1e-6
    assert_allclose(pv1, 0.05, rtol=rtol)
    assert_allclose(pv2, 0.05, rtol=rtol)

    # check one-sided, note all confint are central
    pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[0],
                                  method=method,
                                  compare="diff",
                                  alternative="larger").pvalue
    pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[1],
                                  method=method,
                                  compare="diff",
                                  alternative="smaller").pvalue

    assert_allclose(pv1, 0.025, rtol=rtol)
    assert_allclose(pv2, 0.025, rtol=rtol)


methods_ratio = ["wald-log", "score-log",  # "waldccv",
                 ]


@pytest.mark.parametrize("method", methods_ratio)
def test_rate_poisson_ratio_consistency(method):
    compare = "ratio"
    # check consistency between test and confint for ratio of poisson rates
    count1, n1, count2, n2 = 30, 400 / 10, 7, 300 / 10   # avoid low=0
    ci = confint_poisson_2indep(count1, n1, count2, n2, method=method,
                                compare=compare)
    pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[0],
                                  method=method, compare=compare).pvalue
    pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[1],
                                  method=method, compare=compare).pvalue

    rtol = 1e-10
    if method in ["score", "score-log"]:
        # numerical root finding, lower precision
        rtol = 1e-6
    assert_allclose(pv1, 0.05, rtol=rtol)
    assert_allclose(pv2, 0.05, rtol=rtol)

    # check one-sided, note all confint are central
    pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[0],
                                  method=method,
                                  compare=compare,
                                  alternative="larger").pvalue
    pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, value=ci[1],
                                  method=method,
                                  compare=compare,
                                  alternative="smaller").pvalue

    assert_allclose(pv1, 0.025, rtol=rtol)
    assert_allclose(pv2, 0.025, rtol=rtol)


methods_diff_ratio = ["wald", "score", "etest", "etest-wald",
                      ]


@pytest.mark.parametrize("method", methods_diff_ratio)
def test_rate_poisson_diff_ratio_consistency(method):
    # check consistency between test for rate and diff for equality null
    count1, n1, count2, n2 = 30, 400 / 10, 7, 300 / 10   # avoid low=0
    t1 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                 method=method, compare="ratio")
    t2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                 method=method, compare="diff")

    assert_allclose(t1.tuple, t2.tuple, rtol=1e-13)

    t1 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                 method=method, compare="ratio",
                                 alternative="larger")
    t2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                 method=method, compare="diff",
                                 alternative="larger")

    assert_allclose(t1.tuple, t2.tuple, rtol=1e-13)

    t1 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                 method=method, compare="ratio",
                                 alternative="smaller")
    t2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                 method=method, compare="diff",
                                 alternative="smaller")

    assert_allclose(t1.tuple, t2.tuple, rtol=1e-13)


def test_twosample_poisson():
    # testing against two examples in Gu et al

    # example 1
    count1, n1, count2, n2 = 60, 51477.5, 30, 54308.7

    s1, pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, method="wald")
    pv1r = 0.000356
    assert_allclose(pv1, pv1r*2, rtol=0, atol=5e-6)
    assert_allclose(s1, 3.384913, atol=0, rtol=5e-6)  # regression test

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="score")
    pv2r = 0.000316
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-6)
    assert_allclose(s2, 3.417402, atol=0, rtol=5e-6)  # regression test

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                      method="wald-log")
    pv2r = 0.000420
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-6)
    assert_allclose(s2, 3.3393, atol=0, rtol=5e-6)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                      method="score-log")
    pv2r = 0.000200
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-6)
    assert_allclose(s2, 3.5406, atol=0, rtol=5e-5)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="sqrt")
    pv2r = 0.000285
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-6)
    assert_allclose(s2, 3.445485, atol=0, rtol=5e-6)  # regression test

    # two-sided
    # example2
    # I don't know why it's only 2.5 decimal agreement, rounding?
    count1, n1, count2, n2 = 41, 28010, 15, 19017
    s1, pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, method="wald",
                                      value=1.5)
    pv1r = 0.2309
    assert_allclose(pv1, pv1r*2, rtol=0, atol=5e-4)
    assert_allclose(s1, 0.735447, atol=0, rtol=5e-6)  # regression test

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="score",
                                      value=1.5)
    pv2r = 0.2398
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-4)
    assert_allclose(s2, 0.706631, atol=0, rtol=5e-6)  # regression test

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                      method="wald-log", value=1.5)
    pv2r = 0.2402
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-4)
    assert_allclose(s2, 0.7056, atol=0, rtol=5e-4)

    with pytest.warns(FutureWarning):
        s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                          method="score-log", ratio_null=1.5)
    pv2r = 0.2303
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-4)
    assert_allclose(s2, 0.7380, atol=0, rtol=5e-4)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="sqrt",
                                      value=1.5)
    pv2r = 0.2499
    assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-3)
    assert_allclose(s2, 0.674401, atol=0, rtol=5e-6)  # regression test

    # one-sided
    # example 1 onesided
    count1, n1, count2, n2 = 60, 51477.5, 30, 54308.7

    s1, pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, method="wald",
                                      alternative="larger")
    pv1r = 0.000356
    assert_allclose(pv1, pv1r, rtol=0, atol=5e-6)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="score",
                                      alternative="larger")
    pv2r = 0.000316
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-6)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="sqrt",
                                      alternative="larger")
    pv2r = 0.000285
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-6)

    # 'exact-cond', 'cond-midp'

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                      method="exact-cond",
                                      value=1, alternative="larger")
    pv2r = 0.000428  # typo in Gu et al, switched pvalues between C and M
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                      method="cond-midp",
                                      value=1, alternative="larger")
    pv2r = 0.000310
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)

    _, pve1 = etest_poisson_2indep(count1, n1, count2, n2,
                                   method="score",
                                   alternative="larger")
    pve1r = 0.000298
    assert_allclose(pve1, pve1r, rtol=0, atol=5e-4)

    _, pve1 = etest_poisson_2indep(count1, n1, count2, n2,
                                   method="wald",
                                   alternative="larger")
    pve1r = 0.000298
    assert_allclose(pve1, pve1r, rtol=0, atol=5e-4)

    # example2 onesided
    # I don't know why it's only 2.5 decimal agreement, rounding?
    count1, n1, count2, n2 = 41, 28010, 15, 19017
    s1, pv1 = smr.test_poisson_2indep(count1, n1, count2, n2, method="wald",
                                      value=1.5, alternative="larger")
    pv1r = 0.2309
    assert_allclose(pv1, pv1r, rtol=0, atol=5e-4)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="score",
                                      value=1.5, alternative="larger")
    pv2r = 0.2398
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2, method="sqrt",
                                      value=1.5, alternative="larger")
    pv2r = 0.2499
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)

    # 'exact-cond', 'cond-midp'

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                      method="exact-cond",
                                      value=1.5, alternative="larger")
    pv2r = 0.2913
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)

    s2, pv2 = smr.test_poisson_2indep(count1, n1, count2, n2,
                                      method="cond-midp",
                                      value=1.5, alternative="larger")
    pv2r = 0.2450
    assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)

    _, pve2 = etest_poisson_2indep(count1, n1, count2, n2,
                                   method="score",
                                   value=1.5, alternative="larger")
    pve2r = 0.2453
    assert_allclose(pve2, pve2r, rtol=0, atol=5e-4)

    _, pve2 = etest_poisson_2indep(count1, n1, count2, n2,
                                   method="wald",
                                   value=1.5, alternative="larger")
    pve2r = 0.2453
    assert_allclose(pve2, pve2r, rtol=0, atol=5e-4)


cases_diff_ng = [
    ("wald",  (2.2047, 0.0137), (1.5514, 0.06040)),
    ("score", (2.0818, 0.0187), (1.5023, 0.06651)),
    ("etest-wald",  (2.2047, 0.0184), (1.5514, 0.06626)),
    # Ng reports 0.07088 at value=0.0002, looks strange
    ("etest-score", (2.0818, 0.0179), (1.5023, 0.06626)),  # 0.07088)),
    ]


@pytest.mark.parametrize("case", cases_diff_ng)
def test_twosample_poisson_diff(case):
    # testing against two examples Ng et al 2007
    # methods_diff = ["wald", "score", "etest-wald", "etest"]

    meth, res1, res2 = case
    count1, exposure1, count2, exposure2 = 41, 28010, 15, 19017

    value = 0
    t = smr.test_poisson_2indep(count1, exposure1, count2, exposure2,
                                value=value,
                                method=meth, compare="diff",
                                alternative="larger", etest_kwds=None)
    assert_allclose((t.statistic, t.pvalue), res1, atol=0.0006)

    value = 0.0002
    t = smr.test_poisson_2indep(count1, exposure1, count2, exposure2,
                                value=value,
                                method=meth, compare="diff",
                                alternative="larger", etest_kwds=None)
    assert_allclose((t.statistic, t.pvalue), res2, atol=0.0007)


def test_twosample_poisson_r():
    # testing against R package `exactci
    from .results.results_rates import res_pexact_cond, res_pexact_cond_midp

    # example 1 from Gu
    count1, n1, count2, n2 = 60, 51477.5, 30, 54308.7

    res2 = res_pexact_cond
    res1 = smr.test_poisson_2indep(count1, n1, count2, n2, method="exact-cond")
    assert_allclose(res1.pvalue, res2.p_value, rtol=1e-13)
    assert_allclose(res1.ratio, res2.estimate, rtol=1e-13)
    assert_equal(res1.ratio_null, res2.null_value)

    res2 = res_pexact_cond_midp
    res1 = smr.test_poisson_2indep(count1, n1, count2, n2, method="cond-midp")
    assert_allclose(res1.pvalue, res2.p_value, rtol=0, atol=5e-6)
    assert_allclose(res1.ratio, res2.estimate, rtol=1e-13)
    assert_equal(res1.ratio_null, res2.null_value)

    # one-sided
    # > pe = poisson.exact(c(60, 30), c(51477.5, 54308.7), r=1.2,
    #                      alternative="less", tsmethod="minlike", midp=TRUE)
    # > pe$p.value
    pv2 = 0.9949053964701466
    rest = smr.test_poisson_2indep(count1, n1, count2, n2, method="cond-midp",
                                   value=1.2, alternative="smaller")
    assert_allclose(rest.pvalue, pv2, rtol=1e-12)
    # > pe = poisson.exact(c(60, 30), c(51477.5, 54308.7), r=1.2,
    #           alternative="greater", tsmethod="minlike", midp=TRUE)
    # > pe$p.value
    pv2 = 0.005094603529853279
    rest = smr.test_poisson_2indep(count1, n1, count2, n2, method="cond-midp",
                                   value=1.2, alternative="larger")
    assert_allclose(rest.pvalue, pv2, rtol=1e-12)
    # > pe = poisson.exact(c(60, 30), c(51477.5, 54308.7), r=1.2,
    #           alternative="greater", tsmethod="minlike", midp=FALSE)
    # > pe$p.value
    pv2 = 0.006651774552714537
    rest = smr.test_poisson_2indep(count1, n1, count2, n2, method="exact-cond",
                                   value=1.2, alternative="larger")
    assert_allclose(rest.pvalue, pv2, rtol=1e-12)
    # > pe = poisson.exact(c(60, 30), c(51477.5, 54308.7), r=1.2,
    #                      alternative="less", tsmethod="minlike", midp=FALSE)
    # > pe$p.value
    pv2 = 0.9964625674930079
    rest = smr.test_poisson_2indep(count1, n1, count2, n2, method="exact-cond",
                                   value=1.2, alternative="smaller")
    assert_allclose(rest.pvalue, pv2, rtol=1e-12)


def test_confint_poisson_2indep():
    # test other confint methods

    count1, exposure1, count2, exposure2 = 60, 51477.5, 30, 54308.7

    ci = confint_poisson_2indep(count1, exposure1, count2, exposure2,
                                method="mover", compare="ratio", alpha=0.1,
                                method_mover="jeff",
                                )
    ci1 = (1.4667, 3.0608)  # LTW 2014
    assert_allclose(ci, ci1, atol=0.05)
    # moverci(60, 51477.5, 30, 54308.7, type = "jeff", distrib="poi",
    # contrast="RR", cc=0, level=0.9)
    #   Lower Estimate    Upper
    # est 1.466768 2.104135 3.058634
    # ratesci uses different center for mover-jeffreys
    ci1 = (1.466768, 3.058634)
    assert_allclose(ci, ci1, rtol=0.001)

    ci = confint_poisson_2indep(count1, exposure1, count2, exposure2,
                                method="mover", compare="ratio", alpha=0.1,
                                method_mover="score",
                                )
    ci1 = (1.4611, 3.0424)  # LTW 2014
    assert_allclose(ci, ci1, atol=0.05)

    ci = confint_poisson_2indep(count1, exposure1, count2, exposure2,
                                method="waldcc", compare="ratio", alpha=0.1,
                                )
    ci1 = (1.4523, 3.0154)  # LTW 2014
    assert_allclose(ci, ci1, atol=0.0005)

    ci = confint_poisson_2indep(count1, exposure1, count2, exposure2,
                                method="score", compare="ratio", alpha=0.05,
                                )
    ci1 = (1.365962, 3.259306)
    assert_allclose(ci, ci1, atol=5e-6)

    # with diff
    # use larger rates for diff
    exposure1 /= 1000
    exposure2 /= 1000
    ci = confint_poisson_2indep(count1, exposure1, count2, exposure2,
                                method="mover", compare="diff", alpha=0.05,
                                method_mover="jeff",
                                )
    # moverci(60, 51477.5/1000, 30, 54308.7/1000, type = "jeff", distrib="poi",
    # contrast="RD", cc=0, level=0.95)
    ci1 = (0.2629322, 0.9786493)
    assert_allclose(ci, ci1, atol=0.005)
    ci = confint_poisson_2indep(count1, exposure1, count2, exposure2,
                                method="score", compare="diff", alpha=0.05,
                                )
    # scoreci(60, 51477.5/1000, 30, 54308.7/1000, distrib="poi", contrast="RD",
    # skew=FALSE, cc=0, level=0.95)
    ci1 = (0.265796, 0.989192)
    assert_allclose(ci, ci1, atol=5e-6)

    # example in Li et al 2011, scaled
    ci = confint_poisson_2indep(count2, exposure2, count1, exposure1,
                                method="mover", compare="diff", alpha=0.1,
                                method_mover="jeff",
                                )
    # moverci(30, 54308.7/1000, 60, 51477.5/1000, type = "jeff", distrib="poi",
    # contrast="RD", cc=0, level=0.90)
    ci1 = (-0.9183272231752, -0.3188611692202)
    assert_allclose(ci, ci1, atol=0.005)
    ci1 = (-0.9195, -0.3193)  # from Li et al 2011
    assert_allclose(ci, ci1, atol=0.005)

    ci = confint_poisson_2indep(count2, exposure2, count1, exposure1,
                                method="mover", compare="diff", alpha=0.1,
                                method_mover="jeff",
                                )
    # scoreci(30, 54308.7/1000, 60, 51477.5/1000, distrib="poi", contrast="RD",
    # skew=FALSE, cc=0, level=0.90)
    ci1 = (-0.9232, -0.3188)  # from Li et al 2011
    assert_allclose(ci, ci1, atol=0.006)


def test_tost_poisson():

    count1, n1, count2, n2 = 60, 51477.5, 30, 54308.7
    # # central conf_int from R exactci
    low, upp = 1.339735721772650, 3.388365573616252

    res = smr.tost_poisson_2indep(count1, n1, count2, n2, low, upp,
                                  method="exact-cond")

    assert_allclose(res.pvalue, 0.025, rtol=1e-12)
    methods = ["wald", "score", "sqrt", "exact-cond", "cond-midp"]

    # test that we are in the correct range for other methods
    for meth in methods:
        res = smr.tost_poisson_2indep(count1, n1, count2, n2, low, upp,
                                      method=meth)
        assert_allclose(res.pvalue, 0.025, atol=0.01)


cases_alt = {
    ("two-sided", "wald"): 0.07136366497984171,
    ("two-sided", "score"): 0.0840167525117227,
    ("two-sided", "sqrt"): 0.0804675114297235,
    ("two-sided", "exact-cond"): 0.1301269270479679,
    ("two-sided", "cond-midp"): 0.09324590196774807,
    ("two-sided", "etest"): 0.09054824785458056,
    ("two-sided", "etest-wald"): 0.06895289560607239,

    ("larger", "wald"): 0.03568183248992086,
    ("larger", "score"): 0.04200837625586135,
    ("larger", "sqrt"): 0.04023375571486175,
    ("larger", "exact-cond"): 0.08570447732927276,
    ("larger", "cond-midp"): 0.04882345224905293,
    ("larger", "etest"): 0.043751060642682936,
    ("larger", "etest-wald"): 0.043751050280207024,

    ("smaller", "wald"): 0.9643181675100791,
    ("smaller", "score"): 0.9579916237441386,
    ("smaller", "sqrt"): 0.9597662442851382,
    ("smaller", "exact-cond"): 0.9880575728311669,
    ("smaller", "cond-midp"): 0.9511765477509471,
    ("smaller", "etest"): 0.9672396898656999,
    ("smaller", "etest-wald"): 0.9672397002281757
    }


@pytest.mark.parametrize("case", list(cases_alt.keys()))
def test_alternative(case):
    # regression test numbers, but those are close to each other
    alt, meth = case
    count1, n1, count2, n2 = 6, 51., 1, 54.
    _, pv = smr.test_poisson_2indep(count1, n1, count2, n2, method=meth,
                                    value=1.2, alternative=alt)
    assert_allclose(pv, cases_alt[case], rtol=1e-13)


class TestMethodsCompare2indep():
    # check that all methods for test and confint produce similar results

    @pytest.mark.parametrize(
        "compare, meth",
        [("ratio", meth) for meth
         in method_names_poisson_2indep["test"]["ratio"]] +
        [("diff", meth) for meth
         in method_names_poisson_2indep["test"]["diff"]]
        )
    def test_test(self, meth, compare):
        count1, n1, count2, n2 = 60, 514.775, 40, 543.0870000
        tst = smr.test_poisson_2indep(count1, n1, count2, n2, method=meth,
                                      compare=compare,
                                      value=None, alternative="two-sided")

        assert_allclose(tst.pvalue, 0.0245, rtol=0.2)

        # compare with degenerate interval for nonequivalence
        if compare == "ratio":
            f = 1.00
            low, upp = 1 / f, f
        else:
            v = 0.00
            low, upp = -v, v

        tst2 = nonequivalence_poisson_2indep(count1, n1, count2, n2, low, upp,
                                             method=meth,
                                             compare=compare)
        if "cond" in meth or "etest" in meth:
            # comparison is not exact for noncentral methods
            rtol = 0.1
        else:
            rtol = 1e-12
        assert_allclose(tst2.pvalue, tst.pvalue, rtol=rtol)

        # check corner case count2 = 0, see issue #8313
        if not PYTHON_IMPL_WASM:  # No fp exception support in WASM
            with pytest.warns(RuntimeWarning):
                smr.test_poisson_2indep(
                    count1, n1, 0, n2, method=meth,
                    compare=compare,
                    value=None, alternative="two-sided"
                )

    @pytest.mark.parametrize(
        "compare, meth",
        [("ratio", meth) for meth
         in method_names_poisson_2indep["confint"]["ratio"]] +
        [("diff", meth) for meth
         in method_names_poisson_2indep["confint"]["diff"]]
        )
    def test_confint(self, meth, compare):
        count1, n1, count2, n2 = 60, 514.775, 40, 543.0870000
        if compare == "ratio":
            ci_val = [1.04, 2.34]
        else:
            ci_val = [0.0057, 0.081]

        ci = confint_poisson_2indep(count1, n1, count2, n2, method=meth,
                                    compare=compare, alpha=0.05)
        assert_allclose(ci, ci_val, rtol=0.1)

    @pytest.mark.parametrize(
        "compare, meth",
        [("ratio", meth) for meth
         in method_names_poisson_2indep["test"]["ratio"]] +
        [("diff", meth) for meth
         in method_names_poisson_2indep["test"]["diff"]]
        )
    def test_test_vectorized(self, meth, compare):
        if "etest" in meth:
            pytest.skip("nonequivalence etest not vectorized")

        count1, n1, count2, n2 = 60, 514.775, 40, 543.0870000
        count1v = np.array([count1, count2])
        n1v = np.array([n1, n2])
        nfact = 1.
        # scipy binom test requires int dtype
        count2v = np.array([count2, count1 * nfact], dtype=int)
        n2v = np.array([n2, n1 * nfact])
        count1, n1, count2, n2 = count1v, n1v, count2v, n2v

        # compare with degenerate interval for nonequivalence
        if compare == "ratio":
            f = 1.00
            low, upp = 1 / f, f
        else:
            v = 0.00
            low, upp = -v, v

        tst2 = nonequivalence_poisson_2indep(count1, n1, count2, n2, low, upp,
                                             method=meth,
                                             compare=compare)
        assert tst2.statistic.shape == (2,)
        assert tst2.pvalue.shape == (2,)

        if not ("cond" in meth or "etest" in meth):
            # not all methods are vectorized for smr.test_poisson_2indep
            tst = smr.test_poisson_2indep(count1, n1, count2, n2, method=meth,
                                          compare=compare,
                                          value=None, alternative="two-sided")

            assert_allclose(tst2.pvalue, tst.pvalue, rtol=1e-12)

        # check vectorized equivalence test
        if compare == "ratio":
            f = 1.5
            low, upp = 1 / f, f
        else:
            v = 0.5
            low, upp = -v, v

        tst0 = smr.tost_poisson_2indep(count1[0], n1[0], count2[0], n2[0],
                                       low, upp,
                                       method=meth,
                                       compare=compare)
        tst1 = smr.tost_poisson_2indep(count1[1], n1[1], count2[1], n2[1],
                                       low, upp,
                                       method=meth,
                                       compare=compare)

        tst2 = smr.tost_poisson_2indep(count1, n1, count2, n2, low, upp,
                                       method=meth,
                                       compare=compare)
        assert tst2.statistic.shape == (2,)
        assert tst2.pvalue.shape == (2,)

        assert_allclose(tst2.statistic[0], tst0.statistic, rtol=1e-12)
        assert_allclose(tst2.pvalue[0], tst0.pvalue, rtol=1e-12)
        assert_allclose(tst2.statistic[1], tst1.statistic, rtol=1e-12)
        assert_allclose(tst2.pvalue[1], tst1.pvalue, rtol=1e-12)


