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

import matplotlib.pyplot as plt
import random
import numpy as np
import matplotlib.animation as animation


def save_fungi_ca_gif():  # save the gif file to the path
    target_gif_path = "E:/engineering space/figure/gif format/"
    target_gif_name = "aggressive_ca.gif"
    target_gif_full = target_gif_path + target_gif_name
    anim_1.save(target_gif_full, writer='pillow')  # save the animation
    print('Saved')


def calculate_random_neighbour(stay, left, up, right, down, no_reproduce):  # get the relative location of the newly
    total_rate = stay + left + up + right + down + no_reproduce  # reproduced fungi cell
    if total_rate == 0.0:
        total_rate = 1
    cap1 = stay / total_rate  # cap 1 to 5 is the boundary of probability area
    cap2 = cap1 + left / total_rate
    cap3 = cap2 + up / total_rate
    cap4 = cap3 + right / total_rate
    cap5 = cap4 + down / total_rate
    rnd = random.random()
    if 0 <= rnd < cap1:
        return 0
    elif cap1 <= rnd < cap2:  # divide the range into 6 parts
        return 1
    elif cap2 <= rnd < cap3:  # randomly select one part
        return 2
    elif cap3 <= rnd < cap4:
        return 3
    elif cap4 <= rnd < cap5:
        return 4
    else:
        return 5


def get_extension_parameters(i_get_extension_parameters):  # get the extension parameters
    extension_parameters = extension_rate_list[fungi_list[i_get_extension_parameters][3]]
    return extension_parameters


def get_random_neighbour(i_get_random_neighbour):  # calculate the density of neighbour mesh as well as self mesh
    extension_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_max
    stay_rate = extension_rate_get_random_neighbour - stay_density
    if 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_max
        left_rate = extension_rate_get_random_neighbour * (
                    1 - left_density + 0.1)  # further calculate the transferring rate of
        # the new fungi
    else:
        left_density = 0.5
        left_rate = 0

    if 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_max
        up_rate = extension_rate_get_random_neighbour * (1 - up_density + 0.1)
    else:
        up_density = 0.5
        up_rate = 0

    if 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_max
        right_rate = extension_rate_get_random_neighbour * (1 - right_density + 0.1)
    else:
        right_density = 0.5
        right_rate = 0

    if 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_max
        down_rate = extension_rate_get_random_neighbour * (1 - down_density + 0.1)
    else:
        down_density = 0.5
        down_rate = 0

    total_density = stay_density + left_density + up_density + right_density + down_density
    no_reproduce_rate = total_density / 8
    neighbour_location = calculate_random_neighbour(stay_rate, left_rate, up_rate, right_rate, down_rate,
                                                    no_reproduce_rate)
    return neighbour_location


def reproduce_one_cell(i_reproduce_one_cell):  # reproduce a single fungi cell
    iterations = 0
    while True:  # loop until the random neighbour meets the boundary conditions
        iterations += 1
        # print('iteration:', iterations)
        reproduce_location = get_random_neighbour(i_reproduce_one_cell)
        # the followings are boundary conditions
        if iterations < 3:
            if reproduce_location == 0:  # stay in one mesh
                current_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 mesh
                    break  # this is the location wanted
                else:
                    # print(iterations)  # how many loops have been taken
                    continue  # random again

            if reproduce_location == 1:  # go to left mesh
                if fungi_list[i_reproduce_one_cell][1] == 0:  # reaches the left boundary
                    # print(iterations)
                    continue  # random again
                elif 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 wanted
                else:
                    continue

            if reproduce_location == 2:  # go to up mesh
                if fungi_list[i_reproduce_one_cell][0] == 0:  # reaches the up boundary
                    # print(iterations)
                    continue  # random again
                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 wanted
                else:
                    continue

            if reproduce_location == 3:  # go to right mesh
                if fungi_list[i_reproduce_one_cell][1] == patch_size - 1:  # reaches the right boundary
                    # print(iterations)
                    continue  # random again
                elif 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 wanted
                else:
                    continue

            if reproduce_location == 4:  # go to down mesh
                if fungi_list[i_reproduce_one_cell][0] == patch_size - 1:  # reaches the down boundary
                    continue
                    # 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 wanted
                else:
                    continue
        else:
            reproduce_location = 5
            break

    global fungi_num
    # take reproducing activity
    if 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 fungi
        if 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 die
            my_board[fungi_list[i_reproduce_one_cell][0]][fungi_list[i_reproduce_one_cell][1]] -= 1  # remove fungi
            # in number board
            del fungi_list[i_reproduce_one_cell]  # the death of a fungi
        fungi_num = len(fungi_list)
        update_board(fungi_num)  # get a board according to the new fungi reproduced
    else:
        fungi_list[i_reproduce_one_cell][2] += 1  # the life span of a fungi


