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:
1. def MCMC_step(T: int, lattice: np.array):
2. [rows, cols] = lattice.shape
3. for r in range(1,rows1):
4. for c in range(1,cols1):
5. sum_NN = (lattice[r1, c]+
lattice[r+1, c]+
lattice[r, c+1]+
lattice[r, c1])
6. E_a = 0.5*lattice[r,c]*sum_NN
7. E_b = 1*E_a
8. if E_b < E_a:
9. lattice[r, c] *= 1
10. elif np.exp((E_b  E_a)*T) > np.random.rand():
11. lattice[r, c] *= 1
12. 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 511. 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 511, and the total number of commands \(T\) for lines 511 is $\( \Large T_{5:11} = C_1 \)$Now, because there is a
for
loop on line 4, we know that the commands on lines 511 will be repeated for \(x = (cols1)  1\) times, because that is the number of iterations of the line 4for
loop, where the index variablec
goes through the different column values of the lattice. Furthermore, we can tell from line 2 that the value ofcols
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 thisfor
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 thefor
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 forloop. The nested forloop 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 \)$
In general, good coding practices means you will repeat a command as few times as possible. DON’T put anything inside a forloop (i.e. a variable initialization) unless it absolutely must be there, and be careful using forloops 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 humanfriendly, not computerfriendly.The big Python data science libraries like NumPy, SciPy, SciKitLearn, 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 preoptimized 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 
Measures how long a code snippet takes to execute 

Speed 
Measures how long a loop takes to execute with a loading bar!!! 

Memory 
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 outofmemory (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 noninteracting 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 noninteracting 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,rows1):
for c in range(1,cols1):
# sum over the nearest neighbors
sum_NN =
# calculate the energy
E_a =
# recalculate 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:
Rerun 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/ortqdm
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 netmagnetization curves for each system size.
The function main_exercise1b
is provided for you below; it calls the two functions you just completed.
Note: This computation may take some time to complete. Discuss with your instructor how long you should wait for a simulation to run.
def main_exercise1b(sqrt_N: int):
"""
Main function to complete a simple MetropolisHastings 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#
Copypaste 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 functionMCMC_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.
# Copypaste 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 MetropolisHastings 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?
Answer:
What things are different in the Timing plot?
Answer:
How effective is using the
numba
package at improving your code?
Answer:
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 tomain_exercise1c
:create a new
list
variable:lattice_at_T
. This variable will hold a set of lattices at different temperaturesupdate the
lattice_at_T
variable after the net magnetismM
is calculated by appending thespin_lattice
to the list.return the
lattice_at_T
variable along with theinverse_temperature
andM
variablesA 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 MetropolisHastings 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("Netmagnetism 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?
Answer:
How does this compare to the plotted values?
Answer:
How would you describe the shapes of the domains that form at \(\beta \in [0.75, 0.85]\)?
Answer:
Do you have a lattice that is completely aligned? For what \(\beta\) value? Is there one that is close?
Answer:
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 710t\ dt} \) 
\(t\) is time; \(t = [0, 1]\) s 
Gives position at 
2. 
\({\Delta S}\) \({= \int_{300}^{400}\frac{mc}{T}\ dT }\) 
\(m\) is mass; \(m=1\) kg 
Change in entropy for 
3. 
\(\Phi = \int_1^2 \frac{Q}{4 \pi \epsilon_o r^2} dr\) 
\(r\) is distance; \(r = [1, 2]\) m 
\(\Phi\) is the electric potential energy 
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 
Planck’s radiation 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.
For example, if I want to integrate the expression
\[ \Large > F = \int 3x^2 + 17\ dx > \]then my integrand is
\[ \Large > f(x) = 3x^2 + 17 > \]and I would write the following function:
def integrand(x): f_x = 3*np.power(x, 2) + 17 return f_xThis function takes
x
as my function argument, and returns the calculated valuef_x
. Note that I am not yet evaluating the limits of my integrand.
# Helpful constants:
pi = np.pi #unitless
c = 2.99E8 #m/s
h = 6.62607015E34 #J
k = 1.380649E23 #J/K
epsilon = 8.854187817E12 #F/m
sigma = 5.6704E8 #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.
Hint: if one of your limits is \(\infty\), it is okay to approximate it with a large number. Another way to do it is to plot [x, f(x)] and visually estimate the most important region of your integration.
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 forloop 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 redo the experiment to calculate this analysis; if you do it will slow down your forloop 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 they[:N]
functionality will give you the firstN
values in this list. The firstN
values can then be used to calculate a \(\langle F^N \rangle\) value for the firstN
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.
Answer:
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.
Answer: