Введение
Вероятно, лучший способ понять это - создать итератор в Python и поэкспериментировать с ним. Это будет медленно, но это подтвердит, что вы делаете правильно. Затем вы используете NpyIter_AdvancedNew, используя параметры по умолчанию везде, где это возможно.
Боюсь, я сам не перевел это на C-код - это заняло у меня слишком много времени. Поэтому я предлагаю вам не принимать этот ответ, поскольку он действительно дает только отправную точку.
Я полагаю, что любые улучшения производительности будут разочаровывающими, учитывая количество усилий, которые написание кода на C (особенно если учесть, что написание быстрого кода требует более глубокого понимания). В конце ответа я предлагаю несколько более простых альтернатив, которые я бы порекомендовал вместо использования C API.
Примеры
Я перевел пару строк из вашего кода в качестве примеров:
d = p[:, b:b+1] - p[:, :b]
Становится
with np.nditer([p[:,b],p[:,:b],None],
op_axes=[[0,-1],[0,1],[0,1]]) as it:
for x,y,z in it:
z[...] = x - y
d = it.operands[2]
Обратите внимание, что вам нужно предварительно нарезать массив p. Я передал один из массивов как None. Это переводит указатель NULL в C и означает, что массив создается с соответствующим размером (с использованием стандартных правил вещания).
В терминах op_axes первый массив только 1D, поэтому я сказал: «сначала переберите ось 0; оси 1 нет». Второй и третий массивы - это два D, поэтому я сказал «итерация по оси 0, затем по оси 1».
В Python он выводит op_flags автоматически. Я не знаю, будет ли это делать в C. Если нет, то они должны быть:
npy_uint32 op_flags[] = { NPY_ITER_READONLY,
NPY_ITER_READONLY,
NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE };
Самый важный бит в том, что третья ось выделена.
На мой взгляд, вы хотите указать op_dtypes в C как
{ PyArray_DescrFromType(NPY_DOUBLE),
PyArray_DescrFromType(NPY_DOUBLE),
NULL }
, чтобы заставить массивы иметь правильный тип (тип третьего выделенного массива может быть получен из двух входов). Таким образом, вы сможете использовать указатели данных на double* в C.
Строка:
d2 = (d*d).sum(axis=0)
переводится как
with np.nditer([d,None],
flags=['reduce_ok'],
op_flags=[['readonly'],['readwrite','allocate']],
op_axes=[[1,0],[0,-1]]) as it:
it.operands[1][...] = 0
for x,z in it:
z[...] += x*x
d2 = it.operands[1]
Самым важным отличием является то, что это уменьшение (второй выходной массив меньше входного, потому что одна из осей является суммой). Поэтому мы передаем 'limit_ok' как флаг.
Второй массив имеет только одну ось, поэтому op_axes равен [0, -1]. Ось второго массива соответствует оси 1 первого массива, поэтому op_axes для первого массива установлено на [1, 0].
При переводе на C строка it.operands[1][...] = 0 усложняется:
Обратите внимание, что если вы хотите сделать сокращение для автоматически распределенного вывода, вы должны использовать NpyIter_GetOperandArray, чтобы получить его ссылку, а затем установить каждое значение в единицу сокращения перед выполнением цикла итерации.
В C я, вероятно, сначала выделю d2 как нулевой массив и передам его итератору.
Альтернативы
Написание кода C API для этого требует большого количества кода, проверки ошибок, подсчета ссылок и т. Д. Хотя это должен быть «простой» перевод (API nditer в основном одинаков как в C, так и в Python), это не так. не легко.
Если бы вы были знакомы с использованием некоторых стандартных инструментов для ускорения Python, например Numba, NumExpr или Cython. Numba и NumExpr - это компиляторы, работающие точно по времени, которые могут делать такие вещи, как избегание размещения промежуточных массивов Cython - это "Python-подобный" язык, в котором вы можете указывать типы. Чтобы показать первые несколько частей, переведенных на Cython:
def grav3(double[:,:] p, M):
G = 6.67408e-2 # m³/s²T
cdef int l = p.shape[1]
a = np.empty(shape=(2, l))
a[:, 0] = 0
cdef double[:,::1] d
cdef double[::1] d2
cdef Py_ssize_t b, i, j
for b in range(1, l):
# computing the distance between body #b and all previous
d = np.empty((p.shape[0],b))
for i in range(d.shape[0]):
for j in range(b):
d[i,j] = p[i,b] - p[i,j]
d2 = np.zeros((d.shape[1]))
for j in range(d.shape[1]):
for i in range(d.shape[0]):
d2[j] += d[i,j]*d[i,j]
if d2[j] == 0:
d2[j] = 1
Здесь я указал некоторые массивы - 1D или 2D двойные массивы double[:] или double[:,:]. Затем я написал циклы явно, что позволяет избежать создания промежуточных звеньев.
Cython генерирует код C, где получает PyArray_DATA, а затем использует PyArray_STRIDES, чтобы определить, куда обращаться в двумерном массиве. Вы можете найти это проще, чем использовать итераторы. Вы можете изучить код, сгенерированный Cython, чтобы увидеть, как он это делает. В Numpy также доступны PyArray_GetPtr функции , которые предоставляют такой доступ, который может оказаться проще, чем использование итераторов.