c-api vs. ctypes

[brennsuppa]

Hallo zusammen!

Ich wurde berufsbedingt in den "Python"-Topf geworfen und bin mich gerade am einlernen, allerdings gibts jetzt ein paar Probleme.

Ich will aus Performancegründen eine aufwändige Berechnung in C ausführen und möchte nun wissen, welche der obigen Methode ihr mir Vorschlagen würdet...

Mit ctypes bin ich immerhin schon soweit, dass ich auf die mehrdimensionalen Arrays in C zugreifen/herumrechnen kann. Gefällt mir also recht gut, aber bei den flags bin ich etwas kritisch...

Wie ich mittels C-API das Array parsen kann weiß ich noch nicht, da brauche ich etwas Unterstützung, falls dieser Ansatz der saubere ist...

Grüße

[CM]

Hallo und willkommen im Forum,

es wird schlecht möglich sein auf eine Frage zu antworten, wenn man nicht weiß, was die Frage ist und worum es geht (welche Art von Daten / Berechnungen). Magst Du uns mehr verraten?

Gruß,
Christian

PS Neben der C-API und ctypes gibt es noch Cython, SWIG und ggf. auch Python-Libraries für recht effektives Handling von Daten - hängt halt vom Typ der Daten ab ...

[brennsuppa]

Gerne verrate ich mehr

Als Minimalbeispiel will ich eigentlich nur ein ndarray dem C-File übergeben und beispielsweise in einer verschachtelten Schleife jeden Wert um 1 erhöhen...

Im Endeffekt übergebe ich 6 Arrays mit unterschiedlichen Parametern um Temperaturberechnungen zu erstellen

Das Array besteht aus Gleitkommazahlen, in ctypes habe ich float64 dafür verwendet.
Die Doku zur C-API ist halt sehr theoretisch gehalten, was es erschwert "rein zu kommen"

[CM]

Klingt für mich reichlich umständlich:

Code: Alles auswählen

>>> import numpy as np
>>> a = np.zeros((3,3))
>>> a
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
>>> a += 1
>>> a
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])
>>> type(a)
<type 'numpy.ndarray'>

Oder was ist für Dich ein ndarray?

[brennsuppa]

schon klar...
das ganze als minimalbeispiel in C ist nicht ganz so trivial, und da haperts bei mir noch ein wenig bzgl. C-API

Ich habe die Berechnungen in Python schon versucht, dauert ca. 14h - in C dasselbe in etwa 1h... deswegen die umstände

Hab aber gerade ein Beispiel dazu gefunden: http://www.scipy.org/Cookbook/C_Extensions/NumPy_arrays

[CM]

Wie Du meinst. Ich möchte Dich noch auf dies aufmerksam machen -- macht das Leben vielleicht etwas einfacher.

HTH
Christian

[HWK]

Vielleicht wäre das auch etwas für Dich?
MfG
HWK

[HerrHagen]

Ich versuch mal ganz plakativ die Vor- und Nachteile der beiden Varianten darzustellen:

C-API:

Code: Alles auswählen

+ volle Kontrolle über alle Kleinigkeiten der Implementierung
- evtl. Probleme mit neuen Python-Versionen
o gar nicht so kompliziert wie man anfangs denkt

c-types:

Code: Alles auswählen

+ wenig bis keine Probleme bei bei neuer Python Version (abgeschirmt vor Änderungen der C-API, kein Neukompilieren)
+ saubere Trennung von C und Python
- hässlicher Wrapper-Code

cython:

Code: Alles auswählen

+ sehr einfach Geschwindigkeitverbesserung erreichbar, minimaler Aufwand
+ geringster "wrapping"-Aufwand
+ wenig bis keine Probleme bei bei neuer Python Version (abgeschirmt vor Änderungen der C-API, kein Neukompilieren)
- zusätzliche Sprache
- langfristige Existenz des Projekts (???, zumindest in Vgl. zu den beiden anderen Varianten)

f2py

Code: Alles auswählen

+ einfache Handhabung, schnelle Ergebnisse
+ Fortran (ideal für numerische Probleme)
- Fortran (sterbende Sprache)

MFG HerrHagen

[Darii]

Dazu kommt bei ctypes noch, dass die Wahrscheinlichkeit größer ist, dass es mit alternativen Python-Implementierungen läuft(pypy und ironpython haben ctypes). Ich würde im Zweifelsfall eigentlich immer zu ctypes greifen.

[brennsuppa]

