compnew.py

import numpy as np
import matplotlib.pyplot as plt
import time

# Grid Structure = x * y * 2(Type, Number)
# Producer - 1
# Resistant - 2
# Sensitive - 3
# Killed - 4

def seed_random_grid(grid, num):
    c = np.random.randint(grid.shape[0], size=(num,2))
    c = tuple(map(tuple, c))
    for i in c:
        grid[i] = 1
    return grid

Round = 0
Final_Round = 30

# gridsize = 10**5

# inhibitionradius = gridsize/100
# grazingradius = gridsize/100

# producers_start = 10
# resistant_start = 10**6
# sensitive_start = 10**8

gridsize = 100

inhibitionradius = gridsize/10
grazingradius = gridsize/10

producers_start = 10
resistant_start = 100
sensitive_start = 10000


producers = producers_start
resistant = resistant_start
sensitive = sensitive_start

# Initialize Grids
producers_grid = np.zeros((gridsize,gridsize))
resistant_grid = np.zeros((gridsize,gridsize))
sensitive_grid = np.zeros((gridsize,gridsize))
resources_grid = np.zeros((gridsize,gridsize))

# Seed grid
# 1 is present. 0 is absent. -1 is killed
print("Seeding producers grid")
producers_grid = seed_random_grid(producers_grid, producers_start)
resistant_grid = seed_random_grid(resistant_grid, resistant_start)
sensitive_grid = seed_random_grid(sensitive_grid, sensitive_start)
resources = [0, 0, 0]              # Producers, Sensitive, Ressitant
grids = [producers_grid, sensitive_grid, resistant_grid]

while Round < Final_Round and len(sensitive_grid[sensitive_grid == 1]) > 0:
    s = time.time()
    # Get all grid positions of non empty cells in producers_grid
    c = np.dstack(np.where(producers_grid == 1))[0]
    for cx,cy in c:
        y,x = np.ogrid[-cx:gridsize-cx, -cy:gridsize-cy]
        mask = (x**2 + y**2 <= inhibitionradius ** 2) & (sensitive_grid[y+cx,x+cy] == 1)
        sensitive_grid[mask] = -1
    print("Killed sensitive cells: "+str(time.time() - s))
    # Distribution of resources
    idx = list(np.ndindex(producers_grid.shape))
    for i in idx:
        _s = time.time()
        y,x = np.ogrid[-i[0]:gridsize-i[0], -i[1]:gridsize-i[1]]
        cx = i[0]
        cy = i[1]
        mask = (x**2 + y**2 <= grazingradius ** 2) & ((producers_grid[y+cx,x+cy] == 1) | (sensitive_grid[y+cx,x+cy] == 1) | (resistant_grid[y+cx,x+cy] == 1))
        total_resource = 1
        resources_grid[mask] += total_resource/np.sum(mask)
        for j in range(0, len(grids)):
            _ = (x**2 + y**2 <= grazingradius ** 2) & (grids[j][y+cx,x+cy] == 1)
            resources[j] += (np.sum(_) * total_resource/np.sum(mask))
    e = time.time()
    Round += 1
    print(Round)
    print("Elapsed: "+str(e - s))