code: def _bvnu_taylor(
      r: Tensor, h: Tensor, k: Tensor, num_points: int = 20
  ) -> Tensor:
      r"""Solves for `P(x > h, y > k)` via Taylor expansion.
  
      Per Section 2.3 of [Genz2004bvnt]_, the bvnu equation (1) may be rewritten as
  
          `L(h, k, r) = L(h, k, s) - s/(2\pi) \int_{0}^{a} f(x) dx \
          f(x) = (1 - x^2){-1/2} e^{-0.5 ((h - sk)/ x)^2} e^{-shk/(1 + (1 - x^2)^{1/2})},`
  
      where `s = sign(r)` and `a = sqrt(1 - r^{2})`. The term `L(h, k, s)` is analytic.
      The second integral is approximated via Taylor expansion. See Sections 2.3 and
      2.4 of [Genz2004bvnt]_.
      """
      _0, _1, _ni2, _i2pi, _sq2pi = get_constants_like(
          values=(0, 1, -0.5, _inv_2pi, _sqrt_2pi), ref=r
      )
  
      x, w = leggauss(num_points, dtype=r.dtype, device=r.device)
      x = x + _1
  
      s = (
          get_constants_like(2, r) * (r > _0).to(r) - _1
      )  # sign of `r` where sign(0) := 1
      sk = s * k
      skh = sk * h
      comp_r2 = _1 - r.square()
  
      a = comp_r2.clip(min=0).sqrt()
      b = safe_sub(h, sk)
      b2 = b.square()
      c = get_constants_like(1 / 8, r) * (
          get_constants_like(4, r) - skh
      )
      d = get_constants_like(1 / 80, r) * (
          get_constants_like(12, r) - skh
      )
  
      # ---- Solve for `L(h, k, s)`
      int_from_0_to_s = case_dispatcher(
          out=torch.empty_like(r),
          cases=[
              (
                  lambda: r > _0,
                  lambda mask: Phi(-torch.maximum(h[mask], k[mask])),
              )
          ],
          default=lambda mask: (Phi(sk[mask]) - Phi(h[mask])).clip(
              min=_0
          ),
      )
  
      # ---- Solve for `s/(2\pi) \int_{0}^{a} f(x) dx`
      # Analytic part
      _a0 = _ni2 * (safe_div(b2, comp_r2) + skh)
      _a1 = c * get_constants_like(1 / 3, r) * (_1 - d * b2)
      _a2 = _1 - b2 * _a1
      abs_b = b.abs()
      analytic_part = torch.subtract(  # analytic part of solution
          a
          * (_a2 + comp_r2 * _a1 + c * d * comp_r2.square())
          * safe_exp(_a0),
          _sq2pi
          * Phi(safe_div(-abs_b, a))
          * abs_b
          * _a2
          * safe_exp(_ni2 * skh),
      )
  
      # Quadrature part
      _b2 = b2.unsqueeze(-1)
      _skh = skh.unsqueeze(-1)
      _q0 = (
          get_constants_like(0.25, r)
          * comp_r2.unsqueeze(-1)
          * x.square()
      )
      _q1 = (_1 - _q0).sqrt()
      _q2 = _ni2 * (_b2 / _q0 + _skh)
  
      _b2 = b2.unsqueeze(-1)
      _c = c.unsqueeze(-1)
      _d = d.unsqueeze(-1)
      vals = (_ni2 * (_b2 / _q0 + _skh)).exp() * torch.subtract(
          _1 + _c * _q0 * (_1 + get_constants_like(5, r) * _d * _q0),
          safe_exp(_ni2 * _q0 / (_1 + _q1).square() * _skh) / _q1,
      )
      mask = _q2 > get_constants_like(-100, r)
      if not mask.all():
          vals[~mask] = _0
      quadrature_part = _ni2 * a * (vals @ w)
  
      # Return `P(x > h, y > k)`
      int_from_0_to_a = _i2pi * s * (analytic_part + quadrature_part)
      return (int_from_0_to_s - int_from_0_to_a).clip(min=_0, max=_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 4 for calling a global function with example usage: _bvnu_taylor(r, h, k, num_points) and returns: unknown