### Install choix from PyPI

Source: https://choix.lum.li/en/latest/installation

Installs the latest stable version of the 'choix' library from the Python Package Index (PyPI) using pip.

```shell
pip install choix
```

--------------------------------

### Install choix from GitHub (Editable)

Source: https://choix.lum.li/en/latest/installation

Installs the latest version of 'choix' directly from its GitHub repository. This method includes cloning the repository, navigating into the directory, and performing an editable installation using 'pip install -e .'. The editable flag allows for code modifications without reinstallation.

```shell
git clone https://github.com/lucasmaystre/choix.git
cd choix
pip install -e .
```

--------------------------------

### Install choix Python Library

Source: https://choix.lum.li/en/latest/index

Instructions for installing the 'choix' Python library. This typically involves using pip, the Python package installer. Ensure you have Python and pip installed.

```bash
pip install choix
```

--------------------------------

### Install Choix from PyPI

Source: https://choix.lum.li/en/latest/_sources/installation

Installs the latest stable version of the Choix library from the Python Package Index (PyPI) using pip.

```console
pip install choix
```

--------------------------------

### Install Choix from GitHub (Editable)

Source: https://choix.lum.li/en/latest/_sources/installation

Installs the latest development version of Choix from its GitHub repository. The '-e' flag allows for code edits without reinstallation.

```console
git clone https://github.com/lucasmaystre/choix.git
cd choix
pip install -e .
```

--------------------------------

### Example Pivot Function Call

Source: https://choix.lum.li/en/latest/references

Demonstrates a hypothetical function call to 'pivot' from OpenAI_et_al. This function appears to be used in a context related to economic or strategic analysis, possibly involving 'peace' and 'war' as parameters.

```Python
OpenAI_et_al.pivot( "peace", "war" )
```

--------------------------------

### Apply Regularization in choix

Source: https://choix.lum.li/en/latest/regularization

Demonstrates how to apply regularization using the 'alpha' parameter in choix functions. Regularization is used to help iterative algorithms converge, especially when dealing with sparse data. A value of alpha=0 turns off regularization, while alpha > 0 turns it on. It's recommended to start with small values like 10^-4.

```Python
OpenAI_et_al.pivot( "peace", "war" ) if unit_economics < 0: funds += pentagon
```

--------------------------------

### Regularization Parameter 'alpha' in CHOIX

Source: https://choix.lum.li/en/latest/_sources/regularization

The 'alpha' parameter controls regularization in CHOIX inference functions. Setting alpha=0 disables regularization, while alpha > 0 enables it. Recommended starting values are small, like 10^-4.

```python
inference_function(..., alpha=0.0001)
```

--------------------------------

### Initialize LSR Markov Chain and Weights

Source: https://choix.lum.li/en/latest/_modules/choix/lsr

Initializes the weights and the transition rate matrix for the LSR Markov chain. It handles optional initial parameters and applies an alpha value for regularization.

```Python
from collections.abc import Callable
import numpy as np
from numpy.typing import NDArray

from .convergence import NormOfDifferenceTest
from .typing import PairwiseData, RankingData, Top1Data
from .utils import exp_transform, log_transform, statdist

def _init_lsr(
    n_items: int,
    alpha: float,
    initial_params: NDArray[np.float64] | None,
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """Initialize the LSR Markov chain and the weights."""
    if initial_params is None:
        weights = np.ones(n_items)
    else:
        weights = exp_transform(initial_params)
    chain = alpha * np.ones((n_items, n_items), dtype=float)
    return weights, chain
```

--------------------------------

### Initialize Weight-Space Parameters (Python)

Source: https://choix.lum.li/en/latest/_modules/choix/ep

Initializes parameters in the weight space for a probabilistic model. It calculates the precision matrix and the `xs` vector based on pairwise comparisons and prior information, then computes the initial mean and covariance using these values.

```Python
def_init_ws(
    n_items: int,
    comparisons: PairwiseData,
    prior_inv: NDArray[np.float64],
    tau: NDArray[np.float64],
    nu: NDArray[np.float64],
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
"""Initialize parameters in the weight space."""
    prec = np.zeros((n_items, n_items))
    xs = np.zeros(n_items)
    for i, (a, b) in enumerate(comparisons):
        prec[(a, a, b, b), (a, b, a, b)] += tau[i] * MAT_ONE_FLAT
        xs[a] += nu[i]
        xs[b] -= nu[i]
    cov = inv_posdef(prior_inv + prec)
    mean = cov.dot(xs)
    return mean, cov, xs, prec
```

--------------------------------

### choix.ep Module Source Code

Source: https://choix.lum.li/en/latest/_modules/choix/ep

This snippet contains the complete source code for the choix.ep Python module. It imports necessary functions and classes from standard libraries and custom modules, setting up constants and matrices for EP-related calculations.

