Add argument for spline boundary conditions

doc/changes/2114.misc

@@ -0,0 +1,2 @@
+Added new argument bc_type to take boundary conditions when creating QobjEvo
+

qutip/core/coefficient.py

@@ -53,7 +53,8 @@ def _return(base, **kwargs):
 
 
 def coefficient(base, *, tlist=None, args={}, args_ctypes={},
-                order=3, compile_opt=None, function_style=None, **kwargs):
+                order=3, compile_opt=None, function_style=None,
+                boundary_conditions=None, **kwargs):
     """Build ``Coefficient`` for time dependent systems:
 
     ```
@@ -149,6 +150,9 @@ def f2_t(t, args):
     compile_opt : CompilationOptions, optional
         Sets of options for the compilation of string based coefficients.
 
+    bc_type: 2-tupule or None, optional
+        Specify boundary conditions for spline interpolation.
+
     **kwargs
         Extra arguments to pass the the coefficients.
     """
@@ -159,6 +163,7 @@ def f2_t(t, args):
         'order': order,
         'compile_opt': compile_opt,
         'function_style': function_style,
+        'bc_type': bc_type,
     })
 
     for type_ in coefficient_builders:

qutip/core/cy/coefficient.pyx

@@ -387,14 +387,23 @@ cdef class InterCoefficient(Coefficient):
     order : int
         Order of the interpolation. Order ``0`` uses the previous (i.e. left)
         value. The order will be reduced to ``len(tlist) - 1`` if it is larger.
+
+    bc_type : 2-Tuple or None, optional
+        Boundary conditions for spline evaluation. Default value is `None` to
+        enable automatically choosing boundary conditions. Otherwise, a two-length
+        tuple of form (`deriv_l`, `deriv_r`) must be passed. Here, `deriv_l` denotes
+        the boundary conditions at ```t = 0``` and deriv_r, the conditions at time
+        ```t```. The condition specified must be an iterable pair of form (`order`,
+        `value`). Refer to Scipy's documentation for further details:
+        https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.make_interp_spline.html
     """
     cdef int order
     cdef double dt
     cdef double[::1] tlist
     cdef complex[:, :] poly
     cdef object np_arrays
 
-    def __init__(self, coeff_arr, tlist, int order, **_):
+    def __init__(self, coeff_arr, tlist, int order, bc_type, **_):
         tlist = np.array(tlist, dtype=np.float64)
         coeff_arr = np.array(coeff_arr, dtype=np.complex128)
 
@@ -418,7 +427,7 @@ cdef class InterCoefficient(Coefficient):
         elif order >= 2:
             # Use scipy to compute the spline and transform it to polynomes
             # as used in scipy's PPoly which is easier for us to use.
-            spline = make_interp_spline(tlist, coeff_arr, k=order)
+            spline = make_interp_spline(tlist, coeff_arr, k=order, bc_type=bc_type)
             # Scipy can move knots, we add them to tlist
             tlist = np.sort(np.unique(np.concatenate([spline.t, tlist])))
             a = np.arange(spline.k+1)

qutip/core/cy/qobjevo.pyx

@@ -99,6 +99,17 @@ cdef class QobjEvo:
         ``qutip.settings.core["function_coefficient_style"]``
         is used. Otherwise the supplied value overrides the global setting.
 
+
+    bc_type : 2-Tuple or None, optional
+        Boundary conditions for spline evaluation. Default value is `None` to
+        enable automatically choosing boundary conditions. Otherwise, a two-length
+        tuple of form (`deriv_l`, `deriv_r`) must be passed. Here, `deriv_l` denotes
+        the boundary conditions at ```t = 0``` and deriv_r, the conditions at time
+        ```t```. The condition specified must be an iterable pair of form (`order`,
+        `value`). Refer to Scipy's documentation for further details:
+        https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.make_interp_spline.html
+
+
     Attributes
     ----------
     dims : list
@@ -182,7 +193,7 @@ cdef class QobjEvo:
     """
     def __init__(QobjEvo self, Q_object, args=None, tlist=None,
                  order=3, copy=True, compress=True,
-                 function_style=None):
+                 function_style=None, bc_type=None):
         if isinstance(Q_object, QobjEvo):
             self.dims = Q_object.dims.copy()
             self.shape = Q_object.shape
@@ -217,29 +228,29 @@ cdef class QobjEvo:
                 self.elements.append(
                     self._read_element(
                         op, copy=copy, tlist=tlist, args=args, order=order,
-                        function_style=function_style
+                        function_style=function_style, bc_type=bc_type
                     )
                 )
         else:
             self.elements.append(
                 self._read_element(
                     Q_object, copy=copy, tlist=tlist, args=args, order=order,
-                    function_style=function_style
+                    function_style=function_style, bc_type=bc_type
                 )
             )
 
         if compress:
             self.compress()
 
-    def _read_element(self, op, copy, tlist, args, order, function_style):
+    def _read_element(self, op, copy, tlist, args, order, function_style, bc_type):
         """ Read a Q_object item and return an element for that item. """
         if isinstance(op, Qobj):
             out = _ConstantElement(op.copy() if copy else op)
             qobj = op
         elif isinstance(op, list):
             out = _EvoElement(
                 op[0].copy() if copy else op[0],
-                coefficient(op[1], tlist=tlist, args=args, order=order)
+                coefficient(op[1], tlist=tlist, args=args, order=order, bc_type=bc_type)
             )
             qobj = op[0]
         elif isinstance(op, _BaseElement):

qutip/tests/core/test_coefficient.py

@@ -372,6 +372,10 @@ def test_CoeffFromScipy():
     from_scipy = coefficient(interp.make_interp_spline(tlist, y, k=3))
     _assert_eq_over_interval(coeff, from_scipy, rtol=1e-8, inside=True)
 
+    coeff = coefficient(y, tlist=tlist, order=3, bc_type=([(2, 0.0)],[(2, 0.0)]))
+    from_scipy = coefficient(interp.make_interp_spline(tlist, y, k=3, bc_type="natural"))
+    _assert_eq_over_interval(coeff, from_scipy, rtol=1e-8, inside=True)
+
 
 @pytest.mark.parametrize('map_func', [
     pytest.param(qutip.solver.parallel.parallel_map, id='parallel_map'),