LORD: Levels based On Recent Discovery¶

class LordThree(AbstractSequentialTest):
    """LORD 3: Online FDR control based on recent discovery with wealth dynamics.

    LORD 3 is a variant of the LORD (significance Levels based On Recent Discovery)
    procedure for online FDR control. The test levels depend on the past only through
    the time of the last discovery and the wealth accumulated at that time.

    LORD procedures have an intuitive interpretation: they start with an error budget
    (alpha-wealth), pay a price each time a hypothesis is tested, and earn back wealth
    when discoveries are made. LORD 3 sets thresholds based on the time since the last
    discovery and the wealth at that time.

    Note:
        This method was superseded by LORD++ and is implemented for demonstrative
        purposes and comparison studies. For practical applications, consider using
        LORD++ or more recent methods like ADDIS or SAFFRON.

    Args:
        alpha: Target FDR level (e.g., 0.05 for 5% FDR). Must be in (0, 1).
        wealth: Initial alpha-wealth for purchasing rejection thresholds.
                Must satisfy 0 < wealth < alpha.
        reward: Reward earned back for each discovery. Must be positive.
                Typical choice is reward = alpha - wealth.

    Attributes:
        wealth: Current alpha-wealth available for testing.
        reward: Fixed reward earned per discovery.
        last_reject: Index of the most recent rejection (0 if none).
        wealth_reject: Alpha-wealth at the time of the last rejection.

    Examples:
        >>> # Basic usage with recommended parameters
        >>> lord3 = LordThree(alpha=0.05, wealth=0.025, reward=0.025)
        >>> decision = lord3.test_one(0.01)  # Test a small p-value
        >>> print(f"Rejected: {decision}")

        >>> # Sequential testing
        >>> p_values = [0.001, 0.3, 0.02, 0.8, 0.005]
        >>> decisions = [lord3.test_one(p) for p in p_values]
        >>> discoveries = sum(decisions)

    References:
        Javanmard, A., and A. Montanari (2018). "Online rules for control of false
        discovery rate and false discovery exceedance." Annals of Statistics,
        46(2):526-554.

        Project Euclid: https://projecteuclid.org/journals/annals-of-statistics/volume-46/issue-2/Online-rules-for-control-of-false-discovery-rate-and-false/10.1214/17-AOS1559.full
    """

    def __init__(
        self,
        alpha: float,
        wealth: float,
        reward: float,
    ):  # fmt: skip
        super().__init__(alpha)
        self.wealth: float = wealth
        self.reward: float = reward
        validity.check_initial_wealth(wealth, alpha)
        validity.check_reward_budget(wealth=wealth, reward=reward, alpha=alpha)

        self.seq = DefaultLordGammaSequence(c=0.07720838)

        self.last_reject: int = 0  # reject index
        self.wealth_reject: float = wealth  # reject wealth
        # Matches onlineFDR's LORD(version = "3") state recursion sentinel.
        self._decision_history: list[bool] = [True]

    def test_one(self, p_val: float) -> bool:
        """Test a single p-value using the LORD 3 procedure.

        The LORD 3 algorithm processes p-values sequentially:
        1. Calculate threshold based on time since last discovery and wealth at that time
        2. Spend alpha-wealth equal to the threshold
        3. Earn back reward if discovery is made
        4. Update last rejection time and wealth if discovery is made

        Args:
            p_val: P-value to test. Must be in [0, 1].

        Returns:
            True if the null hypothesis is rejected (discovery), False otherwise.

        Raises:
            ValueError: If p_val is not in [0, 1].

        Examples:
            >>> lord3 = LordThree(alpha=0.05, wealth=0.025, reward=0.025)
            >>> lord3.test_one(0.01)  # Small p-value, likely rejected
            True
            >>> lord3.test_one(0.8)   # Large p-value, not rejected
            False
        """
        validity.check_p_val(p_val)
        self.num_test += 1

        self.alpha = (
            self.seq.calc_gamma(self.num_test - self.last_reject)  # fmt: skip
            * self.wealth_reject
        )

        is_rejected = p_val <= self.alpha

        # onlineFDR applies the spend as min(alpha_t, W_{t-1}).
        spend = min(self.alpha, self.wealth)

        # onlineFDR LORD-3 recursion credits reward with one-step lag and
        # sentinel R_0 = TRUE on t=2.
        reward_credited = (
            is_rejected
            if self.num_test == 1
            else self._decision_history[self.num_test - 2]
        )

        self.wealth -= spend
        self.wealth += self.reward if reward_credited else 0.0
        self._decision_history.append(is_rejected)

        self.last_reject = self.num_test if is_rejected else self.last_reject
        self.wealth_reject = self.wealth if is_rejected else self.wealth_reject

        return is_rejected