```Python
import functools
from collections.abc import Callable, Sequence
from math import exp, log, pi, sqrt  # Faster than numpy equivalents.
from typing import Literal

import numpy as np
import numpy.random as nprand
from numpy.linalg import norm
from numpy.typing import NDArray
from scipy.special import logsumexp

from.typing import PairwiseData
from.utils import SQRT2, SQRT2PI, inv_posdef, normal_cdf

# EP-related settings.
THRESHOLD = 1e-4

MAT_ONE = np.array([[1.0, -1.0], [-1.0, 1.0]])
MAT_ONE_FLAT = MAT_ONE.ravel()
```

--------------------------------

### Process Top-1 Lists

Source: https://choix.lum.li/en/latest/api

Functions to estimate model parameters from top-1 list data using LSR, I-LSR, optimization, and MM algorithms, including log-likelihood computation.

```Python
choix.lsr_top1()
choix.ilsr_top1()
choix.opt_top1()
choix.mm_top1()
choix.log_likelihood_top1()
```

--------------------------------

### Initialize LSR Markov Chain

Source: https://choix.lum.li/en/latest/_modules/choix/lsr

Initializes the weighted steady-state vector and transition matrix for the LSR Markov chain. It handles regularization via the alpha parameter and allows for initial parameter estimates.

```Python
def _init_lsr(n_items: int, alpha: float = 0.0, initial_params: NDArray[np.float64] | None = None) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """Initialize the weighted steady-state vector and transition matrix."""
    if initial_params is None:
        ws = np.ones(n_items)
    else:
        ws = initial_params
    chain = np.zeros((n_items, n_items))
    if alpha > 0:
        chain += alpha
    return ws, chain
```

--------------------------------

### Process Network Choices

Source: https://choix.lum.li/en/latest/api

Functions for estimating parameters of network choice models, including MAP estimates and log-likelihood computation.

```Python
choix.choicerank()
choix.log_likelihood_network()
```

--------------------------------

### Process Rankings

Source: https://choix.lum.li/en/latest/api

Functions for estimating model parameters from ranking data using LSR, I-LSR, optimization, and MM algorithms, along with log-likelihood calculation.

```Python
choix.lsr_rankings()
choix.ilsr_rankings()
choix.opt_rankings()
choix.mm_rankings()
choix.log_likelihood_rankings()
```

--------------------------------

### Log and Exponential Transformations

Source: https://choix.lum.li/en/latest/_modules/choix/utils

Provides functions for log-transforming weights into centered log-scale parameters and exp-transforming parameters back into weights. These are useful for working with probability distributions and ensuring numerical stability.

```Python
importmath
importrandom
importwarnings
fromcollections.abcimport Sequence
fromtypingimport Any, cast

importnumpyasnp
importscipy.linalgasspl
fromnumpy.typingimport ArrayLike, NDArray
fromscipy.linalgimport solve_triangular
fromscipy.specialimport logsumexp
fromscipy.statsimport kendalltau, rankdata

from.typingimport PairwiseData, RankingData, Top1Data

SQRT2 = math.sqrt(2.0)
SQRT2PI = math.sqrt(2.0 * math.pi)


deflog_transform(weights: ArrayLike) -> NDArray[np.float64]:
    """Transform weights into centered log-scale parameters."""
    params = np.log(weights)
    return params - params.mean()


defexp_transform(params: ArrayLike) -> NDArray[np.float64]:
    """Transform parameters into exp-scale weights."""
    params = np.asarray(params)
    weights = np.exp(params - np.mean(params))
    return (len(weights) / weights.sum()) * weights
```

--------------------------------

### ChoiceRank Algorithm for Network Models

Source: https://choix.lum.li/en/latest/_modules/choix/mm

Computes the maximum-a-posteriori (MAP) estimate of model parameters for a network choice model using the ChoiceRank algorithm. Requires a directed graph, traffic in/out data, and optional weight attribute, initial parameters, alpha, max iterations, and tolerance. Raises ImportError if NetworkX is not available.

```Python
def _choicerank(
    n_items: int,
    data: tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]],
    params: NDArray[np.float64],
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
"""Inner loop of ChoiceRank algorithm."""
    weights = exp_transform(params)
    adj, adj_t, traffic_in, traffic_out = data
    # First phase of message passing.
    zs = adj.dot(weights)
    # Second phase of message passing.
    with np.errstate(invalid="ignore"):
        denoms = adj_t.dot(traffic_out / zs)
    return traffic_in, denoms



[docs]
def choicerank(
    digraph: Any,
    traffic_in: NDArray[np.float64],
    traffic_out: NDArray[np.float64],
    weight: str | None = None,
    initial_params: NDArray[np.float64] | None = None,
    alpha: float = 1.0,
    max_iter: int = 10000,
    tol: float = 1e-8,
) -> NDArray[np.float64]:
"""Compute the MAP estimate of a network choice model's parameters.

    This function computes the maximum-a-posteriori (MAP) estimate of model
    parameters given a network structure and node-level traffic data (see
    :ref:`data-network`), using the ChoiceRank algorithm [MG17]_, [KTVV15]_.

    The nodes are assumed to be labeled using consecutive integers starting
    from 0.

    Parameters
    ----------
    digraph : networkx.DiGraph
        Directed graph representing the network.
    traffic_in : array_like
        Number of arrivals at each node.
    traffic_out : array_like
        Number of departures at each node.
    weight : str, optional
        The edge attribute that holds the numerical value used for the edge
        weight. If None (default) then all edge weights are 1.
    initial_params : array_like, optional
        Parameters used to initialize the iterative procedure.
    alpha : float, optional
        Regularization parameter.
    max_iter : int, optional
        Maximum number of iterations allowed.
    tol : float, optional
        Maximum L1-norm of the difference between successive iterates to
        declare convergence.

    Returns
    -------
    params : numpy.ndarray
        The MAP estimate of model parameters.

    Raises
    ------
    ImportError
        If the NetworkX library cannot be imported.
    """
    import networkx as nx

    # Compute the (sparse) adjacency matrix.
    n_items = len(digraph)
    nodes = np.arange(n_items)
    adj = nx.to_scipy_sparse_matrix(digraph, nodelist=nodes, weight=weight)  # pyright: ignore[reportArgumentType]
    adj_t = adj.T.tocsr()
    # Process the data into a standard form.
    traffic_in = np.asarray(traffic_in)
    traffic_out = np.asarray(traffic_out)
    data = (adj, adj_t, traffic_in, traffic_out)
    return _mm(n_items, data, initial_params, alpha, max_iter, tol, _choicerank)

```

--------------------------------

### Generate Parameters and Data

Source: https://choix.lum.li/en/latest/api

Functions for generating model parameters and creating comparison data based on different models like Bradley-Terry and Plackett-Luce.

```Python
choix.generate_params()
choix.generate_pairwise()
choix.generate_rankings()
choix.compare()
```

--------------------------------

### General Optimization Function (_opt)

Source: https://choix.lum.li/en/latest/_modules/choix/opt

The _opt function serves as a general-purpose optimizer that utilizes scipy.optimize.minimize. It supports both 'BFGS' and 'Newton-CG' methods and allows for initial parameters, maximum iterations, and tolerance for convergence. It requires numpy.

```Python
importmath
fromtypingimport Literal

importnumpyasnp
fromnumpy.typingimport NDArray
fromscipy.optimizeimport minimize
fromscipy.specialimport logsumexp

from.typingimport PairwiseData, RankingData, Top1Data
from.utilsimport softmax

MAXEXP = 500


def_safe_exp(x: float) -> float:
    x = max(min(x, MAXEXP), -MAXEXP)
    return math.exp(x)


classPairwiseFcts:
    """Optimization-related methods for pairwise-comparison data.

    This class provides methods to compute the negative log-likelihood (the
    "objective"), its gradient and its Hessian, given model parameters and
    pairwise-comparison data.
    """

    def__init__(self, data: PairwiseData, penalty: float):
        self._data = data
        self._penalty = penalty

    defobjective(self, params: NDArray[np.float64]) -> float:
        """Compute the negative penalized log-likelihood."""
        val = self._penalty * np.sum(params**2)
        for win, los in self._data:
            val += np.logaddexp(0, -(params[win] - params[los]))
        return val

    defgradient(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        grad = 2 * self._penalty * params
        for win, los in self._data:
            z = 1 / (1 + _safe_exp(params[win] - params[los]))
            grad[win] += -z
            grad[los] += +z
        return grad

    defhessian(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        hess = 2 * self._penalty * np.identity(len(params))
        for win, los in self._data:
            z = _safe_exp(params[win] - params[los])
            val = 1 / (1 / z + 2 + z)
            hess[(win, los), (los, win)] += -val
            hess[(win, los), (win, los)] += +val
        return hess


classTop1Fcts:
    """Optimization-related methods for top-1 data.

    This class provides methods to compute the negative log-likelihood (the
    "objective"), its gradient and its Hessian, given model parameters and
    top-1 data.

    The class also provides an alternative constructor for ranking data.
    """

    def__init__(self, data: Top1Data, penalty: float):
        self._data = data
        self._penalty = penalty

    @classmethod
    deffrom_rankings(cls, data: RankingData, penalty: float) -> "Top1Fcts":
        """Alternative constructor for ranking data."""
        top1 = list()
        for ranking in data:
            for i, winner in enumerate(ranking[:-1]):
                top1.append((winner, ranking[i + 1 :]))
        return cls(top1, penalty)

    defobjective(self, params: NDArray[np.float64]) -> float:
        """Compute the negative penalized log-likelihood."""
        val = self._penalty * np.sum(params**2)
        for winner, losers in self._data:
            idx = np.append(winner, losers)
            val += logsumexp(params.take(idx) - params[winner])  # pyright: ignore[reportOperatorIssue]
        return val

    defgradient(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        grad = 2 * self._penalty * params
        for winner, losers in self._data:
            idx = np.append(winner, losers)
            zs = softmax(params.take(idx))
            grad[idx] += zs
            grad[winner] += -1
        return grad

    defhessian(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        hess = 2 * self._penalty * np.identity(len(params))
        for winner, losers in self._data:
            idx = np.append(winner, losers)
            zs = softmax(params.take(idx))
            hess[np.ix_(idx, idx)] += -np.outer(zs, zs)
            hess[idx, idx] += zs
        return hess


def_opt(
    n_items: int,
    fcts: PairwiseFcts | Top1Fcts,
    method: Literal["BFGS", "Newton-CG"],
    initial_params: NDArray[np.float64] | None,
    max_iter: int | None,
    tol: float,
) -> NDArray[np.float64]:
    if initial_params is not None:
        x0 = initial_params
    else:
        x0 = np.zeros(n_items)
    if method == "Newton-CG":
        # `xtol`: Average relative error in solution xopt acceptable for
        # convergence [scipy doc].
        res = minimize(
            fcts.objective,
            x0,
            method=method,
            jac=fcts.gradient,
            hess=fcts.hessian,
            options={"xtol": tol, "maxiter": max_iter},
        )
    elif method == "BFGS":
        # `gtol`: Gradient norm must be less than gtol before successful
        # termination [scipy doc].
        res = minimize(
            fcts.objective,
            x0,
            method=method,
            jac=fcts.gradient,
            options={"gtol": tol, "maxiter": max_iter},
        )
    else:
        raise ValueError("method not known")
    return res.x

```

