code: def _beta_divergence(X, W, H, beta, square_root=False):
      """Compute the beta-divergence of X and dot(W, H).
  
      Parameters
      ----------
      X : float or array-like of shape (n_samples, n_features)
  
      W : float or array-like of shape (n_samples, n_components)
  
      H : float or array-like of shape (n_components, n_features)
  
      beta : float or {'frobenius', 'kullback-leibler', 'itakura-saito'}
          Parameter of the beta-divergence.
          If beta == 2, this is half the Frobenius *squared* norm.
          If beta == 1, this is the generalized Kullback-Leibler divergence.
          If beta == 0, this is the Itakura-Saito divergence.
          Else, this is the general beta-divergence.
  
      square_root : bool, default=False
          If True, return np.sqrt(2 * res)
          For beta == 2, it corresponds to the Frobenius norm.
  
      Returns
      -------
          res : float
              Beta divergence of X and np.dot(X, H).
      """
      beta = _beta_loss_to_float(beta)
  
      # The method can be called with scalars
      if not sp.issparse(X):
          X = np.atleast_2d(X)
      W = np.atleast_2d(W)
      H = np.atleast_2d(H)
  
      # Frobenius norm
      if beta == 2:
          # Avoid the creation of the dense np.dot(W, H) if X is sparse.
          if sp.issparse(X):
              norm_X = np.dot(X.data, X.data)
              norm_WH = trace_dot(np.linalg.multi_dot([W.T, W, H]), H)
              cross_prod = trace_dot((X * H.T), W)
              res = (norm_X + norm_WH - 2.0 * cross_prod) / 2.0
          else:
              res = squared_norm(X - np.dot(W, H)) / 2.0
  
          if square_root:
              return np.sqrt(res * 2)
          else:
              return res
  
      if sp.issparse(X):
          # compute np.dot(W, H) only where X is nonzero
          WH_data = _special_sparse_dot(W, H, X).data
          X_data = X.data
      else:
          WH = np.dot(W, H)
          WH_data = WH.ravel()
          X_data = X.ravel()
  
      # do not affect the zeros: here 0 ** (-1) = 0 and not infinity
      indices = X_data > EPSILON
      WH_data = WH_data[indices]
      X_data = X_data[indices]
  
      # used to avoid division by zero
      WH_data[WH_data < EPSILON] = EPSILON
  
      # generalized Kullback-Leibler divergence
      if beta == 1:
          # fast and memory efficient computation of np.sum(np.dot(W, H))
          sum_WH = np.dot(np.sum(W, axis=0), np.sum(H, axis=1))
          # computes np.sum(X * log(X / WH)) only where X is nonzero
          div = X_data / WH_data
          res = np.dot(X_data, np.log(div))
          # add full np.sum(np.dot(W, H)) - np.sum(X)
          res += sum_WH - X_data.sum()
  
      # Itakura-Saito divergence
      elif beta == 0:
          div = X_data / WH_data
          res = np.sum(div) - np.prod(X.shape) - np.sum(np.log(div))
  
      # beta-divergence, beta not in (0, 1, 2)
      else:
          if sp.issparse(X):
              # slow loop, but memory efficient computation of :
              # np.sum(np.dot(W, H) ** beta)
              sum_WH_beta = 0
              for i in range(X.shape[1]):
                  sum_WH_beta += np.sum(np.dot(W, H[:, i]) ** beta)
  
          else:
              sum_WH_beta = np.sum(WH**beta)
  
          sum_X_WH = np.dot(X_data, WH_data ** (beta - 1))
          res = (X_data**beta).sum() - beta * sum_X_WH
          res += sum_WH_beta * (beta - 1)
          res /= beta * (beta - 1)
  
      if square_root:
          res = max(
              res, 0
          )  # avoid negative number due to rounding errors
          return np.sqrt(2 * res)
      else:
          return res
  
imports: import itertools
  import time
  import warnings
  from abc import ABC
  from math import sqrt
  from numbers import Integral, Real
  import numpy
  import scipy.sparse
  from scipy import linalg
  from .._config import config_context
  from ..base import BaseEstimator, ClassNamePrefixFeaturesOutMixin, TransformerMixin, _fit_context
  from ..exceptions import ConvergenceWarning
  from ..utils import check_array, check_random_state, gen_batches, metadata_routing
  from ..utils._param_validation import Interval, StrOptions, validate_params
  from ..utils.extmath import randomized_svd, safe_sparse_dot, squared_norm
  from ..utils.validation import check_is_fitted, check_non_negative
  from ._cdnmf_fast import _update_cdnmf_fast
label_desc: Use image node_id 4 for calling a global function with example usage: _beta_divergence(X, W, H, beta, square_root) and returns: np, res, np, res