import math
import matplotlib.pyplot as plt
import numpy as np
import random
import time


############################################################
# RANDOM SEED
############################################################


# seed = time.time_ns()
# seed = 12345


def set_seed(new_seed):
    global seed # tells Python to use the variable outside of the function
    seed = new_seed


def random_nb():
    global seed
    a = 1664525
    c = 1013904223
    m = 2**32
    seed = (a * seed + c) % m
    return seed / m


# for i in range(4):
#      print(random_nb())


def random_int_range(start=0, end=1):
    r = random_nb()  # base random number between 0 and 1
    return int(start + r * (end - start))


def random_poisson(lmbda):
    """Generate a random integer following a Poisson distribution (Knuth algorithm)."""
    L = math.exp(-lmbda)
    k = 0
    p = 1.0

    while True:
        k += 1
        p *= random_nb()
        if p <= L:
            break
    return k - 1


# for i in range(4):
#      print(random_int_range(10, 100))
#      print(random_poisson(10))

# for i in range(4):
#     print(random.randint(1, 100))

# random.seed(0)
# for i in range(4):
#     print(random.randint(1, 100))

# for i in range(4):
#     random.seed(0)
#     print(random.randint(1, 100))


# val = 0
# for _ in range(3):
#     random.seed(0)
#     val += random.random()
#     val -= random.random()
# print(val == 0)

# val = 0
# for _ in range(3):
#     random.seed(0)
#     val += random.random()
#     random.seed(0)
#     val -= random.random()
# print(val == 0)


############################################################
# ROLLING DICE
############################################################


def roll_die():
    """Return an int between 1 and 6"""
    return 2


def roll_die():
    """Return a random int between 1 and 6"""
    return random.choice([1, 2, 3, 4, 5, 6])


def test_roll(n):
    result = ''
    for _ in range(n):
        result += ' ' + str(roll_die())
    return result


# for _ in range(10):
#     print(test_roll(5))


def run_sim(goal, num_trials):
    total = 0
    for i in range(num_trials):
        if i != 0 and i % 100_000 == 0:
            print(f"Starting trial {i}")
        result = ''
        for _ in range(len(goal)):
            result += str(roll_die())
        if result == goal:
            total += 1
    print(f"Actual probability of {goal} = {1 / 6**len(goal):.8f}")
    print(f"Estimated Probability of {goal} = {total / num_trials:.8f}")


# run_sim('11111', 1000)
# run_sim('11111', 1_000_000)


############################################################
# Modelling a Business
############################################################


FIXED_COST = 10_000  # fixed overhead cost in dollars
random.seed(3333)


def run_simulation():
    price = np.random.uniform(10, 20)
    cost = np.random.uniform(5, 15)
    demand = max(np.random.normal(10000, 1000), 0)

    revenue = price * demand
    variable_profit = (price - cost) * demand
    total_profit = variable_profit - FIXED_COST

    return revenue, total_profit


def likelihood_of_loss(num_simulations=100_000):
    losses = 0
    for _ in range(num_simulations):
        _, profit = run_simulation()
        if profit < 0:
            losses += 1
    return losses / num_simulations


# print("Likelihood of a loss", likelihood_of_loss())


def plot_distribution():
    # --- Create distribution ---
    N = 100_000  # number of simulation runs
    revenues = []
    profits = []
    for _ in range(N):
        rev, prof = run_simulation()
        revenues.append(rev)
        profits.append(prof)
    revenues = np.array(revenues)
    profits = np.array(profits)

    # --- Print summary statistics ---
    print(f"Expected revenue: ${np.mean(revenues):,.2f}")
    print(f"Revenue std dev:  ${np.std(revenues, ddof=1):,.2f}")
    print(f"Expected profit (after $10k fixed cost): ${np.mean(profits):,.2f}")
    print(f"Profit std dev:   ${np.std(profits, ddof=1):,.2f}")

    # --- Plot distributions ---
    plt.figure(figsize=(8,5))
    plt.hist(revenues, bins=60, alpha=0.8, label="Revenue")
    plt.title("Monte Carlo Simulation: Revenue Distribution")
    plt.xlabel("Revenue ($)")
    plt.ylabel("Frequency")
    plt.legend()
    plt.tight_layout()
    plt.show()

    plt.figure(figsize=(8,5))
    plt.hist(profits, bins=60, alpha=0.8, color="orange", label="Profit (after fixed cost)")
    plt.title("Monte Carlo Simulation: Profit Distribution")
    plt.xlabel("Profit ($)")
    plt.ylabel("Frequency")
    plt.legend()
    plt.tight_layout()
    plt.show()


# plot_distribution()


############################################################
# ESTIMATING PI
############################################################


def throw_needles(num_needles):
    in_circle = 0
    for _ in range(num_needles):
        x = random.random()
        y = random.random()
        if (x*x + y*y)**0.5 <= 1:
            in_circle += 1
    # Counting needles in one quadrant only, so multiply by 4
    return 4 * in_circle / num_needles


# random.seed(3333)
# print("Estimating pi (3.14159) with increasing number of needles:")
# for p in range(1, 8):
#     n = 10**p
#     est = throw_needles(n)
#     print("With ", n, "needles, estimate for pi:", est)


############################################################
# ESTIMATING PI
############################################################


def integrate(fcn, minX, maxX, num_samples = 1_000_000):
    under_curve = 0
    for _ in range(num_samples):
        x = random.uniform(minX, maxX)
        y = fcn(x)
        under_curve += y
    return under_curve / num_samples * (maxX - minX)


def integrate_and_print(fcn, minX, maxX, num_samples = 1_000_000):
    result = integrate(fcn, minX, maxX, num_samples)
    print(f'Integral of {fcn.__name__} from {minX:.0f} to {maxX:2f} =',
          f'{result:2f}')


# random.seed(0)
# integrate_and_print(np.sin, 0, np.pi)
# integrate_and_print(np.sin, 0, 2* np.pi)