code: def _pc_relate_bm(
      call_expr: CallExpression,
      min_individual_maf: float,
      *,
      k: Optional[int] = None,
      scores_expr: Optional[ArrayNumericExpression] = None,
      min_kinship: Optional[float] = None,
      statistics: str = "all",
      block_size: Optional[int] = None,
      include_self_kinship: bool = False,
  ) -> Table:
      assert 0.0 <= min_individual_maf <= 1.0, (
          f"invalid argument: min_individual_maf={min_individual_maf}. "
          f"Must have min_individual_maf on interval [0.0, 1.0]."
      )
      mt = matrix_table_source("pc_relate_bm/call_expr", call_expr)
      if k and scores_expr is None:
          eigens, scores, _ = _hwe_normalized_blanczos(
              call_expr, k, compute_loadings=False, q_iterations=10
          )
          scores_table = (
              scores.select(__scores=scores.scores)
              .key_by()
              .select("__scores")
          )
          compute_S0 = False
      elif not k and scores_expr is not None:
          analyze(
              "pc_relate_bm/scores_expr", scores_expr, mt._col_indices
          )
          eigens = None
          scores_table = (
              mt.select_cols(__scores=scores_expr)
              .key_cols_by()
              .select_cols("__scores")
              .cols()
          )
          compute_S0 = True
      elif k and scores_expr is not None:
          raise ValueError(
              "pc_relate_bm: exactly one of 'k' and 'scores_expr' "
              "must be set, found both"
          )
      else:
          raise ValueError(
              "pc_relate_bm: exactly one of 'k' and 'scores_expr' "
              "must be set, found neither"
          )
  
      n_missing = scores_table.aggregate(
          agg.count_where(hl.is_missing(scores_table.__scores))
      )
      if n_missing > 0:
          raise ValueError(
              f"Found {n_missing} columns with missing scores array."
          )
      pc_scores = hl.nd.array(
          scores_table.collect(_localize=False).map(
              lambda x: x.__scores
          )
      )
  
      # Define NaN for missing values, otherwise cannot convert expr to block matrix
      nan = hl.float64(float("NaN"))
  
      # Create genotype matrix, set missing GT entries to NaN
      mt = mt.select_entries(
          __gt=call_expr.n_alt_alleles()
      ).unfilter_entries()
      gt_with_nan_expr = hl.or_else(hl.float64(mt.__gt), nan)
      if not block_size:
          block_size = BlockMatrix.default_block_size()
      g = BlockMatrix.from_entry_expr(
          gt_with_nan_expr, block_size=block_size
      )
      g = g.checkpoint(new_temp_file("pc_relate_bm/g", "bm"))
      sqrt_n_samples = hl.nd.array([hl.sqrt(g.shape[1])])
  
      # Recover singular values, S0, as vector of column norms of pc_scores if necessary
      if compute_S0:
          S0 = (
              (pc_scores ** hl.int32(2))
              .sum(0)
              .map(lambda x: hl.sqrt(x))
          )
      else:
          S0 = hl.nd.array(eigens).map(lambda x: hl.sqrt(x))
      # Set first entry of S to sqrt(n), for intercept term in beta
      S = hl.nd.hstack((sqrt_n_samples, S0))._persist()
      # Recover V from pc_scores with inv(S0)
      V0 = (pc_scores * (1 / S0))._persist()
      # Set all entries in first column of V to 1/sqrt(n), for intercept term in beta
      ones_normalized = hl.nd.full((V0.shape[0], 1), (1 / S[0]))
      V = hl.nd.hstack((ones_normalized, V0))
  
      # Compute matrix of regression coefficients for PCs (beta), shape (k, m)
      beta = (
          BlockMatrix.from_ndarray(
              ((1 / S) * V).T, block_size=block_size
          )
          @ g.T
      )
      beta = beta.checkpoint(new_temp_file("pc_relate_bm/beta", "bm"))
  
      # Compute matrix of individual-specific AF estimates (mu), shape (m, n)
      mu = (
          0.5
          * (
              BlockMatrix.from_ndarray(V * S, block_size=block_size)
              @ beta
          ).T
      )
      # Replace entries in mu with NaN if invalid or if corresponding GT is missing (no contribution from that variant)
      mu = mu._apply_map2(
          lambda _mu, _g: hl.if_else(
              _bad_mu(_mu, min_individual_maf) | hl.is_nan(_g), nan, _mu
          ),
          g,
          sparsity_strategy="NeedsDense",
      )
      mu = mu.checkpoint(new_temp_file("pc_relate_bm/mu", "bm"))
  
      # Compute kinship matrix (phi), shape (n, n)
      # Where mu is NaN (missing), set variance and centered AF to 0 (no contribution from that variant)
      variance = _replace_nan(mu * (1.0 - mu), 0.0).checkpoint(
          new_temp_file("pc_relate_bm/variance", "bm")
      )
      centered_af = _replace_nan(g - (2.0 * mu), 0.0)
      phi = _gram(centered_af) / (4.0 * _gram(variance.sqrt()))
      phi = phi.checkpoint(new_temp_file("pc_relate_bm/phi", "bm"))
      ht = phi.entries().rename({"entry": "kin"})
      ht = ht.annotate(
          k0=hl.missing(hl.tfloat64),
          k1=hl.missing(hl.tfloat64),
          k2=hl.missing(hl.tfloat64),
      )
  
      if statistics in ["kin2", "kin20", "all"]:
          # Compute inbreeding coefficient and dominance encoding of GT matrix
          f_i = (2.0 * phi.diagonal()) - 1.0
          gd = g._apply_map2(
              lambda _g, _mu: _dominance_encoding(_g, _mu),
              mu,
              sparsity_strategy="NeedsDense",
          )
          normalized_gd = gd - (variance * (1.0 + f_i))
  
          # Compute IBD2 (k2) estimate
          k2 = _gram(normalized_gd) / _gram(variance)
          ht = ht.annotate(k2=k2.entries()[ht.i, ht.j].entry)
  
          if statistics in ["kin20", "all"]:
              # Get the numerator used in IBD0 (k0) computation (IBS0), compute indicator matrices for homozygotes
              hom_alt = g._apply_map2(
                  lambda _g, _mu: hl.if_else(
                      (_g != 2.0) | hl.is_nan(_mu), 0.0, 1.0
                  ),
                  mu,
                  sparsity_strategy="NeedsDense",
              )
              hom_ref = g._apply_map2(
                  lambda _g, _mu: hl.if_else(
                      (_g != 0.0) | hl.is_nan(_mu), 0.0, 1.0
                  ),
                  mu,
                  sparsity_strategy="NeedsDense",
              )
              ibs0 = _AtB_plus_BtA(hom_alt, hom_ref)
  
              # Get the denominator used in IBD0 (k0) computation
              mu2 = _replace_nan(mu**2.0, 0.0)
              one_minus_mu2 = _replace_nan((1.0 - mu) ** 2.0, 0.0)
              k0_denom = _AtB_plus_BtA(mu2, one_minus_mu2)
  
              # Compute IBD0 (k0) estimates, correct the estimates where phi <= k0_cutoff
              k0 = ibs0 / k0_denom
              k0_cutoff = 2.0 ** (-5.0 / 2.0)
              ht = ht.annotate(k0=k0.entries()[ht.i, ht.j].entry)
              ht = ht.annotate(
                  k0=hl.if_else(
                      ht.kin <= k0_cutoff,
                      1.0 - (4.0 * ht.kin) + ht.k2,
                      ht.k0,
                  )
              )
  
              if statistics == "all":
                  # Finally, compute IBD1 (k1) estimate
                  ht = ht.annotate(k1=1.0 - (ht.k2 + ht.k0))
  
      # Filter table to only have one row for each distinct pair of samples
      ht = ht.filter(ht.i <= ht.j)
      ht = ht.rename({"k0": "ibd0", "k1": "ibd1", "k2": "ibd2"})
  
      if min_kinship is not None:
          ht = ht.filter(ht.kin >= min_kinship)
  
      if statistics != "all":
          fields_to_drop = {
              "kin": ["ibd0", "ibd1", "ibd2"],
              "kin2": ["ibd0", "ibd1"],
              "kin20": ["ibd1"],
          }
          ht = ht.drop(*fields_to_drop[statistics])
  
      if not include_self_kinship:
          ht = ht.filter(ht.i == ht.j, keep=False)
  
      col_keys = hl.literal(
          mt.select_cols().key_cols_by().cols().collect(),
          dtype=hl.tarray(mt.col_key.dtype),
      )
      return ht.key_by(
          i=col_keys[hl.int32(ht.i)], j=col_keys[hl.int32(ht.j)]
      )
  
imports: from typing import Optional
  import hail
  import hail.expr.aggregators
  from hail import ir
  from hail.backend.spark_backend import SparkBackend
  from hail.expr import ArrayNumericExpression, BooleanExpression, CallExpression, Float64Expression, analyze, expr_array, expr_call, expr_float64, matrix_table_source
  from hail.expr.types import tarray
  from hail.linalg import BlockMatrix
  from hail.table import Table
  from hail.typecheck import enumeration, nullable, numeric, typecheck
  from hail.utils import new_temp_file
  from hail.utils.java import Env
  from ..pca import _hwe_normalized_blanczos, hwe_normalized_pca
label_desc: Use image node_id 7 for calling a global function with example usage: _pc_relate_bm(call_expr, min_individual_maf) and returns: ht