def score(
        self,
        val_lower_bounds: Array,
        val_upper_bounds: Array,
        val_targets: Array,
    ) -> jnp.ndarray:
        """
        Compute score function.

        Parameters
        ----------
        val_lower_bounds: Array
            Interval lower bounds computed on a validation set.
        val_upper_bounds: Array
            Interval upper bounds computed on a validation set.
        val_targets: Array
            A two-dimensional array of validation target variables.

        Returns
        -------
        jnp.ndarray
            The conformal scores.
        """
        if val_lower_bounds.shape != val_upper_bounds.shape:
            raise ValueError(
                """The shapes of `val_lower_bounds` and `val_upper_bounds` must be the same, but shapes {} and {}
                were given, respectively.""".format(
                    val_lower_bounds.shape, val_upper_bounds.shape
                )
            )
        if val_lower_bounds.ndim != 1:
            raise ValueError(
                "`val_lower_bounds` and `val_upper_bounds` must be one-dimensional arrays."
            )
        if val_targets.shape[1] != 1:
            raise ValueError(
                """The second dimension of `val_targets` must have only one component."""
            )
        val_targets = val_targets.squeeze(1)
        return jnp.maximum(
            val_lower_bounds - val_targets, val_targets - val_upper_bounds
        )

        Compute score function.

        Parameters
        ----------
        val_lower_bounds: Array
            Interval lower bounds computed on a validation set.
        val_upper_bounds: Array
            Interval upper bounds computed on a validation set.
        val_targets: Array
            A two-dimensional array of validation target variables.

        Returns
        -------
        jnp.ndarray
            The conformal scores.
        
fortuna/conformal/regression/quantile.py