everything should be working but it does not

This commit is contained in:
Daniel Knüttel 2019-07-12 19:53:13 +02:00
parent 757efb2e5b
commit dae8093bca
4 changed files with 79 additions and 49 deletions

View File

@ -13,6 +13,7 @@ interaction_force_function
int i; int i;
for(i = 0; i < 7; i++) for(i = 0; i < 7; i++)
{ {
printf("%x\n", coefficients);
result += coefficients[i] * powf(r, i); result += coefficients[i] * powf(r, i);
} }
result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9])); result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9]));
@ -27,7 +28,7 @@ interaction_ufunc_force
( char ** args ( char ** args
, npy_intp * dimensions , npy_intp * dimensions
, npy_intp * steps , npy_intp * steps
, void * data) , void ** data)
{ {
char * in = args[0] char * in = args[0]
, * out = args[1]; , * out = args[1];
@ -37,7 +38,7 @@ interaction_ufunc_force
npy_intp i; npy_intp i;
float * coefficients = (float *) data; float * coefficients = (float *) data[0];
for(i = 0; i < n; i++) 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[] = static char interaction_types[] =
{ NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT}; { NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT};
static char force_types[] = static char force_types[] =
{ NPY_FLOAT, NPY_FLOAT}; { 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 static int
interaction_UFuncWrapper_init interaction_UFuncWrapper_init
@ -161,21 +152,31 @@ interaction_UFuncWrapper_init
int i; int i;
PyObject * this_coefficient; 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; 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; return -1;
} }
// copy the coefficients. // copy the coefficients.
for(i = 0; i < 19; i++) 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, // XXX: PyFloat_AsDouble might call python code,
// so make sure that nothing bad can happen. // so make sure that nothing bad can happen.
Py_INCREF(this_coefficient); Py_INCREF(this_coefficient);
@ -193,10 +194,10 @@ interaction_UFuncWrapper_init
{ {
self->ufunc = PyUFunc_FromFuncAndData( self->ufunc = PyUFunc_FromFuncAndData(
force_funcs force_funcs
, self->coefficients , (void **) &(self->coefficients)
, force_types , force_types
, 1 , 1
, 2 , 1
, 1 , 1
, PyUFunc_None , PyUFunc_None
, "force_function" , "force_function"
@ -208,10 +209,10 @@ interaction_UFuncWrapper_init
{ {
self->ufunc = PyUFunc_FromFuncAndData( self->ufunc = PyUFunc_FromFuncAndData(
interaction_funcs interaction_funcs
, interaction_data , (void **) &(self->coefficients)
, interaction_types , interaction_types
, 1 , 1
, 5 , 4
, 2 , 2
, PyUFunc_None , PyUFunc_None
, "interaction2D" , "interaction2D"
@ -222,10 +223,12 @@ interaction_UFuncWrapper_init
default: default:
{ {
PyErr_SetString(PyExc_ValueError, "unknown ufunc type, must be 0 or 1"); 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); 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[] = { static PyMethodDef InteractionMethods[] = {
{NULL, NULL, 0, NULL} {NULL, NULL, 0, NULL}
}; };
@ -258,10 +283,12 @@ static struct PyModuleDef moduledef = {
PyMODINIT_FUNC PyMODINIT_FUNC
PyInit_interaction(void) PyInit_interaction(void)
{ {
PyObject * module PyObject * module;
, * ufunc_interaction
, * ufunc_force if(PyType_Ready(&interaction_UFuncWrapperType) < 0)
, * dct; {
return NULL;
}
module = PyModule_Create(&moduledef); module = PyModule_Create(&moduledef);
if(!module) if(!module)
@ -271,14 +298,8 @@ PyInit_interaction(void)
import_array(); import_array();
import_ufunc(); import_ufunc();
ufunc_interaction = PyUFunc_FromFuncAndDataAndSignature( Py_INCREF(&interaction_UFuncWrapperType);
ufunc_force = PyUFunc_FromFuncAndDataAndSignature( PyModule_AddObject(module, "UFuncWrapper", (PyObject *) &interaction_UFuncWrapperType);
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);
return module; return module;
} }

View File

@ -2,8 +2,12 @@
#define interaction_h #define interaction_h
#include <Python.h> #include <Python.h>
#include "structmember.h"
// XXX
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/ndarraytypes.h> #include <numpy/ndarraytypes.h>
#include <numpy/ufuncobject.h> #include <numpy/ufuncobject.h>
#include <stddef.h>
/* /*
* This is a quite generic force function mapping a * This is a quite generic force function mapping a
@ -33,6 +37,6 @@ interaction_ufunc_force
( char ** args ( char ** args
, npy_intp * dimensions , npy_intp * dimensions
, npy_intp * steps , npy_intp * steps
, void * data) , void ** data);
#endif #endif

View File

@ -22,7 +22,7 @@ def configuration(parent_package='', top_path=None):
config = Configuration('brown' config = Configuration('brown'
, parent_package , parent_package
, top_path) , top_path)
config.add_extension('interaction', ['c/interaction/interaction.c']) config.add_extension('interaction', ['c/interaction/interaction.c'], extra_compile_args=["-g"])
return config return config
if __name__ == "__main__": if __name__ == "__main__":

View File

@ -1,4 +1,4 @@
from brown.interaction import interaction2D,force_function from brown.interaction import UFuncWrapper
import numpy as np import numpy as np
from collections import deque from collections import deque
from copy import copy 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) 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(c)
print(len(c))
#x_coords = np.array([0, 1], dtype=np.float16) #x_coords = np.array([0, 1], dtype=np.float16)
#y_coords = np.array([0, 0], dtype=np.float16) #y_coords = np.array([0, 0], dtype=np.float16)
@ -30,11 +31,15 @@ print(c)
#particle_two_evolution = time_evolution[:, 1] #particle_two_evolution = time_evolution[:, 1]
#plt.subplot(2, 1, 1) #plt.subplot(2, 1, 1)
force_function = UFuncWrapper(0, c)
r = np.arange(0, 0.5, 0.01, dtype=np.float16) r = np.arange(0, 0.5, 0.01, dtype=np.float16)
plt.title("Force") plt.title("Force")
plt.xlabel("particle distance") plt.xlabel("particle distance")
plt.ylabel("scalar force") 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.subplot(2, 1, 2)
#plt.title("Particle x-positions over time") #plt.title("Particle x-positions over time")