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

@@ -421,26 +421,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