Merge pull request #12 from bjodah/refactor

Deprecated populate_weights, introduced calculate_weights

.gitignore

@@ -26,3 +26,5 @@ finitediff/_finitediff_templated.cpp
 *.egg-info
 tests/test_finitediff_templated
 examples/demo
+tests/catch.hpp
+examples/a.out

CHANGES.rst

@@ -1,3 +1,25 @@
+v0.3.0
+======
+- Refactored ``finitediff::generate_weights`` (part of official API from v0.3.0)
+- Deprecated ``finitediff::populate_weights``
+- Introduced ``finitediff::calculate_weights`` (to replace populate_weights)
+
+v0.2.5
+======
+- Fixes to Python distribution (include C++ header)
+
+v0.2.4
+======
+- C++11 function ``generate_weights`` provisionally provided
+
+v0.2.3
+======
+- conda recipe related fixes
+
+v0.2.2
+======
+- fixes to setup.py
+
 v0.2.1
 ======
 - More robust setup.py, source distributions better tested.

README.rst

@@ -30,11 +30,11 @@ arbitrary derivative order. Python_ bindings are provided.
 Capabilities
 ============
 finitediff currently provides callbacks for estimation of derivatives
-or interpolation either at a single point or over an array (available 
-from the Python bindings). 
+or interpolation either at a single point or over an array (available
+from the Python bindings).
 
 The user may also manually generate the corresponding weights. (see
-``populate_weights``) 
+``populate_weights``)
 
 
 Documentation
@@ -54,29 +54,30 @@ Generating finite difference weights is simple using C++11:
    #include <vector>
    #include <string>
    #include <iostream>
-   
+
    int main(){
-       const unsigned maxorder = 2;
-       std::vector<std::string> labels {"Zeroth order", "First order", "Second order"};
+       const unsigned max_deriv = 2;
+       std::vector<std::string> labels {"Zeroth derivative (interpolation)", "First derivative", "Second derivative"};
        std::vector<double> x {0, 1, -1, 2, -2};  // Fourth order of accuracy
-       auto coeffs = finitediff::generate_weights(0.0, x, maxorder);
-       for (unsigned order = 0; order <= maxorder; order++){
-           std::cout << labels[order] << ": ";
+       auto coeffs = finitediff::generate_weights(x, max_deriv);
+       for (unsigned deriv_i = 0; deriv_i <= max_deriv; deriv_i++){
+           std::cout << labels[deriv_i] << ": ";
            for (unsigned idx = 0; idx < x.size(); idx++){
-               std::cout << coeffs[order*x.size() + idx] << " ";
+               std::cout << coeffs[deriv_i*x.size() + idx] << " ";
            }
            std::cout << std::endl;
        }
    }
-
+
+
 ::
 
    $ cd examples/
    $ g++ -std=c++11 demo.cpp -I../include
    $ ./a.out
-   Zeroth order: 1 -0 0 0 -0 
-   First order: -0 0.666667 -0.666667 -0.0833333 0.0833333 
-   Second order: -2.5 1.33333 1.33333 -0.0833333 -0.0833333 
+   Zeroth derivative (interpolation): 1 -0 0 0 -0
+   First derivative: -0 0.666667 -0.666667 -0.0833333 0.0833333
+   Second derivative: -2.5 1.33333 1.33333 -0.0833333 -0.0833333
 
 
 and of course using the python bindings:
@@ -88,7 +89,7 @@ and of course using the python bindings:
    >>> c = get_weights(np.array([-1., 0, 1]), 0, maxorder=1)
    >>> np.allclose(c[:, 1], [-.5, 0, .5])
    True
-   
+
 
 see the ``examples/`` directory for more examples.
 
@@ -168,7 +169,7 @@ http://dx.doi.org/10.1137/S0036144596322507
       publisher={SIAM}
       doi={10.1137/S0036144596322507}
     }
-    
+
 
 Which is based on an article of the same author:
 

examples/demo.cpp

