[Cython] Generalized interface to Func1

It is now possible to use any callable Python object when the Cantera C++ core
requires an instance of Func1.

doc/sphinx/cython/index.rst

@@ -9,3 +9,4 @@ Contents:
    thermo
    kinetics
    transport
+   zerodim

doc/sphinx/cython/zerodim.rst

@@ -0,0 +1,9 @@
+.. py:currentmodule:: cantera
+
+Zero-Dimensional Reactor Networks
+=================================
+
+Defining Functions
+------------------
+
+.. autoclass:: Func1

interfaces/cython/cantera/_cantera.pxd

@@ -27,6 +27,16 @@ cdef extern from "cantera/thermo/mix_defs.h":
     cdef int kinetics_type_interface "Cantera::cInterfaceKinetics"
     cdef int kinetics_type_edge "Cantera::cEdgeKinetics"
 
+
+cdef extern from "funcWrapper.h":
+    ctypedef double (*callback_wrapper)(double, void*, void**)
+    cdef int translate_exception()
+
+    cdef cppclass CxxFunc1 "Func1Py":
+        CxxFunc1(callback_wrapper, void*)
+        double eval(double) except +translate_exception
+
+
 cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera":
     cdef cppclass CxxThermoPhase "Cantera::ThermoPhase":
         CxxThermoPhase()
@@ -118,6 +128,7 @@ cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera":
 cdef extern from "cantera/thermo/IdealGasPhase.h":
     cdef cppclass CxxIdealGasPhase "Cantera::IdealGasPhase"
 
+
 cdef extern from "cantera/thermo/SurfPhase.h":
     cdef cppclass CxxSurfPhase "Cantera::SurfPhase":
         CxxSurfPhase()
@@ -192,17 +203,10 @@ cdef extern from "cantera/equil/equil.h" namespace "Cantera":
 cdef extern from "cantera/equil/vcs_MultiPhaseEquil.h" namespace "Cantera":
     int vcs_equilibrate(CxxMultiPhase&, char*, int, int, int, double, int, int, int)
 
-cdef extern from "cantera/numerics/Func1.h":
-    cdef cppclass CxxFunc1 "Cantera::Func1":
-        CxxFunc1()
-        double eval(double)
-        string write(string)
-
-    cdef cppclass CxxSin1 "Cantera::Sin1":
-        CxxSin1(double)
 
 cdef extern from "cantera/zeroD/ReactorBase.h" namespace "Cantera":
     cdef cppclass CxxWall "Cantera::Wall"
+
     cdef cppclass CxxReactorBase "Cantera::ReactorBase":
         CxxReactorBase()
         void setThermoMgr(CxxThermoPhase&)
@@ -354,6 +358,8 @@ cdef class Mixture:
 
 cdef class Func1:
     cdef CxxFunc1* func
+    cdef object callable
+    cdef object exception
 
 cdef class ReactorBase:
     cdef CxxReactorBase* rbase

interfaces/cython/cantera/func1.pyx

@@ -1,20 +1,73 @@
+import sys
+
+cdef double func_callback(double t, void* obj, void** err):
+    """
+    This function is called from C/C++ to evaluate a `Func1` object *obj*,
+    returning the value of the function at *t*. If an exception occurs while
+    evaluating the function, the Python exception info is saved in the
+    two-element array *err*.
+    """
+    try:
+        return (<Func1>obj).callable(t)
+    except Exception as e:
+        exc_type, exc_value = sys.exc_info()[:2]
+
+        # Stash the exception info to prevent it from being garbage collected
+        (<Func1>obj).exception = exc_type, exc_value
+        err[0] = <void*>exc_type
+        err[1] = <void*>exc_value
+        return 0.0
+
+
 cdef class Func1:
-    def __cinit__(self, *args, **kwargs):
-        self.func = NULL # derived classes should allocate this object
+    """
+    This class is used as a wrapper for a function of one variable, i.e.
+    :math:`y = f(t)`, that is defined in Python and can be called by the
+    Cantera C++ core. `Func1` objects are constructed from callable Python
+    objects, e.g. functions or classes which implement the `__call__` method::
+
+        >>> f1 = Func1(math.sin)
+        >>> f1(math.pi/4)
+        0.7071067811865475
+
+        >>> f2 = Func1(lambda t: t**2 + 1)
+        >>> f2(3)
+        10
+
+        >>> class Multiplier(object):
+        ...     def __init__(self, factor):
+        ...         self.factor = factor
+        ...     def __call__(self, t):
+        ...         return self.factor * t
+        >>> f3 = Func1(Multiplier(5))
+        >>> f3(6)
+        30.0
+
+    For simplicity, constant functions can be defined by passing the constant
+    value directly::
+
+        >>> f4 = Func1(2.5)
+        >>> f4(0.1)
+        2.5
+
+    Note that all methods which accept `Func1` objects will also accept the
+    callable object and create the wrapper on their own, so it is generally
+    unnecessary to explicitly create a `Func1` object.
+    """
+    def __cinit__(self, c):
+        self.exception = None
+        if isinstance(c, (float, int)):
+            self.callable = lambda t: c
+        elif hasattr(c, '__call__'):
+            self.callable = c
+        else:
+            raise TypeError('Func1 must be constructed from a number or a '
+                            'callable object')
+
+        self.func = new CxxFunc1(func_callback, <void*>self)
 
     def __dealloc__(self):
         del self.func
 
-    def __init__(self, *args, **kwargs):
-        assert self.func != NULL
-
     def __call__(self, t):
         return self.func.eval(t)
-
-    def __str__(self):
-        cdef bytes s = self.func.write(stringify("t")).c_str()
-        return s.decode()
-
-cdef class Sin1(Func1):
-    def __cinit__(self, freq):
-        self.func = <CxxFunc1*>(new CxxSin1(freq))

interfaces/cython/cantera/funcWrapper.h

@@ -0,0 +1,72 @@
+#ifndef CT_CYTHON_FUNC_WRAPPER
+#define CT_CYTHON_FUNC_WRAPPER
+
+#include "cantera/numerics/Func1.h"
+
+typedef double (*callback_wrapper)(double, void*, void**);
+
+// A C++ exception that holds a Python exception so that it can be re-raised
+// by translate_exception()
+class CallbackError : public Cantera::CanteraError
+{
+public:
+    CallbackError(void* type, void* value) :
+        m_type((PyObject*) type),
+        m_value((PyObject*) value)
+    {}
+    PyObject* m_type;
+    PyObject* m_value;
+};
+
+
+// A function of one variable implemented as a callable Python object
+class Func1Py : public Cantera::Func1
+{
+public:
+    Func1Py(callback_wrapper callback, void* pyobj) :
+        m_callback(callback),
+        m_pyobj(pyobj)
+    {
+    }
+
+    double eval(double t) const {
+        void* err[2] = {0, 0};
+        double y = m_callback(t, m_pyobj, err);
+        if (err[0]) {
+            throw CallbackError(err[0], err[1]);
+        }
+        return y;
+    }
+
+private:
+    callback_wrapper m_callback;
+    void* m_pyobj;
+};
+
+
+// Translate C++ Exceptions generated by Cantera to appropriate Python
+// exceptions. Used with Cython function declarations, e.g:
+//     cdef double eval(double) except +translate_exception
+inline int translate_exception()
+{
+    try
+    {
+        if (!PyErr_Occurred()) {
+            // Let the latest Python exception pass through and ignore the
+            // current one.
+            throw;
+        }
+    } catch (CallbackError& exn) {
+        // Re-raise a Python exception generated in a callback
+        PyErr_SetObject(exn.m_type, exn.m_value);
+    } catch (const std::out_of_range& exn) {
+        PyErr_SetString(PyExc_IndexError, exn.what());
+    } catch (const std::exception& exn) {
+        PyErr_SetString(PyExc_Exception, exn.what());
+    } catch (...) {
+        PyErr_SetString(PyExc_Exception, "Unknown exception");
+    }
+    return 0;
+}
+
+#endif

interfaces/cython/cantera/test/test_func1.py

@@ -5,8 +5,37 @@ import cantera as ct
 from . import utilities
 
 class TestFunc1(utilities.CanteraTest):
-    def test_sin(self):
-        f = ct.Sin1(3)
+    def test_function(self):
+        f = ct.Func1(np.sin)
         self.assertNear(f(0), np.sin(0))
-        self.assertNear(f(0.1), np.sin(3*0.1))
-        self.assertNear(f(0.7), np.sin(3*0.7))
+        self.assertNear(f(0.1), np.sin(0.1))
+        self.assertNear(f(0.7), np.sin(0.7))
+
+    def test_lambda(self):
+        f = ct.Func1(lambda t: np.sin(t)*np.sqrt(t))
+        for t in [0.1, 0.7, 4.5]:
+            self.assertNear(f(t), np.sin(t)*np.sqrt(t))
+
+    def test_callable(self):
+        class Multiplier(object):
+            def __init__(self, factor):
+                self.factor = factor
+            def __call__(self, t):
+                return self.factor * t
+
+        m = Multiplier(8.1)
+        f = ct.Func1(m)
+        for t in [0.1, 0.7, 4.5]:
+            self.assertNear(f(t), 8.1*t)
+
+    def test_constant(self):
+        f = ct.Func1(5)
+        for t in [0.1, 0.7, 4.5]:
+            self.assertNear(f(t), 5)
+
+    def test_failure(self):
+        def fails(t):
+            raise ValueError('bad')
+
+        f = ct.Func1(fails)
+        self.assertRaises(ValueError, f, 0.1)