[solidago] feat: Asymetric uncertainty

solidago/src/solidago/comparisons_to_scores/continuous_bradley_terry.py

@@ -4,13 +4,14 @@
 using coordinate descent.
 """
 import random
+from typing import Tuple, Callable
 
 import numpy as np
 import pandas as pd
 from numba import njit
 
 from solidago.comparisons_to_scores.base import ComparisonsToScoresAlgorithm
-from solidago.solvers.optimize import brentq
+from solidago.solvers.optimize import brentq, SignChangeIntervalNotFoundError
 
 
 DEFAULT_ALPHA = 0.20  # Signal-to-noise hyperparameter
@@ -29,6 +30,19 @@ def contributor_loss_partial_derivative(theta_a, theta_b, r_ab, alpha):
         + r_ab
     )
 
+@njit
+def continuous_bradley_terry_log_likelihood(theta_a, theta_b, r_ab, r_max) -> float:
+    theta_ab = theta_a - theta_b
+    normalized_r_ab = r_ab / r_max
+    positive_exponential_term = np.exp((normalized_r_ab + 1) * theta_ab)
+    negative_exponential_term = np.exp((normalized_r_ab - 1) * theta_ab)
+    return np.where(
+        np.abs(theta_ab) < EPSILON,
+        1 / 2,
+        np.log(theta_ab / (positive_exponential_term - negative_exponential_term)),
+    ).sum()
+
+
 
 @njit
 def Delta_theta(theta_ab):
@@ -39,6 +53,64 @@ def Delta_theta(theta_ab):
     ).sum() ** (-0.5)
 
 
+HIGH_LIKELIHOOD_RANGE_THRESHOLD = 1
+
+
+@njit
+def translated_function(x, f, translation, args=()):
+    """Returns the function x => f(x) - translation"""
+    return f(x, *args) - translation
+
+def get_high_likelihood_range(
+    log_likelihood,
+    maximum_a_posteriori: float,
+    threshold: float = HIGH_LIKELIHOOD_RANGE_THRESHOLD,
+    args=(),
+) -> Tuple[float, float]:
+    """
+    Find a root of a function in a bracketing interval using Brent's method
+    adapted from Scipy's brentq.
+    Uses the classic Brent's method to find a zero of the function `f` on
+    the sign changing interval [a , b].
+
+    Parameters
+    ----------
+    likelihood_function:
+        Python function computing a log likelihood.
+        `f` must be continuous and concave.
+        `f` must be jitted via numba.
+    maximum_a_posteriori:
+        The high liklihood position selected as most likely based on the prior
+        distribution and the observed likelihood
+    threshold:
+        The threshold used to compute the high likelihood range. The range will
+        be the interval with where we have
+        log_likelihood > log_likelihood(maximum_a_posteriori) - threshold
+        The threshold must be strictly positive.
+
+    Returns
+    -------
+    interval:
+        A tuple of float representing an interval containing the
+        maximum_a_posteriori.
+    """
+    if threshold <= 0:
+        raise ValueError("`threshold` must be strictly positive")
+    log_likelihood_at_maximum_a_posteriori = log_likelihood(maximum_a_posteriori, *args)
+    min_log_likelihood = log_likelihood_at_maximum_a_posteriori - threshold
+
+    try:
+        lower_bound = brentq(translated_function, a=maximum_a_posteriori-1, b=maximum_a_posteriori, search_b=False, args=(log_likelihood, min_log_likelihood, args))
+    except SignChangeIntervalNotFoundError:
+        lower_bound = -np.inf
+    try: 
+        upper_bound = brentq(translated_function, a=maximum_a_posteriori, b=maximum_a_posteriori+1, search_a=False, args=(log_likelihood, min_log_likelihood, args))
+    except SignChangeIntervalNotFoundError:
+        upper_bound = np.inf
+
+    return lower_bound, upper_bound
+
+
 @njit
 def coordinate_optimize(r_ab, theta_b, precision, alpha):
     return brentq(
@@ -91,17 +163,17 @@ def pick_next_coordinate():
                 unchanged.clear()
         return theta
 
-    def compute_individual_scores(self, scores: pd.DataFrame, initial_entity_scores=None):
-        scores = scores[["entity_a", "entity_b", "score"]]
+    def compute_individual_scores(self, comparison_scores: pd.DataFrame, initial_entity_scores=None):
+        comparison_scores = comparison_scores[["entity_a", "entity_b", "score"]]
         scores_sym = (
             pd.concat(
                 [
-                    scores,
+                    comparison_scores,
                     pd.DataFrame(
                         {
-                            "entity_a": scores.entity_b,
-                            "entity_b": scores.entity_a,
-                            "score": -1 * scores.score,
+                            "entity_a": comparison_scores.entity_b,
+                            "entity_b": comparison_scores.entity_a,
+                            "score": -1 * comparison_scores.score,
                         }
                     ),
                 ]
@@ -124,14 +196,25 @@ def compute_individual_scores(self, scores: pd.DataFrame, initial_entity_scores=
             initial_scores = initial_scores.to_numpy()
         theta_star_numpy = self.coordinate_descent(coord_to_subset, initial_scores=initial_scores)
         delta_star_numpy = np.zeros(len(theta_star_numpy))
+        raw_score_lower_bound = np.zeros(len(theta_star_numpy))
+        raw_score_upper_bound = np.zeros(len(theta_star_numpy))
         for idx_a in range(len(theta_star_numpy)):
-            indices_b, _r_ab = coord_to_subset[idx_a]
+            indices_b, r_ab = coord_to_subset[idx_a]
+            lower_bound, upper_bound = get_high_likelihood_range(
+                continuous_bradley_terry_log_likelihood,
+                theta_star_numpy[idx_a],
+                args=(theta_star_numpy[indices_b], r_ab, self.r_max),
+            )
+            raw_score_lower_bound[idx_a] = lower_bound
+            raw_score_upper_bound[idx_a] = upper_bound
             delta_star_numpy[idx_a] = Delta_theta(theta_star_numpy[idx_a] - theta_star_numpy[indices_b])
 
         result = pd.DataFrame(
             {
                 "raw_score": theta_star_numpy,
                 "raw_uncertainty": delta_star_numpy,
+                "raw_score_lower_bound": raw_score_lower_bound,
+                "raw_score_upper_bound": raw_score_upper_bound,
             },
             index=entities_index,
         )

solidago/src/solidago/solvers/optimize.py

@@ -7,7 +7,7 @@
 All rights reserved.
 """
 # pylint: skip-file
