"""
Tests for discrete models

Notes
-----
DECIMAL_3 is used because it seems that there is a loss of precision
in the Stata *.dta -> *.csv output, NOT the estimator for the Poisson
tests.
"""

# pylint: disable-msg=E1101
from statsmodels.compat.pandas import assert_index_equal

import os
import warnings

import numpy as np
from numpy.testing import (
    assert_,
    assert_allclose,
    assert_almost_equal,
    assert_array_equal,
    assert_array_less,
    assert_equal,
)
import pandas as pd
import pytest
from scipy import stats
from scipy.stats import nbinom

import statsmodels.api as sm
from statsmodels.discrete.discrete_margins import _iscount, _isdummy
from statsmodels.discrete.discrete_model import (
    CountModel,
    GeneralizedPoisson,
    Logit,
    MNLogit,
    NegativeBinomial,
    NegativeBinomialP,
    Poisson,
    Probit,
)
import statsmodels.formula.api as smf
from statsmodels.tools.sm_exceptions import (
    ConvergenceWarning,
    PerfectSeparationError,
    SpecificationWarning,
    ValueWarning,
)

from .results.results_discrete import Anes, DiscreteL1, RandHIE, Spector

try:
    import cvxopt  # noqa:F401

    has_cvxopt = True
except ImportError:
    has_cvxopt = False


DECIMAL_14 = 14
DECIMAL_10 = 10
DECIMAL_9 = 9
DECIMAL_4 = 4
DECIMAL_3 = 3
DECIMAL_2 = 2
DECIMAL_1 = 1
DECIMAL_0 = 0


def load_anes96():
    data = sm.datasets.anes96.load()
    data.endog = np.asarray(data.endog)
    data.exog = np.asarray(data.exog)
    return data


def load_spector():
    data = sm.datasets.spector.load()
    data.endog = np.asarray(data.endog)
    data.exog = np.asarray(data.exog)
    return data


def load_randhie():
    data = sm.datasets.randhie.load()
    data.endog = np.asarray(data.endog)
    data.exog = np.asarray(data.exog, dtype=float)
    return data


def check_jac(self, res=None):
    # moved from CheckModelResults
    if res is None:
        res1 = self.res1
    else:
        res1 = res

    exog = res1.model.exog
    # basic cross check
    jacsum = res1.model.score_obs(res1.params).sum(0)
    score = res1.model.score(res1.params)
    assert_almost_equal(jacsum, score, DECIMAL_9)  # Poisson has low precision ?

    if isinstance(res1.model, (NegativeBinomial, MNLogit)):
        # skip the rest
        return

    # check score_factor
    # TODO: change when score_obs uses score_factor for DRYing
    s1 = res1.model.score_obs(res1.params)
    sf = res1.model.score_factor(res1.params)
    if not isinstance(sf, tuple):
        s2 = sf[:, None] * exog
    else:
        sf0, sf1 = sf
        s2 = np.column_stack((sf0[:, None] * exog, sf1))

    assert_allclose(s2, s1, rtol=1e-10)

    # check hessian_factor
    h1 = res1.model.hessian(res1.params)
    hf = res1.model.hessian_factor(res1.params)
    if not isinstance(hf, tuple):
        h2 = (hf * exog.T).dot(exog)
    else:
        hf0, hf1, hf2 = hf
        h00 = (hf0 * exog.T).dot(exog)
        h10 = np.atleast_2d(hf1.T.dot(exog))
        h11 = np.atleast_2d(hf2.sum(0))
        h2 = np.vstack((np.column_stack((h00, h10.T)), np.column_stack((h10, h11))))

    assert_allclose(h2, h1, rtol=1e-10)


def check_distr(res):
    distr = res.get_distribution()
    distr1 = res.model.get_distribution(res.params)
    m = res.predict()
    m2 = distr.mean()
    assert_allclose(m, np.squeeze(m2), rtol=1e-10)
    m2 = distr1.mean()
    assert_allclose(m, np.squeeze(m2), rtol=1e-10)

    v = res.predict(which="var")
    v2 = distr.var()
    assert_allclose(v, np.squeeze(v2), rtol=1e-10)


