Fix interpolate_by_finite_diff and native impl. thereof (#25)

* Native interpolate_by_finite_diff
* Updates to CI scripts

.drone.yml

@@ -1,3 +1,7 @@
+clone:
+  depth: 2
+  recursive: true
+
 build:
   image: bjodah/bjodahimg18dev:v1.0
   environment:
@@ -10,6 +14,7 @@ build:
     - (cd examples/; make -B CXX=g++ EXTRA_COMPILE_ARGS="-DNDEBUG -O2")
     - (cd tests/; make -B CXX=clang++-6.0 EXTRA_COMPILE_ARGS="-fsanitize=address -O0 -g")
     - (cd tests/; make -B CXX=g++ EXTRA_COMPILE_ARGS="-DNDEBUG -O2")
+    - (cd tests/; make -B CXX=g++ EXTRA_COMPILE_ARGS="-DNDEBUG -O2 -DFINITEDIFF_OPENMP -fopenmp")
     - (cd tests/; make -f fortran_tests.mk)
     - ./scripts/ci.sh finitediff
     - ./scripts/generate_docs.sh

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "external/newton_interval"]
+	path = external/newton_interval
+	url = git://github.com/bjodah/newton_interval
+[submodule "finitediff/external/newton_interval"]
+	path = finitediff/external/newton_interval
+	url = git://github.com/bjodah/newton_interval

CHANGES.rst

@@ -1,3 +1,8 @@
+v0.6.0
+======
+- New C API (prefixed function names with ``finitediff_``)
+- New C function: finitediff_interpolate_by_finite_diff (optionally OpenMP parallelized)
+
 v0.5.4
 ======
 - Re-release for Zenodo

MANIFEST.in

@@ -1,4 +1,4 @@
-graft finitediff/external
+graft finitediff/external/newton_interval
 include finitediff/include/finitediff_c.h
 include finitediff/include/finitediff_c.pxd
 include finitediff/include/finitediff_templated.hpp

README.rst

@@ -44,6 +44,10 @@ from the Python bindings).
 The user may also manually generate the corresponding weights. (see
 ``calculate_weights``)
 
+Finitediff can be conditionally compiled to make ``finitediff_interpolate_by_finite_diff``
+multithreaded. Then the number of threads used is set through the environment variable
+``FINITEDIFF_NUM_THREADS`` (or ``OMP_NUM_THREADS``).
+
 
 Documentation
 -------------
@@ -141,19 +145,6 @@ You need either a C, C++ or a Fortran 90 compiler. On debian based linux systems
 
 See `setup.py <setup.py>`_ for optional (Python) dependencies.
 
-Notes
-=====
-There is a git subtree under finitediff, update through::
-
-    git subtree pull --prefix finitediff/external/newton_interval newton_interval master --squash
-
-
-where the repo "newton_interval" is https://github.com/bjodah/newton_interval.git
-
-First time you need to add it::
-
-    git subtree add --prefix finitediff/external/newton_interval git://github.com/bjodah/newton_interval master
-
 
 References
 ==========

appveyor.yml

@@ -21,6 +21,7 @@ install:
     - cmd: conda update --yes --quiet conda
 
     - cmd: set PYTHONUNBUFFERED=1
+    - cmd: git submodule update --init --recursive
 
 # Skip .NET project specific build phase.
 build: off

finitediff/_finitediff_c.pyx

@@ -5,7 +5,7 @@ cimport numpy as cnp
 import numpy as np
 
 from newton_interval cimport get_interval, get_interval_from_guess
-from finitediff_c cimport apply_fd, calculate_weights
+from finitediff_c cimport finitediff_calc_and_apply_fd, finitediff_calculate_weights, finitediff_interpolate_by_finite_diff
 
 
 def get_weights(grid, double xtgt, int n=-1, int maxorder=0):
@@ -35,7 +35,7 @@ def get_weights(grid, double xtgt, int n=-1, int maxorder=0):
         n = xarr.size
     cdef cnp.ndarray[cnp.float64_t, ndim=2, mode='fortran'] c = \
         np.empty((n, maxorder+1), order='F')
-    calculate_weights(&c[0, 0], n, &xarr[0], n, maxorder, xtgt)
+    finitediff_calculate_weights(&c[0, 0], n, &xarr[0], n, maxorder, xtgt)
     return c
 
 
@@ -94,7 +94,7 @@ def derivatives_at_point_by_finite_diff(
     cdef cnp.ndarray[cnp.float64_t, ndim=1] yout = np.empty((maxorder+1)*nsets)
     if xarr.size < maxorder+1:
         raise ValueError("xdata too short for requested derivative order")
-    apply_fd(&yout[0], xarr.size, nsets, maxorder, xarr.size, &xarr[0], &yarr[0], xarr.size, xtgt)
+    finitediff_calc_and_apply_fd(&yout[0], xarr.size, nsets, maxorder, xarr.size, &xarr[0], &yarr[0], xarr.size, xtgt)
     if reshape is None:
         reshape = ydata.ndim != 1
     if reshape:
@@ -163,38 +163,32 @@ def interpolate_by_finite_diff(
     """
     ydata = np.asarray(ydata)
     xtgts = np.asarray(xtgts)
-    cdef int nin = ntail+nhead
-    cdef int nout = xtgts.size
-    cdef cnp.ndarray[cnp.float64_t, ndim=1] xarr = np.ascontiguousarray(grid, dtype=np.float64)
-    cdef cnp.ndarray[cnp.float64_t, ndim=1] tgts = np.ascontiguousarray(xtgts, dtype=np.float64)
-    cdef cnp.ndarray[cnp.float64_t, ndim=1] yarr = np.ascontiguousarray(np.ravel(ydata, order=yorder), dtype=np.float64)
-    if yarr.size % xarr.size:
+    cdef:
+        int flag
+        int nin = ntail+nhead
+        int nout = xtgts.size
+        cnp.ndarray[cnp.float64_t, ndim=1] xgrd = np.ascontiguousarray(grid, dtype=np.float64)
+        cnp.ndarray[cnp.float64_t, ndim=1] tgts = np.ascontiguousarray(xtgts, dtype=np.float64)
+        cnp.ndarray[cnp.float64_t, ndim=1] yarr = np.ascontiguousarray(np.ravel(ydata, order=yorder), dtype=np.float64)
+        int nsets = yarr.size // xgrd.size
+        cnp.ndarray[cnp.float64_t, ndim=1] yout = np.zeros(
+            (nout*nsets*(maxorder+1)), order='C', dtype=np.float64)
+
+    if yarr.size % xgrd.size:
         raise ValueError("Incompatible shapes: grid & ydata")
-    cdef int nsets = yarr.size // xarr.size
-    cdef cnp.ndarray[cnp.float64_t, ndim=1] yout = np.zeros(
-        (nout*nsets*(maxorder+1)), order='C', dtype=np.float64)
-    cdef int i, j=0
 
-    if xarr.size < min(ntail, nhead) + 1:
+    flag = finitediff_interpolate_by_finite_diff(
+        <double*>yout.data, nout, nsets, maxorder, nsets*(maxorder+1), maxorder+1,
+        ntail, nhead, <double*>xgrd.data, xgrd.size, <double*>yarr.data, xgrd.size,
+        <double*>tgts.data
+    )
+    if flag == 1:
+        raise ValueError("Bad alloc")
+    elif flag == 2:
         raise ValueError("grid is too small")
-    if nhead+ntail < maxorder+1:
-        raise ValueError("nhead+ntail < maxorder+1")
-
-    for i in range(nout):
-        j = max(0, get_interval_from_guess(
-            &xarr[0], xarr.size, tgts[i], j))
-        j = min(j, xarr.size - nin)
-        apply_fd(
-            &yout[i*nsets*(maxorder+1)],
-            maxorder+1,
-            nsets,
-            maxorder,
-            nin,
-            &xarr[j],
-            &yarr[j],
-            xarr.size,
-            tgts[i]
-        )
+    if flag == 4:
+        raise ValueError("too few points")
+
     if reshape is None:
         reshape = ydata.ndim != 1
     if reshape: