Regressors

online_cp.regressors.ConformalRidgeRegressor

Bases: ConformalRegressor

Conformal ridge regression (Algorithm 2.4 in Algorithmic Learning in a Random World)

Let's create a dataset with noisy evaluations of the function f(x1,x2) = x1+x2:

import numpy as np np.random.seed(31337) # only needed for doctests N = 30 X = np.random.uniform(0, 1, (N, 2)) y = X.sum(axis=1) + np.random.normal(0, 0.1, N)

Import the library and create a regressor:

cp = ConformalRidgeRegressor()

Learn the whole dataset:

cp.learn_initial_training_set(X, y)

Predict an object (output may not be exactly the same, as the dataset depends on the random seed):

interval = cp.predict(np.array([0.5, 0.5]), bounds="both") print("(%.2f, %.2f)" % (interval.lower, interval.upper)) (0.73, 1.23)

You can of course learn a new data point online:

cp.learn_one(np.array([0.5, 0.5]), 1.0)

The prediction set is the closed interval whose boundaries are indicated by the output.

We can then predict again:

interval = cp.predict(np.array([2, 4]), bounds="both") print("(%.2f, %.2f)" % (interval.lower, interval.upper)) (5.39, 6.33)

Source code in src/online_cp/regressors.py

class ConformalRidgeRegressor(ConformalRegressor):
    """
    Conformal ridge regression (Algorithm 2.4 in Algorithmic Learning in a Random World)

    Let's create a dataset with noisy evaluations of the function f(x1,x2) = x1+x2:

    >>> import numpy as np
    >>> np.random.seed(31337)  # only needed for doctests
    >>> N = 30
    >>> X = np.random.uniform(0, 1, (N, 2))
    >>> y = X.sum(axis=1) + np.random.normal(0, 0.1, N)

    Import the library and create a regressor:

    >>> cp = ConformalRidgeRegressor()

    Learn the whole dataset:

    >>> cp.learn_initial_training_set(X, y)

    Predict an object (output may not be exactly the same, as the dataset
    depends on the random seed):
    >>> interval = cp.predict(np.array([0.5, 0.5]), bounds="both")
    >>> print("(%.2f, %.2f)" % (interval.lower, interval.upper))
    (0.73, 1.23)

    You can of course learn a new data point online:

    >>> cp.learn_one(np.array([0.5, 0.5]), 1.0)

    The prediction set is the closed interval whose boundaries are indicated by the output.

    We can then predict again:

    >>> interval = cp.predict(np.array([2, 4]), bounds="both")
    >>> print("(%.2f, %.2f)" % (interval.lower, interval.upper))
    (5.39, 6.33)
    """

    # TODO: Fix gracefull error handling when the matrix is singular. It should raise an exception, but we could
    #       specify that it can be handled by changing the ridge parameter.

    def __init__(
        self, a=0, warnings=True, autotune=False, verbose=0, rnd_state=None, studentised=False, epsilon=default_epsilon
    ) -> None:
        """
        The ridge parameter (L2 regularisation) is a.
        Setting autotune=True automatically tunes the ridge parameter using generalized cross validation when learning initial training set.
        """
        super().__init__(epsilon=epsilon)

        self.a = a
        self.X = None
        self.y = None
        self.p = None
        self.Id = None
        self.XTXinv = None

        # Should we raise warnings
        self.warnings = warnings
        # Do we autotune ridge prarmeter on warning
        self.autotune = autotune

        self.verbose = verbose
        self.rnd_gen = np.random.default_rng(rnd_state)

        # Do we use the studentised residuals
        self.studentised = studentised

    def learn_initial_training_set(self, X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None:
        self.X = X
        self.y = y
        self.p = X.shape[1]
        self.Id = np.identity(self.p)
        if self.autotune:
            self._tune_ridge_parameter()
        else:
            self.XTXinv = np.linalg.inv(self.X.T @ self.X + self.a * self.Id)

    def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
        """
        Learn a single example. If we have already computed X and XTXinv, use them for update. Then the last row of X is the object with label y.
        >>> cp = ConformalRidgeRegressor()
        >>> cp.learn_one(np.array([1, 0]), 1)
        >>> cp.X
        array([[1, 0]])
        >>> cp.y
        array([1])
        """
        # Learn label y
        if self.y is None:
            self.y = np.array([y])
        else:
            if hasattr(self, "h"):
                self.y = np.append(self.y, y.reshape(1, self.h), axis=0)
            else:
                self.y = np.append(self.y, y)

        if precomputed is not None:
            X = precomputed["X"]
            XTXinv = precomputed["XTXinv"]

            if X is not None:
                self.X = X
                self.p = self.X.shape[1]
                self.Id = np.identity(self.p)

            if XTXinv is not None:
                self.XTXinv = XTXinv

            else:
                if self.X.shape[0] == 1:
                    # print(self.X)
                    # print(self.Id)
                    # print(self.a)
                    self.XTXinv = np.linalg.inv(self.X.T @ self.X + self.a * self.Id)
                else:
                    # Update XTX_inv (inverse of Kernel matrix plus regularisation) Use the Sherman-Morrison formula to update the hat matrix
                    # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
                    self.XTXinv -= (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)

                    # Check the rank
                    if self.warnings:
                        self.check_matrix_rank(self.XTXinv)

        else:
            # Learn object x
            if self.X is None:
                self.X = x.reshape(1, -1)
                self.p = self.X.shape[1]
                self.Id = np.identity(self.p)
            elif self.X.shape[0] == 1:
                self.X = np.append(self.X, x.reshape(1, -1), axis=0)
                self.XTXinv = np.linalg.inv(self.X.T @ self.X + self.a * self.Id)
            else:
                self.X = np.append(self.X, x.reshape(1, -1), axis=0)
                # Update XTX_inv (inverse of Kernel matrix plus regularisation) Use the Sherman-Morrison formula to update the hat matrix
                # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
                self.XTXinv -= (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)

                # Check the rank
                if self.warnings:
                    self.check_matrix_rank(self.XTXinv)

    def compute_A_and_B_OLD(self, X, XTXinv, y):
        n = X.shape[0]
        # Hat matrix (This block is the time consuming one...)
        H = X @ XTXinv @ X.T
        C = np.identity(n) - H
        A = C @ np.append(y, 0)  # Elements of this vector are denoted ai
        B = C @ np.append(np.zeros((n - 1,)), 1)  # Elements of this vector are denoted bi
        if self.studentised:
            h = H.diagonal()
            A = A / np.sqrt(1 - h)
            B = B / np.sqrt(1 - h)
        # Nonconformity scores are A + yB = y - yhat
        return A, B

    def compute_A_and_B(self, X, XTXinv, y):
        """
        Efficient and correct computation of A and B for conformal ridge regression.
        X: (n, d) augmented matrix (last row is test object)
        XTXinv: (d, d) inverse of X.T @ X + a*I for augmented X
        y: (n-1,) training labels (no test label)
        """
        y_ext = np.append(y, 0)  # y with test point (last row) as 0

        # Compute beta using the augmented X and y_ext (just like the old code)
        beta = XTXinv @ X.T @ y_ext  # (d, d) @ (d, n) @ (n,) -> (d,)

        # Fitted values for all points (including test)
        y_hat = X @ beta  # (n, d) @ (d,) -> (n,)

        # Compute hat matrix diagonal for all points using XTXinv (augmented)
        H_diag = np.sum(X @ XTXinv * X, axis=1)  # (n,)

        # Compute last column of H efficiently
        h_col = X @ XTXinv @ X[-1]  # (n, d) @ (d,) -> (n,)

        # A and B for each point
        A = y_ext - y_hat
        B = -h_col
        B[-1] += 1  # e_{-1}[-1] = 1

        if self.studentised:
            A = A / np.sqrt(1 - H_diag + 1e-12)
            B = B / np.sqrt(1 - H_diag + 1e-12)

        return A, B

    def predict(
        self,
        x: NDArray[np.floating[Any]],
        epsilon: float | NDArray[np.floating[Any]] | None = None,
        bounds: str = "both",
        return_update: bool = False,
        debug_time: bool = False,
    ) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
        """
        This function makes a prediction.

        If you start with no training,
        you get a null prediciton between
        -infinity and +infinity.

        >>> cp = ConformalRidgeRegressor()
        >>> cp.predict(np.array([0.506, 0.22, -0.45]), bounds="both")
        (-inf, inf)
        """

        def build_precomputed(X, XTXinv, A, B):
            computed = {
                "X": X,  # The updated matrix of objects
                "XTXinv": XTXinv,  # The updated kernel matrix
                "A": A,
                "B": B,
            }
            return computed

        if epsilon is None:
            epsilon = self.epsilon

        if self._safe_size_check(self.X) > 0:
            tic = time.time()
            # Add row to X matrix
            X = np.append(self.X, x.reshape(1, -1), axis=0)
            toc_add_row = time.time() - tic
            n = X.shape[0]
            XTXinv = None

            # Check that the significance level is not too small. If it is, return infinite prediction interval
            # For multi-level: only bail out if even the largest epsilon is too small
            eps_check = max(epsilon) if hasattr(epsilon, '__iter__') else epsilon
            if bounds == "both":
                if not (eps_check >= 2 / n):
                    if self.warnings:
                        eps_warn = min(epsilon) if hasattr(epsilon, '__iter__') else epsilon
                        warnings.warn(
                            f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(2 / eps_warn))} examples. Increase or add more examples",
                            stacklevel=2,
                        )
                    if hasattr(epsilon, '__iter__'):
                        predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                        result = MultiLevelPredictionInterval(predictions)
                    else:
                        result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                    if return_update:
                        return result, build_precomputed(X, XTXinv, None, None)
                    else:
                        return result
            else:
                if not (eps_check >= 1 / n):
                    if self.warnings:
                        eps_warn = min(epsilon) if hasattr(epsilon, '__iter__') else epsilon
                        warnings.warn(
                            f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(1 / eps_warn))} examples. Increase or add more examples",
                            stacklevel=2,
                        )
                    if hasattr(epsilon, '__iter__'):
                        predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                        result = MultiLevelPredictionInterval(predictions)
                    else:
                        result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                    if return_update:
                        return result, build_precomputed(X, XTXinv, None, None)
                    else:
                        return result

            tic = time.time()
            # Update XTX_inv (inverse of Kernel matrix plus regularisation) Use the Sherman-Morrison formula to update the hat matrix
            # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula

            XTXinv = self.XTXinv - (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)
            toc_update_XTXinv = time.time() - tic

            tic = time.time()
            A, B = self.compute_A_and_B(X, XTXinv, self.y)
            toc_nc = time.time() - tic

            if self.studentised:
                tic = time.time()
                t = (A[:-1] - A[-1]) / (B[-1] - B[:-1])
                t.sort()
                l_dic = {i + 1: val for i, val in enumerate(t)}
                u_dic = {i + 1: val for i, val in enumerate(t)}
                toc_dics = time.time() - tic
            else:
                tic = time.time()
                l_dic, u_dic = self._vectorised_l_and_u(A, B)
                toc_dics = time.time() - tic

            if bounds == "both":
                if hasattr(epsilon, '__iter__'):
                    predictions = {}
                    for eps in epsilon:
                        lo = self._get_lower(l_dic=l_dic, epsilon=eps / 2, n=n)
                        up = self._get_upper(u_dic=u_dic, epsilon=eps / 2, n=n)
                        predictions[eps] = self._construct_Gamma(lo, up, eps)
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    lower = self._get_lower(l_dic=l_dic, epsilon=epsilon / 2, n=n)
                    upper = self._get_upper(u_dic=u_dic, epsilon=epsilon / 2, n=n)
                    result = self._construct_Gamma(lower, upper, epsilon)
            elif bounds == "lower":
                if hasattr(epsilon, '__iter__'):
                    predictions = {}
                    for eps in epsilon:
                        lo = self._get_lower(l_dic=l_dic, epsilon=eps, n=n)
                        predictions[eps] = self._construct_Gamma(lo, np.inf, eps)
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    lower = self._get_lower(l_dic=l_dic, epsilon=epsilon, n=n)
                    result = self._construct_Gamma(lower, np.inf, epsilon)
            elif bounds == "upper":
                if hasattr(epsilon, '__iter__'):
                    predictions = {}
                    for eps in epsilon:
                        up = self._get_upper(u_dic=u_dic, epsilon=eps, n=n)
                        predictions[eps] = self._construct_Gamma(-np.inf, up, eps)
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    upper = self._get_upper(u_dic=u_dic, epsilon=epsilon, n=n)
                    result = self._construct_Gamma(-np.inf, upper, epsilon)
            else:
                raise Exception

            if debug_time:
                print(f"Add row: {toc_add_row}")
                print(f"Update kernel: {toc_update_XTXinv}")
                print(f"NC scores: {toc_nc}")
                print(f"l and u: {toc_dics}")
                print()
        else:
            # With just one object, and no label, we cannot predict any meaningful interval
            X = x.reshape(1, -1)
            XTXinv = None
            A = None
            B = None

            if hasattr(epsilon, '__iter__'):
                predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                result = MultiLevelPredictionInterval(predictions)
            else:
                result = self._construct_Gamma(-np.inf, np.inf, epsilon)

        if return_update:
            return result, build_precomputed(X, XTXinv, A, B)
        else:
            return result

    def compute_p_value(self, x, y, bounds="both", precomputed=None, tau=None, smoothed=True):
        """
        Computes the smoothed p-value of the example (x, y).
        """
        if tau is None and smoothed:
            tau = self.rnd_gen.uniform(0, 1)
        if precomputed is not None:
            assert np.allclose(x, precomputed["X"][-1])
            A = precomputed["A"]
            B = precomputed["B"]
        else:
            if self.XTXinv is not None:
                X = np.append(self.X, x.reshape(1, -1), axis=0)
                XTXinv = self.XTXinv - (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)
                A, B = self.compute_A_and_B(X, XTXinv, self.y)
            else:
                A, B = None, None

        if A is not None and B is not None:
            if bounds == "both":
                E = A + y * B
                Alpha = np.zeros_like(A)
                for i, e in enumerate(E):
                    alpha = min((E >= e).sum(), (E <= e).sum())
                    Alpha[i] = alpha
                c_type = "conformity"
            elif bounds == "lower":
                Alpha = -(A + y * B)
                c_type = "nonconformity"
            elif bounds == "upper":
                Alpha = A + y * B
                c_type = "nonconformity"
            else:
                raise Exception('bounds must be one of "both", "lower", "upper"')

            if smoothed:
                p = self._compute_p_value(Alpha, tau, c_type=c_type)
            else:
                p = self._compute_p_value(Alpha, c_type=c_type)
        else:
            if smoothed:
                p = tau
            else:
                p = 1

        return p

    def change_ridge_parameter(self, a):
        """
        Change the ridge parameter
        >>> cp = ConformalRidgeRegressor()
        >>> cp.learn_one(np.array([1, 0]), 1)
        >>> cp.change_ridge_parameter(1)
        >>> cp.a
        1
        """
        self.a = a
        if self.X is not None:
            self.XTXinv = np.linalg.inv(self.X.T @ self.X + self.a * self.Id)

    def _tune_ridge_parameter(self, a0=None):
        """
        Tune ridge parameter with Generalized cross validation https://pages.stat.wisc.edu/~wahba/stat860public/pdf1/golub.heath.wahba.pdf
        """
        XTX = self.X.T @ self.X
        n = self.X.shape[0]
        In = np.identity(n)

        def GCV(a):
            try:
                A = self.X @ np.linalg.inv(XTX + a * self.Id) @ self.X.T
                max_diag_H = np.max(np.diag(A))
                if max_diag_H > 1:
                    return np.inf
                return (1 / n) * np.linalg.norm((In - A) @ self.y) ** 2 / ((1 / n) * np.trace(In - A)) ** 2
            except (np.linalg.LinAlgError, ZeroDivisionError):
                return np.inf

        # Initial guess
        if a0 is None:
            a0 = 1e-6  # Just a small pertubation to avoid numerical issues

        # Bounds to ensure a >= 0
        res = minimize(
            GCV, x0=a0, bounds=Bounds(lb=1e-6, keep_feasible=True)
        )  # May be relevant to pass some arguments here, or even use another minimizer.
        a = res.x[0]

        if self.verbose > 0:
            print(f"New ridge parameter: {a}")
        self.change_ridge_parameter(a)

    # TODO
    def prune_training_set(self):
        """
        Just an idea at the moment, but perhaps we should have some inclusion criteria for examples to only include the informative ones. Could improve accuracy, but also significantly decrease computation time if we have a large dataset.
        """
        raise NotImplementedError

    def check_matrix_rank(self, M):
        """
        Check if a matrix has full rank <==> is invertible
        Returns False if matrix is rank deficient
        NOTE In numerical linear algebra it is a bit more subtle. The condition number can tell us more.

        >>> cp = ConformalRidgeRegressor(warnings=False)
        >>> cp.check_matrix_rank(np.array([[1, 0], [1, 0]]))
        False
        >>> cp.check_matrix_rank(np.array([[1, 0], [0, 1]]))
        True
        """
        if np.linalg.matrix_rank(M) < M.shape[0]:
            if self.warnings:
                warnings.warn(
                    f"The matrix X is rank deficient. Condition number: {np.linalg.cond(M)}. Consider changing the ridge parameter",
                    stacklevel=2,
                )
            return False
        else:
            return True

