Update quantum code from matgwinc, and vectorize

Before this change and !31 (merged):

%timeit gwinc.gwinc(freq, ifo)

519 ms ± 15.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

After both:

%timeit gwinc.gwinc(freq, ifo)

50.4 ms ± 1.97 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%prun gwinc.gwinc(freq, ifo)

     5304 function calls (4870 primitive calls) in 0.102 seconds

Ordered by: internal time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    74    0.033    0.000    0.034    0.000 coatingthermal.py:743(getCoatRefl)
     1    0.014    0.014    0.016    0.016 quantum.py:257(shotradSignalRecycled)
    10    0.008    0.001    0.008    0.001 suspension.py:30(tst_force_to_tst_displ)
     2    0.007    0.004    0.007    0.004 coatingthermal.py:446(getCoatTOPhase)
     1    0.006    0.006    0.016    0.016 suspension.py:51(suspQuad)
     7    0.006    0.001    0.006    0.001 {built-in method numpy.core.multiarray.matmul}

I think the biggest remaining chunk of run time goes to the coating design step in precompIFO (which can be skipped by setting CoatLayerOpticalThickness in the ifo struct).

gwinc/noise/quantum.py

@@ -322,34 +322,41 @@ def shotradSignalRecycled(f, ifo):
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     # Coefficients [BnC, Equations 5.8 to 5.12]
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+    exp_1jbeta = exp(1j*beta)
+    cos_beta = exp_1jbeta.real
+    invexp_1jbeta = 1/exp_1jbeta
+    exp_2jbeta = exp_1jbeta**2
+    cos_2beta = exp_2jbeta.real
+    invexp_2jbeta = 1/exp_2jbeta
+    exp_4jbeta = exp_2jbeta**2
     C11_L   = ( (1+rho**2) * ( cos(2*phi) + Kappa/2 * sin(2*phi) ) -
-                2*rho*cos(2*beta) - 1/4*epsilon * ( -2 * (1+exp(2j*beta))**2 * rho + 4 * (1+rho**2) *
-                                                    cos(beta)**2*cos(2*phi) + ( 3+exp(1j*2*beta) ) *
+                2*rho*cos_2beta - 1/4*epsilon * ( -2 * (1+exp_2jbeta)**2 * rho + 4 * (1+rho**2) *
+                                                    cos_beta**2*cos(2*phi) + ( 3+exp_2jbeta ) *
                                                     Kappa * (1+rho**2) * sin(2*phi) ) +
-                lambda_SR * ( exp(2j*beta)*rho-1/2 * (1+rho**2) * ( cos(2*phi)+Kappa/2 * sin(2*phi) ) ) )
+                lambda_SR * ( exp_2jbeta*rho-1/2 * (1+rho**2) * ( cos(2*phi)+Kappa/2 * sin(2*phi) ) ) )
 
     C22_L   = C11_L
 
     C12_L   = tau**2 * ( - ( sin(2*phi) + Kappa*sin(phi)**2 )+
-                         1/2*epsilon*sin(phi) * ( (3+exp(2j*beta)) * Kappa * sin(phi) + 4*cos(beta)**2 * cos(phi)) +
+                         1/2*epsilon*sin(phi) * ( (3+exp_2jbeta) * Kappa * sin(phi) + 4*cos_beta**2 * cos(phi)) +
                          1/2*lambda_SR * ( sin(2*phi)+Kappa*sin(phi)**2) )
 
     C21_L   = tau**2 * ( (sin(2*phi)-Kappa*cos(phi)**2 ) +
-                         1/2*epsilon*cos(phi) * ( (3+exp(2j*beta) )*Kappa*sin(phi) - 4*cos(beta)**2*sin(phi) ) +
+                         1/2*epsilon*cos(phi) * ( (3+exp_2jbeta )*Kappa*sin(phi) - 4*cos_beta**2*sin(phi) ) +
                          1/2*lambda_SR * ( -sin(2*phi) + Kappa*cos(phi)**2) )
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-    D1_L    = ( - (1+rho*exp(2j*beta) ) * sin(phi) +
-                1/4*epsilon * ( 3+rho+2*rho*exp(4*1j*beta) + exp(2j*beta)*(1+5*rho) ) * sin(phi)+
-                1/2*lambda_SR * exp(2j*beta) * rho * sin(phi) )
+    D1_L    = ( - (1+rho*exp_2jbeta ) * sin(phi) +
+                1/4*epsilon * ( 3+rho+2*rho*exp_4jbeta + exp_2jbeta*(1+5*rho) ) * sin(phi)+
+                1/2*lambda_SR * exp_2jbeta * rho * sin(phi) )
 
-    D2_L    = ( - (-1+rho*exp(2j*beta) ) * cos(phi) +
-                1/4*epsilon * ( -3+rho+2*rho*exp(4*1j*beta) + exp(2j*beta) * (-1+5*rho) ) * cos(phi)+
-                1/2*lambda_SR * exp(2j*beta) * rho * cos(phi) )
+    D2_L    = ( - (-1+rho*exp_2jbeta ) * cos(phi) +
+                1/4*epsilon * ( -3+rho+2*rho*exp_4jbeta + exp_2jbeta * (-1+5*rho) ) * cos(phi)+
+                1/2*lambda_SR * exp_2jbeta * rho * cos(phi) )
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     P11     = 0.5 * sqrt(lambda_SR) * tau * \
-              ( -2*rho*exp(2j*beta)+2*cos(2*phi)+Kappa*sin(2*phi) )
+              ( -2*rho*exp_2jbeta+2*cos(2*phi)+Kappa*sin(2*phi) )
     P22     = P11
     P12     = -sqrt(lambda_SR)*tau*sin(phi)*(2*cos(phi)+Kappa*sin(phi) )
     P21     =  sqrt(lambda_SR)*tau*cos(phi)*(2*sin(phi)-Kappa*cos(phi) )
