Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[solidago] feat: Asymetric uncertainty #1781

Closed
wants to merge 2 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably log of 1/2 instead

np.log(theta_ab / (positive_exponential_term - negative_exponential_term)),
).sum()



@njit
def Delta_theta(theta_ab):
Expand All @@ -39,6 +53,64 @@ def Delta_theta(theta_ab):
).sum() ** (-0.5)


HIGH_LIKELIHOOD_RANGE_THRESHOLD = 1
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lenhoanglnh Do you confirm the idea of using a Likelihood lower bound to compute the uncertainty interval, rather than a more standard 90% confidence interval?



@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(
Expand Down Expand Up @@ -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,
}
),
]
Expand All @@ -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,
Copy link
Member Author

@lfaucon lfaucon Sep 24, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lenhoanglnh When calculating the high likelihood range, should we include only the likelihood of observed comparisons, or also the prior/regularization (something like alpha * theta^2)?

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,
)
Expand Down
76 changes: 70 additions & 6 deletions solidago/src/solidago/solvers/optimize.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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`")
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

optimize by not calling f(a) and f(b) twice


if xtol <= 0:
raise ValueError("xtol is too small (<= 0)")
Expand Down
98 changes: 97 additions & 1 deletion solidago/tests/test_comparisons_to_scores.py
Original file line number Diff line number Diff line change
@@ -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)
Expand Down Expand Up @@ -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
4 changes: 2 additions & 2 deletions solidago/tests/test_resilient_primitives.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Loading
Loading