__init__(a=0, warnings=True, autotune=False, verbose=0, rnd_state=None, studentised=False, epsilon=default_epsilon) -> None

The ridge parameter (L2 regularisation) is a. Setting autotune=True automatically tunes the ridge parameter using generalized cross validation when learning initial training set.

Source code in src/online_cp/regressors.py

def __init__(
    self, a=0, warnings=True, autotune=False, verbose=0, rnd_state=None, studentised=False, epsilon=default_epsilon
) -> None:
    """
    The ridge parameter (L2 regularisation) is a.
    Setting autotune=True automatically tunes the ridge parameter using generalized cross validation when learning initial training set.
    """
    super().__init__(epsilon=epsilon)

    self.a = a
    self.X = None
    self.y = None
    self.p = None
    self.Id = None
    self.XTXinv = None

    # Should we raise warnings
    self.warnings = warnings
    # Do we autotune ridge prarmeter on warning
    self.autotune = autotune

    self.verbose = verbose
    self.rnd_gen = np.random.default_rng(rnd_state)

    # Do we use the studentised residuals
    self.studentised = studentised

learn_one(x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None

Learn a single example. If we have already computed X and XTXinv, use them for update. Then the last row of X is the object with label y.

cp = ConformalRidgeRegressor() cp.learn_one(np.array([1, 0]), 1) cp.X array([[1, 0]]) cp.y array([1])

Source code in src/online_cp/regressors.py

def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
    """
    Learn a single example. If we have already computed X and XTXinv, use them for update. Then the last row of X is the object with label y.
    >>> cp = ConformalRidgeRegressor()
    >>> cp.learn_one(np.array([1, 0]), 1)
    >>> cp.X
    array([[1, 0]])
    >>> cp.y
    array([1])
    """
    # Learn label y
    if self.y is None:
        self.y = np.array([y])
    else:
        if hasattr(self, "h"):
            self.y = np.append(self.y, y.reshape(1, self.h), axis=0)
        else:
            self.y = np.append(self.y, y)

    if precomputed is not None:
        X = precomputed["X"]
        XTXinv = precomputed["XTXinv"]

        if X is not None:
            self.X = X
            self.p = self.X.shape[1]
            self.Id = np.identity(self.p)

        if XTXinv is not None:
            self.XTXinv = XTXinv

        else:
            if self.X.shape[0] == 1:
                # print(self.X)
                # print(self.Id)
                # print(self.a)
                self.XTXinv = np.linalg.inv(self.X.T @ self.X + self.a * self.Id)
            else:
                # Update XTX_inv (inverse of Kernel matrix plus regularisation) Use the Sherman-Morrison formula to update the hat matrix
                # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
                self.XTXinv -= (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)

                # Check the rank
                if self.warnings:
                    self.check_matrix_rank(self.XTXinv)

    else:
        # Learn object x
        if self.X is None:
            self.X = x.reshape(1, -1)
            self.p = self.X.shape[1]
            self.Id = np.identity(self.p)
        elif self.X.shape[0] == 1:
            self.X = np.append(self.X, x.reshape(1, -1), axis=0)
            self.XTXinv = np.linalg.inv(self.X.T @ self.X + self.a * self.Id)
        else:
            self.X = np.append(self.X, x.reshape(1, -1), axis=0)
            # Update XTX_inv (inverse of Kernel matrix plus regularisation) Use the Sherman-Morrison formula to update the hat matrix
            # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
            self.XTXinv -= (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)

            # Check the rank
            if self.warnings:
                self.check_matrix_rank(self.XTXinv)

compute_A_and_B(X, XTXinv, y)

Efficient and correct computation of A and B for conformal ridge regression. X: (n, d) augmented matrix (last row is test object) XTXinv: (d, d) inverse of X.T @ X + a*I for augmented X y: (n-1,) training labels (no test label)

Source code in src/online_cp/regressors.py

def compute_A_and_B(self, X, XTXinv, y):
    """
    Efficient and correct computation of A and B for conformal ridge regression.
    X: (n, d) augmented matrix (last row is test object)
    XTXinv: (d, d) inverse of X.T @ X + a*I for augmented X
    y: (n-1,) training labels (no test label)
    """
    y_ext = np.append(y, 0)  # y with test point (last row) as 0

    # Compute beta using the augmented X and y_ext (just like the old code)
    beta = XTXinv @ X.T @ y_ext  # (d, d) @ (d, n) @ (n,) -> (d,)

    # Fitted values for all points (including test)
    y_hat = X @ beta  # (n, d) @ (d,) -> (n,)

    # Compute hat matrix diagonal for all points using XTXinv (augmented)
    H_diag = np.sum(X @ XTXinv * X, axis=1)  # (n,)

    # Compute last column of H efficiently
    h_col = X @ XTXinv @ X[-1]  # (n, d) @ (d,) -> (n,)

    # A and B for each point
    A = y_ext - y_hat
    B = -h_col
    B[-1] += 1  # e_{-1}[-1] = 1

    if self.studentised:
        A = A / np.sqrt(1 - H_diag + 1e-12)
        B = B / np.sqrt(1 - H_diag + 1e-12)

    return A, B

predict(x: NDArray[np.floating[Any]], epsilon: float | NDArray[np.floating[Any]] | None = None, bounds: str = 'both', return_update: bool = False, debug_time: bool = False) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]

This function makes a prediction.

If you start with no training, you get a null prediciton between -infinity and +infinity.

cp = ConformalRidgeRegressor() cp.predict(np.array([0.506, 0.22, -0.45]), bounds="both") (-inf, inf)

Source code in src/online_cp/regressors.py

def predict(
    self,
    x: NDArray[np.floating[Any]],
    epsilon: float | NDArray[np.floating[Any]] | None = None,
    bounds: str = "both",
    return_update: bool = False,
    debug_time: bool = False,
) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
    """
    This function makes a prediction.

    If you start with no training,
    you get a null prediciton between
    -infinity and +infinity.

    >>> cp = ConformalRidgeRegressor()
    >>> cp.predict(np.array([0.506, 0.22, -0.45]), bounds="both")
    (-inf, inf)
    """

    def build_precomputed(X, XTXinv, A, B):
        computed = {
            "X": X,  # The updated matrix of objects
            "XTXinv": XTXinv,  # The updated kernel matrix
            "A": A,
            "B": B,
        }
        return computed

    if epsilon is None:
        epsilon = self.epsilon

    if self._safe_size_check(self.X) > 0:
        tic = time.time()
        # Add row to X matrix
        X = np.append(self.X, x.reshape(1, -1), axis=0)
        toc_add_row = time.time() - tic
        n = X.shape[0]
        XTXinv = None

        # Check that the significance level is not too small. If it is, return infinite prediction interval
        # For multi-level: only bail out if even the largest epsilon is too small
        eps_check = max(epsilon) if hasattr(epsilon, '__iter__') else epsilon
        if bounds == "both":
            if not (eps_check >= 2 / n):
                if self.warnings:
                    eps_warn = min(epsilon) if hasattr(epsilon, '__iter__') else epsilon
                    warnings.warn(
                        f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(2 / eps_warn))} examples. Increase or add more examples",
                        stacklevel=2,
                    )
                if hasattr(epsilon, '__iter__'):
                    predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                if return_update:
                    return result, build_precomputed(X, XTXinv, None, None)
                else:
                    return result
        else:
            if not (eps_check >= 1 / n):
                if self.warnings:
                    eps_warn = min(epsilon) if hasattr(epsilon, '__iter__') else epsilon
                    warnings.warn(
                        f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(1 / eps_warn))} examples. Increase or add more examples",
                        stacklevel=2,
                    )
                if hasattr(epsilon, '__iter__'):
                    predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                if return_update:
                    return result, build_precomputed(X, XTXinv, None, None)
                else:
                    return result

        tic = time.time()
        # Update XTX_inv (inverse of Kernel matrix plus regularisation) Use the Sherman-Morrison formula to update the hat matrix
        # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula

        XTXinv = self.XTXinv - (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)
        toc_update_XTXinv = time.time() - tic

        tic = time.time()
        A, B = self.compute_A_and_B(X, XTXinv, self.y)
        toc_nc = time.time() - tic

        if self.studentised:
            tic = time.time()
            t = (A[:-1] - A[-1]) / (B[-1] - B[:-1])
            t.sort()
            l_dic = {i + 1: val for i, val in enumerate(t)}
            u_dic = {i + 1: val for i, val in enumerate(t)}
            toc_dics = time.time() - tic
        else:
            tic = time.time()
            l_dic, u_dic = self._vectorised_l_and_u(A, B)
            toc_dics = time.time() - tic

        if bounds == "both":
            if hasattr(epsilon, '__iter__'):
                predictions = {}
                for eps in epsilon:
                    lo = self._get_lower(l_dic=l_dic, epsilon=eps / 2, n=n)
                    up = self._get_upper(u_dic=u_dic, epsilon=eps / 2, n=n)
                    predictions[eps] = self._construct_Gamma(lo, up, eps)
                result = MultiLevelPredictionInterval(predictions)
            else:
                lower = self._get_lower(l_dic=l_dic, epsilon=epsilon / 2, n=n)
                upper = self._get_upper(u_dic=u_dic, epsilon=epsilon / 2, n=n)
                result = self._construct_Gamma(lower, upper, epsilon)
        elif bounds == "lower":
            if hasattr(epsilon, '__iter__'):
                predictions = {}
                for eps in epsilon:
                    lo = self._get_lower(l_dic=l_dic, epsilon=eps, n=n)
                    predictions[eps] = self._construct_Gamma(lo, np.inf, eps)
                result = MultiLevelPredictionInterval(predictions)
            else:
                lower = self._get_lower(l_dic=l_dic, epsilon=epsilon, n=n)
                result = self._construct_Gamma(lower, np.inf, epsilon)
        elif bounds == "upper":
            if hasattr(epsilon, '__iter__'):
                predictions = {}
                for eps in epsilon:
                    up = self._get_upper(u_dic=u_dic, epsilon=eps, n=n)
                    predictions[eps] = self._construct_Gamma(-np.inf, up, eps)
                result = MultiLevelPredictionInterval(predictions)
            else:
                upper = self._get_upper(u_dic=u_dic, epsilon=epsilon, n=n)
                result = self._construct_Gamma(-np.inf, upper, epsilon)
        else:
            raise Exception

        if debug_time:
            print(f"Add row: {toc_add_row}")
            print(f"Update kernel: {toc_update_XTXinv}")
            print(f"NC scores: {toc_nc}")
            print(f"l and u: {toc_dics}")
            print()
    else:
        # With just one object, and no label, we cannot predict any meaningful interval
        X = x.reshape(1, -1)
        XTXinv = None
        A = None
        B = None

        if hasattr(epsilon, '__iter__'):
            predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
            result = MultiLevelPredictionInterval(predictions)
        else:
            result = self._construct_Gamma(-np.inf, np.inf, epsilon)

    if return_update:
        return result, build_precomputed(X, XTXinv, A, B)
    else:
        return result

compute_p_value(x, y, bounds='both', precomputed=None, tau=None, smoothed=True)

Computes the smoothed p-value of the example (x, y).

