Ugly Gaussian symbolic integration, how to extract coefficients and exponent argument?

Question

I have a very ugly function with 157 terms that I want to integrate over three variables. Each term is a Gaussian integral in those variables multiplied by a polynomial in the variables. However, in each term the argument of the exponential is different (still Gaussian) and each polynomial is different. I want to I will give an example term at the bottom (I've also attached a picture in case this is too tough to parse).

Is there a way to extract the coefficient in front of the exponent in such an expression, and to extract the argument of the exponent? I'm asking because I saw a post that deals with improving the speed of these types of integrals. That post provides a module to do the integration given the exponent's argument, the coefficient expression, and the variables. So if I can extract the relevant argument and exponent, I should in principle be able to use that module.

(I E^((3 I m \[Omega] (-2 q q1 + (q^2 + q1^2) Cos[t1 \[Omega]]) Csc[
    t1 \[Omega]])/(2 \[HBar]) + (
  I m \[Omega] (-2 q1 q2 + (q1^2 + 
        q2^2) Cos[(-t1 + t2) \[Omega]]) Csc[(-t1 + t2) \[Omega]])/(
  2 \[HBar]) + (
  I m \[Omega] (-2 q q2 + (q^2 + 
        q2^2) Cos[(-t2 + t3 - I \[Beta]) \[Omega]]) Csc[(-t2 + t3 - 
       I \[Beta]) \[Omega]])/(2 \[HBar])) m q1^5 q2 \[Lambda] Cot[
  t1 \[Omega]] Csc[t1 \[Omega]]^3 Sqrt[-((
  I m \[Omega] Csc[t1 \[Omega]])/\[HBar])] Sqrt[-((
  I m \[Omega] Csc[(-t1 + t2) \[Omega]])/\[HBar])] Sqrt[-((
  I m \[Omega] Csc[(-t2 + t3 - 
       I \[Beta]) \[Omega]])/\[HBar])])/(64 Sqrt[2] \[Pi]^(
 5/2) \[HBar])

Answer 1

I think you're referring to my answer in the linked thread. Using that, you could do the following:

gaussMoment[fPre_, fExp_, vars_] := 
 Module[{coeff, dist, ai, μ, norm}, 
  coeff = CoefficientArrays[fExp, vars, "Symmetric" -> True];
  ai = Inverse[2 coeff[[3]]];
  μ = -ai.coeff[[2]];
  dist = MultinormalDistribution[μ, -ai];
  norm = 1/PDF[dist, vars] /. Thread[vars -> μ];
  Simplify[
   norm Exp[1/2 coeff[[2]].μ + coeff[[1]]] Distribute@
     Expectation[fPre, vars \[Distributed] dist]]]


expr = (I E^((3 I m ω (-2 q q1 + (q^2 + q1^2) Cos[
              t1 ω]) Csc[t1 ω])/(2 ℏ) + (I m ω (-2 q1 q2 + \
(q1^2 + q2^2) Cos[(-t1 + t2) ω]) Csc[(-t1 + 
              t2) ω])/(2 ℏ) + (I m ω (-2 q q2 + \
(q^2 + q2^2) Cos[(-t2 + t3 - I β) ω]) Csc[(-t2 + t3 - 
              I β) ω])/(2 ℏ)) m q1^5 q2 λ \
Cot[t1 ω] Csc[t1 ω]^3 Sqrt[-((I m ω Csc[
           t1 ω])/ℏ)] Sqrt[-((I m ω Csc[(-t1 + 
              t2) ω])/ℏ)] Sqrt[-((I m ω Csc[(-t2 +
               t3 - I β) ω])/ℏ)])/(64 Sqrt[2] Pi^(5/2) ℏ);

expr /. Times[f___, Power[E, g_]] :> gaussMoment[Times[f], g, {q, q2, q1}]

$$\frac{15 \lambda \hbar ^2 \cot (\text{t1} \omega ) \csc ^3(\text{t1} \omega ) \sqrt{-\frac{i m \omega \csc (\text{t1} \omega )}{\hbar }} (\sin (\omega (i \beta +2 \text{t1}-\text{t3}))-2 \sin (\omega (\text{t3}-i \beta )))^2 \sqrt{\frac{i m \omega \csc (\omega (\text{t1}-\text{t2}))}{\hbar }} (-2 \sin (\omega (-i \beta +\text{t1}-\text{t2}+\text{t3}))+\sin (\omega (i \beta +\text{t1}+\text{t2}-\text{t3}))+3 \sin (\omega (\text{t1}-\text{t2}))) \sqrt{\frac{i m \omega \csc (\omega (i \beta +\text{t2}-\text{t3}))}{\hbar }} \sqrt{\frac{i \hbar ^3 \sin (\text{t1} \omega ) \sin (\omega (\text{t2}-\text{t1})) \sin (\omega (-i \beta -\text{t2}+\text{t3}))}{m^3 \omega ^3 (\cos (\omega (i \beta +2 \text{t1}-\text{t3}))-4 \cos (\omega (\text{t3}-i \beta ))+3)}}}{256 \sqrt{2} \pi m^2 \omega ^3 (\cos (\omega (i \beta +2 \text{t1}-\text{t3}))-4 \cos (\omega (\text{t3}-i \beta ))+3)^3}$$

What I did here is to replace the combination of E^(something) times any prefactor by the Gaussian integral, as calculated with the function gaussMoment from the answer linked here. This assumes that the expression is simply a sum of such exponentials, so that the integration over q, q1 and q2 can be distributed over it. That's what the replacement rule essentially does.