def conformal_interval(
        self,
        bootstrap_indices: Array,
        bootstrap_train_preds: Array,
        bootstrap_test_preds: Array,
        train_targets: Array,
        error: float,
        return_residuals: bool = False,
    ) -> 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.

        Parameters
        ----------
        bootstrap_indices: Array
            The indices, randomly sampled with replacement, of the training data in the time series used to train the
            model. The first dimension is over the different samples of indices. The second dimension contain the
            data points for each sample, which must be as many as thnumber of data points in the time series used for
            training. A simple way of obtaining indices randomly sampled with replacement is
            :code:`numpy.random.choice(T, size=(B, T))`,
            where :code:`T` is the number of training points in the time series,
            and :code:`B` is the number of bootstrap samples.
            It is the user job to make sure that the models are trained upon the data corresponding to the random
            indices.
        bootstrap_train_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at each of the training data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the training inputs.
            There may be a third dimension, corresponding to the dimensionality of the predictions,
            but if so this must be one.
        bootstrap_test_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at each of the test data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the test inputs. There may be a third dimension,
            corresponding to the dimensionality of the predictions, but if so this must be one.
        train_targets: Array
            The target variables of the training data points in the time series.
        error: float
            The desired coverage error. This must be a scalar between 0 and 1, extremes included.
        return_residuals: bool
            If True, return the residual errors computed over the training data. These are used in
            :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval_from_residuals`.

        Returns
        -------
        Union[Array, Tuple[Array, Array]]
            The conformal intervals. The two components of the second dimension correspond to the left and right
            interval bounds. If :code:`return_residuals` is set to True, then it returns also the residuals computed on
            the training set.
        """
        n_bootstraps, n_train_times = bootstrap_indices.shape
        if (
            jnp.min(bootstrap_indices) < 0
            or jnp.max(bootstrap_indices) >= n_train_times
        ):
            raise ValueError(
                f"All elements f `bootstrap_indices` must be integers from 0 to {n_train_times - 1} "
                f"corresponding to the indices of the data used for training in each of the bootstrap "
                f"samples."
            )
        if bootstrap_train_preds.shape[0] != n_bootstraps:
            raise ValueError(
                "The first dimension of `bootstrap_train_preds` and `bootstrap_indices` corresponds to "
                "the number of bootstrap samples, and must have the same size. However, "
                f"{bootstrap_train_preds.shape[0]} and {n_bootstraps} were found, respectively."
            )
        if bootstrap_train_preds.shape[1] != n_train_times:
            raise ValueError(
                "The second dimension of `bootstrap_train_preds` and `bootstrap_indices` corresponds to "
                "the number of data points in the time series used for training, "
                "and must have the same size. "
                f"However, {bootstrap_train_preds.shape[1]} and {n_train_times} were found, respectively."
            )
        if bootstrap_test_preds.shape[0] != n_bootstraps:
            raise ValueError(
                "The first dimension of `bootstrap_test_preds` and `bootstrap_indices` corresponds to the "
                "number of bootstrap samples, and must have the same size. However, "
                f"{bootstrap_test_preds.shape[0]} and {n_bootstraps} were found, respectively."
            )
        if bootstrap_train_preds.ndim == 3:
            if bootstrap_train_preds.shape[2] == 1:
                bootstrap_train_preds = bootstrap_train_preds.squeeze(2)
            else:
                raise ValueError(
                    "This method is supported only for scalar model predictions. However, `bootstrap_train_preds` has "
                    "third dimension greater than 1."
                )
        if bootstrap_test_preds.ndim == 3:
            if bootstrap_test_preds.shape[2] == 1:
                bootstrap_test_preds = bootstrap_test_preds.squeeze(2)
            else:
                raise ValueError(
                    "This method is supported only for scalar model predictions. However, `bootstrap_test_preds` has "
                    "third dimension greater than 1."
                )
        if train_targets.shape[0] != bootstrap_train_preds.shape[1]:
            raise ValueError(
                "The first dimension of `train_targets` and the second dimension of "
                "`bootstrap_train_preds` correspond to the number of data points in the time series used "
                "for training, and must have the same size. However, "
                f"{train_targets.shape[0]} and {bootstrap_train_preds.shape[1]} were found, respectively."
            )
        if train_targets.ndim == 2:
            if train_targets.shape[1] == 1:
                train_targets = train_targets.squeeze(1)
            else:
                raise ValueError(
                    "This method is supported only for scalar target variables. However, `train_targets` has "
                    "second dimension greater than 1."
                )

        in_bootstrap_indices = np.zeros((n_bootstraps, n_train_times), dtype=bool)
        np.put_along_axis(in_bootstrap_indices, bootstrap_indices, values=1, axis=1)
        aggr_bootstrap_test_preds = np.zeros(
            (n_train_times,) + bootstrap_test_preds.shape[1:]
        )
        train_residuals = np.zeros((n_train_times,) + train_targets.shape[1:])

        for t in range(n_train_times):
            which_bootstraps = np.where(~(in_bootstrap_indices[:, t]))[0]
            if len(which_bootstraps) > 0:
                aggr_bootstrap_train_pred = self.aggregation_fun(
                    bootstrap_train_preds[which_bootstraps, t]
                )
                train_residuals[t] = np.abs(
                    train_targets[t] - aggr_bootstrap_train_pred
                )
                aggr_bootstrap_test_preds[t] = self.aggregation_fun(
                    bootstrap_test_preds[which_bootstraps]
                )
            else:
                train_residuals[t] = np.abs(train_targets[t])

        test_quantiles = jnp.quantile(aggr_bootstrap_test_preds, q=1 - error, axis=0)
        residuals_quantile = jnp.quantile(train_residuals, q=1 - error, axis=0)

        left = test_quantiles - residuals_quantile
        right = test_quantiles + residuals_quantile

        conformal_intervals = jnp.array(list(zip(left, right)))
        if not return_residuals:
            return conformal_intervals
        return conformal_intervals, 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.

        Parameters
        ----------
        bootstrap_indices: Array
            The indices, randomly sampled with replacement, of the training data in the time series used to train the
            model. The first dimension is over the different samples of indices. The second dimension contain the
            data points for each sample, which must be as many as thnumber of data points in the time series used for
            training. A simple way of obtaining indices randomly sampled with replacement is
            :code:`numpy.random.choice(T, size=(B, T))`,
            where :code:`T` is the number of training points in the time series,
            and :code:`B` is the number of bootstrap samples.
            It is the user job to make sure that the models are trained upon the data corresponding to the random
            indices.
        bootstrap_train_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at each of the training data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the training inputs.
            There may be a third dimension, corresponding to the dimensionality of the predictions,
            but if so this must be one.
        bootstrap_test_preds: Array
            Model predictions for each of the bootstrap samples of data of the time series used for training the model,
            evaluated at each of the test data inputs of the time series. The first dimension is over the
            different bootstrap samples. The second dimensions is over the test inputs. There may be a third dimension,
            corresponding to the dimensionality of the predictions, but if so this must be one.
        train_targets: Array
            The target variables of the training data points in the time series.
        error: float
            The desired coverage error. This must be a scalar between 0 and 1, extremes included.
        return_residuals: bool
            If True, return the residual errors computed over the training data. These are used in
            :meth:`~fortuna.conformal.regression.enbpi.EnbPI.conformal_interval_from_residuals`.

        Returns
        -------
        Union[Array, Tuple[Array, Array]]
            The conformal intervals. The two components of the second dimension correspond to the left and right
            interval bounds. If :code:`return_residuals` is set to True, then it returns also the residuals computed on
            the training set.
        
fortuna/conformal/regression/enbpi.py