Улучшение времени выполнения симуляции MonteCarlo с использованием numpy - PullRequest
Новая
       17

Улучшение времени выполнения симуляции MonteCarlo с использованием numpy

0 голосов
/ 24 марта 2019

Я давно работаю над симуляцией Монте-Карло на двумерной модели Изинга. У меня проблема чисто вычислительная и не требует каких-либо знаний о модели Изинга.

Мне удалось полностью реализовать модель. Тем не менее, код работает недостаточно быстро, так как я хочу моделировать довольно большие системы. Запустив подходящее количество циклов Монте-Карло, время выполнения составляет около получаса с использованием самой большой системы, которую я пробовал.

Основная причина высокой продолжительности выполнения состоит в том, что я не могу поддерживать выполнение функций в "пустой форме", мне приходится использовать обычный цикл for Python, который, на мой взгляд, является самой большой проблемой. Две функции, с которыми мне нужна помощь, следующие: 

def helix(S, i, N, L, beta, r, B = 0, J_v = 1, J_p = 1):
    nn = i + 1
    if nn >= N:
        nn -= N
    sum_nn = S[nn]
    nn = i - 1
    if nn < 0:
        nn += N
    sum_nn += S[nn]
    nn = i + L
    if nn >= N:
        nn -= N
    sum_nn += S[nn]
    nn = i - L
    if nn < 0:
        nn += N
    sum_nn += S[nn]
    del_E = 2 * S[i] * sum_nn
    if del_E <= 0:
        S[i] = - S[i]
        return del_E, 2 * S[i] / N
    elif np.exp(-beta * sum_nn) > r:
        S[i] = - S[i]
        return del_E, 2 * S[i] / N
    else:
        return 0, 0


def random_sweep(S, N, L, beta, B = 0, J_v = 1, J_p = 1):
    del_m_sweep = 0
    del_He_sweep = 0
    rand_list_index = np.random.randint(N, size=N)
    rand_numbers = np.random.rand(N)
    for i in range(N):
        del_E, del_m = helix(S, rand_list_index[i], N, L, beta, 
                              rand_numbers[i], B, J_v, J_p)
        del_He_sweep += del_E
        del_m_sweep += del_m
    return del_He_sweep, del_m_sweep

Я хочу, чтобы все было в форме, так как я понимаю, что это очень поможет с вычислительной скоростью. В основном я хочу заменить цикл for в random_sweep, просто вызвав функцию helix(S, rand_list_index, N, L, beta, rand_numbers, B, J_v, J_p), и с использованием гибкой гибкости он автоматически вызывает функцию спирали N раз (на основе аргументов rand_list_index и rand_numbers вместо rand_list_index[i] и rand_numbers[i]). Причина, по которой я сейчас не могу этого сделать, заключается в том, что мне нужно обновить массив S с индексом i. Есть ли способ обойти это обновление?

Я уже давно пытаюсь это реализовать, поэтому буду очень признателен за любую помощь!

...