Computational Modeling
Fall 2008

For today, you should:

1) work on Chapter 7 

2) prepare an evaluation packet:

   ONE of your chapters 4, 5, 6
   one copy of an evaluation rubric, filled out (just leave some room)

Today:

1) quiz 5 solutions

2) percolation

3) work on Chapter 7

For next time you should:

1) work on Chapter 7

Quiz 5 Solutions
----------------

1) In a 2-state totalistic CA with 9 neighbors, there are 10 possible
   states for the neighborhood (0 through 9), so there are 2^10
   possible rules.

   In a GoL-like CA, there are 9 possible states for the neighborhood,
   but there are two tables (one for each state of the center cell),
   so the total number of rules is 2^18.

2) Einstein seemed to believe that the Copenhagen Interpretation makes
   claims about how the world really is (deterministic or
   non-deterministic, etc.), and he had a strong opinion about that
   implication.  That suggests that he held a realist stance toward
   this theory.

   An instrumentalist would probably take the position that quantum
   mechanics is a useful model with no implications about how the
   world is.

3) The "microcanonical ensemble" is the set of values Q_n
   for 0 <= n <= N, where

Q_n = Prob { percolation | n porous cells } 

The "canonical ensemble" is a function Q(p) for 0 <= p <= 1 where

Q(p) = Prob { percolation | p = Prob {any cell is porous} }

Given the Q_n, you can compute Q(p):

Q(p) = sum_n p_n Q_n

Where p_n = Prob { n porous cells | p }

So p_n comes from the binomial distribution:
en.wikipedia.org/wiki/Binomial_distribution


Here is the simplest implementation of canonical_ensemble:

def choose(N, k):
    return factorial(N) / factorial(k) / factorial(N-k)

def canonical_ensemble1(Q, p):
    N = len(Q)-1
    sum = 0

    for n, qn in enumerate(Q):
        pn = choose(N, n) * p**n * (1-p)**(N-n)
        sum += pn * qn
        
    return sum

The only problem is that many implementations of factorial
are linear in N, which makes canonical_ensemble quadratic.

The easiest fix is to accelerate factorial using a memo:

factorial_memo = {0:1}

def factorial(n):
    try:
        return factorial_memo[n]
    except KeyError:
        return n * factorial(n-1)

Now factorial is amortized constant time and canonical_ensemble is O(n).

Other solutions:

1) precompute all the factorials and store them in a list.

2) compute the choose function using a recurrence relation:


def canonical_ensemble2(Q, p):
    N = len(Q)-1
    choose = 1
    sum = 0
    prod1 = 1
    prod2 = (1-p) ** N

    for n, qn in enumerate(Q):
        pn = choose * prod1 * prod2
        sum += pn * qn
        
        # update choose for the next iteration
        prod1 *= p
        prod2 /= (1-p)
        choose *= (N-n)
        choose /= (n+1)

    return sum

This solution has low overhead, but it is harder to prove correctness.


Percolation
-----------

My chapter glosses the difference between site and bond percolation,
but they are different enough that I should not have done that.

Bond percolation: sites represent locations, bonds represent
porous flow between sites.

1) Start with a lattice of sites, each a cluster of 1.

2) Add bonds between sites, joining clusters as you go.

3) Stop when the lattice starts percolating and record the number
   of bonds, x.

Site percolation: sites represent locations that may be porous.

1) Start with a grid of sites, all non-porous, and no clusters.

2) Switch a cell to porous, make it a cluster of 1, and join
   it to adjacent cluster(s), by some definition of "adjacent".

3) Stop when the grid starts percolating and record the number of
   porous cells, x.


How do Newman and Ziff join clusters efficiently?

What is their other major source of efficiency?


Suppose we run the N&Z algorithm 10000 times and record
the value of x where the lattice/grid starts to percolate.

How do we get from 10000 values of x to Q_n?

Q_n = Prob { percolation | n porous cells } 

And how do we get from the Q_n to Q(p)?