class LordPlusPlus(AbstractSequentialTest):
    """Implements LORD++, an improved variant that superseded LORD1 and LORD2.

    LORD++ uses a wealth-based approach where alpha levels are determined by
    accumulated wealth and the gamma sequence. The method tracks rejections
    and spends wealth accordingly.

    References
    ----------
    [1] Ramdas, A., Zrnic, T., Wainwright, M.J. and Jordan, M.I. (2017).
    "SAFFRON: an adaptive algorithm for online control of the false discovery rate."
    arXiv preprint arXiv:1802.09098.

    [2] Javanmard, A., and Montanari, A. (2018).
    "Online Rules for Control of False Discovery Rate and False Discovery Exceedance."
    Annals of Statistics, 46(2):526-554.
    """

    def __init__(self, alpha: float, wealth: float, reward: float | None = None):
        super().__init__(alpha)
        self.alpha0: float = alpha
        self.wealth0: float = wealth
        self.wealth: float = wealth

        validity.check_initial_wealth(wealth, alpha)
        if reward is not None and not math.isclose(reward, alpha):
            raise ValueError(
                "LordPlusPlus guarantee regime requires reward == alpha."
            )
        self.reward: float = alpha

        self.seq = DefaultLordGammaSequence(c=0.07720838)

        self.first_reject: int | None = None  # first rejection index
        self.last_reject: list = []  # rejection indices without first
        self.wealth_at_first_reject: float | None = None

    def test_one(self, p_val: float) -> bool:
        validity.check_p_val(p_val)
        self.num_test += 1

        # Calculate alpha based on LORD++ formula
        self.alpha = self.wealth0 * self.seq.calc_gamma(self.num_test)

        if self.first_reject is not None:
            # Add contribution from first rejection
            self.alpha += (self.alpha0 - self.wealth0) * self.seq.calc_gamma(
                self.num_test - self.first_reject
            )

            # Add contributions from subsequent rejections
            self.alpha += self.alpha0 * sum(
                self.seq.calc_gamma(self.num_test - reject_idx)
                for reject_idx in self.last_reject
            )

        # Ensure we don't spend more than available wealth
        self.alpha = min(self.alpha, self.wealth)

        is_rejected = p_val <= self.alpha

        # Update wealth: spend alpha, gain reward if rejected
        self.wealth -= self.alpha
        if is_rejected:
            self.wealth += self.reward

            if self.first_reject is None:
                # First rejection
                self.first_reject = self.num_test
                self.wealth_at_first_reject = self.wealth
            else:
                # Subsequent rejection
                self.last_reject.append(self.num_test)

        return is_rejected

class LordDependent(AbstractSequentialTest):
    """Implements a variant of LORD for dependent p-values[1]_.

    References
    ----------
    [1] Javanmard, A., and A. Montanari.
    Online rules for control of false discovery rate
    and false discovery exceedance.
    Annals of Statistics, 46(2):526-554, 2018."""

    def __init__(
        self,
        alpha: float,
        wealth: float,
        reward: float,
    ):  # fmt: skip
        super().__init__(alpha)
        self.alpha0: float = alpha
        self.wealth: float = wealth
        self.reward: float = reward
        validity.check_initial_wealth(wealth, alpha)
        validity.check_reward_budget(wealth=wealth, reward=reward, alpha=alpha)

        # onlineFDR LORDdep default xi_i:
        # xi_i = 0.139307 * alpha / (b0 * i * log(max(i,2))^3)
        self.seq = DependentLordGammaSequence(
            c=0.139307 * self.alpha0 / self.reward,
            b0=self.reward,
        )

        self.last_reject: int = 0  # tau
        self.wealth_reject: float = self.wealth  # wealth at tau

    def test_one(self, p_val: float) -> bool:
        validity.check_p_val(p_val)
        self.num_test += 1

        self.alpha = (  # fmt: skip
            self.seq.calc_gamma(self.num_test)  # fmt: skip
            * self.wealth_reject
        )

        is_rejected = p_val <= self.alpha

        self.wealth -= self.alpha
        self.wealth += self.reward if is_rejected else 0
        self.last_reject = self.num_test if is_rejected else self.last_reject
        self.wealth_reject = self.wealth if is_rejected else self.wealth_reject

        return is_rejected