-from typing import Tuple
+from typing import Tuple, Callable, Optional
 
 import numpy as np
 from numba import njit
@@ -39,9 +39,64 @@ def _bisect_interval(a, b, fa, fb) -> Tuple[float, int]:
 
     return root, status
 
+class SignChangeIntervalNotFoundError(RuntimeError):
+    pass
 
 @njit
-def brentq(f, args=(), xtol=_xtol, rtol=_rtol, maxiter=_iter, disp=True, a: float=-1.0, b: float=1.0) -> float:
+def search_sign_change_interval(
+    f: Callable,
+    a: float,
+    b: float,
+    args: Tuple = (),
+    max_iterations: int = 32,
+    search_a: bool = True,
+    search_b: bool = True,
+):
+    """
+        Searches bounds a and b of interval where `f` changes sign. This is
+        achieved by increasing the size of the interval iteratively.
+        Note that the method is not guaranteed to succeed for most functions
+        and highly depends on the initial bounds.
+
+        Parameters
+        ----------
+        f : jitted and callable
+            Python function returning a number.  `f` must be continuous.
+        a : number
+            One end of the bracketing interval [a,b].
+        b : number
+            The other end of the bracketing interval [a,b].
+        args : tuple, optional(default=())
+            Extra arguments to be used in the function call.
+        max_iterations:
+            The maximum number of iteration in the search. /!\ When using a
+            large number of iterations, bounds would become very large and
+            functions may not be well behaved.
+        search_a:
+            If true, the value of `a` provided will be updated to search for an
+            interval where `f` changes sign
+        search_b:
+            If true, the value of `b` provided will be updated to search for an
+            interval where `f` changes sign
+
+        Returns
+        -------
+        a, b:
+            An interval on which the continuous function `f` changes sign
+    """
+    if a >= b:
+        raise ValueError(f"Initial interval bounds should be such that a < b. Found a={a} and b={b}")
+    iteration_count = 0
+    while f(a, *args) * f(b, *args) > 0:
+        if iteration_count > max_iterations:
+            raise SignChangeIntervalNotFoundError("Could not find a sign changing interval")
+        iteration_count+=1
+        a = a-(b-a) if search_a else a
+        b = b+(b-a) if search_b else b
+    return a, b
+
+@njit
+def brentq(f, args=(), xtol=_xtol, rtol=_rtol, maxiter=_iter, disp=True, a: float=-1.0, b: float=1.0, search_a: bool=True, search_b: bool = True) -> float:
     """
     Find a root of a function in a bracketing interval using Brent's method
     adapted from Scipy's brentq.
@@ -69,14 +124,23 @@ def brentq(f, args=(), xtol=_xtol, rtol=_rtol, maxiter=_iter, disp=True, a: floa
         Maximum number of iterations.
     disp : bool, optional(default=True)
         If True, raise a RuntimeError if the algorithm didn't converge.
+    search_a:
+        If true, the value of `a` provided will be updated to search for an
+        interval where `f` changes sign
+    search_b:
+        If true, the value of `b` provided will be updated to search for an
+        interval where `f` changes sign
     Returns
     -------
     root : float
     """
