code: def mean_riemann(
      X=None,
      tol=10e-9,
      maxiter=50,
      init=None,
      sample_weight=None,
      covmats=None,
  ):
      r"""Mean of SPD/HPD matrices according to the Riemannian metric.
  
      The affine-invariant Riemannian mean minimizes the sum of squared
      affine-invariant Riemannian distances :math:`d_R` to all SPD/HPD matrices
      [1]_ [2]_:
  
      .. math::
           \arg \min_{\mathbf{M}} \sum_i w_i \ d_R (\mathbf{M}, \mathbf{X}_i)^2
  
      Parameters
      ----------
      X : ndarray, shape (n_matrices, n, n)
          Set of SPD/HPD matrices.
      tol : float, default=10e-9
          The tolerance to stop the gradient descent.
      maxiter : int, default=50
          The maximum number of iterations.
      init : None | ndarray, shape (n, n), default=None
          A SPD/HPD matrix used to initialize the gradient descent.
          If None, the weighted Euclidean mean is used.
      sample_weight : None | ndarray, shape (n_matrices,), default=None
          Weights for each matrix. If None, it uses equal weights.
  
      Returns
      -------
      M : ndarray, shape (n, n)
          Affine-invariant Riemannian mean.
  
      See Also
      --------
      mean_covariance
  
      References
      ----------
      .. [1] `Principal geodesic analysis for the study of nonlinear statistics
          of shape
          <https://ieeexplore.ieee.org/document/1318725>`_
          P.T. Fletcher, C. Lu, S. M. Pizer, S. Joshi.
          IEEE Trans Med Imaging, 2004, 23(8), pp. 995-1005
      .. [2] `A differential geometric approach to the geometric mean of
          symmetric positive-definite matrices
          <https://epubs.siam.org/doi/10.1137/S0895479803436937>`_
          M. Moakher. SIAM J Matrix Anal Appl, 2005, 26 (3), pp. 735-747
      """
      X = _deprecate_covmats(covmats, X)
      n_matrices, _, _ = X.shape
      sample_weight = check_weights(sample_weight, n_matrices)
      if init is None:
          M = mean_euclid(X, sample_weight=sample_weight)
      else:
          M = init
  
      nu = 1.0
      tau = np.finfo(np.float64).max
      crit = np.finfo(np.float64).max
      for _ in range(maxiter):
          M12, Mm12 = sqrtm(M), invsqrtm(M)
          J = np.einsum(
              "a,abc->bc", sample_weight, logm(Mm12 @ X @ Mm12)
          )
          M = M12 @ expm(nu * J) @ M12
  
          crit = np.linalg.norm(J, ord="fro")
          h = nu * crit
          if h < tau:
              nu = 0.95 * nu
              tau = h
          else:
              nu = 0.5 * nu
          if crit <= tol or nu <= tol:
              break
      else:
          warnings.warn("Convergence not reached")
  
      return M
  
imports: from copy import deepcopy
  import numpy
  import warnings
  from .ajd import ajd_pham
  from .base import sqrtm, invsqrtm, logm, expm, powm
  from .distance import distance_riemann
  from .geodesic import geodesic_riemann
  from .utils import check_weights
label_desc: Use image node_id 11 for calling a global function with example usage: mean_riemann(X, tol, maxiter, init, sample_weight, covmats) and returns: M