code: def mvhermgauss(H: int, D: int):
      
      gh_x, gh_w = hermgauss(H)
      x = np.array(list(itertools.product(*(gh_x,) * D)))  # H**DxD
      w = np.prod(np.array(list(itertools.product(*(gh_w,) * D))), 1)  # H**D
      return x, w
docs: Return the evaluation locations 'xn', and weights 'wn' for a multivariate
      Gauss-Hermite quadrature.
  
      The outputs can be used to approximate the following type of integral:
      int exp(-x)*f(x) dx ~ sum_i w[i,:]*f(x[i,:])
  
      :param H: Number of Gauss-Hermite evaluation points.
      :param D: Number of input dimensions. Needs to be known at call-time.
      :return: eval_locations 'x' (H**DxD), weights 'w' (H**D)