El otro día me planteé la cuestión de cómo crear un módulo para Python que fuera más eficiente que el propio Python, es decir, sin escribirlo en Python. Me puse a mirar la documentación oficial y lo encontré fácilmente, la verdad es que Python está muy bien documentado (cambiaría algunas cosillas, pero más por gustos personales que por que crea que están mal). La principal forma de crear módulos eficientes para Python es hacerlos en C o C++, y ésta es la forma que seguí yo. Antes de continuar, os preguntaréis… ¿Por qué querías hacer un módulo eficiente para Python? Pues… ni idea, jeje, por experimentar, y por ver si con mis conocimientos recién adquiridos (sí, realmente recién adquiridos) de métodos numéricos soy capaz de hacer una buena librería de cálculo científico (en particular para aplicar técnicas de interpolación) para Python. Sólo he visto la PyGMP, pero es posible que sea demasiado general, y no incluye algoritmos de interpolación (que es, en definitiva, lo que yo quiero hacer).
Por ahora sólo he incluido dos algoritmos muy sencillos, el de Neville para evaluar polinomios en un punto (de momento sólo real), y el método de las diferencias divididas de Newton, que nos da las diferencias divididas con las que podemos construir un polinomio interpolatorio. De momento he hecho sólo éstos dos porque he estado haciendo pruebas, y entre ellas han habido pruebas de rendimiento… la verdad es que me quedé pasmado al ver la diferencia entre C y Python (sí, ya sabía que Python era más lento, pero no tanto). El algoritmo de Neville llega a ser un 862% más rápido escrito en C que en Python (tarda sólo un 11.5% de lo que tardaría si estuviera escrito en Python), además, aunque no lo he comprobado, creo que a medida que aumentamos los puntos de soporte, el margen de tiempo aumenta todavía más (supongo que habrá un límite).
Ahora pasaré a explicar brevemente como lo he hecho, aunque para el que quiera profundizar, tendría que leer la documentación de Python directamente, esto debería ser sólo un ejemplo. Lo primero que debemos hacer es crear el código que queremos transformar en un módulo python (ya desde el principio tenemos que empezar utilizando los tipos de la API de Python):
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
| #include <python .h>
/* Neville interpolation on real value */
static PyObject *interpolation_iaval(PyObject *self, PyObject *args)
{
double x;
PyObject *a, *fa;
PyObject *tmp;
int n;
unsigned i, j;
double *da, *dfa;
if (!PyArg_ParseTuple(args, "dOO", &x, &a, &fa))
return PyExc_AttributeError;
if(!(PyList_Check(a) || PyList_Check(fa)))
return PyExc_AttributeError;
n = PyList_Size(a);
if(n!=PyList_Size(fa))
return PyExc_IndexError;
da = (double *)malloc(n*sizeof(double));
dfa = (double *)malloc(n*sizeof(double));
if(da==NULL || dfa==NULL)
{
if(da) free(da);
if(dfa) free(dfa);
return PyErr_NoMemory();
}
for(i=0; i<n ; i++)
{
tmp = PyList_GetItem(a, i);
if(!PyFloat_Check(tmp))
{
free(da);
free(dfa);
return PyExc_AttributeError;
}
da[i] = PyFloat_AsDouble(tmp);
tmp = PyList_GetItem(fa, i);
if(!PyFloat_Check(tmp))
{
free(da);
free(dfa);
return PyExc_AttributeError;
}
dfa[i] = PyFloat_AsDouble(tmp);
}
n--;
for(i=1; i<=n; i++)
for(j=n; j>=i; j--)
dfa[j] = ((x-da[j-i])*dfa[j]-(x-da[j])*dfa[j-1])/(da[j]-da[j-i]);
x = dfa[n];
free(da);
free(dfa);
return Py_BuildValue("d", x);
}
static PyObject *interpolation_divdif(PyObject *self, PyObject *args)
{
PyObject *a, *fa;
PyObject *tmp;
int n;
unsigned i, j;
double *da, *dfa;
if (!PyArg_ParseTuple(args, "OO", &a, &fa))
return PyExc_AttributeError;
if(!(PyList_Check(a) || PyList_Check(fa)))
return PyExc_AttributeError;
n = PyList_Size(a);
if(n!=PyList_Size(fa))
return PyExc_IndexError;
da = (double *)malloc(n*sizeof(double));
dfa = (double *)malloc(n*sizeof(double));
if(da==NULL || dfa==NULL)
{
if(da) free(da);
if(dfa) free(dfa);
return PyErr_NoMemory();
}
for(i=0; i</n><n ; i++)
{
tmp = PyList_GetItem(a, i);
if(!PyFloat_Check(tmp))
{
free(da);
free(dfa);
return PyExc_AttributeError;
}
da[i] = PyFloat_AsDouble(tmp);
tmp = PyList_GetItem(fa, i);
if(!PyFloat_Check(tmp))
{
free(da);
free(dfa);
return PyExc_AttributeError;
}
dfa[i] = PyFloat_AsDouble(tmp);
}
n--;
for(i=1; i<=n; i++)
for(j=n; j>=i; j--)
dfa[j] = (dfa[j]-dfa[j-1])/(da[j]-da[j-i]);
n++;
tmp = PyList_New((Py_ssize_t)(n));
for(i=0; i</n><n ; i++)
PyList_SetItem(tmp, (Py_ssize_t)i, PyFloat_FromDouble(dfa[i]));
return tmp;
}
static PyMethodDef interpolationmethods[] = {
{"iaval", interpolation_iaval, METH_VARARGS,
"float iaval(x,[x1,x2,...,xn], [f(x1),f(x2),..,f(xn)])\n"
"Neville's method for Lagrange's interpolation\n"
"This function returns the value of the interpolation polynomial on the given point (x)." },
{"divdif", interpolation_divdif, METH_VARARGS,
"divdif([x1,x2,...,xn], [f(x1),f(x2),...,f(xn)])\n"
"Newton's Divided differences method for Lagrange's Interpolation\n"
"This function returns the divided differences of the Newton's divided differences interpolation polynomial."},
{NULL, NULL} /* sentinel */
};
DL_EXPORT(void) initinterpolation()
{
Py_InitModule("interpolation", interpolationmethods);
} |
Una vez escrito el código tendría que explicarlo… pero en serio, da pereza explicarlo, :p, mirad la referencia de Python online http://docs.python.org/ext/ext.html (Extendiendo Python, en inglés), http://docs.python.org/api/api.html (Python/C API Reference Manual). Después de entender lo que hemos escrito, tenemos que hacer que sea accesible desde el intérprete Python. Desde al versión 2.3 tenemos un método muy sencillo. Este se basa en crear una especie de “instalador” en un pequeño archivo Python que cuando es interpretado construye e instala el módulo para que pueda ser utilizado sin ningún problema.
1
2
3
4
5
6
7
8
9
| from distutils.core import setup, Extension
module1 = Extension('interpolation',
sources = ['interpolationmodule.c'])
setup (name = 'interpolation',
version = '0.1',
description = 'Interpolation Package',
ext_modules = [module1]) |
Por cierto, el fichero C lo llamé interpolationmodule.c, y a éste fichero Python lo llamaremos setup.py. Para construir el módulo escribimos python setup.py build . Para instalarlo debemos entrar en modo administrador y escribir lo siguiente: python setup.py install . Y con esto finaliza la instalación del módulo. Si queréis probarlo sólo hace falta escribir import interpolation en el intérprete, las funciones disponibles son iaval (método de Neville) y divdif (método de las diferencias divididas). Para los que queráis probar las diferecias de rendimiento, podéis importar el módulo timeit, que se utiliza como sigue.
Asignamos a una variable un objeto Timer:
t = timeit.Timer("interpolation.iaval(3.0, [2.0, 4.0], [5.0, 6.0])”, “import interpolation”)
Luego ejecutamos la función timeit() del objeto que hemos creado. Ésta función lo que hará será llamar a la función que le hemos indicado un millón de veces dentro del entorno que le indicamos en el segundo parámetro. Una vez acabe su trabajo, nos indicará el tiempo transcurrido. Si creamos una función en Python que haga el mismo trabajo y le damos los mismos datos, veremos con timeit la gran diferencia de rendimiento que hay entre ambas.
Bueno, acabo este artículo nada reflexivo, aunque espero, entretenido para frikis como yo.
.
.