@@ -4,14 +4,14 @@
 #include <iostream>
 
 int main(){
-    const unsigned maxorder = 2;
-    std::vector<std::string> labels {"Zeroth order", "First order", "Second order"};
+    const unsigned max_deriv = 2;
+    std::vector<std::string> labels {"Zeroth derivative (interpolation)", "First derivative", "Second derivative"};
     std::vector<double> x {0, 1, -1, 2, -2};  // Fourth order of accuracy
-    auto coeffs = finitediff::generate_weights(0.0, x, maxorder);
-    for (unsigned order = 0; order <= maxorder; order++){
-        std::cout << labels[order] << ": ";
+    auto coeffs = finitediff::generate_weights(x, max_deriv);
+    for (unsigned deriv_i = 0; deriv_i <= max_deriv; deriv_i++){
+        std::cout << labels[deriv_i] << ": ";
         for (unsigned idx = 0; idx < x.size(); idx++){
-            std::cout << coeffs[order*x.size() + idx] << " ";
+            std::cout << coeffs[deriv_i*x.size() + idx] << " ";
         }
         std::cout << std::endl;
     }

include/finitediff_templated.hpp

@@ -8,6 +8,60 @@
 
 namespace finitediff {
 
+    template <typename Real_t>
+    void calculate_weights(const Real_t * const __restrict__ grid, const int len_g,
+                           const int max_deriv, Real_t * const __restrict__ c, const Real_t around=0) {
+        // Parameters
+        // ----------
+        // grid[len_g]: array with grid point locations
+        // len_g: length of grid
+        // c[len_g, max_deriv+1]: weights of order 0 to max_deriv (output argument, contiguous memory, column major order)
+        // max_deriv: highest derivative.
+        // around: location where approximations are to be accurate
+        //
+        // References
+        // ----------
+        // Generation of Finite Difference Formulas on Arbitrarily
+        // Spaced Grids, Bengt Fornberg,
+        // Mathematics of compuation, 51, 184, 1988, 699-706
+
+        Real_t c1, c4, c5;
+        c1 = 1;
+        c4 = grid[0] - around;
+        for (int i=0; i < len_g*(max_deriv+1); ++i)
+            c[i] = 0;
+        c[0] = 1;
+        for (int i=1; i < len_g; ++i){
+            const int mn = (i < max_deriv) ? i : max_deriv; // min(i, max_deriv)
+            Real_t c2 = 1;
+            c5 = c4;
+            c4 = grid[i] - around;
+            for (int j=0; j<i; ++j){
+                const Real_t c3 = grid[i] - grid[j];
+                const Real_t c3_r = 1/c3;
+                c2 = c2*c3;
+                if (j == i-1){
+                    const Real_t c2_r = 1/c2;
+                    for (int k=mn; k>=1; --k){
+                        const Real_t tmp1 = c[i - 1 + (k-1)*len_g];
+                        const Real_t tmp2 = c[i - 1 + k*len_g];
+                        c[i + k*len_g] = c1*(k*tmp1 - c5*tmp2)*c2_r;
+                    }
+                    c[i] = -c1*c5*c[i-1]*c2_r;
+                }
+                for (int k=mn; k>=1; --k){
+                    const Real_t tmp1 = c[j + k*len_g];
+                    const Real_t tmp2 = c[j + (k-1)*len_g];
+                    c[j + k*len_g] = (c4*tmp1 - k*tmp2)*c3_r;
+                }
+                c[j] = c4*c[j]*c3_r;
+            }
+            c1 = c2;
+        }
+    }
+
+
+    // populate_weights is deprecated due to counter-intuitive parameter "nd"
     template <typename Real_t>
     void populate_weights(const Real_t z, const Real_t * const __restrict__ x, const int nd,
                           const int m, Real_t * const __restrict__ c) {
@@ -25,34 +79,36 @@ namespace finitediff {
         // Spaced Grids, Bengt Fornberg,
         // Mathematics of compuation, 51, 184, 1988, 699-706
 
-        Real_t c1, c2, c3, c4, c5;
+        Real_t c1, c4, c5;
         c1 = 1;
         c4 = x[0] - z;
         for (int i=0; i < (nd+1)*(m+1); ++i)
             c[i] = 0;
         c[0] = 1;
         for (int i=1; i <= nd; ++i){
             const int mn = (i < m) ? i : m; // min(i, m)
-            c2 = 1;
+            Real_t c2 = 1;
             c5 = c4;
             c4 = x[i] - z;
             for (int j=0; j<i; ++j){
-                c3 = x[i] - x[j];
+                const Real_t c3 = x[i] - x[j];
+                const Real_t c3_r = 1/c3;
                 c2 = c2*c3;
                 if (j == i-1){
+                    const Real_t c2_r = 1/c2;
                     for (int k=mn; k>=1; --k){
                         const Real_t tmp1 = c[i - 1 + (k-1)*(nd+1)];
                         const Real_t tmp2 = c[i - 1 + k*(nd+1)];
-                        c[i + k*(nd+1)] = c1*(k*tmp1 - c5*tmp2)/c2;
+                        c[i + k*(nd+1)] = c1*(k*tmp1 - c5*tmp2)*c2_r;
                     }
-                    c[i] = -c1*c5*c[i-1]/c2;
+                    c[i] = -c1*c5*c[i-1]*c2_r;
                 }
                 for (int k=mn; k>=1; --k){
                     const Real_t tmp1 = c[j + k*(nd+1)];
                     const Real_t tmp2 = c[j + (k-1)*(nd+1)];
-                    c[j + k*(nd+1)] = (c4*tmp1 - k*tmp2)/c3;
+                    c[j + k*(nd+1)] = (c4*tmp1 - k*tmp2)*c3_r;
                 }
-                c[j] = c4*c[j]/c3;
+                c[j] = c4*c[j]*c3_r;
             }
             c1 = c2;
         }
@@ -65,7 +121,7 @@ namespace finitediff {
                   const Real_t xtgt,
                   Real_t * const __restrict__ out){
         Real_t * const c = new Real_t[nin * (maxorder+1)];
-        populate_weights<Real_t>(xtgt, xdata, nin-1, maxorder, c);
+        calculate_weights<Real_t>(xdata, nin, maxorder, c, xtgt);
         for (int j=0; j <= maxorder; ++j){
             out[j] = 0;
             for (int i=0; i<nin; ++i)
@@ -78,12 +134,14 @@ namespace finitediff {
 // Pre-processor macro __cplusplus == 201103L in ISO C++11 compliant compilers. (e.g. GCC >= 4.7.0)
 #if __cplusplus > 199711L
     template<typename Real_t, template<typename, typename...> class Cont, typename... Args>
-    Cont<Real_t, Args...> generate_weights(const Real_t z, const Cont<Real_t, Args...>& x, const unsigned maxorder){
-        Cont<Real_t, Args...> coeffs(x.size()*(maxorder+1));
-        if (x.size() < maxorder + 1){
+    Cont<Real_t, Args...> generate_weights(const Cont<Real_t, Args...>& x, int maxorder=-1, const Real_t around=0){
+        // Cont<Real_t, Args...> must have contiguous memory storage (e.g. std::vector)
+        const unsigned maxorder_ = (maxorder < 0) ? (x.size()+1)/2 : maxorder;
+        Cont<Real_t, Args...> coeffs(x.size()*(maxorder_+1));
+        if (x.size() < maxorder_ + 1){
             throw std::logic_error("size of x insufficient");
         }
-        populate_weights<Real_t>(z, &x[0], x.size() - 1, maxorder, &coeffs[0]);
+        calculate_weights<Real_t>(&x[0], x.size(), maxorder_, &coeffs[0], around);
         return coeffs;
     }
 #endif

tests/test_finitediff_templated.cpp

@@ -9,7 +9,7 @@ T inline abs_(T v) { return v < 0 ? -v : v; }
 TEST_CASE( "weight correctness", "finitediff::generate_weights") {
 
     std::vector<double> x3 {-1, 0, 1};
-    auto coeffs3 = finitediff::generate_weights(0.0, x3, 2);
+    auto coeffs3 = finitediff::generate_weights(x3);
     REQUIRE( coeffs3.size() == 3*3);
 
     // Zeroth order
@@ -29,7 +29,7 @@ TEST_CASE( "weight correctness", "finitediff::generate_weights") {
 
 
     std::vector<double> x5 {-2, -1, 0, 1, 2};
-    auto coeffs5 = finitediff::generate_weights(0.0, x5, 2);
+    auto coeffs5 = finitediff::generate_weights(x5, 2);
     REQUIRE( coeffs5.size() == 5*3);
 
     // Zeroth order