from math import gcd
import numpy as np
import matplotlib.pyplot as plt

## Cantor staircase

def cantor(n):
    return [0.0] + cant(0.0, 1.0, n) + [1.0]


def cant(x, y, n):
    if n == 0:
        return []

    new_pts = [2.0 * x / 3.0 + y / 3.0, x / 3.0 + 2.0 * y / 3.0]
    return cant(x, new_pts[0], n - 1) + new_pts + cant(new_pts[1], y, n - 1)


xcantor = np.array(cantor(5))
ycantor = np.cumsum(np.ones(len(xcantor)) / (len(xcantor) - 2)) - 1.0 / (len(xcantor) - 2)
ycantor[-1] = 1

# Unit circle
tcircle = np.linspace(0, 2 * np.pi, 100)

## Brownian motion
nbrow = 300
dt = 1.0 / (nbrow - 1)
brow = np.cumsum(np.sqrt(dt) * np.random.randn(nbrow))

## Weierstrass function
def weirstrass(x, a, b):
    n = 200  # number of terms in the series
    result = 0.0
    for k in range(n):
        result += a**k * np.cos(b**k * np.pi * x)
    return result

## Thomae's function
xthom = []
ythom = []
for q in range(1, 101):
    for p in range(1, q):
        xthom.append(-p * 1.0 / q * 1.0)
        ythom.append(-1 / (q / gcd(p, q)))

fig, axs = plt.subplots(1, 4, figsize=(20, 3))

axs[0].plot(np.arange(nbrow)/nbrow, brow, color="k")
#axs[0].set_title("Realization of 1D brownian motion")
axs[1].plot(xcantor, ycantor, color="k")
#axs[1].set_title("Cantor staircase graph")
axs[2].plot(np.linspace(0, 1, 300), weirstrass(np.linspace(-1, 1, 300), 0.5, 3.0), color="k")
#axs[2].set_title("Weierstrass function")
axs[3].scatter(xthom, ythom, color="k", marker='x')
#axs[3].set_title("Negative Thomae's popcorn function")
fig.tight_layout()

#plt.subplots_adjust(wspace=0.6)
plt.savefig("lower_semicontinuous_functions.png")