code: def mean_alm(
      X=None, tol=1e-14, maxiter=100, sample_weight=None, covmats=None
  ):
      r"""Ando-Li-Mathias (ALM) mean of SPD/HPD matrices.
  
      Return the geometric mean recursively [1]_, generalizing from:
  
      .. math::
          \mathbf{M} = X_1^{\frac{1}{2}} (X_1^{-\frac{1}{2}}X_2^{\frac{1}{2}}
                       X_1^{-\frac{1}{2}})^{\frac{1}{2}} X_1^{\frac{1}{2}}
  
      and requiring a high number of iterations.
      Extremely slow, due to the recursive formulation.
  
      Parameters
      ----------
      X : ndarray, shape (n_matrices, n, n)
          Set of SPD/HPD matrices.
      tol : float, default=10e-14
          The tolerance to stop the gradient descent.
      maxiter : int, default=100
          The maximum number of iterations.
      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)
          ALM mean.
  
      Notes
      -----
      .. versionadded:: 0.3
  
      See Also
      --------
      mean_covariance
  
      References
      ----------
      .. [1] `Geometric Means
          <https://www.sciencedirect.com/science/article/pii/S0024379503008693>`_
          T. Ando, C.-K. Li, and R. Mathias. Linear Algebra and its Applications.
          Volume 385, July 2004, Pages 305-334.
      """
      X = _deprecate_covmats(covmats, X)
      n_matrices, _, _ = X.shape
      sample_weight = check_weights(sample_weight, n_matrices)
  
      if n_matrices == 2:
          alpha = sample_weight[1] / sample_weight[0] / 2
          M = geodesic_riemann(X[0], X[1], alpha=alpha)
          return M
  
      M = X
      M_iter = np.zeros_like(M)
      for _ in range(maxiter):
          for h in range(n_matrices):
              s = np.mod(
                  np.arange(h, h + n_matrices - 1) + 1, n_matrices
              )
              M_iter[h] = mean_alm(M[s], sample_weight=sample_weight[s])
  
          norm_iter = np.linalg.norm(M_iter[0] - M[0], 2)
          norm_c = np.linalg.norm(M[0], 2)
          if (norm_iter / norm_c) < tol:
              break
          M = deepcopy(M_iter)
      else:
          warnings.warn("Convergence not reached")
  
      return M_iter.mean(axis=0)
  
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 3 for calling a global function with example usage: mean_alm(X, tol, maxiter, sample_weight, covmats) and returns: M_iter, M