class CheckModelMixin:
    # Assertions about the Model object, as opposed to the Results
    # Assumes that mixed-in class implements:
    #   res1

    def test_fit_regularized_invalid_method(self):
        # GH#5224 check we get ValueError when passing invalid "method" arg
        model = self.res1.model

        with pytest.raises(ValueError, match=r"is not supported, use either"):
            model.fit_regularized(method="foo")


class CheckModelResults(CheckModelMixin):
    """
    res2 should be the test results from RModelWrap
    or the results as defined in model_results_data
    """

    def test_params(self):
        assert_almost_equal(self.res1.params, self.res2.params, DECIMAL_4)

    def test_conf_int(self):
        assert_allclose(self.res1.conf_int(), self.res2.conf_int, rtol=8e-5)

    def test_zstat(self):
        assert_almost_equal(self.res1.tvalues, self.res2.z, DECIMAL_4)

    def test_pvalues(self):
        assert_almost_equal(self.res1.pvalues, self.res2.pvalues, DECIMAL_4)

    def test_cov_params(self):
        if not hasattr(self.res2, "cov_params"):
            raise pytest.skip("TODO: implement res2.cov_params")
        assert_almost_equal(self.res1.cov_params(), self.res2.cov_params, DECIMAL_4)

    def test_llf(self):
        assert_almost_equal(self.res1.llf, self.res2.llf, DECIMAL_4)

    def test_llnull(self):
        assert_almost_equal(self.res1.llnull, self.res2.llnull, DECIMAL_4)

    def test_llr(self):
        assert_almost_equal(self.res1.llr, self.res2.llr, DECIMAL_3)

    def test_llr_pvalue(self):
        assert_almost_equal(self.res1.llr_pvalue, self.res2.llr_pvalue, DECIMAL_4)

    @pytest.mark.xfail(
        reason="Test has not been implemented for this class.",
        strict=True,
        raises=NotImplementedError,
    )
    def test_normalized_cov_params(self):
        raise NotImplementedError

    def test_bse(self):
        assert_almost_equal(self.res1.bse, self.res2.bse, DECIMAL_4)

    def test_dof(self):
        assert_equal(self.res1.df_model, self.res2.df_model)
        assert_equal(self.res1.df_resid, self.res2.df_resid)

    def test_aic(self):
        assert_almost_equal(self.res1.aic, self.res2.aic, DECIMAL_3)

    def test_bic(self):
        assert_almost_equal(self.res1.bic, self.res2.bic, DECIMAL_3)

    def test_predict(self):
        assert_almost_equal(
            self.res1.model.predict(self.res1.params), self.res2.phat, DECIMAL_4
        )

    def test_predict_xb(self):
        assert_almost_equal(
            self.res1.model.predict(self.res1.params, which="linear"),
            self.res2.yhat,
            DECIMAL_4,
        )

    def test_loglikeobs(self):
        # basic cross check
        llobssum = self.res1.model.loglikeobs(self.res1.params).sum()
        assert_almost_equal(llobssum, self.res1.llf, DECIMAL_14)

    def test_jac(self):
        check_jac(self)

    def test_summary_latex(self):
        # see #7747, last line of top table was dropped
        summ = self.res1.summary()
        ltx = summ.as_latex()
        n_lines = len(ltx.splitlines())
        if not isinstance(self.res1.model, MNLogit):
            # skip MNLogit which creates several params tables
            assert n_lines == 19 + np.size(self.res1.params)
        assert "Covariance Type:" in ltx

    def test_distr(self):
        check_distr(self.res1)


