code: def eager_integrate_gaussian_gaussian(
      log_measure, integrand, reduced_vars
  ):
      assert log_measure.is_full_rank
      reduced_names = frozenset(v.name for v in reduced_vars)
      real_vars = frozenset(
          v.name for v in reduced_vars if v.dtype == "real"
      )
      if real_vars:
          lhs_reals = frozenset(
              k
              for k, d in log_measure.inputs.items()
              if d.dtype == "real"
          )
          rhs_reals = frozenset(
              k
              for k, d in integrand.inputs.items()
              if d.dtype == "real"
          )
          if lhs_reals == real_vars and rhs_reals <= real_vars:
              inputs = OrderedDict(
                  (k, d)
                  for t in (log_measure, integrand)
                  for k, d in t.inputs.items()
              )
              lhs_white_vec, lhs_prec_sqrt = align_gaussian(
                  inputs, log_measure
              )
              rhs_white_vec, rhs_prec_sqrt = align_gaussian(
                  inputs, integrand
              )
              lhs = Gaussian(
                  white_vec=lhs_white_vec,
                  prec_sqrt=lhs_prec_sqrt,
                  inputs=inputs,
              )
  
              # Compute the expectation of a non-normalized quadratic form.
              # See "The Matrix Cookbook" (November 15, 2012) ss. 8.2.2 eq. 380.
              # http://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf
              # If x ~ N(mean,cov) then
              #   E[(x-m)' A (x-m)] = (m-mean)'A(m-mean) + Tr(A cov)     # eq. 380
              # To perform this computation in rhs's internal space, we first transform
              # lhs to rhs's whitened space
              mean = _vm(lhs._mean, rhs_prec_sqrt)
              norm = ops.exp(lhs._log_normalizer)
              # Then in rhs's whitened space, A = I so Tr(A cov) = Tr(cov).
              vmv_term = _norm2(rhs_white_vec - mean)
              trace_term = (
                  (
                      ops.triangular_solve(
                          rhs_prec_sqrt, lhs._precision_chol
                      )
                      ** 2
                  )
                  .sum(-1)
                  .sum(-1)
              )
              data = (-0.5) * norm * (vmv_term + trace_term)
  
              inputs = OrderedDict(
                  (k, d)
                  for k, d in inputs.items()
                  if k not in reduced_names
              )
              result = Tensor(data, inputs)
              return result.reduce(ops.add, reduced_names - real_vars)
  
          raise NotImplementedError(
              "TODO implement partial integration"
          )
  
      return None
  
imports: from collections import OrderedDict
  from functools import reduce
  from typing import Tuple, Union
  import funsor.ops
  from funsor.cnf import Contraction, GaussianMixture
  from funsor.constant import Constant
  from funsor.delta import Delta
  from funsor.gaussian import Gaussian, _norm2, _vm, align_gaussian
  from funsor.interpretations import eager, normalize
  from funsor.tensor import Tensor
  from funsor.terms import Funsor, FunsorMeta, Number, Subs, Unary, Variable, _convert_reduced_vars, substitute, to_funsor
label_desc: Use image node_id 7 for calling a global function with example usage: eager_integrate_gaussian_gaussian(log_measure, integrand, reduced_vars) and returns: None, result