# 2-D HOT Forest Model¶

In [519]:
import numpy as np
import random
import matplotlib.pyplot as plt
from scipy import stats
#from scipy.stats import expon
#from scipy.stats import boltzmann
from matplotlib import rc
rc('text', usetex=True)
from scipy.ndimage import measurements
import math
from itertools import groupby
from operator import itemgetter
import matplotlib.animation as animation


## Defining our Spark Probabilities¶

In [527]:
L = 32
l = L/10

probs = np.zeros([L,L])
def p_row(p,L):
return math.exp(-p/l)

row = np.array([p_row(i,L) for i in range(L)])
prow = row/sum(row)
probs = np.array([prow[i]*prow[j] for i in range(L) for j in range(L)])
probs = probs.reshape(L,L)


#### Visualizing the Lightening Strike Probabilities¶

In [528]:
plt.imshow(probs, cmap='viridis',interpolation='nearest')
plt.show()


## Lightning Strike Forest Madness¶

In [529]:
D = [L**2]

In [530]:
for d in D:
land = np.zeros([L,L]) #makes empty land for new D trial
EMPTY = [[i,j]for i in range(L) for j in range(L)] #fills EMPTY array holding tree-less spaces on land
#FULL = []
YIELD = [] #starts array to record yield for plot
k=0
xs = []
ys = []
c = 0
while len(EMPTY) != 0:
COST = {} #dict to hold the location and yield for each spot in D to compare
for spot in range(d):
em = np.random.randint(0,len(EMPTY)) # returns random int in EMPTY
x = EMPTY[em][0]
y = EMPTY[em][1]
land[x][y] =1  #puts a tree in empty spot
forests, num = measurements.label(land)
area = measurements.sum(land, forests, index=np.arange(forests.max() + 1))
LAND = np.array([land[i][j]*area[forests[i][j]] for i in range(L) for j in range(L)])
LAND = LAND.reshape(L,L)
fire = [LAND*probs]
COST[np.sum(fire)] = [x,y,sum(sum(land))-np.sum(fire)] #records/calc.s the loc and Yield for spot
land[x][y] = 0  #takes tree away
xs.append(COST[min(COST)][0]) #records x coord. spot of minimum avg damage
ys.append(COST[min(COST)][1]) #records y coord. spot of minimum avg damage
YIELD.append([k+1,COST[min(COST)][2]]) #records tree placement number and yield
land[xs[k]][ys[k]] = 1  #puts a tree in place
#FULL.append([xs[k],ys[k]])
full = [xs[k],ys[k]]
EMPTY = [o for o in EMPTY if o != full] #takes tree location out of the EMPTY array
if YIELD[k][1] < YIELD[k-1][1]-25 and c == 0: #change diff as L increases 25, 65, 315
c = YIELD[k-1][0]
cyield = YIELD[k-1][1]
plt.figure(figsize = (10,10))
plt.axes(xlim=(-1, L), ylim=(-1, L))
plt.scatter(xs,ys,color='g',s=150) #s=5:L=128, s=22:L=64, s=150:L=32
plt.title('Tree Layout at Peak Yield '+str(cyield)+', for L = '+str(L)+' and for D = '+str(d),Fontsize=18)
plt.savefig('/Users/sarahhowerter/Google Drive/UVM/PoCS/Assignments/Assignment #7/PeakPlot-D'+str(d)+'L'+str(L)+'.png',bbox_inches='tight')
AREA = np.sort(area)
forestsize, numforests = zip(*Counter(AREA).items())
forests = [str(int(forestsize[i])) for i in range(len(forestsize))]
size = [int(numforests[i]) for i in range(len(numforests))]
plt.figure(figsize = (10,3))
plt.xticks(np.arange(len(forests)),forests)
plt.bar(np.arange(len(forests)),size,align='center',color='olive',edgecolor='black',linewidth=0.5)
plt.title('Forest Size Distribution at Peak Yield '+str(cyield)+', for L = '+str(L)+' and for D = '+str(d),Fontsize=14)
plt.savefig('/Users/sarahhowerter/Google Drive/UVM/PoCS/Assignments/Assignment #7/ForDist-D'+str(d)+'L'+str(L)+'.png',bbox_inches='tight')
plt.show()
plt.figure(figsize = (10,10))
plt.loglog(forestsize,numforests,'.',color='r',ms=15)
plt.title('Forest Size Distribution on a LogLog Plot at Peak Yield '+str(cyield)+', for L = '+str(L)+' and for D = '+str(d),Fontsize=14)
plt.savefig('/Users/sarahhowerter/Google Drive/UVM/PoCS/Assignments/Assignment #7/ForDistPlot-D'+str(d)+'L'+str(L)+'.png',bbox_inches='tight')
print('The peak Yield, '+str(cyield)+' occurs at the '+str(c)+'th tree placement' )
if k % 100 == 0:
print('The '+str(k+1)+'th tree has been planted.')
k += 1
trees = [x[0]/(L**2) for x in YIELD]
Yield = [y[1]/(L**2) for y in YIELD]
plt.figure(figsize = (10,10))
plt.axes(xlim=(0,1), ylim=(0,1))
plt.scatter(trees, Yield, color ='g',s=10)
plt.axvline((c/L**2),color = 'y')
plt.title('Yield Curve for L = '+str(L)+'and for D = '+str(d),Fontsize=18)
plt.xlabel('Density',Fontsize=14)
plt.ylabel('Yield',Fontsize=14)
plt.text(0.1,0.9,'Density at Peak Yield = '+str(c/(L**2)),Fontsize=12,bbox=dict(facecolor='none', edgecolor='olive', boxstyle='round,pad=1'))
plt.savefig('/Users/sarahhowerter/Google Drive/UVM/PoCS/Assignments/Assignment #7/YieldPlot-D'+str(d)+'L'+str(L)+'.png',bbox_inches='tight')
plt.show()


