Assignment 13: Monte Carlo Integration and Sampling Methods

Assignment 13: Monte Carlo Integration and Sampling Methods#

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm, trange
import timeit
from numba import jit
import copy
from ipywidgets import interact, interactive
import ipywidgets as widgets

Problem 1: Optimizing Speed of the 2D Ising Model Simulations#

In this problem, you will work on optimizing the execution time for the 2D Ising model simulations from the MonteCarloSamplingMethods lecture notebook. That code takes a few minutes to run for even a fairly small system of atoms. We want to optimize the code to make it as fast as possible so that we can repeat our results many times to make sure they are reliable, study larger systems, and use computational resources wisely.

Optimizing code requires two different types of analysis:

Complexity analysis. This approach considers the number of commands required to finish your simulation, and can be used to estimate the amount of time it takes to run your code before you run it. Complexity analysis leverages mathematical cleverness to solve a problem using the fewest number of steps possible. It does not depend on what programming language you use, or what computer your code will run on, which makes it a very powerful and flexible technique. For a good introduction, see here.
Execution optimization. This approach considers how to best use the hardware resources available to you, so it requires an understanding of both the different code libraries and languages you are using, and how the hardware uses these libraries to return your result. Significant improvements can be made, but it can be difficult to reduplicate the optimization for your algorithm on another machine or in a different programming language.

Execution optimization is more challenging than complexity analysis and is hardware specific, so in general you will want to start with complexity analysis, and then add execution optimization as the second optimization step.

Complexity Analysis#

To gain an understanding of complexity analysis, let’s start by looking at a function that defines just the MCMC step. Prepare to dust off your counting skills:

def MCMC_step(T: int, lattice: np.array):
    [rows, cols] = lattice.shape
    for r in range(1,rows-1):
        for c in range(1,cols-1):
            sum_NN = (lattice[r-1, c]+
                         lattice[r+1, c]+
                         lattice[r, c+1]+
                         lattice[r, c-1])
            E_a = -0.5*lattice[r,c]*sum_NN           
            E_b = -1*E_a
            if E_b < E_a:
                lattice[r, c] *= -1
            elif np.exp(-(E_b - E_a)*T) > np.random.rand():
               lattice[r, c] *= -1
   return lattice

We want to try and estimate how many commands happen inside this function. The actual number of commands is not as important as gaining intuition about what parts of the code take the most time to run. We want to find the slowest things first, so that we can fix them the fastest. Use the following guide lines to help with the analysis:

Work from inside nested code to outside nested code.
Corollary: Find the for loops first.

There are two for loops defined, on lines 3 and 4, and the for loop on line 4 is inside the for loop on line 3.

The for loop on line 4 repeatedly executes commands on lines 5-11. Since these commands have a constant number of steps, where the number of steps does not depend on the input arguments to the function, we can say that there are $C_1$ commands for lines 5-11, and the total number of commands $T$ for lines 5-11 is $$ \Large T_{5:11} = C_1 $$
Now, because there is a for loop on line 4, we know that the commands on lines 5-11 will be repeated for $x = (cols-1) - 1$ times, because that is the number of iterations of the line 4 for loop, where the index variable c goes through the different column values of the lattice. Furthermore, we can tell from line 2 that the value of cols depends on the input argument to this function. So when including the command on line 4, we will want to use a variable $x$ for the number of times this for loop executes, and the total number of commands becomes $$ \Large T_{4:11} = C_1 x $$
Next, we need to consider the for loop on line 3. For this notebook, we will only use square matrices, which means that the number of rows and columns will be the same, and the for loop on line 3 will execute $x$ times, the same as the loop on line 4. The total number of commands increases to $$ \Large T_{3:11} = C_1x^2 $$
Finally, we need to account for the command on line 2. This command is expected to have a constant run time because it is returning the shape of the input, and returning the shape requires the same calculations regardless of what the shape actually is. Therefore, we can say that the number of commands for line 2 is $C_2$. We can omit lines 1 and 12 because they initialize and exit the function, so the total time is $$ \Large T = C_1x^2 + C_2 $$

Here’s what we know from this analysis

The primary parameter controlling the speed of this algorithm is the size of the lattice that we are simulating. The $x$ value is the number of rows/columns in the system. The value of the other parameter, $T$, does not change how fast or slow the code is.
The slowest part of the code will be the nested for-loop. The nested for-loop structure is why we have the highest power term in $T$ is the $x^2$. This is a polynomial run time, which means that if we change from a $10\times10$ lattice to a $20\times20$ lattice, we should expect the time it takes to run the code to approximately quadruple. This can be shown mathematically as follows: $$ \Large T_1(x=10) = 100C_1 + C_2 $$

\[ \Large T_2(x=20) = 400C_1 + C_2 \]

\[ \Large T_2 \approx 4 T_1 \]

In general, good coding practices means you will repeat a command as few times as possible. DON’T put anything inside a for-loop (i.e. a variable initialization) unless it absolutely must be there, and be careful using for-loops in the first place.

This function is already very minimal, which means that to optimize it further we need to consider how it can be executed on the hardware.

Execution Optimization#

Before we get started optimizing, here are some things to know:

Computers convert every programming language to machine code in the end. However, the compiler a computer uses to do this conversion can make a huge difference in how quickly your code runs.
Python is an interpreted language, which means that every time a line of code is executed the computer pretends it has never compiled it before (even if it has, as in the case of a for loop). This means Python is easy to read and easy to change as you work, but the flexibility slows the computer down. Python was designed to be human-friendly, not computer-friendly.
The big Python data science libraries like NumPy, SciPy, SciKit-Learn, TensorFlow/Keras, and PyTorch have typically undergone execution optimization to make them as fast as possible. For example, NumPy uses the C programming language to speed up calculations; calculations actually aren’t executed in Python because C has a very good compiler. As a scientist, sometimes you will need decide whether to choose the speed gained by using a pre-optimized function, or the transparency that comes from coding your own version.
Python libraries for visualization like Matplotlib and Pandas will always be slower than libraries used for calculation, even if they have some calculation functionality included (cough Pandas cough). Keep this in mind when deciding where in your code to generate plots.

Execution Optimization Toolkit#

Here are some common tools useful for analyzing how your code runs. These tools focus on optimizing speed and memory; there are other tools for optimizing other aspects of programming like documentation and code consistency.

What you are optimizing	What you use to optimize it	How it works
Speed	timeit	Measures how long a code snippet takes to execute
Speed	tqdm	Measures how long a loop takes to execute with a loading bar!!!
Memory	memory-profiler	See which variables and programs use what memory resource

Note that these are diagnostic tools that provide a way to measure what the computer does with your code. It is still up to you to find solutions to slow parts or out-of-memory (OOM) errors and implement them.

Problem 1a) Complete the functions `intialize_lattice` and `MCMC_step`#

Use the information provided above and the commented hints in the function templates provided below. Feel free to copy code from above as you fill in the blanks, but be careful! Some variable names have been changed…

# Complete the initialize_lattice function below

def initialize_lattice(sqrt_n: int):
    """
    Function to initialize lattice.  Adds a border of zeros
    to represent non-interacting atoms and make the neighbor
    calculation easier

    sqrt_n: The square root of the number of atoms in the lattice
    returns padded_lattice: A lattice of size (sqrt_n+2, sqrt_n+2)
    """

    # initialize the lattice using np.random.random for a lattice of size
    # sqrt_n x sqrt_n

    new_lattice =

    # create a lattice of zeros that has an extra row on the top and bottom,
    # and an extra column on the left and the right

    padded_lattice =

    #mask lattice by setting values above 0.5 to 1, and everything else to -1
    new_lattice[new_lattice>0.5] =
    new_lattice[new_lattice!=1] =

    # added step to create non-interacting atoms
    padded_lattice[1:sqrt_n+1, 1:sqrt_n+1] =

    return np.array(padded_lattice)

# Complete the MCMC_step function below

def MCMC_step(beta: float, lattice: np.array):
    """
    Function to repeat the Monte Carlo Markov Chain for this system.
    beta: the inverse temperature value for the MCMC step
    lattice: the system of spins that will be simulated
    returns: an updated version of the input lattice
    """

    # Figure out the size of the lattice
    [rows, cols] = lattice.shape

    # keep the neighbors inside the region
    for r in range(1,rows-1):
        for c in range(1,cols-1):

            # sum over the nearest neighbors
            sum_NN =

            # calculate the energy
            E_a =

            # re-calculate the energy for a spin state change
            E_b = -1*E_a

            # choose whether to keep the new state or not
            if #<ENTER LOGIC STATEMENT HERE>
                lattice[r, c] *= -1

    return lattice

Problem 1b) Perform the following experiment#

By calling the main_exercise1b function using different arguments:

Re-run the simulation for lattice sizes of $\sqrt{N}$ = [10, 15, 20, 25, 40].
For each simulation of $\sqrt{N}$ values, record how long it took to run. Use timeit and/or tqdm to measure the how long the code takes to run. (What is the difference between the reported values if you use both?)
Make a plot showing how the time it took to run depends on the size of the simulation.
Make a plot showing the net-magnetization curves for each system size.

The function main_exercise1b is provided for you below; it calls the two functions you just completed.

def main_exercise1b(sqrt_N: int):
    """
    Main function to complete a simple Metropolis-Hastings Ising model
    experiment.

    sqrt_N: integer that is the square root of the number of atoms in the system

    returns [inverse_temperature, M]: where M is the net magnetization at each
            temperature
    """

    spin_lattice = initialize_lattice(sqrt_N)

    inverse_temperature = np.linspace(0, 2, 1000)

    M = [] #empty variable to hold net magnetism

    # For each temperature
    for beta in tqdm(inverse_temperature):

        # Repeat the MCMC step 100 times to make sure the system is stable
        for n in range(100):

            spin_lattice = MCMC_step(beta, spin_lattice)

        M.append(np.abs(np.sum(np.sum(spin_lattice)))/(sqrt_N*sqrt_N))

    return inverse_temperature, M

Call the main function in the code cell below like this.

start = timeit.default_timer()
beta_N5, M_N5 = main_exercise1b(5)
print("\nTime to completion for N = 5: "+str(timeit.default_timer() - start))

Make sure to keep the beta and M values for plotting.

# Enter your code to call the main function here.  You can either
# call the function multiple times in one cell, or add code cells
# for each call

# Manually enter your results for each run here:
num_atoms = []
measured_time = []

# Put your plotting code here

Problem 1c) Optimize the code to see if you can reduce the time it takes to run#

Copy-paste the MCMC_step function from above into the empty code cell below. We are going to modify it to make it faster, but we also want to keep the original for reference. Rename the function MCMC_step_optimized.
Add a numba “decorator” to the line just above the function declaration. It should look like this:

@jit()
def MCMC_step_optimized(beta: float, lattice: np.array):
...

The @jit() decorator comes from the numba python package and, when placed before the function, tells the Python compiler to treat this code like a statically compiled C object. The first time the function is called, it won’t show any improvement because it isn’t compiled yet. However, the second time the function is called, the compiled version should still be in memory, and we will see some improvements.

# Copy-paste the MCMC_step function from above here

A slightly altered main function main_exercise1c is provided for you here. No need to change anything. Just call it the same way you did before.

start = timeit.default_timer()
beta_N5, M_N5 = main_exercise1c(5)
print("\nTime to completion for N = 5: "+str(timeit.default_timer() - start))

#from tqdm.gui import tqdm_gui
# Here is a slightly altered main function for you

def main_exercise1c(sqrt_N: int):
    """
    Main function to complete a simple Metropolis-Hastings Ising model
    experiment.

    sqrt_N: integer that is the square root of the number of atoms in the system

    returns [inverse_temperature, M]: where M is the net magnetization at each
            temperature
    """

    spin_lattice = initialize_lattice(sqrt_N)

    inverse_temperature = np.linspace(0, 2, 1000)

    M = np.zeros(1000) #empty variable to hold net magnetism

    # For each temperature
    for t in range(len(inverse_temperature)):
        beta = inverse_temperature[t]

        # Repeat the MCMC step 500 times to make sure the system is stable
        for n in range(500):

            spin_lattice = MCMC_step_optimized(beta, spin_lattice)

        M[t] = (np.abs(np.sum(np.sum(spin_lattice)))/(sqrt_N*sqrt_N))

    return inverse_temperature, M

Repeat the experiment with the same values for $\sqrt{N}$ from above, and record the new time values. Create a two more plots to compare with the ones created above (Time vs $\sqrt{N}$, and Net Magnetization vs Temp), and compare them with the plot from earlier.

# Call the optimized version of the functions here:

# Create plotting code here

Answer the following questions:

What things are the same in the Net Magnetization plot?

What things are different in the Timing plot?

How effective is using the numba package at improving your code?

Problem 2: Study of Magnetic Domain Formation#

Now that we have code up and running for a large system efficiently, let’s watch some magnetic domains form. Complete the following experiment:

Write a new main function main_exercise2 by making the following changes to main_exercise1c:
- create a new list variable: lattice_at_T. This variable will hold a set of lattices at different temperatures
- update the lattice_at_T variable after the net magnetism M is calculated by appending the spin_lattice to the list.
- return the lattice_at_T variable along with the inverse_temperature and M variables
- A template is started for you below.
Choose the largest value of N that you can simulate within a reasonable amount of time. Use the code in Exercise 1C to figure out how large to go.
Execute the new main_exercise2 function. Save all three outputs.
Plotting code is created for you below.

def main_exercise2(sqrt_N: int):
    """
    Main function to complete a simple Metropolis-Hastings Ising model
    experiment.

    sqrt_N: integer that is the square root of the number of atoms in the system

    returns [inverse_temperature, M]: where M is the net magnetization at each
            temperature
    """

    # YOUR CODE HERE
    raise NotImplementedError()

beta, mag, lattices = main_exercise2(100)

fig, ax = plt.subplots(1, 1)
ax.plot(beta, mag)
ax.set_xlabel(r"Inverse Temperature $\beta$")
ax.set_ylabel("Net-magnetism M")
plt.show()

def plot_mag_domains(x):
    plt.figure(2)
    plt.imshow(lattices[x])
    plt.title(r"$\beta$ = "+str(beta[x]))
    plt.show()

interactive_plot = interactive(plot_mag_domains, x=(0,999))
output = interactive_plot.children[-1]
output.layout.height = '500px'
interactive_plot

Answer the following questions:

Explore the set of simulated lattices above using the slider, and see if you can detect changes in the magnetic domains generated. Edit this cell with your answers to the questions.

At what inverse temperature $\beta < 1$ do you notice a significant change in the lattice domain states?

How does this compare to the plotted values?

How would you describe the shapes of the domains that form at $\beta \in [0.75, 0.85]$?

Do you have a lattice that is completely aligned? For what $\beta$ value? Is there one that is close?

Problem 3: Monte Carlo Integration#

In this problem, you will study the accuracy of Monte Carlo integration in each of four different expressions, each with some physical significance, shown in the table below:

Expression #	Function	Interval	Notes
1.	${x = \int_0^1 7-10t\ dt} $	$t$ is time; $t = [0, 1]$ s	Gives position at time $t$ for this system
2.	${\Delta S}$ ${= \int_{300}^{400}\frac{mc}{T}\ dT }$	$m$ is mass; $m=1$ kg $c$ is specific heat capacity; $c = 4190$ J/kg K $T$ is temperature; $T = [300, 400]$ K	Change in entropy for thermal processes
3.	$\Phi = \int_1^2 \frac{Q}{4 \pi \epsilon_o r^2} dr$	$r$ is distance; $r = [1, 2]$ m $\epsilon_o$ is the Permittivity of Free Space $Q$ is the charge; $Q = 1$ C	$\Phi$ is the electric potential energy gained by moving along line $r$
4.	$I = \int_0^\infty \frac{2 \pi h c^2}{\lambda^5(e^{hc/\lambda k T} - 1)}\ d\lambda$	$h$ is Planck’s constant $c$ is speed of light $k$ is Boltzmann’s Constant $T$ is the absolute temperature; T = 400K $\lambda$ is wavelength; $\lambda = [0, \infty]$ m	Planck’s radiation law; Integrating gives Stefan Boltzmann Law

Analytically integrate each for the region and values provided, and record your answer in the analytical_result variables below:

analytical_result_expr1 = None # replace the None's with your results
analytical_result_expr2 = None
analytical_result_expr3 = None
analytical_result_expr4 = None

Show your work in the cell below, either in a picture file for written derivations or in Latex

Write the each expression to be integrated as a python function.

# Helpful constants:
pi = np.pi #unitless
c = 2.99E8 #m/s
h = 6.62607015E-34 #J
k = 1.380649E-23 #J/K
epsilon = 8.854187817E-12 #F/m
sigma = 5.6704E-8 #W/(m^2 K^4)

def integrand(x):
    # YOUR CODE HERE
    raise NotImplementedError()

Randomly choose values for x from a distribution between the limits of the definite integral.

lower_limit = # this can be a float
upper_limit = # this can be a float
num_x_vals  = # this must be an integer value less than or equal to 10^8
x_vals = np.random.uniform(lower_limit, upper_limit, num_x_vals)

Calculate the f_x values:

y = integrand(x_vals)

approx_int = ((upper_limit - lower_limit)*np.sum(y))/(num_x_vals - 1)
print(f"The Monte Carlo approximation is {approx_int}")

Calculate the error between the approx_int and the analytical_result variables using one or more of the metrics discussed above

mse = None # replace with your calculation
print(f"The Mean Squared Error is {mse}")

pe = None # replace with your calculation
print(f"The Percent Error is {pe}%")

Finally, we want to visualize how the error decreases as the number of random trials num_x_vals increases. Write code to the do the following:

Using the error metric you decided on above, write a for-loop that calculates the error as a function of the number of points you sampled. For example, calculate the error when you summed two values of $\langle F^N \rangle$, then calculate the error for three summed values of $\langle F^N \rangle$, and so on until you have calculated the errors for the full range of $\langle F^N \rangle$.
IMPORTANT: You do not need to re-do the experiment to calculate this analysis; if you do it will slow down your for-loop and potentially crash your notebook kernel. Instead, you will want to reuse all of the integrand values are stored in the y variable. Python indexing into this list using the y[:N] functionality will give you the first N values in this list. The first N values can then be used to calculate a $\langle F^N \rangle$ value for the first N samples.
Make a figure showing how the error changes with the number of values in the sum.

error_data = [] 
# Write code here to fill error_data with the percentage error corresponding to each of the number of points you sampled in the MC integration

Finally, plot it

plt.plot(np.linspace(2, 2000000, 1999998, endpoint=True), error_data)
plt.xlabel("Number of Values in Sum")
plt.ylabel("Percent Error")

Answer the following questions:

Model vs Simulation: In your own words, describe the difference between a model and a simulation. Give your own example of a model, and how you would simulate it.

Markov Chain: In your own words, describe a Markov Chain and its properties. Give your own example of a stochastic system and how you would implement a Markov Chain for it.

Expression #	Function	Interval	Notes
1.	\({x = \int_0^1 7-10t\ dt} \)	\(t\) is time; \(t = [0, 1]\) s	Gives position at time \(t\) for this system
2.	\({\Delta S}\) \({= \int_{300}^{400}\frac{mc}{T}\ dT }\)	\(m\) is mass; \(m=1\) kg \(c\) is specific heat capacity; \(c = 4190\) J/kg K \(T\) is temperature; \(T = [300, 400]\) K	Change in entropy for thermal processes
3.	\(\Phi = \int_1^2 \frac{Q}{4 \pi \epsilon_o r^2} dr\)	\(r\) is distance; \(r = [1, 2]\) m \(\epsilon_o\) is the Permittivity of Free Space \(Q\) is the charge; \(Q = 1\) C	\(\Phi\) is the electric potential energy gained by moving along line \(r\)
4.	\(I = \int_0^\infty \frac{2 \pi h c^2}{\lambda^5(e^{hc/\lambda k T} - 1)}\ d\lambda\)	\(h\) is Planck’s constant \(c\) is speed of light \(k\) is Boltzmann’s Constant \(T\) is the absolute temperature; T = 400K \(\lambda\) is wavelength; \(\lambda = [0, \infty]\) m	Planck’s radiation law; Integrating gives Stefan Boltzmann Law

Assignment 13: Monte Carlo Integration and Sampling Methods

Contents

Assignment 13: Monte Carlo Integration and Sampling Methods#

Problem 1: Optimizing Speed of the 2D Ising Model Simulations#

Complexity Analysis#

Execution Optimization#

Execution Optimization Toolkit#

Problem 1a) Complete the functions intialize_lattice and MCMC_step#

Problem 1b) Perform the following experiment#

Problem 1c) Optimize the code to see if you can reduce the time it takes to run#

Problem 2: Study of Magnetic Domain Formation#

Problem 3: Monte Carlo Integration#

Problem 1a) Complete the functions `intialize_lattice` and `MCMC_step`#