"""
    The file provides a class for logistic regression utilities with polynomial
    logit (log odds) function:

        p(x|theta) = 1/(1 + exp(-F(x|beta(theta))))

    with log odds F of polynomial form:

        F(x|beta) = sum_{j=0}^degree beta_j x^j

    with decision function (aka logit) F and parameters

        beta(theta) = [beta(theta)_j]_{j=0}^degree

    where theta are regression parameters and len(theta) = degree + 1.

    Coefficients beta(theta) can be constrained to be monotonic
    function of x by using monotonic cubic transformation:

        beta(theta) = mc.forward_map(theta)

    where mc is module mono_cubic2.

    NOTES: The model defines the conditional probability

        Prob(Y = y|x, theta) = 1/(1 + exp(-s(y) F(x| beta(theta))))

    where

        s(y) = 2*y - 1

    with x in R and y in {0,1}. For degree = 1 this is standard logistic regression

        beta(theta) = (theta[0], theta[1]).

    and generally without monotonicity condition

        beta(theta) = [theta_i]_{i=0}^degree

    For degree = 3 and mono = True the coefficients beta(theta) are constrained to be
    monotonic by using monotonic cubic transformation.

    We have data

        {(x_i, y_i) : x_i in R, y_i in {0,1}, i = 0, ..., n-1 }

    and the model is fit to data by minimizing negative log-likelihood function:

        nlff = -sum_i log(Prob(Y = y_i| x_i, theta))  : neg. log likelihood

    with respect to parameters theta. We can have regularization term in cost function
    and this case we minimize cost function:

        cost(theta) = nllf(theta) + lambda_0*|theta| + lambda_1*|theta|_2^2

    Author: Martin Horvat, January 2026
"""


import numpy as np
import scipy
import scipy.optimize

from . import monotonic as mc

def resize_with_const(v, n, val=0):
    """
    Resize a 1D vector to a specified length `n`.

    If the input vector `v` is longer than `n`, it is truncated.
    If it is shorter, it is padded with the constant value `val`.
    If it is already of length `n`, it is returned unchanged.

    Parameters:
        v (array-like): Input 1D vector (list or NumPy array).
        n (int): Target length of the output vector.
        val (scalar, optional): Value used to pad if `v` is shorter than `n`. Default is 0.

    Returns:
        np.ndarray: Resized 1D NumPy array of length `n`.
    """
    v = np.asarray(v)

    if len(v) == n: return v

    if len(v) > n: return v[:n]

    return np.concatenate([v, np.full(n - len(v), val)])