class CheckBinaryResults(CheckModelResults):
    def test_pred_table(self):
        assert_array_equal(self.res1.pred_table(), self.res2.pred_table)

    def test_resid_dev(self):
        assert_almost_equal(self.res1.resid_dev, self.res2.resid_dev, DECIMAL_4)

    def test_resid_generalized(self):
        assert_almost_equal(
            self.res1.resid_generalized, self.res2.resid_generalized, DECIMAL_4
        )

    @pytest.mark.smoke
    def test_resid_response(self):
        assert isinstance(self.res1.resid_response, np.ndarray)


class CheckMargEff:
    """
    Test marginal effects (margeff) and its options
    """

    def test_nodummy_dydxoverall(self):
        me = self.res1.get_margeff()
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dydx, DECIMAL_4)
        assert_almost_equal(me.margeff_se, self.res2.margeff_nodummy_dydx_se, DECIMAL_4)

        me_frame = me.summary_frame()
        eff = me_frame["dy/dx"].values
        assert_allclose(eff, me.margeff, rtol=1e-13)
        assert_equal(me_frame.shape, (me.margeff.size, 6))

    def test_nodummy_dydxmean(self):
        me = self.res1.get_margeff(at="mean")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dydxmean, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_dydxmean_se, DECIMAL_4
        )

    def test_nodummy_dydxmedian(self):
        me = self.res1.get_margeff(at="median")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dydxmedian, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_dydxmedian_se, DECIMAL_4
        )

    def test_nodummy_dydxzero(self):
        me = self.res1.get_margeff(at="zero")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dydxzero, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_dydxzero, DECIMAL_4
        )

    def test_nodummy_dyexoverall(self):
        me = self.res1.get_margeff(method="dyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dyex, DECIMAL_4)
        assert_almost_equal(me.margeff_se, self.res2.margeff_nodummy_dyex_se, DECIMAL_4)

    def test_nodummy_dyexmean(self):
        me = self.res1.get_margeff(at="mean", method="dyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dyexmean, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_dyexmean_se, DECIMAL_4
        )

    def test_nodummy_dyexmedian(self):
        me = self.res1.get_margeff(at="median", method="dyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dyexmedian, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_dyexmedian_se, DECIMAL_4
        )

    def test_nodummy_dyexzero(self):
        me = self.res1.get_margeff(at="zero", method="dyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_dyexzero, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_dyexzero_se, DECIMAL_4
        )

    def test_nodummy_eydxoverall(self):
        me = self.res1.get_margeff(method="eydx")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eydx, DECIMAL_4)
        assert_almost_equal(me.margeff_se, self.res2.margeff_nodummy_eydx_se, DECIMAL_4)

    def test_nodummy_eydxmean(self):
        me = self.res1.get_margeff(at="mean", method="eydx")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eydxmean, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_eydxmean_se, DECIMAL_4
        )

    def test_nodummy_eydxmedian(self):
        me = self.res1.get_margeff(at="median", method="eydx")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eydxmedian, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_eydxmedian_se, DECIMAL_4
        )

    def test_nodummy_eydxzero(self):
        me = self.res1.get_margeff(at="zero", method="eydx")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eydxzero, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_eydxzero_se, DECIMAL_4
        )

    def test_nodummy_eyexoverall(self):
        me = self.res1.get_margeff(method="eyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eyex, DECIMAL_4)
        assert_almost_equal(me.margeff_se, self.res2.margeff_nodummy_eyex_se, DECIMAL_4)

    def test_nodummy_eyexmean(self):
        me = self.res1.get_margeff(at="mean", method="eyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eyexmean, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_eyexmean_se, DECIMAL_4
        )

    def test_nodummy_eyexmedian(self):
        me = self.res1.get_margeff(at="median", method="eyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eyexmedian, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_eyexmedian_se, DECIMAL_4
        )

    def test_nodummy_eyexzero(self):
        me = self.res1.get_margeff(at="zero", method="eyex")
        assert_almost_equal(me.margeff, self.res2.margeff_nodummy_eyexzero, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_nodummy_eyexzero_se, DECIMAL_4
        )

    def test_dummy_dydxoverall(self):
        me = self.res1.get_margeff(dummy=True)
        assert_almost_equal(me.margeff, self.res2.margeff_dummy_dydx, DECIMAL_4)
        assert_almost_equal(me.margeff_se, self.res2.margeff_dummy_dydx_se, DECIMAL_4)

    def test_dummy_dydxmean(self):
        me = self.res1.get_margeff(at="mean", dummy=True)
        assert_almost_equal(me.margeff, self.res2.margeff_dummy_dydxmean, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_dummy_dydxmean_se, DECIMAL_4
        )

    def test_dummy_eydxoverall(self):
        me = self.res1.get_margeff(method="eydx", dummy=True)
        assert_almost_equal(me.margeff, self.res2.margeff_dummy_eydx, DECIMAL_4)
        assert_almost_equal(me.margeff_se, self.res2.margeff_dummy_eydx_se, DECIMAL_4)

    def test_dummy_eydxmean(self):
        me = self.res1.get_margeff(at="mean", method="eydx", dummy=True)
        assert_almost_equal(me.margeff, self.res2.margeff_dummy_eydxmean, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_dummy_eydxmean_se, DECIMAL_4
        )

    def test_count_dydxoverall(self):
        me = self.res1.get_margeff(count=True)
        assert_almost_equal(me.margeff, self.res2.margeff_count_dydx, DECIMAL_4)
        assert_almost_equal(me.margeff_se, self.res2.margeff_count_dydx_se, DECIMAL_4)

    def test_count_dydxmean(self):
        me = self.res1.get_margeff(count=True, at="mean")
        assert_almost_equal(me.margeff, self.res2.margeff_count_dydxmean, DECIMAL_4)
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_count_dydxmean_se, DECIMAL_4
        )

    def test_count_dummy_dydxoverall(self):
        me = self.res1.get_margeff(count=True, dummy=True)
        assert_almost_equal(
            me.margeff, self.res2.margeff_count_dummy_dydxoverall, DECIMAL_4
        )
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_count_dummy_dydxoverall_se, DECIMAL_4
        )

    def test_count_dummy_dydxmean(self):
        me = self.res1.get_margeff(count=True, dummy=True, at="mean")
        assert_almost_equal(
            me.margeff, self.res2.margeff_count_dummy_dydxmean, DECIMAL_4
        )
        assert_almost_equal(
            me.margeff_se, self.res2.margeff_count_dummy_dydxmean_se, DECIMAL_4
        )


