Posts: 5
Threads: 1
Joined: Feb 2019
Feb-12-2019, 02:06 PM
(This post was last modified: Feb-12-2019, 02:08 PM by buran.)
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 |
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 ===")
L = 20.0
numCells = 100
totalTime = (40 * 12)
print_every_n_steps = 480
print_time_every = 200
source_x = 0.2
source_y = 0.2
city_x = 12
city_y = 16
cell_city_x = int(numCells * (city_x / L))
cell_city_y = int(numCells * (city_y / L))
verbosity = 0
Kmax = 3.5
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;
iteration = str(n+1)
print("Iteration:" + iteration)
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) :
if(verbosity > 0):
print("Cell " + str(i) + " , " + str(j))
sij = 0
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):
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):
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):
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):
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) :
sij = C[0][1] + C[1][0] + C[1][1]
elif(i==0 and j==numCells-1):
sij = C[1][numCells-1] + C[0][numCells-2] + C[1][numCells-2]
elif(i==numCells-1 and j==0):
sij = C[numCells-1][1] + C[numCells-2][0] + C[numCells-2][1]
else:
sij = C[numCells-1][numCells-2] + C[numCells-2][numCells-1] + C[numCells-2][numCells-2]
p = (sij/8)*((K[i][j])/Kmax)
n = random.uniform(0,1)
if p>n:
Cnew[i][j]=1
else :
Cnew[i][j] = 0
C = Cnew;
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)
|
Posts: 4,220
Threads: 97
Joined: Sep 2016
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.
Posts: 5
Threads: 1
Joined: Feb 2019
(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.
Posts: 4,220
Threads: 97
Joined: Sep 2016
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.
Posts: 5
Threads: 1
Joined: Feb 2019
(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.
Posts: 4,220
Threads: 97
Joined: Sep 2016
Does this work? I don't have the files necessary to test it.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 |
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 ===")
L = 20.0
numCells = 100
totalTime = (40 * 12)
print_every_n_steps = 480
print_time_every = 200
source_x = 0.2
source_y = 0.2
city_x = 12
city_y = 16
cell_city_x = int(numCells * (city_x / L))
cell_city_y = int(numCells * (city_y / L))
verbosity = 0
Kmax = 3.5
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;
contaminated = set()
contaminated.add((cell_source_x, cell_source_y))
iteration = str(n+1)
print("Iteration:" + iteration)
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
p = (count/8)*((K[possible_x][possible_y])/Kmax)
n = random.uniform(0,1)
if p>n:
Cnew[possible_x][possible_y]=1
contaminated.add(coordinates)
C = Cnew;
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)
|
Posts: 5
Threads: 1
Joined: Feb 2019
Feb-19-2019, 09:25 PM
(This post was last modified: Feb-19-2019, 09:25 PM by benbrown03.)
(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
Posts: 4,220
Threads: 97
Joined: Sep 2016
Whoops, change line 81 to this:
1 |
if 0 <= neighbor_x < numCells and 0 <= neighbor_y < numCells:
|
Posts: 5
Threads: 1
Joined: Feb 2019
(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:
1 |
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.
Posts: 4,220
Threads: 97
Joined: Sep 2016
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.
|
