code: def sample_polytope(
      A: Tensor,
      b: Tensor,
      x0: Tensor,
      n: int = 10000,
      n0: int = 100,
      seed: Optional[int] = None,
  ) -> Tensor:
      r"""
      Hit and run sampler from uniform sampling points from a polytope,
      described via inequality constraints A*x<=b.
  
      Args:
          A: A Tensor describing inequality constraints
              so that all samples satisfy Ax<=b.
          b: A Tensor describing the inequality constraints
              so that all samples satisfy Ax<=b.
          x0: A `d`-dim Tensor representing a starting point of the chain
              satisfying the constraints.
          n: The number of resulting samples kept in the output.
          n0: The number of burn-in samples. The chain will produce
              n+n0 samples but the first n0 samples are not saved.
          seed: The seed for the sampler. If omitted, use a random seed.
  
      Returns:
          (n, d) dim Tensor containing the resulting samples.
      """
      n_tot = n + n0
      seed = (
          seed
          if seed is not None
          else torch.randint(0, 1000000, (1,)).item()
      )
      with manual_seed(seed=seed):
          rands = torch.rand(n_tot, dtype=A.dtype, device=A.device)
  
      # pre-sample samples from hypersphere
      d = x0.size(0)
      # uniform samples from unit ball in d dims
      Rs = sample_hypersphere(
          d=d, n=n_tot, dtype=A.dtype, device=A.device, seed=seed
      ).unsqueeze(-1)
  
      # compute matprods in batch
      ARs = (A @ Rs).squeeze(-1)
      out = torch.empty(n, A.size(-1), dtype=A.dtype, device=A.device)
      x = x0.clone()
      for i, (ar, r, rnd) in enumerate(zip(ARs, Rs, rands)):
          # given x, the next point in the chain is x+alpha*r
          # it also satisfies A(x+alpha*r)<=b which implies A*alpha*r<=b-Ax
          # so alpha<=(b-Ax)/ar for ar>0, and alpha>=(b-Ax)/ar for ar<0.
          # b - A @ x is always >= 0, clamping for numerical tolerances
          w = (b - A @ x).squeeze().clamp(min=0.0) / ar
          pos = w >= 0
          alpha_max = w[pos].min()
          # important to include equality here in cases x is at the boundary
          # of the polytope
          neg = w <= 0
          alpha_min = w[neg].max()
          # alpha~Unif[alpha_min, alpha_max]
          alpha = alpha_min + rnd * (alpha_max - alpha_min)
          x = x + alpha * r
          if i >= n0:  # save samples after burn-in period
              out[i - n0] = x.squeeze()
      return out
  
imports: from __future__ import annotations
  from abc import ABC, abstractmethod
  from contextlib import contextmanager
  from typing import Any, Generator, Iterable, List, Optional, Tuple, TYPE_CHECKING
  import numpy
  import scipy
  import torch
  from botorch.exceptions.errors import BotorchError
  from botorch.sampling.qmc import NormalQMCEngine
  from botorch.utils.transforms import unnormalize
  from scipy.spatial import Delaunay, HalfspaceIntersection
  from torch import LongTensor, Tensor
  from torch.quasirandom import SobolEngine
label_desc: Use image node_id 6 for calling a global function with example usage: sample_polytope(A, b, x0, n, n0, seed) and returns: out