Source code in src/online_cp/regressors.py

def compute_p_value(self, x, y, bounds="both", precomputed=None, tau=None, smoothed=True):
    """
    Computes the smoothed p-value of the example (x, y).
    """
    if tau is None and smoothed:
        tau = self.rnd_gen.uniform(0, 1)
    if precomputed is not None:
        assert np.allclose(x, precomputed["X"][-1])
        A = precomputed["A"]
        B = precomputed["B"]
    else:
        if self.XTXinv is not None:
            X = np.append(self.X, x.reshape(1, -1), axis=0)
            XTXinv = self.XTXinv - (self.XTXinv @ np.outer(x, x) @ self.XTXinv) / (1 + x.T @ self.XTXinv @ x)
            A, B = self.compute_A_and_B(X, XTXinv, self.y)
        else:
            A, B = None, None

    if A is not None and B is not None:
        if bounds == "both":
            E = A + y * B
            Alpha = np.zeros_like(A)
            for i, e in enumerate(E):
                alpha = min((E >= e).sum(), (E <= e).sum())
                Alpha[i] = alpha
            c_type = "conformity"
        elif bounds == "lower":
            Alpha = -(A + y * B)
            c_type = "nonconformity"
        elif bounds == "upper":
            Alpha = A + y * B
            c_type = "nonconformity"
        else:
            raise Exception('bounds must be one of "both", "lower", "upper"')

        if smoothed:
            p = self._compute_p_value(Alpha, tau, c_type=c_type)
        else:
            p = self._compute_p_value(Alpha, c_type=c_type)
    else:
        if smoothed:
            p = tau
        else:
            p = 1

    return p

change_ridge_parameter(a)

Change the ridge parameter

cp = ConformalRidgeRegressor() cp.learn_one(np.array([1, 0]), 1) cp.change_ridge_parameter(1) cp.a 1

Source code in src/online_cp/regressors.py

def change_ridge_parameter(self, a):
    """
    Change the ridge parameter
    >>> cp = ConformalRidgeRegressor()
    >>> cp.learn_one(np.array([1, 0]), 1)
    >>> cp.change_ridge_parameter(1)
    >>> cp.a
    1
    """
    self.a = a
    if self.X is not None:
        self.XTXinv = np.linalg.inv(self.X.T @ self.X + self.a * self.Id)

prune_training_set()

Just an idea at the moment, but perhaps we should have some inclusion criteria for examples to only include the informative ones. Could improve accuracy, but also significantly decrease computation time if we have a large dataset.

Source code in src/online_cp/regressors.py

def prune_training_set(self):
    """
    Just an idea at the moment, but perhaps we should have some inclusion criteria for examples to only include the informative ones. Could improve accuracy, but also significantly decrease computation time if we have a large dataset.
    """
    raise NotImplementedError

check_matrix_rank(M)

Check if a matrix has full rank <==> is invertible Returns False if matrix is rank deficient NOTE In numerical linear algebra it is a bit more subtle. The condition number can tell us more.

cp = ConformalRidgeRegressor(warnings=False) cp.check_matrix_rank(np.array([[1, 0], [1, 0]])) False cp.check_matrix_rank(np.array([[1, 0], [0, 1]])) True

Source code in src/online_cp/regressors.py

def check_matrix_rank(self, M):
    """
    Check if a matrix has full rank <==> is invertible
    Returns False if matrix is rank deficient
    NOTE In numerical linear algebra it is a bit more subtle. The condition number can tell us more.

    >>> cp = ConformalRidgeRegressor(warnings=False)
    >>> cp.check_matrix_rank(np.array([[1, 0], [1, 0]]))
    False
    >>> cp.check_matrix_rank(np.array([[1, 0], [0, 1]]))
    True
    """
    if np.linalg.matrix_rank(M) < M.shape[0]:
        if self.warnings:
            warnings.warn(
                f"The matrix X is rank deficient. Condition number: {np.linalg.cond(M)}. Consider changing the ridge parameter",
                stacklevel=2,
            )
        return False
    else:
        return True

online_cp.regressors.ConformalNearestNeighboursRegressor

Bases: ConformalRegressor

Conformal k-nearest neighbours regressor (ALRW2 §2.4).

Produces prediction intervals with guaranteed coverage using the nonconformity measure |y - ŷ_kNN| where ŷ_kNN is the leave-one-out k-nearest neighbours prediction.

import numpy as np np.random.seed(42) N = 30 X = np.random.uniform(0, 1, (N, 1)) y = np.sin(2 * np.pi * X[:, 0]) + np.random.normal(0, 0.1, N) cp = ConformalNearestNeighboursRegressor(k=3) cp.learn_initial_training_set(X, y) interval = cp.predict(np.array([0.25])) bool(interval.lower < np.sin(2 * np.pi * 0.25) < interval.upper) True

Source code in src/online_cp/regressors.py

