Python Forum

if anyone could help me improve the efficiency of this code it would be greatly appreciated. I would like it to only loop over the newly contaminated cells at each time step instead of looping over the whole grid. Any other ideas welcome.

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import colors
import random
import os


print("=== Cellular Automata for Radioactive Waste Disposal ===")

# ======= User Inputs

L = 20.0 # km

numCells = 100 # This needs to remain the same for!

totalTime = (40 * 12) # number of time steps (in months)

print_every_n_steps = 480 # don't make this too small!
print_time_every = 200

source_x = 0.2 # x - coordinate location of source of radio active waste
source_y = 0.2 # y - coordinate location of course  of radio active waste

city_x = 12
city_y = 16
cell_city_x = int(numCells * (city_x / L))
cell_city_y = int(numCells * (city_y / L))

verbosity = 0

#Permeability Maximum Variable

# *** Model Developement 3 ***
Kmax = 3.5#Permeability Maximum Variable
#K = np.loadtxt("perm.txt", dtype='f', delimiter=' ') 
Ktmp = np.loadtxt("perm.txt")

K=Ktmp.astype(np.float)

cmap = mpl.cm.jet
norm = mpl.colors.Normalize(vmin=0, vmax = 1)
fig, ax = plt.subplots()
ax.imshow(K, cmap=cmap, norm=norm)

ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
ax.axes.get_xaxis().set_visible(False)
ax.axes.get_yaxis().set_visible(False)
results_dir = os.path.join(os.getcwd(),'solutions/')
if not os.path.isdir(results_dir):
  os.makedirs(results_dir)
    
plt.savefig(results_dir + "perm.png")

cell_source_x = int(numCells * (source_x / L))
cell_source_y = int(numCells * (source_y / L))
count = 0

for n in range(0, 20):
  C = np.zeros((numCells	,numCells ))

  C[cell_source_x,	cell_source_y] = 1; #1 is contaminated
  iteration = str(n+1)
  print("Iteration:" + iteration)
  #C[cell_city_x][cell_city_y] = 1 #will show city 

  for t in range(0,totalTime+1):
    
    if (t % print_time_every == 0) : 
      print("Time = " + str(t))

    Cnew = np.zeros((numCells	,numCells ))

    for i in range(0,numCells):

      for j in range(0,numCells):

        if (C[i][j] == 0) :#and i+1 j+1 ==0 etc

          if(verbosity > 0):
            print("Cell " + str(i) + " , " + str(j))

          # *** Model Developement 1 **
          #sij = contaminated neighbours
          
          sij = 0
          
  #using 'if' and 'elif' to assess the 9 different cases for a given cell (i,j)

          if(i!=0 and i!=numCells-1 and j!=0 and j!=numCells-1):
            sij= C[i][j+1] + C[i][j-1] + C[i-1][j] + C[i-1][j-1] + C[i-1][j+1] + C[i+1][j] + C[i+1][j-1] + C[i+1][j+1]
          
          elif(i==0 and j!=0 and j!=numCells-1):#mid left
            sij= C[i][j+1] + C[i][j-1] + C[i+1][j+1] + C[i+1][j] + C[i+1][j-1]

          elif(i==numCells-1 and j!=0 and j!=numCells-1):#mid right
            sij= C[i][j+1] + C[i][j-1] + C[i-1][j+1] + C[i-1][j] + C[numCells-1][numCells-1]
          
          elif(i!=0 and i!=numCells-1 and j==numCells-1):#mid top
            sij= C[i-1][j] + C[i-1][j-1] + C[i][j-1] + C[i+1][j-1] + C[i+1][j]
            
          elif(i!=0 and i!=numCells-1 and j==0):#mid bottom
            sij= C[i-1][j] + C[i-1][j+1] + C[i][j+1] + C[i+1][j+1] + C[i+1][j]

          elif(i == 0 and j==0) : #bottom left
            sij = C[0][1] + C[1][0] + C[1][1]
          
          elif(i==0 and j==numCells-1):#top left
            sij = C[1][numCells-1] + C[0][numCells-2] + C[1][numCells-2]
          
          elif(i==numCells-1 and j==0):#bottom right
            sij = C[numCells-1][1] + C[numCells-2][0] + C[numCells-2][1]
          
          else:#top right
            sij = C[numCells-1][numCells-2] + C[numCells-2][numCells-1] + C[numCells-2][numCells-2]
          
          
          # *** Model Developement 2 ***
          p = (sij/8)*((K[i][j])/Kmax) #gives probability
          n = random.uniform(0,1)
          if p>n:
            Cnew[i][j]=1#make this =p to give range of intensity?


        else :

          Cnew[i][j] = 0
            
    
    C = Cnew; # Set new contamination Matrix to old C.
    #this is how it goes through, adding on each time
    # Plot grid and save to file

    if (t % print_every_n_steps == 0) : 
      cmap = mpl.cm.jet
      norm = mpl.colors.Normalize(vmin=0, vmax = 1)
      fig, ax = plt.subplots()
      ax.imshow(C, cmap=cmap, norm=norm)

      ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
      ax.axes.get_xaxis().set_visible(False)
      ax.axes.get_yaxis().set_visible(False)
      results_dir = os.path.join(os.getcwd(),'solutions/')
      if not os.path.isdir(results_dir):
        os.makedirs(results_dir)
      
      plt.savefig(results_dir + "Run:" + (iteration) + "_solution_time_" + str(t) + ".png")

  city = C[cell_city_x][cell_city_y]    
  if city == 1:
    print("City Dead")
    count += 1
  else:
    print("City Alive")
  print(count)

