Cython: самый быстрый способ присвоить точечные индексы грубой сетке - PullRequest
0 голосов
/ 12 февраля 2019

Я пытаюсь разработать симуляцию NBody для бедного человека, и я бы хотел ускорить процесс взаимодействия, рассматривая только моих ближайших соседей.Кажется, что самый простой способ сделать это состоит в том, чтобы связать все координаты точек с грубой сеткой, а затем зациклить точки, которые содержатся в ближайших соседних точках сетки.

Изначально я собирался сделать quadTreeдля запросов к ближайшим соседям, но данные довольно динамичны от итерации к итерации, и я ожидаю обработки симуляций, которые имеют более миллиона координатных точек, поэтому мне показалось, что решение на основе сетки было бы лучшим местом для начала.

В настоящее время у меня есть следующий рабочий код:

import numpy as np
cimport numpy as np
from cython.view cimport array as cvarray
cimport cython
from libc.math cimport sinh, cosh, sin, cos, acos, exp, sqrt, fabs, M_PI, floor
from collections import defaultdict

DTYPE = np.float64

ctypedef np.float64_t DTYPE_t
cdef DTYPE_t pi = 3.141592653589793


# Define Grid Structure
@cython.cdivision(True)
@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False)  # turn off negative index wrapping for entire function
def gridStructure(DTYPE_t[:,:,:] X, Py_ssize_t ref, int xSize, int ySize):
    Grid = defaultdict(int)
    cdef Py_ssize_t ii,jj
    cdef Py_ssize_t N = X.shape[1]
    cdef DTYPE_t[:] maxes = np.max(X[ref,:,:],axis=0)
    cdef DTYPE_t[:] mines = np.min(X[ref,:,:],axis=0)
    cdef Py_ssize_t pointX, pointY, index
    for ii in range(0,N):
        pointX = int(floor((X[ref,ii,0]-mines[0])/maxes[0]*xSize))
        pointY = int(floor((X[ref,ii,1]-mines[1])/maxes[1]*ySize))
        index = pointX + xSize*pointY
        if index in Grid.keys():
            Grid[index].append(ii)
        else:
            Grid[index] = [ii]
    return Grid

Я использую стандартную структуру dict для хранения ключа (целочисленного индекса сетки) и набора целых чисел, которые указывают на координаты впросмотр памяти.Можно с уверенностью предположить, что сетка редкая;координаты обычно имеют неоднородное распределение.Я думаю, что моя структура достаточно проста, поэтому я должен использовать структуру на основе C вместо dict, чтобы я мог избавиться от оболочек Python.

В конце концов, эта функция не будет вызываться пользователем напрямую, поэтому яна этом уровне не требуется никаких взаимодействий с питоном.

В настоящее время код может скопировать 1000000 точек на сетку 256x256 примерно за секунду, возможно ли выполнить эту операцию быстрее для потенциальной симуляции бедняков в реальном времени?Любой совет для правильной структуры данных для ускорения этой функции будет принята с благодарностью.Я в основном хочу вернуть некоторую структуру данных, где я могу вызвать целочисленный ключ и вернуть набор точек, связанных с этим ключом.Спасибо!

1 Ответ

0 голосов
/ 15 февраля 2019

ОК, у меня есть лучшее решение.Я нашел этот пост здесь: Эффективная (и хорошо объясненная) реализация Quadtree для двумерного обнаружения столкновений

, которая говорит о том, как выполнить логику на сетке в начале его второго ответа.,Поскольку сетка имеет фиксированный размер и количество элементов в моем симе фиксировано, я получаю ускорение в 20 раз, если я просто предварительно выделяю память.Вот новая функция сетки:

import numpy as np
cimport numpy as np

cimport cython
from libc.math cimport floor

DTYPE = np.float64
ctypedef np.float64_t DTYPE_t

@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
def gridStructure(DTYPE_t[:,:,:] X, np.int_t ref, np.int_t xSize, np.int_t[:] grid, np.int_t[:] eList):
    cdef Py_ssize_t N = X.shape[1]
    cdef np.int_t[:] lastEl = np.zeros(xSize*xSize, dtype=np.int)
    cdef DTYPE_t[:] maxes = np.max(X[ref,:,:],axis=0)
    cdef DTYPE_t[:] mines = np.min(X[ref,:,:],axis=0)
    cdef DTYPE_t widthX = max(maxes[0]-mines[0],maxes[1]-mines[1])
    cdef DTYPE_t ratio = xSize/widthX
    cdef np.int_t index,pointX,pointY
    cdef np.int_t ii,zz

    for ii in range(0,N):
        pointX = int(floor((X[ref,ii,0]-mines[0])*ratio*0.99))
        pointY = int(floor((X[ref,ii,1]-mines[1])*ratio*0.99))
        index = int(pointX + xSize*pointY)
        if grid[index] == -1:
            grid[index] = ii
            lastEl[index] = ii
        else:
            zz = lastEl[index]
            eList[zz] = ii
            lastEl[index] = ii

    return 0

1 000 000 точек в сетке 256x256 за 50 мс.Я могу жить с этим сейчас.

Код предполагает, что элементы для сетки и списка элементов (eList) заполняются с -1 при вводе.

РЕДАКТИРОВАТЬ: я удалил оригинальный ответпотому что предыдущее решение, которое я разместил, было неверным.Мне пришлось добавить дополнительный буфер памяти для хранения последней позиции указателя для каждого индекса.

...