How to write to a matrix within a while loop in a monte carlo simulation in python 3

Question

Okay so I am trying to complete the infection simulation code for a monte carlo simulation in python. We were given the thread shell and just need to complete this. I will also have to add a vaccine function but it should be pretty similar to the infect function so Im trying to make sure that it works first. I know that I need to read from my original matrix but write to a new matrix for each day/iteration, but I can't figure out how to write to the new matrix. I tried using np.append but it forced me to redefine my original matrix A, or it said that list indices must be integers, not floats. I've tried working through examples on here and other places, but they all seem to be using lists or they don't involve calling a function within the while loop. Any help would be greatly appreciated.

import random
import math
import numpy as np



def infect(Pop,i,j,n,m,tau):
    t = 0
    if (i > 1) and (i < n) and (j > 1) and (j < m):
        if (Pop[i-1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j+1]>0):
            t = (np.random.rand() < tau)
        if (Pop[i+1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j-1]>0):
            t = (np.random.rand() < tau)
    if (i == 1) and (j == 1):
        if (Pop[i,j+1]>0):
            t = (np.random.rand() < tau)
        if (Pop[i+1,j]):
            t = (np.random.rand() < tau)
    if (i == 1) and (j != m) and (j > 1):
        if (Pop[i,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i+1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j-1]>0):
            t = (np.random.rand() < tau)
    if (i == 1) and (j == m):
        if (Pop[i+1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j-1]>0):
            t = (np.random.rand() < tau)
    if (i == n) and (j == 1):
        if (Pop[i-1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j+1]>0):
            t = (np.random.rand() < tau)
    if (i < n) and (i > 1) and (j == 1):
        if (Pop[i-1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j+1]>0):
            t = (np.random.rand() < tau)
        if (Pop[i+1,j]>0):
            t = (np.random.rand() < tau)
    if (i < n) and (i > 1) and (j == m):
        if (Pop[i-1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i+1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j-1]>0):
            t = (np.random.rand() < tau)
    if (i == n) and (j > 1) and (j < m):
        if (Pop[i-1,j]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j+1]>0):
            t = (np.random.rand() < tau)
        if (Pop[i,j-1]>0):
            t = (np.random.rand() < tau)
    if (i == n) and (j == m):
        if (Pop[i,j-1]>0):
            t = (np.random.rand() < tau)
        if (Pop[i-1,j]>0):
            t = (np.random.rand() < tau)       

    p = 0
    if (t==True):
        p = 1
    return p


i = 1
j = 1
n = 10
m = 10
k = int(input("Number of Days to Recover from Illness?"))
d = 0.0
tau = 0.5
mu = 0.2
A = np.zeros((n,m))
if d == 0:
    n1 = random.sample(range(n),1)
    m1 = random.sample(range(m),1)
    A[n1,m1] = 1
    print(A)

while d < 100:
    while True:

        if (A[i,j]==0):
            x = infect(A,i,j,n,m,tau)
        print(x)
                #A_new.append(x)

Answer 1

You create an extra matrix and every time through the loop you swap the two references. For example,

A1 = np.zeros((m, n))
A2 = np.zeros((m, n))
Anow = A1 # reference, not copy
Aafter = A2
while d < 100:
    x = infectAll(Anow, Aafter, n,m,tau)
    Anow, Aafter = Aafter, Anow

and infectAll() sweeps over the whole matrix so it would be something like

def infectAll(Ain, Aout, n, m, tau):
    for i in range(m):
        for j in range(n):
             if Anow[i,j] == 0:
                  Aafter[i,j] = infect(Anow, i, j, n, m, tau)

I like Andrei's more compact code but there is no need to create a new copy of A every time, so the best would be to combine the above technique with Andrei's approach.

Answer 2

As I understand you need dynamic of your infection matrix and tau is probability of neighbor infection. You can use 3-dimensional array for this and optimize your code like this:

from copy import copy
import numpy as np
def infect(A, tau):
    B = copy(A)
    for i in range(m):
        for j in range(n):
            is_infected = False
            for neighbor in [A[i-1,j], A[i+1,j], A[i,j-1], A[i,j+1]]:
                if neighbor: 
                    B[i,j] = int(A[i,j] or (np.random.rand() < tau))
    return B

D = np.zeros((T + 1, m, n))
A = np.zeros((m, n))
A[i,j] = 1
for t in range(T):
    D[t,:,:] = A
    A = infect(A, tau)
D[T,:,:] = A