${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)
* 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
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;
* 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;
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
%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
%End
%RaiseCode
- printf("lALLA\n");
const char *detail = sipExceptionRef.what();
SIP_BLOCK_THREADS
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))