class TestProbitNewton(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        cls.res1 = Probit(data.endog, data.exog).fit(method="newton", disp=0)
        res2 = Spector.probit
        cls.res2 = res2

    def test_init_kwargs(self):
        endog = self.res1.model.endog
        exog = self.res1.model.exog
        z = np.ones(len(endog))
        with pytest.warns(ValueWarning, match="unknown kwargs"):
            # unsupported keyword
            Probit(endog, exog, weights=z)


class TestProbitBFGS(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        cls.res1 = Probit(data.endog, data.exog).fit(method="bfgs", disp=0)
        res2 = Spector.probit
        cls.res2 = res2


class TestProbitNM(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2
        cls.res1 = Probit(data.endog, data.exog).fit(method="nm", disp=0, maxiter=500)


class TestProbitPowell(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2
        cls.res1 = Probit(data.endog, data.exog).fit(method="powell", disp=0, ftol=1e-8)


class TestProbitCG(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2

        # fmin_cg fails to converge on some machines - reparameterize
        from statsmodels.tools.transform_model import StandardizeTransform

        transf = StandardizeTransform(data.exog)
        exog_st = transf(data.exog)
        res1_st = Probit(data.endog, exog_st).fit(
            method="cg", disp=0, maxiter=1000, gtol=1e-08
        )
        start_params = transf.transform_params(res1_st.params)
        assert_allclose(start_params, res2.params, rtol=1e-5, atol=1e-6)

        cls.res1 = Probit(data.endog, data.exog).fit(
            start_params=start_params, method="cg", maxiter=1000, gtol=1e-05, disp=0
        )

        assert_array_less(cls.res1.mle_retvals["fcalls"], 100)


class TestProbitNCG(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2
        cls.res1 = Probit(data.endog, data.exog).fit(
            method="ncg", disp=0, avextol=1e-8, warn_convergence=False
        )
        # converges close enough but warnflag is 2 for precision loss


class TestProbitBasinhopping(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2
        fit = Probit(data.endog, data.exog).fit
        np.random.seed(1)
        cls.res1 = fit(
            method="basinhopping",
            disp=0,
            niter=5,
            minimizer={"method": "L-BFGS-B", "tol": 1e-8},
        )


class TestProbitMinimizeDefault(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2
        fit = Probit(data.endog, data.exog).fit
        cls.res1 = fit(method="minimize", disp=0, niter=5, tol=1e-8)


class TestProbitMinimizeDogleg(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2
        fit = Probit(data.endog, data.exog).fit
        cls.res1 = fit(
            method="minimize", disp=0, niter=5, tol=1e-8, min_method="dogleg"
        )


class TestProbitMinimizeAdditionalOptions(CheckBinaryResults):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=False)
        res2 = Spector.probit
        cls.res2 = res2
        cls.res1 = Probit(data.endog, data.exog).fit(
            method="minimize",
            disp=0,
            maxiter=500,
            min_method="Nelder-Mead",
            xatol=1e-4,
            fatol=1e-4,
        )


class CheckLikelihoodModelL1:
    """
    For testing results generated with L1 regularization
    """

    def test_params(self):
        assert_almost_equal(self.res1.params, self.res2.params, DECIMAL_4)

    def test_conf_int(self):
        assert_almost_equal(self.res1.conf_int(), self.res2.conf_int, DECIMAL_4)

    def test_bse(self):
        assert_almost_equal(self.res1.bse, self.res2.bse, DECIMAL_4)

    def test_nnz_params(self):
        assert_almost_equal(self.res1.nnz_params, self.res2.nnz_params, DECIMAL_4)

    def test_aic(self):
        assert_almost_equal(self.res1.aic, self.res2.aic, DECIMAL_3)

    def test_bic(self):
        assert_almost_equal(self.res1.bic, self.res2.bic, DECIMAL_3)


class TestProbitL1(CheckLikelihoodModelL1):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=True)
        alpha = np.array([0.1, 0.2, 0.3, 10])  # / data.exog.shape[0]
        cls.res1 = Probit(data.endog, data.exog).fit_regularized(
            method="l1",
            alpha=alpha,
            disp=0,
            trim_mode="auto",
            auto_trim_tol=0.02,
            acc=1e-10,
            maxiter=1000,
        )
        res2 = DiscreteL1.probit
        cls.res2 = res2

    def test_cov_params(self):
        assert_almost_equal(self.res1.cov_params(), self.res2.cov_params, DECIMAL_4)


class TestMNLogitL1(CheckLikelihoodModelL1):

    @classmethod
    def setup_class(cls):
        anes_data = load_anes96()
        anes_exog = anes_data.exog
        anes_exog = sm.add_constant(anes_exog, prepend=False)
        mlogit_mod = sm.MNLogit(anes_data.endog, anes_exog)
        alpha = 10.0 * np.ones((mlogit_mod.J - 1, mlogit_mod.K))  # / anes_exog.shape[0]
        alpha[-1, :] = 0
        cls.res1 = mlogit_mod.fit_regularized(
            method="l1",
            alpha=alpha,
            trim_mode="auto",
            auto_trim_tol=0.02,
            acc=1e-10,
            disp=0,
        )
        res2 = DiscreteL1.mnlogit
        cls.res2 = res2


class TestLogitL1(CheckLikelihoodModelL1):

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=True)
        cls.alpha = 3 * np.array([0.0, 1.0, 1.0, 1.0])  # / data.exog.shape[0]
        cls.res1 = Logit(data.endog, data.exog).fit_regularized(
            method="l1",
            alpha=cls.alpha,
            disp=0,
            trim_mode="size",
            size_trim_tol=1e-5,
            acc=1e-10,
            maxiter=1000,
        )
        res2 = DiscreteL1.logit
        cls.res2 = res2

    def test_cov_params(self):
        assert_almost_equal(self.res1.cov_params(), self.res2.cov_params, DECIMAL_4)


@pytest.mark.skipif(
    not has_cvxopt, reason="Skipped test_cvxopt since cvxopt is not available"
)
class TestCVXOPT:

    @classmethod
    def setup_class(cls):
        if not has_cvxopt:
            pytest.skip("Skipped test_cvxopt since cvxopt is not available")
        cls.data = sm.datasets.spector.load()
        cls.data.endog = np.asarray(cls.data.endog)
        cls.data.exog = np.asarray(cls.data.exog)
        cls.data.exog = sm.add_constant(cls.data.exog, prepend=True)

    def test_cvxopt_versus_slsqp(self):
        # Compares results from cvxopt to the standard slsqp
        self.alpha = 3.0 * np.array([0, 1, 1, 1.0])  # / self.data.endog.shape[0]
        res_slsqp = Logit(self.data.endog, self.data.exog).fit_regularized(
            method="l1",
            alpha=self.alpha,
            disp=0,
            acc=1e-10,
            maxiter=1000,
            trim_mode="auto",
        )
        res_cvxopt = Logit(self.data.endog, self.data.exog).fit_regularized(
            method="l1_cvxopt_cp",
            alpha=self.alpha,
            disp=0,
            abstol=1e-10,
            trim_mode="auto",
            auto_trim_tol=0.01,
            maxiter=1000,
        )
        assert_almost_equal(res_slsqp.params, res_cvxopt.params, DECIMAL_4)


class TestSweepAlphaL1:

    @classmethod
    def setup_class(cls):
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=True)
        cls.model = Logit(data.endog, data.exog)
        cls.alphas = np.array(
            [[0.1, 0.1, 0.1, 0.1], [0.4, 0.4, 0.5, 0.5], [0.5, 0.5, 1, 1]]
        )  # / data.exog.shape[0]
        cls.res1 = DiscreteL1.sweep

    def test_sweep_alpha(self):
        for i in range(3):
            alpha = self.alphas[i, :]
            res2 = self.model.fit_regularized(
                method="l1",
                alpha=alpha,
                disp=0,
                acc=1e-10,
                trim_mode="off",
                maxiter=1000,
            )
            assert_almost_equal(res2.params, self.res1.params[i], DECIMAL_4)


class CheckL1Compatability:
    """
    Tests compatability between l1 and unregularized by setting alpha such
    that certain parameters should be effectively unregularized, and others
    should be ignored by the model.
    """

    def test_params(self):
        m = self.m
        assert_almost_equal(
            self.res_unreg.params[:m], self.res_reg.params[:m], DECIMAL_4
        )
        # The last entry should be close to zero
        # handle extra parameter of NegativeBinomial
        kvars = self.res_reg.model.exog.shape[1]
        assert_almost_equal(0, self.res_reg.params[m:kvars], DECIMAL_4)

    def test_cov_params(self):
        m = self.m
        # The restricted cov_params should be equal
        assert_almost_equal(
            self.res_unreg.cov_params()[:m, :m],
            self.res_reg.cov_params()[:m, :m],
            DECIMAL_1,
        )

    def test_df(self):
        assert_equal(self.res_unreg.df_model, self.res_reg.df_model)
        assert_equal(self.res_unreg.df_resid, self.res_reg.df_resid)

    def test_t_test(self):
        m = self.m
        kvars = self.kvars
        # handle extra parameter of NegativeBinomial
        extra = getattr(self, "k_extra", 0)
        t_unreg = self.res_unreg.t_test(np.eye(len(self.res_unreg.params)))
        t_reg = self.res_reg.t_test(np.eye(kvars + extra))
        assert_almost_equal(t_unreg.effect[:m], t_reg.effect[:m], DECIMAL_3)
        assert_almost_equal(t_unreg.sd[:m], t_reg.sd[:m], DECIMAL_3)
        assert_almost_equal(np.nan, t_reg.sd[m])
        assert_allclose(t_unreg.tvalue[:m], t_reg.tvalue[:m], atol=3e-3)
        assert_almost_equal(np.nan, t_reg.tvalue[m])

    def test_f_test(self):
        m = self.m
        kvars = self.kvars
        # handle extra parameter of NegativeBinomial
        extra = getattr(self, "k_extra", 0)
        f_unreg = self.res_unreg.f_test(np.eye(len(self.res_unreg.params))[:m])
        f_reg = self.res_reg.f_test(np.eye(kvars + extra)[:m])
        assert_allclose(f_unreg.fvalue, f_reg.fvalue, rtol=3e-5, atol=1e-3)
        assert_almost_equal(f_unreg.pvalue, f_reg.pvalue, DECIMAL_3)

    def test_bad_r_matrix(self):
        kvars = self.kvars
        with pytest.raises(ValueError):
            self.res_reg.f_test(np.eye(kvars))


class TestPoissonL1Compatability(CheckL1Compatability):

    @classmethod
    def setup_class(cls):
        cls.kvars = 10  # Number of variables
        cls.m = 7  # Number of unregularized parameters
        rand_data = load_randhie()
        rand_exog = rand_data.exog.view(float).reshape(len(rand_data.exog), -1)
        rand_exog = sm.add_constant(rand_exog, prepend=True)
        # Drop some columns and do an unregularized fit
        exog_no_PSI = rand_exog[:, : cls.m]
        mod_unreg = sm.Poisson(rand_data.endog, exog_no_PSI)
        cls.res_unreg = mod_unreg.fit(method="newton", disp=False)
        # Do a regularized fit with alpha, effectively dropping the last column
        alpha = 10 * len(rand_data.endog) * np.ones(cls.kvars)
        alpha[: cls.m] = 0
        cls.res_reg = sm.Poisson(rand_data.endog, rand_exog).fit_regularized(
            method="l1",
            alpha=alpha,
            disp=False,
            acc=1e-10,
            maxiter=2000,
            trim_mode="auto",
        )


class TestNegativeBinomialL1Compatability(CheckL1Compatability):

    @classmethod
    def setup_class(cls):
        cls.kvars = 10  # Number of variables
        cls.m = 7  # Number of unregularized parameters
        rand_data = load_randhie()
        rand_exog = rand_data.exog.view(float).reshape(len(rand_data.exog), -1)
        rand_exog_st = (rand_exog - rand_exog.mean(0)) / rand_exog.std(0)
        rand_exog = sm.add_constant(rand_exog_st, prepend=True)
        # Drop some columns and do an unregularized fit
        exog_no_PSI = rand_exog[:, : cls.m]
        mod_unreg = sm.NegativeBinomial(rand_data.endog, exog_no_PSI)
        cls.res_unreg = mod_unreg.fit(method="newton", disp=False)
        # Do a regularized fit with alpha, effectively dropping the last column
        alpha = 10 * len(rand_data.endog) * np.ones(cls.kvars + 1)
        alpha[: cls.m] = 0
        alpha[-1] = 0  # do not penalize alpha

        mod_reg = sm.NegativeBinomial(rand_data.endog, rand_exog)
        cls.res_reg = mod_reg.fit_regularized(
            method="l1",
            alpha=alpha,
            disp=False,
            acc=1e-10,
            maxiter=2000,
            trim_mode="auto",
        )
        cls.k_extra = 1  # 1 extra parameter in nb2


class TestNegativeBinomialGeoL1Compatability(CheckL1Compatability):

    @classmethod
    def setup_class(cls):
        cls.kvars = 10  # Number of variables
        cls.m = 7  # Number of unregularized parameters
        rand_data = load_randhie()
        rand_exog = rand_data.exog.view(float).reshape(len(rand_data.exog), -1)
        rand_exog = sm.add_constant(rand_exog, prepend=True)
        # Drop some columns and do an unregularized fit
        exog_no_PSI = rand_exog[:, : cls.m]
        mod_unreg = sm.NegativeBinomial(
            rand_data.endog, exog_no_PSI, loglike_method="geometric"
        )
        cls.res_unreg = mod_unreg.fit(method="newton", disp=False)
        # Do a regularized fit with alpha, effectively dropping the last columns
        alpha = 10 * len(rand_data.endog) * np.ones(cls.kvars)
        alpha[: cls.m] = 0
        mod_reg = sm.NegativeBinomial(
            rand_data.endog, rand_exog, loglike_method="geometric"
        )
        cls.res_reg = mod_reg.fit_regularized(
            method="l1",
            alpha=alpha,
            disp=False,
            acc=1e-10,
            maxiter=2000,
            trim_mode="auto",
        )

        assert_equal(mod_reg.loglike_method, "geometric")


class TestLogitL1Compatability(CheckL1Compatability):

    @classmethod
    def setup_class(cls):
        cls.kvars = 4  # Number of variables
        cls.m = 3  # Number of unregularized parameters
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=True)
        # Do a regularized fit with alpha, effectively dropping the last column
        alpha = np.array([0, 0, 0, 10])
        cls.res_reg = Logit(data.endog, data.exog).fit_regularized(
            method="l1", alpha=alpha, disp=0, acc=1e-15, maxiter=2000, trim_mode="auto"
        )
        # Actually drop the last columnand do an unregularized fit
        exog_no_PSI = data.exog[:, : cls.m]
        cls.res_unreg = Logit(data.endog, exog_no_PSI).fit(disp=0, tol=1e-15)


class TestMNLogitL1Compatability(CheckL1Compatability):

    @classmethod
    def setup_class(cls):
        cls.kvars = 4  # Number of variables
        cls.m = 3  # Number of unregularized parameters
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=True)
        alpha = np.array([0, 0, 0, 10])
        cls.res_reg = MNLogit(data.endog, data.exog).fit_regularized(
            method="l1", alpha=alpha, disp=0, acc=1e-15, maxiter=2000, trim_mode="auto"
        )
        # Actually drop the last columnand do an unregularized fit
        exog_no_PSI = data.exog[:, : cls.m]
        cls.res_unreg = MNLogit(data.endog, exog_no_PSI).fit(
            disp=0, gtol=1e-15, method="bfgs", maxiter=1000
        )

    def test_t_test(self):
        m = self.m
        kvars = self.kvars
        t_unreg = self.res_unreg.t_test(np.eye(m))
        t_reg = self.res_reg.t_test(np.eye(kvars))
        assert_almost_equal(t_unreg.effect, t_reg.effect[:m], DECIMAL_3)
        assert_almost_equal(t_unreg.sd, t_reg.sd[:m], DECIMAL_3)
        assert_almost_equal(np.nan, t_reg.sd[m])
        assert_almost_equal(t_unreg.tvalue, t_reg.tvalue[:m], DECIMAL_3)

    @pytest.mark.skip("Skipped test_f_test for MNLogit")
    def test_f_test(self):
        pass


class TestProbitL1Compatability(CheckL1Compatability):

    @classmethod
    def setup_class(cls):
        cls.kvars = 4  # Number of variables
        cls.m = 3  # Number of unregularized parameters
        data = load_spector()
        data.exog = sm.add_constant(data.exog, prepend=True)
        alpha = np.array([0, 0, 0, 10])
        cls.res_reg = Probit(data.endog, data.exog).fit_regularized(
            method="l1", alpha=alpha, disp=0, acc=1e-15, maxiter=2000, trim_mode="auto"
        )
        # Actually drop the last columnand do an unregularized fit
        exog_no_PSI = data.exog[:, : cls.m]
        cls.res_unreg = Probit(data.endog, exog_no_PSI).fit(disp=0, tol=1e-15)


class CompareL1:
    """
    For checking results for l1 regularization.
