docs: def tensor_jacobian_product(fun, argnum=0):
      """Builds a function that returns the exact tensor-Jacobian product, that
      is the Jacobian matrix left-multiplied by tensor. The returned function
      has arguments (*args, tensor, **kwargs)."""
      def vector_dot_fun(*args, **kwargs):
          args, vector = args[:-1], args[-1]
          return np.tensordot(vector, fun(*args, **kwargs), axes=np.ndim(vector))
      return jacobian(vector_dot_fun, argnum)