Python-Versionen sind kein Problem - bin auf python 2.6 und da bleib ich mit diesem Projekt auch

Ich danke euch allen mal für die neuen Infos!

[brennsuppa]

ich habe jetzt mal diese beiden Fälle auscodiert... vl interessierts wen!
als referenz fürs c_api habe ich http://www.scipy.org/Cookbook/C_Extensions/NumPy_arrays verwendet und als referenz für ctypes http://docs.python.org/extending/extend ... le-example

ctypes-python:

Code: Alles auswählen

__all__ = ['arr']

import numpy as np
import ctypes
import os

_path = os.path.dirname('__file__')
lib = np.ctypeslib.load_library('pyctypesversion', _path)

lib.arr.restype=None
lib.arr.argtypes = [np.ctypeslib.ndpointer(dtype=np.float64,ndim=2,flags='aligned, contiguous, writeable'),
                    np.ctypeslib.ndpointer(dtype=np.float64,ndim=2,flags='aligned'),
                    ctypes.POINTER(np.ctypeslib.c_intp)]

def arr(values):
    values = np.require(values, np.float64, ['ALIGNED'])
    result = np.zeros_like(values)

    lib.arr(result, values, values.ctypes.shape)
    return result

ctypes C-File:

Code: Alles auswählen

void arr(double *res, double *values, int *dim)
{
  int M,N,i,j;
  M = dim[0]; N = dim[1];

  for (i=0; i<N; i++) 
    for (j=0; j<M; j++) 
      res[i*N+j] = values[i*N+j]+3;	   
}

c-api python:

Code: Alles auswählen

import c_api_module
import numpy as n
import sys

def arr(values):
    # .... Check arguments, double NumPy matrices?
    test=n.zeros((2,2)) # create a NumPy matrix as a test object to check values
    typetest= type(test) # get its type

    datatest=test.dtype   # get data type  
    
    if type(values) != typetest:
       raise 'In matsq, matrix argument is not *NumPy* array'
    if values.dtype != datatest:
       raise 'In matsq, matrix argument is not *Float* NumPy matrix'
    # .... Call C extension function
    return c_api_module.arr(values)

c-api C-File:

Code: Alles auswählen

#include "Python.h"
#include "arrayobject.h"
#include "c_api_module.h"
#include <math.h>

/* define methods */
static PyMethodDef c_api_modulemethods[] = {
  {"arr", arr, METH_VARARGS},
  {NULL, NULL}     /* Sentinel - marks the end of this structure */
};

/* initialize test functions */ 
void initc_api_module()  {
  (void) Py_InitModule("c_api_module", c_api_modulemethods);
   import_array();
}


/* array computing
   returns a NEW NumPy array                                   */
static PyObject *arr(PyObject *self, PyObject *args)
{
    PyArrayObject *values, *result;
    double **cin, **cout;
    int i,j,n,m, dims[2];
    
    /* Parse tuples */
    if (!PyArg_ParseTuple(args, "O!", 
        &PyArray_Type, &values))  return NULL;
    if (NULL == values)  return NULL;
    
    /* Check that object input is 'double' type and a matrix */
    if (not_doublematrix(values)) return NULL;
    
    /* Get the dimensions of the input */
    n=dims[0]=values->dimensions[0];
    m=dims[1]=values->dimensions[1];
    
    /* Make a new double matrix of same dims */
    result=(PyArrayObject *) PyArray_SimpleNew(2,dims,NPY_DOUBLE);
        
    /* Change contiguous arrays into C ** arrays (Memory is Allocated!) */
    cin=pymatrix_to_Carrayptrs(values);
    cout=pymatrix_to_Carrayptrs(result);
    
    /* Do the calculation. */
    for ( i=0; i<n; i++)  {
        for ( j=0; j<m; j++)  {
            cout[i][j]= cin[i][j]+3;
    }  }
        
    /* Free memory, close file and return */
    free_Carrayptrs(cin);
    free_Carrayptrs(cout);
    return PyArray_Return(result);
}

double **pymatrix_to_Carrayptrs(PyArrayObject *arrayin) { 
   double **c, *a;
   int i,n,m;

   n=arrayin->dimensions[0];
   m=arrayin->dimensions[1];
   c=ptrvector(n);
   a=(double *) arrayin->data; /* pointer to arrayin data as double */
   for ( i=0; i<n; i++) {
      c[i]=a+i*m; }
   return c;
}
double **ptrvector(long n) { 
   double **v;
   v=(double **)malloc((size_t) (n*sizeof(double)));
   if (!v)   {
      printf("In **ptrvector. Allocation of memory for double array failed.");
      exit(0); }
   return v;
}

