code: def _comp_sums_meg(beta, ctheta, lut_fun, n_fact, volume_integral):
      """Lead field dot products using Legendre polynomial (P_n) series.
  
      Parameters
      ----------
      beta : array, shape (n_points * n_points, 1)
          Coefficients of the integration.
      ctheta : array, shape (n_points * n_points, 1)
          Cosine of the angle between the sensor integration points.
      lut_fun : callable
          Look-up table for evaluating Legendre polynomials.
      n_fact : array
          Coefficients in the integration sum.
      volume_integral : bool
          If True, compute volume integral.
  
      Returns
      -------
      sums : array, shape (4, n_points * n_points)
          The results.
      """
      # Compute the sums occurring in the evaluation.
      # Two point magnetometers on the xz plane are assumed.
      # The four sums are:
      #  * sums[:, 0]    n(n+1)/(2n+1) beta^(n+1) P_n
      #  * sums[:, 1]    n/(2n+1) beta^(n+1) P_n'
      #  * sums[:, 2]    n/((2n+1)(n+1)) beta^(n+1) P_n'
      #  * sums[:, 3]    n/((2n+1)(n+1)) beta^(n+1) P_n''
  
      # This is equivalent, but slower:
      # sums = np.sum(bbeta[:, :, np.newaxis].T * n_fact * coeffs, axis=1)
      # sums = np.rollaxis(sums, 2)
      # or
      # sums = np.einsum('ji,jk,ijk->ki', bbeta, n_fact, lut_fun(ctheta)))
      sums = np.empty((n_fact.shape[1], len(beta)))
      # beta can be e.g. 3 million elements, which ends up using lots of memory
      # so we split up the computations into ~50 MB blocks
      n_chunk = 50000000 // (8 * max(n_fact.shape) * 2)
      lims = np.concatenate(
          [np.arange(0, beta.size, n_chunk), [beta.size]]
      )
      for start, stop in zip(lims[:-1], lims[1:]):
          bbeta = np.tile(
              beta[start:stop][np.newaxis], (n_fact.shape[0], 1)
          )
          bbeta[0] *= beta[start:stop]
          np.cumprod(bbeta, axis=0, out=bbeta)  # run inplace
          np.einsum(
              "ji,jk,ijk->ki",
              bbeta,
              n_fact,
              lut_fun(ctheta[start:stop]),
              out=sums[:, start:stop],
          )
      return sums
  
file_path: mne-python/mne/forward/_lead_dots.py
imports: import os
  import os.path
  import numpy
  from numpy.polynomial import legendre
  from ..parallel import parallel_func
  from ..utils import _get_extra_data_path, fill_doc, logger, verbose
label_desc: Use image node_id 6 for calling a global function with example usage: _comp_sums_meg(beta, ctheta, lut_fun, n_fact, volume_integral) and returns: sums