Just keep a list of the contaminated cells and loop over that. They seem to get contaminated at only one point, so you can just append them at that point. If you only want the newly contaminated cells, reset the list each generation.

(Feb-12-2019, 03:38 PM)ichabod801 Wrote: [ -> ]Thread

I understand what you mean i'm just struggling to implement it but will continue trying. any extra help would be perfect.

When you do Cnew[i][j] = 1, just append a tuple of the coordinates to you contaminated list: contaminated.append((i, j)). Then you loop through that list with for cx, xy in contaminated:.

Now, you'll have to do the check for the next set of cells in reverse. I would suggest a dictionary of tuples, reset each generation to empty. Ideally this would be a defaultdict, but you could just use get with a 0 default. So go through the neighbors of each contaminated space, and add one to the dictionary for the neighbor's coordinates. After all the contaminated squares are checked, loop through the dictionary and use the totals to check for new contaminations. Add those to your contaminated list.

(Feb-12-2019, 08:22 PM)ichabod801 Wrote: [ -> ]When you do Cnew[i][j] = 1, just append a tuple of the coordinates to you contaminated list: contaminated.append((i, j)). Then you loop through that list with for cx, xy in contaminated:. Now, you'll have to do the check for the next set of cells in reverse. I would suggest a dictionary of tuples, reset each generation to empty. Ideally this would be a defaultdict, but you could just use get with a 0 default. So go through the neighbors of each contaminated space, and add one to the dictionary for the neighbor's coordinates. After all the contaminated squares are checked, loop through the dictionary and use the totals to check for new contaminations. Add those to your contaminated list.

you make it sound so easy!! i have tried to implement this but i have had no luck.

Does this work? I don't have the files necessary to test it.

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import colors
import random
import os
import collections


print("=== Cellular Automata for Radioactive Waste Disposal ===")

# ======= User Inputs

L = 20.0 # km

numCells = 100 # This needs to remain the same for!

totalTime = (40 * 12) # number of time steps (in months)

print_every_n_steps = 480 # don't make this too small!
print_time_every = 200

source_x = 0.2 # x - coordinate location of source of radio active waste
source_y = 0.2 # y - coordinate location of course  of radio active waste

city_x = 12
city_y = 16
cell_city_x = int(numCells * (city_x / L))
cell_city_y = int(numCells * (city_y / L))

verbosity = 0

#Permeability Maximum Variable

# *** Model Developement 3 ***
Kmax = 3.5#Permeability Maximum Variable
#K = np.loadtxt("perm.txt", dtype='f', delimiter=' ')
Ktmp = np.loadtxt("perm.txt")

