Issue# 1762 - Add an option to build SWIG python bindings with threads enabled

.github/workflows/build-ci.yml

@@ -151,12 +151,12 @@ jobs:
       if: ${{ matrix.enable-mpi == false && matrix.python-version == 3.6 }}
       run: ./configure --enable-maintainer-mode --with-coverage --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF}
 
-    - name: Run configure with single-precision floating point
+    - name: Run configure with single-precision floating point and swig threads
       if: ${{ matrix.enable-mpi == true && matrix.python-version == 3.9 }}
       run: |
         mkdir -p build &&
         pushd build &&
-        ../configure --enable-maintainer-mode --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF} --enable-single &&
+        ../configure --enable-maintainer-mode --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF} --enable-single --enable-swig-python-threads &&
         popd
 
     - name: Run make

configure.ac

@@ -31,6 +31,33 @@ else
 fi
 AC_SUBST(MEEP_SINGLE)
 
+##############################################################################
+# build with SWIG threads for Python
+
+AC_ARG_ENABLE(swig-python-threads,
+              [AS_HELP_STRING([--enable-swig-python-threads],[enable SWIG threads for Python bindings])],
+	            enable_swig_python_threads=$enableval, enable_swig_python_threads=no)
+if test "x$enable_swig_python_threads" = "xyes"; then
+    MEEP_SWIG_PYTHON_THREADS=1
+else
+    MEEP_SWIG_PYTHON_THREADS=0
+fi
+AC_SUBST(MEEP_SWIG_PYTHON_THREADS)
+AM_CONDITIONAL(MEEP_SWIG_PYTHON_THREADS, test "x$enable_swig_python_threads" = "xyes")
+
+##############################################################################
+# build with SWIG debugging enabled for Python
+
+AC_ARG_ENABLE(swig-python-debug,
+              [AS_HELP_STRING([--enable-swig-python-debug],[enable SWIG debug for Python bindings])],
+	            enable_swig_python_debug=$enableval, enable_swig_python_debug=no)
+if test "x$enable_swig_python_debug" = "xyes"; then
+    MEEP_SWIG_PYTHON_DEBUG=1
+else
+    MEEP_SWIG_PYTHON_DEBUG=0
+fi
+AC_SUBST(MEEP_SWIG_PYTHON_DEBUG)
+
 ##############################################################################
 
 # Check for mpiCC immediately after getting C++ compiler...

doc/docs/Python_Developer_Information.md

