code: def bvn(
      r: Tensor, xl: Tensor, yl: Tensor, xu: Tensor, yu: Tensor
  ) -> Tensor:
      r"""A function for computing bivariate normal probabilities.
  
      Calculates `P(xl < x < xu, yl < y < yu)` where `x` and `y` are bivariate normal with
      unit variance and correlation coefficient `r`. See Section 2.4 of [Genz2004bvnt]_.
  
      This method uses a sign flip trick to improve numerical performance. Many of `bvnu`s
      internal branches rely on evaluations `Phi(-bound)`. For `a < b < 0`, the term
      `Phi(-a) - Phi(-b)` goes to zero faster than `Phi(b) - Phi(a)` because
      `finfo(dtype).epsneg` is typically much larger than `finfo(dtype).tiny`. In these
      cases, flipping the sign can prevent situations where `bvnu(...) - bvnu(...)` would
      otherwise be zero due to round-off error.
  
      Args:
          r: Tensor of correlation coefficients.
          xl: Tensor of lower bounds for `x`, same shape as `r`.
          yl: Tensor of lower bounds for `y`, same shape as `r`.
          xu: Tensor of upper bounds for `x`, same shape as `r`.
          yu: Tensor of upper bounds for `y`, same shape as `r`.
  
      Returns:
          Tensor of probabilities `P(xl < x < xu, yl < y < yu)`.
  
      """
      if not (r.shape == xl.shape == xu.shape == yl.shape == yu.shape):
          raise UnsupportedError(
              "Arguments to `bvn` must have the same shape."
          )
  
      # Sign flip trick
      _0, _1, _2 = get_constants_like(values=(0, 1, 2), ref=r)
      flip_x = xl.abs() > xu  # is xl more negative than xu is positive?
      flip_y = yl.abs() > yu
      flip = (flip_x & (~flip_y | yu.isinf())) | (
          flip_y & (~flip_x | xu.isinf())
      )
      if (
          flip.any()
      ):  # symmetric calls to `bvnu` below makes swapping bounds unnecessary
          sign = _1 - _2 * flip.to(dtype=r.dtype)
          xl = sign * xl  # becomes `-xu` if flipped
          xu = sign * xu  # becomes `-xl`
          yl = sign * yl  # becomes `-yu`
          yu = sign * yu  # becomes `-yl`
  
      p = (
          bvnu(r, xl, yl)
          - bvnu(r, xu, yl)
          - bvnu(r, xl, yu)
          + bvnu(r, xu, yu)
      )
      return p.clip(_0, _1)
  
file_path: botorch/botorch/utils/probability/bvn.py
imports: from __future__ import annotations
  from math import pi
  from typing import Optional, Tuple
  import torch
  from botorch.exceptions import UnsupportedError
  from botorch.utils.probability.utils import case_dispatcher, get_constants_like, leggauss, ndtr, phi, STANDARDIZED_RANGE
  from botorch.utils.safe_math import div, exp, mul, sub
  from torch import Tensor
label_desc: Use image node_id 1 for calling a global function with example usage: bvn(r, xl, yl, xu, yu) and returns: p