void free_Carrayptrs(double **v) { 
   free((char*) v);
}

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

c-api H-File:

Code: Alles auswählen

static PyObject *arr(PyObject *self, PyObject *args);
PyArrayObject *pymatrix(PyObject *objin);
double **pymatrix_to_Carrayptrs(PyArrayObject *arrayin);
double **ptrvector(long n);
void free_Carrayptrs(double **v);
int  not_doublematrix(PyArrayObject *mat);

Ich weiß, das ist jetzt einiges an Code... aber man sieht halt wie unterschiedlich die Ansätze sind, ich werde jetzt die Laufzeit beider Versionen testen und mich danach entscheiden was ich verwenden will![/code]

[Leonidas]

Wo du am vergleichen bist solltest du auch Cython ausprobieren.

[brennsuppa]

ja, will ich noch machen

[HerrHagen]

@brennsuppa: Dein C-Beispiel scheint mir etwas zu kompliziert. Der gängige Weg um auf den Inhalt eines numpy-array zuzugreifen ist PyArray_GETPTR. Das funktioniert dann auch mit nicht-contiguous-arrays (also sowas a[::2, ::3]. Hier mal wie ich es implementiert hätte:

Code: Alles auswählen

#include <Python.h>
#include "numpy/arrayobject.h"

#define GET2_DOUBLE(a, y, x)		(*((double *) PyArray_GETPTR2((a), (y), (x))))

static PyObject *_test_add3(PyObject *self, PyObject *args)
{ 
    PyArrayObject *a;		//in

	if(!PyArg_ParseTuple(args, "O!", &PyArray_Type, &a)) 
		return NULL; 

	int Y = PyArray_DIM(a, 0);
    int X = PyArray_DIM(a, 1);

    PyArrayObject *out;     //out

	npy_intp dims[2] = {Y, X};
	out = (PyArrayObject *) PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);

    int x, y;

	for (y=0; y<Y; y++){
		for (x=0; x<X; x++){
            GET2_DOUBLE(out, y, x) = GET2_DOUBLE(a, y, x) + 3.;
		}
	}

	return PyArray_Return(out);
}


// Methodenliste
static PyMethodDef _testMethods[] = 
    { 
		{"add3", _test_add3, METH_VARARGS, "Addiere 3"}, 
		{NULL, NULL, 0, NULL} 
    };

// Modul-Initialisierung
PyMODINIT_FUNC init_test(void) 
    { 
    Py_InitModule("_test", _testMethods);
	import_array();
    }

MFG HerrHagen

[brennsuppa]

das problem mit den contiguous arrays habe ich bisher "elegant" ignoriert!
Ich danke für dein Beispiel @HerrHagen, werds gleich testen

bzgl. Performance: da liegt ctypes hinten, werds also hart codiert lassen und versuche mittels openMP einen Multicore auszunutzen!

[HerrHagen]

Wäre nett wenn du etwas Code posten würdest wenn du ein Beispiel mit OpenMP zum laufen gekriegt hast. Die Thematik interessiert mich schon seit längeren...

[HerrHagen]

Vieleicht interessant zu diesem Thema:
http://conference.scipy.org/proceedings ... l_text.pdf

[HWK]

Hier einmal eine Lösung mit f2py:
File test.f95

Code: Alles auswählen

subroutine test(arrayin, arrayout, m, n)
  implicit none
  REAL, dimension(m, n), intent(in) :: arrayin
  REAL, dimension(m, n), intent(out) :: arrayout
  integer, intent(in) :: m, n
  arrayout = arrayin + 3
end subroutine test

File testpy.py

Code: Alles auswählen

import numpy
import test

a = numpy.arange(9).reshape(3, 3)
b = test.test(a)
print a
print b

Ausgabe

Code: Alles auswählen

>C:\Programme\Python26\pythonw -u "testpy.py"
[[0 1 2]
 [3 4 5]
 [6 7 8]]
[[  3.   4.   5.]
 [  6.   7.   8.]
 [  9.  10.  11.]]
>Exit code: 0

MfG
HWK

[brennsuppa]

@openMP: ich konnte noch keine Performancesteigerungen feststellen, auch wenn beide Cores jetzt ausgelastet sind...muss den Code wahrsch. noch optimieren!

wenn ich den Code beim compilieren optimiere (-O2) rennt es schneller als auf meinem Dual Core!