@@ -91,7 +91,20 @@ Classes and functions related to the high-level Python interface to `MPB`. Addit
 
 Definition of `MPBData`, a Python class useful for `MPB` data analysis (documented [here](https://mpb.readthedocs.io/en/latest/Python_Data_Analysis_Tutorial)). This is is a Python port of the functionality available in the [`mpb-data` command line program](https://github.com/NanoComp/mpb/blob/master/utils/mpb-data.c) originally written in C.
 
-## Development and Testing
+## Development
+
+By default, the SWIG Python bindings are built with `threads` disabled (GIL is
+held for all SWIG wrapped python calls by default). You can optionally build the
+Python bindings with `threads` enabled (releasing the GIL for all SWIG wrapped
+Python calls) by passing the `--enable-swig-python-threads`
+option to the configure script.
+
+Since the bindings could be built with `threads` enabled, one needs to be
+careful to protect (acquire the GIL) code that calls back into Python or custom
+python wrapper code that uses PyAPI. Look for `SWIG_PYTHON_THREAD_SCOPED_BLOCK`
+in the SWIG interface files or the custom wrapper code for how this is done.
+
+## Testing
 
 The tests for the Python interface are located in `python/tests`. To run the whole test suite, run `make check` in the `python` build tree. During development it is more convenient to run individual tests. This can be accomplished by running `python <path_to_test>/test.py MyTestCase.test_method`. See the [Python unittest framework documentation](https://docs.python.org/3/library/unittest.html) for more info.
 

python/Makefile.am

@@ -164,8 +164,12 @@ HPPFILES=                                    \
  $(top_srcdir)/src/adjust_verbosity.hpp      \
  meep-python.hpp
 
+if MEEP_SWIG_PYTHON_THREADS
+SWIG_MEEP_FLAGS = -threads
+endif # MEEP_SWIG_PYTHON_THREADS
+
 meep-python.cxx: $(MEEP_SWIG_SRC) $(HPPFILES)
-	$(SWIG) -Wextra $(AM_CPPFLAGS) -outdir $(builddir) -c++ -nofastunpack -python -o $@ $(srcdir)/meep.i
+	$(SWIG) -Wextra $(SWIG_MEEP_FLAGS) $(AM_CPPFLAGS) -outdir $(builddir) -c++ -nofastunpack -python -o $@ $(srcdir)/meep.i
 
 if WITH_MPB
 MPB_SWIG_SRC = mpb.i

python/meep-python.hpp

@@ -1,14 +1,22 @@
 namespace meep {
 
+#ifndef SWIG_PYTHON_THREAD_SCOPED_BLOCK
+#define SWIG_PYTHON_THREAD_SCOPED_BLOCK   SWIG_PYTHON_THREAD_BEGIN_BLOCK
+#endif
+
 // like custom_src_time, but using Python function object, with proper reference counting
 class custom_py_src_time : public src_time {
 public:
   custom_py_src_time(PyObject *fun, double st = -infinity, double et = infinity,
                      std::complex<double> f = 0)
       : func(fun), freq(f), start_time(float(st)), end_time(float(et)) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     Py_INCREF(func);
   }
-  virtual ~custom_py_src_time() { Py_DECREF(func); }
+  virtual ~custom_py_src_time() {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
+    Py_DECREF(func);
+  }
 
   virtual std::complex<double> current(double time, double dt) const {
     if (is_integrated)
@@ -17,6 +25,7 @@ class custom_py_src_time : public src_time {
       return dipole(time);
   }
   virtual std::complex<double> dipole(double time) const {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     float rtime = float(time);
     if (rtime >= start_time && rtime <= end_time) {
       PyObject *py_t = PyFloat_FromDouble(time);
@@ -33,6 +42,7 @@ class custom_py_src_time : public src_time {
   }
   virtual double last_time() const { return end_time; };
   virtual src_time *clone() const {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     Py_INCREF(func); // default copy constructor doesn't incref
     return new custom_py_src_time(*this);
   }

python/meep.i

@@ -23,6 +23,23 @@
 #define SWIG_FILE_WITH_INIT
 #define SWIG_PYTHON_2_UNICODE
 
+/*
+ * In C++ we can use a scoped variable to acquire the GIL and then auto release
+ * on leaving scope, making our code a bit cleaner.
+ *
+ * SWIG_PYTHON_THREAD_SCOPED_BLOCK is a macro that SWIG automatically generates
+ * wrapping a class using an RAII pattern to automatically acquire/release
+ * the GIL. See the generated meep-python.cxx for details.
+ *
+ * We could instead just explicitly call SWIG_PYTHON_THREAD_BEGIN_BLOCK and
+ * SWIG_PYTHON_THREAD_END_BLOCK everywhere - but this is error prone since we
+ * have to ensure that SWIG_PYTHON_THREAD_END_BLOCK is called before every
+ * return statement in a method.
+ *
+ * NOTE: This wont work with plain-old C.
+ */
+#define SWIG_PYTHON_THREAD_SCOPED_BLOCK   SWIG_PYTHON_THREAD_BEGIN_BLOCK
+
 #include <complex>
 #include <string>
 
@@ -112,6 +129,7 @@ static PyObject *py_meep_src_time_object() {
 }
 
 static double py_callback_wrap(const meep::vec &v) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     PyObject *pyv = vec2py(v);
     PyObject *pyret = PyObject_CallFunctionObjArgs(py_callback, pyv, NULL);
     double ret = PyFloat_AsDouble(pyret);
@@ -121,6 +139,7 @@ static double py_callback_wrap(const meep::vec &v) {
 }
 
 static std::complex<double> py_amp_func_wrap(const meep::vec &v) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     PyObject *pyv = vec2py(v);
     PyObject *pyret = PyObject_CallFunctionObjArgs(py_amp_func, pyv, NULL);
     double real = PyComplex_RealAsDouble(pyret);
@@ -134,6 +153,7 @@ static std::complex<double> py_amp_func_wrap(const meep::vec &v) {
 static std::complex<double> py_field_func_wrap(const std::complex<double> *fields,
                                                const meep::vec &loc,
                                                void *data_) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     PyObject *pyv = vec2py(loc);
 
     py_field_func_data *data = (py_field_func_data *)data_;
@@ -163,37 +183,35 @@ static std::complex<double> py_field_func_wrap(const std::complex<double> *field
 }
 
 static meep::vec py_kpoint_func_wrap(double freq, int mode, void *user_data) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     PyObject *py_freq = PyFloat_FromDouble(freq);
     PyObject *py_mode = PyInteger_FromLong(mode);
 
     PyObject *py_result = PyObject_CallFunctionObjArgs((PyObject*)user_data, py_freq, py_mode, NULL);
 
-    if (!py_result) {
-        PyErr_PrintEx(0);
-        Py_DECREF(py_freq);
-        Py_DECREF(py_mode);
-        return meep::vec(0, 0, 0);
-    }
+    meep::vec result;
 
-    vector3 v3;
-    if (!pyv3_to_v3(py_result, &v3)) {
+    if (!py_result) {
         PyErr_PrintEx(0);
-        Py_DECREF(py_freq);
-        Py_DECREF(py_mode);
+        result = meep::vec(0, 0, 0);
+    } else {
+        vector3 v3;
+        if (!pyv3_to_v3(py_result, &v3)) {
+            PyErr_PrintEx(0);
+            result = meep::vec(0, 0, 0);
+        } else {
+            result = meep::vec(v3.x, v3.y, v3.z);
+        }
         Py_XDECREF(py_result);
-        return meep::vec(0, 0, 0);
     }
 
-    meep::vec result(v3.x, v3.y, v3.z);
-
     Py_DECREF(py_freq);
     Py_DECREF(py_mode);
-    Py_DECREF(py_result);
-
     return result;
 }
 
 static void _do_master_printf(const char* stream_name, const char* text) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     PyObject *py_stream = PySys_GetObject((char*)stream_name); // arg is non-const on Python2
 
     Py_XDECREF(PyObject_CallMethod(py_stream, "write", "(s)", text));
@@ -248,6 +266,7 @@ static int pyabsorber_to_absorber(PyObject *py_absorber, meep_geom::absorber *a)
 
 // Wrapper for Python PML profile function
 double py_pml_profile(double u, void *f) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     PyObject *func = (PyObject *)f;
     PyObject *d = PyFloat_FromDouble(u);
 
@@ -573,6 +592,7 @@ meep::eigenmode_data *_get_eigenmode(meep::fields *f, double frequency, meep::di
 }
 
 PyObject *_get_eigenmode_Gk(meep::eigenmode_data *emdata) {
+    SWIG_PYTHON_THREAD_SCOPED_BLOCK;
     // Return value: New reference
     PyObject *v3_class = py_vector3_object();
     PyObject *args = Py_BuildValue("(ddd)", emdata->Gk[0], emdata->Gk[1], emdata->Gk[2]);
@@ -594,6 +614,24 @@ void _get_eigenmode(meep::fields *f, double frequency, meep::direction d, const
 #endif
 %}
 
+/*
+ * These methods extensively use the Python C api (especially allocation) and
+ * hence need to hold the GIL (acquire/release) for key parts of their
+ * implementaion/code. Instead, disable threading for these methods by default.
+ *
+ * TODO: If any of these methods are expensive, we can explicitly allow threads
+ * for the expensive blocks of code in these methods.
+ */
+%feature("nothreadallow") _dft_ldos_J;
+%feature("nothreadallow") _dft_ldos_F;
+%feature("nothreadallow") _dft_ldos_ldos;
+%feature("nothreadallow") _get_farfields_array;
+%feature("nothreadallow") _get_farfield;
+%feature("nothreadallow") py_do_harminv;
+%feature("nothreadallow") _get_array_slice_dimensions;
+%feature("nothreadallow") _get_gradient;
+%feature("nothreadallow") _get_dft_array;
+
 %numpy_typemaps(std::complex<double>, NPY_CDOUBLE, int);
 %numpy_typemaps(std::complex<double>, NPY_CDOUBLE, size_t);
 
@@ -867,6 +905,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
     Py_DECREF(swigobj);
 }
 %}
+
 //--------------------------------------------------
 // end typemaps needed for material grid
 //--------------------------------------------------
@@ -1410,6 +1449,11 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
 }
 
 %exception {
+#ifdef MEEP_SWIG_PYTHON_DEBUG
+  // NOTE: You can do fancier things like timing the calls and using that
+  // to track the most expensive calls etc.
+  master_printf("**SWIG**: $symname\n");
+#endif
   try {
     $action
   } catch (std::runtime_error &e) {
@@ -1522,6 +1566,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
 %template(get_dft_fields_array) _get_dft_array<meep::dft_fields>;
 %template(get_dft_force_array) _get_dft_array<meep::dft_force>;
 %template(get_dft_near2far_array) _get_dft_array<meep::dft_near2far>;
+
 %template(FragmentStatsVector) std::vector<meep_geom::fragment_stats>;
 %template(DftDataVector) std::vector<meep_geom::dft_data>;
 %template(VolumeVector) std::vector<meep::volume>;