def reproduce_round():  # reproduce all the current fungi cells
    for i_reproduce_round in range(fungi_num):
        reproduce_one_cell(i_reproduce_round)  # reproduce a single fungi cell


def update_board(fungi_num_update_board):  # get a number board according to the new fungi
    my_board[fungi_list[fungi_num_update_board - 1][0]][fungi_list[fungi_num_update_board - 1][1]] += 1


def 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)]  # initialize
    global red_fungi_num
    global green_fungi_num
    global blue_fungi_num
    red_fungi_num = 0
    green_fungi_num = 0
    blue_fungi_num = 0
    # color board
    for i_paint_rgb_board in range(fungi_num):
        if fungi_list[i_paint_rgb_board][3] == 0:
            red_fungi_num += 1
            paint_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 += 1
            paint_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 += 1
            paint_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_board


def attack_round():
    global fungi_num
    for i_attack_round in range(fungi_num):
        if i_attack_round == fungi_num - 1:
            break
        if fungi_list[i_attack_round][2] >= reproduce_time_to_die_list[fungi_list[i_attack_round][3]]:
            # when reproduced an exact times to die
            my_board[fungi_list[i_attack_round][0]][fungi_list[i_attack_round][1]] -= 1  # remove fungi in number board
            del fungi_list[i_attack_round]  # the death of a fungi
            fungi_num = len(fungi_list)
            # update_board(fungi_num)  # get a board according to the new fungi reproduced
        else:
            continue


def update():  # update the board information and fungi list in a round
    reproduce_round()  # reproduce all the current fungi cells
    rgb_board_update = get_rgb_board()  # get a color board after a round of reproduction
    attack_round()
    return rgb_board_update


def calculate_initial_position():  # calculate the initial position of cells for fair competition
    radius_calculate_initial_position = patch_size / 6  # distance between initial cells
    middle_point_initial_position = patch_size / 2
    for i_calculate_initial_position in range(initial_fungi_num):
        rotation_degree = i_calculate_initial_position / initial_fungi_num * 2 * np.pi
        x_initial_position = middle_point_initial_position + np.cos(rotation_degree) * radius_calculate_initial_position
        y_initial_position = middle_point_initial_position + np.sin(rotation_degree) * radius_calculate_initial_position
        x_initial_list.append(int(round(x_initial_position)))
        y_initial_list.append(int(round(y_initial_position)))


def animate_func(frame_sequence):  # depict the animation
    ax_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_list
    red_fungi_rate = red_fungi_num / fungi_num
    green_fungi_rate = green_fungi_num / fungi_num
    blue_fungi_rate = blue_fungi_num / fungi_num
    if 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_rate
    elif 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_rate
    elif 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_rate


if __name__ == '__main__':
    patch_size = 25  # size of our patch
    # frame_number = 30  # number of frames
    one_mesh_max = 7  # the maximum number of fungi cells in a single mesh
    initial_fungi_num = 3  # the initial number of fungi cells as well as species
    fungi_list = []  # list of cells
    my_board = np.zeros((patch_size, patch_size))  # initialize number board
    x_initial_list = []  # the x, y position of initial cells
    y_initial_list = []
    # calculate_initial_position()  # calculate the initial position of cells for fair competition
    red_fungi_num = 1  # number of each species of fungi
    green_fungi_num = 1
    blue_fungi_num = 1
    rate_limit = 0.01
    for i_initial in range(initial_fungi_num):  # initialize fungi list
        x_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 num
        update_board(i_initial + 1)
    fungi_num = initial_fungi_num
    extension_rate_list = [0.6, 0.2, 0.2]  #
    extension_rate_list = np.array(extension_rate_list)
    # different species of fungi have different extension rate
    attack_parameter_list = [0.3, 0.1, 0.5]  # density parameter of different species of fungi
    reproduce_time_to_die_list = [2, 2, 2]  # each cell reproduce two times until death
    rgb_board = [[[1.0 for k in range(3)] for j in range(patch_size)] for i in range(patch_size)]  # initialize color
    # board
    fig_1 = plt.figure(1)  # create a plot
    ax_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 file
    anim_2 = animation.FuncAnimation(fig_1, animate_func)
    plt.show()
    # print(my_board)