code: def _fit_single(
      X,
      y=None,
      n_clusters=2,
      init="random",
      random_state=None,
      metric="riemann",
      max_iter=100,
      tol=1e-4,
      n_jobs=1,
  ):
      """helper to fit a single run of centroid."""
      # init random state if provided
      mdm = MDM(metric=metric, n_jobs=n_jobs)
      mdm.metric_mean, mdm.metric_dist = _check_metric(metric)
      squared_norms = [np.linalg.norm(x, ord="fro") ** 2 for x in X]
      mdm.covmeans_ = _init_centroids(
          X,
          n_clusters,
          init,
          random_state=random_state,
          x_squared_norms=squared_norms,
      )
      if y is not None:
          mdm.classes_ = np.unique(y)
      else:
          mdm.classes_ = np.arange(n_clusters)
  
      labels = mdm.predict(X)
      k = 0
      while True:
          old_labels = labels.copy()
          mdm.fit(X, old_labels)
          dist = mdm._predict_distances(X)
          labels = mdm.classes_[dist.argmin(axis=1)]
          k += 1
          if (k > max_iter) | (
              np.mean(labels == old_labels) > (1 - tol)
          ):
              break
      inertia = sum(
          [
              sum(dist[labels == mdm.classes_[i], i])
              for i in range(len(mdm.classes_))
          ]
      )
      return labels, inertia, mdm
  
imports: import numpy
  from joblib import Parallel, delayed
  from scipy.stats import norm, chi2
  import sklearn
  from sklearn.base import BaseEstimator, ClassifierMixin, TransformerMixin, ClusterMixin, clone
  from sklearn.cluster import KMeans
  from .classification import MDM, _check_metric
  from .utils.mean import mean_covariance
  from .utils.geodesic import geodesic
label_desc: Use image node_id 2 for calling a global function with example usage: _fit_single(X, y, n_clusters, init, random_state, metric, max_iter, tol, n_jobs) and returns: labels,  inertia,  mdm