Може бути багато причин, щоб використовувати код на С в інтерпретаторі Python, наприклад:
- наявність бібліотек, написаних лише на С;
- робота з hardware, де існують виклики з нативними для С структурами даних;
- пришвидшення роботи коду.
Зупинимось на останньому пункті та розглянемо бібліотеку ctypes і python/c API, які дають можливість запускати код на С в інтерпретаторі Python. Порівняння здійснювалось на найпростішому алгоритмі пошуку власних векторів матриці — степеневому методі (power iteration).
Розділи статті:
- Power iteration
- Реалізація на чистому Python
- Реалізація на чистому С
- ctypes
- python/c api
- benchmarks
- Висновки
Power iteration
Power iteration є надзвичайно легким у реалізації алгоритмом пошуку власних векторів. Зазвичай він потрібен для роботи з розрідженими матрицями. Хоча алгоритм і простий, але використовується, наприклад, пошуковою системою Google (для розрахунку рангу сторінок) та в рекомендаційній системі Twitter. Наша стаття оглядова, тому цей метод добре нам підходить: через легкість реалізації та неоптимальність роботи (наявність великої кількості циклів, що може досить сповільнити роботу скриптів Python).
Детальне математичне викладення методу та доведення його збіжності є на вікіпедії. У статті наведемо лише кроки його імплементації. У спрощеному варіанті вони будуть такими:
- Робимо припущення про значення власного вектора.
- Знаходимо добуток матриці та пробного власного вектора — на виході отримуємо новий вектор.
- Виконується нормалізація отриманого вектора та перевірка на збіжність (різниця між значеннями, отриманими на поточній та попередній ітераціях).
- Якщо на попередньому кроці збіжність досягнена, отриманий вектор — власний вектор матриці.
Реалізація на Python
Синтаксис Python робить перше знайомство з реалізацією алгоритму надзвичайно приємним. Зауважимо, що більшість конструкцій імплементовані в найочевидніший спосіб і без додаткових бібліотек — для легшого сприйняття та чистоти порівняння бенчмарків наприкінці статті.
pure_python.py
def dot(matrix, vector):
""" Calculates the dot product of matrix and vector. Result is new vector. """
product = []
for row in matrix:
row_product = 0
for a_element, b_element in zip(row, vector):
row_product += a_element * b_element
product.append(row_product)
return product
def is_equal(a, b, precision=0.000001):
""" Checks equality between two float numbers with some precision"""
for a_element, b_element in zip(a, b):
if abs(a_element - b_element) > precision:
return False
return True
eigenvector = INITIAL_GUESS # Make initial guess for eigenvector
while not is_equal(eigenvector, prev_value):
prev_value = eigenvector
eigenvector = dot(TEST_MATRIX, eigenvector)
norm = max(eigenvector, key=abs)
eigenvector = [element / norm for element in eigenvector]
Реалізація на С
Далі представлена імплементація алгоритму на С, яка не особливо відрізняється від попереднього варіанту. Єдине, тут потрібен додатковий метод copy_array для копіювання масиву та специфіки роботи з пам'яттю.
pure_c.c
float *dot(float matrix[N][N], float* vector) {
float *result_vector = (float*) malloc(N * sizeof(float));
float row_column_product;
float element_value;
for (int row_index=0; row_index < N; row_index++){
row_column_product = 0;
for (int col_index=0; col_index < N; col_index++) {
element_value = *(vector + col_index);
row_column_product += matrix[row_index][col_index] * element_value;
}
result_vector[row_index] = row_column_product;
}
return result_vector;
}
float *copy_array(float *array) {
float *result_vector = (float*) malloc(N * sizeof(float));
for (int i=0; i < N; i++) {
result_vector[i] = *(array+i);
}
return result_vector;
}
float *normalize(float* vector) {
float *result_vector = (float*) malloc(N * sizeof(float));
float max_value = 0;
for (int i=0; i < N; i++) {
if ( fabs(*(vector + i)) > fabs(max_value) ) {
max_value = *(vector + i);
}
}
for (int i=0; i < N; i++) {
result_vector[i] = *(vector + i) / max_value;
}
return result_vector;
}
bool is_equal(float* a, float* b) {
for (int i=0; i < N; i++) {
if ( fabs(*(a+i) - *(b+i)) > PRECISION ) {
return false;
}
}
return true;
}
int main(int argc, char** argv) {
double* prev_value;
double* vector = INITIAL_GUESS;
do {
prev_value = copy_array(vector);
vector = dot(TEST_MATRIX, vector);
vector = normalize(vector);
} while ( !is_equal(vector, prev_value) );
return 0;
}
Описані функції використовуватимуться в подальших реалізаціях через їхні виклики за допомогою ctypes та python/c API.
ctypes
ctypes — бібліотека для Python, що забезпечує сумісні з С типи даних. Також з нею можна робити виклики функцій динамічних бібліотек. Ось цікавий туторіал про використання цієї бібліотеки на простих прикладах і з варіантами ООП.
Основна складність роботи з бібліотекою полягає у наведенні складних типів даних — через неочевидний для Python синтаксис, наприклад (ctypes.c_float * len(list_to_pointer))(*list_to_pointer). Для цих операцій було введено дві додаткові функції: cast — для списків, та case2d — для тестової матриці.
Щоб приєднати динамічну бібліотеку, код на С спочатку компілюється у певний спосіб (варіант для Linux з використанням gcc):
gcc -shared -Wl,-install_name,testlib.so -o testlib.so -fPIC testlib.c
А потім потрібно ініціалізувати в Python-скрипті pure_c = ctypes.CDLL('ctype_c.so'). Далі можуть бути вказані типи вхідних та вихідних параметрів, що частково спрощує життя надалі.
ctypes_python.py
def case2d(matrix):
""" Cast python represented matrix to ctypes """
return (ctypes.c_float * N * N)(*(tuple(matrix_el) for matrix_el in matrix))
def cast(list_to_pointer: list):
""" Cast python represented list/vector to ctypes """
return (ctypes.c_float * len(list_to_pointer))(*list_to_pointer)
pure_c = ctypes.CDLL('ctype_c.so')
# defines return type of C functions, so ctypes can convert their result in builtin python datatype.
pure_c.dot.restype = ctypes.POINTER(ctypes.c_float * N)
pure_c.normalize.restype = ctypes.POINTER(ctypes.c_float * N)
pure_c.is_equal.restype = ctypes.POINTER(ctypes.c_bool)
eigenvector = INITIAL_GUESS
while not pure_c.is_equal(cast(eigenvector), cast(prev_value)):
prev_value = eigenvector
eigenvector = pure_c.dot(case2d(TEST_MATRIX), cast(eigenvector)).contents
eigenvector = pure_c.normalize(cast(eigenvector)).contents
На цьому кроці процес інтеграції видається досить простим, окрім моментів конвертації складних структур Python в C. Але складність останніх не варто недооцінювати. Для ґрунтовніших задач це стане досить нетривіальною проблемою, особливо якщо врахувати читабельність функції case2d для інших розробників. А от зручність роботи в тому, що всі кроки з інтеграції мов відбуваються в Python-скрипті, тож не потрібно зайвий раз відкривати код, написаний на С.
python/c api
С extension є передбаченою (на рівні інтерпретатора) можливістю розширювати Python-код модулями на С. Це робиться приведенням структур даних останнього до базових типів самого інтерпретатора Python. Альтернативний туторіал про використання C extention можна глянути тут.
Щоб так реалізувати power iteration, перш за все необхідно імпортувати С-модуль (що міститься в самому інтерпретаторі python), модуль для виконання математичних операцій та інші допоміжні бібліотеки:
power_iteration.c
#include <Python.h>
#include <math.h>
#include <stdbool.h
Далі потрібно сформувати декоратори для раніше створених функцій, написаних на С, зокрема normalize, isequal, dot, а також ввести допоміжні функції для конвертації матриці parseMatrix та списків makelist та parseList.
power_iteration.c
// Convert C array to Python list
static PyObject *makelist(float array[], size_t size) {
PyObject *lst = PyList_New(size);
for (size_t i = 0; i != size; ++i) {
PyList_SET_ITEM(lst, i, PyFloat_FromDouble(array[i]));
}
return lst;
}
// Convert Python list to C array
static float *parseList(PyObject* self, PyObject *initial_vector) {
float *vector = (float*) malloc(N * sizeof(float));
if (PyObject_Length(initial_vector) != N) {
return NULL;
}
for (int i = 0; i < N; i++) {
PyFloatObject *item = PyList_GetItem(initial_vector, i);
float value = PyFloat_AS_DOUBLE(item);
vector[i] = value;
}
return vector;
}
// Convert Python matrix (list with list as elements) to appropriate C array
static float *parseMatrix(PyObject* self, PyObject *initial_matrix) {
float *matrix = (float*) malloc(N * N * sizeof(float));
if (PyObject_Length(initial_matrix) != N) {
printf("error, invalid length of matrix\
");
return NULL;
}
for (int i = 0; i < N; i++) {
PyFloatObject *row = PyList_GetItem(initial_matrix, i);
for (int j = 0; j < N; j++) {
PyFloatObject *item = PyList_GetItem(row, j);
float value = PyFloat_AS_DOUBLE(item);
*(matrix + i*N + j) = value;
}
}
return matrix;
}
static PyObject* normalize(PyObject* self, PyObject* args) {
PyObject *i_vector;
if (!PyArg_ParseTuple(args, "O", &i_vector)) {
printf("error, when parsing arguments\
");
return NULL;
}
float *vector = parseList(self, i_vector);
float *normalized_vector = c_normalize(vector);
free(vector);
return makelist(normalized_vector, N);
}
static PyObject* is_equal(PyObject* self, PyObject* args) {
PyObject *i_first_vector, *i_second_vector;
if (!PyArg_ParseTuple(args, "OO", &i_first_vector, &i_second_vector)) {
printf("error, when parsing arguments\
");
return NULL;
}
float *first_vector = parseList(self, i_first_vector);
float *second_vector = parseList(self, i_second_vector);
bool _is_equal = c_is_equal(first_vector, second_vector);
free(first_vector);
free(second_vector);
if (_is_equal) {
Py_RETURN_TRUE;
} else {
Py_RETURN_FALSE;
}
}
static PyObject* dot(PyObject* self, PyObject* args) {
PyObject *i_matrix, *i_vector;
if (!PyArg_ParseTuple(args, "OO", &i_matrix, &i_vector)) {
printf("error, when parsing arguments\
");
return NULL;
}
float *matrix = parseMatrix(self, i_matrix);
float *vector = parseList(self, i_vector);
float *r_vector = c_dot(matrix, vector);
free(matrix);
free(vector);
return makelist(r_vector, N);
}
Залишається лише формальна частина. Робимо мапінг функцій та їхніх назв у майбутній бібліотеці Python та додаємо певний опис модуля:
power_iteration.c
static PyMethodDef myMethods[] = {
{ "normalize", normalize, METH_VARARGS, "Normalize vector" },
{ "is_equal", is_equal, METH_VARARGS, "Check if two lists are equal" },
{ "dot", dot, METH_VARARGS, "Performs multiplication matrix and array" },
{ NULL, NULL, 0, NULL }
};
static struct PyModuleDef power_iteration = {
PyModuleDef_HEAD_INIT, "power_iteration", "Test C extension for power iteration", -1, myMethods
};
PyMODINIT_FUNC PyInit_power_iteration(void) {
return PyModule_Create(&power_iteration);
}
Тепер необхідно сформувати файл setup.py, для того щоб новостворену бібліотеку можна було встановити:
setup.py
from distutils.core import setup, Extension
setup(name='power_iteration', version='1.0', ext_modules=[Extension('power_iteration', ['power_iteration.c'])])
Після встановлення залишається імпортувати модуль power_iteration та використати його методи, вже без жодних конвертацій типів даних і додаткових залежностей.
cextension_python.py
import power_iteration
eigenvector = INITIAL_GUESS
while not power_iteration.is_equal(eigenvector, prev_value):
prev_value = eigenvector
eigenvector = power_iteration.dot(TEST_MATRIX, eigenvector)
eigenvector = power_iteration.normalize(eigenvector)
Тож вся складність інтеграції зводиться до формування пакету power_iteration, який далі можна імпортувати у скрипт.
Тут необхідно зробити певний відступ. Дійсно, використання C extension вимагає значно більш додаткових рядків коду (на С, та ще й з використанням неочевидних, на перший погляд, макросів і типів даних з Python.h). Та робити це не так складно, оскільки все досить однотипно і реалізуються шляхом копіювання та вставки вже написаного. Нічого нового вигадувати не доводиться. В результаті реалізація підходу з C extension була значно швидшою, ніж аналогічний підхід з використанням ctypes.
benchmarks
Перейдемо до порівняння часу виконання алгоритму power iteration, використовуючи різні підходи.
На цій діаграмі наведений час пошуку власного значення (в секундах) для матриці розмірності 1000x1000. Розрахунки виконувались на Intel Pentium B960 (2.2 GHZ). Для основних випадків назви стовпчиків відповідають назві файлу, в якому він реалізовувався. Але для зручності наведу опис та назви підходів тут:
- pure_python — реалізація на чистому Python, без додаткових модулів та розширень.
- numpy — реалізація з використанням numpy. Крім типу даних array, використовувались також методи модуля: dot, allclose та поелементне ділення для нормалізації.
- pure_c — реалізація на С, без інтеграції з Python.
- ctypes — реалізація з інтеграцією коду написаного на С в Python, використовуючи бібліотеку ctypes.
- c_extension — реалізація з інтеграцією коду написаного на С в Python із застосуванням С extension.
- ctypes_full — реалізація з інтеграцією коду написаного на С в Python, використовуючи бібліотеку ctypes. Основний цикл був винесений в код на С, а в Python отримується лише результат — власний вектор матриці.
- c_extension_full — аналогічне до ctypes_full, але з інтеграцією за допомогою С extension.
Важливе зауваження! pure_c містить дуже неоптимальний імпорт тестової матриці, тому не рекомендується приймати час виконання pure_c за абсолютний та придатний для порівняння з іншими підходами.
Висновки
Зазначимо, що усі подальші твердження є суб'єктивною думкою автора і базуються лише на даних, вказаних вище. Отже:
- Складність інтеграції коду на С в Python за допомогою ctypes полягає у нетривіальній конвертації складних типів даних.
- Щоб інтегрувати код на С в Python з використанням C extension, потрібно вносити зміни мовою С, а це не дуже зручно для Python-розробників.
- Отриманий в результаті підходу C extension пакет є зручним для повторних використань і не вимагає особливих знань про структуру коду С.
- Отриманий алгоритм на C extension працює майже в 38 разів швидше, ніж його аналог на ctypes. Швидкість зумовлюється часом конвертації типів даних. В першому випадку це відбувається засобами С, а в другому — Python.
- Значно ефективніше винести цілий модуль в С та не звертатись багаторазово до окремих його методів. Тоді час не витрачатиметься на постійну конвертацію між типами даних.
Ще немає коментарів