"""
    Fitting data

        {(x_i, y_i) in R x {0,1} : i = 0, ..., n-1}

    to model function

        f(x|theta) = 1/(1 + exp(-F(x|beta(theta))))

    with log odds of polynomial form:

        log(f(x|theta)/(1 - f(x|theta))) = F(x| beta(theta))

    where F is decision function (aka logit)

        F(x|beta) = sum_{i=0}^degree beta_i x^i

    and coefficients

        beta(theta) = [beta_i(theta)]_{i=0}^degrees

    Conditional probability

        Prob(Y = y_i|x, theta) = 1/(1 + exp(-s(y) F(x|beta(theta)))

    with

        s(y) = 2*y -1
"""
class LogisticPolyRegression:

    """
        Class constructor.

        Input:
            degree: int, degree of polynomial
            mono: boolean, default False
            lambda: None or tuple float, L2 regularization

    """
    def __init__(self, degree = 1, mono = False, lam = None):

        self.degree = degree
        self.mono = mono
        self.big = 1e3
        self.small = 1e-8
        self.lam = lam

        if mono and self.degree not in (1, 3):
            raise ValueError(
                f"Monotonic regression supports only degree 1 or 3, "
                f"not degree {self.degree}."
            )

        self.mono1 = self.mono and (self.degree == 1)
        self.mono3 = self.mono and (self.degree == 3)

    def _get_bounds(self):
        """Return optimizer bounds for the regression parameters."""
        bounds = [(-self.big, self.big)] * (self.degree + 1)

        # A positive linear coefficient makes both the logit and probability
        # strictly increasing in x.
        if self.mono1:
            bounds[1] = (self.small, self.big)

        return bounds

    def _validate_data(self, x, y):
        """Validate and return one-dimensional predictor and response arrays."""
        x = np.asarray(x, dtype=float)
        y = np.asarray(y)

        if x.ndim != 1 or y.ndim != 1:
            raise ValueError("x and y must be one-dimensional arrays.")
        if len(x) != len(y):
            raise ValueError(
                f"x and y must have the same length; got {len(x)} and {len(y)}."
            )
        if len(x) == 0:
            raise ValueError("x and y must not be empty.")
        if not np.all(np.isfinite(x)):
            raise ValueError("x must contain only finite values.")
        if not np.all(np.isfinite(y)):
            raise ValueError("y must contain only finite values.")
        if not np.all(np.isin(y, (0, 1))):
            raise ValueError("y must contain only the binary values 0 and 1.")

        return x, y.astype(int, copy=False)

    @staticmethod
    def _validate_sample_count(m):
        """Validate a requested number of generated parameter samples."""
        if isinstance(m, (bool, np.bool_)) or not isinstance(m, (int, np.integer)):
            raise TypeError("m must be a positive integer.")
        if m < 1:
            raise ValueError("m must be at least 1.")
        return int(m)

    def _validate_bootstrap_request(self, m, max_attempts):
        """Validate bootstrap size and return an explicit attempt limit."""
        m = self._validate_sample_count(m)

        if max_attempts is None:
            max_attempts = max(100, 10 * m)
        elif (
            isinstance(max_attempts, (bool, np.bool_))
            or not isinstance(max_attempts, (int, np.integer))
        ):
            raise TypeError("max_attempts must be an integer or None.")

        max_attempts = int(max_attempts)
        if max_attempts < m:
            raise ValueError("max_attempts must be at least m.")

        return m, max_attempts

    def _collect_bootstrap_theta(
        self,
        sampler,
        theta0,
        m,
        max_attempts,
    ):
        """Collect exactly ``m`` successful bootstrap parameter estimates."""
        samples = []
        attempts = 0

        while len(samples) < m and attempts < max_attempts:
            attempts += 1
            xb, yb = sampler()

            # A binary logistic-regression fit requires both classes.
            if np.unique(yb).size < 2:
                continue

            result = self.fit(
                xb,
                yb,
                theta0=theta0,
                method="local",
            )

            theta = np.asarray(result["theta"], dtype=float)
            if result["success"] and np.all(np.isfinite(theta)):
                samples.append(theta.copy())

        if len(samples) < m:
            raise RuntimeError(
                f"Generated only {len(samples)} successful bootstrap fits "
                f"from {attempts} attempts; requested {m}."
            )

        return np.stack(samples)

    def _fit_bootstrap_reference(self, x, y, seed):
        """Fit the original data before generating bootstrap samples."""
        if np.unique(y).size < 2:
            raise ValueError("Bootstrap data must contain both response classes.")

        result = self.fit(x, y, method="diff_evol", seed=seed)
        if not result["success"]:
            raise RuntimeError(
                f"Initial fit for bootstrap sampling failed: {result['message']}"
            )

        return np.asarray(result["theta"], dtype=float)

    """
        Mapping regression parameters theta to coefficients beta

            beta = beta(theta)

        used in decision function:

            F(x|beta) = sum_{i=0}^degree beta_i x^i

        Input:
            theta

        Return:
            beta
    """
    def get_beta(self, theta):

        beta = mc.forward_map(theta) if self.mono3 else theta
        return np.array(beta)

    """
        Mapping beta to regression parameters used in decision function.

        Input:
            beta: coefficient beta

        Return:
            theta
    """
    def get_theta(self, beta):

        theta = mc.backward_map(beta) if self.mono3 else beta
        return np.array(theta)

    """
        Calculate jacobian between beta and  regression parameters

            J = d(beta)/d(theta)
              = [d(beta_i)/d(theta_a)]_{i,a}
        and

            H = [d^2 beta_i/(d(theta_a) d(theta_b))]_{i,a,b}

        Input:
            theta
            hess: boolean, False

        Return:
            J       if hess = True
            (J, H)  if hess = False
    """
    def get_jac_beta(self, theta, hess = False):

        n = len(theta)

        J = mc.forward_map_jacobian(theta) if self.mono3 else np.eye(n)

        if not hess: return J

        H = mc.forward_map_hessian(theta) if self.mono3 else np.zeros(shape = (n, n, n))

        return (J, H)


    """
        Model function

            f(x|theta) = 1/(1 + exp(-F(x|beta)))  beta = beta(theta)

        Input:
            x: scalar value or a array of values
            theta: array of r = degree + 1 floats, model parameters array of floats

        Return:
            model function values
    """
    def model(self, x, theta):

        beta = self.get_beta(theta)  # beta
        X = np.column_stack([x**i for i in range(len(beta))])
        F = X @ beta              # decision function, X beta

        return scipy.special.expit(F)

    def get_x50(self, theta):
        """Return the first predictor value at which the model equals 0.5.

        Since ``expit(F) = 0.5`` exactly when ``F = 0``, ``x50`` is the
        smallest real root of the polynomial logit

            F(x) = sum_i beta_i(theta) x**i.

        For a monotonic fitted model the root is unique.  Defining "first" as
        the smallest real root also makes the result unambiguous for an
        unconstrained polynomial with several 0.5 crossings.

        Parameters
        ----------
        theta
            Model parameter vector of length ``degree + 1``.

        Returns
        -------
        float
            The smallest real solution of ``model(x, theta) = 0.5``.

        Raises
        ------
        ValueError
            If ``theta`` is invalid, the model never reaches 0.5 at a real
            predictor value, or the model is identically 0.5 and hence has no
            unique first crossing.
        """
        theta = np.asarray(theta, dtype=float)
        expected_shape = (self.degree + 1,)
        if theta.shape != expected_shape:
            raise ValueError(
                f"theta must have shape {expected_shape}; got {theta.shape}."
            )
        if not np.all(np.isfinite(theta)):
            raise ValueError("theta must contain only finite values.")

        beta = np.asarray(self.get_beta(theta), dtype=float)
        nonzero = np.flatnonzero(beta != 0.0)

        if nonzero.size == 0:
            raise ValueError(
                "x50 is not uniquely defined because the model equals 0.5 "
                "for every x."
            )

        polynomial_degree = int(nonzero[-1])
        if polynomial_degree == 0:
            raise ValueError("The model does not reach 0.5 for any real x.")

        coefficients = beta[:polynomial_degree + 1]

        # Avoid the unnecessary loss of precision of a general polynomial
        # root solver in the common linear-logistic case.
        if polynomial_degree == 1:
            return float(-coefficients[0] / coefficients[1])

        roots = np.roots(coefficients[::-1])
        real_roots = []
        machine_tolerance = 100 * np.finfo(float).eps

        for root in roots:
            candidate = float(root.real)
            root_scale = max(1.0, abs(candidate))

            # Repeated real roots may acquire a small imaginary part in a
            # numerical polynomial-root calculation.  In that case, also
            # accept the real component when its scaled polynomial residual
            # is negligible.
            residual = abs(
                np.polynomial.polynomial.polyval(candidate, coefficients)
            )
            coefficient_scale = np.polynomial.polynomial.polyval(
                abs(candidate), np.abs(coefficients)
            )
            small_imaginary_part = (
                abs(root.imag) <= machine_tolerance * root_scale
            )
            small_residual = residual <= 1e-10 * max(
                coefficient_scale, np.finfo(float).tiny
            )

            if small_imaginary_part or small_residual:
                real_roots.append(candidate)

        if not real_roots:
            raise ValueError("The model does not reach 0.5 for any real x.")

        return float(min(real_roots))

    def get_s50(self, theta):
        """Return the probability slope at the model's first 0.5 crossing.

        If ``p(x) = expit(F(x))``, then

            dp/dx = p(x) * (1 - p(x)) * F'(x).

        At ``x50``, ``p(x50) = 0.5``, so the reported midpoint slope is
        ``F'(x50) / 4``.  The derivative is with respect to the predictor on
        the scale supplied to the model.

        Parameters
        ----------
        theta
            Model parameter vector of length ``degree + 1``.

        Returns
        -------
        float
            ``d model(x, theta) / dx`` evaluated at ``x = get_x50(theta)``.
        """
        x50 = self.get_x50(theta)
        beta = np.asarray(self.get_beta(theta), dtype=float)
        derivative_coefficients = np.arange(1, len(beta)) * beta[1:]
        logit_slope = np.polynomial.polynomial.polyval(
            x50, derivative_coefficients
        )
        return float(logit_slope / 4.0)

    """
        Calculate negative log-likelihood function

            nllf = -sum_i log(Prob(Y = y_i| x_i, theta))   : neg. log likelihood

            grad = [d(nllf)/d(theta_a)]_a                  : jacobian

        where

            Prob(Y = y_i| x_i, theta) = 1/(1 + exp(-s_i F(x_i| beta)))  beta=beta(theta)
            s_i = 2*y_i  - 1

        Input:
            x: array of n floats
            y: array of n int in {0,1}
            theta: array of r = degree + 1 floats, model parameters array of floats
            jac: boolean, default False, if jacobian is needed

        Return:
            nllf           : if jac is false
            (nllf, grad)   : if jac is true
    """
    def get_nllf(self, x, y, theta, jac=False):
        beta = self.get_beta(theta)              # shape (p,)
        J = self.get_jac_beta(theta) if jac else None   # shape (p, q)

        X = np.column_stack([x**i for i in range(len(beta))])   # shape (n, p)
        s = 2.0 * y - 1.0

        eta = X @ beta
        z = s * eta

        # stable negative log-likelihood
        nllf = -np.sum(scipy.special.log_expit(z))

        if not jac:
            return nllf

        grad = -(s * (1.0 - scipy.special.expit(z))) @ (X @ J)
        return nllf, grad

    """
        Penalty function

         Input:
            theta: array of r = degree + 1 floats, model parameters array of floats
            jac: boolean, default False, if jacobian is needed

        Return:
            val           : if jac is false
            (val, grad)   : if jac is true
    """
    def penalty(self, theta, jac = False):

        val = self.lam[0]*np.sum(np.abs(theta)) + self.lam[1]*np.sum(theta**2)

        if jac:
            grad = self.lam[0]*np.sign(theta) + 2*self.lam[1]*theta
            return (val, grad)

        return val

    """
        Cost function
    """
    def get_cost(self, x, y, theta, jac = False):

        val = self.get_nllf(x, y, theta, jac)

        if self.lam is not None:
            pen = self.penalty(theta, jac)
            return (val[0] + pen[0], val[1] + pen[1]) if jac else val + pen

        return val

    """
        Estimate parameters.

        Input:
            x: array of n floats
            y: array of n int in {0,1}

        Return:
            theta0
    """
    def get_est_theta(self, x, y):

        L = np.log(2*len(x) + 1)

        if self.mono3:
            z = L*np.polyfit(x, 2.0*y - 1, 1)[::-1]
            beta = resize_with_const(z, self.degree + 1, 1e-8)
        else:
            beta = L*np.polyfit(x, 2.0*y - 1, self.degree)[::-1]

        return self.get_theta(beta)

    """
        Performing logistic regression with log odds of polynomial form:

            log(f(x|theta)/(1 - f(x|theta))) = F(x|beta)  beta = beta(theta)

        and this gives

            f(x|theta) = 1/(1 + exp(-F(x|beta)))

        where coefficient beta = [beta_i(theta)]_{i=0}^degree, with decision
        function (aka logit)

            F(x|beta) = sum_{i=0}^degree beta_i x^i

        Input:
            x: array of n floats
            y: array of n int in {0, 1}

        Return:
            dict {"theta", "cost", "success"}
    """
    def fit(self, x, y, theta0 = None, method = "local", seed = None):
        """Fit the model using a local or global optimization method.

        ``theta0`` is used as the initial point for local optimization.  If it is
        omitted, an initial estimate is calculated from the data.  ``seed`` is
        passed to stochastic global optimizers for reproducible fits.
        """
        x, y = self._validate_data(x, y)
        bnds = self._get_bounds()

        if method == "local":

            if theta0 is None:
                theta0 = self.get_est_theta(x, y)

            theta0 = np.asarray(theta0, dtype=float)
            expected_shape = (self.degree + 1,)
            if theta0.shape != expected_shape:
                raise ValueError(
                    f"theta0 must have shape {expected_shape}; got {theta0.shape}."
                )
            if not np.all(np.isfinite(theta0)):
                raise ValueError("theta0 must contain only finite values.")

            # L-BFGS-B requires its initial point to satisfy the bounds.
            lower = np.array([bound[0] for bound in bnds])
            upper = np.array([bound[1] for bound in bnds])
            theta0 = np.clip(theta0, lower, upper)

            cf = lambda theta: self.get_cost(x, y, theta, jac = True)
            res = scipy.optimize.minimize(
                cf,
                x0 = theta0,
                method = "L-BFGS-B",
                jac = True,
                bounds = bnds,
                tol = 1e-12,
            )

        elif method == "diff_evol":

            cf = lambda theta: self.get_cost(x, y, theta, jac = False)
            res_global = scipy.optimize.differential_evolution(
                cf,
                bounds = bnds,
                tol = 1e-8,
                polish = False,
                seed = seed,
            )

            cf = lambda theta: self.get_cost(x, y, theta, jac = True)
            res = scipy.optimize.minimize(
                cf,
                x0 = res_global.x,
                method = "L-BFGS-B",
                jac = True,
                bounds = bnds,
                tol = 1e-12,
            )

        elif method == "anneal":

            cf = lambda theta: self.get_cost(x, y, theta, jac = False)
            res = scipy.optimize.dual_annealing(
                cf,
                bounds = bnds,
                seed = seed,
            )

        else:
            raise ValueError(
                f"Unsupported fitting method {method!r}; expected "
                "'local', 'diff_evol', or 'anneal'."
            )

        return {
            "theta": res.x,
            "cost": res.fun,
            "success": res.success,
            "message": res.message,
            "nit": res.nit,
        }

    def _parametric_bootstrap_deviance(
        self,
        x,
        y,
        theta,
        m,
        seed,
        max_attempts,
    ):
        """Return observed deviance and its refitted bootstrap p-value."""
        m, max_attempts = self._validate_bootstrap_request(m, max_attempts)
        observed = 2.0 * self.get_nllf(x, y, theta)
        fitted_probabilities = self.model(x, theta)
        rng = np.random.default_rng(seed)

        simulated = []
        attempts = 0
        while len(simulated) < m and attempts < max_attempts:
            attempts += 1
            y_sim = rng.binomial(1, fitted_probabilities)

            # Degenerate samples do not support the fitted binary model.
            if np.unique(y_sim).size < 2:
                continue

            result = self.fit(
                x,
                y_sim,
                theta0=theta,
                method="local",
            )
            theta_sim = np.asarray(result["theta"], dtype=float)
            if not result["success"] or not np.all(np.isfinite(theta_sim)):
                continue

            simulated.append(2.0 * self.get_nllf(x, y_sim, theta_sim))

        if len(simulated) < m:
            raise RuntimeError(
                f"Generated only {len(simulated)} successful goodness-of-fit "
                f"bootstrap fits from {attempts} attempts; requested {m}."
            )

        simulated = np.asarray(simulated)
        p_value = (1 + np.count_nonzero(simulated >= observed)) / (m + 1)
        return observed, p_value

    def goodness_of_fit(
        self,
        x,
        y,
        theta,
        thresh = 0.5,
        regularization = False,
        bootstrap_samples = 1000,
        bootstrap_seed = 1977,
        bootstrap_max_attempts = None,
    ):
        """Calculate fit summaries and a bootstrap goodness-of-fit test.

        By default, AIC and BIC are based on the unpenalized log-likelihood,
        even when the model was fitted with a regularization penalty.  This
        keeps their objective common when comparing models fitted with
        different penalties.  Set ``regularization=True`` to include the
        configured penalty in the objective used for AIC and BIC.

        Goodness of fit is assessed using the unpenalized logistic deviance

            D = -2 sum_i [y_i log(p_i) + (1-y_i) log(1-p_i)].

        Its p-value is calibrated by a parametric bootstrap.  Each bootstrap
        response is sampled independently from Bernoulli(p_i) and the model is
        refitted with the same constraints and configured penalty before its
        deviance is calculated.  The returned p-value is

            (1 + number of simulated deviances >= observed deviance)
            / (bootstrap_samples + 1).

        This replaces the Pearson chi-square approximation, which is not valid
        when continuous predictors give approximately one Bernoulli observation
        per covariate pattern.  A large p-value means that the observed
        discrepancy is not unusual under the fitted model; it does not prove
        that the model is correct.
        """
        x, y = self._validate_data(x, y)
        theta = np.asarray(theta, dtype=float)

        # model probabilities
        p = self.model(x, theta)

        # Log likelihood, or the negative penalized objective when explicitly
        # requested.  With regularization=False this remains comparable across
        # models fitted using different penalty strengths.
        llf = -(
            self.get_cost(x, y, theta)
            if regularization
            else self.get_nllf(x, y, theta)
        )

        # information criteria
        k, n = len(theta), len(x)
        AIC = 2*k - 2*llf
        BIC = k*np.log(n) - 2*llf

        dof = n - k
        deviance, deviance_p_value = self._parametric_bootstrap_deviance(
            x,
            y,
            theta,
            m=bootstrap_samples,
            seed=bootstrap_seed,
            max_attempts=bootstrap_max_attempts,
        )

        # using model as classifier
        matches = y == np.heaviside(p - thresh, 1)

        # accuracy A
        A = np.count_nonzero(matches)/n

        return {"LLF": llf,
                "AIC": AIC,
                "BIC": BIC,
                "A" : A,
                "deviance": deviance,
                "p-value(deviance_bootstrap)": deviance_p_value,
                "deviance_bootstrap_samples": bootstrap_samples,
                "n": n, "k": k, "dof": dof}

    """
        Calculate hessian of cost function with respect to parameters theta

        Input:
            x: array of n floats
            y: array of n int in {0,1}
            theta: array of r = degree+1 floats, model parameters
        Return:
            H matrix of shape (r, r)
    """
    def get_cost_hessian(self, x, y, theta):

        # coefficients
        beta = self.get_beta(theta)

        # design matrix -- add column of 1's at the beginning of your X_train matrix
        X = np.column_stack([x**i for i in range(len(beta))])

        # Jacobian J = [dbeta_i/dtheta_a]_{ia}
        # Hessian H = [d^2 beta_i/(d(theta_a) d(theta_b))]_{i,a,b}
        J, H = self.get_jac_beta(theta, hess = True)

        # signs
        s = 2.0*y - 1

        # decision function for conditional probability Prob(Y = y| x)
        V = s[:,None]*X
        F = V @ beta

        # probabilities p_i = P(Y=y_i | x_i)
        p = scipy.special.expit(F)
        q = 1 - p

        # calculate hessian
        L = V @ J
        H = (L.T*(q*p))@L - np.tensordot(q@V, H, axes = ([0], [0]))

        if self.lam is not None:
            return H + 2*self.lam[1]*np.eye(len(theta))   # H' = H + lambda id

        return H

    def get_cov(self, x, y, theta, method="model_sandwich"):
        """Estimate the covariance matrix of fitted parameters.

        Parameters
        ----------
        x, y
            Predictor and binary response arrays.
        theta
            Fitted regression parameters.
        method
            ``"model_sandwich"`` (default) uses

                A^+ [sum_i p_i(1-p_i) g_i g_i.T] A^+,

            where ``A`` is the observed Hessian of the penalized objective and
            ``g_i = d eta_i / d theta``.  This accounts for the fact that the
            deterministic penalty changes the Hessian but does not contribute
            sampling variability.

            ``"robust_sandwich"`` replaces the middle matrix by the empirical
            outer product of unpenalized scores,

                sum_i [(y_i-p_i) g_i][(y_i-p_i) g_i].T.

            ``"inverse_hessian"`` returns the previous approximation ``A^+``
            for comparison.  It should not be treated as the frequentist
            covariance of a penalized estimator.

        Notes
        -----
        ``+`` denotes a symmetric Moore-Penrose pseudoinverse.  These are local
        normal approximations.  For the monotonic cubic model, they can be
        unreliable near ``epsilon=0``, where the estimate is on the boundary
        and the parameter map loses rank.  Bootstrap inference is preferred in
        that case.  A nonzero L1 penalty is not supported because its Hessian is
        undefined at zero.
        """
        x, y = self._validate_data(x, y)
        theta = np.asarray(theta, dtype=float)

        valid_methods = {
            "model_sandwich",
            "robust_sandwich",
            "inverse_hessian",
        }
        if method not in valid_methods:
            raise ValueError(
                f"Unknown covariance method {method!r}; expected one of "
                f"{sorted(valid_methods)}."
            )

        if self.lam is not None and self.lam[0] != 0:
            raise ValueError(
                "Covariance estimation does not support a nonzero L1 penalty."
            )

        # Bread: exact observed Hessian of the smooth penalized objective.
        bread = self.get_cost_hessian(x, y, theta)
        bread = 0.5 * (bread + bread.T)
        bread_inv = np.linalg.pinv(bread, hermitian=True)

        if method == "inverse_hessian":
            return bread_inv

        beta = self.get_beta(theta)
        X = np.column_stack([x**i for i in range(len(beta))])
        J = self.get_jac_beta(theta)
        gradient_eta = X @ J
        probabilities = self.model(x, theta)

        if method == "model_sandwich":
            weights = probabilities * (1.0 - probabilities)
            meat = (gradient_eta.T * weights) @ gradient_eta
        else:
            scores = (y - probabilities)[:, None] * gradient_eta
            meat = scores.T @ scores

        cov = bread_inv @ meat @ bread_inv
        return 0.5 * (cov + cov.T)

    """
        Calculating standard errors fo model parameters for normal distribution of parameters:

            theta ~ N(mean_theta, cov_theta)

        Input:

            cov_theta: array of rxr floats, variance-covariance matrix of parameters

        Return:
            array of r floats, standard errors of parameters
    """
    def get_SE_theta_normal(self, cov_theta):

        # computing standard errors of parameters
        return np.sqrt(np.clip(np.diag(cov_theta),a_min=0, a_max = None))

    """
        Calculating quantiles of the model parameters at given probabilities p
        for normal distribution of parameters:

            theta ~ N(mean_theta, cov_theta)

        Input:
            probs: array of m floats, probabilities
            mean_theta: array of r = degree+1 floats, mean model parameters
            cov_theta: array of rxr floats, variance-covariance matrix of parameters

        Return:
            array of mxr floats
    """
    def get_theta_quantiles_normal(self, probs, mean_theta, cov_theta):

        # mean and standard variance parameters
        locs = mean_theta
        scales = np.sqrt(np.clip(np.diag(cov_theta), a_min=0, a_max=None))

        # computing quantiles of parameters
        return locs + np.outer(scipy.stats.norm.ppf(probs), scales)

    """
        Calculating quantiles of the model values

            f(x|theta) = 1/(1 + exp(-F(x|beta)))     beta = beta(theta)

        with

            F(x|beta) = sum_{i=0}^degree x^i beta_i

        at given probabilities p and values x assuming
        normal distribution of parameters:

            theta ~ N(mean_theta, cov_theta)

        This distribution is asymptotic MLE distribution of parameters.

        Input:
            x: array of n float
            probs: array of m floats, probabilities
            mean_theta: array of r = degree+1 floats, mean model parameters
            cov_theta: array of rxr floats, variance-covariance matrix of parameters

        Return:
            array of mxn floats
    """

    def get_model_quantiles_normal(self, x, probs, mean_theta, cov_theta,
                                   exact = True, seed = 1977, m = 10**5):

        mean_beta = self.get_beta(mean_theta)
        X = np.column_stack([x**i for i in range(len(mean_beta))])

        if exact and self.mono3:

            # init random generator
            rng = np.random.default_rng(seed)
            theta = rng.multivariate_normal(mean_theta, cov_theta, size = m)

            # get betas
            beta = np.apply_along_axis(self.get_beta, 1, theta)

            # quantiles of decision function
            Q = np.quantile(X@beta.T, probs, axis = 1)

            return np.apply_along_axis(scipy.special.expit, 1, Q)

        # J = d(beta)/d(theta)
        J = self.get_jac_beta(mean_theta)

        # transform data
        S = X@J

        # mean and standard variance of logit (aka log of odds)
        locs = X@mean_beta
        scales = np.sqrt(np.diag(S@cov_theta@S.T))

        # computing quantiles of logit
        Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)

        # convert logit to expit
        return scipy.special.expit(Q)

    """
        Calculating quantiles using delta method of the model values

            f(x|theta) = 1/(1 + exp(-F(x|beta))) beta = beta(theta)

        with

            F(x|beta) = sum_{i=0}^degree x^i beta_i

        at given probabilities p and values x assuming
        normal distribution of parameters :

            theta ~ N(mean_theta, cov_theta)

        This distribution is asymptotic MLE distribution of parameters.
        We approximate exact model with linear expansion

            f(x|theta) = p(x|theta_mean) + dp/db(x|theta_mean) (theta - theta_mean)

        and the last term is normally distributed.

        Input:
            x: array of n float
            probs: array of m floats, probabilities
            mean_theta: array of r = degree+1 floats, mean model parameters
            cov_theta: array of rxr floats, variance-covariance matrix of parameters

        Return:
            array of mxn floats
    """
    def get_model_quantiles_delta(self, x, probs, mean_theta, cov_theta):

        mean_beta = self.get_beta(mean_theta)
        X = np.column_stack([x**i for i in range(len(mean_beta))])
        F = X@mean_beta

        # J = d(beta)/d(theta)
        J = self.get_jac_beta(mean_theta)

        # S = d(F)/d(theta)
        S = X@J

        # attributes of normal distribution of model values
        locs = scipy.special.expit(F)

        derivative = locs * (1 - locs)
        scales = np.sqrt(np.diag(S @ cov_theta @ S.T)) * derivative

        # computing quantiles of logit
        Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)

        return np.clip(Q, a_min = 0, a_max = 1)

    """
        Generate parameters assuming normal distribution.

        Input:
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator

        Return:
            array of mx(degree + 1)
    """
    def get_normal_theta(self, x, y, m, seed = 1977):
        x, y = self._validate_data(x, y)
        m = self._validate_sample_count(m)

        res = self.fit(x, y, method = "diff_evol", seed = seed)

        if not res["success"]:
            raise RuntimeError(f"Fit did not succeed: {res['message']}")

        # optimal parameters
        theta = res["theta"]

        # covariance matrix of parameters
        cov = self.get_cov(x, y, theta)

        # init random generator
        rng = np.random.default_rng(seed)

        return rng.multivariate_normal(theta, cov, size = m)

    """
        Generate m parameters via non-parametric bootstrapping with a minimal constraint
        that both groups should be present in the sampled data:

        boostrapped sample = (xb, yb)  by sampling with replacement pairs (x_i, y_i)
                                    with condition that yb can not be just 0 or just 1

        Input:
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator

        Return:
            array of mx(degree + 1)
    """
    def get_nonparam_boots_theta(
        self,
        x,
        y,
        m,
        seed = 1977,
        max_attempts = None,
    ):
        x, y = self._validate_data(x, y)
        m, max_attempts = self._validate_bootstrap_request(m, max_attempts)
        theta0 = self._fit_bootstrap_reference(x, y, seed)

        n = len(x)
        rng = np.random.default_rng(seed)

        def sampler():
            idx = rng.choice(n, n, replace=True)
            return x[idx], y[idx]

        return self._collect_bootstrap_theta(
            sampler,
            theta0,
            m,
            max_attempts,
        )


    """
        Generate m parameters via non-parametric stratified bootstrapping:

            boostrapped sample = (xb, yb)

                xb = (sampled with replacement from x0, sampled with replacement from x1)
                yb = (0...0, 1...1)

        Input:
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator

        Return:
            array of mx(degree + 1)
    """
    def get_nonparam_stratified_boots_theta(
        self,
        x,
        y,
        m,
        seed = 1977,
        max_attempts = None,
    ):
        x, y = self._validate_data(x, y)
        m, max_attempts = self._validate_bootstrap_request(m, max_attempts)
        theta0 = self._fit_bootstrap_reference(x, y, seed)

        xs = [x[y == i] for i in range(2)]
        ns = [len(group) for group in xs]
        yb = np.concatenate(
            [np.full(ns[i], i, dtype=int) for i in range(2)]
        )

        rng = np.random.default_rng(seed)

        def sampler():
            xb = np.concatenate(
                [rng.choice(xs[i], ns[i], replace=True) for i in range(2)]
            )
            return xb, yb

        return self._collect_bootstrap_theta(
            sampler,
            theta0,
            m,
            max_attempts,
        )

    """
        Generate m parameters via parametric bootstrapping:

            boostrapped sample = (x, yb) yb ~ B(model(x, fitted theta))

        where B is Bernoulli distribution

        Input:
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator

        Return:
            array of mx(degree + 1)

        Ref:

        * https://www.scirp.org/journal/paperinformation?paperid=70962
        * https://en.wikipedia.org/wiki/Bernoulli_distribution
    """
    def get_parametric_boots_theta(
        self,
        x,
        y,
        m,
        seed = 1977,
        max_attempts = None,
    ):
        x, y = self._validate_data(x, y)
        m, max_attempts = self._validate_bootstrap_request(m, max_attempts)
        theta0 = self._fit_bootstrap_reference(x, y, seed)

        p = self.model(x, theta0)
        rng = np.random.default_rng(seed)

        def sampler():
            y_sim = rng.binomial(n = 1, p = p)
            return x, y_sim

        return self._collect_bootstrap_theta(
            sampler,
            theta0,
            m,
            max_attempts,
        )