Move to NumPy in Python API
authorMaxim Koltsov <maks@omrb.pnpi.spb.ru>
Sat, 22 Nov 2014 16:35:58 +0000 (19:35 +0300)
committerMaxim Koltsov <maks@omrb.pnpi.spb.ru>
Sat, 22 Nov 2014 16:35:58 +0000 (19:35 +0300)
Drop custom implementation of python array interface and use numpy C API
for this task. This makes code cleaner and easier to use on python side.

src/python/CMakeLists.txt
src/python/include/numpy_conv.h [moved from src/python/include/vec.h with 53% similarity]
src/python/sip/definitions.sip
src/python/sip/trajectoryanalysis/analysismodule.sip
src/python/test.py

index 1aec631f198420afd1d040f0f983af6b801c7b16..70f980001782ebe90d3d095b8acec1d0d7048216 100644 (file)
@@ -44,19 +44,26 @@ include_directories(
     ${SIP_INCLUDE_DIR}
     ${GROMACS_INCLUDE_DIRS}
     ${CMAKE_SOURCE_DIR}/src/python/include
+    ${PYTHON_SITE_PACKAGES_INSTALL_DIR}/numpy/core/include
 )
 
+add_definitions(-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION)
+
 set(SIP_INCLUDES ${CMAKE_BINARY_DIR} sip)
 set(SIP_EXTRA_OPTIONS ${SIP_EXTRA_OPTIONS} -e -o)
 
+# This is needed for NumPy, which does not like code in many files
+set(SIP_CONCAT_PARTS 1)
+
 file(GLOB common_files_sip sip/*.sip)
 
 file(GLOB options_files_sip sip/options/*.sip)
 set(SIP_EXTRA_FILES_DEPEND ${options_files_sip} ${common_files_sip})
-add_sip_python_module(gromacs.Options sip/options/Options.sip ${GROMACS_LIBRARIES})
+add_sip_python_module(gromacs.Options sip/options/Options.sip libgromacs)
 
 file(GLOB trajectoryanalysis_files_sip sip/trajectoryanalysis/*.sip)
 set(SIP_EXTRA_FILES_DEPEND ${trajectoryanalysis_files_sip} ${common_files_sip})
-add_sip_python_module(gromacs.TrajectoryAnalysis sip/trajectoryanalysis/TrajectoryAnalysis.sip ${GROMACS_LIBRARIES})
+add_sip_python_module(gromacs.TrajectoryAnalysis
+    sip/trajectoryanalysis/TrajectoryAnalysis.sip libgromacs)
 
 python_install(__init__.py ${PYTHON_SITE_PACKAGES_INSTALL_DIR}/gromacs)
similarity index 53%
rename from src/python/include/vec.h
rename to src/python/include/numpy_conv.h
index 087f27f3defa9311fa9b4407a0c6fefe404f4265..70d9dc519091b94e1c95b6c187afe74cf0594497 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
-#ifndef VEC_H
-#define VEC_H
+#ifndef NUMPY_CONV_H
+#define NUMPY_CONV_H
 
-#include <cstdlib>
-#include <stdexcept>
-#include <gromacs/math/vectypes.h>
+#include <Python.h>
+#include <numpy/ndarrayobject.h>
 
-template<typename T> class RealVector {
-public:
-    T x, y, z;
-    RealVector(T v[3]) : x(v[0]), y(v[1]), z(v[2]) {};
-    T operator[] (size_t i) {
-        switch (i) {
-        case 0:
-            return x;
-        case 1:
-            return y;
-        case 2:
-            return z;
-        default:
-            throw std::out_of_range("Inexistent component requested");
-        };
-    };
-};
+PyObject* array2dToNumpy(int dim1, int dim2, void *data) {
+    npy_intp dims[] = {dim1, dim2};
+    #ifdef GMX_DOUBLE
+        return PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, (double*) data);
+    #else
+        return PyArray_SimpleNewFromData(2, dims, NPY_FLOAT, (float*) data);
+    #endif
+}
 
-template<typename T, typename Adaptor, typename Inner> class Array {
-private:
-    Inner *data;
-    size_t length;
-public:
-    Array(Inner *data, size_t length) : data(data), length(length) {};
-    Adaptor operator[] (size_t i) {
-        if (i >= length)
-            throw std::out_of_range("Array index out of range");
-
-        return Adaptor(data[i]);
-    };
-    size_t len() { return length; };
-};
-
-class Matrix {
-private:
-    real data[3][3];
-public:
-    Matrix(real data[3][3]) {
-        std::copy((real*) data, ((real*) data) + 9, (real*) this->data);
-    };
-    real get(size_t i, size_t j) {
-        if (i > 2 || j > 2)
-            throw std::out_of_range("Matrix index out of range");
-        return data[i][j];
-    }
-    RealVector<real> get(size_t i) {
-        if (i > 2)
-            throw std::out_of_range("Matrix index out of range");
-        return RealVector<real>(data[i]);
-    }
-};
-
-typedef RealVector<real> py_rvec;
-typedef Array<real, py_rvec, rvec> rvecs;
+PyObject* array1dToNumpy(int dim, void *data) {
+    npy_intp n_dim = dim;
+    #ifdef GMX_DOUBLE
+        return PyArray_SimpleNewFromData(1, &n_dim, NPY_DOUBLE, (double*) data);
+    #else
+        return PyArray_SimpleNewFromData(1, &n_dim, NPY_FLOAT, (float*) data);
+    #endif
+}
 #endif
index b2c5b5d183707979b55a2c2c40079d9e94cbabca..35e9631c6d0ef14191d81bd57805b0700dad222c 100644 (file)
     SIP_UNBLOCK_THREADS
 %End
 };
-
-template<T> class RealVector /NoDefaultCtors/ {
-%TypeHeaderCode
-#include "vec.h"
-%End
-
-public:
-    T x;
-    T y;
-    T z;
-    T operator[] (unsigned int i) throw (std::out_of_range);
-    PyObject* __str__();
-    %MethodCode
-        size_t size = PyOS_snprintf(NULL, 0, "Vector [%f, %f, %f]", sipCpp->x, sipCpp->y, sipCpp->z);
-        char str[size + 1];
-        PyOS_snprintf(str, size + 1, "Vector [%f, %f, %f]", sipCpp->x, sipCpp->y, sipCpp->z);
-
-        return PyUnicode_FromString(str);
-    %End
-};
-
-template<T, Adaptor, Inner> class Array /NoDefaultCtors/ {
-%TypeHeaderCode
-#include "vec.h"
-%End
-
-public:
-    Adaptor operator[] (unsigned int i) throw (std::out_of_range);
-    SIP_SSIZE_T __len__();
-    %MethodCode
-        sipRes = sipCpp->len();
-    %End
-};
-
-class Matrix /NoDefaultCtors/ {
-%TypeHeaderCode
-#include "vec.h"
-%End
-
-public:
-    py_rvec __getitem__(int) throw (std::out_of_range);
-    %MethodCode
-        try {
-            sipRes = new py_rvec(sipCpp->get(a0));
-        } catch (const std::out_of_range &e) {
-            sipIsErr = 1;
-            SIP_BLOCK_THREADS
-            PyErr_SetString(PyExc_IndexError, e.what());
-            SIP_UNBLOCK_THREADS
-        }
-    %End
-    double __getitem__(SIP_PYTUPLE) throw (std::out_of_range);
-    %MethodCode
-        size_t i, j;
-        if (PyArg_ParseTuple(a0, "nn", &i, &j))
-            try {
-                sipRes = sipCpp->get(i, j);
-            } catch (const std::out_of_range &e) {
-                sipIsErr = 1;
-                SIP_BLOCK_THREADS
-                PyErr_SetString(PyExc_IndexError, e.what());
-                SIP_UNBLOCK_THREADS
-            }
-        else
-            sipIsErr = 1;
-    %End
-};
-
-typedef double* rvec;
-typedef double* dvec;
-typedef double** matrix;
-typedef double** tensor;
-typedef int* ivec;
-typedef int** imatrix;
-
-typedef RealVector<double> py_rvec;
-typedef Array<double, py_rvec, rvec> rvecs;
index c482ba316676c8bc8645b0ca5b4a893d8b7fe454..c5ef94e94896044433f1a2e0cd1e97b8c2f46965 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
+%PostInitialisationCode
+    import_array();
+%End
+
 struct t_trxframe {
 
 %TypeHeaderCode
 #include <gromacs/fileio/trx.h>
-#include "vec.h"
+#include "numpy_conv.h"
 %End
 
     int flags;
@@ -64,40 +68,36 @@ struct t_trxframe {
     bool bPrec;
     double prec;
     bool bX;
-    rvec *x {
+    SIP_PYOBJECT x {
     %GetCode
-        rvecs *arr = new rvecs(sipCpp->x, sipCpp->natoms);
-        sipPy = sipConvertFromNewType(arr, sipType_rvecs, NULL);
+        sipPy = array2dToNumpy(sipCpp->natoms, 3, sipCpp->x);
     %End
     %SetCode
 
     %End
     };
     bool bV;
-    rvec *v {
+    SIP_PYOBJECT v {
     %GetCode
-        rvecs *arr = new rvecs(sipCpp->v, sipCpp->natoms);
-        sipPy = sipConvertFromNewType(arr, sipType_rvecs, NULL);
+        sipPy = array2dToNumpy(sipCpp->natoms, 3, sipCpp->v);
     %End
     %SetCode
 
     %End
     };
     bool bF;
-    rvec *f {
+    SIP_PYOBJECT f {
     %GetCode
-        rvecs *arr = new rvecs(sipCpp->f, sipCpp->natoms);
-        sipPy = sipConvertFromNewType(arr, sipType_rvecs, NULL);
+        sipPy = array2dToNumpy(sipCpp->natoms, 3, sipCpp->f);
     %End
     %SetCode
 
     %End
     };
     bool bBox;
-    matrix box {
+    SIP_PYOBJECT box {
     %GetCode
-        Matrix *mat = new Matrix(sipCpp->box);
-        sipPy = sipConvertFromNewType(mat, sipType_Matrix, NULL);
+        sipPy = array2dToNumpy(3, 3, sipCpp->box);
     %End
     %SetCode
 
@@ -113,43 +113,39 @@ struct t_pbc {
 
 %TypeHeaderCode
 #include <gromacs/pbcutil/pbc.h>
-#include "vec.h"
+#include "numpy_conv.h"
 %End
     int ePBC;
     int ndim_ePBC;
     int ePBCDX;
     int dim;
-    matrix box {
+    SIP_PYOBJECT box {
     %GetCode
-        Matrix *mat = new Matrix(sipCpp->box);
-        sipPy = sipConvertFromNewType(mat, sipType_Matrix, NULL);
+        sipPy = array2dToNumpy(3, 3, sipCpp->box);
     %End
     %SetCode
 
     %End
     };
-    rvec fbox_diag {
+    SIP_PYOBJECT fbox_diag {
     %GetCode
-        py_rvec *vec = new py_rvec(sipCpp->fbox_diag);
-        sipPy = sipConvertFromNewType(vec, sipType_py_rvec, NULL);
+        sipPy = array1dToNumpy(3, sipCpp->fbox_diag);
     %End
     %SetCode
 
     %End
     };
-    rvec hbox_diag {
+    SIP_PYOBJECT hbox_diag {
     %GetCode
-        py_rvec *vec = new py_rvec(sipCpp->fbox_diag);
-        sipPy = sipConvertFromNewType(vec, sipType_py_rvec, NULL);
+        sipPy = array1dToNumpy(3, sipCpp->hbox_diag);
     %End
     %SetCode
 
     %End
     };
-    rvec mhbox_diag {
+    SIP_PYOBJECT mhbox_diag {
     %GetCode
-        py_rvec *vec = new py_rvec(sipCpp->mhbox_diag);
-        sipPy = sipConvertFromNewType(vec, sipType_py_rvec, NULL);
+        sipPy = array1dToNumpy(3, sipCpp->mhbox_diag);
     %End
     %SetCode
 
@@ -207,7 +203,6 @@ using namespace std;
 %End
 
 %RaiseCode
-    printf("lALLA\n");
     const char *detail = sipExceptionRef.what();
 
     SIP_BLOCK_THREADS
index 5c01936423f86ae9f99eaf3117adcfb8486eb1b2..e2b77c7626229304294fda41f210ef77004d69d7 100644 (file)
@@ -16,8 +16,8 @@ class M(TrajectoryAnalysis.TrajectoryAnalysisModule):
 
     def analyzeFrame(self, frnr, frame, pbc, data):
         print('python: Analyzing frame {}, {} atoms'.format(frnr, frame.natoms))
-        #print(frame.box[0,0], frame.box[0,1], frame.box[0,2])
-        print(pbc.box[0], pbc.box[1], pbc.box[2])
+        print(frame.box[0], frame.box[1], frame.box[2])
+        print(frame.x[0], frame.v[0], frame.f[0])
 
     def finishAnalysis(self, nframes):
         print('python: Analyzed {} frames'.format(nframes))