class ConformalNearestNeighboursRegressor(ConformalRegressor):
    """
    Conformal k-nearest neighbours regressor (ALRW2 §2.4).

    Produces prediction intervals with guaranteed coverage using the
    nonconformity measure |y - ŷ_kNN| where ŷ_kNN is the leave-one-out
    k-nearest neighbours prediction.

    >>> import numpy as np
    >>> np.random.seed(42)
    >>> N = 30
    >>> X = np.random.uniform(0, 1, (N, 1))
    >>> y = np.sin(2 * np.pi * X[:, 0]) + np.random.normal(0, 0.1, N)
    >>> cp = ConformalNearestNeighboursRegressor(k=3)
    >>> cp.learn_initial_training_set(X, y)
    >>> interval = cp.predict(np.array([0.25]))
    >>> bool(interval.lower < np.sin(2 * np.pi * 0.25) < interval.upper)
    True
    """

    def __init__(
        self,
        k=1,
        distance="euclidean",
        distance_func=None,
        aggregation="mean",
        verbose=0,
        rnd_state=None,
        epsilon=default_epsilon,
    ):
        """
        Parameters
        ----------
        k : int
            Number of nearest neighbours.
        distance : str
            Distance metric (passed to scipy.spatial.distance).
        distance_func : callable, optional
            Custom distance function. If provided, `distance` is ignored.
            Signature: distance_func(X, y=None) where X is (n, d) and y is
            (m, d) or None. Returns (n, n) if y is None, else (n, m).
        aggregation : {'mean', 'median'}
            How to aggregate k nearest neighbour labels.
        epsilon : float
            Default significance level.
        """
        super().__init__(epsilon=epsilon)
        self.k = k
        if aggregation not in ("mean", "median"):
            raise ValueError(f"aggregation must be 'mean' or 'median', got '{aggregation}'")
        self.aggregation = aggregation
        self.distance = distance
        if distance_func is None:
            self.distance_func = self._standard_distance_func
        else:
            self.distance_func = distance_func
            self.distance = "custom"
        self.X = None
        self.y = None
        self.D = None
        self.verbose = verbose
        self.rnd_gen = np.random.default_rng(rnd_state)

    def _standard_distance_func(self, X, y=None):
        X = np.atleast_2d(X)
        if y is None:
            return squareform(pdist(X, metric=self.distance))
        else:
            y = np.atleast_2d(y)
            return cdist(X, y, metric=self.distance)

    @staticmethod
    def _update_distance_matrix(D, d):
        """Extend (n, n) distance matrix to (n+1, n+1) given new distances d."""
        d = np.asarray(d).ravel()
        n = D.shape[0]
        D_new = np.empty((n + 1, n + 1), dtype=np.result_type(D.dtype, d.dtype))
        D_new[:n, :n] = D
        D_new[:n, n] = d
        D_new[n, :n] = d
        D_new[n, n] = 0.0
        return D_new

    def _knn_predictions(self, D, y):
        """
        Compute leave-one-out k-NN predictions for all points.

        For point i, the prediction is the aggregation of the labels of the
        k nearest neighbours of i (excluding i itself).

        Parameters
        ----------
        D : ndarray, shape (n, n)
            Distance matrix.
        y : ndarray, shape (n,)
            Labels.

        Returns
        -------
        y_hat : ndarray, shape (n,)
            Leave-one-out k-NN predictions.
        """
        n = D.shape[0]
        k = min(self.k, n - 1)  # can't have more neighbours than n-1
        if k == 0:
            return np.zeros(n)

        agg_func = np.mean if self.aggregation == "mean" else np.median

        # Set diagonal to inf so a point is never its own neighbour
        D_work = D.copy()
        np.fill_diagonal(D_work, np.inf)

        # Find k nearest neighbour indices for all points at once
        # np.argpartition is O(n) per row
        knn_idx = np.argpartition(D_work, k - 1, axis=1)[:, :k]

        # Gather neighbour labels and aggregate
        y_neighbours = y[knn_idx]  # shape (n, k)
        y_hat = agg_func(y_neighbours, axis=1)
        return y_hat

    def learn_initial_training_set(self, X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None:
        """Batch-initialize with training data.

        Parameters
        ----------
        X : ndarray, shape (n, d)
            Training objects.
        y : ndarray, shape (n,)
            Training labels.
        """
        self.X = np.atleast_2d(X)
        self.y = np.asarray(y, dtype=float)
        self.D = self.distance_func(self.X)

    def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
        """Incrementally add one observation.

        Parameters
        ----------
        x : array-like, shape (d,)
            New object.
        y : float
            New label.
        precomputed : dict, optional
            If provided, should contain 'D' (the already-updated distance matrix).
        """
        x = np.asarray(x).ravel()

        if self.X is None:
            self.X = x.reshape(1, -1)
            self.y = np.array([y], dtype=float)
            self.D = self.distance_func(self.X)
        else:
            if precomputed is not None and "D" in precomputed:
                self.D = precomputed["D"]
            else:
                d = self.distance_func(self.X, x.reshape(1, -1)).ravel()
                self.D = self._update_distance_matrix(self.D, d)
            self.X = np.vstack([self.X, x.reshape(1, -1)])
            self.y = np.append(self.y, y)

    def predict(
        self,
        x: NDArray[np.floating[Any]],
        epsilon: float | NDArray[np.floating[Any]] | None = None,
        bounds: str = "both",
        return_update: bool = False,
    ) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
        """Predict a conformal prediction interval for test object x.

        Parameters
        ----------
        x : array-like, shape (d,)
            Test object.
        epsilon : float or array-like, optional
            Significance level(s). If None, uses self.epsilon.
        bounds : {'both', 'lower', 'upper'}
            Which bounds to compute.
        return_update : bool
            If True, also return precomputed dict for subsequent learn_one.

        Returns
        -------
        result : ConformalPredictionInterval or MultiLevelPredictionInterval
        precomputed : dict, optional
            Returned if return_update=True. Contains 'D'.
        """
        x = np.asarray(x).ravel()

        if epsilon is None:
            epsilon = self.epsilon

        n = self._safe_size_check(self.X)

        if n == 0:
            # No training data
            if hasattr(epsilon, '__iter__'):
                predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                result = MultiLevelPredictionInterval(predictions)
            else:
                result = self._construct_Gamma(-np.inf, np.inf, epsilon)
            if return_update:
                return result, {"D": None}
            return result

        # Augment distance matrix with test point
        d = self.distance_func(self.X, x.reshape(1, -1)).ravel()
        D_aug = self._update_distance_matrix(self.D, d)
        n_aug = n + 1  # augmented size

        # Check significance level is feasible
        eps_check = max(epsilon) if hasattr(epsilon, '__iter__') else epsilon
        min_needed = 2 if bounds == "both" else 1
        if not (eps_check >= min_needed / n_aug):
            if hasattr(epsilon, '__iter__'):
                predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                result = MultiLevelPredictionInterval(predictions)
            else:
                result = self._construct_Gamma(-np.inf, np.inf, epsilon)
            if return_update:
                return result, {"D": D_aug}
            return result

        # Compute leave-one-out k-NN predictions for all training points
        # in the augmented set. The test point's label is unknown, so we
        # need to determine the interval of y values that would be included.
        #
        # Key insight: For the n training points, their k-NN predictions
        # may or may not change when the test point is added (only if the
        # test point enters their k-neighbourhood). For the test point,
        # its k-NN prediction is always computed from training labels only
        # (since we exclude self).
        #
        # Strategy: compute all training nonconformity scores under the
        # augmented distance matrix, then find the threshold.

        k = min(self.k, n_aug - 1)
        agg_func = np.mean if self.aggregation == "mean" else np.median

        # For each training point i (0..n-1), compute its leave-one-out
        # k-NN prediction in the augmented set (which includes the test point)
        D_work = D_aug.copy()
        np.fill_diagonal(D_work, np.inf)

        # k-NN predictions for all n+1 points (leave-one-out)
        knn_idx = np.argpartition(D_work, k - 1, axis=1)[:, :k]

        # For training points (0..n-1): their predictions don't depend on
        # the test label since we only use labels of their k neighbours.
        # BUT if the test point (index n) is among a training point's k-NN,
        # its label y_test would appear in the average. We handle this by
        # noting that y_test is what we're searching over.
        #
        # However, for the standard conformal regressor approach, we compute
        # the nonconformity scores as |y_i - hat_y_i| where hat_y_i is
        # computed with ALL labels known (including the hypothesized test label).
        # This makes the prediction interval depend on the test label y in a
        # complex way when the test point enters a training point's k-neighbourhood.
        #
        # Simplification (valid and common): Use the "deleted" approach where
        # for each training point i, the k-NN prediction uses the bag without i.
        # The test point's label only affects training points that have it as
        # a k-NN. For the test point itself, its prediction only uses training
        # labels.
        #
        # Even simpler (and what most practical k-NN conformal regressors do):
        # Compute training residuals on the ORIGINAL training set (without the
        # test point), then compare the test point's residual against them.
        # This is valid because the nonconformity scores are exchangeable.

        # Training residuals (on original training set, leave-one-out)
        y_hat_train = self._knn_predictions(self.D, self.y)
        alpha_train = np.abs(self.y - y_hat_train)

        # Test point k-NN prediction (using training data only)
        # Its k nearest neighbours among training points:
        k_test = min(self.k, n)
        test_knn_idx = np.argpartition(d, k_test - 1)[:k_test]
        y_hat_test = agg_func(self.y[test_knn_idx])

        # The test nonconformity score is |y - y_hat_test| and we need
        # p(y) = |{i: alpha_i >= |y - y_hat_test|}| / (n+1) > epsilon
        # The prediction interval is: y_hat_test +/- alpha_(ceil((1-eps)(n+1)))

        # Sort training residuals
        alpha_sorted = np.sort(alpha_train)

        # Build the interval
        if hasattr(epsilon, '__iter__'):
            predictions = {}
            for eps in epsilon:
                lo, up = self._compute_interval(alpha_sorted, y_hat_test, eps, n_aug, bounds)
                predictions[eps] = self._construct_Gamma(lo, up, eps)
            result = MultiLevelPredictionInterval(predictions)
        else:
            lo, up = self._compute_interval(alpha_sorted, y_hat_test, epsilon, n_aug, bounds)
            result = self._construct_Gamma(lo, up, epsilon)

        if return_update:
            return result, {"D": D_aug}
        return result

    @staticmethod
    def _compute_interval(alpha_sorted, y_hat, epsilon, n, bounds):
        """
        Compute prediction interval from sorted training residuals.

        The prediction set is {y : |y - y_hat| <= alpha_(j)} where
        j = ceil((1 - epsilon) * n) - 1 (0-indexed into sorted alphas).
        For two-sided: split epsilon equally.
        """
        if bounds == "both":
            # Two-sided: the interval is symmetric around y_hat
            # We need rank such that p-value > epsilon
            # threshold index: we want the (ceil((1-eps)*n) - 1)-th sorted alpha
            # (0-indexed), but there are only n-1 training alphas.
            idx = int(np.ceil((1 - epsilon) * n)) - 1
            if idx >= len(alpha_sorted):
                # All training residuals are included
                threshold = alpha_sorted[-1] if len(alpha_sorted) > 0 else np.inf
            elif idx < 0:
                return -np.inf, np.inf
            else:
                threshold = alpha_sorted[idx]
            return y_hat - threshold, y_hat + threshold
        elif bounds == "lower":
            idx = int(np.ceil((1 - epsilon) * n)) - 1
            if idx >= len(alpha_sorted):
                threshold = alpha_sorted[-1] if len(alpha_sorted) > 0 else np.inf
            elif idx < 0:
                return -np.inf, np.inf
            else:
                threshold = alpha_sorted[idx]
            return y_hat - threshold, np.inf
        elif bounds == "upper":
            idx = int(np.ceil((1 - epsilon) * n)) - 1
            if idx >= len(alpha_sorted):
                threshold = alpha_sorted[-1] if len(alpha_sorted) > 0 else np.inf
            elif idx < 0:
                return -np.inf, np.inf
            else:
                threshold = alpha_sorted[idx]
            return -np.inf, y_hat + threshold
        else:
            raise ValueError(f"bounds must be 'both', 'lower', or 'upper', got '{bounds}'")

    def compute_p_value(self, x: NDArray[np.floating[Any]], y: float, tau: float | None = None, smoothed: bool = True) -> float:
        """Compute conformal p-value for a test pair (x, y).

        Parameters
        ----------
        x : array-like, shape (d,)
            Test object.
        y : float
            Hypothesized label.
        tau : float, optional
            Randomization value in [0, 1] for smoothed p-value.
            If None and smoothed=True, drawn uniformly at random.
        smoothed : bool
            Whether to compute the smoothed p-value (default True).

        Returns
        -------
        float
            The conformal p-value.
        """
        x = np.asarray(x).ravel()
        n = self._safe_size_check(self.X)
        if n == 0:
            return 1.0

        if smoothed and tau is None:
            tau = self.rnd_gen.uniform(0, 1)

        # Training residuals (leave-one-out k-NN)
        y_hat_train = self._knn_predictions(self.D, self.y)
        alpha_train = np.abs(self.y - y_hat_train)

        # Test point k-NN prediction
        d = self.distance_func(self.X, x.reshape(1, -1)).ravel()
        k_test = min(self.k, n)
        agg_func = np.mean if self.aggregation == "mean" else np.median
        test_knn_idx = np.argpartition(d, k_test - 1)[:k_test]
        y_hat_test = agg_func(self.y[test_knn_idx])

        # Test nonconformity score
        alpha_test = np.abs(y - y_hat_test)

        # Compute p-value: fraction of training alphas >= test alpha
        Alpha = np.append(alpha_train, alpha_test)
        if smoothed:
            return self._compute_p_value(Alpha, tau, c_type="nonconformity")
        else:
            return self._compute_p_value(Alpha, c_type="nonconformity")

__init__(k=1, distance='euclidean', distance_func=None, aggregation='mean', verbose=0, rnd_state=None, epsilon=default_epsilon)

Parameters:

Name	Type	Description	Default
k	int	Number of nearest neighbours.	1
distance	str	Distance metric (passed to scipy.spatial.distance).	'euclidean'
distance_func	callable	Custom distance function. If provided, distance is ignored. Signature: distance_func(X, y=None) where X is (n, d) and y is (m, d) or None. Returns (n, n) if y is None, else (n, m).	None
aggregation	('mean', 'median')	How to aggregate k nearest neighbour labels.	'mean'
epsilon	float	Default significance level.	default_epsilon

Source code in src/online_cp/regressors.py

def __init__(
    self,
    k=1,
    distance="euclidean",
    distance_func=None,
    aggregation="mean",
    verbose=0,
    rnd_state=None,
    epsilon=default_epsilon,
):
    """
    Parameters
    ----------
    k : int
        Number of nearest neighbours.
    distance : str
        Distance metric (passed to scipy.spatial.distance).
    distance_func : callable, optional
        Custom distance function. If provided, `distance` is ignored.
        Signature: distance_func(X, y=None) where X is (n, d) and y is
        (m, d) or None. Returns (n, n) if y is None, else (n, m).
    aggregation : {'mean', 'median'}
        How to aggregate k nearest neighbour labels.
    epsilon : float
        Default significance level.
    """
    super().__init__(epsilon=epsilon)
    self.k = k
    if aggregation not in ("mean", "median"):
        raise ValueError(f"aggregation must be 'mean' or 'median', got '{aggregation}'")
    self.aggregation = aggregation
    self.distance = distance
    if distance_func is None:
        self.distance_func = self._standard_distance_func
    else:
        self.distance_func = distance_func
        self.distance = "custom"
    self.X = None
    self.y = None
    self.D = None
    self.verbose = verbose
    self.rnd_gen = np.random.default_rng(rnd_state)

learn_initial_training_set(X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None

Batch-initialize with training data.

Parameters:

Name	Type	Description	Default
X	(ndarray, shape(n, d))	Training objects.	required
y	(ndarray, shape(n))	Training labels.	required

Source code in src/online_cp/regressors.py

def learn_initial_training_set(self, X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None:
    """Batch-initialize with training data.

    Parameters
    ----------
    X : ndarray, shape (n, d)
        Training objects.
    y : ndarray, shape (n,)
        Training labels.
    """
    self.X = np.atleast_2d(X)
    self.y = np.asarray(y, dtype=float)
    self.D = self.distance_func(self.X)

learn_one(x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None

Incrementally add one observation.

Parameters:

Name	Type	Description	Default
x	(array - like, shape(d))	New object.	required
y	float	New label.	required
precomputed	dict	If provided, should contain 'D' (the already-updated distance matrix).	None

Source code in src/online_cp/regressors.py

def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
    """Incrementally add one observation.

    Parameters
    ----------
    x : array-like, shape (d,)
        New object.
    y : float
        New label.
    precomputed : dict, optional
        If provided, should contain 'D' (the already-updated distance matrix).
    """
    x = np.asarray(x).ravel()

    if self.X is None:
        self.X = x.reshape(1, -1)
        self.y = np.array([y], dtype=float)
        self.D = self.distance_func(self.X)
    else:
        if precomputed is not None and "D" in precomputed:
            self.D = precomputed["D"]
        else:
            d = self.distance_func(self.X, x.reshape(1, -1)).ravel()
            self.D = self._update_distance_matrix(self.D, d)
        self.X = np.vstack([self.X, x.reshape(1, -1)])
        self.y = np.append(self.y, y)

predict(x: NDArray[np.floating[Any]], epsilon: float | NDArray[np.floating[Any]] | None = None, bounds: str = 'both', return_update: bool = False) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]

Predict a conformal prediction interval for test object x.

Parameters:

Name	Type	Description	Default
x	(array - like, shape(d))	Test object.	required
epsilon	float or array - like	Significance level(s). If None, uses self.epsilon.	None
bounds	('both', 'lower', 'upper')	Which bounds to compute.	'both'
return_update	bool	If True, also return precomputed dict for subsequent learn_one.	False

Returns:

Name	Type	Description
result	ConformalPredictionInterval or MultiLevelPredictionInterval
precomputed	(dict, optional)	Returned if return_update=True. Contains 'D'.

Source code in src/online_cp/regressors.py

def predict(
    self,
    x: NDArray[np.floating[Any]],
    epsilon: float | NDArray[np.floating[Any]] | None = None,
    bounds: str = "both",
    return_update: bool = False,
) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
    """Predict a conformal prediction interval for test object x.

    Parameters
    ----------
    x : array-like, shape (d,)
        Test object.
    epsilon : float or array-like, optional
        Significance level(s). If None, uses self.epsilon.
    bounds : {'both', 'lower', 'upper'}
        Which bounds to compute.
    return_update : bool
        If True, also return precomputed dict for subsequent learn_one.

    Returns
    -------
    result : ConformalPredictionInterval or MultiLevelPredictionInterval
    precomputed : dict, optional
        Returned if return_update=True. Contains 'D'.
    """
    x = np.asarray(x).ravel()

    if epsilon is None:
        epsilon = self.epsilon

    n = self._safe_size_check(self.X)

    if n == 0:
        # No training data
        if hasattr(epsilon, '__iter__'):
            predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
            result = MultiLevelPredictionInterval(predictions)
        else:
            result = self._construct_Gamma(-np.inf, np.inf, epsilon)
        if return_update:
            return result, {"D": None}
        return result

    # Augment distance matrix with test point
    d = self.distance_func(self.X, x.reshape(1, -1)).ravel()
    D_aug = self._update_distance_matrix(self.D, d)
    n_aug = n + 1  # augmented size

    # Check significance level is feasible
    eps_check = max(epsilon) if hasattr(epsilon, '__iter__') else epsilon
    min_needed = 2 if bounds == "both" else 1
    if not (eps_check >= min_needed / n_aug):
        if hasattr(epsilon, '__iter__'):
            predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
            result = MultiLevelPredictionInterval(predictions)
        else:
            result = self._construct_Gamma(-np.inf, np.inf, epsilon)
        if return_update:
            return result, {"D": D_aug}
        return result

    # Compute leave-one-out k-NN predictions for all training points
    # in the augmented set. The test point's label is unknown, so we
    # need to determine the interval of y values that would be included.
    #
    # Key insight: For the n training points, their k-NN predictions
    # may or may not change when the test point is added (only if the
    # test point enters their k-neighbourhood). For the test point,
    # its k-NN prediction is always computed from training labels only
    # (since we exclude self).
    #
    # Strategy: compute all training nonconformity scores under the
    # augmented distance matrix, then find the threshold.

    k = min(self.k, n_aug - 1)
    agg_func = np.mean if self.aggregation == "mean" else np.median

    # For each training point i (0..n-1), compute its leave-one-out
    # k-NN prediction in the augmented set (which includes the test point)
    D_work = D_aug.copy()
    np.fill_diagonal(D_work, np.inf)

    # k-NN predictions for all n+1 points (leave-one-out)
    knn_idx = np.argpartition(D_work, k - 1, axis=1)[:, :k]

    # For training points (0..n-1): their predictions don't depend on
    # the test label since we only use labels of their k neighbours.
    # BUT if the test point (index n) is among a training point's k-NN,
    # its label y_test would appear in the average. We handle this by
    # noting that y_test is what we're searching over.
    #
    # However, for the standard conformal regressor approach, we compute
    # the nonconformity scores as |y_i - hat_y_i| where hat_y_i is
    # computed with ALL labels known (including the hypothesized test label).
    # This makes the prediction interval depend on the test label y in a
    # complex way when the test point enters a training point's k-neighbourhood.
    #
    # Simplification (valid and common): Use the "deleted" approach where
    # for each training point i, the k-NN prediction uses the bag without i.
    # The test point's label only affects training points that have it as
    # a k-NN. For the test point itself, its prediction only uses training
    # labels.
    #
    # Even simpler (and what most practical k-NN conformal regressors do):
    # Compute training residuals on the ORIGINAL training set (without the
    # test point), then compare the test point's residual against them.
    # This is valid because the nonconformity scores are exchangeable.

    # Training residuals (on original training set, leave-one-out)
    y_hat_train = self._knn_predictions(self.D, self.y)
    alpha_train = np.abs(self.y - y_hat_train)

    # Test point k-NN prediction (using training data only)
    # Its k nearest neighbours among training points:
    k_test = min(self.k, n)
    test_knn_idx = np.argpartition(d, k_test - 1)[:k_test]
    y_hat_test = agg_func(self.y[test_knn_idx])

    # The test nonconformity score is |y - y_hat_test| and we need
    # p(y) = |{i: alpha_i >= |y - y_hat_test|}| / (n+1) > epsilon
    # The prediction interval is: y_hat_test +/- alpha_(ceil((1-eps)(n+1)))

    # Sort training residuals
    alpha_sorted = np.sort(alpha_train)

    # Build the interval
    if hasattr(epsilon, '__iter__'):
        predictions = {}
        for eps in epsilon:
            lo, up = self._compute_interval(alpha_sorted, y_hat_test, eps, n_aug, bounds)
            predictions[eps] = self._construct_Gamma(lo, up, eps)
        result = MultiLevelPredictionInterval(predictions)
    else:
        lo, up = self._compute_interval(alpha_sorted, y_hat_test, epsilon, n_aug, bounds)
        result = self._construct_Gamma(lo, up, epsilon)

    if return_update:
        return result, {"D": D_aug}
    return result

compute_p_value(x: NDArray[np.floating[Any]], y: float, tau: float | None = None, smoothed: bool = True) -> float

Compute conformal p-value for a test pair (x, y).

Parameters:

Name	Type	Description	Default
x	(array - like, shape(d))	Test object.	required
y	float	Hypothesized label.	required
tau	float	Randomization value in [0, 1] for smoothed p-value. If None and smoothed=True, drawn uniformly at random.	None
smoothed	bool	Whether to compute the smoothed p-value (default True).	True

Returns:

Type	Description
float	The conformal p-value.

Source code in src/online_cp/regressors.py

def compute_p_value(self, x: NDArray[np.floating[Any]], y: float, tau: float | None = None, smoothed: bool = True) -> float:
    """Compute conformal p-value for a test pair (x, y).

    Parameters
    ----------
    x : array-like, shape (d,)
        Test object.
    y : float
        Hypothesized label.
    tau : float, optional
        Randomization value in [0, 1] for smoothed p-value.
        If None and smoothed=True, drawn uniformly at random.
    smoothed : bool
        Whether to compute the smoothed p-value (default True).

    Returns
    -------
    float
        The conformal p-value.
    """
    x = np.asarray(x).ravel()
    n = self._safe_size_check(self.X)
    if n == 0:
        return 1.0

    if smoothed and tau is None:
        tau = self.rnd_gen.uniform(0, 1)

    # Training residuals (leave-one-out k-NN)
    y_hat_train = self._knn_predictions(self.D, self.y)
    alpha_train = np.abs(self.y - y_hat_train)

    # Test point k-NN prediction
    d = self.distance_func(self.X, x.reshape(1, -1)).ravel()
    k_test = min(self.k, n)
    agg_func = np.mean if self.aggregation == "mean" else np.median
    test_knn_idx = np.argpartition(d, k_test - 1)[:k_test]
    y_hat_test = agg_func(self.y[test_knn_idx])

    # Test nonconformity score
    alpha_test = np.abs(y - y_hat_test)

    # Compute p-value: fraction of training alphas >= test alpha
    Alpha = np.append(alpha_train, alpha_test)
    if smoothed:
        return self._compute_p_value(Alpha, tau, c_type="nonconformity")
    else:
        return self._compute_p_value(Alpha, c_type="nonconformity")

online_cp.regressors.KernelConformalRidgeRegressor

Bases: ConformalRegressor

Source code in src/online_cp/regressors.py

class KernelConformalRidgeRegressor(ConformalRegressor):
    # TODO Add doctests to methods where applicable

    def __init__(self, kernel: Any, a: float = 0, warnings: bool = True, verbose: int = 0, rnd_state: int | None = None, epsilon: float = default_epsilon) -> None:
        """
        KernelConformalRidgeRegressor requires a kernel. Some common kernels are found in kernels.py, but it is
        also compatible with (most) kernels from e.g. scikit-learn.
        Custom kernels can also be passed as callable functions.
        """
        super().__init__(epsilon=epsilon)

        self.a = a
        self.X = None
        self.y = None
        self.p = None
        self.Id = None
        self.K = None
        self.Kinv = None

        self.kernel = kernel

        # Should we raise warnings
        self.warnings = warnings

        self.verbose = verbose

        self.rnd_gen = np.random.default_rng(rnd_state)

    def learn_initial_training_set(self, X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None:
        self.X = X
        self.y = y
        Id = np.identity(self.X.shape[0])

        self.K = self.kernel(self.X)
        self.Kinv = np.linalg.inv(self.K + self.a * Id)

    @staticmethod
    def _update_Kinv(Kinv, k, kappa):
        # print(f'K: {K}')
        # print(f'k: {k}')
        # print(f'kappa: {kappa}')
        d = 1 / (kappa - k.T @ Kinv @ k)
        return np.block([[Kinv + d * Kinv @ k @ k.T @ Kinv, -d * Kinv @ k], [-d * k.T @ Kinv, d]])

    @staticmethod
    def _update_K(K, k, kappa):
        # print(f'K: {K}')
        # print(f'k: {k}')
        # print(f'kappa: {kappa}')
        return np.block([[K, k], [k.T, kappa]])

    def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
        """
        Learn a single example
        """
        # Learn label y
        if self.y is None:
            self.y = np.array([y])
        else:
            self.y = np.append(self.y, y)

        if precomputed is not None:
            X = precomputed["X"]
            K = precomputed["K"]
            Kinv = precomputed["Kinv"]

            if X is not None:
                self.X = X

            if K is not None and Kinv is not None:
                self.K = K
                self.Kinv = Kinv

            else:
                Id = np.identity(self.X.shape[0])
                self.K = self.kernel(self.X)
                self.Kinv = np.linalg.inv(self.K + self.a * Id)

        else:
            # Learn object x
            if self.X is None:
                self.X = x.reshape(1, -1)
                Id = np.identity(self.X.shape[0])
                self.K = self.kernel(self.X)
                self.Kinv = np.linalg.inv(self.K + self.a * Id)
            elif self.X.shape[0] == 1:
                self.X = np.append(self.X, x.reshape(1, -1), axis=0)
                Id = np.identity(self.X.shape[0])
                self.K = self.kernel(self.X)
                self.Kinv = np.linalg.inv(self.K + self.a * Id)
            else:
                k = self.kernel(self.X, x).reshape(-1, 1)
                kappa = self.kernel(x, x)
                self.K = self._update_K(self.K, k, kappa)
                self.Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)
                self.X = np.append(self.X, x.reshape(1, -1), axis=0)

    @staticmethod
    def compute_A_and_B(X, K, Kinv, y):
        # print(f'X: {X}')
        # print(f'K: {K}')
        # print(f'Kinv: {Kinv}')
        # print(f'y: {y}')
        n = X.shape[0]
        H = Kinv @ K
        C = np.identity(n) - H
        A = C @ np.append(y, 0)  # Elements of this vector are denoted ai
        B = C @ np.append(np.zeros((n - 1,)), 1)  # Elements of this vector are denoted bi
        # Nonconformity scores are A + yB = y - yhat
        return A, B

    def predict(
        self,
        x: NDArray[np.floating[Any]],
        epsilon: float | NDArray[np.floating[Any]] | None = None,
        bounds: str = "both",
        return_update: bool = False,
        debug_time: bool = False,
    ) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
        """
        This function makes a prediction.

        If you start with no training,
        you get a null prediciton between
        -infinity and +infinity.

        TODO Add possibility to learn object to save time

        >>> cp = ConformalRidgeRegressor()
        >>> cp.predict(np.array([0.506, 0.22, -0.45]), bounds="both")
        (-inf, inf)
        """

        def build_precomputed(X, K, Kinv, A, B):
            computed = {
                "X": X,  # The updated matrix of objects
                "K": K,  # The updated kernel matrix
                "Kinv": Kinv,
                "A": A,
                "B": B,
            }
            return computed

        if epsilon is None:
            epsilon = self.epsilon

        if self.X is not None:
            tic = time.time()

            # Temporarily update kernel matrix
            k = self.kernel(self.X, x).reshape(-1, 1)
            kappa = self.kernel(x, x)
            K = self._update_K(self.K, k, kappa)
            Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)

            toc_update_kernel = time.time() - tic

            tic = time.time()
            # Add row to X matrix
            X = np.append(self.X, x.reshape(1, -1), axis=0)
            toc_add_row = time.time() - tic
            n = X.shape[0]

            # Check that the significance level is not too small. If it is, return infinite prediction interval
            eps_check = min(epsilon) if hasattr(epsilon, '__iter__') else epsilon
            if bounds == "both":
                if not (eps_check >= 2 / n):
                    if self.warnings:
                        warnings.warn(
                            f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(2 / eps_check))} examples. Increase or add more examples",
                            stacklevel=2,
                        )
                    if hasattr(epsilon, '__iter__'):
                        predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                        result = MultiLevelPredictionInterval(predictions)
                    else:
                        result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                    if return_update:
                        return result, build_precomputed(
                            X, K, Kinv, None, None
                        )
                    else:
                        return result
            else:
                if not (eps_check >= 1 / n):
                    if self.warnings:
                        warnings.warn(
                            f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(1 / eps_check))} examples. Increase or add more examples",
                            stacklevel=2,
                        )
                    if hasattr(epsilon, '__iter__'):
                        predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                        result = MultiLevelPredictionInterval(predictions)
                    else:
                        result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                    if return_update:
                        return result, build_precomputed(
                            X, K, Kinv, None, None
                        )
                    else:
                        return result

            tic = time.time()
            A, B = self.compute_A_and_B(X, K, Kinv, self.y)
            toc_nc = time.time() - tic

            tic = time.time()
            l_dic, u_dic = self._vectorised_l_and_u(A, B)
            toc_dics = time.time() - tic

            if bounds == "both":
                if hasattr(epsilon, '__iter__'):
                    predictions = {}
                    for eps in epsilon:
                        lo = self._get_lower(l_dic=l_dic, epsilon=eps / 2, n=n)
                        up = self._get_upper(u_dic=u_dic, epsilon=eps / 2, n=n)
                        predictions[eps] = self._construct_Gamma(lo, up, eps)
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    lower = self._get_lower(l_dic=l_dic, epsilon=epsilon / 2, n=n)
                    upper = self._get_upper(u_dic=u_dic, epsilon=epsilon / 2, n=n)
                    result = self._construct_Gamma(lower, upper, epsilon)
            elif bounds == "lower":
                if hasattr(epsilon, '__iter__'):
                    predictions = {}
                    for eps in epsilon:
                        lo = self._get_lower(l_dic=l_dic, epsilon=eps, n=n)
                        predictions[eps] = self._construct_Gamma(lo, np.inf, eps)
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    lower = self._get_lower(l_dic=l_dic, epsilon=epsilon, n=n)
                    result = self._construct_Gamma(lower, np.inf, epsilon)
            elif bounds == "upper":
                if hasattr(epsilon, '__iter__'):
                    predictions = {}
                    for eps in epsilon:
                        up = self._get_upper(u_dic=u_dic, epsilon=eps, n=n)
                        predictions[eps] = self._construct_Gamma(-np.inf, up, eps)
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    upper = self._get_upper(u_dic=u_dic, epsilon=epsilon, n=n)
                    result = self._construct_Gamma(-np.inf, upper, epsilon)
            else:
                raise Exception

            if debug_time:
                print(f"Add row: {toc_add_row}")
                print(f"Update kernel: {toc_update_kernel}")
                print(f"NC scores: {toc_nc}")
                print(f"l and u: {toc_dics}")
                print()
        else:
            # With just one object, and no label, we cannot predict any meaningful interval
            X = x.reshape(1, -1)
            K = None
            Kinv = None
            A = None
            B = None

            if hasattr(epsilon, '__iter__'):
                predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                result = MultiLevelPredictionInterval(predictions)
            else:
                result = self._construct_Gamma(-np.inf, np.inf, epsilon)

        if return_update:
            return result, build_precomputed(X, K, Kinv, A, B)
        else:
            return result

    def compute_p_value(self, x, y, bounds="both", precomputed=None, tau=None, smoothed=True):
        """
        Computes the smoothed p-value of the example (x, y).
        """
        if tau is None and smoothed:
            tau = self.rnd_gen.uniform(0, 1)
        if precomputed is not None:
            assert np.allclose(x, precomputed["X"][-1])
            A = precomputed["A"]
            B = precomputed["B"]

        else:
            if self.Kinv is not None:
                k = self.kernel(self.X, x).reshape(-1, 1)
                kappa = self.kernel(x, x)
                K = self._update_K(self.K, k, kappa)
                Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)
                X = np.append(self.X, x.reshape(1, -1), axis=0)
                A, B = self.compute_A_and_B(X, K, Kinv, self.y)
            else:
                A, B = None, None

        if A is not None and B is not None:
            if bounds == "both":
                E = A + y * B
                Alpha = np.zeros_like(A)
                for i, e in enumerate(E):
                    alpha = min((E >= e).sum(), (E <= e).sum())
                    Alpha[i] = alpha
                c_type = "conformity"
            elif bounds == "lower":
                Alpha = -(A + y * B)
                c_type = "nonconformity"
            elif bounds == "upper":
                Alpha = A + y * B
                c_type = "nonconformity"
            else:
                raise Exception('bounds must be one of "both", "lower", "upper"')

            if smoothed:
                p = self._compute_p_value(Alpha, tau, c_type=c_type)
            else:
                p = self._compute_p_value(Alpha, c_type=c_type)
        else:
            if smoothed:
                p = tau
            else:
                p = 1

        return p

    def compute_smoothed_p_value(self, x, y, precomputed=None):
        """
        Computes the smoothed p-value of the example (x, y).
        Smoothed p-values can be used to test the exchangeability assumption.
        """

        # Inner method to compute the p-value from NC scores
        def calc_p(A, B, y):
            # Nonconformity scores are A + yB = y - yhat
            Alpha = A + y * B
            alpha_y = Alpha[-1]
            gt = np.where(Alpha > alpha_y)[0].shape[0]
            eq = np.where(Alpha == alpha_y)[0].shape[0]
            tau = self.rnd_gen.uniform(0, 1)
            p_y = (gt + tau * eq) / Alpha.shape[0]
            return p_y

        if precomputed is not None:
            A = precomputed["A"]
            B = precomputed["B"]
            X = precomputed["X"]
            K = precomputed["K"]
            Kinv = precomputed["Kinv"]

            if A is not None and B is not None:
                p_y = calc_p(A, B, y)
            else:
                if Kinv is not None and X is not None and K is not None:
                    A, B = self.compute_A_and_B(X, K, Kinv, self.y)
                    p_y = calc_p(A, B, y)

                else:
                    p_y = self.rnd_gen.uniform(0, 1)

        else:
            if self.Kinv is not None:
                k = self.kernel(self.X, x).reshape(-1, 1)
                kappa = self.kernel(x, x)
                K = self._update_K(self.K, k, kappa)
                Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)
                X = np.append(self.X, x.reshape(1, -1), axis=0)

                A, B = self.compute_A_and_B(X, K, Kinv, self.y)
                p_y = calc_p(A, B, y)

            else:
                p_y = self.rnd_gen.uniform(0, 1)
        return p_y

