import matplotlib.pyplot as plt
import numpy as np

# Define the step function
def step_function(x):
    return 0 if x < 0 else 1

# Define the sigmoid function
def sigmoid_function(x):
    return 1 / (1 + np.exp(-x))

# Generate x values
x_high_res = np.linspace(-10, 10, 1000)  # High resolution for the functions
x_low_res = np.linspace(-10, 10, 71)  # Low resolution for the bars

# Calculate y values for both functions
y_step = [step_function(i) for i in x_high_res]
y_sigmoid = [sigmoid_function(i) for i in x_high_res]

# Plot both functions
plt.plot(x_high_res, y_step, label=r'$\mathbb{1}_{x>y}$')
plt.plot(x_high_res, y_sigmoid, label='F')

# Create a series of vertical bars to represent the area between the two functions
for i in range(len(x_low_res) - 1):
    bar_height = abs(sigmoid_function(x_low_res[i]) - step_function(x_low_res[i]))
    bar_width = x_low_res[i+1]-x_low_res[i] if i != len(x_low_res) - 2 else x_low_res[-1]-x_low_res[-2]
    plt.bar(x_low_res[i], bar_height, bottom=min(step_function(x_low_res[i]), sigmoid_function(x_low_res[i])), width=bar_width, color='grey', align='edge', alpha=0.5)

# Add an annotation for the grey area
plt.annotate('CRPS', xy=(0.75, 0.75), xytext=(-3, 0.7),
             arrowprops=dict(facecolor='black', shrink=0.05))

# Add labels and title
plt.xlabel('x')
#plt.title('Step Function vs Sigmoid Function')
plt.legend()

# Display the plot
plt.show()