class LordDiscard(AbstractSequentialTest):
    """Implements LORD++ with discarding as
    described in [1]_.

    References
    ----------
    [1] Tian, J., and A. Ramdas.
    ADDIS: an adaptive discarding algorithm for
    online FDR control with conservative nulls.
    In Advances in Neural Information Processing Systems
    (NeurIPS 2019), vol. 32. Curran Associates, Inc., 2019."""

    def __init__(self, alpha: float, wealth: float, tau: float):
        super().__init__(alpha)
        self.alpha0: float = alpha
        self.wealth0: float = wealth
        self.tau: float = tau
        validity.check_initial_wealth(wealth, alpha)
        validity.check_tau(tau)

        self.seq = DefaultLordGammaSequence(c=0.07720838)

        self.first_reject: int | None = None  # first rejection index
        self.last_reject: list = []  # without first rejection

    def _compute_alpha(self, tested_index: int) -> float:
        alpha = self.wealth0 * self.seq.calc_gamma(tested_index)

        if self.first_reject is not None:
            alpha += (
                (self.tau * self.alpha0 - self.wealth0)
                * self.seq.calc_gamma(tested_index - self.first_reject)
            )

        if self.last_reject:
            alpha += (
                self.tau
                * self.alpha0
                * sum(
                    self.seq.calc_gamma(tested_index - reject_idx)
                    for reject_idx in self.last_reject
                )
            )

        return float(alpha)

    def test_one(self, p_val: float) -> bool:
        validity.check_p_val(p_val)
        next_tested_index = self.num_test + 1
        # Expose the same per-step threshold semantics as onlineFDR,
        # including discarded p-values.
        self.alpha = self._compute_alpha(next_tested_index)

        if p_val > self.tau:
            return False  # discard

        self.num_test = next_tested_index

        is_rejected = p_val <= min(self.tau, self.alpha)

        if is_rejected:
            if self.first_reject is None:
                self.first_reject = self.num_test
            else:
                self.last_reject.append(self.num_test)

        return is_rejected

class LORDMemoryDecay(AbstractSequentialTest):
    """LORD variant with memory decay for time series anomaly detection.

    This variant is designed for non-stationary time series where recent
    discoveries are more relevant than older ones. Unlike standard LORD
    variants, it does NOT track wealth. Instead, it uses a decay factor
    to down-weight older rejections and a smoothing parameter to control
    the base detection threshold.

    The algorithm spends:
        alpha_t = alpha * eta * max(gamma(t), 1-delta)
                  + alpha * sum_r decay(t,r) * gamma(t-r-l)

    where the sum is over past rejections r, and decay(t,r) = delta^(t-r-l).

    Key differences from standard LORD:
    - No wealth tracking or accumulation
    - Uses decay to forget old discoveries
    - Ensures minimum spending via max(gamma(t), 1-delta)

    References
    ----------
    [1] Rebjock, Q., B. Kurt, T. Januschowski, and L. Callot.
    Online false discovery rate control for anomaly detection in time series.
    In Advances in Neural Information Processing Systems (NeurIPS 2021),
    vol. 34, pp. 26487-26498. Curran Associates, Inc., 2021.
    """

    def __init__(self, alpha: float, delta: float = 0.99, eta: float = 0.5, l: int = 0):
        """
        Parameters
        ----------
        alpha : float
            Overall significance level in (0, 1). This is the target FDR level.
        delta : float, optional
            Decay factor in (0, 1) for down-weighting older rejections.
            Default is 0.99 (1% decay per time step). Lower values mean
            faster forgetting of old discoveries.
        eta : float, optional
            Smoothing/scaling factor that controls the base detection threshold.
            - eta * alpha * max(gamma(t), 1-delta) is spent at each step
            - eta=0.1: Conservative (10% of budget)
            - eta=0.5: Moderate (50% of budget) [default]
            - eta=1.0: Aggressive (full budget)
        l : int, optional
            Dependency lag parameter. Set l>0 if p-values have serial dependence.
            Default is 0 (assumes independence).
        """
        super().__init__(alpha)
        self.alpha0: float = alpha
        self.delta: float = delta
        self.eta: float = eta
        self.l: int = l

        validity.check_decay_factor(delta)
        if not 0 < eta <= 1:
            raise ValueError(f"eta must be in (0, 1], got {eta}")

        self.seq = DefaultLordGammaSequence(c=0.07720838)

        self.rejection_times: list[int] = []  # all rejection times
        self._gamma_cache: dict[int, float] = {}  # Cache for efficiency

    def test_one(self, p_val: float) -> bool:
        validity.check_p_val(p_val)
        self.num_test += 1

        # Base component with smoothing and minimum threshold
        if self.num_test not in self._gamma_cache:
            self._gamma_cache[self.num_test] = self.seq.calc_gamma(self.num_test)

        gamma_t = self._gamma_cache[self.num_test]
        self.alpha = self.alpha0 * self.eta * max(gamma_t, 1 - self.delta)

        # Add decayed contributions from past rejections
        for reject_idx in self.rejection_times:
            time_diff = self.num_test - reject_idx - self.l
            if time_diff > 0:
                # Cache gamma values for efficiency
                if time_diff not in self._gamma_cache:
                    self._gamma_cache[time_diff] = self.seq.calc_gamma(time_diff)

                decay_weight = self.delta**time_diff
                gamma_val = self._gamma_cache[time_diff]
                self.alpha += self.alpha0 * decay_weight * gamma_val

        is_rejected = p_val <= self.alpha

        if is_rejected:
            self.rejection_times.append(self.num_test)

        return is_rejected