def conformal_interval_from_residuals(
        self,
        train_residuals: Array,
        bootstrap_new_train_preds: Array,
        bootstrap_new_test_preds: Array,
        new_train_targets: Array,
        error: float,
    ) -> Union[Array, Tuple[Array, Array]]:
        """
        Compute a coverage interval for each of the test inputs of the time series, at the desired coverage error.
        This is supported only for one-dimensional target variables, and for one times series at the time.
        This method assumes the residuals over the training sets have already been computed via
        :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval`
        using the flag :code:`return_residuals=True`.
        If so, this method takes in the predictions for each bootstrap samples over the new incoming training and test
        data points, and computes conformal intervals without having to retrain the model. Compared to predicting
        conformal intervals over all test data points at once with
        :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval`,
        this method exploits incoming data points as a form of online feedback,
        and improves the conformal intervals accordingly.

        Parameters
        ----------
        train_residuals: Array
            Residuals over the training data points computed via
            :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval` using the flag
            :code:`return_residuals=True`.
        bootstrap_new_train_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at the new incoming training data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the new training inputs.
            There may be a third dimension, corresponding to the dimensionality of the predictions,
            but if so this must be one. As an example, suppose that the bootstrap has been done over data points from
            1 to T, and the model was trained over each bootstrap sample. A prediction was then made on a batch of
            data points from T+1 to T+T1. After the data from T+1 to T+T1 is observed, this will be used as the new
            batch of training data.
        bootstrap_new_test_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at the new incoming test data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the new test inputs.
            There may be a third dimension, corresponding to the dimensionality of the predictions,
            but if so this must be one. As an example, suppose that the bootstrap has been done over data points from
            1 to T, and the model was trained over each bootstrap sample. A prediction was then made on a batch of
            data points from T+1 to T+T1. After the data from T+1 to T+T1 is observed, data from T+T1+1 to T+2*T1 may
            be taken as the new batch of test data.
        new_train_targets: Array
            The target variables of the new batch of training data points in the time series.
        error: float
            The desired coverage error. This must be a scalar between 0 and 1, extremes included.
        Returns
        -------
        Union[Array, Tuple[Array, Array]]
            - The conformal intervals. The two components of the second dimension correspond to the left and right
            interval bounds.
            - A new set of residuals, which includes the residuals computed on the new batch of training data points.
            The number of stored training residuals is kept constant by eliminating the oldest ones.
        """
        if bootstrap_new_train_preds.shape[0] != bootstrap_new_test_preds.shape[0]:
            raise ValueError(
                "The first dimensions of `bootstrap_new_train_preds` and `bootstrap_new_test_preds` "
                "correspond to the number of bootstrap samples, "
                "and must have the same size. However, "
                f"{bootstrap_new_train_preds.shape[0]} and "
                f"{bootstrap_new_test_preds.shape[0]} were found, respectively."
            )
        if bootstrap_new_train_preds.ndim == 3:
            if bootstrap_new_train_preds.shape[2] == 1:
                bootstrap_new_train_preds = bootstrap_new_train_preds.squeeze(2)
            else:
                raise ValueError(
                    "This method is supported only for scalar model predictions. However, "
                    "`bootstrap_new_train_preds` has third dimension greater than 1."
                )
        if bootstrap_new_test_preds.ndim == 3:
            if bootstrap_new_test_preds.shape[2] == 1:
                bootstrap_new_test_preds = bootstrap_new_test_preds.squeeze(2)
            else:
                raise ValueError(
                    "This method is supported only for scalar model predictions. "
                    "However, `bootstrap_new_test_preds` has third dimension greater than 1."
                )
        if new_train_targets.shape[0] != bootstrap_new_train_preds.shape[1]:
            raise ValueError(
                "The first dimension of `new_train_targets` and the second dimension of "
                "`bootstrap_new_train_preds` correspond to the number of data points in the time series "
                "used for training, and must have the same size. However, "
                f"{new_train_targets.shape[0]} and {bootstrap_new_train_preds.shape[1]} were found, "
                f"respectively."
            )
        if new_train_targets.ndim == 2:
            if new_train_targets.shape[1] == 1:
                new_train_targets = new_train_targets.squeeze(1)
            else:
                raise ValueError(
                    "This method is supported only for scalar target variables. However, `new_train_targets` has "
                    "second dimension greater than 1."
                )

        train_residuals[:-1] = train_residuals[1:]
        train_residuals[-1] = np.abs(
            new_train_targets - self.aggregation_fun(bootstrap_new_train_preds)
        )
        aggr_test_preds = self.aggregation_fun(bootstrap_new_test_preds)
        residuals_quantile = jnp.quantile(train_residuals, q=1 - error, axis=0)

        left = aggr_test_preds - residuals_quantile
        right = aggr_test_preds + residuals_quantile

        return jnp.array(list(zip(left, right))), train_residuals

        Compute a coverage interval for each of the test inputs of the time series, at the desired coverage error.
        This is supported only for one-dimensional target variables, and for one times series at the time.
        This method assumes the residuals over the training sets have already been computed via
        :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval`
        using the flag :code:`return_residuals=True`.
        If so, this method takes in the predictions for each bootstrap samples over the new incoming training and test
        data points, and computes conformal intervals without having to retrain the model. Compared to predicting
        conformal intervals over all test data points at once with
        :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval`,
        this method exploits incoming data points as a form of online feedback,
        and improves the conformal intervals accordingly.

        Parameters
        ----------
        train_residuals: Array
            Residuals over the training data points computed via
            :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval` using the flag
            :code:`return_residuals=True`.
        bootstrap_new_train_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at the new incoming training data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the new training inputs.
            There may be a third dimension, corresponding to the dimensionality of the predictions,
            but if so this must be one. As an example, suppose that the bootstrap has been done over data points from
            1 to T, and the model was trained over each bootstrap sample. A prediction was then made on a batch of
            data points from T+1 to T+T1. After the data from T+1 to T+T1 is observed, this will be used as the new
            batch of training data.
        bootstrap_new_test_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at the new incoming test data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the new test inputs.
            There may be a third dimension, corresponding to the dimensionality of the predictions,
            but if so this must be one. As an example, suppose that the bootstrap has been done over data points from
            1 to T, and the model was trained over each bootstrap sample. A prediction was then made on a batch of
            data points from T+1 to T+T1. After the data from T+1 to T+T1 is observed, data from T+T1+1 to T+2*T1 may
            be taken as the new batch of test data.
        new_train_targets: Array
            The target variables of the new batch of training data points in the time series.
        error: float
            The desired coverage error. This must be a scalar between 0 and 1, extremes included.
        Returns
        -------
        Union[Array, Tuple[Array, Array]]
            - The conformal intervals. The two components of the second dimension correspond to the left and right
            interval bounds.
            - A new set of residuals, which includes the residuals computed on the new batch of training data points.
            The number of stored training residuals is kept constant by eliminating the oldest ones.
        
fortuna/conformal/regression/enbpi.py