#!/usr/bin/python3
import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt

C = 0.9652  # Normalizing constant for f

def f(x):
    y = np.cos(50 * x) + np.sin(20 * x)
    return y * y 

def mutate(x):
    z = rnd.rand()
    if z < 0.5:
        return z * 2
    else:
        y = x + rnd.normal(0, 1)
        # Wrap around
        y = y - int(y)
        if y < 0:
            y = 1 + y
        return y

def metropolis():
    x = 0.5           # Initial state
    N = 100000          # Number of samples
    X = np.empty(N)   # Array of samples

    for i in range(0, N):
        # Propose a new sample
        y = mutate(x)

        # Acceptance probability
        a = min(1, f(y) / f(x)) # Symmetric proposal, proposal densities cancel out

        if rnd.rand() < a:
            x = y

        # Store the sample
        X[i] = x

    # Plot the normalized version of f
    D = np.arange(0, 1, 0.005)
    Y = np.vectorize(lambda x: f(x) / C)(D)
    # Plot the histogram of the samples we got from the algorithm
    plt.plot(D, Y, "r-")
    plt.hist(X, bins = 200, normed=True)

def main():
    metropolis()
    plt.show()

if __name__ == "__main__":
    main()