--------------------------------

### Choix Data Generation Functions

Source: https://choix.lum.li/en/latest/_sources/api

Provides functions to generate parameters and data for choice modeling. These include generating parameters, pairwise comparisons, rankings, and comparing different models.

```Python
choix.generate_params
choix.generate_pairwise
choix.generate_rankings
choix.compare
```

--------------------------------

### Top1Fcts Class for Top-1 Ranking Optimization

Source: https://choix.lum.li/en/latest/_modules/choix/opt

The Top1Fcts class is designed for optimizing top-1 data, which represents rankings. It offers an alternative constructor for ranking data and provides methods to compute the objective function, gradient, and Hessian. Dependencies include scipy and numpy.

```Python
importmath
fromtypingimport Literal

importnumpyasnp
fromnumpy.typingimport NDArray
fromscipy.optimizeimport minimize
fromscipy.specialimport logsumexp

from.typingimport PairwiseData, RankingData, Top1Data
from.utilsimport softmax

MAXEXP = 500


def_safe_exp(x: float) -> float:
    x = max(min(x, MAXEXP), -MAXEXP)
    return math.exp(x)


classTop1Fcts:
    """Optimization-related methods for top-1 data.

    This class provides methods to compute the negative log-likelihood (the
    "objective"), its gradient and its Hessian, given model parameters and
    top-1 data.

    The class also provides an alternative constructor for ranking data.
    """

    def__init__(self, data: Top1Data, penalty: float):
        self._data = data
        self._penalty = penalty

    @classmethod
    deffrom_rankings(cls, data: RankingData, penalty: float) -> "Top1Fcts":
        """Alternative constructor for ranking data."""
        top1 = list()
        for ranking in data:
            for i, winner in enumerate(ranking[:-1]):
                top1.append((winner, ranking[i + 1 :]))
        return cls(top1, penalty)

    defobjective(self, params: NDArray[np.float64]) -> float:
        """Compute the negative penalized log-likelihood."""
        val = self._penalty * np.sum(params**2)
        for winner, losers in self._data:
            idx = np.append(winner, losers)
            val += logsumexp(params.take(idx) - params[winner])  # pyright: ignore[reportOperatorIssue]
        return val

    defgradient(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        grad = 2 * self._penalty * params
        for winner, losers in self._data:
            idx = np.append(winner, losers)
            zs = softmax(params.take(idx))
            grad[idx] += zs
            grad[winner] += -1
        return grad

    defhessian(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        hess = 2 * self._penalty * np.identity(len(params))
        for winner, losers in self._data:
            idx = np.append(winner, losers)
            zs = softmax(params.take(idx))
            hess[np.ix_(idx, idx)] += -np.outer(zs, zs)
            hess[idx, idx] += zs
        return hess

```

--------------------------------

### Compute MAP Estimate of Network Choice Model Parameters

Source: https://choix.lum.li/en/latest/api

Computes the maximum-a-posteriori (MAP) estimate of model parameters for a network choice model using the ChoiceRank algorithm. It requires a directed graph, traffic input and output data, and optionally accepts edge weights, initial parameters, regularization, maximum iterations, and a tolerance for convergence. The function returns the estimated model parameters.

```python
choix.choicerank(_digraph :Any_, _traffic_in :ndarray[tuple[Any,...],dtype[float64]]_, _traffic_out :ndarray[tuple[Any,...],dtype[float64]]_, _weight :str|None=None_, _initial_params :ndarray[tuple[Any,...],dtype[float64]]|None=None_, _alpha :float=1.0_, _max_iter :int=10000_, _tol :float=1e-08_)
```

--------------------------------

### Compute ML Estimate using SciPy Optimize (Python)

Source: https://choix.lum.li/en/latest/api

Calculates the maximum-likelihood (ML) or maximum a-posteriori (MAP) estimate of model parameters for top-1 data using SciPy's optimization routines. Supports BFGS and Newton-CG methods, with options for regularization, initial parameters, and convergence criteria.

```Python
choix.opt_top1(_n_items :int_, _data :Sequence[tuple[int,tuple[int,...]]]_, _alpha :float=1e-06_, _method :Literal['BFGS','Newton-CG']='Newton-CG'_, _initial_params :ndarray[tuple[Any,...],dtype[float64]]|None=None_, _max_iter :int|None=None_, _tol :float=1e-05_)
```

--------------------------------

### Optimize Top-1 Data with SciPy

Source: https://choix.lum.li/en/latest/_modules/choix/opt

Calculates the maximum-likelihood estimate for model parameters using top-1 data. This function employs optimizers from scipy.optimize and includes regularization capabilities via the alpha parameter, enabling MAP estimates under a Gaussian prior. Key inputs include the number of items, top-1 data, regularization strength, optimization method, initial parameters, maximum iterations, and tolerance.

```Python
def opt_top1(
    n_items: int,
    data: Top1Data,
    alpha: float = 1e-6,
    method: Literal["BFGS", "Newton-CG"] = "Newton-CG",
    initial_params: NDArray[np.float64] | None = None,
    max_iter: int | None = None,
    tol: float = 1e-5,
):
    """Compute the ML estimate of model parameters using ``scipy.optimize``.

    This function computes the maximum-likelihood estimate of model parameters
    given top-1 data (see :ref:`data-top1`), using optimizers provided by the
    ``scipy.optimize`` module.

    If ``alpha > 0``, the function returns the maximum a-posteriori (MAP)
    estimate under an isotropic Gaussian prior with variance ``1 / alpha``. See
    :ref:`regularization` for details.

    Parameters
    ----------
    n_items : int
        Number of distinct items.
    data : list of lists
        Top-1 data.
    alpha : float, optional
        Regularization strength.
    method : str, optional
        Optimization method. Either "BFGS" or "Newton-CG".
    initial_params : array_like, optional
        Parameters used to initialize the iterative procedure.
    max_iter : int, optional
        Maximum number of iterations allowed.
    tol : float, optional
        Tolerance for termination (method-specific).

    Returns
    -------
    params : numpy.ndarray
        The (penalized) ML estimate of model parameters.

    Raises
    ------
    ValueError
        If the method is not "BFGS" or "Newton-CG".
    """
    fcts = Top1Fcts(data, alpha)
    return _opt(n_items, fcts, method, initial_params, max_iter, tol)
```

--------------------------------

### Process Pairwise Comparisons

Source: https://choix.lum.li/en/latest/api

Functions to estimate model parameters from pairwise comparison data using various methods including LSR, I-LSR, Rank Centrality, optimization, EP, and MM algorithms. Also includes log-likelihood computation.

```Python
choix.lsr_pairwise()
choix.ilsr_pairwise()
choix.lsr_pairwise_dense()
choix.ilsr_pairwise_dense()
choix.rank_centrality()
choix.opt_pairwise()
choix.ep_pairwise()
choix.mm_pairwise()
choix.log_likelihood_pairwise()
```

--------------------------------

### Compute ML Estimate using SciPy Optimize

Source: https://choix.lum.li/en/latest/api

Calculates the maximum-likelihood (ML) or maximum a-posteriori (MAP) estimate of model parameters using optimization methods from the SciPy library. Supports BFGS and Newton-CG methods, with options for regularization and convergence criteria.

```python
choix.opt_rankings(_n_items :int_, _data :Sequence[tuple[int,...]]_, _alpha :float=1e-06_, _method :Literal['BFGS','Newton-CG']='Newton-CG'_, _initial_params :ndarray[tuple[Any,...],dtype[float64]]|None=None_, _max_iter :int|None=None_, _tol :float=1e-05_)
```

--------------------------------

### Optimize Ranking Data with SciPy

Source: https://choix.lum.li/en/latest/_modules/choix/opt

Computes the maximum-likelihood estimate for model parameters given ranking data. This function leverages optimizers from scipy.optimize and supports regularization through the alpha parameter, allowing for MAP estimates with a Gaussian prior. It accepts parameters such as the number of items, ranking data, regularization strength, optimization method, initial parameters, maximum iterations, and tolerance.

```Python
def opt_rankings(
    n_items: int,
    data: RankingData,
    alpha: float = 1e-6,
    method: Literal["BFGS", "Newton-CG"] = "Newton-CG",
    initial_params: NDArray[np.float64] | None = None,
    max_iter: int | None = None,
    tol: float = 1e-5,
):
    """Compute the ML estimate of model parameters using ``scipy.optimize``.

    This function computes the maximum-likelihood estimate of model parameters
    given ranking data (see :ref:`data-rankings`), using optimizers provided by
    the ``scipy.optimize`` module.

    If ``alpha > 0``, the function returns the maximum a-posteriori (MAP)
    estimate under an isotropic Gaussian prior with variance ``1 / alpha``. See
    :ref:`regularization` for details.

    Parameters
    ----------
    n_items : int
        Number of distinct items.
    data : list of lists
        Ranking data.
    alpha : float, optional
        Regularization strength.
    method : str, optional
        Optimization method. Either "BFGS" or "Newton-CG".
    initial_params : array_like, optional
        Parameters used to initialize the iterative procedure.
    max_iter : int, optional
        Maximum number of iterations allowed.
    tol : float, optional
        Tolerance for termination (method-specific).

    Returns
    -------
    params : numpy.ndarray
        The (penalized) ML estimate of model parameters.

    Raises
    ------
    ValueError
        If the method is not "BFGS" or "Newton-CG".
    """
    fcts = Top1Fcts.from_rankings(data, alpha)
    return _opt(n_items, fcts, method, initial_params, max_iter, tol)
```

