diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index efb50a0..d9b78a0 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -13,6 +13,7 @@ interaction_force_function int i; for(i = 0; i < 7; i++) { + printf("%x\n", coefficients); result += coefficients[i] * powf(r, i); } result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9])); @@ -27,7 +28,7 @@ interaction_ufunc_force ( char ** args , npy_intp * dimensions , npy_intp * steps - , void * data) + , void ** data) { char * in = args[0] , * out = args[1]; @@ -37,7 +38,7 @@ interaction_ufunc_force npy_intp i; - float * coefficients = (float *) data; + float * coefficients = (float *) data[0]; for(i = 0; i < n; i++) { @@ -122,31 +123,21 @@ interaction_ufunc_float2D } } } - -typedef struct -{ - PyObject_HEAD - float coefficients[16]; - PyObject * ufunc; -} interaction_UFuncWrapper; - -static PyTypeObject interaction_UFuncWrapperType -{ - PyTypeObject_HEAD_INIT(NULL, 0) - .tp_name = "brown.interaction.UFuncWrapper", - .tp_doc = "A wrapper that wraps the ufuncs for interaction and force, storing the coefficients", - .tp_basicsize = sizeof(interaction_UFuncWrapper), - .tp_itemsize = 0, - .tp_flags = Py_TPFLAGS_DEFAULT, - .tp_new = PyType_GenericNew, - -}; - - static char interaction_types[] = { NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT}; static char force_types[] = { NPY_FLOAT, NPY_FLOAT}; +static PyUFuncGenericFunction force_funcs[1] = + { interaction_ufunc_force}; +static PyUFuncGenericFunction interaction_funcs[1] = + { interaction_ufunc_float2D}; + +typedef struct +{ + PyObject_HEAD + float coefficients[19]; + PyObject * ufunc; +} interaction_UFuncWrapper; static int interaction_UFuncWrapper_init @@ -161,21 +152,31 @@ interaction_UFuncWrapper_init int i; PyObject * this_coefficient; - if(!PyArgs_ParseTupleAndKeywords(args, kwds, "BO!", &type, &PyList_Type, &coefficients)) + char *kwords[] = { "type_", "coefficients", NULL}; + + if(!PyArg_ParseTupleAndKeywords(args, kwds, "bO", kwords, &type, &coefficients)) + { + return -1; + } + if(!PySequence_Check(coefficients)) { return -1; } - if(PyList_Size(coefficients) != 19) + if(PySequence_Size(coefficients) != 19) { - PyErr_SetString("coefficients must have length 19", PyExc_ValueError); + PyErr_SetString(PyExc_ValueError, "coefficients must have length 19"); return -1; } // copy the coefficients. for(i = 0; i < 19; i++) { - this_coefficient = PyList_GetItem(coefficients, i); + this_coefficient = PySequence_GetItem(coefficients, i); + if(!this_coefficient) + { + return -1; + } // XXX: PyFloat_AsDouble might call python code, // so make sure that nothing bad can happen. Py_INCREF(this_coefficient); @@ -193,10 +194,10 @@ interaction_UFuncWrapper_init { self->ufunc = PyUFunc_FromFuncAndData( force_funcs - , self->coefficients + , (void **) &(self->coefficients) , force_types , 1 - , 2 + , 1 , 1 , PyUFunc_None , "force_function" @@ -208,10 +209,10 @@ interaction_UFuncWrapper_init { self->ufunc = PyUFunc_FromFuncAndData( interaction_funcs - , interaction_data + , (void **) &(self->coefficients) , interaction_types , 1 - , 5 + , 4 , 2 , PyUFunc_None , "interaction2D" @@ -222,10 +223,12 @@ interaction_UFuncWrapper_init default: { PyErr_SetString(PyExc_ValueError, "unknown ufunc type, must be 0 or 1"); - return -1 + return -1; } - } + Py_INCREF(self->ufunc); + printf("%x\n", self->coefficients); + return 0; } @@ -238,6 +241,28 @@ interaction_UFuncWrapper_call return PyObject_Call(self->ufunc, args, kwargs); } +static PyMemberDef interaction_UFuncWrapper_members[] = +{ + {"ufunc", T_OBJECT_EX, offsetof(interaction_UFuncWrapper, ufunc), 0, "ufunc"}, + {NULL} +}; + +static PyTypeObject interaction_UFuncWrapperType = +{ + PyVarObject_HEAD_INIT(NULL, 0) + .tp_name = "brown.interaction.UFuncWrapper", + .tp_doc = "A wrapper that wraps the ufuncs for interaction and force, storing the coefficients", + .tp_basicsize = sizeof(interaction_UFuncWrapper), + .tp_itemsize = 0, + .tp_flags = Py_TPFLAGS_DEFAULT, + .tp_new = PyType_GenericNew, + .tp_init = (initproc) interaction_UFuncWrapper_init, + .tp_call = interaction_UFuncWrapper_call, + .tp_members = interaction_UFuncWrapper_members +}; + + + static PyMethodDef InteractionMethods[] = { {NULL, NULL, 0, NULL} }; @@ -258,10 +283,12 @@ static struct PyModuleDef moduledef = { PyMODINIT_FUNC PyInit_interaction(void) { - PyObject * module - , * ufunc_interaction - , * ufunc_force - , * dct; + PyObject * module; + + if(PyType_Ready(&interaction_UFuncWrapperType) < 0) + { + return NULL; + } module = PyModule_Create(&moduledef); if(!module) @@ -271,14 +298,8 @@ PyInit_interaction(void) import_array(); import_ufunc(); - ufunc_interaction = PyUFunc_FromFuncAndDataAndSignature( - ufunc_force = PyUFunc_FromFuncAndDataAndSignature( - - dct = PyModule_GetDict(module); - PyDict_SetItemString(dct, "interaction2D", ufunc_interaction); - PyDict_SetItemString(dct, "force_function", ufunc_force); - Py_DECREF(ufunc_interaction); - Py_DECREF(ufunc_force); + Py_INCREF(&interaction_UFuncWrapperType); + PyModule_AddObject(module, "UFuncWrapper", (PyObject *) &interaction_UFuncWrapperType); return module; } diff --git a/c/interaction/interaction.h b/c/interaction/interaction.h index e152077..a47bf49 100644 --- a/c/interaction/interaction.h +++ b/c/interaction/interaction.h @@ -2,8 +2,12 @@ #define interaction_h #include +#include "structmember.h" +// XXX +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include +#include /* * This is a quite generic force function mapping a @@ -33,6 +37,6 @@ interaction_ufunc_force ( char ** args , npy_intp * dimensions , npy_intp * steps - , void * data) + , void ** data); #endif diff --git a/setup.py b/setup.py index 621c2f0..ec09553 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ def configuration(parent_package='', top_path=None): config = Configuration('brown' , parent_package , top_path) - config.add_extension('interaction', ['c/interaction/interaction.c']) + config.add_extension('interaction', ['c/interaction/interaction.c'], extra_compile_args=["-g"]) return config if __name__ == "__main__": diff --git a/test.py b/test.py index 9bdce85..ceff3fb 100644 --- a/test.py +++ b/test.py @@ -1,4 +1,4 @@ -from brown.interaction import interaction2D,force_function +from brown.interaction import UFuncWrapper import numpy as np from collections import deque from copy import copy @@ -6,6 +6,7 @@ import matplotlib.pyplot as plt c = np.array([5, 10, 20, 0, 0, 0, 0, 1, -2, 0, -2, -0.1, 2, 0, 0, 0, 0, 0, 0], dtype=np.float16) print(c) +print(len(c)) #x_coords = np.array([0, 1], dtype=np.float16) #y_coords = np.array([0, 0], dtype=np.float16) @@ -30,11 +31,15 @@ print(c) #particle_two_evolution = time_evolution[:, 1] #plt.subplot(2, 1, 1) + +force_function = UFuncWrapper(0, c) r = np.arange(0, 0.5, 0.01, dtype=np.float16) plt.title("Force") plt.xlabel("particle distance") plt.ylabel("scalar force") -plt.plot(r, force_function(r, c)) +print(r) +print(force_function.ufunc(r)) +plt.plot(r, force_function.ufunc(r)) #plt.subplot(2, 1, 2) #plt.title("Particle x-positions over time")