Skip to content

Instantly share code, notes, and snippets.

@zsylvester
Created May 19, 2020 22:20
Show Gist options
  • Select an option

  • Save zsylvester/5d6410909fc3733754e5e83d4191f64c to your computer and use it in GitHub Desktop.

Select an option

Save zsylvester/5d6410909fc3733754e5e83d4191f64c to your computer and use it in GitHub Desktop.
attempt to convert dpcore_py.c to Python 3
/*
* dpcore_py.c - calculate dynamic programming inner loop
* Python extension version
* 2014-05-30 Dan Ellis [email protected]
*/
/* see http://wiki.scipy.org/Cookbook/C_Extensions/NumPy_arrays */
#include <Python.h>
#include <numpy/arrayobject.h>
#include <math.h>
#include <numpy/npy_math.h>
static void calc_dpcore(
double *pM, /* input scores */
int rows, /* size of arrays */
int cols,
double pen, /* nondiagonal penalty score */
double *pD, /* output best-cost matrix */
int *pP /* output traceback matrix */
)
{
/* Data is passed in as pointers to contiguous, row-first memory
blocks for each array, that we interpret appropriately here */
int i, j, k, tb;
double d1, d2, d3, v;
double *costs;
int *steps;
int ncosts;
/* setup cost matrix */
int ii;
ncosts = 3;
costs = (double *)malloc(ncosts*sizeof(double));
steps = (int *)malloc(ncosts*2*sizeof(int));
steps[0] = 1; steps[1] = 1; costs[0] = 0.0;
steps[2] = 0; steps[3] = 1; costs[1] = pen;
steps[4] = 1; steps[5] = 0; costs[2] = pen;
/* do dp */
v = pM[0];
tb = 0; /* value to use for 0,0 */
for (j = 0; j < cols; ++j) {
for (i = 0; i < rows; ++i) {
d1 = pM[i*cols + j];
for (k = 0; k < ncosts; ++k) {
if ( i >= steps[2*k] && j >= steps[2*k+1] ) {
d2 = d1 + costs[k] + pD[(i-steps[2*k])*cols + (j-steps[2*k+1])];
if (d2 < v) {
v = d2;
tb = k;
}
}
}
pD[i*cols + j] = v;
pP[i*cols + j] = tb;
v = NPY_INFINITY;
}
}
free((void *)costs);
free((void *)steps);
}
//////////////////////////// python extension wrapper /////////////////////
/* ==== Check that PyArrayObject is a double (Float) type and a matrix ==========
return 1 if an error and raise exception */
int not_doublematrix(PyArrayObject *mat) {
if (mat->descr->type_num != NPY_DOUBLE || mat->nd != 2) {
PyErr_SetString(PyExc_ValueError,
"In not_doublematrix: array must be of type Float and 2 dimensional (n x m).");
return 1; }
return 0;
}
static PyObject *
dpcore_py_dpcore(PyObject *self, PyObject *args)
{
PyArrayObject *Sin, *Dout, *Pout;
double pen;
double *Sptr, *Dptr;
int *Pptr;
int nrows, ncols;
npy_intp dims[2];
/* parse input args */
if (!PyArg_ParseTuple(args, "O!d",
&PyArray_Type, &Sin, &pen))
return NULL;
if (Sin == NULL) return NULL;
/* Check that object input is 'double' type and a matrix
Not needed if python wrapper function checks before call to this routine */
if (not_doublematrix(Sin)) return NULL;
/* Get the dimension of the input */
nrows = Sin->dimensions[0];
ncols = Sin->dimensions[1];
/* Set up output matrices */
dims[0] = nrows;
dims[1] = ncols;
Dout = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_DOUBLE);
Pout = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_INT);
/* Change contiguous arrays into C *arrays */
Sptr = (double *)Sin->data;
Dptr = (double *)Dout->data;
Pptr = (int *)Pout->data;
/* run calculation */
calc_dpcore(Sptr, nrows, ncols, pen, Dptr, Pptr);
/* return the result */
PyObject *tupleresult = PyTuple_New(2);
PyTuple_SetItem(tupleresult, 0, PyArray_Return(Dout));
PyTuple_SetItem(tupleresult, 1, PyArray_Return(Pout));
return tupleresult;
}
/* standard hooks to Python, per http://docs.python.org/2/extending/extending.html */
static PyMethodDef DpcoreMethods[] = {
{"dpcore", dpcore_py_dpcore, METH_VARARGS},
{NULL, NULL} /* Sentinel */
};
/* ==== Initialize the C_test functions ====================== */
// Module name must be _dpcoremodule in compile and linked
static struct PyModuleDef dpcore_definition = {
PyModuleDef_HEAD_INIT,
"_dpcore_py", /* m_name */
"A Python module for dynamic time warping", /* m_doc */
-1, /* m_size */
DpcoreMethods, /* m_methods */
};
PyMODINIT_FUNC PyInit__dpcore_py(void) {
Py_Initialize();
return PyModule_Create(&dpcore_definition);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment