2021年美赛A题真菌元胞自动机Python实现

import matplotlib.pyplot as plt
import random
import numpy as np
import matplotlib.animation as animationdef save_fungi_ca_gif():  # save the gif file to the pathtarget_gif_path = "E:/engineering space/figure/gif format/"target_gif_name = "aggressive_ca.gif"target_gif_full = target_gif_path + target_gif_nameanim_1.save(target_gif_full, writer='pillow')  # save the animationprint('Saved')def calculate_random_neighbour(stay, left, up, right, down, no_reproduce):  # get the relative location of the newlytotal_rate = stay + left + up + right + down + no_reproduce  # reproduced fungi cellif total_rate == 0.0:total_rate = 1cap1 = stay / total_rate  # cap 1 to 5 is the boundary of probability areacap2 = cap1 + left / total_ratecap3 = cap2 + up / total_ratecap4 = cap3 + right / total_ratecap5 = cap4 + down / total_raternd = random.random()if 0 <= rnd < cap1:return 0elif cap1 <= rnd < cap2:  # divide the range into 6 partsreturn 1elif cap2 <= rnd < cap3:  # randomly select one partreturn 2elif cap3 <= rnd < cap4:return 3elif cap4 <= rnd < cap5:return 4else:return 5def get_extension_parameters(i_get_extension_parameters):  # get the extension parametersextension_parameters = extension_rate_list[fungi_list[i_get_extension_parameters][3]]return extension_parametersdef get_random_neighbour(i_get_random_neighbour):  # calculate the density of neighbour mesh as well as self meshextension_rate_get_random_neighbour = get_extension_parameters(i_get_random_neighbour)stay_density = my_board[fungi_list[i_get_random_neighbour][0], fungi_list[i_get_random_neighbour][1]] / one_mesh_maxstay_rate = extension_rate_get_random_neighbour - stay_densityif fungi_list[i_get_random_neighbour][1] != 0:left_density = my_board[fungi_list[i_get_random_neighbour][0], fungi_list[i_get_random_neighbour][1] - 1] \/ one_mesh_maxleft_rate = extension_rate_get_random_neighbour * (1 - left_density + 0.1)  # further calculate the transferring rate of# the new fungielse:left_density = 0.5left_rate = 0if fungi_list[i_get_random_neighbour][0] != 0:up_density = my_board[fungi_list[i_get_random_neighbour][0] - 1, fungi_list[i_get_random_neighbour][1]] \/ one_mesh_maxup_rate = extension_rate_get_random_neighbour * (1 - up_density + 0.1)else:up_density = 0.5up_rate = 0if fungi_list[i_get_random_neighbour][1] != patch_size - 1:right_density = my_board[fungi_list[i_get_random_neighbour][0], fungi_list[i_get_random_neighbour][1] + 1] \/ one_mesh_maxright_rate = extension_rate_get_random_neighbour * (1 - right_density + 0.1)else:right_density = 0.5right_rate = 0if fungi_list[i_get_random_neighbour][0] != patch_size - 1:down_density = my_board[fungi_list[i_get_random_neighbour][0] + 1, fungi_list[i_get_random_neighbour][1]] \/ one_mesh_maxdown_rate = extension_rate_get_random_neighbour * (1 - down_density + 0.1)else:down_density = 0.5down_rate = 0total_density = stay_density + left_density + up_density + right_density + down_densityno_reproduce_rate = total_density / 8neighbour_location = calculate_random_neighbour(stay_rate, left_rate, up_rate, right_rate, down_rate,no_reproduce_rate)return neighbour_locationdef reproduce_one_cell(i_reproduce_one_cell):  # reproduce a single fungi celliterations = 0while True:  # loop until the random neighbour meets the boundary conditionsiterations += 1# print('iteration:', iterations)reproduce_location = get_random_neighbour(i_reproduce_one_cell)# the followings are boundary conditionsif iterations < 3:if reproduce_location == 0:  # stay in one meshcurrent_mesh_num = my_board[fungi_list[i_reproduce_one_cell][0], fungi_list[i_reproduce_one_cell][1]]if current_mesh_num < one_mesh_max:# haven't arrived the maximum of one meshbreak  # this is the location wantedelse:# print(iterations)  # how many loops have been takencontinue  # random againif reproduce_location == 1:  # go to left meshif fungi_list[i_reproduce_one_cell][1] == 0:  # reaches the left boundary# print(iterations)continue  # random againelif my_board[fungi_list[i_reproduce_one_cell][0], fungi_list[i_reproduce_one_cell][1] - 1] \< one_mesh_max:break  # this is the location wantedelse:continueif reproduce_location == 2:  # go to up meshif fungi_list[i_reproduce_one_cell][0] == 0:  # reaches the up boundary# print(iterations)continue  # random againelif my_board[fungi_list[i_reproduce_one_cell][0] - 1, fungi_list[i_reproduce_one_cell][1]] \< one_mesh_max:break  # this is the location wantedelse:continueif reproduce_location == 3:  # go to right meshif fungi_list[i_reproduce_one_cell][1] == patch_size - 1:  # reaches the right boundary# print(iterations)continue  # random againelif my_board[fungi_list[i_reproduce_one_cell][0], fungi_list[i_reproduce_one_cell][1] + 1] \< one_mesh_max:break  # this is the location wantedelse:continueif reproduce_location == 4:  # go to down meshif fungi_list[i_reproduce_one_cell][0] == patch_size - 1:  # reaches the down boundarycontinue# print(iterations)elif my_board[fungi_list[i_reproduce_one_cell][0] + 1, fungi_list[i_reproduce_one_cell][1]] \< one_mesh_max:break  # this is the location wantedelse:continueelse:reproduce_location = 5breakglobal fungi_num# take reproducing activityif reproduce_location != 5:if reproduce_location == 0:fungi_list.append([fungi_list[i_reproduce_one_cell][0], fungi_list[i_reproduce_one_cell][1], 0, fungi_list[i_reproduce_one_cell][3]])# print('Stayed')if reproduce_location == 1:fungi_list.append([fungi_list[i_reproduce_one_cell][0], fungi_list[i_reproduce_one_cell][1] - 1, 0,fungi_list[i_reproduce_one_cell][3]])if reproduce_location == 2:fungi_list.append([fungi_list[i_reproduce_one_cell][0] - 1, fungi_list[i_reproduce_one_cell][1], 0,fungi_list[i_reproduce_one_cell][3]])if reproduce_location == 3:fungi_list.append([fungi_list[i_reproduce_one_cell][0], fungi_list[i_reproduce_one_cell][1] + 1, 0,fungi_list[i_reproduce_one_cell][3]])if reproduce_location == 4:fungi_list.append([fungi_list[i_reproduce_one_cell][0] + 1, fungi_list[i_reproduce_one_cell][1], 0,fungi_list[i_reproduce_one_cell][3]])fungi_list[i_reproduce_one_cell][2] += 1  # the reproduced time of a fungiif fungi_list[i_reproduce_one_cell][2] == reproduce_time_to_die_list[fungi_list[i_reproduce_one_cell][3]]:# when reproduced an exact times to diemy_board[fungi_list[i_reproduce_one_cell][0]][fungi_list[i_reproduce_one_cell][1]] -= 1  # remove fungi# in number boarddel fungi_list[i_reproduce_one_cell]  # the death of a fungifungi_num = len(fungi_list)update_board(fungi_num)  # get a board according to the new fungi reproducedelse:fungi_list[i_reproduce_one_cell][2] += 1  # the life span of a fungidef reproduce_round():  # reproduce all the current fungi cellsfor i_reproduce_round in range(fungi_num):reproduce_one_cell(i_reproduce_round)  # reproduce a single fungi celldef update_board(fungi_num_update_board):  # get a number board according to the new fungimy_board[fungi_list[fungi_num_update_board - 1][0]][fungi_list[fungi_num_update_board - 1][1]] += 1def get_rgb_board():paint_rgb_board = [[[1.0 for k in range(3)] for j in range(patch_size)] for i in range(patch_size)]  # initializeglobal red_fungi_numglobal green_fungi_numglobal blue_fungi_numred_fungi_num = 0green_fungi_num = 0blue_fungi_num = 0# color boardfor i_paint_rgb_board in range(fungi_num):if fungi_list[i_paint_rgb_board][3] == 0:red_fungi_num += 1paint_rgb_board[fungi_list[i_paint_rgb_board][0]][fungi_list[i_paint_rgb_board][1]][1] -= (0.99/ one_mesh_max)paint_rgb_board[fungi_list[i_paint_rgb_board][0]][fungi_list[i_paint_rgb_board][1]][2] -= (0.99/ one_mesh_max)elif fungi_list[i_paint_rgb_board][3] == 1:green_fungi_num += 1paint_rgb_board[fungi_list[i_paint_rgb_board][0]][fungi_list[i_paint_rgb_board][1]][0] -= (0.99/ one_mesh_max)paint_rgb_board[fungi_list[i_paint_rgb_board][0]][fungi_list[i_paint_rgb_board][1]][2] -= (0.99/ one_mesh_max)else:blue_fungi_num += 1paint_rgb_board[fungi_list[i_paint_rgb_board][0]][fungi_list[i_paint_rgb_board][1]][0] -= (0.99/ one_mesh_max)paint_rgb_board[fungi_list[i_paint_rgb_board][0]][fungi_list[i_paint_rgb_board][1]][1] -= (0.99/ one_mesh_max)return paint_rgb_boarddef attack_round():global fungi_numfor i_attack_round in range(fungi_num):if i_attack_round == fungi_num - 1:breakif fungi_list[i_attack_round][2] >= reproduce_time_to_die_list[fungi_list[i_attack_round][3]]:# when reproduced an exact times to diemy_board[fungi_list[i_attack_round][0]][fungi_list[i_attack_round][1]] -= 1  # remove fungi in number boarddel fungi_list[i_attack_round]  # the death of a fungifungi_num = len(fungi_list)# update_board(fungi_num)  # get a board according to the new fungi reproducedelse:continuedef update():  # update the board information and fungi list in a roundreproduce_round()  # reproduce all the current fungi cellsrgb_board_update = get_rgb_board()  # get a color board after a round of reproductionattack_round()return rgb_board_updatedef calculate_initial_position():  # calculate the initial position of cells for fair competitionradius_calculate_initial_position = patch_size / 6  # distance between initial cellsmiddle_point_initial_position = patch_size / 2for i_calculate_initial_position in range(initial_fungi_num):rotation_degree = i_calculate_initial_position / initial_fungi_num * 2 * np.pix_initial_position = middle_point_initial_position + np.cos(rotation_degree) * radius_calculate_initial_positiony_initial_position = middle_point_initial_position + np.sin(rotation_degree) * radius_calculate_initial_positionx_initial_list.append(int(round(x_initial_position)))y_initial_list.append(int(round(y_initial_position)))def animate_func(frame_sequence):  # depict the animationax_1.imshow(update())update_extension_rate()# ax_1.text(patch_size / 12, patch_size / 12, frame_sequence, color='white', size=20)print('total:', fungi_num, 'red:', red_fungi_num, 'green:', green_fungi_num, 'blue:', blue_fungi_num)print('frame number:', frame_sequence + 1)def update_extension_rate():global extension_rate_listred_fungi_rate = red_fungi_num / fungi_numgreen_fungi_rate = green_fungi_num / fungi_numblue_fungi_rate = blue_fungi_num / fungi_numif red_fungi_rate <= rate_limit:extension_rate_list[0] = fungi_num / red_fungi_num# total_extension_rate = sum(extension_rate_list)# extension_rate_list = extension_rate_list / total_extension_rateelif green_fungi_rate <= rate_limit:extension_rate_list[1] = fungi_num / green_fungi_num# total_extension_rate = sum(extension_rate_list)# extension_rate_list = extension_rate_list / total_extension_rateelif blue_fungi_rate <= rate_limit:extension_rate_list[2] = fungi_num / blue_fungi_num# total_extension_rate = sum(extension_rate_list)# extension_rate_list = extension_rate_list / total_extension_rateif __name__ == '__main__':patch_size = 25  # size of our patch# frame_number = 30  # number of framesone_mesh_max = 7  # the maximum number of fungi cells in a single meshinitial_fungi_num = 3  # the initial number of fungi cells as well as speciesfungi_list = []  # list of cellsmy_board = np.zeros((patch_size, patch_size))  # initialize number boardx_initial_list = []  # the x, y position of initial cellsy_initial_list = []# calculate_initial_position()  # calculate the initial position of cells for fair competitionred_fungi_num = 1  # number of each species of fungigreen_fungi_num = 1blue_fungi_num = 1rate_limit = 0.01for i_initial in range(initial_fungi_num):  # initialize fungi listx_initial = np.random.randint(patch_size)  # int(patch_size / 10) + i_initial * int(8 * patch_size / 10)  #y_initial = np.random.randint(patch_size)  # x_initial  ## fungi_list.append([x_initial_list[i_initial], y_initial_list[i_initial], 0, i_initial])  # the indexes are# x, y,fungi_list.append([x_initial, y_initial, 0, i_initial])# reproduce time and species numupdate_board(i_initial + 1)fungi_num = initial_fungi_numextension_rate_list = [0.6, 0.2, 0.2]  #extension_rate_list = np.array(extension_rate_list)# different species of fungi have different extension rateattack_parameter_list = [0.3, 0.1, 0.5]  # density parameter of different species of fungireproduce_time_to_die_list = [2, 2, 2]  # each cell reproduce two times until deathrgb_board = [[[1.0 for k in range(3)] for j in range(patch_size)] for i in range(patch_size)]  # initialize color# boardfig_1 = plt.figure(1)  # create a plotax_1 = fig_1.add_subplot(1, 1, 1)  # create an axes# ims = []  # prepared for the frames to put in# anim_1 = animation.ArtistAnimation(fig_1, ims, repeat=False)  # create an animation# save_fungi_ca_gif()  # save the gif fileanim_2 = animation.FuncAnimation(fig_1, animate_func)plt.show()# print(my_board)