-    while f(a, *args) > 0:
-        a = a - 2 * (b-a)
-    while f(b, *args) < 0:
-        b = b + 2 * (b-a)
+    a, b = search_sign_change_interval(f, a, b, args=args, search_a=search_a, search_b=search_b)
+    if f(a, *args) == 0:
+        return a
+    if f(b, *args) == 0:
+        return b
+    if f(a, *args) * f(b, *args) > 0:
+        raise ValueError("Function `f` should have opposite sign on bounds `a` and `b`")
 
     if xtol <= 0:
         raise ValueError("xtol is too small (<= 0)")

solidago/tests/test_comparisons_to_scores.py

@@ -1,10 +1,13 @@
 import random
 
+from numba import njit
+import numpy as np
 import pandas as pd
 import pytest
 
-from solidago.comparisons_to_scores import ContinuousBradleyTerry, HookeIndividualScores
 
+from solidago.comparisons_to_scores import ContinuousBradleyTerry, HookeIndividualScores
+from solidago.comparisons_to_scores.continuous_bradley_terry import get_high_likelihood_range
 
 matrix_inversion = HookeIndividualScores(r_max=10)
 continuous_bradley_terry = ContinuousBradleyTerry(r_max=10)
@@ -145,3 +148,96 @@ def test_comparisons_non_connex(self, algo):
         assert scores.loc["B"].raw_score == pytest.approx(scores.loc["F"].raw_score, abs=1e-4)
         assert scores.loc["C"].raw_score == pytest.approx(scores.loc["G"].raw_score, abs=1e-4)
         assert scores.loc["D"].raw_score == pytest.approx(scores.loc["H"].raw_score, abs=1e-4)
