Я изучаю, как кодировать на Python, и я пытаюсь воссоздать код, который я делал в колледже.
Код основан на двухмерной модели Изинга, применяемой к эпидемиологии.Что он делает:
- , он строит массив 2D 100x100, используя numpy, и присваивает значение -1 каждому элементу.
- Энергия рассчитывается на основе функции
calc_h в приведенном ниже сценарии. - Затем код случайным образом выбирает ячейку из решетки, изменяет значение на 1, а затем снова вычисляет энергию системы.
- Затем код сравнивается, если энергия системы меньше или равна предыдущей конфигурации.Если это так, он «принимает» изменение.Если это не так, вероятность сравнивается со случайным числом, чтобы определить, принято ли изменение.Эта часть выполняется в функции
metropolis. - Код повторяет процесс, используя цикл while до максимальной указанной итерации,
max_iterations.-Код подсчитывает количество элементов со значением -1 (которое является переменной s) и количество элементов со значением 1 (которое является переменной i) в функции countSI.Сценарий добавляет к текстовому файлу каждые 500 итераций.
ПРОБЛЕМА
Я запустил сценарий и, кроме того, что для его выполнения потребовалось слишком много времени, подсчет останавливается на15500. Код не выдает никакой ошибки, он просто продолжает работать.Я ждал около 3 часов, чтобы он закончился, но он все еще идет только до 15500 итераций.
Я попытался закомментировать запись в блок csv и вместо этого сначала распечатать значения, чтобы увидеть, как это происходит,и там я вижу, он снова останавливается на 15500.
Я понятия не имею, что не так, поскольку он не выдает никакой ошибки, и код не останавливается.
Вот и всескрипт.Я помещаю описание того, что часть делает ниже каждого блока:
import numpy as np
import random as r
import math as m
import csv
init_size = input("Input array size: ")
size = int(init_size)
эта часть инициализирует размер двумерного массива.В целях наблюдения я выбрал латекс 100 на 100.
def check_overflow(indx, size):
if indx == size - 1:
return -indx
else:
return 1
Я использую эту функцию для функции calc_h, чтобы инициализировать условие круговой границы.Проще говоря, края решетки связаны друг с другом.
def calc_h(pop, J1, size):
h_sum = 0
r = 0
c = 0
while r < size:
buffr = check_overflow(r, size)
while c < size:
buffc = check_overflow(c, size)
h_sum = h_sum + J1*pop[r,c] * pop[r,c-1] * pop[r-1,c] * pop[r+buffr,c] * pop[r,c+buffc]
c = c + 1
c = 0
r = r + 1
return h_sum
эта функция вычисляет энергию системы, беря сумму произведений на значение ячейки, его верхние, нижние, левые и правые соседи, умноженные на константу J.
def metropolis(h, h0, T_):
if h <= h0:
return 1
else:
rand_h = r.random()
p = m.exp(-(h - h0)/T_)
if rand_h <= p:
return 1
else:
return 0
Это определяет, будет ли принято изменение от -1 до 1 в зависимости от того, что calc_h получает.
def countSI(pop, sz, iter):
s = np.count_nonzero(pop == -1)
i = np.count_nonzero(pop == 1)
row = [iter, s, i]
with open('data.txt', 'a') as si_csv:
tally_data = csv.writer(si_csv)
tally_data.writerow(row)
si_csv.seek(0)
Эта часть подсчитывает количество -1 и 1 в решетке.
def main():
J = 1
T = 4.0
max_iterations = 150000
population = np.full((size, size), -1, np.int8) #initialize population array
2Dмассив инициализируется в population.
h_0 = calc_h(population, J, size)
turn = 1
while turn <= max_iterations:
inf_x = r.randint(1,size) - 1
inf_y = r.randint(1,size) - 1
while population[inf_x,inf_y] == 1:
inf_x = r.randint(1,size) - 1
inf_y = r.randint(1,size) - 1
population[inf_x, inf_y] = 1
h = calc_h(population, J, size)
accept_i = metropolis(h,h_0,T)
Это основной цикл, в котором выбирается случайная ячейка, и определяется ли изменение или нет, определяется функцией metropolis.
if (accept_i == 0):
population[inf_x, inf_y] = -1
if turn%500 == 0 :
countSI(population, size, turn)
Сценарий рассчитывает каждую 500-ю итерацию.
turn = turn + 1
h_0 = h
main()
Ожидаемый вывод - текстовый файл с подсчетом числаs и i каждую 500-ю итерацию.что-то похожее на это:
500,9736,264
1000,9472,528
1500,9197,803
2000,8913,1087
2500,8611,1389
3000,8292,1708
3500,7968,2032
4000,7643,2357
4500,7312,2688
5000,6960,3040
5500,6613,3387
6000,6257,3743
6500,5913,4087
7000,5570,4430
7500,5212,4788
Я понятия не имею, с чего начать.Сначала я подумал, что проблема связана с записью в csv, но исследование функции печати доказывает обратное.Я постарался сделать его максимально лаконичным.
Надеюсь, вы, ребята, сможете помочь!Я действительно хочу выучить этот язык и начать симулировать множество вещей, и я думаю, что этот мини-проект - отличный стартовый шаг для меня.
Большое спасибо!