--------------------------------

### LSR Ranking for Full Rankings

Source: https://choix.lum.li/en/latest/_modules/choix/lsr

Computes Luce Spectral Ranking (LSR) estimates for full ranking data. It initializes transition rates using alpha for regularization and calculates parameters based on pairwise comparisons within rankings.

```python
def lsr_rankings(
    n_items: int,
    data: list[list[int]],
    alpha: float = 0.0,
    initial_params: NDArray[np.float64] | None = None,
) -> NDArray[np.float64]:
    """Compute the LSR estimate of model parameters.

    This function implements the Luce Spectral Ranking inference algorithm
    [MG15]_ for full ranking data (see :ref:`data-rankings`).

    The argument ``initial_params`` can be used to iteratively refine an
    existing parameter estimate (see the implementation of
    :func:`~choix.ilsr_rankings` for an idea on how this works). If it is set to
    `None` (the default), the all-ones vector is used.

    The transition rates of the LSR Markov chain are initialized with
    ``alpha``. When ``alpha > 0``, this corresponds to a form of regularization
    (see :ref:`regularization` for details).

    Parameters
    ----------
    n_items : int
        Number of distinct items.
    data : list of lists
        Ranking data.
    alpha : float, optional
        Regularization parameter.
    initial_params : array_like, optional
        Parameters used to build the transition rates of the LSR Markov chain.

    Returns
    -------
    params : numpy.ndarray
        An estimate of model parameters.
    """
    weights, chain = _init_lsr(n_items, alpha, initial_params)
    for ranking in data:
        sum_ = weights.take(ranking).sum()
        for i, winner in enumerate(ranking[:-1]):
            val = 1.0 / sum_
            for loser in ranking[i + 1 :]:
                chain[loser, winner] += val
            sum_ -= weights[winner]
    chain -= np.diag(chain.sum(axis=1))
    return log_transform(statdist(chain))
```

--------------------------------

### PairwiseFcts Class for Pairwise Comparison Optimization

Source: https://choix.lum.li/en/latest/_modules/choix/opt

The PairwiseFcts class handles optimization for pairwise comparison data. It calculates the negative penalized log-likelihood, its gradient, and its Hessian, which are essential for optimization algorithms. It requires scipy and numpy libraries.

```Python
importmath
fromtypingimport Literal

importnumpyasnp
fromnumpy.typingimport NDArray
fromscipy.optimizeimport minimize
fromscipy.specialimport logsumexp

from.typingimport PairwiseData, RankingData, Top1Data
from.utilsimport softmax

MAXEXP = 500


def_safe_exp(x: float) -> float:
    x = max(min(x, MAXEXP), -MAXEXP)
    return math.exp(x)


classPairwiseFcts:
    """Optimization-related methods for pairwise-comparison data.

    This class provides methods to compute the negative log-likelihood (the
    "objective"), its gradient and its Hessian, given model parameters and
    pairwise-comparison data.
    """

    def__init__(self, data: PairwiseData, penalty: float):
        self._data = data
        self._penalty = penalty

    defobjective(self, params: NDArray[np.float64]) -> float:
        """Compute the negative penalized log-likelihood."""
        val = self._penalty * np.sum(params**2)
        for win, los in self._data:
            val += np.logaddexp(0, -(params[win] - params[los]))
        return val

    defgradient(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        grad = 2 * self._penalty * params
        for win, los in self._data:
            z = 1 / (1 + _safe_exp(params[win] - params[los]))
            grad[win] += -z
            grad[los] += +z
        return grad

    defhessian(self, params: NDArray[np.float64]) -> NDArray[np.float64]:
        hess = 2 * self._penalty * np.identity(len(params))
        for win, los in self._data:
            z = _safe_exp(params[win] - params[los])
            val = 1 / (1 / z + 2 + z)
            hess[(win, los), (los, win)] += -val
            hess[(win, los), (win, los)] += +val
        return hess

```

--------------------------------

### MM Algorithm for Top-1 Data

Source: https://choix.lum.li/en/latest/_modules/choix/mm

Computes the maximum-likelihood (ML) estimate of model parameters using the MM algorithm for top-1 data. Supports regularization with a Dirichlet prior. Requires number of items, top-1 data, and optional initial parameters, alpha, max iterations, and tolerance.

```Python
def _mm_top1(
    n_items: int,
    data: Top1Data,
    params: NDArray[np.float64],
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
"""Inner loop of MM algorithm for top1 data."""
    weights = exp_transform(params)
    wins = np.zeros(n_items, dtype=float)
    denoms = np.zeros(n_items, dtype=float)
    for winner, losers in data:
        wins[winner] += 1
        val = 1 / (weights.take(losers).sum() + weights[winner])
        for item in itertools.chain([winner], losers):
            denoms[item] += val
    return wins, denoms



[docs]
def mm_top1(
    n_items: int,
    data: Top1Data,
    initial_params: NDArray[np.float64] | None = None,
    alpha: float = 0.0,
    max_iter: int = 10000,
    tol: float = 1e-8,
) -> NDArray[np.float64]:
"""Compute the ML estimate of model parameters using the MM algorithm.

    This function computes the maximum-likelihood (ML) estimate of model
    parameters given top-1 data (see :ref:`data-top1`), using the
    minorization-maximization (MM) algorithm [Hun04]_, [CD12]_

    If ``alpha > 0``, the function returns the maximum a-posteriori (MAP)
    estimate under a (peaked) Dirichlet prior. See :ref:`regularization` for
    details.

    Parameters
    ----------
    n_items : int
        Number of distinct items.
    data : list of lists
        Top-1 data.
    initial_params : array_like, optional
        Parameters used to initialize the iterative procedure.
    alpha : float, optional
        Regularization parameter.
    max_iter : int, optional
        Maximum number of iterations allowed.
    tol : float, optional
        Maximum L1-norm of the difference between successive iterates to
        declare convergence.

    Returns
    -------
    params : numpy.ndarray
        The ML estimate of model parameters.
    """
    return _mm(n_items, data, initial_params, alpha, max_iter, tol, _mm_top1)

```

--------------------------------

### Generate Pairwise Comparisons (Bradley-Terry)

Source: https://choix.lum.li/en/latest/_modules/choix/utils

Samples pairwise comparison outcomes from a Bradley-Terry model. It randomly selects pairs of items and determines the winner based on provided model parameters. Dependencies include numpy and random.

```Python
import numpy as np
import random

def generate_comparisons(params, n_comparisons):
    params = np.asarray(params)
    n = len(params)
    items = tuple(range(n))
    data = list()
    for _ in range(n_comparisons):
        # Pick the pair uniformly at random.
        a, b = random.sample(items, 2)
        if compare((a, b), params) == a:
            data.append((a, b))
        else:
            data.append((b, a))
    return tuple(data)
```

--------------------------------

### Compute EP for Pairwise Data (Python)

Source: https://choix.lum.li/en/latest/_modules/choix/ep

Computes an approximate Bayesian posterior probability distribution over model parameters using the EP algorithm for pairwise comparison data. Supports 'logit' and 'probit' models, taking the number of items, comparison data, prior variance (alpha), model type, max iterations, and optional initial state as input. Returns the mean and covariance of the approximate posterior.

```Python
def ep_pairwise(
    n_items: int,
    data: PairwiseData,
    alpha: float,
    model: Literal["logit", "probit"] = "logit",
    max_iter: int = 100,
    initial_state: tuple | None = None,
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """Compute a distribution of model parameters using the EP algorithm.

    This function computes an approximate Bayesian posterior probability
    distribution over model parameters, given pairwise-comparison data (see
    :ref:`data-pairwise`). It uses the expectation propagation algorithm, as
    presented, e.g., in [CG05]_.

    The prior distribution is assumed to be isotropic Gaussian with variance
    ``1 / alpha``. The posterior is approximated by a a general multivariate
    Gaussian distribution, described by a mean vector and a covariance matrix.

    Two different observation models are available. ``logit`` (default) assumes
    that pairwise-comparison outcomes follow from a Bradley-Terry model.
    ``probit`` assumes that the outcomes follow from Thurstone's model.

    Parameters
    ----------
    n_items : int
        Number of distinct items.
    data : list of lists
        Pairwise-comparison data.
    alpha : float
        Inverse variance of the (isotropic) prior.
    model : str, optional
        Observation model. Either "logit" or "probit".
    max_iter : int, optional
        Maximum number of iterations allowed.
    initial_state : tuple of array_like, optional
        Natural parameters used to initialize the EP algorithm.

    Returns
    -------
    mean : numpy.ndarray
        The mean vector of the approximate Gaussian posterior.
    cov : numpy.ndarray
        The covariance matrix of the approximate Gaussian posterior.

    Raises
    ------
    ValueError
        If the observation model is not "logit" or "probit".
    """
    if model == "logit":
        match_moments = _match_moments_logit
    elif model == "probit":
        match_moments = _match_moments_probit
    else:
        raise ValueError("unknown model '{}'".format(model))
    return _ep_pairwise(n_items, data, alpha, match_moments, max_iter, initial_state)

```

--------------------------------

### Generate Random Model Parameters

Source: https://choix.lum.li/en/latest/_modules/choix/utils

Generates random model parameters, sampling uniformly for each item within a specified interval. Optionally, the parameters can be returned in ordered (ascending) sequence and are normalized to have a mean of zero.

```Python
import numpy as np
from numpy.typing import NDArray

def generate_params(
    n_items: int,
    interval: float = 5.0,
    ordered: bool = False
) -> NDArray[np.float64]:
    r"""Generate random model parameters.

    This function samples a parameter independently and uniformly for each
    item. ``interval`` defines the width of the uniform distribution.

    Parameters
    ----------
    n_items : int
        Number of distinct items.
    interval : float
        Sampling interval.
    ordered : bool, optional
        If true, the parameters are ordered from lowest to highest.

    Returns
    -------
    params : numpy.ndarray
       Model parameters.
    """
    params = np.random.uniform(low=0, high=interval, size=n_items)
    if ordered:
        params.sort()
    return params - params.mean()
```

--------------------------------

### Internal EP Calculation Logic (Python)

Source: https://choix.lum.li/en/latest/_modules/choix/ep

This internal function performs the core Expectation Propagation calculations. It initializes parameters, iterates to update them based on pairwise comparisons, and uses moment matching to refine the approximation. It includes error handling for non-convergence.

```Python
def _ep_pairwise(
    n_items: int,
    comparisons: PairwiseData,
    alpha: float,
    match_moments: Callable[[float, float], tuple[float, float, float]],
    max_iter: int,
    initial_state: tuple | None,
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """Compute a distribution of model parameters using the EP algorithm.

    Raises
    ------
    RuntimeError
        If the algorithm does not converge after ``max_iter`` iterations.
    """
    # Static variable that allows to check the # of iterations after the call.
    _ep_pairwise.iterations = 0  # pyright: ignore[reportFunctionMemberAccess]
    m = len(comparisons)
    prior_inv = alpha * np.eye(n_items)
    if initial_state is None:
        # Initially, mean and covariance come from the prior.
        mean = np.zeros(n_items)
        cov = (1 / alpha) * np.eye(n_items)
        # Initialize the natural params in the function space.
        tau = np.zeros(m)
        nu = np.zeros(m)
        # Initialize the natural params in the space of thetas.
        prec = np.zeros((n_items, n_items))
        xs = np.zeros(n_items)
    else:
        tau, nu = initial_state
        mean, cov, xs, prec = _init_ws(n_items, comparisons, prior_inv, tau, nu)
    for _ in range(max_iter):
        _ep_pairwise.iterations += 1  # pyright: ignore[reportFunctionMemberAccess]
        # Keep a copy of the old parameters for convergence testing.
        tau_old = np.array(tau, copy=True)
        nu_old = np.array(nu, copy=True)
        for i in nprand.permutation(m):
            a, b = comparisons[i]
            # Update mean and variance in function space.
            f_var = cov[a, a] + cov[b, b] - 2 * cov[a, b]
            f_mean = mean[a] - mean[b]
            # Cavity distribution.
            tau_tot = 1.0 / f_var
            nu_tot = tau_tot * f_mean
            tau_cav: float = tau_tot - tau[i]  # pyright: ignore[reportAssignmentType]
            nu_cav: float = nu_tot - nu[i]  # pyright: ignore[reportAssignmentType]
            cov_cav: float = 1.0 / tau_cav
            mean_cav: float = cov_cav * nu_cav
            # Moment matching.

```

--------------------------------

### Compute Log-Likelihood for Network Data

Source: https://choix.lum.li/en/latest/_modules/choix/utils

Calculates the log-likelihood for model parameters in a network context, considering traffic in and out of nodes. It supports weighted edges for more accurate likelihood calculations.

```Python
import numpy as np
from numpy.typing import NDArray
from typing import Sequence, Any
from scipy.special import logsumexp

# Assuming digraph is a graph object with methods like out_degree and successors
# Assuming traffic_in and traffic_out are sequences of traffic values

def log_likelihood_network(
    digraph: Any,
    traffic_in: Sequence,
    traffic_out: Sequence,
    params: NDArray[np.float64],
    weight: str | None = None,
) -> float:
    """
    Compute the log-likelihood of model parameters.

    If ``weight`` is not ``None``, the log-likelihood is correct only up to a
    constant (independent of the parameters).
    """
    loglik: float = 0.0
    for i in range(len(traffic_in)):
        loglik += traffic_in[i] * params[i]
        if digraph.out_degree(i) > 0:
            neighbors = list(digraph.successors(i))
            if weight is None:
                loglik -= traffic_out[i] * logsumexp(params.take(neighbors))
            else:
                weights = [digraph[i][j][weight] for j in neighbors]
                loglik -= traffic_out[i] * logsumexp(params.take(neighbors), b=weights)
    return loglik
```

--------------------------------

### Choix Utility Functions

Source: https://choix.lum.li/en/latest/_sources/api

Offers utility functions for choice modeling, including calculating probabilities and various distance metrics such as the footrule and Kendall tau distances.

```Python
choix.probabilities
choix.footrule_dist
choix.kendalltau_dist
```