The 1th tree has been planted.
The 101th tree has been planted.
The 201th tree has been planted.
The 301th tree has been planted.
The 401th tree has been planted.
The 501th tree has been planted.
The 601th tree has been planted.
The 701th tree has been planted.
The 801th tree has been planted.
The 901th tree has been planted.

/Applications/anaconda/lib/python3.6/site-packages/matplotlib/axes/_base.py:3193: UserWarning: Attempting to set identical bottom==top results
in singular transformations; automatically expanding.
bottom=1.0, top=1.0
'bottom=%s, top=%s') % (bottom, top))

The peak Yield, 921.673632302 occurs at the 990th tree placement
The 1001th tree has been planted.


# Animating the Madness¶

In [335]:
L = 64
D = L
l = L/10

probs = np.zeros([L,L])
def p_row(p,L):
return math.exp(-p/l)

row = np.array([p_row(i,L) for i in range(L)])
prow = row/sum(row)
probs = np.array([prow[i]*prow[j] for i in range(L) for j in range(L)])
probs = probs.reshape(L,L)

In [336]:
land = np.zeros([L,L]) #makes empty land for new D trial
EMPTY = [[i,j]for i in range(L) for j in range(L)] #fills EMPTY array holding tree-less spaces on land
FULL = []
YIELD = [] #starts array to record yield for plot
k = 0
xs = []
ys = []
c = 0
ims = []
fig = plt.figure(figsize = (6,6))

In [337]:
while len(EMPTY) != 0:
COST = {} #dict to hold the location and yield for each spot in D to compare
for spot in range(D):
em = np.random.randint(0,len(EMPTY)) # returns random int in EMPTY
x = EMPTY[em][0]
y = EMPTY[em][1]
land[x][y] =1  #puts a tree in empty spot
forests, num = measurements.label(land)
area = measurements.sum(land, forests, index=np.arange(forests.max() + 1))
LAND = np.array([land[i][j]*area[forests[i][j]] for i in range(L) for j in range(L)])
LAND = LAND.reshape(L,L)
fire = [LAND*probs]
COST[np.sum(fire)] = [x,y,sum(sum(land))-np.sum(fire)] #records/calc.s the loc and Yield for spot
land[x][y] = 0  #takes tree away
xs.append(COST[min(COST)][0]) #records x coord. spot of minimum avg damage
ys.append(COST[min(COST)][1]) #records y coord. spot of minimum avg damage
YIELD.append([k+1,COST[min(COST)][2]])
land[xs[k]][ys[k]] = 1  #puts a tree in place
#FULL.append([xs[k],ys[k]])
full = [xs[k],ys[k]]
EMPTY = [o for o in EMPTY if o != full] #takes tree location out of the EMPTY array
if (k + 1) % int((L**2)/35) == 0:
plt.figure(figsize = (6,6))
plt.axes(xlim=(-1, L), ylim=(-1, L))
plt.scatter(xs,ys,color='g',s=14)
plt.title('Forest Layout Animation for L = '+str(L)+' & for D = '+str(D),Fontsize=14)
plt.savefig('/Users/sarahhowerter/Google Drive/UVM/PoCS/Assignments/Assignment #7/Forest-D'+str(D)+'L'+str(L)+'/ForD'+str(D)+'L'+str(L)+'-'+str(k+1)+'.png')
#im = plt.scatter(xs,ys,color='g',s=70)
#ims.append(im)
plt.show()
if YIELD[k][1] < YIELD[k-1][1]-25 and c == 0:
c = YIELD[k-1][0]
plt.figure(figsize = (6,6))
plt.axes(xlim=(-1, L), ylim=(-1, L))
plt.scatter(xs,ys,color='olive',s=14)
plt.title('Forest Layout Animation for L = '+str(L)+' and for D = '+str(D), Fontsize=14)
plt.savefig('/Users/sarahhowerter/Google Drive/UVM/PoCS/Assignments/Assignment #7/Forest-D'+str(D)+'L'+str(L)+'/ForD'+str(D)+'L'+str(L)+'-'+str(k+1)+'.png')
#im = plt.scatter(xs,ys,color='olive',s=70)
#ims.append(im)
plt.show()
k += 1

#ani = animation.ArtistAnimation(fig, ims, interval=50, blit=True, repeat_delay=1000)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-337-e135897d3660> in <module>()
8         forests, num = measurements.label(land)
9         area = measurements.sum(land, forests, index=np.arange(forests.max() + 1))
---> 10         LAND = np.array([land[i][j]*area[forests[i][j]] for i in range(L) for j in range(L)])
11         LAND = LAND.reshape(L,L)
12         fire = [LAND*probs]

<ipython-input-337-e135897d3660> in <listcomp>(.0)
8         forests, num = measurements.label(land)
9         area = measurements.sum(land, forests, index=np.arange(forests.max() + 1))
---> 10         LAND = np.array([land[i][j]*area[forests[i][j]] for i in range(L) for j in range(L)])
11         LAND = LAND.reshape(L,L)
12         fire = [LAND*probs]

KeyboardInterrupt: