code: def mean_power(
      X=None,
      p=None,
      *,
      sample_weight=None,
      zeta=10e-10,
      maxiter=100,
      covmats=None,
  ):
      r"""Power mean of SPD/HPD matrices.
  
      Power mean of order p is the solution of [1]_ [2]_:
  
      .. math::
          \mathbf{M} = \sum_i w_i \ \mathbf{M} \sharp_p \mathbf{X}_i
  
      where :math:`\mathbf{A} \sharp_p \mathbf{B}` is the geodesic between
      matrices :math:`\mathbf{A}` and :math:`\mathbf{B}`.
  
      Parameters
      ----------
      X : ndarray, shape (n_matrices, n, n)
          Set of SPD/HPD matrices.
      p : float
          Exponent, in [-1,+1]. For p=0, it returns
          :func:`pyriemann.utils.mean.mean_riemann`.
      sample_weight : None | ndarray, shape (n_matrices,), default=None
          Weights for each matrix. If None, it uses equal weights.
      zeta : float, default=10e-10
          Stopping criterion.
      maxiter : int, default=100
          The maximum number of iterations.
  
      Returns
      -------
      M : ndarray, shape (n, n)
          Power mean.
  
      Notes
      -----
      .. versionadded:: 0.3
  
      See Also
      --------
      mean_covariance
  
      References
      ----------
      .. [1] `Matrix Power means and the Karcher mean
          <https://www.sciencedirect.com/science/article/pii/S0022123611004101>`_
          Y. Lim and M. Palfia. Journal of Functional Analysis, Volume 262,
          Issue 4, 15 February 2012, Pages 1498-1514.
      .. [2] `Fixed Point Algorithms for Estimating Power Means of Positive
          Definite Matrices
          <https://hal.archives-ouvertes.fr/hal-01500514>`_
          M. Congedo, A. Barachant, and R. Bhatia. IEEE Transactions on Signal
          Processing, Volume 65, Issue 9, pp.2211-2220, May 2017
      """
      X = _deprecate_covmats(covmats, X)
      if p is None:
          raise ValueError("Input p can not be None")
      if not isinstance(p, (int, float)):
          raise ValueError(
              "Power mean only defined for a scalar exponent"
          )
      if p < -1 or 1 < p:
          raise ValueError("Exponent must be in [-1,+1]")
  
      if p == 1:
          return mean_euclid(X, sample_weight=sample_weight)
      elif p == 0:
          return mean_riemann(X, sample_weight=sample_weight)
      elif p == -1:
          return mean_harmonic(X, sample_weight=sample_weight)
  
      n_matrices, n, _ = X.shape
      sample_weight = check_weights(sample_weight, n_matrices)
      phi = 0.375 / np.abs(p)
  
      G = np.einsum("a,abc->bc", sample_weight, powm(X, p))
      if p > 0:
          K = invsqrtm(G)
      else:
          K = sqrtm(G)
  
      eye_n, sqrt_n = np.eye(n), np.sqrt(n)
      crit = 10 * zeta
      for _ in range(maxiter):
          H = np.einsum(
              "a,abc->bc",
              sample_weight,
              powm(K @ powm(X, np.sign(p)) @ K.conj().T, np.abs(p)),
          )
          K = powm(H, -phi) @ K
  
          crit = np.linalg.norm(H - eye_n) / sqrt_n
          if crit <= zeta:
              break
      else:
          warnings.warn("Convergence not reached")
  
      M = K.conj().T @ K
      if p > 0:
          M = np.linalg.inv(M)
  
      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 10 for calling a global function with example usage: mean_power(X, p) and returns: M, mean_euclid, mean_riemann, mean_harmonic