+
+
+@pytest.mark.parametrize("comparisons",[
+    pd.DataFrame({
+        "entity_a": ["A", "A", "B", "C"],
+        "entity_b": ["B", "C", "C", "D"],
+        "score": [4, 4, 2, 4]
+    }),
+    pd.DataFrame({
+        "entity_a": ["X", "Z", "B", "C", "A", "A"],
+        "entity_b": ["Y", "X", "C", "D", "B", "C"],
+        "score": [4, 4, 2, 0, -1, -11]
+    }),
+    pd.DataFrame({
+        "entity_a": ["X", "Z", "B", "C", "A", "A"],
+        "entity_b": ["Y", "X", "C", "D", "B", "C"],
+        "score": [-11, 11, -11, 11, -11, 11],
+    })
+])
+def test_gbt_uncertainty_uncertainty_bounds_contain_score(comparisons):
+    model = ContinuousBradleyTerry(r_max=11)
+    individual_scores = model.compute_individual_scores(comparisons)
+    rows = [row for _, row in individual_scores.iterrows()]
+    assert all(
+        row["raw_score_lower_bound"] <= row["raw_score"] <= row["raw_score_upper_bound"]
+        for row in rows
+    )
+
+def test_gbt_uncertainty_has_infinite_range_with_max_comparison():
+    model = ContinuousBradleyTerry(r_max=4)
+    comparisons = pd.DataFrame({
+        "entity_a": ["A", "A", "B", "C"],
+        "entity_b": ["B", "C", "C", "D"],
+        "score": [4, 4, 2, 4]
+    })
+    individual_scores = model.compute_individual_scores(comparisons)
+    assert individual_scores.raw_score_lower_bound.loc["A"] == -np.inf
+    assert individual_scores.raw_score_upper_bound.loc["A"] < np.inf
+    assert individual_scores.raw_score_upper_bound.loc["D"] == np.inf
+
+
+def test_gbt_uncertainty_has_finite_range_with_non_max_comparison():
+    model = ContinuousBradleyTerry(r_max=5)
+    comparisons = pd.DataFrame({
+        "entity_a": ["A", "A", "B", "C"],
+        "entity_b": ["B", "C", "C", "D"],
+        "score": [4, 4, 2, 4]
+    })
+    individual_scores = model.compute_individual_scores(comparisons)
+    assert individual_scores.raw_score_lower_bound.loc["A"] != -np.inf
+    assert individual_scores.raw_score_upper_bound.loc["A"] != np.inf
+    assert individual_scores.raw_score_lower_bound.loc["D"] != -np.inf
+    assert individual_scores.raw_score_upper_bound.loc["D"] != np.inf
+
+
+@njit
+def quadratic_log_likelihood_mock(x):
+    return -10 - (x-3)**2
+
+@njit
+def asymetrix_log_likelihood_mock(x):
+    return -10 - (x-3)**2 if x > 3 else -10 - (x-3)**2 / 9
+
+@pytest.mark.parametrize("log_likelihood,max_a_posteriori,threshold,expected_range", [
+    (quadratic_log_likelihood_mock, 3, 1, (2, 4)),
+    (asymetrix_log_likelihood_mock, 3, 1, (0, 4)),
+    (quadratic_log_likelihood_mock, 3, 9, (0, 6)),
+    (quadratic_log_likelihood_mock, 0, 7, (-1, 7)),
+])
+def test_get_high_likelihood_range_has_expected_value(log_likelihood, max_a_posteriori,threshold,expected_range):
+    found_lower, found_upper = get_high_likelihood_range(log_likelihood, max_a_posteriori, threshold)
+    expected_lower, expected_upper =  expected_range
+    assert found_lower == pytest.approx(expected_lower)
+    assert found_upper == pytest.approx(expected_upper)
+
+
+@njit
+def lost_against_magnus_log_likelihood_mock(x):
+    # This represents the likelihood of losing a game against world champion
+    # Magnus Carlsen. This observation would tell you that the elo score `x`
+    # cannot be extermely high, but it may very likely be arbitrarily low.
+    # In other words, with such a likelihood function, we expect infinite
+    # uncertainty on the left. 
+    L = 0.005 
+    MAGNUS_ELO = 2839  # Source https://ratings.fide.com/profile/1503014 as of Sep 23rd
+    return np.log(1-1./(1 + np.exp(L * (MAGNUS_ELO - x))))
+
+def test_get_high_likelihood_range_supports_infinite_uncertainty():
+    found_lower, _ = get_high_likelihood_range(
+        lost_against_magnus_log_likelihood_mock,
+        2000,
+    )
+    assert found_lower == -np.inf

solidago/tests/test_resilient_primitives.py

@@ -154,8 +154,8 @@ def test_qrunc_with_W_equals_zero(self):
         (10, np.array([1] * 1000), np.array([-1] * 500 + [1] * 500), np.array([0.1] * 1000), 0.10, -1.0712),
         (10, np.array([1] * 1000), np.array([-1] * 100 + [1] * 900), np.array([1e-6] * 1000), 0.10, 0.),
         (0.0001, np.array([1] * 1000), np.array([-1] * 102 + [1] * 898), np.array([1e-6] * 1000), 0.01, -1),
-        (1e-12, np.array([1000] * 1000), np.arange(1000, 2000, 1), np.array([1e-6] * 1000), 0.90, 1899.3929),
+        (1e-12, np.array([1000] * 1000), np.arange(1000, 2000, 1), np.array([1e-6] * 1000), 0.90, 1899.39),
     ]
 )
 def test_qr_quantile_returns_expected_results(W,w,x,delta,quantile,expected_result):
-    assert pytest.approx(expected_result, abs=1e-4) == QrQuantile(W,w,x,delta,quantile)
+    assert pytest.approx(expected_result, abs=1e-2) == QrQuantile(W,w,x,delta,quantile)