Merge branch 'feature/numpyvector-performance' into 'master'

[bugfix][NumPyVector] capture raw pointer to avoid performance penalty.

See merge request core/dune-common!1041

dune/python/common/numpyvector.hh

 #ifndef DUNE_PYTHON_COMMON_NUMPYVECTOR_HH
 #define DUNE_PYTHON_COMMON_NUMPYVECTOR_HH
 
+#include <dune/common/exceptions.hh>
 #include <dune/common/densevector.hh>
 #include <dune/common/ftraits.hh>
 
@@ -67,17 +68,21 @@ namespace Dune
         : array_( pybind11::buffer_info( nullptr, sizeof( T ),
                   pybind11::format_descriptor< T >::value, 1, { size }, { sizeof( T ) } )
                 ),
+          dataPtr_( static_cast< value_type * >( array_.request(true).ptr ) ),
           size_(size)
       {}
 
       NumPyVector ( pybind11::buffer buf )
         : array_( buf ),
+          dataPtr_( nullptr ),
           size_( 0 )
       {
         pybind11::buffer_info info = buf.request();
         if (info.ndim != 1)
           DUNE_THROW( InvalidStateException, "NumPyVector can only be created from one-dimensional array" );
         size_ = info.shape[0];
+
+        dataPtr_ = static_cast< value_type * >( array_.request(true).ptr );
       }
 
       NumPyVector ( const This &other ) = delete;
@@ -109,11 +114,13 @@ namespace Dune
 
       inline const value_type *data () const
       {
-        return static_cast< const value_type * >( const_cast<pybind11::array_t< T >&>(array_).request(false).ptr );
+        assert( dataPtr_ );
+        return dataPtr_;
       }
       inline value_type *data ()
       {
-        return static_cast< value_type * >( array_.request(true).ptr );
+        assert( dataPtr_ );
+        return dataPtr_;
       }
       pybind11::array_t< T > &coefficients()
       {
@@ -133,8 +140,9 @@ namespace Dune
         return size_;
       }
 
-    private:
+    protected:
       pybind11::array_t< T > array_;
+      value_type* dataPtr_;
       size_type size_;
     };
 

dune/python/test/pythontests.py

 from io import StringIO
+from numpy import array
+
 classACode="""
 struct MyClassA {
   MyClassA(int a,int b) : a_(a), b_(b) {}
+  double b() { return b_; }
   int a_,b_;
 };
 """
 classBCode="""
+#include <iostream>
 #include <cmath>
+#include <numeric>
+#include <dune/python/common/numpyvector.hh>
 template <class T> struct MyClassB {
-  MyClassB(T &t, int p) : a_(std::pow(t.a_,p)), b_(std::pow(t.b_,p)) {}
-  int a_,b_;
+  // the following fails to work correctly
+  // MyClassB(T &t, int p, const Dune::Python::NumPyVector<double> &b)
+  // - with 'const NumPyVector<double> b_;':  call to deleter ctor
+  // - with 'const NumPyVector<double> &b_;': segfaults (this is expected
+  //        since array_t->NumPyArray in the __init__ function is explicit and
+  //        leads to dangling reference)
+  // The following has to be used:
+  MyClassB(T &t, int p, pybind11::array_t< double >& b)
+  : a_(std::pow(t.a_,p)), b_(b) {}
+  double b() { return std::accumulate(b_.begin(),b_.end(),0.); }
+  int a_;
+  const Dune::Python::NumPyVector<double> b_;
 };
 """
 runCode="""
+#include <iostream>
 template <class T> int run(T &t) {
-  return t.a_ * t.b_;
+  return t.a_ * t.b();
 }
 """
 
+runVec="""
+#include <dune/python/common/numpyvector.hh>
+template <class T>
+void run(pybind11::array_t< T >& a)
+{
+  Dune::Python::NumPyVector< T > x( a );
+  for( size_t i=0; i<x.size(); ++i )
+    x[ i ] += i;
+}
+"""
+
+def test_numpyvector():
+    """
+    Test correct exchange of numpy arrays to C++ side.
+    """
+    from dune.generator.algorithm import run
+    x = array([0.]*100)
+    run("run",StringIO(runVec),x)
+    for i in range(len(x)):
+        assert x[i] == float(i)
+
 def test_class_export():
     from dune.generator.importclass import load
     from dune.generator.algorithm   import run
     from dune.generator import path
     from dune.typeregistry import generateTypeName
+    a = 2.
+    x = array([2.]*10)
     cls = load("MyClassA",StringIO(classACode),10,20)
     assert run("run",StringIO(runCode),cls) == 10*20
     clsName,includes = generateTypeName("MyClassB",cls)
-    cls = load(clsName,StringIO(classBCode),cls,2)
-    assert run("run",StringIO(runCode),cls) == 10**2*20**2
+    cls = load(clsName,StringIO(classBCode),cls,2,x)
+    assert run("run",StringIO(runCode),cls) == 10**2*10*a
+    x[:] = array([3.]*10)[:]
+    assert run("run",StringIO(runCode),cls) == 10**2*10*3
+    # the following does not work
+    x = array([4.]*10)
+    # the 'B' class still keeps the old vector 'x' alive
+    assert run("run",StringIO(runCode),cls) == 10**2*10*3
+
 if __name__ == "__main__":
     from dune.common.module import get_dune_py_dir
     _ = get_dune_py_dir()
     test_class_export()
+    test_numpyvector()