code: def _estimate_imbalance(
      info, data, cals, n_imbalance, origin, ext_order
  ):
      """Estimate gradiometer imbalance parameters."""
      mag_scale = 100.0
      n_iterations = 3
      mag_picks = pick_types(info, meg="mag", exclude=())
      grad_picks = pick_types(info, meg="grad", exclude=())
      data = data.copy()
      data[mag_picks, :] *= mag_scale
      del mag_picks
  
      grad_imb = np.zeros((len(grad_picks), n_imbalance))
      exp = dict(origin=origin, int_order=0, ext_order=ext_order)
      all_coils = _prep_mf_coils(info, ignore_ref=True)
      grad_point_coils = _get_grad_point_coilsets(
          info, n_imbalance, ignore_ref=True
      )
      S_orig = _trans_sss_basis(exp, all_coils, coil_scale=mag_scale)
      S_orig /= cals[:, np.newaxis]
      # Compute point gradiometers for each grad channel
      this_cs = np.array([mag_scale], float)
      S_pt = np.array(
          [
              _trans_sss_basis(exp, coils, None, this_cs)
              for coils in grad_point_coils
          ]
      )
      for k in range(n_iterations):
          S_tot = S_orig.copy()
          # In theory we could zero out the homogeneous components with:
          #     S_tot[grad_picks, :3] = 0
          # But in practice it doesn't seem to matter
          S_recon = S_tot[grad_picks]
  
          # Add influence of point magnetometers
          S_tot[grad_picks, :] += np.einsum(
              "ij,ijk->jk", grad_imb.T, S_pt
          )
  
          # Compute multipolar moments
          mm = np.dot(_col_norm_pinv(S_tot.copy())[0], data)
  
          # Use good channels to recalculate
          prev_imb = grad_imb.copy()
          data_recon = np.dot(S_recon, mm)
          assert S_pt.shape == (
              n_imbalance,
              len(grad_picks),
              S_tot.shape[1],
          )
          khi_pts = (S_pt @ mm).transpose(1, 2, 0)
          assert khi_pts.shape == (
              len(grad_picks),
              data.shape[1],
              n_imbalance,
          )
          residual = data[grad_picks] - data_recon
          assert residual.shape == (len(grad_picks), data.shape[1])
          d = (residual[:, np.newaxis, :] @ khi_pts)[:, 0]
          assert d.shape == (len(grad_picks), n_imbalance)
          dinv, _, _ = _reg_pinv(
              khi_pts.swapaxes(-1, -2) @ khi_pts, rcond=1e-6
          )
          assert dinv.shape == (
              len(grad_picks),
              n_imbalance,
              n_imbalance,
          )
          grad_imb[:] = (d[:, np.newaxis] @ dinv)[:, 0]
          # This code is equivalent but hits a np.linalg.pinv bug on old NumPy:
          # grad_imb[:] = np.sum(  # dot product across the time dim
          #     np.linalg.pinv(khi_pts) * residual[:, np.newaxis], axis=-1)
          deltas = np.linalg.norm(grad_imb - prev_imb) / max(
              np.linalg.norm(grad_imb), np.linalg.norm(prev_imb)
          )
          logger.debug(
              f"        Iteration {k + 1}/{n_iterations}: "
              f"max ∆ = {100 * deltas.max():7.3f}%"
          )
      imb = np.zeros((len(data), n_imbalance))
      imb[grad_picks] = grad_imb
      return imb
  
imports: from collections import defaultdict
  from functools import partial
  import numpy
  from scipy.optimize import fmin_cobyla
  from .._fiff.pick import pick_info, pick_types
  from .._fiff.tag import _coil_trans_to_loc, _loc_to_coil_trans
  from ..bem import _check_origin
  from ..io import BaseRaw
  from ..transforms import _find_vector_rotation
  from ..utils import _check_fname, _check_option, _ensure_int, _pl, _reg_pinv, _validate_type, check_fname, logger, verbose
  from .maxwell import _col_norm_pinv, _get_grad_point_coilsets, _prep_fine_cal, _prep_mf_coils, _read_cross_talk, _trans_sss_basis
label_desc: Use image node_id 9 for calling a global function with example usage: _estimate_imbalance(info, data, cals, n_imbalance, origin, ext_order) and returns: imb