K=Ktmp.astype(np.float)

cmap = mpl.cm.jet
norm = mpl.colors.Normalize(vmin=0, vmax = 1)
fig, ax = plt.subplots()
ax.imshow(K, cmap=cmap, norm=norm)

ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
ax.axes.get_xaxis().set_visible(False)
ax.axes.get_yaxis().set_visible(False)
results_dir = os.path.join(os.getcwd(),'solutions/')
if not os.path.isdir(results_dir):
  os.makedirs(results_dir)

plt.savefig(results_dir + "perm.png")

cell_source_x = int(numCells * (source_x / L))
cell_source_y = int(numCells * (source_y / L))
count = 0

for n in range(0, 20):
  C = np.zeros((numCells    ,numCells ))

  C[cell_source_x,  cell_source_y] = 1; #1 is contaminated
  contaminated = set()
  contaminated.add((cell_source_x, cell_source_y))
  iteration = str(n+1)
  print("Iteration:" + iteration)
  #C[cell_city_x][cell_city_y] = 1 #will show city

  for t in range(0,totalTime+1):

    if (t % print_time_every == 0) :
      print("Time = " + str(t))

    Cnew = np.zeros((numCells   ,numCells ))

    neighbors = collections.defaultdict(int)
    for source_x, source_y in contaminated:
      for off_x, off_y in [(0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1), (-1, 0), (-1, 1)]:
        neighbor_x, neighbor_y = source_x + off_x, source_y + off_y
        if 0 <= neighbor_x < numCells and 0 <= neighbor_y <= numCells:
          neighbors[(neighbor_x, neighbor_y)] += 1

    for coordinates, count in neighbors.items():
      possible_x, possible_y = coordinates
      # *** Model Developement 2 ***
      p = (count/8)*((K[possible_x][possible_y])/Kmax) #gives probability
      n = random.uniform(0,1)
      if p>n:
        Cnew[possible_x][possible_y]=1#make this =p to give range of intensity?
        contaminated.add(coordinates)



    C = Cnew; # Set new contamination Matrix to old C.
    #this is how it goes through, adding on each time
    # Plot grid and save to file

    if (t % print_every_n_steps == 0) :
      cmap = mpl.cm.jet
      norm = mpl.colors.Normalize(vmin=0, vmax = 1)
      fig, ax = plt.subplots()
      ax.imshow(C, cmap=cmap, norm=norm)

      ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
      ax.axes.get_xaxis().set_visible(False)
      ax.axes.get_yaxis().set_visible(False)
      results_dir = os.path.join(os.getcwd(),'solutions/')
      if not os.path.isdir(results_dir):
        os.makedirs(results_dir)

      plt.savefig(results_dir + "Run:" + (iteration) + "_solution_time_" + str(t) + ".png")

  city = C[cell_city_x][cell_city_y]
  if city == 1:
    print("City Dead")
    count += 1
  else:
    print("City Alive")
  print(count)

(Feb-19-2019, 08:04 PM)ichabod801 Wrote: [ -> ]import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import colors
import random
import os
import collections
 
 
print("=== Cellular Automata for Radioactive Waste Disposal ===")
 
# ======= User Inputs
 
L = 20.0 # km
 
numCells = 100 # This needs to remain the same for!
 
totalTime = (40 * 12) # number of time steps (in months)
 
print_every_n_steps = 480 # don't make this too small!
print_time_every = 200
 
source_x = 0.2 # x - coordinate location of source of radio active waste
source_y = 0.2 # y - coordinate location of course  of radio active waste
 
city_x = 12
city_y = 16
cell_city_x = int(numCells * (city_x / L))
cell_city_y = int(numCells * (city_y / L))
 
verbosity = 0
 
#Permeability Maximum Variable
 
# *** Model Developement 3 ***
Kmax = 3.5#Permeability Maximum Variable
#K = np.loadtxt("perm.txt", dtype='f', delimiter=' ')
Ktmp = np.loadtxt("perm.txt")
 
K=Ktmp.astype(np.float)
 
