app.py

import numpy as np
from pyDOE2.doe_lhs import lhs
import matplotlib.pyplot as plt
import streamlit as st
import mpmath

# Set the precision to 50 decimal places
mpmath.mp.dps = 50

# Compute the value of pi with 50 decimal places of precision
pi = mpmath.pi

st.title("Compute Pi with Monte Carlo")

st.markdown(
    r"""
    We know that area of a circle is $\pi r^2$. If we enclose the circle in a square of side $2r$, then the area of the square is $4r^2$. The ratio of the area of the circle to the area of the square is $\frac{\pi r^2}{4r^2} = \frac{\pi}{4}$. If we randomly sample points in the square, then the ratio of the number of points in the circle to the total number of points is also $\frac{\pi}{4}$.
            """
)


# draw a circle with matplotlib
def draw_the_plot():
    fig, ax = plt.subplots()
    ax.set_aspect("equal")
    # draw a square with matplotlib
    square = plt.Rectangle(
        (-1, -1), 2, 2, facecolor="lightblue", edgecolor="r", linewidth=2
    )
    ax.add_artist(square)
    # add a two sides arrow outside the square to indicate the length of the side
    upper_margin = 1.05
    ax.arrow(
        -1,
        upper_margin,
        2,
        0,
        head_width=0.05,
        color="k",
        length_includes_head=True,
        clip_on=False,
    )
    ax.arrow(
        1,
        upper_margin,
        -2,
        0,
        head_width=0.05,
        color="k",
        length_includes_head=True,
        clip_on=False,
    )
    # annotate the arrow
    ax.annotate(
        "$2r$",
        xy=(0, upper_margin),
        xytext=(0, upper_margin),
        horizontalalignment="center",
        verticalalignment="bottom",
    )
    circle = plt.Circle((0, 0), 1, facecolor="lightgreen", edgecolor="b", linewidth=2)
    ax.add_artist(circle)
    # draw a horizontal arrow from the center to the edge
    ax.arrow(0, 0, 1, 0, head_width=0.05, color="k", length_includes_head=True)
    # annotate the arrow
    ax.annotate(
        "$r$",
        xy=(0.5, 0),
        xytext=(0.5, -0.1),
        horizontalalignment="center",
        verticalalignment="bottom",
    )

    ax.set_xlim(-1.1, 1.1)
    ax.set_ylim(-1.1, 1.1)
    ax.set_xticks([])
    ax.set_yticks([])
    # remove the x and y axis
    ax.spines["left"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_visible(False)
    ax.spines["bottom"].set_visible(False)

    return fig, ax


fig, ax = draw_the_plot()
st.pyplot(fig)

st.subheader("Monte Carlo Simulation")

seed = st.number_input("Random seed", min_value=0, max_value=10000000, value=0, step=1)
N = st.number_input(
    "Number of points", min_value=1, max_value=10000000, value=10000, step=1
)
type = st.selectbox("Type of sampling", ["Random", "Grid"])

if type == "Random":
    samples = lhs(2, samples=N, random_state=seed) * 2 - 1
else:
    samples = np.linspace(-1, 1, int(np.sqrt(N)))
    samples = np.array(np.meshgrid(samples, samples)).T.reshape(-1, 2)

# compute the distance of each point from the origin
distances = np.linalg.norm(samples, axis=1)

# compute the number of points inside the circle
inside = samples[distances < 1, :]
print(samples.shape, inside.shape)

st.markdown(
    f"""
    Points inside the circle: = {len(inside)}\n
    Ratio: = {len(inside)} / {len(samples)} = {len(inside) / len(samples)}\n
    $\pi$ = {len(inside) / len(samples)} * 4\n
    | $\pi$ | value |
    | --- | --- |
    | estimated $\pi$ | {len(inside) / len(samples) * 4:.50f} |
    | real $\pi$ | {pi} |
    """
)

fig, ax = draw_the_plot()
ax.scatter(samples[:, 0], samples[:, 1], s=0.1, color="r")
ax.scatter(inside[:, 0], inside[:, 1], s=0.1, color="b")
st.pyplot(fig)