From aacd1473aa4afbca4b7936e3478b745454e565ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 10:48:49 +0200 Subject: [PATCH 1/9] some work towards objectifying the UFuncs --- c/interaction/interaction.c | 76 +++++++++++++++++++++++-------------- c/interaction/interaction.h | 6 +++ 2 files changed, 54 insertions(+), 28 deletions(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index 54c55bb..c829654 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -30,24 +30,14 @@ interaction_ufunc_force , void * data) { char * in = args[0] - , * raw_coefficients = args[1] - , * out = args[2]; + , * out = args[1]; npy_intp n = dimensions[0]; npy_intp in_step = steps[0] - , raw_coefficients_steps = steps[1] - , out_step = steps[2]; + , out_step = steps[1]; npy_intp i; - // Copy the coefficients to a safe array. - // This is because NumPy arrays are not - // necessarily cast-safe. - float coefficients[19]; - for(i = 0; i < 19; i++) - { - coefficients[i] = *(float *)raw_coefficients; - raw_coefficients += raw_coefficients_steps; - } + float * coefficients = (float *) data; for(i = 0; i < n; i++) { @@ -73,27 +63,17 @@ interaction_ufunc_float2D , * y_old = args[1] , * p_x_old = args[2] , * p_y_old = args[3] - , * raw_coefficients = args[4] - , * p_x_new = args[5] - , * p_y_new = args[6]; + , * p_x_new = args[4] + , * p_y_new = args[5]; npy_intp x_old_steps = steps[0] , y_old_steps = steps[1] , p_x_old_steps = steps[2] , p_y_old_steps = steps[3] - , raw_coefficients_steps = steps[4] - , p_x_new_steps = steps[5] - , p_y_new_steps = steps[6]; + , p_x_new_steps = steps[4] + , p_y_new_steps = steps[5]; - // Copy the coefficients to a safe array. - // This is because NumPy arrays are not - // necessarily cast-safe. - float coefficients[19]; - for(i = 0; i < 19; i++) - { - coefficients[i] = *(float *)raw_coefficients; - raw_coefficients += raw_coefficients_steps; - } + float *coefficients = (float *) data; // Compute the new momenta: @@ -143,6 +123,46 @@ 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 coeficients", + .tp_basicsize = sizeof(interaction_UFuncWrapper), + .tp_itemsize = 0, + .tp_flags = Py_TPFLAGS_DEFAULT, + .tp_new = PyType_GenericNew, + +}; + +static int +interaction_UFuncWrapper_init + ( interaction_UFuncWrapper * self + , PyObject * args + , PyObject * kwds) +{ + char type; + PyObject * coefficients; + + if(!PyArgs_ParseTupleAndKeywords(args, kwds, "" // FIXME +} + + +static PyObject * +interaction_UFuncWrapper_call + (interaction_UFuncWrapper * self, + PyObject * data) +{ + return +} + static PyMethodDef InteractionMethods[] = { {NULL, NULL, 0, NULL} }; diff --git a/c/interaction/interaction.h b/c/interaction/interaction.h index b42833a..e152077 100644 --- a/c/interaction/interaction.h +++ b/c/interaction/interaction.h @@ -28,5 +28,11 @@ interaction_ufunc_float2D , npy_intp * steps , void * data); +static void +interaction_ufunc_force + ( char ** args + , npy_intp * dimensions + , npy_intp * steps + , void * data) #endif -- 2.45.1 From 9220314b681bdda24d61301e9c1b30369088d78a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 11:07:24 +0200 Subject: [PATCH 2/9] a little more work --- c/interaction/interaction.c | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index c829654..9048dee 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -148,10 +148,23 @@ interaction_UFuncWrapper_init , PyObject * args , PyObject * kwds) { + // 0: ufunc_force + // 1: interaction2D char type; PyObject * coefficients; - if(!PyArgs_ParseTupleAndKeywords(args, kwds, "" // FIXME + if(!PyArgs_ParseTupleAndKeywords(args, kwds, "BO!", &type, &PyList_Type, &coefficients)) + { + return -1; + } + + if(PyList_Size(coefficients) != 19) + { + PyErr_SetString("coefficients must have length 19", PyExc_ValueError); + return -1; + } + + } @@ -167,6 +180,7 @@ static PyMethodDef InteractionMethods[] = { {NULL, NULL, 0, NULL} }; +// FIXME PyUFuncGenericFunction interaction_funcs[] = { &interaction_ufunc_float2D}; PyUFuncGenericFunction force_funcs[] = @@ -206,7 +220,7 @@ PyInit_interaction(void) return NULL; } import_array(); - import_umath(); + import_ufunc(); ufunc_interaction = PyUFunc_FromFuncAndDataAndSignature( interaction_funcs -- 2.45.1 From 757efb2e5bbc95c5200821ac4ff2bf0595462f54 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 13:51:43 +0200 Subject: [PATCH 3/9] some further work --- c/interaction/interaction.c | 105 ++++++++++++++++++++++-------------- 1 file changed, 66 insertions(+), 39 deletions(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index 9048dee..efb50a0 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -134,7 +134,7 @@ 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 coeficients", + .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, @@ -142,6 +142,12 @@ static PyTypeObject interaction_UFuncWrapperType }; + +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 int interaction_UFuncWrapper_init ( interaction_UFuncWrapper * self @@ -152,6 +158,8 @@ interaction_UFuncWrapper_init // 1: interaction2D char type; PyObject * coefficients; + int i; + PyObject * this_coefficient; if(!PyArgs_ParseTupleAndKeywords(args, kwds, "BO!", &type, &PyList_Type, &coefficients)) { @@ -164,35 +172,76 @@ interaction_UFuncWrapper_init return -1; } + // copy the coefficients. + for(i = 0; i < 19; i++) + { + this_coefficient = PyList_GetItem(coefficients, i); + // XXX: PyFloat_AsDouble might call python code, + // so make sure that nothing bad can happen. + Py_INCREF(this_coefficient); + self->coefficients[i] = PyFloat_AsDouble(this_coefficient); + Py_DECREF(this_coefficient); + if(PyErr_Occurred()) + { + return -1; + } + } + switch(type) + { + case 0: + { + self->ufunc = PyUFunc_FromFuncAndData( + force_funcs + , self->coefficients + , force_types + , 1 + , 2 + , 1 + , PyUFunc_None + , "force_function" + , "computes the scalar force between two particles with given coefficients" + , 0); + break; + } + case 1: + { + self->ufunc = PyUFunc_FromFuncAndData( + interaction_funcs + , interaction_data + , interaction_types + , 1 + , 5 + , 2 + , PyUFunc_None + , "interaction2D" + , "Update the momenta according to the given coefficients and positions" + , 0); + break; + } + default: + { + PyErr_SetString(PyExc_ValueError, "unknown ufunc type, must be 0 or 1"); + return -1 + } + + } } static PyObject * interaction_UFuncWrapper_call - (interaction_UFuncWrapper * self, - PyObject * data) + (interaction_UFuncWrapper * self + , PyObject * args + , PyObject * kwargs) { - return + return PyObject_Call(self->ufunc, args, kwargs); } static PyMethodDef InteractionMethods[] = { {NULL, NULL, 0, NULL} }; -// FIXME -PyUFuncGenericFunction interaction_funcs[] = - { &interaction_ufunc_float2D}; -PyUFuncGenericFunction force_funcs[] = - { &interaction_ufunc_force}; - -static char interaction_types[] = - { NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT}; -static char force_types[] = - { NPY_FLOAT, NPY_FLOAT, NPY_FLOAT}; - -static void *interaction_data[] = {NULL}; -static void *force_data[] = {NULL}; static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT @@ -223,29 +272,7 @@ PyInit_interaction(void) import_ufunc(); ufunc_interaction = PyUFunc_FromFuncAndDataAndSignature( - interaction_funcs - , interaction_data - , interaction_types - , 1 - , 5 - , 2 - , PyUFunc_None - , "interaction2D" - , "Update the momenta according to the given coefficients and positions" - , 0 - , "(n),(n),(n),(n),(19)->(n),(n)"); ufunc_force = PyUFunc_FromFuncAndDataAndSignature( - force_funcs - , force_data - , force_types - , 1 - , 2 - , 1 - , PyUFunc_None - , "force_function" - , "computes the scalar force between two particles with given coefficients" - , 0 - , "(n),(19)->(n)"); dct = PyModule_GetDict(module); PyDict_SetItemString(dct, "interaction2D", ufunc_interaction); -- 2.45.1 From dae8093bcad7fd7be0f8e6057973dd42da419b8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 19:53:13 +0200 Subject: [PATCH 4/9] everything should be working but it does not --- c/interaction/interaction.c | 111 +++++++++++++++++++++--------------- c/interaction/interaction.h | 6 +- setup.py | 2 +- test.py | 9 ++- 4 files changed, 79 insertions(+), 49 deletions(-) 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") -- 2.45.1 From d3a227ac192fc5bacd18441a5c086122be4c749c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 20:04:40 +0200 Subject: [PATCH 5/9] added some printfs for SO --- c/interaction/interaction.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index d9b78a0..ddf5368 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -38,6 +38,8 @@ interaction_ufunc_force npy_intp i; + printf("DATA IS: %x\n", data); + printf("DATA[0] IS: %x\n", data[0]); float * coefficients = (float *) data[0]; for(i = 0; i < n; i++) @@ -227,7 +229,8 @@ interaction_UFuncWrapper_init } } Py_INCREF(self->ufunc); - printf("%x\n", self->coefficients); + printf("PASSED = %x\n", &(self->coefficients)); + printf("PASSED[0] = %x\n", self->coefficients); return 0; } -- 2.45.1 From c52fbe796ddcc6e5c519d5be8d9ce61d2735e621 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 20:24:02 +0200 Subject: [PATCH 6/9] everything is technically working --- c/interaction/interaction.c | 33 ++++++++-------- c/interaction/interaction.h | 2 +- test.py | 75 ++++++++++++++++++------------------- 3 files changed, 52 insertions(+), 58 deletions(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index ddf5368..70adbbc 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -13,7 +13,6 @@ 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])); @@ -28,7 +27,7 @@ interaction_ufunc_force ( char ** args , npy_intp * dimensions , npy_intp * steps - , void ** data) + , void * data) { char * in = args[0] , * out = args[1]; @@ -38,9 +37,7 @@ interaction_ufunc_force npy_intp i; - printf("DATA IS: %x\n", data); - printf("DATA[0] IS: %x\n", data[0]); - float * coefficients = (float *) data[0]; + float * coefficients = (float *) data; for(i = 0; i < n; i++) { @@ -139,6 +136,7 @@ typedef struct PyObject_HEAD float coefficients[19]; PyObject * ufunc; + void *data[1]; } interaction_UFuncWrapper; static int @@ -189,29 +187,30 @@ interaction_UFuncWrapper_init return -1; } } + self->data[0] = (void *)self->coefficients; switch(type) { case 0: { self->ufunc = PyUFunc_FromFuncAndData( - force_funcs - , (void **) &(self->coefficients) - , force_types - , 1 - , 1 - , 1 - , PyUFunc_None - , "force_function" - , "computes the scalar force between two particles with given coefficients" - , 0); + force_funcs // func + , self->data // data + , force_types //types + , 1 // ntypes + , 1 // nin + , 1 // nout + , PyUFunc_None // identity + , "force_function" // name + , "computes the scalar force between two particles with given coefficients" // doc + , 0); // unused break; } case 1: { self->ufunc = PyUFunc_FromFuncAndData( interaction_funcs - , (void **) &(self->coefficients) + , self->data , interaction_types , 1 , 4 @@ -229,8 +228,6 @@ interaction_UFuncWrapper_init } } Py_INCREF(self->ufunc); - printf("PASSED = %x\n", &(self->coefficients)); - printf("PASSED[0] = %x\n", self->coefficients); return 0; } diff --git a/c/interaction/interaction.h b/c/interaction/interaction.h index a47bf49..dbafc86 100644 --- a/c/interaction/interaction.h +++ b/c/interaction/interaction.h @@ -37,6 +37,6 @@ interaction_ufunc_force ( char ** args , npy_intp * dimensions , npy_intp * steps - , void ** data); + , void * data); #endif diff --git a/test.py b/test.py index ceff3fb..25b8cd5 100644 --- a/test.py +++ b/test.py @@ -4,48 +4,45 @@ from collections import deque from copy import copy 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) -# -#x_momenta = np.array([0, 0], dtype=np.float16) -#y_momenta = np.array([0, 0], dtype=np.float16) -# -#time = np.arange(0, 5, 1, dtype=np.float16) -# -#time_evolution = deque() -# -#for t in time: -# x_momenta, y_momenta = interaction2D(x_coords, y_coords, x_momenta, y_momenta, c) -# x_coords += x_momenta -# y_coords += y_momenta -# -# time_evolution.append(copy(x_coords)) -# -#time_evolution = np.array(time_evolution) -# -#particle_one_evolution = time_evolution[:, 0] -#particle_two_evolution = time_evolution[:, 1] - -#plt.subplot(2, 1, 1) - +c = np.array([5, 10, 20, 30, 0, 0, 0, 1, -20, 0, -2, -0.1, 2, 0, 0, 0, 0, 0, 0], dtype=np.float16) force_function = UFuncWrapper(0, c) -r = np.arange(0, 0.5, 0.01, dtype=np.float16) +interaction2D = UFuncWrapper(1, c) + +x_coords = np.array([0, 1], dtype=np.float16) +y_coords = np.array([0, 0], dtype=np.float16) + +x_momenta = np.array([0, 0], dtype=np.float16) +y_momenta = np.array([0, 0], dtype=np.float16) + +time = np.arange(0, 5, 1, dtype=np.float16) + +time_evolution = deque() + +for t in time: + x_momenta, y_momenta = interaction2D(x_coords, y_coords, x_momenta, y_momenta) + x_coords += x_momenta + y_coords += y_momenta + + time_evolution.append(copy(x_coords)) + +time_evolution = np.array(time_evolution) + +particle_one_evolution = time_evolution[:, 0] +particle_two_evolution = time_evolution[:, 1] + +plt.subplot(2, 1, 1) + +r = np.arange(0, 50, 0.1, dtype=np.float16) plt.title("Force") plt.xlabel("particle distance") plt.ylabel("scalar force") -print(r) -print(force_function.ufunc(r)) -plt.plot(r, force_function.ufunc(r)) +plt.plot(r, force_function(r)) -#plt.subplot(2, 1, 2) -#plt.title("Particle x-positions over time") -#plt.xlabel("time") -#plt.ylabel("particle x-position") -#h0, = plt.plot(time, particle_one_evolution, label="particle 1") -#h1, = plt.plot(time, particle_two_evolution, label="particle 2") -#plt.legend(handles=[h0, h1]) +plt.subplot(2, 1, 2) +plt.title("Particle x-positions over time") +plt.xlabel("time") +plt.ylabel("particle x-position") +h0, = plt.plot(time, particle_one_evolution, label="particle 1") +h1, = plt.plot(time, particle_two_evolution, label="particle 2") +plt.legend(handles=[h0, h1]) plt.show() -- 2.45.1 From dba5664cf3d19e79e8e363de3475bd83e8773f12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 21:51:24 +0200 Subject: [PATCH 7/9] technically everything works, but CoM is not preserved --- c/interaction/interaction.c | 9 +++--- setup.py | 55 +++++++++++++++++++------------------ test.py | 8 ++++-- 3 files changed, 38 insertions(+), 34 deletions(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index 70adbbc..e178b9b 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -115,10 +115,11 @@ interaction_ufunc_float2D // Update the momenta. delta_p_x = delta_x_e * interaction_force_function(r, coefficients); delta_p_y = delty_y_e * interaction_force_function(r, coefficients); - *(float *)(p_x_new + i*p_x_new_steps) -= delta_p_x; - *(float *)(p_y_new + i*p_y_new_steps) -= delta_p_y; - *(float *)(p_x_new + j*p_x_new_steps) += delta_p_x; - *(float *)(p_y_new + j*p_y_new_steps) += delta_p_y; + //printf("%d, %d: %f, %f\n", i, j, delta_p_x, delta_p_y); + *(float *)(p_x_new + i*p_x_new_steps) += delta_p_x; + *(float *)(p_y_new + i*p_y_new_steps) += delta_p_y; + *(float *)(p_x_new + j*p_x_new_steps) -= delta_p_x; + *(float *)(p_y_new + j*p_y_new_steps) -= delta_p_y; } } } diff --git a/setup.py b/setup.py index ec09553..20b4715 100644 --- a/setup.py +++ b/setup.py @@ -1,30 +1,31 @@ -#from distutils.core import setup, Extension -# -# -#interaction = Extension("brown.interaction", -# sources = ["c/interaction/interaction.c"]) -# -#setup(name = "brown", -# version = "0.0.1", -# description = "Me playing around with single-atom classical gases", -# ext_modules = [interaction], -# packages = [ -# "brown" -# ], -# package_dir = {"brown": "py/brown"}, -# #url="https://github.com/daknuett/python3-nf", -# author = "Daniel Knüttel", -# author_email = "daniel.knuettel@daknuett.eu") +from distutils.core import setup, Extension -def configuration(parent_package='', top_path=None): - from numpy.distutils.misc_util import Configuration - config = Configuration('brown' - , parent_package - , top_path) - config.add_extension('interaction', ['c/interaction/interaction.c'], extra_compile_args=["-g"]) - return config +interaction = Extension("brown.interaction", + sources = ["c/interaction/interaction.c"]) -if __name__ == "__main__": - from numpy.distutils.core import setup - setup(configuration=configuration) +setup(name = "brown", + version = "0.0.1", + description = "Me playing around with single-atom classical gases", + ext_modules = [interaction], + packages = [ + "brown" + ], + package_dir = {"brown": "py/brown"}, + #url="https://github.com/daknuett/python3-nf", + author = "Daniel Knüttel", + author_email = "daniel.knuettel@daknuett.eu") + +#def configuration(parent_package='', top_path=None): +# from numpy.distutils.misc_util import Configuration +# +# config = Configuration('brown' +# , parent_package +# , top_path="py/") +# config.add_extension('interaction', ['c/interaction/interaction.c'], extra_compile_args=["-g"]) +# #config.add_subpackage("brown.brown", "py/brown") +# return config +# +#if __name__ == "__main__": +# from numpy.distutils.core import setup +# setup(configuration=configuration) diff --git a/test.py b/test.py index 25b8cd5..beb8dd3 100644 --- a/test.py +++ b/test.py @@ -14,14 +14,16 @@ y_coords = np.array([0, 0], dtype=np.float16) x_momenta = np.array([0, 0], dtype=np.float16) y_momenta = np.array([0, 0], dtype=np.float16) -time = np.arange(0, 5, 1, dtype=np.float16) +time = np.arange(0, 50, 1, dtype=np.float16) time_evolution = deque() +dt = 0.1 + for t in time: x_momenta, y_momenta = interaction2D(x_coords, y_coords, x_momenta, y_momenta) - x_coords += x_momenta - y_coords += y_momenta + x_coords += dt * x_momenta + y_coords += dt * y_momenta time_evolution.append(copy(x_coords)) -- 2.45.1 From c01a3233769b52f16e6f881374b215ca7900de58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 21:51:57 +0200 Subject: [PATCH 8/9] added brown interator --- py/brown/brown.py | 66 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 py/brown/brown.py diff --git a/py/brown/brown.py b/py/brown/brown.py new file mode 100644 index 0000000..320b128 --- /dev/null +++ b/py/brown/brown.py @@ -0,0 +1,66 @@ +from brown.interaction import UFuncWrapper +import numpy as np + +class BrownIterator(object): + def __init__(self + , m + , c + , x + , y + , px + , py + , borders_x + , borders_y + , speed_of_light=1e3 # The value that will replace NaN momenta + , border_dampening=1 + , dt=0.1): + self._max_m = m + self._i = 0 + self.c = c + self.x = x + self.y = y + self.px = px + self.py = py + self.borders_x = borders_x + self.borders_y = borders_y + self.dt = dt + self.speed_of_light = speed_of_light + self.border_dampening = border_dampening + self._interaction = UFuncWrapper(1, c) + def __iter__(self): + self._i = 0 + return self + def __next__(self): + self._i += 1 + if(self._i > self._max_m and self._max_m > 0): + raise StopIteration() + if(self._i == 1): + return self.x, self.y + + self.px, self.py = self._interaction(self.x, self.y, self.px, self.py) + self.px[np.isnan(self.px)] = self.speed_of_light + self.py[np.isnan(self.py)] = self.speed_of_light + self._reflect_at_borders() + + self.x += self.dt * self.px + self.y += self.dt * self.py + + return self.x, self.y + + + def _reflect_at_borders(self): + if(len(self.borders_x) > 0): + self.px[self.x <= self.borders_x[0]] *= -self.border_dampening + self.x[self.x <= self.borders_x[0]] = self.borders_x[0] + if(len(self.borders_x) > 1): + self.px[self.x >= self.borders_x[1]] *= -self.border_dampening + self.x[self.x >= self.borders_x[1]] = self.borders_x[1] + if(len(self.borders_y) > 0): + self.py[self.y <= self.borders_y[0]] *= -self.border_dampening + self.y[self.y <= self.borders_y[0]] = self.borders_y[0] + if(len(self.borders_y) > 1): + self.py[self.y >= self.borders_y[1]] *= -self.border_dampening + self.y[self.y >= self.borders_y[1]] = self.borders_y[1] + + + -- 2.45.1 From f8f0a3818c859e6a7bfc0b53e9bfbcbf2c3e2c72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Fri, 12 Jul 2019 21:52:12 +0200 Subject: [PATCH 9/9] added animation --- particles.py | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 particles.py diff --git a/particles.py b/particles.py new file mode 100644 index 0000000..969bdc3 --- /dev/null +++ b/particles.py @@ -0,0 +1,57 @@ +from brown.interaction import UFuncWrapper +from brown.brown import BrownIterator +import numpy as np +from collections import deque +from copy import copy +import matplotlib.pyplot as plt +import matplotlib.animation as ani + +c = np.array([5, 10, 20, 30, 0, 0, 0, 1, -20, 0, -2, -0.1, 2, 0, 0, 0, 0, 0, 0], dtype=np.float16) +#force_function = UFuncWrapper(0, c) +#interaction2D = UFuncWrapper(1, c) + +borders_x = [-100, 100] +borders_y = [-100, 100] +n_particles = 6 +frames = 1000 + +x_coords = np.random.uniform(borders_x[0] / 2, borders_x[1] / 2, n_particles).astype(np.float16) +y_coords = np.random.uniform(borders_y[0] / 2, borders_y[1] / 2, n_particles).astype(np.float16) + + +x_momenta = np.zeros(n_particles, dtype=np.float16) +y_momenta = np.zeros(n_particles, dtype=np.float16) + + + +fig = plt.figure(figsize=(7, 7)) +ax = fig.add_axes([0, 0, 1, 1], frameon=False) +ax.set_xlim(*borders_x) +ax.set_xticks([]) +ax.set_ylim(*borders_y) +ax.set_yticks([]) + +plot, = ax.plot(x_coords, y_coords, "b.") +center_of_mass, = ax.plot(x_coords.mean(), y_coords.mean(), "r-") +center_of_mass_history_x = deque([x_coords.mean()]) +center_of_mass_history_y = deque([y_coords.mean()]) + +brown = BrownIterator(-1, c + , x_coords, y_coords + , y_momenta, y_momenta + , borders_x, borders_y + , border_dampening=1 + , dt=0.001) + +u = iter(brown) + +def update(i): + data = next(u) + center_of_mass_history_x.append(x_coords.mean()) + center_of_mass_history_y.append(y_coords.mean()) + + plot.set_data(*data) + center_of_mass.set_data(center_of_mass_history_x, center_of_mass_history_y) + +animation = ani.FuncAnimation(fig, update, range(frames), interval=1) +plt.show() -- 2.45.1