cmap = mpl.cm.jet
norm = mpl.colors.Normalize(vmin=0, vmax = 1)
fig, ax = plt.subplots()
ax.imshow(K, cmap=cmap, norm=norm)
 
ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
ax.axes.get_xaxis().set_visible(False)
ax.axes.get_yaxis().set_visible(False)
results_dir = os.path.join(os.getcwd(),'solutions/')
if not os.path.isdir(results_dir):
  os.makedirs(results_dir)
 
plt.savefig(results_dir + "perm.png")
 
cell_source_x = int(numCells * (source_x / L))
cell_source_y = int(numCells * (source_y / L))
count = 0
 
for n in range(0, 20):
  C = np.zeros((numCells    ,numCells ))
 
  C[cell_source_x,  cell_source_y] = 1; #1 is contaminated
  contaminated = set()
  contaminated.add((cell_source_x, cell_source_y))
  iteration = str(n+1)
  print("Iteration:" + iteration)
  #C[cell_city_x][cell_city_y] = 1 #will show city
 
  for t in range(0,totalTime+1):
 
    if (t % print_time_every == 0) :
      print("Time = " + str(t))
 
    Cnew = np.zeros((numCells   ,numCells ))
 
    neighbors = collections.defaultdict(int)
    for source_x, source_y in contaminated:
      for off_x, off_y in [(0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1), (-1, 0), (-1, 1)]:
        neighbor_x, neighbor_y = source_x + off_x, source_y + off_y
        if 0 <= neighbor_x < numCells and 0 <= neighbor_y <= numCells:
          neighbors[(neighbor_x, neighbor_y)] += 1
 
    for coordinates, count in neighbors.items():
      possible_x, possible_y = coordinates
      # *** Model Developement 2 ***
      p = (count/8)*((K[possible_x][possible_y])/Kmax) #gives probability
      n = random.uniform(0,1)
      if p>n:
        Cnew[possible_x][possible_y]=1#make this =p to give range of intensity?
        contaminated.add(coordinates)
 
 
 
    C = Cnew; # Set new contamination Matrix to old C.
    #this is how it goes through, adding on each time
    # Plot grid and save to file
 
    if (t % print_every_n_steps == 0) :
      cmap = mpl.cm.jet
      norm = mpl.colors.Normalize(vmin=0, vmax = 1)
      fig, ax = plt.subplots()
      ax.imshow(C, cmap=cmap, norm=norm)
 
      ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
      ax.axes.get_xaxis().set_visible(False)
      ax.axes.get_yaxis().set_visible(False)
      results_dir = os.path.join(os.getcwd(),'solutions/')
      if not os.path.isdir(results_dir):
        os.makedirs(results_dir)
 
      plt.savefig(results_dir + "Run:" + (iteration) + "_solution_time_" + str(t) + ".png")
 
  city = C[cell_city_x][cell_city_y]
  if city == 1:
    print("City Dead")
    count += 1
  else:
    print("City Alive")
  print(count)

that looks to be loads faster thank you so much for your help
i get this error come up after running it a few times

Traceback (most recent call last):
File "main.py", line 90, in <module>
Cnew[possible_x][possible_y]=1#make this =p to give range of intensity?
IndexError: index 100 is out of bounds for axis 0 with size 100

Whoops, change line 81 to this:

        if 0 <= neighbor_x < numCells and 0 <= neighbor_y < numCells:

(Feb-19-2019, 08:04 PM)ichabod801 Wrote: [ -> ]if 0 <= neighbor_x < numCells and 0 <= neighbor_y < numCells:

(Feb-19-2019, 10:08 PM)ichabod801 Wrote: [ -> ]Whoops, change line 81 to this:

 if 0 <= neighbor_x < numCells and 0 <= neighbor_y < numCells: 

It seems to be running well now thank you the solutions we get from it seems to be different to the previous code the spread of waste seems to be different. almost like it skips loads.

I've gone over the code, I don't see anything that should be making a functional difference. The only thing I explicitly left out was setting Cnew cells to 0, I figured that was taken care of by the initialization using np.zeros. I'm not sure what you mean by "skipping load" though.