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: