diff --git a/c++/cpp2py/CMakeLists.txt b/c++/cpp2py/CMakeLists.txt index 4012ad3..662e03e 100644 --- a/c++/cpp2py/CMakeLists.txt +++ b/c++/cpp2py/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(cpp2py signal_handler.cpp exceptions.cpp) +add_library(cpp2py signal_handler.cpp exceptions.cpp numpy_proxy.cpp) add_library(cpp2py::cpp2py ALIAS cpp2py) target_compile_options(cpp2py PRIVATE -std=c++14 -fPIC) diff --git a/c++/cpp2py/converters/vector.hpp b/c++/cpp2py/converters/vector.hpp index 8b77265..7a6f7c3 100644 --- a/c++/cpp2py/converters/vector.hpp +++ b/c++/cpp2py/converters/vector.hpp @@ -1,47 +1,123 @@ #pragma once -//#include -//#include +#include +#include +#include + +#include "../macros.hpp" +#include "../numpy_proxy.hpp" namespace cpp2py { + template static void delete_pycapsule(PyObject *capsule) { + auto *ptr = static_cast *>(PyCapsule_GetPointer(capsule, "guard")); + delete ptr; + } + + // Convert vector to numpy_proxy, WARNING: Deep Copy + template numpy_proxy make_numpy_proxy_from_vector(std::vector const &v) { + + auto *data_ptr = new std::unique_ptr{new T[v.size()]}; + std::copy(begin(v), end(v), data_ptr->get()); + auto capsule = PyCapsule_New(data_ptr, "guard", &delete_pycapsule); + + return {1, // rank + npy_type>, + (void *)data_ptr->get(), + std::is_const_v, + v_t{static_cast(v.size())}, // extents + v_t{sizeof(T)}, // strides + capsule}; + } + + // Make a new vector from numpy view + template std::vector make_vector_from_numpy_proxy(numpy_proxy const &p) { + EXPECTS(p.extents.size() == 1); + EXPECTS(p.strides == v_t{sizeof(T)}); + + T *data = static_cast(p.data); + long size = p.extents[0]; + + std::vector v(size); + std::copy(data, data + size, begin(v)); + return v; + } + + // -------------------------------------- + template struct py_converter> { - - // -------------------------------------- - - static PyObject *c2py(std::vector const &v) { - PyObject *list = PyList_New(0); - for (auto const &x : v) { - pyref y = py_converter::c2py(x); - if (y.is_null() or (PyList_Append(list, y) == -1)) { - Py_DECREF(list); - return NULL; - } // error + + static PyObject *c2py(std::vector const &v) { + + if constexpr (has_npy_type) { + return make_numpy_proxy_from_vector(v).to_python(); + } else { // Convert to Python List + PyObject *list = PyList_New(0); + for (auto const &x : v) { + pyref y = py_converter::c2py(x); + if (y.is_null() or (PyList_Append(list, y) == -1)) { + Py_DECREF(list); + return NULL; + } // error + } + return list; } - return list; } - // -------------------------------------- + // -------------------------------------- + + static bool is_convertible(PyObject *ob, bool raise_exception) { + _import_array(); - static bool is_convertible(PyObject *ob, bool raise_exception) { - if (!PySequence_Check(ob)) goto _false; - { - pyref seq = PySequence_Fast(ob, "expected a sequence"); - int len = PySequence_Size(ob); - for (int i = 0; i < len; i++) - if (!py_converter::is_convertible(PySequence_Fast_GET_ITEM((PyObject *)seq, i), raise_exception)) goto _false; //borrowed ref - return true; + // Special case: 1-d ndarray of builtin type + if (PyArray_Check(ob)) { + PyArrayObject *arr = (PyArrayObject *)(ob); +#ifdef PYTHON_NUMPY_VERSION_LT_17 + int rank = arr->nd; +#else + int rank = PyArray_NDIM(arr); +#endif + if (PyArray_TYPE(arr) == npy_type and rank == 1) return true; + } + + if (!PySequence_Check(ob)) { + if (raise_exception) { PyErr_SetString(PyExc_TypeError, "Cannot convert a non-sequence to std::vector"); } + return false; } - _false: - if (raise_exception) { PyErr_SetString(PyExc_TypeError, "Cannot convert to std::vector"); } - return false; - } - // -------------------------------------- - - static std::vector py2c(PyObject *ob) { pyref seq = PySequence_Fast(ob, "expected a sequence"); + int len = PySequence_Size(ob); + for (int i = 0; i < len; i++) { + if (!py_converter::is_convertible(PySequence_Fast_GET_ITEM((PyObject *)seq, i), false)) { // borrowed ref + if (raise_exception) { + auto err = std::string{"Cannot convert sequence to std::vector due to element at position "} + std::to_string(i); + PyErr_SetString(PyExc_TypeError, err.c_str()); + } + return false; + } + } + return true; + } + + // -------------------------------------- + + static std::vector py2c(PyObject *ob) { + _import_array(); + + // Special case: 1-d ndarray of builtin type + if (PyArray_Check(ob)) { + PyArrayObject *arr = (PyArrayObject *)(ob); +#ifdef PYTHON_NUMPY_VERSION_LT_17 + int rank = arr->nd; +#else + int rank = PyArray_NDIM(arr); +#endif + if (rank == 1) return make_vector_from_numpy_proxy(make_numpy_proxy(ob)); + } + + ASSERT(PySequence_Check(ob)); std::vector res; - int len = PySequence_Size(ob); + pyref seq = PySequence_Fast(ob, "expected a sequence"); + int len = PySequence_Size(ob); for (int i = 0; i < len; i++) res.push_back(py_converter::py2c(PySequence_Fast_GET_ITEM((PyObject *)seq, i))); //borrowed ref return res; } diff --git a/c++/cpp2py/macros.hpp b/c++/cpp2py/macros.hpp new file mode 100644 index 0000000..24188ba --- /dev/null +++ b/c++/cpp2py/macros.hpp @@ -0,0 +1,71 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011 by O. Parcollet + * + * TRIQS is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. + * + * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along with + * TRIQS. If not, see . + * + ******************************************************************************/ +#pragma once + +#include + +#define AS_STRING(...) AS_STRING2(__VA_ARGS__) +#define AS_STRING2(...) #__VA_ARGS__ + +#ifdef __clang__ +#define REQUIRES(...) __attribute__((enable_if(__VA_ARGS__, AS_STRING(__VA_ARGS__)))) +#elif __GNUC__ +#define REQUIRES(...) requires(__VA_ARGS__) +#endif + +#define PRINT(X) std::cerr << AS_STRING(X) << " = " << X << " at " << __FILE__ << ":" << __LINE__ << '\n' + +#define FORCEINLINE __inline__ __attribute__((always_inline)) + +#define EXPECTS(X) \ + if (!(X)) { \ + std::cerr << "Precondition " << AS_STRING(X) << " violated at " << __FILE__ << ":" << __LINE__ << "\n"; \ + std::terminate(); \ + } +#define ASSERT(X) \ + if (!(X)) { \ + std::cerr << "Assertion " << AS_STRING(X) << " violated at " << __FILE__ << ":" << __LINE__ << "\n"; \ + std::terminate(); \ + } +#define ENSURES(X) \ + if (!(X)) { \ + std::cerr << "Postcondition " << AS_STRING(X) << " violated at " << __FILE__ << ":" << __LINE__ << "\n"; \ + std::terminate(); \ + } + +#define EXPECTS_WITH_MESSAGE(X, ...) \ + if (!(X)) { \ + std::cerr << "Precondition " << AS_STRING(X) << " violated at " << __FILE__ << ":" << __LINE__ << "\n"; \ + std::cerr << "Error message : " << __VA_ARGS__ << std::endl; \ + std::terminate(); \ + } +#define ASSERT_WITH_MESSAGE(X, ...) \ + if (!(X)) { \ + std::cerr << "Assertion " << AS_STRING(X) << " violated at " << __FILE__ << ":" << __LINE__ << "\n"; \ + std::cerr << "Error message : " << __VA_ARGS__ << std::endl; \ + std::terminate(); \ + } +#define ENSURES_WITH_MESSAGE(X, ...) \ + if (!(X)) { \ + std::cerr << "Postcondition " << AS_STRING(X) << " violated at " << __FILE__ << ":" << __LINE__ << "\n"; \ + std::cerr << "Error message : " << __VA_ARGS__ << std::endl; \ + std::terminate(); \ + } diff --git a/c++/cpp2py/numpy_proxy.cpp b/c++/cpp2py/numpy_proxy.cpp new file mode 100644 index 0000000..5d59367 --- /dev/null +++ b/c++/cpp2py/numpy_proxy.cpp @@ -0,0 +1,118 @@ +#include "numpy_proxy.hpp" + +namespace cpp2py { + + // Make a new view_info + PyObject *numpy_proxy::to_python() { + + // Apparently we can not get rid of this + _import_array(); + +#ifdef PYTHON_NUMPY_VERSION_LT_17 + int flags = NPY_BEHAVED & ~NPY_OWNDATA; +#else + int flags = NPY_ARRAY_BEHAVED & ~NPY_ARRAY_OWNDATA; +#endif + // make the array read only + if (is_const) flags &= ~NPY_ARRAY_WRITEABLE; + PyObject *result = + PyArray_NewFromDescr(&PyArray_Type, PyArray_DescrFromType(element_type), rank, extents.data(), strides.data(), data, flags, NULL); + if (not result) return nullptr; // the Python error is set + + if (!PyArray_Check(result)) { + PyErr_SetString(PyExc_RuntimeError, "The python object is not a numpy array"); + return nullptr; + } + + PyArrayObject *arr = (PyArrayObject *)(result); +#ifdef PYTHON_NUMPY_VERSION_LT_17 + arr->base = base; + assert(arr->flags == (arr->flags & ~NPY_OWNDATA)); +#else + int r = PyArray_SetBaseObject(arr, base); + //EXPECTS(r == 0); + //EXPECTS(PyArray_FLAGS(arr) == (PyArray_FLAGS(arr) & ~NPY_ARRAY_OWNDATA)); +#endif + base = nullptr; // ref is stolen by the new object + + return result; + } + + // ---------------------------------------------------------- + + // Extract a view_info from python + numpy_proxy make_numpy_proxy(PyObject *obj) { + + // Apparently we can not get rid of this + _import_array(); + + if (obj == NULL) return {}; + if (not PyArray_Check(obj)) return {}; + + numpy_proxy result; + + // extract strides and lengths + PyArrayObject *arr = (PyArrayObject *)(obj); + +#ifdef PYTHON_NUMPY_VERSION_LT_17 + result.rank = arr->nd; +#else + result.rank = PyArray_NDIM(arr); +#endif + + result.element_type = PyArray_TYPE(arr); + result.extents.resize(result.rank); + result.strides.resize(result.rank); + result.data = PyArray_DATA(arr); + // base is ignored, stays at nullptr + +#ifdef PYTHON_NUMPY_VERSION_LT_17 + for (long i = 0; i < result.rank; ++i) { + result.extents[i] = size_t(arr->dimensions[i]); + result.strides[i] = std::ptrdiff_t(arr->strides[i]); + } +#else + for (size_t i = 0; i < result.rank; ++i) { + result.extents[i] = size_t(PyArray_DIMS(arr)[i]); + result.strides[i] = std::ptrdiff_t(PyArray_STRIDES(arr)[i]); + } +#endif + + //PRINT(result.rank); + //PRINT(result.element_type); + //PRINT(result.data); + + return result; + } + + // ---------------------------------------------------------- + + PyObject *make_numpy_copy(PyObject *obj, int rank, long element_type) { + + if (obj == nullptr) return nullptr; + + // From obj, we ask the numpy library to make a numpy, and of the correct type. + // This handles automatically the cases where : + // - we have list, or list of list/tuple + // - the numpy type is not the one we want. + // - adjust the dimension if needed + // If obj is an array : + // - if Order is same, don't change it + // - else impose it (may provoque a copy). + // if obj is not array : + // - Order = FortranOrder or SameOrder - > Fortran order otherwise C + + int flags = 0; //(ForceCast ? NPY_FORCECAST : 0) ;// do NOT force a copy | (make_copy ? NPY_ENSURECOPY : 0); + //if (!(PyArray_Check(obj) )) + //flags |= ( IndexMapType::traversal_order == indexmaps::mem_layout::c_order(rank) ? NPY_C_CONTIGUOUS : NPY_F_CONTIGUOUS); //impose mem order +#ifdef PYTHON_NUMPY_VERSION_LT_17 + flags |= (NPY_C_CONTIGUOUS); //impose mem order + flags |= (NPY_ENSURECOPY); +#else + flags |= (NPY_ARRAY_C_CONTIGUOUS); // impose mem order + flags |= (NPY_ARRAY_ENSURECOPY); +#endif + return PyArray_FromAny(obj, PyArray_DescrFromType(element_type), rank, rank, flags, NULL); // new ref + } + +} // namespace nda::python diff --git a/c++/cpp2py/numpy_proxy.hpp b/c++/cpp2py/numpy_proxy.hpp new file mode 100644 index 0000000..ed489dc --- /dev/null +++ b/c++/cpp2py/numpy_proxy.hpp @@ -0,0 +1,58 @@ +#pragma once +#include +#include + +#include +#include + +namespace cpp2py { + + using v_t = std::vector; + + // the basic information for a numpy array + struct numpy_proxy { + int rank = 0; + long element_type = 0; + void *data = nullptr; + bool is_const = false; + v_t extents, strides; + PyObject *base = nullptr; // The ref. counting guard typically + + // Returns a new ref (or NULL if failure) with a new numpy. + // If failure, return null with the Python exception set + PyObject *to_python(); + }; + + // From a numpy, extract the info. Better than a constructor, I want to use the aggregate constructor of the struct also. + numpy_proxy make_numpy_proxy(PyObject *); + + // Make a copy in Python with the given rank and element_type + // If failure, return null with the Python exception set + PyObject *make_numpy_copy(PyObject *obj, int rank, long elements_type); + + // + template static constexpr long npy_type = -1; + template static constexpr bool has_npy_type = (npy_type >= 0); + +#define NPY_CONVERT(C, P) template <> static constexpr long npy_type = P; + NPY_CONVERT(bool, NPY_BOOL) + NPY_CONVERT(char, NPY_STRING) + NPY_CONVERT(signed char, NPY_BYTE) + NPY_CONVERT(unsigned char, NPY_UBYTE) + NPY_CONVERT(short, NPY_SHORT) + NPY_CONVERT(unsigned short, NPY_USHORT) + NPY_CONVERT(int, NPY_INT) + NPY_CONVERT(unsigned int, NPY_UINT) + NPY_CONVERT(long, NPY_LONG) + NPY_CONVERT(unsigned long, NPY_ULONG) + NPY_CONVERT(long long, NPY_LONGLONG) + NPY_CONVERT(unsigned long long, NPY_ULONGLONG) + NPY_CONVERT(float, NPY_FLOAT) + NPY_CONVERT(double, NPY_DOUBLE) + NPY_CONVERT(long double, NPY_LONGDOUBLE) + NPY_CONVERT(std::complex, NPY_CFLOAT) + NPY_CONVERT(std::complex, NPY_CDOUBLE) + NPY_CONVERT(std::complex, NPY_CLONGDOUBLE) +#undef NPY_CONVERT + +} // namespace cpp2py