\n", "/* Header to test of C modules for arrays for Python: C_test.c */\n", "\n", "/* ==== Prototypes =================================== */\n", "\n", "// .... Python callable Vector functions ..................\n", "static PyObject *vecfcn1(PyObject *self, PyObject *args);\n", "static PyObject *vecsq(PyObject *self, PyObject *args);\n", "\n", "/* .... C vector utility functions ..................*/\n", "PyArrayObject *pyvector(PyObject *objin);\n", "double *pyvector_to_Carrayptrs(PyArrayObject *arrayin);\n", "int not_doublevector(PyArrayObject *vec);\n", "\n", "\n", "/* .... Python callable Matrix functions ..................*/\n", "static PyObject *rowx2(PyObject *self, PyObject *args);\n", "static PyObject *rowx2_v2(PyObject *self, PyObject *args);\n", "static PyObject *matsq(PyObject *self, PyObject *args);\n", "static PyObject *contigmat(PyObject *self, PyObject *args);\n", "\n", "/* .... C matrix utility functions ..................*/\n", "PyArrayObject *pymatrix(PyObject *objin);\n", "double **pymatrix_to_Carrayptrs(PyArrayObject *arrayin);\n", "double **ptrvector(long n);\n", "void free_Carrayptrs(double **v);\n", "int not_doublematrix(PyArrayObject *mat);\n", "\n", "/* .... Python callable integer 2D array functions ..................*/\n", "static PyObject *intfcn1(PyObject *self, PyObject *args);\n", "\n", "/* .... C 2D int array utility functions ..................*/\n", "PyArrayObject *pyint2Darray(PyObject *objin);\n", "int **pyint2Darray_to_Carrayptrs(PyArrayObject *arrayin);\n", "int **ptrintvector(long n);\n", "void free_Cint2Darrayptrs(int **v);\n", "int not_int2Darray(PyArrayObject *mat);\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Source file, C_arraytest.c:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "/* A file to test imorting C modules for handling arrays to Python */\n", "\n", "#include \"Python.h\"\n", "#include \"arrayobject.h\"\n", "#include \"C_arraytest.h\"\n", "#include" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Now, lets look at the source code in smaller chunks.\n", "\n", "#### Headers\n", "\n", "You must include the following headers with Python.h **always** the first header included.\n", "\n", " #include \"Python.h\"\n", " #include \"arrayobject.h\"\n", "\n", "I also include the header C_arraytest.h which contains the prototype of the matsq function:\n", "\n", " static PyObject *matsq(PyObject *self, PyObject *args);\n", "\n", "The static keyword in front of a function declaration makes this function private to your extension module. The linker just won't see it. This way you can use the same intuitional function names(i.e. sum, check, trace) for all extension modules without having name clashes between them at link time. The type of the function is `PyObject *` because it will always be returning to a Python calling function so you can (must, actually) return a Python object. The arguments are always the same,\n", "\n", " PyObject *self and PyObject *args\n", "\n", "The first one self is never used, but necessary because of how Python passes arguments. The second args is a pointer to a Python tuple that contains all of the arguments (B,i,x) of the function.\n", "\n", "#### Method definitions\n", "\n", "This sets up a table of function names that will be the interface from your Python code to your C extension. The name of the C extension module will be `_C_arraytest` (note the leading underscore). It is important to get the name right each time it is used because there are strict requirements on using the module name in the code. The name appears first in the method definitions table as the first part of the table name:\n", "\n", " static PyMethodDef _C_arraytestMethods[] = {\n", " ...,\n", " {\"matsq\", matsq, METH_VARARGS},\n", " ...,\n", " {NULL, NULL} /* Sentinel - marks the end of this structure */\n", " };\n", "\n", "where I used ellipses (...) to ignore other code not relevant to this function. The `METH_VARARGS` parameter tells the compiler that you will pass the arguments the usual way without keywords as in the example `A=matsq(B,i,x)` above. There are ways to use Python keywords, but I have not tried them out. The table should always end with {NULL, NULL} which is just a \"marker\" to note the end of the table.\n", "\n", "#### Initializations\n", "\n", "These functions tell the Python interpreter what to call when the module is loaded. Note the name of the module (`_C_arraytest`) must come directly after the init in the name of the initialization structure.\n", "\n", " void init_C_arraytest() {\n", " (void) Py_InitModule(\"_C_arraytest\", _C_arraytestMethods);\n", " import_array(); // Must be present for NumPy. Called first after above line.\n", " }\n", "\n", "The order is important and you must call these two initialization functions first.\n", "\n", "#### The matsqfunction code\n", "\n", "Now here is the actual function that you will call from Python code. I will split it up and explain each section.\n", "\n", "The function name and type:\n", "\n", " static PyObject *matsq(PyObject *self, PyObject *args)\n", "\n", "You can see they match the prototype in C_arraytest.h.\n", "\n", "The local variables:\n", "\n", " PyArrayObject *matin, *matout;\n", " double **cin, **cout, dfactor;\n", " int i,j,n,m, dims[2], ifactor;\n", "\n", "The `PyArrayObjects` are structures defined in the NumPy header file and they will be assigned pointers to the actual input and output NumPy arrays (A and B). The C arrays `cin` and `cout` are Cpointers that will point (eventually) to the actual data in the NumPy arrays and allow you to manipulate it. The variable `dfactor` will be the Python float y, `ifactor` will be the Python int i, the variables i,j,n, and m will be loop variables (i and j) and matrix dimensions (n= number of rows, m= number of columns) in A and B. The array dims will be used to access n and m from the NumPy array. All this happens below. First we have to extract the input variables (A, i, y) from the args tuple. This is done by the call,\n", "\n", " /* Parse tuples separately since args will differ between C fcns */\n", " if (!PyArg_ParseTuple(args, \"O!id\",\n", " &PyArray_Type, &matin, &ifactor, &dfactor)) return NULL;\n", "\n", "The `PyArg_ParseTuple` function takes the args tuple and using the format string that appears next (\"O!id\" ) it assigns each member of the tuple to a C variable. Note you must pass all C variables by reference. This is true even if the C variable is a pointer to a string (see code in vecfcn1 routine). The format string tells the parsing function what type of variable to use. The common variables for Python all have letter names (e.g. s for string, i for integer, d for (double - the Python float)). You can find a list of these and many more in Guido's tutorial (http://docs.python.org/ext/ext.html). For data types that are not in standard Python like the NumPy arrays you use the O! notation which tells the parser to look for a type structure (in this case a NumPy structure `PyArray_Type`) to help it convert the tuple member that will be assigned to the local variable ( matin ) pointing to the NumPy array structure. Note these are also passed by reference. The order must be maintained and match the calling interface of the Python function you want. The format string defines the interface and if you do not call the function from Python so the number of arguments match the number in the format string, you will get an error. This is good since it will point to where the problem is.\n", "\n", "If this doesn't work we return NULL which will cause a Python exception.\n", "\n", " if (NULL == matin) return NULL;\n", "\n", "Next we have a check that the input matrix really is a matrix of NumPy type double. This test is also done in the Python wrapper for this C extension. It is better to do it there, but I include the test here to show you that you can do testing in the C extension and you can \"reach into\" the NumPy structure to pick out it's parameters. The utility function `not_doublematrix` is explained later.\n", "\n", " /* Check that object input is 'double' type and a matrix\n", " Not needed if python wrapper function checks before call to this routine */\n", " if (not_doublematrix(matin)) return NULL;\n", "\n", "Here's an example of reaching into the NumPy structure to get the dimensions of the matrix matin and assign them to local variables as mentioned above.\n", "\n", " /* Get the dimensions of the input */\n", " n=dims[0]=matin->dimensions[0];\n", " m=dims[1]=matin->dimensions[1];\n", "\n", "Now we use these matrix parameters to generate a new NumPy matrix matout (our output) right here in our C extension. PyArray_FromDims(2,dims,NPY_DOUBLE) is a utility function provided by NumPy (not me) and its arguments tell NumPy the rank of the NumPy object (2), the size of each dimension (dims), and the data type (NPY_DOUBLE). Other examples of creating different NumPy arrays are in the other C extensions.\n", "\n", " /* Make a new double matrix of same dims */\n", " matout=(PyArrayObject *) PyArray_FromDims(2,dims,NPY_DOUBLE);\n", "\n", "To actually do our calculations we need C structures to handle our data so we generate two C 2-dimensional arrays (cin and cout) which will point to the data in matin and matout, respectively. Note, here memory is allocated since we need to create an array of pointers to C doubles so we can address cin and cout like usual C matrices with two indices. This memory must be released at the end of this C extension. Memory allocation like this is not always necessary. See the routines for NumPy vector manipulation and treating NumPy matrices like contiguous arrays (as they are in NumPy) in the C extension (the routine contigmat).\n", "\n", " /* Change contiguous arrays into C ** arrays (Memory is Allocated!) */\n", " cin=pymatrix_to_Carrayptrs(matin);\n", " cout=pymatrix_to_Carrayptrs(matout);\n", "\n", "Finally, we get to the point where we can manipulate the matrices and do our calculations. Here is the part where the original equation operations $B_{ij}= i y (A_{ij})^2$ are carried out. Note, we are directly manipulating the data in the original NumPy arrays A and B passed to this extension. So anything you do here to the components of cin or cout will be done to the original matrices and will appear there when you return to the Python code.\n", "\n", " /* Do the calculation. */\n", " for ( i=0; i\n", "\n", "/* #### Globals #################################### */\n", "\n", "/* ==== Set up the methods table ====================== */\n", "static PyMethodDef _C_arraytestMethods[] = {\n", "\t{\"vecfcn1\", vecfcn1, METH_VARARGS},\n", "\t{\"vecsq\", vecsq, METH_VARARGS},\n", "\t{\"rowx2\", rowx2, METH_VARARGS},\n", "\t{\"rowx2_v2\", rowx2_v2, METH_VARARGS},\n", "\t{\"matsq\", matsq, METH_VARARGS},\n", "\t{\"contigmat\", contigmat, METH_VARARGS},\n", "\t{\"intfcn1\", intfcn1, METH_VARARGS},\n", "\t{NULL, NULL} /* Sentinel - marks the end of this structure */\n", "};\n", "\n", "/* ==== Initialize the C_test functions ====================== */\n", "// Module name must be _C_arraytest in compile and linked \n", "void init_C_arraytest() {\n", "\t(void) Py_InitModule(\"_C_arraytest\", _C_arraytestMethods);\n", "\timport_array(); // Must be present for NumPy. Called first after above line.\n", "}\n", "\n", "/* #### Vector Extensions ############################## */\n", "\n", "/* ==== vector function - manipulate vector in place ======================\n", " Multiply the input by 2 x dfac and put in output\n", " Interface: vecfcn1(vec1, vec2, str1, d1)\n", " vec1, vec2 are NumPy vectors, \n", " str1 is Python string, d1 is Python float (double)\n", " Returns integer 1 if successful */\n", "static PyObject *vecfcn1(PyObject *self, PyObject *args)\n", "{\n", "\tPyArrayObject *vecin, *vecout; // The python objects to be extracted from the args\n", "\tdouble *cin, *cout; // The C vectors to be created to point to the \n", "\t // python vectors, cin and cout point to the row\n", "\t // of vecin and vecout, respectively\n", "\tint i,j,n;\n", "\tconst char *str;\n", "\tdouble dfac;\n", "\t\n", "\t/* Parse tuples separately since args will differ between C fcns */\n", "\tif (!PyArg_ParseTuple(args, \"O!O!sd\", &PyArray_Type, &vecin,\n", "\t\t&PyArray_Type, &vecout, &str, &dfac)) return NULL;\n", "\tif (NULL == vecin) return NULL;\n", "\tif (NULL == vecout) return NULL;\n", "\t\n", "\t// Print out input string\n", "\tprintf(\"Input string: %s\\n\", str);\n", "\t\n", "\t/* Check that objects are 'double' type and vectors\n", "\t Not needed if python wrapper function checks before call to this routine */\n", "\tif (not_doublevector(vecin)) return NULL;\n", "\tif (not_doublevector(vecout)) return NULL;\n", "\t\n", "\t/* Change contiguous arrays into C * arrays */\n", "\tcin=pyvector_to_Carrayptrs(vecin);\n", "\tcout=pyvector_to_Carrayptrs(vecout);\n", "\t\n", "\t/* Get vector dimension. */\n", "\tn=vecin->dimensions[0];\n", "\t\n", "\t/* Operate on the vectors */\n", "\tfor ( i=0; i dimensions[0];\n", "\t\n", "\t/* Make a new double vector of same dimension */\n", "\tvecout=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);\n", "\t\t\n", "\t/* Change contiguous arrays into C *arrays */\n", "\tcin=pyvector_to_Carrayptrs(vecin);\n", "\tcout=pyvector_to_Carrayptrs(vecout);\n", "\t\n", "\t/* Do the calculation. */\n", "\tfor ( i=0; i dimensions[0];\n", "\treturn (double *) arrayin->data; /* pointer to arrayin data as double */\n", "}\n", "/* ==== Check that PyArrayObject is a double (Float) type and a vector ==============\n", " return 1 if an error and raise exception */ \n", "int not_doublevector(PyArrayObject *vec) {\n", "\tif (vec->descr->type_num != NPY_DOUBLE || vec->nd != 1) {\n", "\t\tPyErr_SetString(PyExc_ValueError,\n", "\t\t\t\"In not_doublevector: array must be of type Float and 1 dimensional (n).\");\n", "\t\treturn 1; }\n", "\treturn 0;\n", "}\n", "\n", "/* #### Matrix Extensions ############################## */\n", "\n", "/* ==== Row x 2 function - manipulate matrix in place ======================\n", " Multiply the 2nd row of the input by 2 and put in output\n", " interface: rowx2(mat1, mat2)\n", " mat1 and mat2 are NumPy matrices\n", " Returns integer 1 if successful */\n", "static PyObject *rowx2(PyObject *self, PyObject *args)\n", "{\n", "\tPyArrayObject *matin, *matout; // The python objects to be extracted from the args\n", "\tdouble **cin, **cout; // The C matrices to be created to point to the \n", "\t // python matrices, cin and cout point to the rows\n", "\t // of matin and matout, respectively\n", "\tint i,j,n,m;\n", "\t\n", "\t/* Parse tuples separately since args will differ between C fcns */\n", "\tif (!PyArg_ParseTuple(args, \"O!O!\", &PyArray_Type, &matin,\n", "\t\t&PyArray_Type, &matout)) return NULL;\n", "\tif (NULL == matin) return NULL;\n", "\tif (NULL == matout) return NULL;\n", "\t\n", "\t/* Check that objects are 'double' type and matrices\n", "\t Not needed if python wrapper function checks before call to this routine */\n", "\tif (not_doublematrix(matin)) return NULL;\n", "\tif (not_doublematrix(matout)) return NULL;\n", "\t\t\n", "\t/* Change contiguous arrays into C ** arrays (Memory is Allocated!) */\n", "\tcin=pymatrix_to_Carrayptrs(matin);\n", "\tcout=pymatrix_to_Carrayptrs(matout);\n", "\t\n", "\t/* Get matrix dimensions. */\n", "\tn=matin->dimensions[0];\n", "\tm=matin->dimensions[1];\n", "\t\n", "\t/* Operate on the matrices */\n", "\tfor ( i=0; i dimensions[0];\n", "\tm=matin->dimensions[1];\n", "\t\n", "\t/* Operate on the matrices */\n", "\tfor ( i=0; i dimensions[0];\n", "\tm=dims[1]=matin->dimensions[1];\n", "\t\n", "\t/* Make a new double matrix of same dims */\n", "\tmatout=(PyArrayObject *) PyArray_FromDims(2,dims,NPY_DOUBLE);\n", "\t\t\n", "\t/* Change contiguous arrays into C ** arrays (Memory is Allocated!) */\n", "\tcin=pymatrix_to_Carrayptrs(matin);\n", "\tcout=pymatrix_to_Carrayptrs(matout);\n", "\t\n", "\t/* Do the calculation. */\n", "\tfor ( i=0; i dimensions[0];\n", "\tm=dims[1]=matin->dimensions[1];\n", "\tncomps=n*m;\n", "\t\n", "\t/* Make a new double matrix of same dims */\n", "\tmatout=(PyArrayObject *) PyArray_FromDims(2,dims,NPY_DOUBLE);\n", "\t\t\n", "\t/* Change contiguous arrays into C * arrays pointers to PyArrayObject data */\n", "\tcin=pyvector_to_Carrayptrs(matin);\n", "\tcout=pyvector_to_Carrayptrs(matout);\n", "\t\n", "\t/* Do the calculation. */\n", "\tprintf(\"In contigmat, cout (as contiguous memory) =\\n\");\n", "\tfor ( i=0; i dimensions[0];\n", "\tm=arrayin->dimensions[1];\n", "\tc=ptrvector(n);\n", "\ta=(double *) arrayin->data; /* pointer to arrayin data as double */\n", "\tfor ( i=0; i descr->type_num != NPY_DOUBLE || mat->nd != 2) {\n", "\t\tPyErr_SetString(PyExc_ValueError,\n", "\t\t\t\"In not_doublematrix: array must be of type Float and 2 dimensional (n x m).\");\n", "\t\treturn 1; }\n", "\treturn 0;\n", "}\n", "\n", "/* #### Integer 2D Array Extensions ############################## */\n", "\n", "/* ==== Integer function - manipulate integer 2D array in place ======================\n", " Replace >=0 integer with 1 and < 0 integer with 0 and put in output\n", " interface: intfcn1(int1, afloat)\n", " int1 is a NumPy integer 2D array, afloat is a Python float\n", " Returns integer 1 if successful */\n", "static PyObject *intfcn1(PyObject *self, PyObject *args)\n", "{\n", "\tPyArrayObject *intin, *intout; // The python objects to be extracted from the args\n", "\tint **cin, **cout; // The C integer 2D arrays to be created to point to the \n", "\t // python integer 2D arrays, cin and cout point to the rows\n", "\t // of intin and intout, respectively\n", "\tint i,j,n,m, dims[2];\n", "\tdouble afloat;\n", "\t\n", "\t/* Parse tuples separately since args will differ between C fcns */\n", "\tif (!PyArg_ParseTuple(args, \"O!d\", \n", "\t\t&PyArray_Type, &intin, &afloat)) return NULL;\n", "\tif (NULL == intin) return NULL;\n", "\t\n", "\tprintf(\"In intfcn1, the input Python float = %e, a C double\\n\",afloat);\n", "\t\n", "\t/* Check that object input is int type and a 2D array\n", "\t Not needed if python wrapper function checks before call to this routine */\n", "\tif (not_int2Darray(intin)) return NULL;\n", "\t\n", "\t/* Get the dimensions of the input */\n", "\tn=dims[0]=intin->dimensions[0];\n", "\tm=dims[1]=intin->dimensions[1];\n", "\t\n", "\t/* Make a new int array of same dims */\n", "\tintout=(PyArrayObject *) PyArray_FromDims(2,dims,NPY_LONG);\n", "\t\t\n", "\t/* Change contiguous arrays into C ** arrays (Memory is Allocated!) */\n", "\tcin=pyint2Darray_to_Carrayptrs(intin);\n", "\tcout=pyint2Darray_to_Carrayptrs(intout);\n", "\t\n", "\t/* Do the calculation. */\n", "\tfor ( i=0; i = 0) {\n", "\t\t\t\tcout[i][j]= 1; }\n", "\t\t\telse {\n", "\t\t\t\tcout[i][j]= 0; }\n", "\t} }\n", "\t\n", "\tprintf(\"In intfcn1, the output array is,\\n\\n\");\n", "\n", "\tfor ( i=0; i dimensions[0];\n", "\tm=arrayin->dimensions[1];\n", "\tc=ptrintvector(n);\n", "\ta=(int *) arrayin->data; /* pointer to arrayin data as int */\n", "\tfor ( i=0; i descr->type_num != NPY_LONG || mat->nd != 2) {\n", "\t\tPyErr_SetString(PyExc_ValueError,\n", "\t\t\t\"In not_int2Darray: array must be of type int and 2 dimensional (n x m).\");\n", "\t\treturn 1; }\n", "\treturn 0;\n", "}\n", "\n", "// EOF\n", "