Changes to multi-threading in bayer helms computations

Only running parallel for-loops if the number of fields is
high (>= 100) for now whilst I experiment with different
scheduling and implementing parallelisation in different
sections of the code.

finesse/cmath/cmath.pxd

@@ -55,7 +55,7 @@ cdef class ComplexMath:
     cdef complex_t u_nm(int n, int m, complex_t qx, complex_t qy, double x, double y, double nr, double lambda0) nogil
 
     @staticmethod
-    cdef complex_t u_single_plane(int n, complex_t q, double x, double nr) nogil
+    cdef complex_t u_single_plane(int n, complex_t q, double x, double nr, double lambda0) nogil
 
 cdef class Gaussian:
     @staticmethod
@@ -65,7 +65,7 @@ cdef class Gaussian:
     cdef double w0_size(complex_t q, double nr, double lambda0) nogil
 
     @staticmethod
-    cdef double w_size(complex_t q, double nr) nogil
+    cdef double w_size(complex_t q, double nr, double lambda0) nogil
 
     @staticmethod
     cdef double z_q(complex_t q) nogil

finesse/knm.pyx

@@ -239,6 +239,44 @@ cdef complex_t k_nmnm(
 
     return z1 * z2
 
+cdef void _bayer_helms_outer_loop(
+    complex_t[:, ::1] out,
+    complex_t qx1, complex_t qy1,
+    complex_t qx2, complex_t qy2,
+    double xgamma, double ygamma,
+    double nr_out,
+    int[:, ::1] all_tem_HG,
+    double lambda0,
+    bint reverse_gouy,
+    Py_ssize_t outer_index,
+    Py_ssize_t N
+) nogil:
+    cdef:
+        Py_ssize_t n1, n2, m1, m2
+        Py_ssize_t i = outer_index
+        Py_ssize_t j
+        complex_t knm_single
+
+    n1 = all_tem_HG[i][0]
+    m1 = all_tem_HG[i][1]
+    for j in range(N):
+        n2 = all_tem_HG[j][0]
+        m2 = all_tem_HG[j][1]
+        knm_single = k_nmnm(
+            n1, m1, n2, m2,
+            qx1, qy1, qx2, qy2,
+            xgamma, ygamma,
+            nr_out,
+            lambda0
+        )
+        if reverse_gouy:
+            knm_single = rev_gouy(
+                knm_single,
+                n1, m1, n2, m2,
+                qx1, qy1, qx2, qy2
+            )
+        out[i][j] = knm_single
+
 cpdef void run_bayer_helms(
     complex_t[:, ::1] out,
     complex_t qx1, complex_t qy1,
@@ -289,30 +327,22 @@ cpdef void run_bayer_helms(
     """
     cdef:
         Py_ssize_t N = all_tem_HG.shape[0]
-        int n1, n2, m1, m2
-        Py_ssize_t i, j
-        complex_t knm_single
-
-    for i in prange(N, nogil=True):
-        n1 = all_tem_HG[i][0]
-        m1 = all_tem_HG[i][1]
-        for j in range(N):
-            n2 = all_tem_HG[j][0]
-            m2 = all_tem_HG[j][1]
-            knm_single = k_nmnm(
-                n1, m1, n2, m2,
-                qx1, qy1, qx2, qy2,
-                xgamma, ygamma,
-                nr_out,
-                lambda0
-            )
-            if reverse_gouy:
-                knm_single = rev_gouy(
-                    knm_single,
-                    n1, m1, n2, m2,
-                    qx1, qy1, qx2, qy2
-                )
-            out[i][j] = knm_single
+        Py_ssize_t i
+
+    # NOTE: only do parallel for loop if nfields >= 100 for now
+    if N >= 100:
+        # TODO do some in-depth profiling of different schedules
+        for i in prange(N, nogil=True, schedule="static"):
+            _bayer_helms_outer_loop(out, qx1, qy1, qx2, qy2,
+                             xgamma, ygamma, nr_out,
+                             all_tem_HG, lambda0, reverse_gouy,
+                             i, N)
+    else:
+        for i in range(N):
+            _bayer_helms_outer_loop(out, qx1, qy1, qx2, qy2,
+                             xgamma, ygamma, nr_out,
+                             all_tem_HG, lambda0, reverse_gouy,
+                             i, N)
 
 cpdef void zero_tem00_phase(complex_t[:, ::1] knm_mat, complex_t[:, ::1] out=None):
     """
@@ -422,7 +452,7 @@ cpdef void rev_all_gouy(
         complex_t gouy_rev_knm
         Py_ssize_t N = knm_mat.shape[0]
 
-    for i in prange(N, nogil=True):
+    for i in range(N):
         n1 = all_tem_HG[i][0]
         m1 = all_tem_HG[i][1]
         for j in range(N):

playground/profiling/mismatching_cavity.py

+import cProfile
+import pstats
+
+import finesse
+
+model = finesse.parse("""
+l L0 P=1
+
+s s0 L0.p1 ITM.p1 L=1
+
+m ITM R=0.9 T=0.1 Rc=-2.5
+s CAV ITM.p2 ETM.p1 L=1
+m ETM R=0.99 T=0.01 Rc=2.5
+
+attr ITM Rc -2.5
+attr ETM Rc 2.5
+
+maxtem 6
+
+cav FP ITM.p2.o ITM.p2.i
+
+xaxis &ETM.Rc 2.5 10 200
+""")
+ifo = model.model
+ifo.create_mismatch(ifo.ITM.p1.i, w0_mm=15, z_mm=7)
+ifo.add_all_ad(ifo.ETM.p2.o, 0)
+
+cProfile.run("model.run()", "mismatched_cav_scan_profile.prof")
+p = pstats.Stats("mismatched_cav_scan_profile.prof")
+p.sort_stats("time")
+p.print_stats(25)
+