@@ -360,22 +367,22 @@ def shotradSignalRecycled(f, ifo):
     #   as well as the input-output relation Mc and the signal matrix Md
 
     Q11     = 1 / \
-              ( exp(-2j*beta)+rho**2*exp(2j*beta)-rho*(2*cos(2*phi)+Kappa*sin(2*phi)) +
-                1/2*epsilon*rho * (exp(-2j*beta)*cos(2*phi)+exp(2j*beta)*
-                                   ( -2*rho-2*rho*cos(2*beta)+cos(2*phi)+Kappa*sin(2*phi) ) +
+              ( invexp_2jbeta+rho**2*exp_2jbeta-rho*(2*cos(2*phi)+Kappa*sin(2*phi)) +
+                1/2*epsilon*rho * (invexp_2jbeta*cos(2*phi)+exp_2jbeta*
+                                   ( -2*rho-2*rho*cos_2beta+cos(2*phi)+Kappa*sin(2*phi) ) +
                                    2*cos(2*phi)+3*Kappa*sin(2*phi))-1/2*lambda_SR*rho *
-                ( 2*rho*exp(2j*beta)-2*cos(2*phi)-Kappa*sin(2*phi) ) )
+                ( 2*rho*exp_2jbeta-2*cos(2*phi)-Kappa*sin(2*phi) ) )
     Q22     = Q11
     Q12     = 0
     Q21     = 0
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-    N11     = sqrt(epsilon/2)*tau *(Kappa*(1+rho*exp(2j*beta))*sin(phi)+
-                                    2*cos(beta)*(exp(-1j*beta)*cos(phi)-rho*exp(1j*beta)*(cos(phi)+Kappa*sin(phi))))
-    N22     = -sqrt(2*epsilon)*tau*(-exp(-1j*beta)+rho*exp(1j*beta))*cos(beta)*cos(phi)
-    N12     = -sqrt(2*epsilon)*tau*(exp(-1j*beta)+rho*exp(1j*beta))*cos(beta)*sin(phi);
+    N11     = sqrt(epsilon/2)*tau *(Kappa*(1+rho*exp_2jbeta)*sin(phi)+
+                                    2*cos_beta*(invexp_1jbeta*cos(phi)-rho*exp_1jbeta*(cos(phi)+Kappa*sin(phi))))
+    N22     = -sqrt(2*epsilon)*tau*(-invexp_1jbeta+rho*exp_1jbeta)*cos_beta*cos(phi)
+    N12     = -sqrt(2*epsilon)*tau*(invexp_1jbeta+rho*exp_1jbeta)*cos_beta*sin(phi);
     N21     = sqrt(2*epsilon)*tau*(-Kappa*(1+rho)*cos(phi)+
-                                   2*cos(beta)*(exp(-1j*beta)+rho*exp(1j*beta))*cos(beta)*sin(phi))
+                                   2*cos_beta*(invexp_1jbeta+rho*exp_1jbeta)*cos_beta*sin(phi))
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     # overall coefficient
@@ -421,26 +428,23 @@ def getProdTF(lhs, rhs):
         raise Exception('Matrix size mismatch size(lhs, 2) = %d != %d = size(rhs, 1)' % (lhs.shape[1], rhs.shape[0]))
     N = lhs.shape[0]
     M = rhs.shape[1]
+    if len(lhs.shape) == 3:
+        lhs = np.transpose(lhs, axes=(2, 0, 1))
+    if len(rhs.shape) == 3:
+        rhs = np.transpose(rhs, axes=(2, 0, 1))
 
     # compute product
-    if len(lhs.shape) < 3 or lhs.shape[2] == 1:
-        Nfreq = rhs.shape[2]
-        rslt = zeros((N, M, Nfreq), dtype=complex)
-        for n in range(Nfreq):
-            rslt[:, :, n] = np.dot(np.squeeze(lhs), rhs[:, :, n])
-    elif len(rhs.shape) < 3 or rhs.shape[2] == 1:
-        Nfreq = lhs.shape[2]
-        rslt = zeros((N, M, Nfreq), dtype=complex)
-        for n in range(Nfreq):
-            rslt[:, :, n] = np.dot(lhs[:, :, n], np.squeeze(rhs))
-    elif lhs.shape[2] == rhs.shape[2]:
-        Nfreq = lhs.shape[2]
-        rslt = zeros((N, M, Nfreq), dtype=complex)
-        for n in range(Nfreq):
-            rslt[:, :, n] = np.dot(lhs[:, :, n], rhs[:, :, n])
+    if len(lhs.shape) < 3 or lhs.shape[0] == 1:
+        rslt = np.matmul(lhs, rhs)
+    elif len(rhs.shape) < 3 or rhs.shape[0] == 1:
+        rslt = np.matmul(lhs, rhs)
+    elif lhs.shape[0] == rhs.shape[0]:
+        rslt = np.matmul(lhs, rhs)
     else:
         raise Exception('Matrix size mismatch lhs.shape[2] = %d != %d = rhs.shape[2]' % (lhs.shape[2], rhs.shape[2]))
 
+    if len(rslt.shape) == 3:
+        rslt = np.transpose(rslt, axes=(1, 2, 0))
     return rslt