__init__(kernel: Any, a: float = 0, warnings: bool = True, verbose: int = 0, rnd_state: int | None = None, epsilon: float = default_epsilon) -> None

KernelConformalRidgeRegressor requires a kernel. Some common kernels are found in kernels.py, but it is also compatible with (most) kernels from e.g. scikit-learn. Custom kernels can also be passed as callable functions.

Source code in src/online_cp/regressors.py

def __init__(self, kernel: Any, a: float = 0, warnings: bool = True, verbose: int = 0, rnd_state: int | None = None, epsilon: float = default_epsilon) -> None:
    """
    KernelConformalRidgeRegressor requires a kernel. Some common kernels are found in kernels.py, but it is
    also compatible with (most) kernels from e.g. scikit-learn.
    Custom kernels can also be passed as callable functions.
    """
    super().__init__(epsilon=epsilon)

    self.a = a
    self.X = None
    self.y = None
    self.p = None
    self.Id = None
    self.K = None
    self.Kinv = None

    self.kernel = kernel

    # Should we raise warnings
    self.warnings = warnings

    self.verbose = verbose

    self.rnd_gen = np.random.default_rng(rnd_state)

learn_one(x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None

Learn a single example

Source code in src/online_cp/regressors.py

def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
    """
    Learn a single example
    """
    # Learn label y
    if self.y is None:
        self.y = np.array([y])
    else:
        self.y = np.append(self.y, y)

    if precomputed is not None:
        X = precomputed["X"]
        K = precomputed["K"]
        Kinv = precomputed["Kinv"]

        if X is not None:
            self.X = X

        if K is not None and Kinv is not None:
            self.K = K
            self.Kinv = Kinv

        else:
            Id = np.identity(self.X.shape[0])
            self.K = self.kernel(self.X)
            self.Kinv = np.linalg.inv(self.K + self.a * Id)

    else:
        # Learn object x
        if self.X is None:
            self.X = x.reshape(1, -1)
            Id = np.identity(self.X.shape[0])
            self.K = self.kernel(self.X)
            self.Kinv = np.linalg.inv(self.K + self.a * Id)
        elif self.X.shape[0] == 1:
            self.X = np.append(self.X, x.reshape(1, -1), axis=0)
            Id = np.identity(self.X.shape[0])
            self.K = self.kernel(self.X)
            self.Kinv = np.linalg.inv(self.K + self.a * Id)
        else:
            k = self.kernel(self.X, x).reshape(-1, 1)
            kappa = self.kernel(x, x)
            self.K = self._update_K(self.K, k, kappa)
            self.Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)
            self.X = np.append(self.X, x.reshape(1, -1), axis=0)

predict(x: NDArray[np.floating[Any]], epsilon: float | NDArray[np.floating[Any]] | None = None, bounds: str = 'both', return_update: bool = False, debug_time: bool = False) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]

This function makes a prediction.

If you start with no training, you get a null prediciton between -infinity and +infinity.

TODO Add possibility to learn object to save time

cp = ConformalRidgeRegressor() cp.predict(np.array([0.506, 0.22, -0.45]), bounds="both") (-inf, inf)

Source code in src/online_cp/regressors.py

def predict(
    self,
    x: NDArray[np.floating[Any]],
    epsilon: float | NDArray[np.floating[Any]] | None = None,
    bounds: str = "both",
    return_update: bool = False,
    debug_time: bool = False,
) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
    """
    This function makes a prediction.

    If you start with no training,
    you get a null prediciton between
    -infinity and +infinity.

    TODO Add possibility to learn object to save time

    >>> cp = ConformalRidgeRegressor()
    >>> cp.predict(np.array([0.506, 0.22, -0.45]), bounds="both")
    (-inf, inf)
    """

    def build_precomputed(X, K, Kinv, A, B):
        computed = {
            "X": X,  # The updated matrix of objects
            "K": K,  # The updated kernel matrix
            "Kinv": Kinv,
            "A": A,
            "B": B,
        }
        return computed

    if epsilon is None:
        epsilon = self.epsilon

    if self.X is not None:
        tic = time.time()

        # Temporarily update kernel matrix
        k = self.kernel(self.X, x).reshape(-1, 1)
        kappa = self.kernel(x, x)
        K = self._update_K(self.K, k, kappa)
        Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)

        toc_update_kernel = time.time() - tic

        tic = time.time()
        # Add row to X matrix
        X = np.append(self.X, x.reshape(1, -1), axis=0)
        toc_add_row = time.time() - tic
        n = X.shape[0]

        # Check that the significance level is not too small. If it is, return infinite prediction interval
        eps_check = min(epsilon) if hasattr(epsilon, '__iter__') else epsilon
        if bounds == "both":
            if not (eps_check >= 2 / n):
                if self.warnings:
                    warnings.warn(
                        f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(2 / eps_check))} examples. Increase or add more examples",
                        stacklevel=2,
                    )
                if hasattr(epsilon, '__iter__'):
                    predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                if return_update:
                    return result, build_precomputed(
                        X, K, Kinv, None, None
                    )
                else:
                    return result
        else:
            if not (eps_check >= 1 / n):
                if self.warnings:
                    warnings.warn(
                        f"Significance level epsilon is too small for training set. Need at least {int(np.ceil(1 / eps_check))} examples. Increase or add more examples",
                        stacklevel=2,
                    )
                if hasattr(epsilon, '__iter__'):
                    predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
                    result = MultiLevelPredictionInterval(predictions)
                else:
                    result = self._construct_Gamma(-np.inf, np.inf, epsilon)
                if return_update:
                    return result, build_precomputed(
                        X, K, Kinv, None, None
                    )
                else:
                    return result

        tic = time.time()
        A, B = self.compute_A_and_B(X, K, Kinv, self.y)
        toc_nc = time.time() - tic

        tic = time.time()
        l_dic, u_dic = self._vectorised_l_and_u(A, B)
        toc_dics = time.time() - tic

        if bounds == "both":
            if hasattr(epsilon, '__iter__'):
                predictions = {}
                for eps in epsilon:
                    lo = self._get_lower(l_dic=l_dic, epsilon=eps / 2, n=n)
                    up = self._get_upper(u_dic=u_dic, epsilon=eps / 2, n=n)
                    predictions[eps] = self._construct_Gamma(lo, up, eps)
                result = MultiLevelPredictionInterval(predictions)
            else:
                lower = self._get_lower(l_dic=l_dic, epsilon=epsilon / 2, n=n)
                upper = self._get_upper(u_dic=u_dic, epsilon=epsilon / 2, n=n)
                result = self._construct_Gamma(lower, upper, epsilon)
        elif bounds == "lower":
            if hasattr(epsilon, '__iter__'):
                predictions = {}
                for eps in epsilon:
                    lo = self._get_lower(l_dic=l_dic, epsilon=eps, n=n)
                    predictions[eps] = self._construct_Gamma(lo, np.inf, eps)
                result = MultiLevelPredictionInterval(predictions)
            else:
                lower = self._get_lower(l_dic=l_dic, epsilon=epsilon, n=n)
                result = self._construct_Gamma(lower, np.inf, epsilon)
        elif bounds == "upper":
            if hasattr(epsilon, '__iter__'):
                predictions = {}
                for eps in epsilon:
                    up = self._get_upper(u_dic=u_dic, epsilon=eps, n=n)
                    predictions[eps] = self._construct_Gamma(-np.inf, up, eps)
                result = MultiLevelPredictionInterval(predictions)
            else:
                upper = self._get_upper(u_dic=u_dic, epsilon=epsilon, n=n)
                result = self._construct_Gamma(-np.inf, upper, epsilon)
        else:
            raise Exception

        if debug_time:
            print(f"Add row: {toc_add_row}")
            print(f"Update kernel: {toc_update_kernel}")
            print(f"NC scores: {toc_nc}")
            print(f"l and u: {toc_dics}")
            print()
    else:
        # With just one object, and no label, we cannot predict any meaningful interval
        X = x.reshape(1, -1)
        K = None
        Kinv = None
        A = None
        B = None

        if hasattr(epsilon, '__iter__'):
            predictions = {eps: self._construct_Gamma(-np.inf, np.inf, eps) for eps in epsilon}
            result = MultiLevelPredictionInterval(predictions)
        else:
            result = self._construct_Gamma(-np.inf, np.inf, epsilon)

    if return_update:
        return result, build_precomputed(X, K, Kinv, A, B)
    else:
        return result

compute_p_value(x, y, bounds='both', precomputed=None, tau=None, smoothed=True)

Computes the smoothed p-value of the example (x, y).

Source code in src/online_cp/regressors.py

def compute_p_value(self, x, y, bounds="both", precomputed=None, tau=None, smoothed=True):
    """
    Computes the smoothed p-value of the example (x, y).
    """
    if tau is None and smoothed:
        tau = self.rnd_gen.uniform(0, 1)
    if precomputed is not None:
        assert np.allclose(x, precomputed["X"][-1])
        A = precomputed["A"]
        B = precomputed["B"]

    else:
        if self.Kinv is not None:
            k = self.kernel(self.X, x).reshape(-1, 1)
            kappa = self.kernel(x, x)
            K = self._update_K(self.K, k, kappa)
            Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)
            X = np.append(self.X, x.reshape(1, -1), axis=0)
            A, B = self.compute_A_and_B(X, K, Kinv, self.y)
        else:
            A, B = None, None

    if A is not None and B is not None:
        if bounds == "both":
            E = A + y * B
            Alpha = np.zeros_like(A)
            for i, e in enumerate(E):
                alpha = min((E >= e).sum(), (E <= e).sum())
                Alpha[i] = alpha
            c_type = "conformity"
        elif bounds == "lower":
            Alpha = -(A + y * B)
            c_type = "nonconformity"
        elif bounds == "upper":
            Alpha = A + y * B
            c_type = "nonconformity"
        else:
            raise Exception('bounds must be one of "both", "lower", "upper"')

        if smoothed:
            p = self._compute_p_value(Alpha, tau, c_type=c_type)
        else:
            p = self._compute_p_value(Alpha, c_type=c_type)
    else:
        if smoothed:
            p = tau
        else:
            p = 1

    return p

compute_smoothed_p_value(x, y, precomputed=None)

Computes the smoothed p-value of the example (x, y). Smoothed p-values can be used to test the exchangeability assumption.

Source code in src/online_cp/regressors.py

def compute_smoothed_p_value(self, x, y, precomputed=None):
    """
    Computes the smoothed p-value of the example (x, y).
    Smoothed p-values can be used to test the exchangeability assumption.
    """

    # Inner method to compute the p-value from NC scores
    def calc_p(A, B, y):
        # Nonconformity scores are A + yB = y - yhat
        Alpha = A + y * B
        alpha_y = Alpha[-1]
        gt = np.where(Alpha > alpha_y)[0].shape[0]
        eq = np.where(Alpha == alpha_y)[0].shape[0]
        tau = self.rnd_gen.uniform(0, 1)
        p_y = (gt + tau * eq) / Alpha.shape[0]
        return p_y

    if precomputed is not None:
        A = precomputed["A"]
        B = precomputed["B"]
        X = precomputed["X"]
        K = precomputed["K"]
        Kinv = precomputed["Kinv"]

        if A is not None and B is not None:
            p_y = calc_p(A, B, y)
        else:
            if Kinv is not None and X is not None and K is not None:
                A, B = self.compute_A_and_B(X, K, Kinv, self.y)
                p_y = calc_p(A, B, y)

            else:
                p_y = self.rnd_gen.uniform(0, 1)

    else:
        if self.Kinv is not None:
            k = self.kernel(self.X, x).reshape(-1, 1)
            kappa = self.kernel(x, x)
            K = self._update_K(self.K, k, kappa)
            Kinv = self._update_Kinv(self.Kinv, k, kappa + self.a)
            X = np.append(self.X, x.reshape(1, -1), axis=0)

            A, B = self.compute_A_and_B(X, K, Kinv, self.y)
            p_y = calc_p(A, B, y)

        else:
            p_y = self.rnd_gen.uniform(0, 1)
    return p_y

online_cp.regressors.ConformalLassoRegressor

Bases: ConformalRegressor

Conformal prediction with Lasso/elastic net using the piecewise linear homotopy from Lei (2019). Computes exact conformal prediction sets without grid search by tracing how residuals evolve as the test label varies.

References

Lei, J. (2019). Fast exact conformalization of the Lasso using piecewise linear homotopy. Biometrika, 106(4), 751–767.

When rho=0 (default), this is pure Lasso. When rho>0, this solves the elastic net: min (1/2)||y - Xβ||² + lam||β||₁ + (rho/2)||β||₂²

Parameters:

Name	Type	Description	Default
lam	float	L1 regularization parameter (lambda). Must be non-negative.	1.0
rho	float	L2 regularization parameter (default 0.0). When rho=0, pure Lasso.	0.0
epsilon	float	Significance level for prediction sets (default 0.1).	default_epsilon
autotune	bool	If True, tune lambda via K-fold cross-validation in learn_initial_training_set.	False
n_folds	int	Number of CV folds for autotuning (default 5).	5
search_range_factor	float	Factor to extend the y search range beyond [y_min, y_max] (default 0.25).	0.25
max_homotopy_steps	int	Maximum number of homotopy breakpoints to trace per direction (default 1000).	1000
verbose	int	Verbosity level.	0
warnings	bool	Whether to emit warnings.	True
rnd_state	int or None	Random seed.	None

Examples:

>>> import numpy as np
>>> np.random.seed(42)
>>> N = 50
>>> X = np.random.normal(size=(N, 10))
>>> beta_true = np.array([3, 1.5, 0, 0, 2, 0, 0, 0, 0, 0])
>>> y = X @ beta_true + np.random.normal(scale=0.5, size=N)
>>> cp = ConformalLassoRegressor(lam=0.5)
>>> cp.learn_initial_training_set(X[:30], y[:30])
>>> interval = cp.predict(X[30], epsilon=0.1)
>>> y[30] in interval
True

Source code in src/online_cp/regressors.py

class ConformalLassoRegressor(ConformalRegressor):
    """
    Conformal prediction with Lasso/elastic net using the piecewise linear
    homotopy from Lei (2019). Computes exact conformal prediction sets
    without grid search by tracing how residuals evolve as the test label varies.

    References
    ----------
    Lei, J. (2019). Fast exact conformalization of the Lasso using piecewise
    linear homotopy. *Biometrika*, 106(4), 751–767.

    When rho=0 (default), this is pure Lasso. When rho>0, this solves the
    elastic net: min (1/2)||y - Xβ||² + lam||β||₁ + (rho/2)||β||₂²

    Parameters
    ----------
    lam : float
        L1 regularization parameter (lambda). Must be non-negative.
    rho : float
        L2 regularization parameter (default 0.0). When rho=0, pure Lasso.
    epsilon : float
        Significance level for prediction sets (default 0.1).
    autotune : bool
        If True, tune lambda via K-fold cross-validation in learn_initial_training_set.
    n_folds : int
        Number of CV folds for autotuning (default 5).
    search_range_factor : float
        Factor to extend the y search range beyond [y_min, y_max] (default 0.25).
    max_homotopy_steps : int
        Maximum number of homotopy breakpoints to trace per direction (default 1000).
    verbose : int
        Verbosity level.
    warnings : bool
        Whether to emit warnings.
    rnd_state : int or None
        Random seed.

    Examples
    --------
    >>> import numpy as np
    >>> np.random.seed(42)
    >>> N = 50
    >>> X = np.random.normal(size=(N, 10))
    >>> beta_true = np.array([3, 1.5, 0, 0, 2, 0, 0, 0, 0, 0])
    >>> y = X @ beta_true + np.random.normal(scale=0.5, size=N)
    >>> cp = ConformalLassoRegressor(lam=0.5)
    >>> cp.learn_initial_training_set(X[:30], y[:30])
    >>> interval = cp.predict(X[30], epsilon=0.1)
    >>> y[30] in interval
    True
    """

    def __init__(
        self,
        lam=1.0,
        rho=0.0,
        epsilon=default_epsilon,
        autotune=False,
        n_folds=5,
        search_range_factor=0.25,
        max_homotopy_steps=1000,
        verbose=0,
        warnings=True,
        rnd_state=None,
    ):
        super().__init__(epsilon=epsilon)
        self.lam = lam
        self.rho = rho
        self.autotune = autotune
        self.n_folds = n_folds
        self.search_range_factor = search_range_factor
        self.max_homotopy_steps = max_homotopy_steps
        self.verbose = verbose
        self.warnings = warnings
        self.rnd_gen = np.random.default_rng(rnd_state)

        self.X = None
        self.y = None
        self.beta = None  # Current Lasso solution
        self.Sigma = None  # X^T X / n (sample covariance)

    def learn_initial_training_set(self, X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None:
        """Fit initial Lasso on the training data."""
        self.X = X.copy()
        self.y = y.copy()
        n, p = X.shape
        self.Sigma = X.T @ X / n

        if self.autotune:
            self._tune_lambda()

        self.beta = _solve_lasso(self.X, self.y, self.lam, rho=self.rho)

    def _tune_lambda(self):
        """Tune lambda via K-fold cross-validation."""
        n = self.X.shape[0]
        indices = np.arange(n)
        self.rnd_gen.shuffle(indices)
        folds = np.array_split(indices, self.n_folds)

        # Lambda grid: geometric sequence
        lam_max = np.max(np.abs(self.X.T @ self.y)) / n
        lam_grid = np.geomspace(lam_max, lam_max * 1e-3, num=50)

        best_lam = self.lam
        best_mse = np.inf

        for lam in lam_grid:
            mse = 0.0
            for k in range(self.n_folds):
                val_idx = folds[k]
                train_idx = np.concatenate([folds[j] for j in range(self.n_folds) if j != k])
                X_tr, y_tr = self.X[train_idx], self.y[train_idx]
                X_val, y_val = self.X[val_idx], self.y[val_idx]
                beta_k = _solve_lasso(X_tr, y_tr, lam, rho=self.rho)
                mse += np.mean((y_val - X_val @ beta_k) ** 2)
            mse /= self.n_folds
            if mse < best_mse:
                best_mse = mse
                best_lam = lam

        self.lam = best_lam
        if self.verbose > 0:
            print(f"Tuned lambda: {self.lam:.6f} (CV MSE: {best_mse:.4f})")

    def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
        """
        Learn a new data point. Updates the training set and refits Lasso.

        If precomputed is provided (from predict with return_update=True),
        uses the cached Lasso solution to avoid refitting.
        """
        x = np.atleast_1d(x).ravel()

        # Update training set
        if self.X is None:
            self.X = x.reshape(1, -1)
            self.y = np.array([y])
        else:
            self.X = np.vstack([self.X, x.reshape(1, -1)])
            self.y = np.append(self.y, y)

        # Update sample covariance incrementally: Sigma_new = (n*Sigma_old + x x^T) / (n+1)
        n = self.X.shape[0]
        self.Sigma = ((n - 1) * self.Sigma + np.outer(x, x)) / n

        if precomputed is not None and precomputed.get("beta") is not None:
            self.beta = precomputed["beta"]
        else:
            # Refit from scratch (warm-started from current beta)
            self.beta = _solve_lasso(self.X, self.y, self.lam, rho=self.rho, warm_start=self.beta)

    def predict(
        self,
        x: NDArray[np.floating[Any]],
        epsilon: float | NDArray[np.floating[Any]] | None = None,
        return_update: bool = False,
    ) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
        """
        Compute the conformal prediction set at x using the homotopy algorithm.

        Returns a ConformalPredictionInterval (if the set is a single interval)
        or a list of (lower, upper) tuples otherwise.

        If epsilon is a list/array, returns a MultiLevelPredictionInterval.
        """
        if epsilon is None:
            epsilon = self.epsilon

        # Handle multi-level epsilon
        if hasattr(epsilon, '__iter__'):
            predictions = {}
            for eps in epsilon:
                result = self.predict(x, epsilon=eps, return_update=False)
                predictions[eps] = result
            result = MultiLevelPredictionInterval(predictions)
            if return_update:
                return result, {"beta": None}
            return result

        x = np.atleast_1d(x).ravel()
        n = self.X.shape[0]

        if n < 2:
            result = self._construct_Gamma(-np.inf, np.inf, epsilon)
            if return_update:
                return result, {"beta": None}
            return result

        # y_{n+1}(0) = x_{n+1}^T beta_hat — the "neutral" label
        y0 = x @ self.beta

        # Compute search range
        y_range = self.y.max() - self.y.min()
        if y_range == 0:
            y_range = 1.0
        t_max = (self.y.max() - y0) + self.search_range_factor * y_range
        t_min = (self.y.min() - y0) - self.search_range_factor * y_range

        # Run homotopy in both directions
        intervals_pos = self._run_homotopy(x, direction=+1, t_bound=t_max)
        intervals_neg = self._run_homotopy(x, direction=-1, t_bound=t_min)

        # Merge intervals (shift from t-space to y-space)
        all_intervals = []
        for a, b in intervals_neg:
            all_intervals.append((y0 + a, y0 + b))
        # Add t=0 point
        # Check if t=0 is in the prediction set
        if self._t_in_prediction_set(x, 0.0, epsilon):
            all_intervals.append((y0, y0))
        for a, b in intervals_pos:
            all_intervals.append((y0 + a, y0 + b))

        # Merge overlapping/adjacent intervals
        merged = self._merge_intervals(all_intervals)

        if not merged:
            result = self._construct_Gamma(np.nan, np.nan, epsilon)
        elif len(merged) == 1:
            result = self._construct_Gamma(merged[0][0], merged[0][1], epsilon)
        else:
            # Return the smallest enclosing interval (conservative)
            result = self._construct_Gamma(merged[0][0], merged[-1][1], epsilon)

        # Build precomputed for learn_one
        precomputed_dict = None
        if return_update:
            precomputed_dict = {"beta": None}  # Will be set if we can extract from homotopy

        if return_update:
            return result, precomputed_dict
        return result

    def compute_p_value(self, x: NDArray[np.floating[Any]], y: float, smoothed: bool = True, tau: float | None = None) -> float:
        """
        Compute the conformal p-value for (x, y) given current training set.
        """
        x = np.atleast_1d(x).ravel()

        if tau is None and smoothed:
            tau = self.rnd_gen.uniform()

        # Augment training set with (x, y)
        X_aug = np.vstack([self.X, x.reshape(1, -1)])
        y_aug = np.append(self.y, y)

        # Fit Lasso on augmented data
        beta_aug = _solve_lasso(X_aug, y_aug, self.lam, rho=self.rho, warm_start=self.beta)

        # Compute residuals
        residuals = np.abs(y_aug - X_aug @ beta_aug)

        # p-value: fraction of residuals >= residual of test point
        r_test = residuals[-1]
        if smoothed and tau is not None:
            gt = np.sum(residuals > r_test)
            eq = np.sum(residuals == r_test)
            p = (gt + tau * eq) / len(residuals)
        else:
            geq = np.sum(residuals >= r_test)
            p = geq / len(residuals)

        return p

    def _t_in_prediction_set(self, x_new, t, epsilon):
        """Check if a specific t value puts y_{n+1}(t) in the prediction set."""
        n = self.X.shape[0]
        threshold = int(np.ceil((n + 1) * (1 - epsilon)))

        # Augment
        X_aug = np.vstack([self.X, x_new.reshape(1, -1)])
        y_aug = np.append(self.y, x_new @ self.beta + t)

        # Fit Lasso
        beta_t = _solve_lasso(X_aug, y_aug, self.lam, rho=self.rho, warm_start=self.beta)

        # Residuals
        abs_res = np.abs(y_aug - X_aug @ beta_t)

        # Rank of |r_{n+1}| (how many are <= it, in increasing order)
        rank = np.sum(abs_res <= abs_res[-1])
        return rank <= threshold

    def _run_homotopy(self, x_new, direction, t_bound):
        """
        Run the piecewise linear homotopy in one direction.

        Returns a list of (t_start, t_end) intervals that are IN the prediction set.
        """
        n = self.X.shape[0]
        p = self.X.shape[1]
        lam = self.lam
        epsilon = self.epsilon
        threshold = int(np.ceil((n + 1) * (1 - epsilon)))

        # sign of direction
        sign = 1 if direction > 0 else -1
        t_bound_abs = abs(t_bound)

        # Initial state
        beta_k = self.beta.copy()
        J_k = np.where(np.abs(beta_k) > 1e-12)[0]  # active set
        signs_k = np.sign(beta_k[J_k])  # signs on active set

        # Dual variable on inactive set: v = X^T(y - X*beta) evaluated at t=0
        # At t=0 the augmented data is (X; x_new) with label y_{n+1}=x_new^T beta
        # The augmented problem's subgradient:
        # v_j = sum_i (y_i - x_i^T beta) x_{i,j} + (y_{n+1}(0) - x_new^T beta) * x_new_j
        # But y_{n+1}(0) = x_new^T beta, so the second term is 0
        # v = X^T (y - X beta) = X^T r (initial residuals)
        residuals = self.y - self.X @ beta_k
        v_full = self.X.T @ residuals  # subgradient (not including test point contribution at t=0)
        # Note: for the augmented problem at t=0, the n+1-th contribution is 0

        J_c_k = np.setdiff1d(np.arange(p), J_k)  # inactive set
        v_inactive = v_full[J_c_k]

        # Residuals at t=0 for training points
        r_train = residuals.copy()  # r_i(0) = y_i - x_i^T beta for i=1..n
        r_test = 0.0  # r_{n+1}(0) = y_{n+1}(0) - x_new^T beta(0) = 0

        # The Sigma hat used in the paper: (1/n) X^T X
        # But actually the formula uses sum_{i=1}^{n+1} x_i x_i^T = n*Sigma + x_new*x_new^T
        # Let's define Sigma_aug = X^T X + x_new x_new^T (unnormalized)
        XtX = self.X.T @ self.X  # n * Sigma
        XtX_aug = XtX + np.outer(x_new, x_new)  # sum_{i=1}^{n+1} x_i x_i^T

        t_accumulated = 0.0
        intervals_in_set = []

        for _step in range(self.max_homotopy_steps):
            if t_accumulated >= t_bound_abs:
                break

            # Compute eta(k) and gamma(k)
            # eta(k) = (Sigma_aug_{J_k})^{-1} x_{n+1,J_k} / (1 + x_{n+1,J_k}^T (Sigma_aug_{J_k})^{-1} x_{n+1,J_k})
            # But the paper uses n^{-1} * Sigma_hat = (1/n) sum x_i x_i^T, so Sigma_hat_{J_k} = XtX_aug[J_k][:,J_k]/n?
            # Actually re-reading: the paper defines Sigma_hat = (1/n) sum_{i=1}^n x_i x_i^T
            # And the formula has (sum_{i=1}^{n+1} x_{i,J} x_{i,J}^T)^{-1} = (n Sigma_hat_J + x_{n+1,J} x_{n+1,J}^T)^{-1}
            # which equals XtX_aug[J,J]^{-1}

            if len(J_k) == 0:
                # All variables inactive — beta(t) = 0 for all t in this piece
                # r_{n+1}(t) grows linearly with slope 1 (since eta=0)
                # This piece extends until some dual variable hits +/-lambda
                if len(J_c_k) == 0:
                    break
                # gamma(k) = x_{n+1,J_c} (since J is empty, no correction term)
                gamma_k = sign * x_new[J_c_k]

                # Breakpoint: dual variable hits boundary
                dt_dual = np.full(len(J_c_k), np.inf)
                for idx, _j in enumerate(J_c_k):
                    g = gamma_k[idx]
                    if g > 1e-15:
                        dt_dual[idx] = (lam - sign * v_inactive[idx]) / g
                    elif g < -1e-15:
                        dt_dual[idx] = (-lam - sign * v_inactive[idx]) / g
                dt_dual = np.where(dt_dual > 1e-12, dt_dual, np.inf)

                dt_k = np.min(dt_dual)
                dt_k = min(dt_k, t_bound_abs - t_accumulated)

                # In this piece: r_i(t) = r_train_i (constant), r_{n+1}(t) = sign*t (grows)
                # Find sub-intervals in prediction set
                sub_intervals = self._find_intervals_in_piece(
                    r_train,
                    r_test,
                    slopes_train=np.zeros(n),
                    slope_test=sign * 1.0,
                    dt_k=dt_k,
                    t_accumulated=t_accumulated,
                    sign=sign,
                    threshold=threshold,
                    n=n,
                )
                intervals_in_set.extend(sub_intervals)

                # Advance
                t_accumulated += dt_k
                r_test += sign * 1.0 * dt_k
                v_inactive += gamma_k * dt_k

                # Update active set
                if dt_k < t_bound_abs - t_accumulated + 1e-12:
                    entering = J_c_k[np.argmin(dt_dual)]
                    J_k = np.append(J_k, entering)
                    signs_k = np.append(signs_k, np.sign(v_inactive[np.argmin(dt_dual)]))
                    J_c_k = np.setdiff1d(np.arange(p), J_k)
                    v_inactive = v_full[J_c_k]  # will be recomputed below
                continue

            # Normal case: |J_k| > 0
            Sigma_J = XtX_aug[np.ix_(J_k, J_k)] + self.rho * np.eye(len(J_k))
            x_J = x_new[J_k]

            try:
                Sigma_J_inv = np.linalg.inv(Sigma_J)
            except np.linalg.LinAlgError:
                # Singular — cannot continue homotopy
                break

            Sigma_J_inv_x = Sigma_J_inv @ x_J

            # eta(k) = (XtX_aug_J)^{-1} x_{n+1,J}  [paper eq. (5), first form]
            # Note: the Sherman-Morrison form with training-only covariance has a
            # 1/(1 + ...) denominator, but that is already baked into the
            # inverse of XtX_aug.  Do NOT divide by denom again.
            eta_k = sign * Sigma_J_inv_x  # slope of beta_{J_k}(t) w.r.t. |t|

            # gamma(k) for inactive variables  [paper eq. (10), first form]
            if len(J_c_k) > 0:
                Sigma_JcJ = XtX_aug[np.ix_(J_c_k, J_k)]
                gamma_k = sign * (x_new[J_c_k] - Sigma_JcJ @ Sigma_J_inv_x)
            else:
                gamma_k = np.array([])

            # Slopes of residuals:
            # dr_i/d(delta) = -x_{i,J}^T eta_k  for training points
            # dr_{n+1}/d(delta) = sign * (1 - x_J^T eta_unsigned)
            slopes_train = -(self.X[:, J_k] @ eta_k)  # (n,) — dr_i/d(delta_t)
            slope_test = sign * (1.0 - x_J @ Sigma_J_inv_x)  # dr_{n+1}/d(delta_t)

            # Find breakpoint t_{k+1}
            # Primal: beta_j(t_k) + eta_j(k) * dt = 0 for j in J_k
            beta_J = beta_k[J_k]
            dt_primal = np.full(len(J_k), np.inf)
            for idx in range(len(J_k)):
                if abs(eta_k[idx]) > 1e-15:
                    dt = -beta_J[idx] / eta_k[idx]
                    if dt > 1e-12:
                        dt_primal[idx] = dt

            # Dual: |v_j(t_k) + gamma_j * dt| = lambda for j in J_c
            dt_dual = np.full(len(J_c_k), np.inf)
            for idx in range(len(J_c_k)):
                g = gamma_k[idx]
                v_j = v_inactive[idx]
                if abs(g) > 1e-15:
                    # v_j + g*dt = +lambda or -lambda
                    dt1 = (lam - v_j) / g
                    dt2 = (-lam - v_j) / g
                    candidates = []
                    if dt1 > 1e-12:
                        candidates.append(dt1)
                    if dt2 > 1e-12:
                        candidates.append(dt2)
                    if candidates:
                        dt_dual[idx] = min(candidates)

            dt_k = min(
                np.min(dt_primal) if len(dt_primal) > 0 else np.inf,
                np.min(dt_dual) if len(dt_dual) > 0 else np.inf,
            )
            dt_k = min(dt_k, t_bound_abs - t_accumulated)

            if dt_k <= 0 or not np.isfinite(dt_k):
                break

            # Find sub-intervals in this piece that are in the prediction set
            sub_intervals = self._find_intervals_in_piece(
                r_train,
                r_test,
                slopes_train=slopes_train,
                slope_test=slope_test,
                dt_k=dt_k,
                t_accumulated=t_accumulated,
                sign=sign,
                threshold=threshold,
                n=n,
            )
            intervals_in_set.extend(sub_intervals)

            # Advance state
            beta_k[J_k] += eta_k * dt_k
            r_train += slopes_train * dt_k
            r_test += slope_test * dt_k
            if len(J_c_k) > 0:
                v_inactive += gamma_k * dt_k
            t_accumulated += dt_k

            # Update active set based on what hit the boundary
            min_primal = np.min(dt_primal) if len(dt_primal) > 0 else np.inf
            min_dual = np.min(dt_dual) if len(dt_dual) > 0 else np.inf

            if dt_k >= t_bound_abs - (t_accumulated - dt_k):
                break  # Reached search boundary

            if min_primal <= min_dual:
                # A variable leaves the active set
                leaving_idx = np.argmin(dt_primal)
                leaving_var = J_k[leaving_idx]
                beta_k[leaving_var] = 0.0
                J_k = np.delete(J_k, leaving_idx)
                signs_k = np.delete(signs_k, leaving_idx)
            else:
                # A variable enters the active set
                entering_idx = np.argmin(dt_dual)
                entering_var = J_c_k[entering_idx]
                entering_sign = np.sign(v_inactive[entering_idx])
                J_k = np.append(J_k, entering_var)
                signs_k = np.append(signs_k, entering_sign)

            J_c_k = np.setdiff1d(np.arange(p), J_k)
            # Recompute v_inactive for new inactive set
            # v = X_aug^T * r_aug where r_aug includes test point
            # Since we track r_train and r_test:
            v_full = self.X.T @ r_train + x_new * r_test
            v_inactive = v_full[J_c_k]

        return intervals_in_set

    def _find_intervals_in_piece(
        self, r_train, r_test, slopes_train, slope_test, dt_k, t_accumulated, sign, threshold, n
    ):
        """
        Within one homotopy piece of length dt_k, find sub-intervals where
        |r_{n+1}(t)| has rank <= threshold among all |r_i(t)|.

        Residuals are linear in delta_t (local parameter within the piece):
            r_i(delta) = r_train[i] + slopes_train[i] * delta   for i=0..n-1
            r_{n+1}(delta) = r_test + slope_test * delta

        Returns intervals in GLOBAL t-space (t_accumulated + sign*delta mapped to t).
        """
        # Find all crossing points using the JIT-compiled function
        raw_crossings = _compute_crossings(r_train, slopes_train, r_test, slope_test, n, dt_k)

        # Build sorted unique crossings with endpoints
        crossings = [0.0]
        for c in raw_crossings:
            crossings.append(c)
        crossings.append(dt_k)
        crossings = sorted(set(crossings))

        # For each sub-interval, check rank at midpoint
        result_intervals = []
        for idx in range(len(crossings) - 1):
            d_start = crossings[idx]
            d_end = crossings[idx + 1]
            if d_end - d_start < 1e-14:
                continue
            d_mid = (d_start + d_end) / 2

            # Compute residuals at midpoint
            r_i_mid = r_train + slopes_train * d_mid
            r_test_mid = r_test + slope_test * d_mid

            abs_r_i = np.abs(r_i_mid)
            abs_r_test = np.abs(r_test_mid)

            # Rank: number of |r_j| <= |r_{n+1}| (including n+1 itself)
            rank = np.sum(abs_r_i <= abs_r_test) + 1  # +1 for itself

            if rank <= threshold:
                # This sub-interval is in the prediction set
                # Map to global t-space
                t_start = sign * (t_accumulated + d_start)
                t_end = sign * (t_accumulated + d_end)
                if t_start > t_end:
                    t_start, t_end = t_end, t_start
                result_intervals.append((t_start, t_end))

        return result_intervals

    @staticmethod
    def _merge_intervals(intervals):
        """Merge overlapping or adjacent intervals."""
        if not intervals:
            return []
        sorted_intervals = sorted(intervals, key=lambda x: x[0])
        merged = [sorted_intervals[0]]
        for a, b in sorted_intervals[1:]:
            if a <= merged[-1][1] + 1e-12:
                merged[-1] = (merged[-1][0], max(merged[-1][1], b))
            else:
                merged.append((a, b))
        return merged

learn_initial_training_set(X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None

Fit initial Lasso on the training data.

Source code in src/online_cp/regressors.py

def learn_initial_training_set(self, X: NDArray[np.floating[Any]], y: NDArray[np.floating[Any]]) -> None:
    """Fit initial Lasso on the training data."""
    self.X = X.copy()
    self.y = y.copy()
    n, p = X.shape
    self.Sigma = X.T @ X / n

    if self.autotune:
        self._tune_lambda()

    self.beta = _solve_lasso(self.X, self.y, self.lam, rho=self.rho)

learn_one(x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None

Learn a new data point. Updates the training set and refits Lasso.

If precomputed is provided (from predict with return_update=True), uses the cached Lasso solution to avoid refitting.

Source code in src/online_cp/regressors.py

def learn_one(self, x: NDArray[np.floating[Any]], y: float, precomputed: dict[str, Any] | None = None) -> None:
    """
    Learn a new data point. Updates the training set and refits Lasso.

    If precomputed is provided (from predict with return_update=True),
    uses the cached Lasso solution to avoid refitting.
    """
    x = np.atleast_1d(x).ravel()

    # Update training set
    if self.X is None:
        self.X = x.reshape(1, -1)
        self.y = np.array([y])
    else:
        self.X = np.vstack([self.X, x.reshape(1, -1)])
        self.y = np.append(self.y, y)

    # Update sample covariance incrementally: Sigma_new = (n*Sigma_old + x x^T) / (n+1)
    n = self.X.shape[0]
    self.Sigma = ((n - 1) * self.Sigma + np.outer(x, x)) / n

    if precomputed is not None and precomputed.get("beta") is not None:
        self.beta = precomputed["beta"]
    else:
        # Refit from scratch (warm-started from current beta)
        self.beta = _solve_lasso(self.X, self.y, self.lam, rho=self.rho, warm_start=self.beta)

predict(x: NDArray[np.floating[Any]], epsilon: float | NDArray[np.floating[Any]] | None = None, return_update: bool = False) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]

Compute the conformal prediction set at x using the homotopy algorithm.

Returns a ConformalPredictionInterval (if the set is a single interval) or a list of (lower, upper) tuples otherwise.

If epsilon is a list/array, returns a MultiLevelPredictionInterval.

Source code in src/online_cp/regressors.py

def predict(
    self,
    x: NDArray[np.floating[Any]],
    epsilon: float | NDArray[np.floating[Any]] | None = None,
    return_update: bool = False,
) -> ConformalPredictionInterval | MultiLevelPredictionInterval | tuple[ConformalPredictionInterval | MultiLevelPredictionInterval, dict[str, Any]]:
    """
    Compute the conformal prediction set at x using the homotopy algorithm.

    Returns a ConformalPredictionInterval (if the set is a single interval)
    or a list of (lower, upper) tuples otherwise.

    If epsilon is a list/array, returns a MultiLevelPredictionInterval.
    """
    if epsilon is None:
        epsilon = self.epsilon

    # Handle multi-level epsilon
    if hasattr(epsilon, '__iter__'):
        predictions = {}
        for eps in epsilon:
            result = self.predict(x, epsilon=eps, return_update=False)
            predictions[eps] = result
        result = MultiLevelPredictionInterval(predictions)
        if return_update:
            return result, {"beta": None}
        return result

    x = np.atleast_1d(x).ravel()
    n = self.X.shape[0]

    if n < 2:
        result = self._construct_Gamma(-np.inf, np.inf, epsilon)
        if return_update:
            return result, {"beta": None}
        return result

    # y_{n+1}(0) = x_{n+1}^T beta_hat — the "neutral" label
    y0 = x @ self.beta

    # Compute search range
    y_range = self.y.max() - self.y.min()
    if y_range == 0:
        y_range = 1.0
    t_max = (self.y.max() - y0) + self.search_range_factor * y_range
    t_min = (self.y.min() - y0) - self.search_range_factor * y_range

    # Run homotopy in both directions
    intervals_pos = self._run_homotopy(x, direction=+1, t_bound=t_max)
    intervals_neg = self._run_homotopy(x, direction=-1, t_bound=t_min)

    # Merge intervals (shift from t-space to y-space)
    all_intervals = []
    for a, b in intervals_neg:
        all_intervals.append((y0 + a, y0 + b))
    # Add t=0 point
    # Check if t=0 is in the prediction set
    if self._t_in_prediction_set(x, 0.0, epsilon):
        all_intervals.append((y0, y0))
    for a, b in intervals_pos:
        all_intervals.append((y0 + a, y0 + b))

    # Merge overlapping/adjacent intervals
    merged = self._merge_intervals(all_intervals)

    if not merged:
        result = self._construct_Gamma(np.nan, np.nan, epsilon)
    elif len(merged) == 1:
        result = self._construct_Gamma(merged[0][0], merged[0][1], epsilon)
    else:
        # Return the smallest enclosing interval (conservative)
        result = self._construct_Gamma(merged[0][0], merged[-1][1], epsilon)

    # Build precomputed for learn_one
    precomputed_dict = None
    if return_update:
        precomputed_dict = {"beta": None}  # Will be set if we can extract from homotopy

    if return_update:
        return result, precomputed_dict
    return result

compute_p_value(x: NDArray[np.floating[Any]], y: float, smoothed: bool = True, tau: float | None = None) -> float

Compute the conformal p-value for (x, y) given current training set.

Source code in src/online_cp/regressors.py

def compute_p_value(self, x: NDArray[np.floating[Any]], y: float, smoothed: bool = True, tau: float | None = None) -> float:
    """
    Compute the conformal p-value for (x, y) given current training set.
    """
    x = np.atleast_1d(x).ravel()

    if tau is None and smoothed:
        tau = self.rnd_gen.uniform()

    # Augment training set with (x, y)
    X_aug = np.vstack([self.X, x.reshape(1, -1)])
    y_aug = np.append(self.y, y)

    # Fit Lasso on augmented data
    beta_aug = _solve_lasso(X_aug, y_aug, self.lam, rho=self.rho, warm_start=self.beta)

    # Compute residuals
    residuals = np.abs(y_aug - X_aug @ beta_aug)

    # p-value: fraction of residuals >= residual of test point
    r_test = residuals[-1]
    if smoothed and tau is not None:
        gt = np.sum(residuals > r_test)
        eq = np.sum(residuals == r_test)
        p = (gt + tau * eq) / len(residuals)
    else:
        geq = np.sum(residuals >= r_test)
        p = geq / len(residuals)

    return p