Reading: Percolation handout II
Today we will consider the case of undirected percolation and a recursive approach to solving this problem. For Part 1 of your homework you constructed a vertical_flow function that took as input a grid of blocked (=0) or vacant (=1) sites and returned a flow matrix of empty (=0) or full (=1) sites. In the vertical percolation case we copied the values of sites to flow for the first row. For each site not on the first row we simply checked to see if the corresponding site was vacant in sites and if the site above it in flow was full. If we now allow the flow to move vertically and horizontally into vacant sites we need another approach.
Suppose you are looking at a grid of blocked and vacant sites. How would you determine which sites get filled up in the undirected percoalation model? One natural way for us humans to approach the problem is to start at a site in the top row and follow all possible paths of vacant cells from there. As we find the paths, we fill up the vacant cells that make them up. Then we move to the next site on the top row and do the same, filling up any new paths of vacant cells that originate there. And so on.
Okay, this sounds pretty good. Let's try and formalize things a little bit. When we say follow all possible paths of vacant cells how exactly do we do that? In general there are two flavors of searching for these kinds of paths. One is to explore every possible path one step at a time. This approach is called Breadth-first search. The other is to go as far as you can one way, then backtrack to an earlier a decision point and follow another path as far as it goes and then backwards and so on. This approach is called Depth-first search. Both are good, we'll use a depth-first approach.
A Recursive Approach: Let's step back and look at the big picture for a moment. We want to write a function that given an nxn matrix of vacant/blocked sites we'll call sites, will produce an nxn matrix of empty/full sites we'll call full. Okay, so first off lets construct an nxn grid of empty sites called full. Now we have to fill it appropriately. Suppose we are the flow and we're starting from position (0,0) in the empty matrix full and we wish to explore all possible paths from here. What do we need to consider:
sites even vacant? If it's blocked we can't do anything so stop.full Great. How do we check to see if we can move down? Well, we look at that site and see if it's empty in sites and if so we fill it in full and check to see where we can go from there. Notice this is similar to what we just did. That is, we're starting at some specific site and we're seeing if we can fill it and if so, where can we go from there. So lets rewrite the steps above more generally beginning a some arbitrary site (i,j).
sites? If so, stop.fullObserve that since we are considering an arbitrary site (i,j) we will pretend we have no idea how we got here and therefore we must check to see if we can exit the site in all four directions. So we can rewrite the above as:
sites?fullFurthermore, if it's really an arbitrary i and j we need to check if they're even in the grid. That is, is (i,j) even within the bounds of the grid indices? So now we have
sites matrix blocked?fullAnd finally one more special case to consider. What if we get to site i,j and we have already visited it and filled it? Well in that case we don't want to do anything either. We're ready to rewrite the above for the final time now.
sites blocked?fullHow is this recursive? Well, notice that steps B through D are really just calls to this same algorithm on slightly different (i,j). Looking down is checking (i+1,j). Looking right is checking (i, j+1). Looking left is checking (i,j-1). Looking up is, you guessed it, checking (i-1,j). The base case is really just steps 1 through 4. Think about why in a finite grid we will always end up in a situation that is treated by one steps 1-4. That is, a situation where we just stop.
So if we call this algorithm on each site in the top row of the empty matrix full, at the end of the last call, full will accurately represent the flow due to undirected percolation. Let's call the algorithm above flow_from. What would the necessary inputs be for flow_from? Well we need i and j and sites and full. To write the more general function undirected_percolaton, we simply have to make the empty full matrix and then call flow_from on each site in the top row. Notice that for this to work, flow_from will modify full. That's okay. This is one of those times when a side effect is okay because it's only visible to the caller of a support function. That is, the side effect is due to the support function flow_from. To the caller of undirected_percolation there is no visible side effect since the matrix full is created inside undirected_percolation. So as far as side effects go, this isn't a bad one.
So I know you may be confused by all of this. To clear things up, stop. Imagine once again that you are figuring out the flow and the way you're doing it is by starting at each site in the top row, following the flow for as long as it can go and filling up the new full matrix as you go. All we have done above is an attempt to articulate this a little more formally. It's now your job to go from here to a Python implementation. Good luck!
Flow Visualization: To visualize the flow (and to ultimately create your plot) you can use a module from the matplotlib library. The module is pyplot. Here's how you use it.
import numpy as np
import matplotlib.pyplot as plt
#lets make a random 20 by 20 matrix
a=np.random.rand(20,20)
#lets make a matrix of zeros and ones from this
b=a<.6
b=b.astype(int) # converts from bools to ints
#lets make another matrix of zeros and ones (with fewer ones)
c=a<.3
c=c.astype(int) # converts from bools to ints
d=b+c # let d be the element-wise sum
# now to visualize
plt.matshow(d)
plt.show()
What just happened?
So using the pyplot function matshow combined with show we can visualize a matrix. Investigate the online documentation for pyplot and figure out how to customize the colors, add titles, and more.
Plotting with matplotlib: We use the pyplot module to generate plots as well. All that is required are two numpy arrays of the same size to plot against each other. Here are some simple examples.
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 3, 20)
y = np.linspace(0, 9, 20)
plt.plot(x, y) # line plot
#plt.plot(x, y, 'o') # dot plot
plt.show() # shows the plot
What just happened? We created two 20 point arrays using np.linspace and then plotted them.
import numpy as np
import matplotlib.pyplot as plt
def f(t):
return t**2*np.exp(-t**2)
t=np.linspace(0,3,51)
y= f(t) # makes a numpy array
plt.plot(t,y)
plt.show()
What just happened? This time we created a function. We used it to populate another array y in a single step!
Summary: For your assignment you need to generate a plot of estimated percolation probability against site vacancy probability. To do this you will generate a random site vacancy matrix using site vacancy probability p (for example p=.1). Then you'll check to see if it percolates. Repeat this many (for example 5000) times for a single value of p. Use the number of times it percolates divided by the number of trials as your estimate for the percoaltion probability. Do this for many values of site vacancy probability p (for example 30 different values between 0 and 1). At the end of this process you have 30 values of estimated percolation probability and the 30 corresponding site vacancy probabilities. Put these values in arrays and plot them!
Consider the function definition
def f(a,b):
print('a is: ',a)
print('b is: ',b)
return
f(2,3)
As expected the first parameter was passed into the parameter a and the second parameter into b. Now observe what happens if we invoke the function in this way:
f(b=10,a=0)
What just happened? This time we specified the parameters by keyword, not by position. Of course we needed to know what the formal parameter names were but if we do, as long as we specifiy them by name, we don't have to worry about order. The only rule is that positional arguments come before keyword arguments in the function call. So for example this is okay:
f(6,b=12)
but this is not okay:
f(b=12,6)
What about f(6,a=12)? Try it! You'll see that this doesn't work because the argument 6 is passed into the parameter a and then we try to reassign it to 12.
Here's another useful quality of keyword parameters:
def g(a,b=12,c=0):
print('a is: ', a)
print('b is: ', b)
print('c is: ', c)
return
g(6,20)
What just happened?: We can provide default values to parameters in a function defintion. If the caller specifies a value either positionally or via keyword it will override the default value (like b in the example above). If the caller does not provide an argument for parameter with a default (like c above) then the default value will be used.
Summary so far: The parameters we've seen so far are called positional-or-keyword parameters. This is because they may be used as either depending on the structure of the function call. In general there are five catagories of function parameters. We'll learn about all of them today.
We've seen these before. In a function defintion these must come after the positional or keyword parameters. So for example:
def f(a,b=0,*args):
print('a is: ',a)
print('b is: ',b)
print('args is: ',args)
return
f(3,2,4,5,6)
The * in front of the parameter args makes it a variable-length-positional parameter. This means it absorbs any extra postional arguments and places them into a tuple called args. Specifying a parameter like this in a function definition gives the function caller the option of providing as many postional arguments as they want without breaking the function.
We can specifiy additional keyword parameters after the variable-length-positional parameter in a function definition. These parameters may only be specified by keyword, never positionally. For this reason we call these keyword-only parameters. Like this:
def g(a,b=0,*args,c=1,d=2):
print('a is: ',a)
print('b is: ',b)
print('args is: ',args)
print('c is: ',c)
print('d is: ',d)
return
g(1,2,3,4,5,d=100)
Play around witht he IPython interpreter try different combinations to see what works. The main rule to remember is that if the parameter comes after the variable length positional paramenter then it is keyword only. Notice that this makes sense since the variable length parameter (typcially but not necessarily called args) will absorb all positional arguments beyond what's been specified before it.
Our fourth catagory of function paramenter absorbs any excess keyword arguments. So just like the variable-length-postional parameter absorbs extra postional arguments into a tuple, this parameter will absorb extra keyword paramenters. Notice If we consider what kind datatype we should use to store an unknown quantity of keyword arguments, it shouldn't take us long to arrive at a dictionary. The keyword argument names will be the keys and their values will be the values. Observe:
def g(a,b=0,*args,c=1,d=2, **kwargs):
print('a is: ',a)
print('b is: ',b)
print('args is: ',args)
print('c is: ',c)
print('d is: ',d)
print('kwargs is: ',kwargs)
return
g(1,2,3,4,5,d=100, e=200, f=300)
What just happened? We used the sequence ** before the name kwargs to tell Python that the parameter kwargs is variable length keyword. It then absorbed any excess keyword arguments into the dictionary named kwargs. Note that there is nothing special about the names args or kwargs. We could have called these parameters anything. It's convention to call them args and kwargs though.
You may have guessed that the fifth and final category of parameter is positional-only. However, we only mention these for completeness sake. We as programmers cannot construct postional-only parameters. Certain built-in functions do have them (for example abs) but that is the extent of their appearance in Python.