brown/c/interaction/interaction.c

#include "interaction.h"
#include "math.h"


#define raise2(x) (x)*(x)

static float 
interaction_force_function
	( float r
	, float * coefficients)
{
	float result = 0;
	int i;
	for(i = 0; i < 7; i++)
	{
		result += coefficients[i] * powf(r, i);
	}
	result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9]));
	result += coefficients[10] * expf(coefficients[11] * (r - coefficients[12]));
	result += coefficients[13] * expf(coefficients[14] * raise2(r - coefficients[15]));
	result += coefficients[16] * expf(coefficients[17] * raise2(r - coefficients[18]));
	return result;
}

static void
interaction_ufunc_force
	( char ** args
	, npy_intp * dimensions
	, npy_intp * steps
	, void * data)
{
	char * in = args[0]
		, * out = args[1]
		, * raw_coefficients = args[2];
	npy_intp n = dimensions[0];
	npy_intp in_step = steps[0]
		, out_step = steps[1]
		, raw_coefficients_steps = steps[2];

	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;
	}

	for(i = 0; i < n; i++)
	{
		*(float *)out = interaction_force_function(*(float *)in, coefficients);
		out += out_step;
		in += in_step;
	}
}
	

static void
interaction_ufunc_float2D
	( char ** args
	, npy_intp * dimensions
	, npy_intp * steps
	, void * data)
{
	npy_intp i;
	npy_intp j;
	npy_intp n = dimensions[0];

	char * x_old = args[0]
		, * 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];

	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];

	// 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;
	}

	// Compute the new momenta: 

	// Stuff we will need:
	float r;         // Distance between interacting particles.
	float delta_x_e; // x-component of unit direction vector.
	float delty_y_e; // y-component of unit direction vector.
	float delta_p_x;
	float delta_p_y;

	// The current x_old[i], y_old[i] coordinates.
	float this_x_i;
	float this_y_i;
	// The current x_old[j], y_old[j] coordinates.
	float this_x_j;
	float this_y_j;


	for(i = 0; i < n; i++)
	{
		this_x_i = *(float *)(x_old + i*x_old_steps);
		this_y_i = *(float *)(y_old + i*y_old_steps);
		
		// copy current momenta
		*(float *)(p_x_new + i*p_x_new_steps) = *(float *)(p_x_old + i*p_x_old_steps);
		*(float *)(p_y_new + i*p_y_new_steps) = *(float *)(p_y_old + i*p_y_old_steps);

		// compute and add the momentum offset
		for(j = 0; j < i; j++)
		{
			this_x_j = *(float *)(x_old + j*x_old_steps);
			this_y_j = *(float *)(y_old + j*y_old_steps);

			// Compute distance and direction between particles i,j.
			r = sqrtf(raise2(this_x_i - this_x_j) + raise2(this_y_i - this_y_j));
			delta_x_e = (this_x_i - this_x_j) / r;
			delty_y_e = (this_y_i - this_y_j) / r;

			// 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;
		}
	}
}

static PyMethodDef InteractionMethods[] = {
	{NULL, NULL, 0, NULL}
};

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
	, "brown.interaction"
	, NULL
	, -1
	, InteractionMethods
	, NULL
	, NULL
	, NULL
	, NULL
};

PyMODINIT_FUNC 
PyInit_interaction(void)
{
	PyObject * module
		, * ufunc_interaction
		, * ufunc_force
		, * dct;

	module = PyModule_Create(&moduledef);
	if(!module)
	{
		return NULL;
	}
	import_array();
	import_umath();

	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);
	PyDict_SetItemString(dct, "force_function", ufunc_force);
	Py_DECREF(ufunc_interaction);
	Py_DECREF(ufunc_force);

	return module;
}