У меня есть строка кода, и я удивлен, что версия, понятая списком, работает медленнее. Я хочу рассчитать «средние расстояния сферы».
Шаги, которые нужно было сделать:
Я использую разреженную матрицу, в которой хранится отношение смежности между узлами. Поскольку это лапласианы, я удалил диагональ, но это не имеет значения. Я использую scipy.shortest_path (или теперь фактически scipy.dijstra) для расчета кратчайшего пути между двумя узлами. Возвращаемое значение на самом деле представляет собой массив для узла «I» для всех элементов до предела (если он добавлен в качестве аргумента для dijkstra). Поскольку я использую подмножество элементов из матрицы, мне понадобится словарь, который может переводить между элементом и индексом в списке кратчайших путей.
Итак, код:
Я импортирую библиотеки:
import numpy as np
import time
import random
from random import randint
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import breadth_first_tree
import scipy as scp
from scipy.stats import uniform
from scipy import io
from scipy import sparse # package for sparse matrices
from scipy import linalg # package for linear algebra
from scipy.sparse import identity
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import breadth_first_order
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import shortest_path,dijkstra
Чтение разреженной матрицы с некоторыми параметрами (nxn):
input_m = "L-4-4.0-0.02-4-2-1.mtx"
L = scp.sparse.csc_matrix(scp.io.mmread(input_m), dtype=int)
ID = identity(np.shape(L)[0], dtype='int8', format='dia')
WA = abs(5*ID - L)
Определите функции, которые возвращают ракушки (все до расстояния n), шар, содержащий все элементы, и еще одну, которая возвращает только ракушку на расстоянии n:
def GetShellBall(WA,n,ind):
p0 = np.zeros(np.shape(L)[0])
p0[ind] = 1
newp = p0
ball = []
shell = []
ball.append(ind)
shell.append([ind])
for it in range(n):
newp = WA@newp
for it2 in np.where(newp)[0]:
if it2 in ball:
newp[it2] = 0
else:
ball.append(it2)
shell.append(np.where(newp)[0])
return ball,shell
def GetShellj(WA,n,ind):
p0 = np.zeros(np.shape(L)[0])
p0[ind] = 1
newp = p0
ball = []
shell = []
ball.append(ind)
shell.append([ind])
for it in range(n):
newp = WA@newp
for it2 in np.where(newp)[0]:
if it2 in ball:
newp[it2] = 0
else:
ball.append(it2)
shell.append(np.where(newp)[0])
return shell[n]
Создайте оболочку и шар вокруг узла на заданном ограниченном расстоянии, используя алгоритм Дейкстры и словарь для сопоставления элементов в шаре и индексов в шпатах (список кратчайших путей).
%%time
it = 0
N = 11
ball,shell = GetShellBall(WA,N,it) #shell and ball around it at distance N
dict_ba = dict(zip(ball,arange(len(ball)))) #dictionary: ball_element - number_i (list in ball)
dict_ab = dict(zip(arange(len(ball)),ball)) #dictionary: num_i (list in ball) - ball element
spaths = np.asarray([dijkstra(csgraph=WA, directed=True,
limit = N,unweighted = True, indices=it,
return_predecessors=False) for it in ball])
Время ЦП: пользователь 7,61 с, система: 11,9 мс, всего: 7,62 с Время стены: 7,63 с
ТЕПЕРЬ начинается «тяжелая часть»: вычислить средние расстояния вокруг всех узлов на расстоянии d от моего исходного узла (так что узел el_i имеет shell_at_dist_i, получить оболочки вокруг всех элементов shell_at_dist_i с радиусом «i» и вычислить среднее расстояние вокруг двух оболочки, затем усредняем по нему, для этого кода циклическая версия быстрее, чем векторизованная:*
%%timeit
#for i in range(len(shell)):
DD = []
for i in range(int(floor(N/2))):
sumd = []
chosen_paths = spaths[[dict_ba[it] for it in shell[i]]]
for eli in shell[i]:
shellj = GetShellj(WA,i,eli)
sumd = [[chp[elj] for elj in shellj] for chp in chosen_paths]
DD.append(mean(sumd))
36,2 мс ± 943 мкс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 10 циклов в каждом)
%%timeit
#for i in range(len(shell)):
DD = []
for i in range(int(floor(N/2))):
sumd = []
chosen_paths = spaths[[dict_ba[it] for it in shell[i]]]
for eli in shell[i]:
sumd = [[chp[elj] for elj in GetShellj(WA,i,eli)] for chp in chosen_paths]
DD.append(mean(sumd))
422 мс ± 7,05 мс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 1 циклу в каждом)
%%timeit
#for i in range(len(shell)):
DD = []
for i in range(int(floor(N/2))):
chosen_paths = spaths[[dict_ba[it] for it in shell[i]]]
sumd = [[[chp[elj] for elj in GetShellj(WA,i,eli)] for chp in chosen_paths] for eli in shell[i]]
DD.append(mean(sumd))
439 мс ± 5,85 мс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 1 циклу в каждом)
%%timeit
#for i in range(len(shell)):
DD = []
for i in range(int(floor(N/2))):
sumd = [
[
[
chp[elj] for elj in GetShellj(WA,i,eli)
] for chp in spaths[[dict_ba[it] for it in shell[i]]]
] for eli in shell[i]]
DD.append(mean(sumd))
445 мс ± 2,68 мс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 1 циклу в каждом)
%%timeit # для i в диапазоне (len (оболочка)):
DD = [mean([
[
[
chp[elj] for elj in GetShellj(WA,i,eli)
] for chp in spaths[[dict_ba[it] for it in shell[i]]]
] for eli in shell[i]]) for i in range(int(floor(N/2)))]
440 мс ± 8,05 мс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 1 циклу в каждом)
Почему циклическая версия намного быстрее? Я хотел бы масштабировать его. Это всего лишь одна отправная точка от матрицы 2300x2300, но мне нужно пойти намного выше (~ полмиллиона х полмиллиона матриц).
Может быть, есть методы, чтобы ускорить его еще больше? (помимо использования C вместо python....)
Обновлено: Решено: ошибка.
Остается вопрос: есть ли более быстрый способ?
Во многих фрагментах кода, которые вы опубликовали, есть ошибка, из-за которой большая часть работы внутреннего цикла игнорируется, а среднее значение вычисляется только с использованием значений из последней итерации.
[[chp[elj] for elj in shellj] for ...]
не эквивалентен [[chp[elj] for elj in GetShellj(WA,i,eli)] for ...]
, последний будет вызывать GetShellj
один раз каждый раз во внешнем цикле, что приводит к неэффективности.
И выигрывают те, кто указывает на ошибку. (спасибо и извините) Но тогда мой вопрос: не должно ли быть быстрее написать цикл for в виде []? (или это полностью эквивалентно?)
@Kregnach Это может быть немного быстрее, но в целом незначительно.
Первая версия быстрее остальных, потому что вы не вызываете GetShellj(WA,i,eli)
столько раз.
Если вы прилепите декоратор @cache или @lru_cache к этой функции, все подходы должны быть примерно такими же быстрыми.
Кэширование потребует хэшируемых входных данных.
Спасибо! Это действительно большая вина. Не могли бы вы уточнить @cache и другие, как я могу украсить функцию? Я никогда не использовал их раньше.
@Kregnach Нет, просто перейдите к документации, чтобы увидеть примеры использования.
Вы правы, справедливо!
"чем векторизованный" - какой векторизованный?