docs: def adaptive_rejection_sampling(logpdf: callable,
                                  a: float, b: float,
                                  domain: Tuple[float, float],
                                  n_samples: int,
                                  random_stream=None):
      """
      Adaptive rejection sampling samples exactly (all samples are i.i.d) and efficiently from any univariate log-concave distribution. The basic idea is to successively determine an envelope of straight-line segments to construct an increasingly accurate approximation of the logarithm.
      It does not require any normalization of the target distribution.
  
      Parameters
      ----------
      logpdf: callable
          Univariate function that computes :math:`log(f(u))`
          for a given :math:`u`, where :math:`f(u)` is proportional
          to the target density to sample from.
  
      a: float
          Lower starting point used to initialize the hulls.
          Must lie in the domain of the logpdf and it
          must hold: :math:`a < b`.
  
      b: float
          Upper starting point used to initialize the hulls.
          Must lie in the domain of the logpdf and it
          must hold: :math:`a < b`.
  
      domain : Tuple[float, float]
          Domain of `logpdf`.
          May be unbounded on either or both sides,
          in which case `(float("-inf"), float("inf"))`
          would be passed.
          If this domain is unbounded to the left,
          the derivative of the logpdf
          for x<= a must be positive.
          If this domain is unbounded to the right the derivative of the logpdf for x>=b
          must be negative.
  
      n_samples: int
          Number of samples to draw.
  
      random_stream : RandomState, optional
          Seeded random number generator object with same interface as a NumPy
          RandomState object. Defaults to `None` in which case a NumPy
          RandomState seeded from `/dev/urandom` if available or the clock if not
          will be used.
  
      Returns
      ----------
      samples : list
          A list of samples drawn from the
          target distribution :math:`f`
          with the given `logpdf`.
  
      Examples
      ----------
      Sampling from a simple gaussian, adaptive rejection sampling style.
      We use the logpdf of a standard gaussian and this small code snippet
      demonstrates that our sample approximation accurately approximates the mean:
  
      >>> from math import isclose
      >>> from numpy import log, exp, mean
      >>> gaussian_logpdf = lambda x, sigma=1: log(exp(-x ** 2 / sigma))
      >>> a, b = -2, 2  # a < b must hold
      >>> domain = (float("-inf"), float("inf"))
      >>> n_samples = 10000
      >>> samples = adaptive_rejection_sampling(logpdf=gaussian_logpdf, a=a, b=b, domain=domain, n_samples=n_samples)
      >>> isclose(mean(samples), 0.0, abs_tol=1e-02)
      True
  
      """
      assert(hasattr(logpdf, "__call__"))
      assert(len(domain) == 2), "Domain must be two-element iterable."
      assert(domain[1] >= domain[0]), "Invalid domain, it must hold: domain[1] >= domain[0]."
      assert(n_samples >= 0), "Number of samples must be >= 0."
  
      if random_stream is None:
          random_stream = RandomState()
  
      if a >= b or isinf(a) or isinf(b) or a < domain[0] or b > domain[1]:
          raise ValueError("invalid a and b")
  
      n_derivative_steps = 1e-3 * (b - a)
  
      S = (a, a + n_derivative_steps, b - n_derivative_steps, b)
  
      if domain[0] == float("-inf"):
          # ensure positive derivative at 'a'
          derivative_sign = sign(logpdf(a + n_derivative_steps) - logpdf(a))
          positive_derivative = derivative_sign > 0
  
          assert(positive_derivative), "derivative at 'a' must be positive, since the domain is unbounded to the left"
  
      if domain[1] == float("inf"):
          # ensure negative derivative at 'b'
          derivative_sign = sign(logpdf(b) - logpdf(b - n_derivative_steps))
          negative_derivative = derivative_sign < 0
  
          assert(negative_derivative), "derivative at 'b' must be negative, since the domain is unbounded to the right"
  
      # initialize a mesh on which to create upper & lower hulls
      n_initial_mesh_points = 3
  
      S = unique(
          (S[0], *(linspace(S[1], S[2], num=n_initial_mesh_points + 2)), S[3])
      )
  
      fS = tuple(logpdf(s) for s in S)
  
      lower_hull, upper_hull = compute_hulls(S=S, fS=fS, domain=domain)
  
      samples = []
  
      while len(samples) < n_samples:
  
          mesh_changed = False
  
          x = sample_upper_hull(upper_hull, random_stream=random_stream)
  
          lh_val, uh_val = evaluate_hulls(x, lower_hull, upper_hull)
  
          U = random_stream.rand()
  
          if log(U) <= lh_val - uh_val:
              # accept u is below lower bound
              samples.append(x)
  
          elif log(U) <= logpdf(x) - uh_val:
              # accept, u is between lower bound and f
              samples.append(x)
  
              mesh_changed = True
          else:
              # reject, u is between f and upper_bound
              mesh_changed = True
  
          if mesh_changed:
              S = sorted([*S, x])
              fS = tuple(logpdf(s) for s in S)
              lower_hull, upper_hull = compute_hulls(S=S, fS=fS, domain=domain)
      return samples