code: def VB_diarization(
      X,
      m,
      invSigma,
      w,
      V,
      pi=None,
      gamma=None,
      maxSpeakers=10,
      maxIters=10,
      epsilon=1e-4,
      loopProb=0.99,
      statScale=1.0,
      alphaQInit=1.0,
      downsample=None,
      VtinvSigmaV=None,
      ref=None,
      plot=False,
      sparsityThr=0.001,
      llScale=1.0,
      minDur=1,
      Fa=1.0,
      Fb=1.0,
  ):
      """
      This a generalized version of speaker diarization described in:
  
      Diez. M., Burget. L., Landini. F., Cernocky. J.
      Analysis of Speaker Diarization based on Bayesian HMM with Eigenvoice Priors
  
      Variable names and equation numbers refer to those used the paper
  
      Inputs:
      X  - T x D array, where columns are D dimensional feature vectors for T frames
      m  - C x D array of GMM component means
      invSigma - C x D array of GMM component inverse covariance matrix diagonals
      w  - C dimensional column vector of GMM component weights
      V  - R x C x D array of eigenvoices
      maxSpeakers - maximum number of speakers expected in the utterance
      maxIters    - maximum number of algorithm iterations
      epsilon     - stop iterating, if obj. fun. improvement is less than epsilon
      loopProb    - probability of not switching speakers between frames
      statScale   - deprecated, use Fa instead
      Fa          - scale sufficient statiscits collected using UBM
      Fb          - speaker regularization coefficient Fb (controls final # of speaker)
      llScale     - scale UBM likelihood (i.e. llScale < 1.0 make atribution of
                    frames to UBM componets more uncertain)
      sparsityThr - set occupations smaller that this threshold to 0.0 (saves memory
                    as the posteriors are represented by sparse matrix)
      alphaQInit  - Dirichlet concentration parameter for initializing gamma
      downsample  - perform diarization on input downsampled by this factor
      VtinvSigmaV       - C x (R**2+R)/2 matrix normally calculated by VB_diarization when
                    VtinvSigmaV is None. However, it can be pre-calculated using function
                    precalculate_VtinvSigmaV(V) and used across calls of VB_diarization.
      minDur      - minimum number of frames between speaker turns imposed by linear
                    chains of HMM states corresponding to each speaker. All the states
                    in a chain share the same output distribution
      ref         - T dim. integer vector with reference speaker ID (0:maxSpeakers)
                    per frame
      plot        - if set to True, plot per-frame speaker posteriors.
  
       Outputs:
       gamma  - S x T matrix of posteriors attribution each frame to one of S possible
            speakers, where S is given by opts.maxSpeakers
       pi - S dimensional column vector of ML learned speaker priors. Ideally, these
            should allow to estimate # of speaker in the utterance as the
            probabilities of the redundant speaker should converge to zero.
       Li - values of auxiliary function (and DER and frame cross-entropy between gamma
            and reference if 'ref' is provided) over iterations.
      """
  
      # The references to equations corresponds to
      # Diez. M., Burget. L., Landini. F., Cernocky. J.
      # Analysis of Speaker Diarization based on Bayesian HMM with Eigenvoice Priors
  
      D = X.shape[1]  # feature dimensionality
      C = len(w)  # number of mixture components
      R = V.shape[0]  # subspace rank
      nframes = X.shape[0]
  
      if VtinvSigmaV is None:
          VtinvSigmaV = precalculate_VtinvSigmaV(V, invSigma)
  
      V = V.reshape(V.shape[0], -1)
  
      if pi is None:
          pi = np.ones(maxSpeakers) / maxSpeakers
      else:
          maxSpeakers = len(pi)
  
      if gamma is None:
          # initialize gamma from flat Dirichlet prior with concentration parameter
          # alphaQInit
          gamma = np.random.gamma(
              alphaQInit, size=(nframes, maxSpeakers)
          )
          gamma = gamma / gamma.sum(1, keepdims=True)
  
      # calculate UBM mixture frame posteriors (i.e. per-frame zero order statistics)
      ll = (
          (X**2).dot(-0.5 * invSigma.T)
          + X.dot(invSigma.T * m.T)
          - 0.5
          * (
              (invSigma * m**2 - np.log(invSigma)).sum(1)
              - 2 * np.log(w)
              + D * np.log(2 * np.pi)
          )
      )
      ll *= llScale
      G = logsumexp_ne(ll, axis=1)
      zeta = exp_ne(ll - G[:, np.newaxis])
      zeta[zeta < sparsityThr] = 0.0
      zeta = zeta * statScale
      G = G * statScale
  
      # Kx = np.sum(zeta * (np.log(w) - np.log(zeta)), 1)
      zeta = coo_matrix(
          zeta
      )  # represent zero-order stats using sparse matrix
      print(
          "Sparsity: ",
          len(zeta.row),
          float(len(zeta.row)) / np.prod(zeta.shape),
      )
      LL = np.sum(G)  # total log-likelihod as calculated using UBM
  
      mixture_sum = coo_matrix(
          (np.ones(C * D), (np.repeat(range(C), D), range(C * D)))
      )
      # Aleready calculated above
      # G = np.sum((zeta.multiply(ll - np.log(w))).toarray(), 1) + Kx  # from eq. (30)
  
      # Calculate per-frame first order statistics projected into the R-dim. subspace
      # V^T \Sigma^{-1} F_m
      F_s = coo_matrix(
          (
              (
                  (X[zeta.row] - m[zeta.col]) * zeta.data[:, np.newaxis]
              ).flat,
              (
                  zeta.row.repeat(D),
                  zeta.col.repeat(D) * D
                  + np.tile(range(D), len(zeta.col)),
              ),
          )
      )
      rho = F_s.tocsr().dot((invSigma.flat * V).T)
      del F_s
      # # The code above is only efficient implementation of the following comented code
      # rho = 0;
      # for ii in range(C):
      #  rho = rho + V[ii*D:(ii+1)*D,:].T.dot(zeta[ii,:] * invSigma[:,[ii]] * \
      #       (X - m[:,[ii]]))
  
      if downsample is not None:
          # Downsample zeta, rho, G and gamma by summing the statistic
          # over 'downsample' frames
          # This speeds-up diarization for the price of lowering its frame resolution
          # downsampler = coo_matrix(
          #   (np.ones(nframes),
          #   (np.ceil(np.arange(nframes)/downsample).astype(int), np.arange(nframes)))
          # )
          downsampler = coo_matrix(
              (
                  np.ones(nframes),
                  (
                      np.ceil(np.arange(nframes) / downsample).astype(
                          int
                      ),
                      np.arange(nframes),
                  ),
              )
          )
          zeta = downsampler.dot(zeta)
          rho = downsampler.dot(rho)
          G = downsampler.dot(G)
          gamma = downsampler.dot(gamma) / downsample
      else:
          downsampler = np.array(1)
  
      Li = [[LL * Fa]]  # for the 0-th iteration,
      if ref is not None:
          Li[-1] += [
              DER(downsampler.T.dot(gamma), ref),
              DER(downsampler.T.dot(gamma), ref, xentropy=True),
          ]
  
      ln_p = np.zeros_like(gamma)
      tr = np.eye(minDur * maxSpeakers, k=1)
      ip = np.zeros(minDur * maxSpeakers)
      for ii in range(maxIters):
          ELBO = 0  # objective function (11) (i.e. VB lower-bound on the evidence)
          sum_gamma_zeta = zeta.T.dot(
              gamma
          ).T  # corresponds to the last sum in eq. (26) for all 's'
          invLnoI_flat = sum_gamma_zeta.astype(VtinvSigmaV.dtype).dot(
              VtinvSigmaV
          )  # eq. (26) except for 'I' and the F_A F_B factors for all 's'
          sum_gamma_rho = gamma.T.dot(rho)  # summation in eq. (17)
          for sid in range(maxSpeakers):
              invL = np.linalg.inv(
                  np.eye(R) + tril_to_sym(invLnoI_flat[sid]) * Fa / Fb
              )  # eq. (18) inverse
              a = invL.dot(sum_gamma_rho[sid]) * Fa / Fb  # eq. (17)
              ln_p[:, sid] = Fa * (
                  G
                  + rho.dot(a)
                  - 0.5
                  * zeta.dot(
                      mixture_sum.dot(
                          (
                              (invL + np.outer(a, a))
                              .astype(V.dtype)
                              .dot(V)
                              * (invSigma.flat * V)
                          ).sum(0)
                      )
                  )
              )  # eq. (23)
              ELBO += (
                  Fb
                  * 0.5
                  * (
                      logdet(invL)
                      - np.sum(np.diag(invL) + a**2, 0)
                      + R
                  )
              )
  
          # Construct transition probability matrix with linear chain of 'minDur'
          # states for each of 'maxSpeaker' speaker. The last state in each chain has
          # self-loop probability 'loopProb' and the transition probabilities to the
          # initial chain states given by vector '(1-loopProb) * pi'. From all other,
          # states, one must move to the next state in the chain with probability one.
          tr[minDur - 1 :: minDur, 0::minDur] = (1 - loopProb) * pi
          tr[
              (np.arange(1, maxSpeakers + 1) * minDur - 1,) * 2
          ] += loopProb
          ip[::minDur] = pi
          # per-frame HMM state posteriors. Note that we can have linear
          # chain of minDur states
          # for each speaker.
          gamma, tll, lf, lb = forward_backward(
              ln_p.repeat(minDur, axis=1), tr, ip
          )  # , np.arange(1,maxSpeakers+1)*minDur-1)
  
          # Right after updating q(Z), tll is E{log p(X|,Y,Z)} - KL{q(Z)||p(Z)}.
          # ELBO now contains -KL{q(Y)||p(Y)}. Therefore,
          # ELBO+ttl is correct value for ELBO.
          ELBO += tll
          Li.append([ELBO])
  
          # ML estimate of speaker prior probabilities, eq. (24)
          pi = gamma[0, ::minDur] + np.exp(
              logsumexp(lf[:-1, minDur - 1 :: minDur], axis=1)[
                  :, np.newaxis
              ]
              + lb[1:, ::minDur]
              + ln_p[1:]
              + np.log((1 - loopProb) * pi)
              - tll
          ).sum(0)
          pi = pi / pi.sum()
  
          # per-frame speaker posteriors (eq. (19)), obtained by summing
          # HMM state posteriors corresponding to each speaker
          gamma = gamma.reshape(len(gamma), maxSpeakers, minDur).sum(
              axis=2
          )
  
          # if reference is provided, report DER, cross-entropy and plot the figures
          if ref is not None:
              Li[-1] += [
                  DER(downsampler.T.dot(gamma), ref),
                  DER(downsampler.T.dot(gamma), ref, xentropy=True),
              ]
  
              if plot:
                  import matplotlib.pyplot
  
                  if ii == 0:
                      matplotlib.pyplot.clf()
                  matplotlib.pyplot.subplot(maxIters, 1, ii + 1)
                  matplotlib.pyplot.plot(downsampler.T.dot(gamma), lw=2)
                  matplotlib.pyplot.imshow(
                      np.atleast_2d(ref),
                      interpolation="none",
                      aspect="auto",
                      cmap=matplotlib.pyplot.cm.Pastel1,
                      extent=(0, len(ref), -0.05, 1.05),
                  )
              print(ii, Li[-2])
  
          if ii > 0 and ELBO - Li[-2][0] < epsilon:
              if ELBO - Li[-1][0] < 0:
                  print(
                      "WARNING: Value of auxiliary function has decreased!"
                  )
              break
  
      if downsample is not None:
          # upsample resulting gamma to match number of frames in the input utterance
          gamma = downsampler.T.dot(gamma)
  
      return gamma, pi, Li
  
file_path: espnet/egs/libri_css/asr1/diarization/VB_diarization.py
imports: import numexpr
  import numpy
  import scipy.linalg
  from scipy.sparse import coo_matrix
label_desc: Use image node_id 1 for calling a global function with example usage: VB_diarization(X, m, invSigma, w, V, pi, gamma, maxSpeakers, maxIters, epsilon, loopProb, statScale, alphaQInit, downsample, VtinvSigmaV, ref, plot, sparsityThr, llScale, minDur, Fa, Fb) and returns: gamma,  pi,  Li